Converting working Fortran code to MATLAB

Hello,

I am trying to convert a Fortran code into a Matlab code. Most of the fortran code consists of matrices / arrays such as T(i,k), u(i,k), x(i), dpdx(x).

I'm wondering what is the best way to replicate those arrays in matlab. For instance the matrix T(i,k) in the for loop, shouldn't I use the zeros command instead of writing T(i,k) = 0 which is part of my initial condtions for my problem??

I'm a programming newbie so some code modifications are not as obvious to me.

One piece of the Fortran code is as following:

"

imin = 0
imax = 100
kmin = 0
kmax = 100
%%%
do i = imin, imax, 1
   x(i) = xmin + real (i - imin, kr) / real (imax - imin, kr) * (xmax - xmin)
   h(i) = hm + x(i) ** 2.0_kr / (2.0_kr * Rad)
   dpdx(i) = 6.0_kr * visc / h(i) ** 3.0_kr * ((U_ls + U_us) * h(i) + cons1)
          do k = kmin, kmax, 1
          z(k) = zmin + real (k - kmin, kr) / real (kmax - kmin, kr) * (h(i) - zmin)
          T(i,k) = 0.0_kr
           A1 = dpdx(i) / (2.0_kr * visc)
          A2 = (- dpdx(i) / (2.0_kr * visc) * h(i) + (U_us - U_ls) / h(i))
          u(i,k) = A1 * z(k) ** 2.0_kr + A2 * z(k) + U_ls
          dudz(i,k) = 2.0_kr * A1 * z(k) + A2
    enddo
enddo

%%%

My matlab equivalent currently looks like this:

for i = 1:imax
    x(i)= xmin + (i - imin) / (imax - imin) * (xmax - xmin);
    h(i) = hmin + x(i).^2 / (2 * Rad);
    dpdx = 6 * visc /(h(i).^3) * ((U_ls * U_us) * h(i) + cons1);
    for k = kmin:kmax
        z(i) = zmin + (k - kmin) / (kmax - kmin) * (h(i) - zmin);
        T(i,k) = 0
        A1 = dpdx / (2 * visc);
        A2 = (-dpdx / (2 * visc)) * h(i) + (U_us - U_ls) / h(i);
        u(i,k) = A1 * z(i).^2 + A2 * z(i) + U_ls;
        dudz(i,k) = 2 * A1 * z(i) + A2;
    end
end

I'm not sure about the i and k arrays in the for loop which in the original code run from 0 to 100 but in matlab they go from 1 to 100, how do I correct this?

Please ignore the other variables I listed like hmin, visc, Rad etc. I would really appreciate it if someone could tell me that the structure of my matlab code is correct and would yield me the same output as the original Fortran code.

[KSSV]

In MATLAB the indices should be posititve integers. You can run the loop from i = 1 to 101.

Answer 1

[Ben Barrowes]

Here is your fortran code cleaned up and compilable:

program loops

integer, parameter :: kr = 8
real x(0:100), h(0:100), dpdx(0:100), z(0:100)
real T(0:100,0:100),u(0:100,0:100),dudz(0:100,0:100)

imin = 0
imax = 100
kmin = 0
kmax = 100

do i = imin, imax, 1
 x(i) = xmin + real (i - imin, kr) / real (imax - imin, kr) * (xmax - xmin)
 h(i) = hm + x(i) ** 2.0_kr / (2.0_kr * Rad)
 dpdx(i) = 6.0_kr * visc / h(i) ** 3.0_kr * ((U_ls + U_us) * h(i) + cons1)
 do k = kmin, kmax, 1
  z(k) = zmin + real (k - kmin, kr) / real (kmax - kmin, kr) * (h(i) - zmin)
  T(i,k) = 0.0_kr
  A1 = dpdx(i) / (2.0_kr * visc)
  A2 = (- dpdx(i) / (2.0_kr * visc) * h(i) + (U_us - U_ls) / h(i))
  u(i,k) = A1 * z(k) ** 2.0_kr + A2 * z(k) + U_ls
  dudz(i,k) = 2.0_kr * A1 * z(k) + A2
 enddo
enddo


end program loops

And here is the result when I run it through f2matlab:

function loops(varargin)
 persistent a1 a2 cons1 dpdx dudz firstCall h hm i imax imin k kmax kmin kr rad t u u_ls u_us visc x xmax xmin z zmin 
 ; 
 if isempty(firstCall),firstCall=1;end; 

 if firstCall;
  kr = 8;
  a1=0;
  a2=0;
  cons1=0;
  dpdx=zeros(1,100+1);
  dudz=zeros(100+1,100+1);
  h=zeros(1,100+1);
  hm=0;
  i=0;
  imax=0;
  imin=0;
  k=0;
  kmax=0;
  kmin=0;
  rad=0;
  t=zeros(100+1,100+1);
  u=zeros(100+1,100+1);
  u_ls=0;
  u_us=0;
  visc=0;
  x=zeros(1,100+1);
  xmax=0;
  xmin=0;
  z=zeros(1,100+1);
  zmin=0;
 end
 firstCall=0;

 imin = 0;
 imax = 100;
 kmin = 0;
 kmax = 100;

 for i = imin: 1: imax
  x(i+1) = xmin + real(i-imin)./real(imax-imin).*(xmax-xmin);
  h(i+1) = hm + x(i+1).^2.0./(2.0.*rad);
  dpdx(i+1) = 6.0.*visc./h(i+1).^3.0.*((u_ls+u_us).*h(i+1)+cons1);
  for k = kmin: 1: kmax
   z(k+1) = zmin + real(k-kmin)./real(kmax-kmin).*(h(i+1)-zmin);
   t(i+1,k+1) = 0.0;
   a1 = dpdx(i+1)./(2.0.*visc);
   a2 =(-dpdx(i+1)./(2.0.*visc).*h(i+1)+(u_us-u_ls)./h(i+1));
   u(i+1,k+1) = a1.*z(k+1).^2.0 + a2.*z(k+1) + u_ls;
   dudz(i+1,k+1) = 2.0.*a1.*z(k+1) + a2;
  end
  k = fix(kmax+ 1);
 end
 i = fix(imax+ 1);

end %program loops

This could be further optimized and cleaned up in matlab, but this gets you working code you can profile and tinker with.

PM me if you have more code or further troubles. Posting larger codes here can get unwieldy.