# [matlab-fortran] translate interpolation

• MATLAB

## Main Question or Discussion Point

hi all

I am new learner on fortran. as a part of practice, i am trying to translate below equation into fortan. I have read about interpolation, but couldnt understand it well. Can anyone here help me to translate interp1 funtion into fortran90....thank you very much.

Code:
  for i=ja:jb
Ez(ia,i) = Ez(ia,i) + cb1(ia,i)*interp1(1:length(Hinc),Hinc,PLT(ia-1,i))*CosA;
Ez(ib,i) = Ez(ib,i) - cb1(ib,i)*interp1(1:length(Hinc),Hinc,PLT(ib,i))*CosA;
end

for i=ia:ib
Ez(i,ja) = Ez(i,ja) + cb1(i,ja)*interp1(1:length(Hinc),Hinc,PLT(i,ja-1))*SinA;
Ez(i,jb) = Ez(i,jb) - cb1(i,ja)*interp1(1:length(Hinc),Hinc,PLT(i,jb))*SinA;
end

Related MATLAB, Maple, Mathematica, LaTeX News on Phys.org
.Scott
Homework Helper
Since you're using only one form of the interp1 function and the first parameter is always 1:length(Hinc), you can save time by only coding you're Fortran function for that case.

So this is what the interp1 is doing in this instance and what your Fortran function will need to do:
* That first argument is simply making Hinc act as a array with 1-based indexing.
* Substitute that first argument with an integer argument that specifies the length of array Hinc. Call it NMAX.
* Call that third argument "xq". It's the "PLT(...)" in each call to interp1.
* Function interp1 is using xq as an index into vector (array) Hinc.
* However, xq is not necessarily an integer value, so compute NLO=FLOOR(xq) and NHI=CEILING(xq).
* If NLO is less than 1 or NHI is greater than NMAX, you have an error condition.
* Else if NLO equals NHI, return Hinc(NLO).
* Else return this linear interpolation: (Hinc(NLO)*(NHI-xq)) + (Hinc(NHI)*(xq-NLO)).

dear scott,

Thank you for your tips, but I am sorry I couldn't digest it well. I put here all the codes that I used for this practice. Here we can see Hinc in 1D array as well as PLT in 2d array. Hope you can help me to translate it... thank you very much for your time.

Code:
clear all
close all

c0 =2.99792458e8;
muz = 4*pi*1E-7;
epsz=1.0/(c0*c0*muz);       %permittivity of free space

f0 = 45.0*10^9;
lambda = c0/f0;

S = 1; %%Courant Time Stability 1=="magic time step
S2D= S/sqrt(2);
Nx = 40; %Number of times the grid cell will be divided by lambda
Ny = Nx;
dx = lambda/Nx; %length of cell in x dir
dy = lambda/Ny; %length of cell in y dir

N=100; %Number of time steps

NI = 40; %Number of cells
NJ = NI;

NIi = max([NI NJ])*2; % This is to prevent reflections from the boundary at NI
Einc = zeros([1 NIi]);
Hinc = zeros([1 NIi]);

%%Connecting Boundary Setup
ia = 5;
ja = ia;
ib = NI-ia-1;
jb = NJ-ja-1;

%***********************************************************************
%     Material parameters
%***********************************************************************
eps=1;
sig=0;
mur=1;
sim=0;

%% Set the whole space to Vacumn
MaterialSpace = ones(NI,NJ);

%**********************************************************************

Ez = zeros(NI,NJ);
Hx = Ez;
Hy = Ez;

dxy = sqrt((dx)^2+ (dy)^2);
dt = S*dxy/(2*c0); %delta time step

%***********************************************************************
%     Updating coefficients
%***********************************************************************
for i=1:NI;
for j=1:NJ;
eaf=dt*sig/(2.0*epsz*eps);
ca(i,j)=(1.0-eaf)/(1.0+eaf);
cb1(i,j)=dt/epsz/eps/dx./(1.0+eaf);
cb2(i,j)=dt/epsz/eps/dy/(1.0+eaf);
haf =dt*sim/(2.0*muz*mur);
da(i,j)=(1.0-haf)/(1.0+haf);
db1(i,j)=dt/muz/mur/dx/(1.0+haf);
db2(i,j)=dt/muz/mur/dy/(1.0+haf);
end
end
E0 = 1;

source_angle = 45;

%% Connecting Region Activation Constants
m0 = 4;

%% Make Sure everythings adds up nicley
CosA = cos(source_angle*pi/180);
SinA = sin(source_angle*pi/180);

% Determine the Courant Stability factor for the 1D Auxillary case to zero
% out the phase velocity difference between a 1D auxiallary wave equation
% which astc as if it is on axis and an off axis 2D update

S3=sqrt((CosA^4+SinA^4));

% grid size for the 1D auxiallary space
dx1=S3*dx;

