COSMOS v7.655  COSMOSv7655
(AirShowerMC)
csPrimAng.f
Go to the documentation of this file.
1 #if defined NEXT486
2 #define IMAG_P dimag
3 #elif defined PCLinux
4 #define IMAG_P dimag
5 #else
6 #define IMAG_P imag
7 #endif
8 
9 ! Test program is in Test/testSPrimAng.f
10 !----------------------samples primary direction cos. at observation place.
11 ! system is det
12 !
13  subroutine csprimang(dir)
14  implicit none
15 !
16 #include "Zglobalc.h"
17 #include "Zcoord.h"
18 #include "Zincidentp.h"
19 #include "Zobs.h"
20 #include "Zobsp.h"
21 
22  type(coord)::dir
23  character*50 msg
24  real*8 u, cs, sn, sint, fai, za
25 
26  if(za1ry .eq. 'is') then
27  call csprimisoang(dir)
28  elseif(za1ry .eq. 'ps') then
29 ! from point source
30  call csprimpsang(dir)
31  elseif(za1ry .eq. 'aps') then
32 ! around point source
33  call csprimapsang(dir)
34  elseif(za1ry .eq. 'cos' ) then
35 ! cos dcos
36  call ksamplin
37  * (1.0d0, 0.d0, real(CosZenith), IMAG(coszenith), za)
38  dir%r(3) = - za ! going down is negative
39  if(obsplane .eq. spherical) then
40 ! actually setting here is dummy.
41  call kcossn(cs, sn)
42  else
43  call rndc(u)
44  fai = (imag_p(azimuth)- real(azimuth)) *u +
45  * real(azimuth)
46  cs = cos(torad*fai)
47  sn = sin(torad*fai)
48  endif
49  sint = sqrt(1.d0-dir%r(3)**2)
50  dir%r(1) = - sint*cs ! - is needed for going down ptcl.
51  dir%r(2) = - sint*sn
52  dir%sys = 'det' ! fof ObsPlane=spherical, this is reset later
53  else
54  write(msg,*) 'strange Za1ry=',za1ry
55  call cerrormsg(msg, 0)
56  endif
57  end
58 #if defined NEXT486
59 #define IMAG_P dimag
60 #elif defined PCLinux
61 #define IMAG_P dimag
62 #else
63 #define IMAG_P imag
64 #endif
65 ! -----------------------------------------
66 
67  subroutine csprimisoang(dir)
68  implicit none
69 !
70 #include "Zglobalc.h"
71 #include "Zcoord.h"
72 #include "Zincidentp.h"
73  type(coord)::dir
74 !
75  real*8 u, cs, sn, sint, fai
76  call rndc(u)
77  dir%r(3) = -( (imag_p(coszenith)- real(CosZenith) )*u +
78  * real(CosZenith) ) ! going down is negative
79  call rndc(u)
80  fai = (imag_p(azimuth)- real(azimuth)) *u +
81  * real(azimuth)
82  cs = cos(torad*fai)
83  sn = sin(torad*fai)
84 !ccc call kcossn(cs, sn)
85  sint = sqrt(1.d0-dir%r(3)**2)
86  dir%r(1) = - sint*cs ! - is needed for going down ptcl.
87  dir%r(2) = - sint*sn
88  dir%sys = 'det'
89  end
90 ! -----------------------------
91 ! sample 1ry angle from point source
92  subroutine csprimpsang(dir)
93  implicit none
94 
95 #include "Ztrack.h"
96 #include "Zincidentp.h"
97 #include "Zincidentv.h"
98  type(coord)::dir
99 !
100  real*8 u, h, w1, w2, w3
101  call rndc(u)
102  h=(imag_p(obsvhour)-real(obsvhour))*u + real(Obsvhour)
103  call rndc(u)
104  if(u .lt. .5) then
105  h=-h
106  endif
107 ! for source SourceDec at hour h, get horizontal vector
108  call kdhtoh(sourcedec, h, w1, w2, w3)
109  w3 = - w3 ! note vector def. is opposit
110  w2 = - w2
111  w1 = - w1
112 ! convert to detector system
113  call khtad(w1, w2, w3, dir%r(1), dir%r(2), dir%r(3))
114  dir%sys = 'det'
115  end
116 ! ---------------------------------------
117 ! arround point sorucec
118  subroutine csprimapsang(dir)
119 
120  implicit none
121 #include "Zglobalc.h"
122 #include "Ztrack.h"
123 #include "Zincidentp.h"
124 #include "Zincidentv.h"
125  type(coord)::dir
126 
127  real*8 u, w3p, decx, h, w1, w2, w3
128 ! sample angle in horizontal system
129  call rndc(u)
130  w3p=(cspsmx -cspsmn)*u + cspsmn
131  decx=acos(w3p)*todeg
132 !c if(decx .gt. 90.) then ! bug found Oct.15. by Kawata.
133  decx=90.-decx
134 !c endif
135  call rndc(u)
136  h=( imag_p(obsvhour)-real(obsvhour))*u + real(Obsvhour)
137  call rndc(u)
138  if(u .lt. .5) then
139  h=-h
140  endif
141  call kdhtoh(decx, h, w1, w2, w3)
142  w3 = - w3 ! note vector direction is opposit.
143  w2 = - w2
144  w1 = - w1
145  call khtad(w1, w2, w3, dir%r(1), dir%r(2), dir%r(3))
146  dir%sys = 'det'
147  end
148 ! ------------------------- init for sampling 1ry angle
149  subroutine cinisprimang
150 
151  implicit none
152 #include "Zglobalc.h"
153 #include "Ztrack.h"
154 #include "Zincidentp.h"
155 #include "Zincidentv.h"
156 
157  real*8 smx, smn, w3min, hmaxd, hmind, hmax, hmin,
158  * wh1, wh2, wh3, w3max
159  integer icon, jcon
160  external cblkincident
161  character*200 msg
162 !
163 ! used if arround point source
164  smx=pi/2 -(sourcedec+ddelta)*torad
165  smn=pi/2 -(sourcedec-ddelta)*torad
166  cspsmx=cos(smx)
167  cspsmn=cos(smn)
168  if(za1ry .eq. 'ps' .or. za1ry .eq. 'aps') then
169  write(msg,*) '1ry from point source is specified'
170  call cerrormsg(msg, 1)
171  w3min=real(coszenith)
172 ! from declination and zenith angle, get time
173 ! from meridian
174  call kdzth2(sourcedec, w3min, hmaxd, icon)
175  hmax=hmaxd
176  if(icon .le. 1) then
177  w3max = imag_p(coszenith)
178  call kdzth2(sourcedec, w3max, hmind, jcon)
179  hmin=hmind
180  if(jcon .eq. 2) then
181  hmin=0.
182  call kdhtoh(sourcedec, 0.d0, wh1, wh2, wh3)
183  w3max=wh3
184  endif
185  obsvhour=cmplx(hmin, hmax, 8)
186  if(hmin .gt. hmax) then
187  write(msg,*) ' invalid CosZenith=',coszenith
188  call cerrormsg(msg, 0)
189  endif
190  if(icon .eq. 1) then
191  write(msg,*) ' max observable zenith is',
192  * ' smaller than the input'
193  call cerrormsg(msg, 1)
194  endif
195  write(msg,*)' range of hour =',hmin,' to ',
196  * hmax,' which corresponds to cos(zenith)=',
197  * w3min, w3max
198  call cerrormsg(msg, 1)
199  elseif(icon .eq. 2) then
200 ! no observation range
201  write(msg,*) ' range of cos(zenith)=', coszenith,
202  * ' but no such angle is realizable'
203  call cerrormsg(msg, 0)
204  endif
205  endif
206  end
207 
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
kdhtoh
subroutine kdhtoh(del, h, w1, w2, w3)
Definition: kceles.f:408
cos
! for length to thickness conversion or v v ! integer maxnodes real Hinf ! This is used when making table for dim simulation ! The slant length for vertical height to km is ! divided by LenStep m steps ! It can cover the slant length of about km for cos
Definition: Zatmos.h:8
rndc
subroutine rndc(u)
Definition: rnd.f:91
kdzth2
subroutine kdzth2(del, w3, h, icon)
Definition: kceles.f:488
csprimang
subroutine csprimang(dir)
Definition: csPrimAng.f:14
csprimapsang
subroutine csprimapsang(dir)
Definition: csPrimAng.f:119
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
csprimpsang
subroutine csprimpsang(dir)
Definition: csPrimAng.f:93
csprimisoang
subroutine csprimisoang(dir)
Definition: csPrimAng.f:68
pi
! constants thru Cosmos real * pi
Definition: Zglobalc.h:2
khtad
subroutine khtad(hx, hy, hz, ax, ay, az)
Definition: kceles.f:392
coszenith
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz *Zfirst pos *Zfirst pos Zfirst pos *Zfirst pos *Zfirst Zfirst vec w *Zfirst vec w Zfirst vec coszenith
Definition: ZavoidUnionMap.h:1
kcossn
subroutine kcossn(cs, sn)
Definition: kcossn.f:13
ksamplin
subroutine ksamplin(a, b, alpha, beta, x)
Definition: ksampLin.f:15
coord
Definition: Zcoord.h:43
cinisprimang
subroutine cinisprimang
Definition: csPrimAng.f:150
spherical
const int spherical
Definition: Zobs.h:19