matlab里如何写fortran,如何在MATLAB中对FORTRAN中输出的数据读入，以及运用

一直使用FORTRAN编译程序，最近想用Matlab里的几个特殊函数，想把这个FORTRAN程序转变成Matlab程序，其实这个程序是由读入数据，梯形公式，格林公式，求导公式，第一类贝塞尔函数(贝塞尔函数)， 第二类贝塞尔函数(诺伊曼函数 )， 第三类贝塞尔函数(汉克尔函数 )组成。想用Matlab，是因为后边的贝塞尔函数想用他们的复数形式，而不知道复数形式的ForTraN怎么编写，因此想直接用MATLaB。求赐教！

Module globdata

implicit none

integer in,jn,kn,im,jm,km,lm,ln

integer unit,key

parameter(im=449,jm=377,km=33,lm=5)

parameter(unit=8)

character*2 num

real*8 re,omer,alfr

real*8,dimension(im,jm) :: x,y,abx,bby,bky,bkx

real*8,allocatable,dimension(:,:,:) :: px,ppy,pp

complex*16,allocatable,dimension(:,:,:,:):: mm

real*8 bes,neu

complex*16 han,grn

complex*16,dimension(im,jm) ::gg

real*8,dimension(im,jm) ::ggy

complex*16 kc

complex*16,parameter:: ii=(0.d0,1.d0)

real*8 coe_4s_1(5),coe_4s_2(5),coe_4c(2),coe_6c(3),coe_5z(6),coe_5f(6)

End Module

program main

use globdata

implicit none

integer i,j,k,l,ys,is,js,xl

real*8 dx,ma

open(unit,file='ini.para')

read(unit,*)omer,alfr,re

read(unit,*)xl,dx

read(unit,*)in,jn,ln,kn

read(unit,*)kc

close(unit)

allocate(mm(in,jn,ln,kn))

allocate(pp(in,jn,kn))

allocate(px(in,jn,kn))

allocate(ppy(in,jn,kn))

open(77,file='pp.plt')

read(77,*)

read(77,*)

do k=2,2

do j=1,jn

do i=1,in

read(77,*)x(i,j),y(i,j),mm(i,j,5,2)

enddo

enddo

enddo

do k=2,2

do j=1,jn

do i=1,in

pp(i,j,k)=mm(i,j,5,2)

enddo

enddo

enddo

call coe_scheme

do k=2,2

do i=1,in

call div_grid_y(pp(i,:,k),ppy(i,:,k),jn)                      call div_grid_y(y(i,:),bky(i,:),jn)

end do

end do

!-------------- --------------------------

ys=270

do k=2,2

do j=ys+1,jn

do i=1,in  !in

!-------------------------------------

do is=1,in

do js=ys,ys

call green(x(i,j)-x(is,ys),y(i,j)-y(is,ys))

gg(is,js)=grn

end do

call div_grid_y(dreal(gg(is,:)),ggy(is,:),jn)

end do

do is=1,in-1

px(i,j,k)=px(i,j,k)+(pp(is,ys,k)*ggy(is,js)/bky(is,js)-ppy(is,ys,k)/bky(is,ys)*gg(is,js)+pp(is+1,ys,k)*ggy(is+1,js)/bky(is+1,js)-ppy(is+1,ys,k)/bky(is+1,ys)*gg(is+1,js))*dx/2.d0

enddo

!-------------------------------------

enddo

enddo

enddo

!----------------------------------------------------------------

open(123,file='px.plt')

write(123,*)'variables="x","y","p"'

write(123,*) 'zone',' i=',in,' j=',jn

do j=1,jn

do i=1,in

write(123,'(3f16.8)') x(i,j),y(i,j),px(i,j,2)

end do

end do

close(123)

end program

!----------the end------------------

!====================法向求导=================================

subroutine div_grid_y(f,df,mn)

use globdata,only: coe_4s_1,coe_4s_2,coe_4c,coe_6c

implicit none

integer*4 mn,i

real*8 f(mn),df(mn)

i=1

