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.
1 个评论
KSSV
2020-6-10
In MATLAB the indices should be posititve integers. You can run the loop from i = 1 to 101.
采纳的回答
Ben Barrowes
2020-6-11
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
