[matlab-fortran] translate interpolation

  • MATLAB
  • Thread starter s_hy
  • Start date
  • #1
61
0

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
 

Answers and Replies

  • #2
.Scott
Homework Helper
2,402
828
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)).
 
  • #3
61
0
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
    pcolor(Ez); shading interp;
    %axis off;
  
end
 
  • #4
.Scott
Homework Helper
2,402
828
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?
 
  • #5
61
0
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
 
  • #6
.Scott
Homework Helper
2,402
828
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?
 

Related Threads for: [matlab-fortran] translate interpolation

  • Last Post
Replies
0
Views
7K
Replies
4
Views
4K
  • Last Post
Replies
9
Views
34K
  • Last Post
Replies
0
Views
5K
Replies
1
Views
2K
Replies
1
Views
12K
Replies
1
Views
6K
Replies
2
Views
4K
Top