How to convert fortran code to matlab?

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

[KSSV]

You have to convert it manually.

Answer 1

[Ben Barrowes]

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

[amir djaouadi]

Ben Barrowes sir please i need i have .for i would to change to .m for trnsys

[amir djaouadi]

Ben Barrowes [e-mail]

[Ben Barrowes]

Hi Amir. I will email you.