一直使用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