[Fortran] 纯文本查看 复制代码
program sasm !model with P32
implicit real*8(a-h,o-z)
integer,parameter::n=15
integer,parameter::i=10
real*8,parameter::a=1d0
real*8,dimension(n)::p12,p21,p23,p32
real*8,dimension(n)::b,c
real*8,dimension(n)::r1,r2
real*8,dimension(n)::c1,c2,c3
real*8,dimension(n)::m1,m2,m3
real*8,dimension(n)::x1,x2,x3
real*8,dimension(n)::yn1,yn2,yn3,yn
open(1,file='D:/program/0809/frequency12.txt')
open(2,file='D:/program/0809/frequency21.txt')
open(3,file='D:/program/0809/frequency23.txt')
open(5,file='D:/program/0809/frequency32.txt')
open(7,file='D:/program/0809/1200.txt')
do ix=1,n
read(1,*)xtemp,ytemp,p12(ix)
read(2,*)xtemp,ytemp,p21(ix)
read(3,*)xtemp,ytemp,p23(ix)
read(5,*)xtemp,ytemp,p32(ix)
end do
close(1)
close(2)
close(3)
close(5)
b(i)=p12(i)+p21(i)+p23(i)+p32(i)
c(i)=p12(i)*p23(i)+p12(i)*p32(i)+p21(i)*p32(i)
r1(i)=(-1d0*b(i)+dsqrt(b(i)**(2d0)-4d0*a*c(i)))/(2d0*a)
r2(i)=(-1d0*b(i)-dsqrt(b(i)**(2d0)-4d0*a*c(i)))/(2d0*a)
c1(i)=p12(i)*(r2(i)+p12(i)+p21(i))/(r1(i)*(r1(i)-r2(i)))
c2(i)=p12(i)*(r1(i)+p12(i)+p21(i))/(r2(i)*(r2(i)-r1(i)))
c3(i)=1d0-c1(i)-c2(i)
m1(i)=(r1(i)+p12(i))/p21(i)
m2(i)=(r2(i)+p12(i))/p21(i)
m3(i)=p12(i)/p21(i)
x1(i)=(r1(i)**(2d0)+r1(i)*(p12(i)+p21(i)+p23(i))+p12(i)*p23(i))/(p21(i)*p32(i))
x2(i)=(r2(i)**(2d0)+r2(i)*(p12(i)+p21(i)+p23(i))+p12(i)*p23(i))/(p21(i)*p32(i))
x3(i)=(p12(i)*p23(i))/(p21(i)*p32(i))
do t=0d0,10000000d0,10000d0
yn1(i)=c1(i)*dexp(r1(i)*t)+c2(i)*dexp(r2(i)*t)+c3(i)
yn2(i)=c1(i)*m1(i)*dexp(r1(i)*t)+c2(i)*m2(i)*dexp(r2(i)*t)+c3(i)*m3(i)
yn3(i)=c1(i)*x1(i)*dexp(r1(i)*t)+c2(i)*x2(i)*dexp(r2(i)*t)+c3(i)*x3(i)
yn(i)=yn1(i)+yn2(i)+yn3(i)
write(*,*)t,yn(i)
write(7,*)yn1(i),yn2(i),yn3(i),yn(i)
end do
stop
end program sasm