df(i  )=coe_4s_1(1)*f(i)+coe_4s_1(2)*f(i+1)+coe_4s_1(3)*f(i+2)+coe_4s_1(4)*f(i+3)+coe_4s_1(5)*f(i+4)

i=1

df(i+1)=coe_4s_2(1)*f(i)+coe_4s_2(2)*f(i+1)+coe_4s_2(3)*f(i+2)+coe_4s_2(4)*f(i+3)+coe_4s_2(5)*f(i+4)

i=3

df(i  )=coe_4c(1)*(f(i+1)-f(i-1))+coe_4c(2)*(f(i+2)-f(i-2))

do I=4,mn-3

df(i) = coe_6c(1) * (f(i+1) - f(i-1)) + coe_6c(2) * (f(i+2) - f(i-2)) + coe_6c(3) * (f(i+3) - f(i-3))

enddo

i=mn-2

df(i  )=coe_4c(1)*(f(i+1)-f(i-1))+coe_4c(2)*(f(i+2)-f(i-2))

i=mn

df(i-1)=-(coe_4s_2(1)*f(i)+coe_4s_2(2)*f(i-1)+coe_4s_2(3)*f(i-2)+coe_4s_2(4)*f(i-3)+coe_4s_2(5)*f(i-4))

i=mn

df(i  )=-(coe_4s_1(1)*f(i)+coe_4s_1(2)*f(i-1)+coe_4s_1(3)*f(i-2)+coe_4s_1(4)*f(i-3)+coe_4s_1(5)*f(i-4))

return

end subroutine

subroutine coe_scheme

use globdata

implicit none

coe_4s_1(1)=-25.d0/12.d0 ;coe_4s_1(2)=48.d0/12.d0 ;coe_4s_1(3)=-36.d0/12.d0

coe_4s_1(4)=16.d0/12.d0 ;coe_4s_1(5)=-3.d0/12.d0

coe_4s_2(1)=-3.d0/12.d0 ;coe_4s_2(2)=-10.d0/12.d0 ;coe_4s_2(3)=18.d0/12.d0

coe_4s_2(4)=-6.d0/12.d0 ;coe_4s_2(5)=1.d0/12.d0

coe_4c(1)=2.d0/3.d0 ;coe_4c(2)=-1.d0/12.d0

coe_6c(1)=3.d0/4.d0 ;coe_6c(2)=-3.d0/20.d0 ;coe_6c(3)=1.d0/60.d0

return

end subroutine

!*************************************************************

!           第一类贝塞尔函数:贝塞尔函数

!*************************************************************/

subroutine Bessel(n,xx)

use globdata

implicit none

integer*4 i,m,n

real*8 xx,yy,zz,t,p,q,s,b0,b1

real*8 a(1:6),b(1:6),c(1:6),d(1:6)

real*8 e(1:5),f(1:5),g(1:5),h(1:5)

a(1)=57568490574.0

a(2)=-13362590354.0

a(3)=651619640.7

a(4)=-11214424.18

a(5)=77392.33017

a(6)=-184.9052456

b(1)=57568490411.0

b(2)=1029532985.0

b(3)=9494680.718

b(4)=59272.64853

b(5)=267.8532712

b(6)=1.0

c(1)=72362614232.0

c(2)=-7895059235.0

c(3)=242396853.1

c(4)=-2972611.439

c(5)=15704.4826

c(6)=-30.16036606

d(1)=144725228443.0

d(2)=2300535178.0

d(3)=18583304.74

d(4)=99447.43394

d(5)=376.9991397

d(6)=1.0

e(1)=1.0

e(2)=-0.1098628627d-2

e(3)=0.2734510407d-4

e(4)=-0.2073370639d-5

e(5)=0.2093887211d-6

f(1)=-0.1562499995d-1

f(2)=0.1430488765d-3

f(3)=-0.6911147651d-5

f(4)=0.7621095161d-6

f(5)=-0.934935152d-7

g(1)=1.0

g(2)=0.183105d-2

g(3)=-0.3516396496d-4

g(4)=0.2457520174d-5

g(5)=-0.240337019d-6

h(1)=0.4687499995d-1

h(2)=-0.2002690873d-3

h(3)=0.8449199096d-5

h(4)=-0.88228987d-6

