How to convert fortran code to matlab?
5 views (last 30 days)
Show older comments
PROGRAM IMPLISIT
C
C *** PROGRAM PENYELESAIAN PERSAMAAN HYPERBOLIK DENGAN METODA
IMPLICIT
C
PARAMETER( IMX=1000 )
DIMENSION U(IMX),UL(IMX),UB(IMX),XP(IMX),UP(IMX,100),TP(IMX),
: A(IMX),B(IMX),C(IMX),D(IMX)
C
C-------------------------------------------------------------------
-----------
C
C UL(I) = HARGA LAMA VARIABLE PADA TITIK I
C UB(I) = HARGA BARU VARIABLE PADA TITIK I
V-12
C UP(I,K) = BEBERAPA HARGA VARIABLE PADA TITIK I DISIMPAN UNTUK PLOTTING
C
C------------------------------------------------------------------------------
C
OPEN(2,FILE='IM1.CON')
OPEN(3,FILE='IM1.OUT')
C
C *** INPUT DATA DAN PARAMETER
C
IM = 801
XL = 0
XR = 20
DX = (XR - XL )/(IM-1)
DT = 0.5
TMAX = 101
CR = 0.1
THETA2 = 0.5
THETA1 = 1.0 - THETA2
IPRINT = 20
CN = CR*DT/DX
C
C *** KOORDINAT NODE
C
DO I=1,IM
XP(I) = XL + (I-1)*DX
ENDDO
C
C *** HARGA AWAL (INITIAL VALUES)
C
DO I=1,IM
IF( XP(I).GT.2.5.AND.XP(I).LE.3.0 )THEN
UB(I) = ( XP(I) - 2.5 )
ELSEIF( XP(I).GT.3.0.AND.XP(I).LE.3.5 )THEN
UB(I) = ( 3.5 - XP(I) )
ELSE
UB(I) = 0.0
ENDIF
ENDDO
C
C *** TETAPKAN HARGA UB(I) PADA SYARAT BATAS EXTERNAL
C
UB(1) = 0
UB(IM) = 0
C
C *** LOOP ITERASI FTCS
C
ITER = 0
IP = IPRINT
TIME = 0
K = 0
100 CONTINUE
V-13
C
C ***** SIMPAN HARGA UB(I) PADA TIME STEP N
C
DO I=1,IM
UL(I) = UB(I)
ENDDO
C
C **** SIMPAN HASIL UNTUK PLOTTING
C
IF( IP.EQ.IPRINT )THEN
IP = 0
K = K + 1
DO I=1,IM
UP(I,K) = UB(I)
TP(K) = TIME
ENDDO
WRITE(*,'(I10,3F10.3)') ITER,TIME,CN
WRITE(2,'(I10,100F10.3)') ITER,TIME,(UB(I),I=1,IM)
ENDIF
ITER = ITER + 1
TIME = TIME + DT
IP = IP + 1
C
C ***** PEMBENTUKAN MATRIK TRIDIAGONAL UNTUK MEMPERBAHARUI HARGA UB(I) PADA TIME STEP N+1
C
DO I=2,IM-1
II = I - 1
IF( I.EQ.2 )THEN
A(II) = 0
B(II) = 1.0
C(II) = THETA1*CN/2
D(II) = UL(I) + THETA2*CN/2*( UL(I+1) - UL(I-1) ) +
: THETA1*CN/2*UL(I-1)
ELSEIF( I.EQ.IM-1 )THEN
A(II) = -THETA1*CN/2
B(II) = 1.0
C(II) = 0
D(II) = UL(I) - THETA2*CN/2*( UL(I+1) - UL(I-1) ) -
: THETA1*CN/2*UL(I+1)
ELSE
A(II) = -THETA1*CN/2
B(II) = 1.0
C(II) = THETA1*CN/2
D(II) = UL(I) - THETA2*CN/2*( UL(I+1) - UL(I-1) )
ENDIF
ENDDO
C
C ***** SELESAIKAN DENGAN TRIDIAGONAL SOLVER
C
CALL TRIDAG(A,B,C,D,U,IM-2)
DO I=1,IM-2
UB(I+1) = U(I)
V-14
ENDDO
IF( TIME.LT.TMAX )GOTO 100
C
C *** PRINT HASIL
C
WRITE(3,'(10X,999F10.3)') (TP(J),J=1,K)
DO I=1,IM
WRITE(3,'(1000F10.3)') XP(I),(UP(I,J),J=1,K)
ENDDO
STOP
END
C
C--------------------------< SUBROUTINE THOMAS >-------------------------------
C
SUBROUTINE tridag(a,b,c,r,u,n)
PARAMETER (NMAX=1000)
REAL a(n),b(n),c(n),r(n),u(n),bet,gam(NMAX)
IF( b(1).eq.0.) PAUSE 'tridag: rewrite equations'
c
c *** Forward Sweep
c
bet=b(1)
u(1)=r(1)/bet
do 11 j=2,n
gam(j)=c(j-1)/bet
bet=b(j)-a(j)*gam(j)
if(bet.eq.0.)pause 'tridag failed'
u(j)=(r(j)-a(j)*u(j-1))/bet
11 continue
c
c *** Backward Sweep
c
do 12 j=n-1,1,-1
u(j)=u(j)-gam(j+1)*u(j+1)
12 continue
RETURN
END
Answers (1)
Ben Barrowes
on 12 Mar 2022
You should try my file exchange program, f2matlab. But you would have to convert it fortran90 style first, including refactoring the "goto" in there.
But I put it through my version of f2matlab after cleaning up the code a bit. I had to clean it up before it would compile, because when you copied and pasted this code, you lost all the indenting. And for example the ":" at the beginning of the line in your code is the continuation character. One thing I didn't get in this code was the "V-12", and "V-13", and "V-14". They almost look like page numbers. I commented them out.
Here is your converted code, which I admit I have not tested... nor can I without input files, etc. But there are support functions as well like writeFmt which I am including as an attachment.
function implisit(varargin)
global unit2fid;
if ~isempty(unit2fid);
unit2fid={};
end
imx = 1000;
a=zeros(1,imx);
b=zeros(1,imx);
c=zeros(1,imx);
cn=0.0;
cr=0.0;
d=zeros(1,imx);
dt=0.0;
dx=0.0;
i=0;
ii=0;
im=0;
ip=0;
iprint=0;
iter=0;
j=0;
k=0;
theta1=0.0;
theta2=0.0;
time_fv=0.0;
tmax=0.0;
tp=zeros(1,imx);
u=zeros(1,imx);
ub=zeros(1,imx);
ul=zeros(1,imx);
up=zeros(imx,100);
xl=0.0;
xp=zeros(1,imx);
xr=0.0;
% *** PROGRAM PENYELESAIAN PERSAMAAN HYPERBOLIK DENGAN METODA IMPLICIT
%------------------------------------------------------------------------------
%
% UL(I) = HARGA LAMA VARIABLE PADA TITIK I
% UB(I) = HARGA BARU VARIABLE PADA TITIK I
% V-12
% UP(I,K) = BEBERAPA HARGA VARIABLE PADA TITIK I DISIMPAN UNTUK PLOTTING
%
%------------------------------------------------------------------------------
%
if exist(strtrim('IM1.CON'))~=2;
fclose(fopen(strtrim('IM1.CON'),'w'));
end
thismlfid=fopen(strtrim('IM1.CON'),'r+');
unit2fid{end+1,1}=2;
unit2fid{end,2}=thismlfid;
unit2fid{end,3}=0;
unit2fid{end,4}=strtrim('IM1.CON');
unit2fid{end,5}=0;
unit2fid{end,6}=0;
% unit2fid maps fortran unit numbers to matlab fid numbers
if exist(strtrim('IM1.OUT'))~=2;
fclose(fopen(strtrim('IM1.OUT'),'w'));
end
thismlfid=fopen(strtrim('IM1.OUT'),'r+');
unit2fid{end+1,1}=3;
unit2fid{end,2}=thismlfid;
unit2fid{end,3}=0;
unit2fid{end,4}=strtrim('IM1.OUT');
unit2fid{end,5}=0;
unit2fid{end,6}=0;
% unit2fid maps fortran unit numbers to matlab fid numbers
%
% *** INPUT DATA DAN PARAMETER
%
im = 801;
xl = 0;
xr = 20;
dx =(xr-xl)./(im-1);
dt = 0.5;
tmax = 101;
cr = 0.1;
theta2 = 0.5;
theta1 = 1.0 - theta2;
iprint = 20;
cn = cr.*dt./dx;
%
% *** KOORDINAT NODE
%
for i = 1: im
xp(i) = xl +(i-1).*dx;
end
i = fix(im+1);
%
% *** HARGA AWAL (INITIAL VALUES)
%
for i = 1: im
if(xp(i) > 2.5 && xp(i) <= 3.0)
ub(i) =(xp(i)-2.5);
elseif(xp(i) > 3.0 && xp(i) <= 3.5)
ub(i) =(3.5-xp(i));
else;
ub(i) = 0.0;
end
end
i = fix(im+1);
%
% *** TETAPKAN HARGA UB(I) PADA SYARAT BATAS EXTERNAL
%
ub(1) = 0;
ub(im) = 0;
%
% *** LOOP ITERASI FTCS
%
iter = 0;
ip = fix(iprint);
time_fv = 0;
k = 0;
while (1);
% V-13
%
% ***** SIMPAN HARGA UB(I) PADA TIME STEP N
%
for i = 1: im
ul(i) = ub(i);
end
i = fix(im+1);
%
% **** SIMPAN HASIL UNTUK PLOTTING
%
if(ip == iprint)
ip = 0;
k = fix(k + 1);
for i = 1: im
up(i,k) = ub(i);
tp(k) = time_fv;
end
i = fix(im+1);
[writeErrFlag,wfso]=writeFmt(1,['%10i',repmat('%10.3f',1,3)],'iter','time_fv','cn');
eval(wfso);
[writeErrFlag,wfso]=writeFmt(2,['%10i',repmat('%10.3f',1,100)],'iter','time_fv',{'ub(i)','i',1,1,im});
eval(wfso);
end
iter = fix(iter + 1);
time_fv = time_fv + dt;
ip = fix(ip + 1);
%
% ***** PEMBENTUKAN MATRIK TRIDIAGONAL UNTUK MEMPERBAHARUI HARGA UB(I) PADA TIME STEP N+1
%
for i = 2: im - 1
ii = fix(i - 1);
if(i == 2)
a(ii) = 0;
b(ii) = 1.0;
c(ii) = theta1.*cn./2;
d(ii) = ul(i) + theta2.*cn./2.*(ul(i+1)-ul(i-1)) + theta1.*cn./2.*ul(i-1);
elseif(i == im-1)
a(ii) = -theta1.*cn./2;
b(ii) = 1.0;
c(ii) = 0;
d(ii) = ul(i) - theta2.*cn./2.*(ul(i+1)-ul(i-1)) - theta1.*cn./2.*ul(i+1);
else;
a(ii) = -theta1.*cn./2;
b(ii) = 1.0;
c(ii) = theta1.*cn./2;
d(ii) = ul(i) - theta2.*cn./2.*(ul(i+1)-ul(i-1));
end
end
i = fix(im - 1+1);
%
% ***** SELESAIKAN DENGAN TRIDIAGONAL SOLVER
%
[a,b,c,d,u,dumvar6]=tridag(a,b,c,d,u,im-2);
for i = 1: im - 2
ub(i+1) = u(i);
% V-14
end
i = fix(im - 2+1);
if(time_fv < tmax)
continue;
end
%
% *** PRINT HASIL
%
[writeErrFlag,wfso]=writeFmt(3,['%10x',repmat('%10.3f',1,999)],{'tp(j)','j',1,1,k});
eval(wfso);
for i = 1: im
[writeErrFlag,wfso]=writeFmt(3,[repmat('%10.3f',1,1000)],'xp(i)',{'up(i,j)','j',1,1,k});
eval(wfso);
end
i = fix(im+1);
tempBreak=1;
break;
end
closeAllOpenFiles(unit2fid);
end %program implisit
3 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!