fortran 笔记

实际上是从test.f中直接拷贝过来的:)
!该文档是用来试验各种功能的。可编译!
include 'gammln.for'
program test1
real*8 gammln,x
x=.8
!怎样判断读一个文件到尾了呢?
! read(123,*,end=124)x即可,到了尾端自动跳到标号为124的语句继续执行.
!如果不用open()语句，将自动生成一个文件放东西
write(33,*)'fort.33'
!怎样在循环中中止一个循环，且不用goto语句
!goto语句用来实现跳出循环的功能, 最好不要在让它表示其它跳转, 否则会晕头转向的.
do i=1,10
if(i>5)then
exit !goto 11 !"exit" is not in the f77 standard
endif
enddo
11 print*,'i',i
!怎样做到注释一块的功能，改成“.false.”即可
!其中有输出单引号的写法
if(.true.) then
print*,'''if(.false.) then'' is the same as /*... */ in c, ',
& 'for notation'
endif
!验证是否支持符号的大于、等于号
if ( x < 1.)then
print *,'good! it can recognize <, not only .lt.'
endif
!调用其它fortran程序文件，见开头的include语句
print*,'include ok! 汉字 ok!^.^'
print*,'gamma:',x,' =',exp(gammln(x))
!比较各精度在运算时的转换，6代表屏幕输出，同"*"(同*同d ^^)
write(6,61),1+1.11111111111d0,1.+1.1111111111111D0,2.1111111111D0
print*,'log(e)=',log(2.7)
!求最大值返回双精度要用dmax1()
print*,(dmax1(5.12345678987654321D0,4.D0)),
& max(5.12345678987654321,4.)
!看看是不是把-当成减号了--不是
print*,1.1d-15
61 format(D30.15)
end

!被一个东东害苦了, 写在下面, 以示反省:
! (nu_m/nu_c)**-0.5D0*(nu_a/nu_m)**(-p_hat/2.D0) 错的
! (nu_m/nu_c)**(-0.5D0)*(nu_a/nu_m)**(-p_hat/2.D0) 对的
!可见括号是多么的重要, 以后有事没事就打个括号, 倍儿有面子&.&
!
!小心:
!print*, 4/5 --> 0
!print*, 4/5. --> 0.800000012
!print*, 4./5. --> 0.800000012
!print*, 4.D0/5 --> 0.8

ifort -r8 foo.f -o foo 这样可以将foo.f中单精度的都转换成双精度计算。比如1e56直接编译不了，如果加上r8就可以编译了。而不用改成1d56.

---

分类内容

库(lib, library): PDA: many useful usage.

---

输入输出:
print*,'Please select the style of the data you scrached. ',
& '(1. x-y(default), 2. logx-logy, 3. logx-y, 4. x-logy) ',
& 'You can just press ENTER to accept the default choice:'
read'(I3)',i_type !要实现print句中的功能的关键语句
if(((i_type-2)*(i_type-3)*(i_type-4)).ne.0)i_type=1
print*,'i_type: ', i_type !test
如果是"read*,i_type", 则运行的时候程序会一直等着, 直到输入了一个值之后才会往下执行, 而如果加了格式控制之后, 就可以直接按enter, 程序继续往下执行. 而i_type被赋值为1.


赋值:

1. integer i/0/

2. integer i

i=0

---

common语句:;
一个common列表中最好只有一种数据类型, 特别是字符型不能和别的类型混在一个公用数据块中.
一个变量不能出现在两个不同的公用块中.

---

一些应用:

用龙格-库塔 (Rounge-Kutta)法解常微分方程, 求解非线性方程的解, 求积分 (Romberg Integration)

---

• 用龙格-库塔(Rounge-Kutta)法解常微分方 程

!参考:Numerical Recipes中的example: xodeinit.for
!用Numerical Recipies上的程序, 采用自适应的变步长方法
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "odeint.for"
include "rkck.for"
include "rkqs.for"
program test_Runge_Kutta
implicit none
integer i,kmax,kount,nbad,nok,nvar,NMAX,KMAXX
parameter (nvar=1,KMAXX=200,NMAX=50)
real dxsav,eps,h1,hmin,x,x1,x2,
& y,ystart(nvar)
COMMON /path/ kmax,kount,dxsav,x(KMAXX),y(NMAX,KMAXX)
! common中的数组还有声明的作用, 所以不能再在前面的声明中声明是数组, 否则就重复了, 会报错.
! (NMAX好像是允许Rounge-Kutta法所求方程组的个数, NMAX=50很高了, 注意这些都只是最大的, 实际比这些要小;)
! NMAX: 微分方程组的维数.
! KMAXX 存储中间量的最大数目, 注意这个值还要和子程序中的相一致;
! nvar 实际Rounge-Kutta法中微分方程组的个数.
! 输出量:
! kount 实际存储的中间量的个数, (并不是实际中间计算的步数)
! nbad 试的过程中失败的次数;
! nok 试的过程中失败的次数;
! x() 自变量
! y(i,KMAXX) 结果
external test_eq,rkqs
open(60,file='testEQ.dat')
x1 = 0. !自变量的最小值
x2 = 10. !自变量的最大值
eps = 1.e-10 !最大误差, 用于决定自适应步长
h1 = 1.e-5 !Intial estimated step.
hmin = 1.e-10 !The minimum step. should be in accord with eps.
x(1) = 0. !The initial value.
ystart(1) = 1. !The initial value.
kmax = KMAXX !the max stored number.
dxsav=(x2-x1)/20.0
! The step of stored x, if dxsav .gt. the step with eps,
! otherwise, the step is greaterthan dxsav, regarding the value of eps.
call odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,test_eq,rkqs)
write(*,'(/1x,a,t30,i4)') 'Successful steps:',nok
write(*,'(1x,a,t30,i3)') 'Bad steps:',nbad
write(*,'(1x,a,t30,i3)') 'Stored intermediate values:',kount
do 20 i=1,kount
write(60,*)x(i),y(1,i)
20 enddo
close(60)
end

subroutine test_eq(x,y,dydx)
real x(*),y(*),dydx(*)
dydx(1) = y(1) !解就是y=e^x
end subroutine test_eq
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

---

• 求解非线性方程的解

!方法:对多项式, 人为取一个大于0一个小于0的两个x1,x2,找中间等于0的x
!参考:Numerical Recipes中的example: xzbrent.for

include 'zbrent.for'
external eq_gamma_3
gamma_3 = zbrent(eq_gamma_3,1.D0,gamma_4,1.D-8)
! eq_gamma_3是一个外部函数作哑元, 第二,三个哑元分别是x1,x2,x在[x1,x2]中, 注意x1和x2的选取要是的中间只有一个根, 否则会漏根; 最后一项是error.

FUNCTION eq_gamma_3(gamma_3)
! see NOTICE: eq_gamma_3
implicit none
real*8 eq_gamma_3
real*8 gamma_3,gamma_4,gamma_34,f
common /C_eq_gamma_3/ gamma_4,f
gamma_34 = gamma_3*gamma_4
& -dsqrt((gamma_3*gamma_3-1.D0)*(gamma_4*gamma_4-1.D0))
eq_gamma_3 = (gamma_34-1.D0)*(4.D0*gamma_34+3.D0)*f
& -(gamma_3-1.D0)*(4.D0*gamma_3+3.D0)
END FUNCTION eq_gamma_3

---

• 求积分

  ---

  1. Romberg Integration

!参考: Numerical Recipes中的example: xqrobm.for

include 'qromb.for'
external fun_D_l
call qromb(fun_D_l,0.D0,z,D)
! fun_D_l是被积函数, 0.D0和z是积分的上下限, D是积分结果
! 误差可以到qromb.for中修改EPS
!
FUNCTION fun_D_l(z)
implicit none
real*8 fun_D_l,omega_l,omega_m,z
! the integral part of the luminosity distance (only for plane universe)
omega_l = 0.7D0
omega_m = 0.3D0
fun_D_l = 1.D0/dsqrt((1.D0+z)*(1.D0+z)*(1.D0+z*omega_m)
& - z*(2.D0+z)*omega_l)
end FUNCTION fun_D_l

但是qromb有局限, 如对跨度太大的积不出来, 例:
fun(x)=1/x**2
call qromb(fun,1.D0,1.D9,result)
结果是:
PAUSE too many steps in qromb statement executed
如果将qromb中的JMAX从20改到30, 在耗时约1分钟之后的结果为23.3017521, 而准确结果是20.73.

---

2. 用龙格-库塔(Rounge-Kutta)法求积分

方法与求微分方程相同: dy/dx=f(x), 只不过这里是f(x), 即x的函数, 所要的结果是y(1,KMAXX)
注意: f(x)这个函数要连续, 否则计算不出来, 会出现步长下溢.
! to do a integration, output is y(1,KMAXX)
! dydx(1) is the function should be integrated.
include "odeint.for"
include "rkck.for"
include "rkqs.for"
PROGRAM test1
C driver for routine odeint
INTEGER KMAXX,NMAX,NVAR
PARAMETER (KMAXX=200,NMAX=50,NVAR=1)
INTEGER i,kmax,kount,nbad,nok,nrhs
REAL dxsav,eps,h1,hmin,x1,x2,x,y,ystart(NVAR)
COMMON /path/ kmax,kount,dxsav,x(KMAXX),y(NMAX,KMAXX)
COMMON nrhs
EXTERNAL derivs,rkqs
nrhs=0
x1=1.0
x2=10.0
ystart(1)=0.0
eps=1.0e-4
h1=.1
hmin=0.0
kmax=100
dxsav=(x2-x1)/20.0
call odeint(ystart,NVAR,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs)
write(*,'(/1x,a,t30,i3)') 'Successful steps:',nok
write(*,'(1x,a,t30,i3)') 'Bad steps:',nbad
write(*,'(1x,a,t30,i3)') 'Function evaluations:',nrhs
write(*,'(1x,a,t30,i3)') 'Stored intermediate values:',kount
write(*,'(/1x,t9,a,t20,a,t33,a)') 'X','Integral'
do 11 i=1,kount
write(*,'(1x,f10.4,2x,2f14.6)') x(i),y(1,i)
11 continue
END

SUBROUTINE derivs(x,y,dydx)
INTEGER nrhs
REAL x,y(*),dydx(*)
real c,h,k,T,nu
COMMON nrhs
c = 2.99792458D10
h = 6.626176D-27
k = 1.380662D-16
T = 1.D4
nrhs=nrhs+1
nu = x
dydx(1) = nu
! dydx(1)= 2.D0*h*nu**3.D0/c/c/(exp(h*nu/k/T)-1.D0)
return
END
C (C) Copr. 1986-92 Numerical Recipes Software ,4-#.

双精度程序
function test_integral(nu)
! test the integration
implicit none
INTEGER KMAXX,NMAX,NVAR
PARAMETER (KMAXX=200,NMAX=50,NVAR=1)
INTEGER i,kmax,kount,nbad,nok,nrhs
REAL*8 dxsav,eps,h1,hmin,x1,x2,x,y,ystart(NVAR),test_integral,nu
COMMON /path/ kmax,kount,dxsav,x(KMAXX),y(NMAX,KMAXX)
COMMON nrhs
EXTERNAL rkqs,ftest
nrhs=0
x1=1.D-6
x2=1.D1
ystart(1)=0.D0
eps=1.0D-4
h1=.1D0
hmin=0.D0
kmax=100
dxsav=(x2-x1)/20.D0
call odeint(ystart,NVAR,x1,x2,eps,h1,hmin,nok,nbad,ftest,rkqs)
write(*,'(/1x,a,t30,i4)') 'Successful steps:',nok
write(*,'(1x,a,t30,i4)') 'Bad steps:',nbad
write(*,'(1x,a,t30,i4)') 'Function evaluations:',nrhs
write(*,'(1x,a,t30,i4)') 'Stored intermediate values:',kount
write(*,'(/1x,t9,a,t20,a,t33,a)') 'X','Integral'
do 11 i=1,kount
write(*,'(1x,D10.4,2x,2D14.6)') x(i),y(1,i)
11 continue
test_integral = y(1,kount)
END function test_integral

subroutine ftest(x,y,dydx)
! test function of integrated for the 'test_integral'
implicit none
integer nrhs
real*8 x,y(*),dydx(*)
common nrhs
nrhs = nrhs+1
dydx(1) = 1.D0
end subroutine
对dydx(1)=1/x**2的测试结果:
eps=1.D-4----- 227320662.
eps=1.D-8----- 1002082.05
eps=1.D-10----- 1000003.2
eps=1.D-14----- 999999.0000240024
eps=1.D-20----- PAUSE too many steps in odeint statement executed
结论: eps不能太大, 太大对积分绝对值也很大的情况, 误差太...大; 也不能太小, 太小会导致不能往下算了. 虽然用起来比Romberg方法麻烦, 但是计算的时候准确度和速度都要好一些.

---

How to print colored output?
怎样做一个待定长度的数组? 答: 在子过程中 如 real A(*)
怎样决定一个文件的长度? 答: end=