h(5)=0.105787412d-6

t=dabs(xx);

if(n.ne.1) then

if(t.lt.8.d0) then

yy=t**2

p=a(6)

q=b(6)

do i=5,1,-1

p=p*yy+a(i)

q=q*yy+b(i)

end do

bes=p/q

else

zz=8.d0/t      !!

yy=zz**2

p=e(5)

q=f(5)

do i=4,1,-1

p=p*yy+e(i)

q=q*yy+f(i)

end do

s=t-0.785398164

p=p*dcos(s)-zz*q*dsin(s)

bes=p*sqrt(0.636619772/t)                        !mark gg

end if

end if

end subroutine

!/*************************************************************

!            第二类贝塞尔函数:诺伊曼函数

!*************************************************************/

subroutine Neumann(n,xx)

use globdata

implicit none

integer*4 i,n

real*8 xx,yy,zz,p,q,s,b0,b1

real*8 a(1:6),b(1:6),c(1:6),d(1:7)

real*8 e(1:5),f(1:5),g(1:5),h(1:5)

a(1)=-2.957821389d9

a(2)=7.062834065d9

a(3)=-5.123598036d8

a(4)=1.087988129d7

a(5)=-8.632792757d4

a(6)=2.284622733d2

b(1)=4.0076544269d10

b(2)=7.452499648d8

b(3)=7.189466438d6

b(4)=4.74472647d4

b(5)=2.261030244d2

b(6)=1.d0

c(1)=-4.900604943d12

c(2)=1.27527439d12

c(3)=-5.153438139d10

c(4)=7.349264551d8

c(5)=-4.237922726d6

c(6)=8.511937935d3

d(1)=2.49958057d13

d(2)=4.244419664d11

d(3)=3.733650367d9

d(4)=2.245904002d7

d(5)=1.02042605d5

d(6)=3.549632885d2

d(7)=1.0

e(1)=1.0

e(2)=-0.1098628627d-2

e(3)=0.2734510407d-4

e(4)=-0.2073370639d-5

e(5)=0.2093887211d-6

f(1)=-0.1562499995d-1

f(2)=0.1430488765d-3

f(3)=-0.6911147651d-5

f(4)=0.7621095161d-6

f(5)=-0.934935152d-7

g(1)=1.0

g(2)=0.183105d-2

g(3)=-0.3516396496d-4

g(4)=0.2457520174d-5

g(5)=-0.240337019d-6

h(1)=0.4687499995d-1

h(2)=-0.2002690873d-3

h(3)=0.8449199096d-5

h(4)=-0.88228987d-6

h(5)=0.105787412d-6

if(xx==0.d0) then

neu=-1.0d70

endif

if(n.ne.1) then

if(xx.lt.8.d0) then

yy=xx*xx

p=a(6)

q=b(6);

do i=5,1,-1

p=p*yy+a(i)

q=q*yy+b(i)

end do

call Bessel(0,xx)

neu=p/q+0.636619772*bes*dlog(xx)

else                             !!

zz=8.0/xx

yy=zz*zz

p=e(5)

q=f(5);

do i=4,1,-1

p=p*yy+e(i)

q=q*yy+f(i)

end do

s=xx-0.785398164

p=p*dsin(s)+zz*q*dcos(s)

neu=p*sqrt(0.636619772/xx)

end if

end if

end subroutine

!/*************************************************************

!            第三类贝塞尔函数:汉克尔函数

!*************************************************************/

subroutine Hankel(n,xx)

use globdata

implicit none

integer*4 n

real*8 xx

call Bessel(n,xx)

call Neumann(n,xx)

han=bes-ii*neu

end subroutine

!====================格林函数================================

subroutine Green(xx,yy)

use globdata

implicit none

real*8 xx,yy

complex*16 p1,p2,ma

ma=1.5

p1=exp((ii*xx*kc*ma)/(1-ma**2))/(4*ii*dsqrt(1-ma**2))

p2=dsqrt(xx**2/((1-ma**2))+yy**2)*dsqrt(1/(1-ma**2))*kc

call Hankel(0,p2)

grn=han*p1

end subroutine