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
| 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
REAL*8 :: lambda_s,sigma_ap,sigma_ep,sigma_es,sigma_as
REAL*8 :: gamma_s,R1,R2,tau_s,h
REAL*8 :: c,vp,vs,A21,pi,N
REAL*8 :: xmin,xmax,Tmin,Tmax,deltaT
REAL*8 :: resultat,Tol
LOGICAL :: fini
REAL*8, DIMENSION(SIZE) :: N1,N2,X,Pp
REAL*8, DIMENSION(SIZE) :: Psplus,Psmoins
REAL*8, DIMENSION(SIZE) :: R12,R21,W12,W21
REAL*8, DIMENSION(TSIZE) :: T,R2_R,Pout
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(30,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
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)
R12=(sigma_ap*PpT(i,j)*gamma_p)/(Aco*h*vp)
R21=(sigma_ep*PpT(i,j)*gamma_p)/(Aco*h*vp)
W12=(sigma_as*(PsplusT(i,j)+PsmoinsT(i,j))*gamma_s)/(Aco*h*vs)
W21=(sigma_es*(PsplusT(i,j)+PsmoinsT(i,j))*gamma_s)/(Aco*h*vs)
N2T(i,j+1)=N2T(i,j)+((R12+W12)*N1T(i,j)-(R21+W21+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
R12=(sigma_ap*PpT(i+1,j)*gamma_p)/(Aco*h*vp)
R21=(sigma_ep*PpT(i+1,j)*gamma_p)/(Aco*h*vp)
W12=(sigma_as*(PsplusT(i+1,j)+PsmoinsT(i+1,j))*gamma_s)/(Aco*h*vs)
W21=(sigma_es*(PsplusT(i+1,j)+PsmoinsT(i+1,j))*gamma_s)/(Aco*h*vs)
N2T(i+1,j+1)=N2T(i+1,j)+((R12+W12)*N1T(i+1,j)-(R21+W21+A21)*N2T(i+1,j))*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)*&
sigma_es*2.D0*nm*h*c*c*(delta_lambda_S/lambda_S**3)*deltaX-alpha_s*Psmoins(i)*deltaX
end do
end do
END program |
