COSMOS v7.655  COSMOSv7655
(AirShowerMC)
crigCut.f
Go to the documentation of this file.
1 ! ***********************************************************
2 ! To know the arrival prob. as a function of direction at a given
3 ! place.
4 !
5  subroutine crigcut(azmin, zen, rig, prob)
6  implicit none
7 #include "ZrigCut.h"
8 ! Before calling this subroutine, you must have called
9 ! crigCut0 to get the table which contains the information
10 ! of cutoff rigidities.
11 
12 !
13  real*8 azmin ! input. value of azimuth angle at a given location.
14  ! deg.
15  ! For format# 5, this is longitude of a place
16  real*8 zen ! input. value of zenith angle at a given location.
17  ! It may be deg, or cos, depending on the table.
18  ! For format# 5, this is latitude of the place.
19  ! It may be deg, or cos, depending on the table.
20  real*8 rig ! input. value of particle rigidity (GV) p/Z
21 
22  real*8 prob ! output. probability that 'r' can enter the earth.
23  ! For format# 5, this is 0 or 1. 1 dose not necessary
24  ! mean that the particle can enter the earth;
25  ! there is possibility so you have to examine.
26  ! if 0, there is no possibility so you may discard it.
27 !
28 ! Note: This system cannot treat multiple locations at a time.
29 !
30 !
31  real*8 rigcut
32 
33  real*8 azm
34 
35  integer i1, i2, i3
36 
37  azm = azmin
38  if(azm .lt. 0.) then
39  azm = azm + 360.
40  else
41  azm = mod(azm, 360.d0)
42  endif
43 
44  if(rdatafmt .eq. 1 ) then
45  call k4ptdi(rigcuttbl, azmsize, zensize, azmsize, minazm,
46  * minzen, dazm, dzen, azm, zen, rigcut)
47  if(rig .lt. rigcut) then
48  prob = 0.
49  else
50  prob = 1.
51  endif
52  elseif(rdatafmt .le. 4) then
53 ! i3 = log10(rig/MinRig)/LogDRig + 1
54 ! i3 <= log10(rig)*20 < i3+1 (i3=0,1..)
55 ! logrig = log10(rig)/LogDRig
56 ! i3 = logrig + 1
57 ! dlogrig=(logrig - (i3-1))*LogDRig
58 
59  i3 = log10(rig/minrig)/logdrig + 1
60  if(i3 .lt. 1) then
61  prob = 0.
62  goto 100
63  endif
64  if(i3 .ge. rigsize) then
65  prob = 1.
66  goto 100
67  endif
68  if(rdatafmt .eq. 4) then
69  i2=1
70  else
71  i2 = azm/dazm + 1
72  i2 =min(azmsize, i2)
73  endif
74  i1 = (zen-zenmax)/dzen +1 ! DZen < 0
75  i1 = min(zensize, i1)
76  call cgetrigprob(
77  * rigtbl2, zensize, azmsize, rigsize, i1, i2, i3,
78  * prob )
79  prob = min(prob, 1.d0)
80  elseif( rdatafmt .eq. 5) then
81 ! essentiall the same as 1 case. but rather should not
82 ! make interpolation
83  call k2dtblv(rigcuttbl, azmsize, zensize, azmsize, minazm,
84  * minzen, dazm, dzen, azm, zen, rigcut)
85  if(rig .lt. rigcut) then
86  prob = 0.
87  else
88  prob = 1.
89  endif
90  else
91  call cerrormsg('format # error in rigidity table',0)
92  endif
93  100 continue
94  end
95  subroutine k2dtblv(tbl, xs, ys, adj, xm,
96  * ym, dx, dy, x, y, ans)
97  integer xs
98  integer ys
99  integer adj
100  real*8 xm
101  real*8 ym
102  real*8 dx
103  real*8 dy
104  real*8 x, y
105  real*8 ans
106  real*8 tbl(adj, ys)
107 
108  integer i, j
109 
110  i = (x-xm)/dx + 1
111  j = (y-ym)/dy + 1
112  ans = tbl(i,j)
113  end
114 
115 
116  subroutine cgetrigcut(azmin, zen, rigcut)
117  implicit none
118 #include "ZrigCut.h"
119 ! get rigidity cut (for old table only); use crigCut now
120 !
121  real*8 azmin ! input. value of azimuth angle at a given location.
122  ! deg.
123  real*8 zen ! input. value of zenith angle at a given location.
124  real*8 rigcut ! output. cutoff rigidity. in GV
125 !
126 
127  real*8 azm
128 
129  azm = azmin
130 
131  if(azm .lt. 0.) then
132  azm = azm + 360.
133  endif
134  azm = mod(azmin, 360.d0)
135  if(rdatafmt .eq. 1) then
136  call k4ptdi(rigcuttbl, azmsize, zensize, azmsize, minazm,
137  * minzen, dazm, dzen, azm, zen, rigcut)
138  else
139  call cerrormsg(
140  * 'only old cutoff table can be used', 0)
141  endif
142  end
143  subroutine cgetrigprob(tbl,
144  * izen, iphi, irig, i1, i2, i3, prob)
145  implicit none
146 #include "ZrigCut.h"
147 
148  integer izen, iphi, irig
149  real tbl(izen, iphi, irig)
150  integer i1, i2, i3
151  real*8 prob
152 
153 ! if( i3 .lt. RigSize ) then
154 ! prob = dlogrig/LogDRig *
155 ! * (tbl(i1,i2,i3+1)-tbl(i1, i2, i3))
156 ! * + tbl(i1,i2,i3)
157 ! else
158 ! prob = tbl(i1,i2,i3)
159 ! endif
160 
161  prob = tbl(i1,i2,i3)
162 
163  end
164 !------------------------- init for crigCut
165  subroutine crigcut0(file)
166  implicit none
167 #include "Zmanagerp.h"
168 #include "ZrigCut.h"
169 
170  character*(*) file
171  character*192 msg
172  integer icon
173 
174  call copenf(tempdev, file, icon)
175  if(icon .ne. 0) then
176  write(msg, *) ' file specification error '
177  call cerrormsg(msg, 0)
178  endif
179 !
180  read(tempdev, '(1x,i1)') rdatafmt ! format 1 or 2 or 3 or 4
181  if(rdatafmt .eq. 1) then
182  call crigcutfmt1
183  elseif(rdatafmt .eq. 2 .or. rdatafmt .eq. 3) then
184  call crigcutfmt2
185  elseif(rdatafmt .eq. 4) then
186  call crigcutfmt4
187  elseif(rdatafmt .eq. 5) then
188  call crigcutfmt5
189  else
190  write(msg, *) 'rigidity cut data format =',rdatafmt,
191  * ' ivalid'
192  call cerrormsg(msg, 0)
193  endif
194  end
195 ! **********************
196  subroutine crigcutfmt1
197  implicit none
198 #include "Zmanagerp.h"
199 #include "ZrigCut.h"
200 
201 
202  character*192 msg
203 
204  integer icon
205 
206  read(tempdev, *) ! skip
207  read(tempdev, *) place, latit, longi, magdec, azmvalue,
208  * dazm, azmsize, zenvalue, dzen, zensize
209 
210 ! skip until -------------- line
211  call cskipcomment(tempdev, icon)
212  if(icon .ne. 0) then
213  call cerrormsg('cut-off table data has no ---- line', 0)
214  endif
215 ! see if azimuthal angle range is from 0 to 360 or from 0 to 360-DAzm
216 
217  if(azmvalue .eq. 'deg') then
218 ! assume min of Azimuthal angle is 0 deg.
219  minazm = 0.
220  if( (azmsize - 1)* dazm .lt. (360.d0- dazm* 0.1d0) ) then
221 ! fill the last col. by the first col.
222 
223  call creadrigcut(tempdev, rigcuttbl, azmsize, zensize,
224  * azmsize+1)
225  azmsize = azmsize + 1
226  call cfillrigcut(rigcuttbl, azmsize, zensize)
227  else
228  call creadrigcut(tempdev, rigcuttbl, azmsize, zensize,
229  * azmsize)
230  endif
231 
232  else
233 
234  write(msg, *)
235  * ' Azimuthal angle unit must be deg for rigidity cut table'
236  call cerrormsg(msg, 1)
237  write(msg, *) ' But it is ', azmvalue
238  call cerrormsg(msg, 0)
239  endif
240 
241  close(tempdev)
242 
243  if(zenvalue .eq. 'cos') then
244 ! DZen should be negative and min is 1.0
245  minzen = 1.0
246  if(dzen .ge. 0.) then
247  write(msg, *)
248  * ' step of Zenith angle for rigidity cut should be < 0'
249  call cerrormsg(msg, 1)
250  write(msg, *) ' because you give it in cos value'
251  call cerrormsg(msg, 0)
252  endif
253  endif
254 !
255 !
256  write(msg,*) 'Rigidity cut-off table has been read:',
257  * ' place=',place,' latitute=',latit, ' longitude=',longi,
258  * ' mag. dec=', magdec
259  call cerrormsg(msg, 1)
260  end
261 ! ******************
262  subroutine creadrigcut(io, tbl, azm, zen, adj)
263  implicit none
264  integer io, azm, zen, adj
265  real*8 tbl(adj, zen)
266 
267 
268  integer i, j, ios
269  character*150 msg
270 
271  do j=1, zen
272  read(io, *, iostat=ios) ( tbl(i, j), i = 1, azm)
273  if(ios .ne. 0) then
274  write(msg, *) 'Unexpected EOF at rigicity table reading'
275  call cerrormsg(msg, 1)
276  write(msg, *) ' line number=', j, ' azm=',azm, ' zen=',
277  * zen, ' adj=',adj
278  call cerrormsg(msg, 0)
279  endif
280  enddo
281  end
282 ! *******************
283  subroutine cfillrigcut( tbl, azm, zen)
284  implicit none
285  integer azm, zen
286  real*8 tbl(azm, zen)
287  integer i
288 
289  do i = 1, zen
290  tbl(azm, i) = tbl(1, i)
291  enddo
292  end
293 ! **********************************
294  subroutine cqrigcutplace(tlt, tlg, mdec)
295  implicit none
296 #include "ZrigCut.h"
297  real*8 tlt, tlg, mdec
298 
299 
300 ! to inform lat, long, magdec
301  tlt = latit
302  tlg = longi
303  mdec = magdec
304  end
305 
306 ! **********************
307  subroutine crigcutfmt2
308  implicit none
309 #include "Zmanagerp.h"
310 #include "ZrigCut.h"
311 
312  character*192 msg
313 
314  integer icon
315 
316 
317  read(tempdev, *) ! skip
318  read(tempdev, *)
319  * place, latit, longi, magdec, zenvalue, zenmax, dzen, zensize,
320  * azmvalue, minazm, dazm, azmsize, minrig, logdrig, rigsize
321 
322 
323 ! skip until -------------- line
324  call cskipcomment(tempdev, icon)
325  if(icon .ne. 0) then
326  call cerrormsg('cut-off table data has no ---- line', 0)
327  endif
328  call crdrigcut2(
329  * rdatafmt, tempdev, rigtbl2, zensize, azmsize, rigsize)
330  close(tempdev)
331 
332  if(zenvalue .eq. 'cos') then
333 ! DZen should be negative
334  if(dzen .ge. 0.) then
335  write(msg, *)
336  * ' step of Zenith angle for rigidity cut should be < 0'
337  call cerrormsg(msg, 1)
338  write(msg, *) ' because you give it in cos value'
339  call cerrormsg(msg, 0)
340  endif
341  endif
342 !
343 !
344  write(msg,*) 'New rigidity cut-off table has been read:',
345  * ' place=',place,' latitute=',latit, ' longitude=',longi,
346  * ' mag. dec=', magdec
347  call cerrormsg(msg, 1)
348  end
349 ! **********************
350  subroutine crigcutfmt4
351  implicit none
352 #include "Zmanagerp.h"
353 #include "ZrigCut.h"
354 
355  character*192 msg
356 
357  integer icon
358 
359 ! real pw
360 
361  read(tempdev, *) ! skip
362  read(tempdev, *)
363  * place, latit, longi, magdec, zenvalue, zenmax, dzen, zensize,
364  * minrig, logdrig, rigsize
365 
366  azmsize=1
367 
368 
369 ! skip until -------------- line
370  call cskipcomment(tempdev, icon)
371  if(icon .ne. 0) then
372  call cerrormsg('cut-off table data has no ---- line', 0)
373  endif
374  call crdrigcut4(
375  * rdatafmt, tempdev, minrig, rigtbl2, zensize, rigsize)
376  close(tempdev)
377 
378  if(zenvalue .eq. 'cos') then
379 ! DZen should be negative
380  if(dzen .ge. 0.) then
381  write(msg, *)
382  * ' step of Zenith angle for rigidity cut should be < 0'
383  call cerrormsg(msg, 1)
384  write(msg, *) ' because you give it in cos value'
385  call cerrormsg(msg, 0)
386  endif
387  endif
388 !
389  write(msg,*) 'New rigidity cut-off table (fmt4) has been read:',
390  * ' place=',place,' latitute=',latit, ' longitude=',longi,
391  * ' mag. dec=', magdec
392  call cerrormsg(msg, 1)
393  end
394 ! ------------------------
395  subroutine crdrigcut2(fmt, io, tbl, izen, iphi, irig)
396 ! read cut off table for format 2
397  implicit none
398  integer fmt ! input. fmt=2 or 3, if 3, data for izen, iphi, irig are
399  ! not given
400  integer io ! input. table dev. number
401  integer izen ! input. table for zenith angle
402  integer iphi ! input. table for azimuthal angle
403  integer irig ! input. table for rigidity
404  real*4 tbl(izen, iphi, irig) ! output. 3D table which shows prob. that
405  ! given ptcl with izen, iphi, irig can enter the earth
406  ! NOTE: this is real*4--------------------
407 
408  integer i1, i2, i3
409  integer j1, j2, j3
410  character*128 msg
411 
412 
413  do i1= 1 , izen
414  do i2 = 1, iphi
415  do i3 = 1, irig
416  if(fmt .eq. 2) then
417  read(io, *, end=100) j1, j2, j3, tbl(i1, i2, i3)
418  if(i1 .ne. j1+1 .or.
419  * i2 .ne. j2+1 .or.
420  * i3 .ne. j3+1) then
421  write(msg, *)
422  * ' data index mismatch in new rigidit cut table',
423  * i1, j1+1, i2, j2+1, i2, j3+1
424  call cerrormsg(msg, 0)
425  endif
426  else
427  read(io, *) tbl(i1, i2, i3)
428  endif
429  enddo
430  enddo
431  enddo
432  return
433  100 continue
434  call cerrormsg(
435  * 'Unexpected EOF in new rigidity cut table',0)
436  end
437 ! ------------------------
438  subroutine crdrigcut4(fmt, io, minval, tbl, izen, irig)
439 ! read cut off table for format 4
440  implicit none
441  integer fmt ! input. fmt=4
442  integer io ! input. table dev. number
443  real*8 minval ! input minimum value of rigidty (for check)
444  integer izen ! input. table for zenith angle
445  integer irig ! input. table for rigidity
446  real*4 tbl(izen, irig) ! output. 2D table which shows prob. that
447  ! given ptcl with izen, irig can enter the earth
448  ! NOTE: this is real*4--------------------
449 
450  integer i1, i3
451  integer j1, idummy
452  real*4 dummy
453  character*128 msg
454 
455 
456  do i1= 1 , izen
457  do i3 = 1, irig
458  read(io, *, end=100)
459  * j1, idummy, idummy, dummy, tbl(i1, i3)
460  if(i3 .eq. 1) then
461  if( abs(dummy - minval)/minval .gt. 1.e-3) then
462  write(msg,*)
463  * 'check min rigidity in headr=', minval,
464  * ' in table=', dummy
465  call cerrormsg(msg, 0)
466  endif
467  endif
468  if(i1 .ne. j1) then
469  write(msg, *)
470  * ' data index mismatch in new rigidit cut table',
471  * i1, j1
472  call cerrormsg(msg, 0)
473  endif
474  enddo
475  enddo
476  return
477  100 continue
478  call cerrormsg(
479  * 'Unexpected EOF in new rigidity cut table',0)
480  end
481 
482 
483 ! **********************
484  subroutine crigcutfmt5
485  implicit none
486 #include "Zmanagerp.h"
487 #include "ZrigCut.h"
488 
489  character*192 msg
490 
491  integer icon
492 
493 
494  read(tempdev, *) ! skip
495  read(tempdev, *)
496  * zenvalue, zenmax, dzen, zensize,
497  * azmvalue, minazm, dazm, azmsize
498 
499 
500 ! skip until -------------- line
501  call cskipcomment(tempdev, icon)
502  if(icon .ne. 0) then
503  call cerrormsg('cut-off table data has no ---- line', 0)
504  endif
505 
506  if(azmvalue .eq. 'deg') then
507 ! assume min of logitudinal angle is 0 deg.
508  minazm = 0.
509  if( (azmsize - 1)* dazm .lt. (360.d0- dazm* 0.1d0) ) then
510 ! fill the last col. by the first col.
511  call creadrigcut(tempdev, rigcuttbl, azmsize, zensize,
512  * azmsize+1)
513  azmsize = azmsize + 1
514  call cfillrigcut(rigcuttbl, azmsize, zensize)
515  else
516  call creadrigcut(tempdev, rigcuttbl, azmsize, zensize,
517  * azmsize)
518  endif
519  close(tempdev)
520  else
521  write(msg, *)
522  * ' Azimuthal angle unit must be deg for rigidity cut table'
523  call cerrormsg(msg, 1)
524  write(msg, *) ' But it is ', azmvalue
525  call cerrormsg(msg, 0)
526  endif
527 
528 
529  if(zenvalue .eq. 'cos') then
530 ! DZen should be negative and min is 1.0
531  minzen = 1.0
532  if(dzen .ge. 0.) then
533  write(msg, *)
534  * ' step of Zenith angle for rigidity cut should be < 0'
535  call cerrormsg(msg, 1)
536  write(msg, *) ' because you give it in cos value'
537  call cerrormsg(msg, 0)
538  endif
539  endif
540 !
541 !
542  write(msg,*) 'Rough rigidity cut-off table has been read:'
543  call cerrormsg(msg, 1)
544  end
545 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
crdrigcut2
subroutine crdrigcut2(fmt, io, tbl, izen, iphi, irig)
Definition: crigCut.f:396
cgetrigcut
subroutine cgetrigcut(azmin, zen, rigcut)
Definition: crigCut.f:117
crigcut
subroutine crigcut(azmin, zen, rig, prob)
Definition: crigCut.f:6
crigcutfmt4
subroutine crigcutfmt4
Definition: crigCut.f:351
crigcut0
subroutine crigcut0(file)
Definition: crigCut.f:166
creadrigcut
subroutine creadrigcut(io, tbl, azm, zen, adj)
Definition: crigCut.f:263
crigcutfmt1
subroutine crigcutfmt1
Definition: crigCut.f:197
crigcutfmt5
subroutine crigcutfmt5
Definition: crigCut.f:485
crigcutfmt2
subroutine crigcutfmt2
Definition: crigCut.f:308
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
copenf
subroutine copenf(io, fnin, icon)
Definition: copenf.f:8
cqrigcutplace
subroutine cqrigcutplace(tlt, tlg, mdec)
Definition: crigCut.f:295
cfillrigcut
subroutine cfillrigcut(tbl, azm, zen)
Definition: crigCut.f:284
crdrigcut4
subroutine crdrigcut4(fmt, io, minval, tbl, izen, irig)
Definition: crigCut.f:439
k4ptdi
subroutine k4ptdi(f, im, jm, iadj, x0, y0, hx, hy, x, y, ans)
Definition: k4ptdi.f:21
cgetrigprob
subroutine cgetrigprob(tbl, izen, iphi, irig, i1, i2, i3, prob)
Definition: crigCut.f:145
cskipcomment
subroutine cskipcomment(io, icon)
Definition: cskipComment.f:19
k2dtblv
subroutine k2dtblv(tbl, xs, ys, adj, xm, ym, dx, dy, x, y, ans)
Definition: crigCut.f:97