This repository has been archived by the owner on Apr 28, 2019. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
convolution_model_Fresnel_xyz.f
409 lines (333 loc) · 11.2 KB
/
convolution_model_Fresnel_xyz.f
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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
c Given Q(z) longitudinal development of excess charge
c computes vector potential at observer position (Fresnel)
c in time domain.
c
c Calculates x,y,z components of vector potential
c Accounts for polarization dependence on shower depth
c
c 2D-model i.e. the lateral distribution is accounted for
c by a parameterisation of the vector potential.
c
c S.I. units
c Based on "Practical and accurate calculations of Askaryan radiation" PRD in press (2012).
c J. Alvarez-Muniz, A. Romero-Wolf, E. Zas
c program model
subroutine model
implicit double precision (a-h,o-z)
parameter(pi=3.1415926)
parameter(c=2.99792458e8) ! speed of light m^s-1
parameter(xmu=12.566370e-7) ! magnetic permeability N A^-2
parameter(e=1.602177e-19) ! electron charge C
character*80 filein
common/filename/filein
common/antennae/xa(100),ya(100),za(100),rhoa(100),Ra(100)
common/antenna_index/iant
common/depthmax/zend
common/ice/xn,x0,rho
common/cherenkov/cher
common/const/factor
common/time/tobs_sec
common/charge/xntot
common/shower/E_TeV
common/units/s_to_ns,xns_to_s
c Vegas routine commons
CHARACTER NAMEOU*12
COMMON/BVEG1/NCALL,ITMX,NPRN,NDEV,XL(18),XU(18),ACC
COMMON/NAME/NAMEOU
external vegas,xintegrand_x,xintegrand_y,xintegrand_z
ncall=10000
itmx=30
c Vegasd
nameou='vegas.out'
c Output file containing waveforms of vector potential vs observer's time
open (unit=22,status='unknown',file='output.dat')
write(*,*) 'File containing long shower profile'
read(*,*) filein
write(*,*) filein
write(*,*) 'Shower energy [TeV]'
read(*,*) E_TeV
write(*,*) E_TeV
c ICE
xn=1.78 ! refractive index
x0=36.08 ! radiation length g cm^-2
rho=0.924 ! density g cm^-3
write(*,*) 'Number of antennas'
read(*,*) nant
write(*,*) nant
c Position of shower maximum in g cm^-2 if antenna
c needs to be referred to zmax
zmax=0.
zmax=zmax/rho/100. ! m
write(*,*) 'Position of antenna i'
write(*,*) '(x,y,z) in g cm^-2'
write(*,*) 'Primary injected at (0,0,0)'
do ia=1,nant
read(*,*) xa(ia),ya(ia),za(ia)
write(*,*) xa(ia),ya(ia),za(ia)
c Convert antenna positions from g cm^-2 to meters
xa(ia)=xa(ia)/rho/100. ! m
ya(ia)=ya(ia)/rho/100. ! m
za(ia)=za(ia)/rho/100. ! m
rhoa(ia)=sqrt(xa(ia)*xa(ia)+ya(ia)*ya(ia)) ! m
Ra(ia)=sqrt(rhoa(ia)*rhoa(ia)+
# (za(ia)-zmax)*(za(ia)-zmax)) ! m
c write(55,*) ia-1,Ra(ia)
end do
c First call to interpolation to initialize zend and xntot
dum=xnep(10.)
write(*,*) 'Integral[N(z)dz] ',xntot
c Limits in integration
c Vegasd
c Integration limits in shower depth z'
xl(1)=0.
xu(1)=zend
c Assumes shower travels at speed of light
beta=1.
cher=acos(1./xn)
c Factor in Eq.(22) PRD paper
factor=-xmu/(4.*pi)
s_to_ns=1.e9
xns_to_s=1.e-9
write(*,*) 'Time resolution [ns]'
read(*,*) dtbin
write(*,*) dtbin
write(*,*) 'Time min. and max. [ns]'
read(*,*) tmin,tmax
itmax=int((tmax-tmin)/dtbin)
c Loop in observer's time
c Compute vector potential at each antenna position
c Electric field can be computed by a simple derivative w.r.t. time.
do iant=1,nant
write(*,*)
#'# Time [ns] Vector potential [Vs] antenna position (x,y,z) m',
#xa(iant),ya(iant),za(iant)
write(22,*)
#'# Time [ns] Vector potential [Vs] antenna position (x,y,z) m',
#xa(iant),ya(iant),za(iant)
vp_x=0.
vp_y=0.
vp_z=0.
c Calculate x,y,z components of vector potential
do it=1,itmax
tobs_plot=tmin+dtbin*(it-1)+dtbin/2.
tobs=tobs_plot+(xn*Ra(iant)/c)*s_to_ns
tobs_sec=tobs*xns_to_s
c x-component of Eq.(22) PRD paper
if(xa(iant).eq.0.) then
vp_x=0.
else
call vegas(1,xintegrand_x,vp_x,sal_x,schi_x)
end if
c y-component of Eq.(22) PRD paper
if(ya(iant).eq.0.) then
vp_y=0.
else
call vegas(1,xintegrand_y,vp_y,sal_y,schi_y)
end if
c z-component of Eq.(22) PRD paper
call vegas(1,xintegrand_z,vp_z,sal_z,schi_z)
vp_x=factor*vp_x
vp_y=factor*vp_y
vp_z=factor*vp_z
write(*,*) tobs,vp_x,vp_y,vp_z
c write(*,*) tobs_plot,vp_x,vp_y,vp_z
write(22,*) tobs,vp_x,vp_y,vp_z
c write(22,*) tobs_plot,vp_x,vp_y,vp_z
end do
write(*,*)
c write(22,*)
end do
end
c ----------------------------------------------------------------------------
c x-component of integral in Eq.(22) PRD paper
function xintegrand_x(qww)
implicit double precision (a-h,o-z)
save
common/x/x(1),xd(17)
common/antennae/xa(100),ya(100),za(100),rhoa(100),Ra(100)
common/antenna_index/iant
common/cherenkov/cher
common/time/tobs_sec
common/dist/R
c Depth in m
z=x(1) ! m
c Argument of F_p function. Eq.(22) PRD.
arg=argument(z,tobs_sec)
c Distance from position in shower depth z' to each antenna.
c Denominator in Eq. (22) PRD paper
R = sqrt(rhoa(iant)*rhoa(iant)+
# (za(iant)-z)*(za(iant)-z)) ! m
c x-component of v_perp
c Account for polarization dependence on shower position
c Eq. (22) PRD paper
u_x=xa(iant)/R
u_z=(za(iant)-z)/R
beta_z=1.
vperp_x=u_x*u_z*beta_z
xintegrand_x=-vperp_x*xnep(z)*F_p(arg)/R
return
end
c ----------------------------------------------------------------------------
c y-component of integral in Eq.(22) PRD paper
function xintegrand_y(qww)
implicit double precision (a-h,o-z)
save
common/x/x(1),xd(17)
common/antennae/xa(100),ya(100),za(100),rhoa(100),Ra(100)
common/antenna_index/iant
common/cherenkov/cher
common/time/tobs_sec
common/dist/R
c Depth in m
z=x(1) ! m
c Argument of F_p function. Eq.(22) PRD.
arg=argument(z,tobs_sec)
c Distance from position in shower depth z' to each antenna.
c Denominator in Eq. (22) PRD paper
R = sqrt(rhoa(iant)*rhoa(iant)+
# (za(iant)-z)*(za(iant)-z)) ! m
c y-component of v_perp
c Account for polarization dependence on shower position
c Eq.(22) PRD
u_y=ya(iant)/R
u_z=(za(iant)-z)/R
beta_z=1.
vperp_y=u_y*u_z*beta_z
xintegrand_y=-vperp_y*xnep(z)*F_p(arg)/R
return
end
c ----------------------------------------------------------------------------
c z-component of integral in Eq.(22) PRD paper
function xintegrand_z(qww)
implicit double precision (a-h,o-z)
save
c Common for vegas routine
common/x/x(1),xd(17)
common/antennae/xa(100),ya(100),za(100),rhoa(100),Ra(100)
common/antenna_index/iant
common/cherenkov/cher
common/time/tobs_sec
common/dist/R
c Depth in m
c Specify first and only variable to be integrated by vegas
z=x(1) ! m
c Argument of F_p function. Eq.(22) PRD.
arg=argument(z,tobs_sec)
c Distance from position in shower depth z' to each antenna
c Denominator in Eq. (22) PRD paper
R = sqrt(rhoa(iant)*rhoa(iant)+
# (za(iant)-z)*(za(iant)-z)) ! m
c z-component of v_perp
c Eq.(22) PRD
u_x=xa(iant)/R
u_y=ya(iant)/R
beta_z=1.
vperp_z=-(u_x*u_x+u_y*u_y)*beta_z
xintegrand_z=-vperp_z*xnep(z)*F_p(arg)/R
return
end
c -----------------------------------------------------------------------------
c Argument of F_p Eq.(22) PRD paper
c -----------------------------------------------------------------------------
function argument(z,t)
implicit double precision (a-h,o-z)
save
parameter(c=2.99792458e8) ! speed of light m^s-1
common/antennae/xa(100),ya(100),za(100),rhoa(100),Ra(100)
common/antenna_index/iant
common/cherenkov/cher
common/ice/xn,x0,rho
common/dist/R
c Assumes shower travels at speed of light
beta=1.
argument=z-(beta*c*t-xn*R) ! m
return
end
c -----------------------------------------------------------------------------
c Function F_p Eq.(15) PRD paper.
c -----------------------------------------------------------------------------
function F_p(arg)
implicit double precision (a-h,o-z)
save
parameter(pi=3.1415926)
parameter(c=2.99792458e8) ! speed of light m^s-1
parameter(xmu=12.566370e-7) ! magnetic permeability N A^-2
parameter(e=1.602177e-19) ! electron charge C
common/cherenkov/cher
common/charge/xntot
common/shower/E_TeV
common/units/s_to_ns,xns_to_s
common/ice/xn,x0,rho
c Factor accompanying the F_p function in Eq.(15) in PRD paper
fc=4.*pi/(xmu*sin(cher))
c Converts arg into time
beta=1.
c Note that Acher peaks at tt=0 which corresponds to the observer time.
c The shift from tobs to tt=0 is done when defining argument
tt=(-arg/(c*beta))*s_to_ns ! Parameterisation of A_Cherenkov with t in ns
c Parameterisation of vector potential at Cherenkov angle
c Eq.(16) PRD paper.
Af=-4.5e-14 ! V s
if (tt.ge.0.) then
Acher=Af*E_TeV*(exp(-abs(tt)/0.057)+(1.+2.87*abs(tt))**(-3))
else
Acher=Af*E_TeV*(exp(-abs(tt)/0.03)+(1.+3.05*abs(tt))**(-3.5))
end if
c Cut fit above +/-5 ns
c if (abs(tt).gt.5.) Acher=1.e-30
c Obtain "shape" of Lambda-function from vp at Cherenkov angle
c xntot = LQ_tot in PRD paper
F_p=Acher*fc/xntot
return
end
C-----------------------------------------------------------
c Simple interpolation in depth of longitudinal profile of excess charge.
C-----------------------------------------------------------
function xnep(z)
implicit double precision (a-h,o-z)
save
common/depthmax/zend
dimension xx(500),yy(500),ylog(500)
character*80 filein
logical firstgo
data firstgo / .true. /
common/filename/filein
common/ice/xn,x0,rho
common/charge/xntot
nnint=1
if (firstgo .eqv. .true.) then
xntot=0.
c Read longitudinal profile
npoints=0
open (unit=11,status='old',file=filein)
read(11,*)
read(11,*)
read(11,*)
do j=1,500
read(11,*,end=11) xx(j),dum1,dum2,yy(j),dum3,dum4,dum5
xx(j)=xx(j)*x0/rho/100. ! m
xntot=xntot+yy(j)
npoints=npoints+1
end do
11 close (unit=21)
xntot=xntot*(xx(2)-xx(1)) ! m
npoints=j-1
firstgo=.false.
zend=xx(npoints)
write(*,*) 'Depth max (m) ',zend
do j=1,npoints
if (yy(j).gt.0.) then
ylog(j)=log(yy(j))
else
ylog(j)=0.
end if
end do
end if
xnep=ddivdif(ylog,xx,npoints,z,nnint)
xnep=-exp(xnep) ! Excess negative charge => - sign
RETURN
END
C ******************************************************************
C ******************************************************************
C ******************************************************************
C ******************************************************************