%update constants for 1D auxillary
Cb = dt/dx1/epsz;        %Cb
Db = dt/dx1/muz;         %Db

% TODO
% UPDATE THIS FOR RECTANGULAR CELLS
% Find the distance from the Plane wave origin to the connectors (eq
% 5.14 Taflove
for i=1:NI
for j=1:NJ
OffSet = round(sqrt((ia*CosA)^2+(ja*SinA)^2));
PLT(i,j)=1/S3*(CosA*(i-ia)+SinA*(j-ja))+OffSet;
end
end

t=0;

for n=1:200
fprintf('Step %d of %d',n,N);

%% Incident Source

source = E0*sin(2*pi*f0*n*dt);

% After a sufficient amount of time zero out the incident wave
if n == NJ+200
Hinc = zeros([1 NIi]);
Einc = zeros([1 NIi]);

end

%%Boundarry Condition at front
rbc1 = Einc(2);
rbc=Einc(NIi-1);

for j=2:NIi
Einc(j) = Einc(j) + Cb*(Hinc(j-1)-Hinc(j));
end

%%boundary Condition at end
Einc(NIi)=rbc;
Einc(1) = rbc1;

%inject the a hard source into the incident field
Einc(m0) = source;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Update all E Fields
for i=2:(NI-1)
for j=2:(NJ-1)
Ez(i,j) = ca(i,j)*Ez(i,j) + cb1(i,j).*...
(Hy(i,j)-Hy(i-1,j))-cb2(i,j)*(Hx(i,j)-Hx(i,j-1));
end
end

%% Update E fields along front and back faces (x==0 and x==jb)
%%this updates the boundary between the total and scattered field
%%region
for i=ja:jb
Ez(ia,i) = Ez(ia,i) + cb1(ia,i)*interp1(1:length(Hinc),Hinc,PLT(ia-1,i))*CosA;          %left
Ez(ib,i) = Ez(ib,i) - cb1(ib,i)*interp1(1:length(Hinc),Hinc,PLT(ib,i))*CosA;              %right
end

for i=ia:ib
Ez(i,ja) = Ez(i,ja) + cb1(i,ja)*interp1(1:length(Hinc),Hinc,PLT(i,ja-1))*SinA;         %front
Ez(i,jb) = Ez(i,jb) - cb1(i,ja)*interp1(1:length(Hinc),Hinc,PLT(i,jb))*SinA;             %back
end

%% Update H Field
t = t+dt/2;

%%Update Incident Field
for j=1:NIi-1
Hinc(j) = Hinc(j) + Db*(Einc(j)-Einc(j+1));
end

%%Update Hx for Space
for i=2:(NI-1)
for j=1:(NJ-1)
Hx(i,j) = da(i,j)*Hx(i,j) + db2(i,j)*((-Ez(i,j+1)+Ez(i,j)));
end
end

%%Update Scattered/Total Field Hx Connectors
for i=ia:ib
Hx(i,ja-1) = Hx(i,ja-1)+ db1(i,ja-1)*interp1(1:length(Einc),Einc,PLT(i,ja));
Hx(i,jb) = Hx(i,jb) - db1(i,jb)*interp1(1:length(Einc),Einc,PLT(i,jb));

end

%Update Hy for Space
for i=1:(NI-1)
for j=2:(NJ-1)
Hy(i,j) = da(i,j)*Hy(i,j) + db1(i,j).*((-Ez(i,j)+Ez(i+1,j)));
end
end

%Update Scattered/Total Field Hy Connectors
for j=ja:jb
Hy(ia-1,j) = Hy(ia-1,j)- db1(ia-1,j)*interp1(1:length(Einc),Einc,PLT(ia,j));
Hy(ib,j) = Hy(ib,j)+ db1(ib,j)*interp1(1:length(Einc),Einc,PLT(ib,j));
end

Trans (n) = sum(Ez(:,jb+1));
Refl(n) = sum(Ez(:,ja-1));
Reference(n) = Einc(m0);

t = t+dt/2;
drawnow

% resolution setting
%axis off;

end

.Scott
Homework Helper
Hinc is a 1D array and you are passing the entire array to interp1.
PLT is a 2D array, but that doesn't matter to interp1 because you are only passing one element of that array to it.
You original question was for assistance in converting interp1 to Fortran. Are you now asking for assistance in converting the whole thing? Have you attempted to write any of the Fortran code yet?

no, i am not asking you to translate the whole codes. sorry for this. i just want to show you the Hinc and PLT and some explaination in simple english (me not a native speaker). by the way, thank you very much

.Scott
Homework Helper
The only reason I asked was that I was hoping you'd be far enough along to show Fortran Code.
I was also thinking that you were more familiar with MatLab than Fortran - but now I'm thinking that you may need assistance with the MatLab part as well. Is this the case?