[matlab-fortran] translate interpolation

[s_hy]

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 couldn't understand it well. Can anyone here help me to translate interp1 funtion into fortran90...thank you very much.

  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

 

[.Scott]

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)).

 

[s_hy]

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.

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
endt=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
    pcolor(Ez); shading interp;
    %axis off;
  
end

 

[.Scott]

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?

 

[s_hy]

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 explanation in simple english (me not a native speaker). by the way, thank you very much

 

[.Scott]

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?