[Fortran] 纯文本查看 复制代码
module phase
use imslf90
implicit none
real*8,parameter,private:: pi=3.141592653589793238462643d0
integer,parameter,private:: N=233, NN=N*2, M=1, tran=2000, dot=2000, npara=101,D=0.001d0
!integer,parameter,private:: mxparm=50, N=30, NN=N*2, M=1, tran=200, dot=200, npara=101 !Fig. 4
real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=144.d0/233.d0, beta=89.d0/144.d0, L=144.d0, Kappa=3.d0*pi
!real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=1.d0, beta=24.d0/30.d0, L=30.d0, Kappa=3.d0*pi !Fig. 4
real*8, private:: changeF
private derivs
public plot
contains
function ran() !returns random number between 0 - 1
implicit none
integer , save :: flag = 0
double precision :: ran
if(flag==0) then
call random_seed()
flag = 1
endif
call random_number(ran) ! built in fortran 90 random number function
end function ran
function normal(mean,sigma)
implicit none
integer :: flag
double precision, parameter :: pi = 3.141592653589793239
double precision :: uuu1, uuu2, yyy1, yyy2, normal, mean, sigma
save flag
data flag /0/
uuu1 = ran(); uuu2 = ran()
if (flag.eq.0) then
yyy1 = sqrt(-2.0d0*log(uuu1))*cos(2.0d0*pi*uuu2)
normal = mean + sigma*yyy1
flag = 1
else
yyy2 = sqrt(-2.0d0*log(uuu1))*sin(2.0d0*pi*uuu2)
normal = mean + sigma*yyy2
flag = 0
endif
end function normal
subroutine plot()
integer i,j,k
real*8 y(NN),dydt(NN)
real*8 hh,th,h,t
real*8 dym(NN),dyt(NN),yt(NN),c(N)
real*8 y0(NN)
real*8 parameters(npara),averageVelocity(npara),parameters2(npara),averageVelocity2(npara),temp
character*38 filenm
! set initial conditions
t = 0.0d0
!call random_seed()
do i=1,N
!call random_number(y(i))
y0(i) =(i-1)*b
y0(i+N)=0.d0
enddo
!绝热增加F
do k=1,npara
y=y0
changeF= 0.01d0*(k-1)
t=0.d0; h=0.001d0
hh=h/2.d0
temp=0.d0
do i=1,dot
call derivs(t,y,dydt)
yt=y+h*dydt
t=t+h
call derivs(t,yt,dyt)
y=y+hh*(dydt+dyt)
do j=1,N
temp=temp+y(j+N)
enddo
enddo
temp=temp/N/dot
print*, k,changeF,temp
averageVelocity(k)=temp
parameters(k)=changeF
y0=y
enddo
!输出结果
filenm="eeee.txt"
open (5555, file = filenm)
write (5555,*) " x= [ "
do i = 1 , npara
write(5555,'(f13.5)') parameters(i)
enddo
write (5555,*) " ]; "
write (5555,*) " y= [ "
do i = 1 , npara
write(5555,'(f13.5)') averageVelocity(i)
enddo
write (5555,*) " ]; "
write (5555,*)
write (5555,*) "figure(8)"
write (5555,*) "plot(x,y,'-b^');"
print*,filenm
close (5555)
endsubroutine plot
subroutine derivs(t,y,dydt)
implicit none
real*8 t
real*8 y(NN),dydt(NN),c(N)
integer i,j
real*8 u(N),v(N),w(-N+1:NN),PForce,f1,f2
do i=1,N
u(i)=y(i); v(i)=y(i+N)
c(i)=normal(5.0d0,2.0d0)
!为了处理周期性边界条件Eq.(3)，引入如下中间变量
w(i)=y(i); w(i+N)=y(i)+N*b; w(i-N)=y(i)-N*b
enddo
do i=1,N
PForce=0.d0;f1=0.d0;f2=0.d0
do j=1,M
f1=f1-Kappa*(1.d0-dexp(b+w(i)-w(i+j))) * dexp(b+w(i)-w(i+j)) ! when x_i<x_j
f2=f2+Kappa*(1.d0-dexp(b-w(i)+w(i-j))) * dexp(b-w(i)+w(i-j))
PForce=f1+f2
enddo
dydt(i) =v(i)+D*c(i)
dydt(i+N)=-gamma*v(i) - dsin(2.d0*pi/a*u(i))/2.d0 - dsin(2.d0*pi*beta/a*u(i))/2.d0 -PForce+changeF+D*c(i)
enddo
end subroutine derivs
endmodule phase
program phase_plot
use phase
implicit none
real*8 t1,t2
call cpu_time(t1)
call plot()
call cpu_time(t2)
print*, "Time used:", t2-t1, "seconds."
endprogram phase_plot