1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
| program dynamic
implicit none
INTEGER, parameter :: Nstep=200,SIZE=1000,TSIZE=200,cst_loop=100
REAL*8, parameter :: W0=30.D0,PsmoinsL=1.D-07,L=500.D0,Lt=100.D0
INTEGER*4 :: i,j,NstepT
REAL*8 :: Rco,Aco,Rfo,Acl,gamma_p,delta_lambda_s
REAL*8 :: lambda_s,sigma_ap,sigma_ep,sigma_es,sigma_as
REAL*8 :: gamma_s,R1,R2,tau_s,h,alpha_s
REAL*8 :: c,vp,vs,A21,pi,N,nm,alpha_p,deltaX
REAL*8 :: xmin,xmax,Tmin,Tmax,deltaT,v
REAL*8 :: resultat,Tol
LOGICAL :: fini
REAL*8, DIMENSION(SIZE) :: N1,N2,X,Pp
REAL*8, DIMENSION(SIZE) :: Psplus,Psmoins
REAL*8, DIMENSION(SIZE,TSIZE) :: R12T,R21T,W12T,W21T
REAL*8, DIMENSION(TSIZE) :: T,R2_R,Pout,R2T
REAL*8, DIMENSION(SIZE,TSIZE) :: PsplusT,PsmoinsT,PpT
REAL*8, DIMENSION(SIZE,TSIZE) :: N1T,N2T
!!!!!!!!!!!!!!!!
OPEN(20,FILE="c:\these_driss\résultats\cloop50\pump.dat",ACTION="read")
OPEN(30,FILE="c:\these_driss\résultats\cloop50\signal+.dat",ACTION="read")
OPEN(40,FILE="c:\these_driss\résultats\cloop50\signal-.dat",ACTION="read")
OPEN(50,FILE="c:\these_driss\résultats\cloop50\population_1.dat",ACTION="read")
OPEN(60,FILE="c:\these_driss\résultats\cloop50\population_2.dat",ACTION="read")
!!!!!!!!! initialisation des tableaux
do j=1,TSIZE
T(j)=0.D0
Pout(j)=0.D0
do i=1,SIZE
PsplusT(i,j)=0.D0
PsmoinsT(i,j)=0.D0
N1T(i,j)=0.D0
N2T(i,j)=0.D0
enddo
enddo
!!!!!!!!! données
DATA sigma_ap,sigma_ep,sigma_es,sigma_as,alpha_p,alpha_s,delta_lambda_s,gamma_s,R1,R2,tau_s,h,vp,nm,c/3.D-25,5.D-26,&
2.5D-25,1.4D-27,3.D-03,5.D-03,2.D-08,0.75D0,0.04D0,0.4D0,3.D-07,6.62D-34,319.D12,1.45D0,3.D08/
Tmin=0.D0
Tmax=300.D-09
deltaT=deltaX/v
NstepT=int((Tmax-Tmin)/deltaT)
write(*,*)'NstepT=',NstepT
!!!!!!!!! conditions initiales
do i=1,Nstep+1
N1T(i,1)=(N1(i)/N)*N
N2T(i,1)=(N2(i)/N)*N
PpT(i,1)=Pp(i)
PsplusT(i,1)=Psplus(i)
do j=1,NstepT+1
PsmoinsT(i,j)=Psmoins(i)
end do
end do
!!!!!!!!! Détérmination de R2(T)
do j=Tmin,Tmax,deltaT
R2_R=R2/300.D-09*j
R2_R=REAL((ifft(fft(R2_R)))) !???????????
end do
do j=0,NstepT
T(j+1)=j*deltaT
IF (j.lt.INT(300.D-09/deltaT)) then
R2T=R2_R(j+1)
else if (j.lt.INT(50.D-06/deltaT)) then
R2T=R2
else
R2T=0.D0
end if
end do
!!!!!!!!! boucles de calcul
do WHILE (j.lt.NstepT+1)
R12T=(sigma_ap*PpT(i,j)*gamma_p)/(Aco*h*vp)
R21T=(sigma_ep*PpT(i,j)*gamma_p)/(Aco*h*vp)
W12T=(sigma_as*(PsplusT(i,j)+PsmoinsT(i,j))*gamma_s)/(Aco*h*vs)
W21T=(sigma_es*(PsplusT(i,j)+PsmoinsT(i,j))*gamma_s)/(Aco*h*vs)
N2T(i,j+1)=N2T(i,j)+((R12T+W12T)*N1T(i,j)-(R21T+W21T+A21)*N2T(i,j))*deltaT
N1T(i,j+1)=N-N2T(i,j+1)
do i=1,Nstep-1
PpT(i+1,j+1)=PpT(i,j)+gamma_p*(sigma_ep*N2T(i,j+1)-sigma_ap*N1T(i,j+1))*PpT(i,j)*deltaX-alpha_p*PpT(i,j)*deltaX
PsplusT(i+1,j+1)=PsplusT(i,j)+gamma_s*(sigma_es*N2T(i,j+1)-sigma_as*N1T(i,j+1))*PsplusT(i,j)*deltaX+gamma_s*N2T(i,j+1)*&
sigma_es*2.D0*nm*h*c*c*(delta_lambda_S/lambda_S**3)*deltaX-alpha_s*PsplusT(i,j)*deltaX
R12T=(sigma_ap*PpT(i,j+1)*gamma_p)/(Aco*h*vp)
R21T=(sigma_ep*PpT(i,j+1)*gamma_p)/(Aco*h*vp)
W12T=(sigma_as*(PsplusT(i,j+1)+PsmoinsT(i,j+1))*gamma_s)/(Aco*h*vs)
W21T=(sigma_es*(PsplusT(i,j+1)+PsmoinsT(i,j+1))*gamma_s)/(Aco*h*vs)
N2T(i+1,j+1)=N2T(i,j)+((R12T+W12T)*N1T(i,j)-(R21T+W21T+A21)*N2T(i,j+1))*deltaT
N1T(i+1,j+1)=N-N2T(i+1,j+1)
end do
PsmoinsT(i,j+1)=R2T(j+1)*PsplusT(i,j+1)
do i=Nstep,2,-1
PsmoinsT(i-1,j+1)=PsmoinsT(i,j)+gamma_s*(sigma_es*N2T(i,j+1)-sigma_as*N1T(i,j+1))*PsmoinsT(i,j)*deltaX+gamma_s*N2T(i,j+1)*&
sigma_es*2.D0*nm*h*c*c*(delta_lambda_S/lambda_S**3)*deltaX-alpha_s*PsmoinsT(i,j)*deltaX
R12T=(sigma_ap*PpT(i-1,j)*gamma_p)/(Aco*h*vp)
R21T=(sigma_ep*PpT(i-1,j)*gamma_p)/(Aco*h*vp)
W12T=(sigma_as*(PsplusT(i-1,j)+PsmoinsT(i-1,j))*gamma_s)/(Aco*h*vs)
W21T=(sigma_es*(PsplusT(i-1,j)+PsmoinsT(i-1,j))*gamma_s)/(Aco*h*vs)
N2T(i-1,j+1)=N2T(i-1,j)+((R12T+W12T)*N1T(i-1,j)-(R21T+W21T+A21)*N2T(i-1,j))*deltaT
N1T(i-1,j+1)=N-N2T(i-1,j+1)
end do
Pout(j+1)=(1-R1)*PsmoinsT(1,j+1)
PsplusT(1,j+2)=R1*PsmoinsT(1,j+1)
i=1
j=j+1
end do
END program |
Partager