COSMOS v7.655  COSMOSv7655
(AirShowerMC)
cgeomag.f
Go to the documentation of this file.
1 cc test crdGeomag
2 c real*8 year
3 c read(*,*) year
4 c call crdGeomag('../../Data/Geomag/wmm', year)
5 c end
6 c **************************************************************
7 c *
8 c * cgeomag: geomagnetic filed strength is obtained
9 c *
10 c **************************************************************
11 c /usage/ call cgeomag(year, llh, h, icon)
12 c year: real*8. input. such as 1990.5
13 c llh: /coord/ input. position around the earth.
14 c in 'llh' form is better. if not 'llh'
15 c conversion is done here.
16 c h: /magfield/. output. magnetic field is set in
17 c the form of 'ned' (north,
18 c east-down). The unit is T.
19 c icon: output. integer*4 0---> o.k
20 c 1---> input parameter wrong.
21 c
22 
23 c
24  subroutine cgeomag(yearin, llh, h, icon)
25  implicit none
26 #include "Zglobalc.h"
27 #include "Zcoord.h"
28 #include "Zmagfield.h"
29 #include "Zearth.h"
30  real*8 yearin
31  record /coord/ llh
32  record /magfield/ h
33  integer icon
34 
35 c ************************************
36  integer nmax
37  parameter(nmax = 12)
38  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
39  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
40  integer nmx
41  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
42 c ************************************
43 
44 c
45 c
46  real*4 r, sumn, sume, sumd, t, cost, sint, x, tlonr, gmnc,
47  * cosml, sinml, hmnc, temp
48  real*8 gn, ge, gd
49  real*4 ssumd, ssumn, ssume
50  integer m, n
51  record /coord/ cdata
52  real*4 kdpmnxsinn, kpmnxn, kpmnxisinn
53 
54 c
55 c check data type
56  if(llh.sys .eq. 'llh') then
57  cdata = llh
58  else ! convet to llh
59  call ctranscoord2('llh', llh, cdata)
60  endif
61 c
62  icon = 0
63  if( abs(yearin-year) .gt. 7.) then
64  icon=1
65  call cerrormsg(
66  * ' Year spec. for geomag data is invalid.', 1)
67 c elseif(abs(cdata.lat) .gt. 90.) then
68  elseif(abs(cdata.r(1)) .gt. 90.5d0) then
69  icon=2
70 c elseif(abs(cdata.long) .gt. 360.) then
71  elseif(abs(cdata.r(2)) .gt. 360.5d0) then
72  icon=2
73 c elseif(cdata.h .gt. 5000.d3) then
74  elseif(cdata.r(3) .gt. 100000.d3) then
75  icon=1
76 c elseif(cdata.h .lt. -3000.d3) then
77  elseif(cdata.r(3) .lt. -3000.d3) then
78  icon=1
79  endif
80  if(icon .ne. 2) then
81 c r=1./( 1.+cdata.h/Eradius )
82  r=1./( 1.+cdata.r(3)/eradius )
83  sumn=0.
84  sume=0.
85  sumd=0.
86 c t=(90.-cdata.lat)*Torad
87  t=(90.-cdata.r(1))*torad
88  cost=cos(t)
89  sint=sin(t)
90  x=cost
91 c tlonr=cdata.long*Torad
92  tlonr=cdata.r(2)*torad
93  do n=1, nmx
94  ssumn = 0.
95  ssume = 0.
96  ssumd = 0.
97  do m=0, n
98  gmnc= gnm(n, m)
99  cosml=cos(m*tlonr)
100  sinml=sin(m*tlonr)
101  if(m .gt. 0) then
102  hmnc=hnm(n,m)
103  ssumn = ssumn+ (gmnc*cosml+hmnc*sinml)*
104  * kdpmnxsinn(m, n, sint, x)
105  ssume = ssume+ m *(-gmnc*sinml + hmnc*cosml)
106  * * kpmnxisinn(m, n, sint, x)
107  ssumd = ssumd + (gmnc*cosml+ hmnc*sinml)
108  * * kpmnxn(m, n, sint, x)
109  else
110  ssumn = ssumn + gmnc*kdpmnxsinn(m, n, sint, x)
111  ssumd = ssumd + gmnc*kpmnxn(m, n, sint, x)
112  endif
113  enddo
114  temp = r**(n+2)
115  sumn = sumn + ssumn*temp
116  sume = sume + ssume*temp
117  sumd = sumd + ssumd*temp *(n + 1)
118  enddo
119 c original formula gives data in nT.
120 c north component
121  gn = -sumn /1.e9 ! to T.
122 c east component
123  ge = -sume /1.e9
124 c down
125  gd = -sumd /1.e9
126  call csetmagfield('ned', gn, ge, gd, h)
127  endif
128  if(icon .eq. 2) then
129  call cerrormsg('Geometrical input data wrong', 0)
130  endif
131  end
132  subroutine crdgeomag(filepath, yearin)
133  implicit none
134 c #include "Zmanagerp.h"
135  character*(*) filepath
136  real*8 yearin
137  integer klena
138  character*71 msg
139 
140 
141 c ************************************
142  integer nmax
143  parameter(nmax = 12)
144  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
145  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
146  integer nmx
147  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
148 c ************************************
149 
150  character*128 path
151  character*1 NULL
152  integer devn, icon, leng, kgetenv
153  save devn
154  data devn /12/
155 
156  null = char(0)
157 
158  if(filepath .eq. ' ') then
159 c assume it is Cosmos/Data/Geomag/igrf
160 c get $COSMOSTOP
161  leng = kgetenv("COSMOSTOP"//null, path)
162  if(leng .eq. 0) then
163  call cerrormsg(
164  * 'Env. variable "COSMOSTOP" cannot be found', 0)
165  endif
166  path = path(1:leng)//"/Data/Geomag/igrf"
167  else
168  path = filepath
169  endif
170  call copenf(devn, path(1:klena(path)), icon)
171  if(icon .eq. 1) then
172  call cerrormsg(path, 1)
173  call cerrormsg(
174  * 'above file cannot be opened for geomagnetic data', 0)
175  endif
176 c first judge if it is IGRF data or WMM data
177 c and branch
178  read(devn, '(a)') msg
179  if( index(msg, 'DGRF') .gt. 0
180  * .or. index(msg, 'IGRF') .gt. 0 ) then
181 c IGRF data
182 c//////////
183  write(*,*) ' aaaaaaa'
184 c///////////
185  call cgmgigrf(devn, path, yearin)
186  elseif(index(msg, 'igrf') .gt. 0 ) then
187 c IGRF data; from 2005
188 c////////
189  write(*,*) ' bbbbbbbbbbb'
190 c/////////////
191  call cgmgigrf(devn, path, yearin)
192  elseif(index(msg, 'WMM') .gt. 0) then
193 c DoD data
194  call cgmgwmm(devn, path, yearin)
195  else
196  call cerrormsg(msg, 1)
197  call cerrormsg('above data for geomag is strange',0)
198  endif
199  end
200 c *************************
201  subroutine cgmgigrf(devn, filepath, yearin)
202  implicit none
203 
204  integer devn ! input file dev. no.
205  character*(*) filepath
206  real*8 yearin ! input. year of geomag.
207 
208  integer i, j, store, m, n
209  character*72 msg
210  character*90 msg2
211  real*4 yearx, dy, coef, g, h, dg, dh
212 
213 
214 c ************************************
215  integer nmax
216  parameter(nmax = 12)
217  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
218  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
219  integer nmx
220  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
221 c ************************************
222 
223 
224  rewind devn
225  store = 0
226  nmx = 0
227 c first find the appropriate year
228 
229  read(devn, '(a)', end=1000) msg
230  do i=1, 10000
231  do while( index(msg, 'DGRF') .eq. 0 .and.
232  * index(msg, 'IGRF') .eq. 0 .and.
233  * index(msg, 'igrf') .eq. 0 )
234  read(devn, '(a)', end=1000) msg
235  enddo
236  if( index(msg, 'igrf') .ne. 0 ) then
237 c >= 2005
238  read(devn, '(a)', end=1000) msg
239  read(msg, '(12x, f6.1)') yearx
240 c////////
241  write(*,*) ' yearx=',yearx
242 c///////////
243  else
244  read(msg, '(12x, f8.2)') yearx
245  endif
246  if(store .eq. 0 .and. yearin - yearx .lt. 5.0
247  * .and. yearin .ge. yearx) then
248  store = store + 1
249  year = yearx
250  elseif(store .eq. 1 .and. yearin .lt. yearx) then
251  store = store + 1
252  dy = yearx - year
253  coef = ( yearin - year)/dy
254  endif
255 c/////////
256  write(*,*) ' store=',store, 'yearin=',yearin
257 c//////////
258  if( store .gt. 0 ) then
259  do j = 1, 1000000
260  read(devn, '(a)', end=1000) msg
261  if( index(msg, 'DGRF') .gt. 0 .or.
262  * index(msg, 'IGRF') .gt. 0 .or.
263  * index(msg, 'igrf') .gt. 0 ) then
264 c new year appeared
265  goto 30 ! loop out
266  else
267  if(yearx .ge. 2005.) then
268  read(msg, *) n, m, g, h
269  dg = 0.
270  dh = 0.
271  else
272  read(msg, '(2i2, 2f8.1,2f8.2)')
273  * m, n, g, h, dg, dh
274  endif
275  if( n .gt. nmax ) then
276  call cerrormsg(msg, 1)
277  write(msg2,*)
278  * ' n=',n, ' for geomag data > ', nmx
279  call cerrormsg(msg2, 0)
280  endif
281  nmx = max(n, nmx)
282  if(n .le. 0) then
283  call cerrormsg(msg2, 1)
284  write(msg2,*)
285  * ' n=',n, ' for geomag data <=0'
286  call cerrormsg(msg2, 0)
287  endif
288  if(m .lt. 0 .or. m .gt. n) then
289  call cerrormsg(msg2, 1)
290  write(msg2,*)
291  * ' m=',m, ' for geomag data < 0 or > n'
292  call cerrormsg(msg2, 0)
293  endif
294  if(store .eq. 1) then
295  gnm(n, m) = g
296  gnmd(n, m) = dg
297  if(m .gt. 0) then
298  hnm(n, m) = h
299  hnmd(n, m) = dh
300  endif
301  elseif( store .eq. 2) then
302  gnm(n, m) =( g-gnm(n, m)) * coef + gnm(n,m)
303  gnmd(n, m) = 0.
304  if(m .gt. 0) then
305  hnm(n, m) = ( h-hnm(n, m)) * coef + hnm(n,m)
306  hnmd(n, m) = 0 .
307  endif
308  else
309  call cerrormsg(
310  * 'store value error in cgeomag', 0)
311  endif
312  endif
313  enddo
314  30 continue
315  if( store .eq. 2) goto 1000
316  else
317  msg = ' '
318  endif
319  enddo
320  1000 continue
321  if(store .eq. 1) then
322 c last year. do extrapolation and clearn gnmd, hnmd
323  dy = yearin - year
324  do n =1, nmx
325  do m=0, n
326  gnm(n, m) = gnm(n, m) + dy*gnmd(n, m)
327  if(m .gt. 0) then
328  hnm(n, m) = hnm(n,m) + dy*hnmd(n, m)
329  endif
330  enddo
331  enddo
332  endif
333 
334  call cerrormsg('Geomagnetic data has been read from', 1)
335  call cerrormsg(filepath,1)
336  write(msg2, *) ' year=',yearin, ' # of expansion terms=', nmx
337  call cerrormsg(msg2, 1)
338  if(store .eq. 0) then
339  call cerrormsg(
340  * 'Geomagnetic data has no appropriate data', 0)
341  endif
342  end
343 
344 c *************************
345  subroutine cgmgwmm(devn, filepath, yearin)
346  implicit none
347 
348  integer devn ! input file dev. no.
349  character*(*) filepath
350  real*8 yearin ! input. year of geomag.
351 
352  integer i, j, store, m, n
353  character*100 msg
354  real*4 yearx, dy, coef, g, h, dg, dh
355 
356 
357 c ************************************
358  integer nmax
359  parameter(nmax = 12)
360  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
361  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
362  integer nmx
363  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
364 c ************************************
365 
366 
367  rewind devn
368  store = 0
369  nmx = 0
370 c first find the appropriate year
371 
372  read(devn, '(a)', end=1000) msg
373  do i=1, 10000
374  do while(index(msg, 'WMM') .eq. 0)
375  read(devn, '(a)', end=1000) msg
376  enddo
377  read(msg, *) yearx
378  if(store .eq. 0 .and. yearin - yearx .lt. 5.0
379  * .and. yearin .ge. yearx) then
380  store = store + 1
381  year = yearx
382  elseif(store .eq. 1 .and. yearin .lt. yearx) then
383  store = store + 1
384  dy = yearx - year
385  coef = ( yearin - year)/dy
386  endif
387  if( store .gt. 0 ) then
388  do j = 1, 1000000
389  read(devn, '(a)', end=1000) msg
390  if(index(msg, 'WMM') .gt. 0) then
391 c new year appeared
392  goto 30 ! loop out
393  else
394  if(msg(1:3) .eq. '999') then
395  n = 999
396  else
397 c not m, n(cf; IGRF)
398  read(msg, *) n, m, g, h, dg, dh
399  endif
400  if(n .lt. 999) then
401  if( n .gt. nmax ) then
402  call cerrormsg(msg, 1)
403  write(msg,*)
404  * ' n=',n, ' for geomag data > ',
405  * nmax
406  call cerrormsg(msg, 0)
407  endif
408  nmx = max(n, nmx)
409  if(n .le. 0) then
410  call cerrormsg(msg, 1)
411  write(msg,*)
412  * ' n=',n, ' for geomag data <=0'
413  call cerrormsg(msg, 0)
414  endif
415  if(m .lt. 0 .or. m .gt. n) then
416  call cerrormsg(msg, 1)
417  write(msg,*)
418  * ' m=',m, ' for geomag data < 0 or > n'
419  call cerrormsg(msg, 0)
420  endif
421  if(store .eq. 1) then
422  gnm(n, m) = g
423  gnmd(n, m) = dg
424  if( m .gt. 0) then
425  hnm(n, m) = h
426  hnmd(n, m) = dh
427  endif
428  elseif( store .eq. 2) then
429  gnm(n, m) =
430  * ( g-gnm(n, m)) * coef + gnm(n,m)
431  gnmd(n, m) = 0.
432  if(m .gt. 0) then
433  hnm(n, m) =
434  * ( h-hnm(n, m)) * coef + hnm(n,m)
435  hnmd(n, m) = 0.
436  endif
437  else
438  call cerrormsg(
439  * 'store value error in cgeomag', 0)
440  endif
441  endif
442  endif
443  enddo
444  30 continue
445  if(store .eq. 2) goto 1000
446  else
447  msg = ' '
448  endif
449  enddo
450  1000 continue
451  if(store .eq. 1) then
452 c last year. do extrapolation and clearn gnmd, hnmd
453  dy = yearin - year
454  do n =1, nmx
455  do m=0, n
456  gnm(n, m) = gnm(n, m) + dy*gnmd(n, m)
457  if(m .gt. 0) then
458  hnm(n, m) = hnm(n,m) + dy*hnmd(n, m)
459  endif
460  enddo
461  enddo
462  endif
463 
464  call cerrormsg('Geomagnetic data has been read from', 1)
465  call cerrormsg(filepath,1)
466  write(msg, *) ' year=',yearin, ' # of expansion terms=', nmx
467  call cerrormsg(msg, 1)
468  if(store .eq. 0) then
469  call cerrormsg(
470  * 'Geomagnetic data has no appropriate data', 0)
471  endif
472  end
473 
474 
475 c *************************
476  subroutine cprgeomag
477  implicit none
478 c print current geomag coeff.
479 
480 
481  integer m, n
482 
483 
484 c ************************************
485  integer nmax
486  parameter(nmax = 12)
487  real*4 gnm(nmax, 0:nmax), hnm(nmax, nmax), year
488  real*4 gnmd(nmax, 0:nmax), hnmd(nmax, nmax)
489  integer nmx
490  common /zmagcoef/ gnm, gnmd, hnm, hnmd, year, nmx
491 c ************************************
492 
493  write(*,*) 'year=',year, ' nmax=', nmx
494  write(*,'(" m n gnm hnm dgnm dhnm")')
495  do n = 1, nmx
496  do m =0, n
497  if(m .eq. 0) then
498  write(*,'(2i3,2f7.0,2f7.1)') m, n,
499  * gnm(n, m), 0., gnmd(n, m), 0.
500  else
501  write(*,'(2i3,2f7.0,2f7.1)') m, n,
502  * gnm(n, m), hnm(n,m), gnmd(n,m), hnmd(n,m)
503  endif
504  enddo
505  enddo
506  end
507 
508 
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
cprgeomag
subroutine cprgeomag
Definition: cgeomag.f:477
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
crdgeomag
subroutine crdgeomag(filepath, yearin)
Definition: cgeomag.f:133
d3
block data cblkIncident data *Za1ry *HeightOfInj d3
Definition: cblkIncident.h:7
llh
latitude latitude this system is used *****************************************************************! type coord sequence union map real z z in m endmap xyz map real ! latitude in deg is to the north ! longitude in deg is to the east *h ! height in m endmap llh map real ! polar angle ! azimuthal angle *radius ! radial distance endmap sph endunion character *sys ! which system llh
Definition: Zcoord.h:25
cgeomag
subroutine cgeomag(yearin, llh, h, icon)
Definition: cgeomag.f:25
sys
*Zfirst p fm *Zfirst p Zfirst p Zfirst p *Zfirst p *Zfirst pos xyz Zfirst pos xyz sys
Definition: ZavoidUnionMap.h:1
cgmgwmm
subroutine cgmgwmm(devn, filepath, yearin)
Definition: cgeomag.f:346
ctranscoord2
subroutine ctranscoord2(sys, a, b)
Definition: ctransCoord2.f:2
cgmgigrf
subroutine cgmgigrf(devn, filepath, yearin)
Definition: cgeomag.f:202
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
h
dE dx *! Nuc Int sampling table h
Definition: cblkMuInt.h:130
csetmagfield
subroutine csetmagfield(sys, b1, b2, b3, b)
brief set Calculated magnetic field to /magfield/ b
Definition: csetMagField.f:6
coord
Definition: Zcoord.h:43