Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[matlab-fortran] translate interpolation

  1. Feb 10, 2014 #1
    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 (Text):

      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

     
     
  2. jcsd
  3. Feb 10, 2014 #2
    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)).
     
  4. Feb 11, 2014 #3
    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 (Text):

    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

     
     
  5. Feb 11, 2014 #4
    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?
     
  6. Feb 11, 2014 #5
    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
     
  7. Feb 11, 2014 #6
    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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [matlab-fortran] translate interpolation
  1. Fortran 90 Vs Matlab (Replies: 5)

Loading...