-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtdpt.f90
149 lines (128 loc) · 4.79 KB
/
tdpt.f90
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
140
141
142
143
144
145
146
147
148
149
program main
use MKL_DFTI
use quad2d_qgaus_mod
use quad2d_qromb_mod
use Const
use Laser
use Functions
implicit none
real(DOUBLE), dimension(7) :: paras1, paras2, paras3, paras4
real(DOUBLE), dimension(:), allocatable :: t, ele
real(DOUBLE) :: omega, tau, dt
real(DOUBLE) :: E1, E2, Ea, Eb, Ei, Ef
complex(DOUBLE) :: Ka, Kb
complex(DOUBLE), dimension(:), allocatable :: Kas, Kbs
integer(SINGLE), parameter :: nE = 150
real(DOUBLE) :: dE, Eg, sigmas(4)
real(DOUBLE), dimension(nE) :: energy, ss(nE)
complex(DOUBLE), dimension(nE) :: KK
real(DOUBLE), dimension(nE, nE) :: spectrum
real(DOUBLE), dimension(nE*nE) :: Ex, Ey
type(cos2field) :: f1, f2, f3
integer(SINGLE) :: iE1, iE2, i, j
character(10) nowtime
! namelist/He/sigma0, E0, yw, ya, P, y0, y1
call CheckDelete("data/ele.bin")
call CheckDelete("data/je.bin")
call CheckDelete("data/sigma.bin")
paras1 = [1.361d1, 9.492d2, 1.469d0, 3.188d0, 2.039d0, 4.434d-1, 2.136d0]
paras2 = [1.720d0, 1.369d4, 3.288d1, 2.963d0, 0d0, 0d0, 0d0]
paras3 = [1.709d1, 2.16d1, 2.645d2, 4.796d0, 4.185d-1, 1.688d0, 8.942d-1]
paras4 = [2.556d0, 4.140d0, 1.337d1, 1.191d1, 1.570d0, 6.634d0, 1.272d-1]
! open (100, file='input.nml', delim='apostrophe')
! read (unit=100, nml=He)
! close (100)
! omega = 53d0/transE
! ! tau = 3*2*PI/omega
tau = 5d-15/2.4189d-17
f1 = cos2field(ItoE0(1.0d15), 1.86, 0.0, tau, tau/2.0)
f2 = cos2field(ItoE0(1.0d15), 1.92, 0.0, tau, tau/2.0)
f3 = cos2field(ItoE0(0.0d15)/30.0, 0.10, 0.0, 4*tau, 2*tau)
dt = 0.05
nt = int((max(f1%tau/2 + f1%tm, f2%tau/2 + f2%tm, f3%tau/2 + f3%tm)+0*tau)/dt)
allocate (t(nt), ele(nt))
allocate (Kas(nt-1), Kbs(nt-1))
forall (i=1:nt) t(i) = i*dt
do i = 1, nt
ele(i) = GetEle(f1, t(i)) + GetEle(f2, t(i)) + GetEle(f3, t(i))
end do
open (100, file='data/ele.bin', access='stream')
write (100) (t(i), ele(i), i=1, nt)
close (100)
! dE = 27.2d0/nE/transE
dE = 3d0/nE
! energy = [1:nE]*dE
energy = [-nE/2:nE/2 - 1]*dE
Ei = -79.0 / transE
Eg = -54.4 / transE
call time(nowtime)
print'("time is:", T20, A)', nowtime
spectrum = 0d0
do i = 1, nE
ss(i) = sigma(i*dE, paras2)
end do
open (100, file="data/sigma.bin", access="stream")
write (100) (i*dE, ss(i), i=1, nE)
close (100)
do iE2 = 1, nE
do iE1 = 1, iE2
! E1 = energy(iE1)
! E2 = energy(iE2)
E1 = (energy(iE1))**2/2
E2 = (energy(iE2))**2/2
Ea = E1 + Eg
Eb = E2 + Eg
Ef = E1 + E2
call quad2d_qromb(t(1), t(nt), Ka, Ea-Ei, Ef-Ea, f1, f2, f3) ! energy in a. u.
call quad2d_qromb(t(1), t(nt), Kb, Eb-Ei, Ef-Eb, f1, f2, f3)
! Ka = Ka*(1+exp(IM*(E1+E2-Ei)*10.0-2*PI))
! Kb = Kb*(1+exp(IM*(E1+E2-Ei)*10.0-2*PI))
sigmas(1) = sigma(Ea - Ei, paras1)
sigmas(2) = sigma(Ef - Ea, paras2)
sigmas(3) = sigma(Eb - Ei, paras1)
sigmas(4) = sigma(Ef - Eb, paras2)
spectrum(iE1, iE2) = abs(energy(iE1)*energy(iE2))*abs(sqrt(sigmas(1)*sigmas(2)/(Ea - Ei)/(Ef - Ea))*Ka &
+ sqrt(sigmas(3)*sigmas(4)/(Eb - Ei)/(Ef - Eb))*Kb)**2
! spectrum(iE1, iE2) = abs(sqrt(sigmas(1)*sigmas(2)/(Ea - Ei)/(Ef - Ea))*Ka &
! + sqrt(sigmas(3)*sigmas(4)/(Eb - Ei)/(Ef - Eb))*Kb)**2
spectrum(iE2, iE1) = spectrum(iE1, iE2)
Ex(nE*(iE2 - 1) + iE1) = (iE1 - nE/2)*dE
Ey(nE*(iE1 - 1) + iE2) = (iE2 - nE/2)*dE
Ex(nE*(iE1 - 1) + iE2) = (iE2 - nE/2)*dE
Ey(nE*(iE2 - 1) + iE1) = (iE1 - nE/2)*dE
! Ex(nE*(iE2 - 1) + iE1) = iE1*dE
! Ey(nE*(iE1 - 1) + iE2) = iE2*dE
! Ex(nE*(iE1 - 1) + iE2) = iE2*dE
! Ey(nE*(iE2 - 1) + iE1) = iE1*dE
end do
if (mod(iE2, 20) == 0) then
print *, "iE2=", iE2
end if
end do
! Ex = Ex * transE
! Ey = Ey * transE
open (100, file='data/je.bin', access='stream')
write (100) ((Ex(nE*(j - 1) + i), Ey(nE*(i - 1) + j), spectrum(i, j), i=1, nE), j=1, nE)
close (100)
call time(nowtime)
print'("time is:", T20, A)', nowtime
deallocate (t, ele)
call system('cd fig & gnuplot ele.gp & gnuplot je.gp & gnuplot sigma.gp')
end program main
function shapefun(x, y, w1, w2, f1, f2, f3) !!! w1, w2 in atomic unit
use Const
use Laser
implicit none
real(DOUBLE), intent(in) :: x, w1, w2
type(cos2field) :: f1, f2, f3
real(DOUBLE) :: Ex
real(DOUBLE), dimension(:), intent(in) :: y
real(DOUBLE), dimension(size(y)) :: Ey
complex(DOUBLE), dimension(size(y)) :: shapefun
integer(SINGLE) :: i
Ex = GetEle(f1, x) + GetEle(f2, x) + GetEle(f3, x)
do i = 1, size(y)
Ey(i) = GetEle(f1, y(i)) + GetEle(f2, y(i)) + GetEle(f3, y(i))
end do
shapefun = Ex*Ey*exp(IM*(w1*x + w2*y))
end function shapefun