COSMOS v7.655  COSMOSv7655
(AirShowerMC)
Functions/Subroutines
cbDefByRK.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine cbdefbyrk (aTrack, leng, newpos, newdir, newmom)
 Compute magnetic deflection (angle and displacement) Use Runge-Kutta-Fehlberg Method. More...
 
subroutine cgeomfunc (L, y, f)
 
subroutine krungekutfds (neq, func, x0, xe, y0, init, relerr, abserr, yn, esterr, nde, ier, h, yl0, yln, err, work)
 
subroutine krungekutfs (neq, func, x0, xe, yh0, init, relerr, abserr, yhn, esterr, nde, ier, h, yl0, yln, err, ak1, ak2, ak3, ak4, ak5, ak6, et, work)
 
subroutine krkflstep (neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, work)
 
subroutine krkfhstep (neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, work)
 
subroutine krkfehl (lindex, neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, w2, w3, w4, w5, w6)
 
subroutine cbdefbyrk2 (aTrack, leng, newpos, newdir, newmom)
 
subroutine krungekutg (func, n, x0, y0, m, h, ny, y, wk, ierr)
 
subroutine cbdefuser (aTrack, leng, newpos, newdir, newmom)
 

Function/Subroutine Documentation

◆ cbdefbyrk()

subroutine cbdefbyrk ( type(track)  aTrack,
real*8  leng,
type(coord)  newpos,
type(coord)  newdir,
real(8), dimension(3), intent(out)  newmom 
)

Compute magnetic deflection (angle and displacement) Use Runge-Kutta-Fehlberg Method.

Parameters
[in]aTracka charged ptcl (datatype /track/)
[in]lenglength travelled in m (datatype real*8)
[out]newposnew position after leng move (datatype /coord/newpos)
[out]newmomchanged mon

Definition at line 12 of file cbDefByRK.f.

References cbdefbyrk2(), cgeomfunc(), d, d0, d3, and krungekutfds().

Referenced by cmagdef().

12 ! ************************************************
13 ! This uses Runge-Kutta method with the automatically adjustable
14 ! step size and takes a long time for execution.
15 !
16 
17 
18  implicit none
19 
20 #include "Ztrack.h"
21 
22  real(8),intent(out):: newmom(3) ! chnaged mom.
23  real*8 geomcnst ! Z/pc/3.3358. Z charge, pc GeV
24  common /zgeomcnst/ geomcnst
25 
26  type(track)::atrack !input. a charged ptcl
27 
28 ! ptcl and magfiled coord is in Exyz
29 
30  real*8 leng !input. length travelled in m
31 
32  type(coord)::newpos ! output. new position after leng move
33  type(coord)::newdir ! // new direction cose. //
34 
35  real*8 cbetat, x0
36 
37  integer first, nde, ierr
38  real*8 relaerr, abserr, h
39  real*8 esterr(6), yn(6), y10(6), y1n(6), err(6), y0(6)
40  real*8 work(6,12)
41  save abserr, relaerr
42 
43  external cgeomfunc
44 
45  data first/1/, relaerr/1.d-2/, abserr/1.d3/
46 
47  cbetat= leng
48 !
49 ! y0(1) = aTrack.pos.xyz.x
50 ! y0(2) = aTrack.pos.xyz.y
51 ! y0(3) = aTrack.pos.xyz.z
52 ! IBM must be like below
53  y0(1) = atrack%pos%xyz%r(1)
54  y0(2) = atrack%pos%xyz%r(2)
55  y0(3) = atrack%pos%xyz%r(3)
56 
57 ! y0(4)= aTrack%vec%w%x
58 ! y0(5)= aTrack%vec%w%y
59 !　 y0(6)= aTrack%vec%w%z
60  y0(4:6) = atrack%vec%w%r(1:3)
61 
62 ! pc in GeV
63  geomcnst = max(sqrt(atrack%p%fm%p(4)**2- atrack%p%mass**2), 1.d-9)
64 
65  geomcnst = atrack%p%charge/geomcnst/3.3358
66 
67  x0= 0.d0 ! this cannot be a literal below; due to bad prog.
68  call krungekutfds(6, cgeomfunc, x0, cbetat, y0, first, relaerr,
69  * abserr, yn, esterr, nde, ierr, h, y10, y1n, err, work)
70  if(ierr .eq. 40000 ) then
71  call cbdefbyrk2(atrack, leng, newpos, newdir, newmom)
72  elseif(ierr .ne. 0) then
73  write(0,*) ' ierr=',ierr
74  call cbdefbyrk2(atrack, leng, newpos, newdir, newmom)
75  else
76 
77 !c first = 1
78 
79 ! newpos.x = yn(1)
80 ! newpos.y = yn(2)
81 ! newpos.z = yn(3)
82 ! for ibm
83  newpos%r(1) = yn(1)
84  newpos%r(2) = yn(2)
85  newpos%r(3) = yn(3)
86  newpos%sys='xyz'
87 ! newdir.x = yn(4)
88 ! newdir.y = yn(5)
89 ! newdir.z = yn(6)
90 ! for ibm
91  newdir%r(1) = yn(4)
92  newdir%r(2) = yn(5)
93  newdir%r(3) = yn(6)
94 ! oh my god next yn(4:6)* sqrt(. ) was
95 ! yn(4:6)** sqrt(...) !!!
96 ! but default UseRungeKutta = 0 did not come here
97  newmom(:) = yn(4:6)*
98  * sqrt(dot_product( atrack%p%fm%p(1:3),atrack%p%fm%p(1:3)))
99  newdir%sys='xyz'
100  endif
cgeomfunc
subroutine cgeomfunc(L, y, f)
Definition: cbDefByRK.f:103
nrl_atmos::first
logical, save first
Definition: cNRLAtmos.f:8
d3
block data cblkIncident data *Za1ry *HeightOfInj d3
Definition: cblkIncident.h:7
track
Definition: Ztrack.h:44
leng
integer leng
Definition: interface2.h:1
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
cbdefbyrk2
subroutine cbdefbyrk2(aTrack, leng, newpos, newdir, newmom)
Definition: cbDefByRK.f:473
coord
Definition: Zcoord.h:43
krungekutfds
subroutine krungekutfds(neq, func, x0, xe, y0, init, relerr, abserr, yn, esterr, nde, ier, h, yl0, yln, err, work)
Definition: cbDefByRK.f:145
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cbdefbyrk2()

subroutine cbdefbyrk2 ( type(track)  aTrack,
real*8  leng,
type(coord)  newpos,
type(coord)  newdir,
real(8), dimension(3), intent(out)  newmom 
)

Definition at line 473 of file cbDefByRK.f.

References cgeomfunc(), d, d0, krungekutg(), and parameter().

Referenced by cbdefbyrk(), and cmagdef().

473 ! ************************************************
474 ! This used Runge-Kutta-Gill method. The step size
475 ! is fixed.
476 !
477 
478  implicit none
479 
480 #include "Ztrack.h"
481  real(8),intent(out):: newmom(3) ! changed mom.
482  real*8 geomcnst ! Z/pc/3.3358. Z charge, pc GeV
483  common /zgeomcnst/ geomcnst
484 
485  type(track)::atrack !input. a charged ptcl
486 
487 ! ptcl and magfiled coord is in Exyz
488 
489  real*8 leng !input. length travelled in m
490  ! since we use Ruuge-Kutta-Gill method with fixed
491  ! step size, we use leng/2 as the step size assuming
492  ! leng is (Larmor radius)/5.
493 
494  type(coord)::newpos ! output. new position after leng move
495  type(coord)::newdir ! // new direction cose. //
496 
497  real*8 cbetat, x0
498 
499  integer ier, m
500 
501  real*8 y0(6)
502  parameter(m=3)
503  real*8 wk(6,4), y(6,m)
504 
505  external cgeomfunc
506 
507 
508  cbetat= leng/(m-1)
509 !
510 ! y0(1) = aTrack.pos.xyz.x
511 ! y0(2) = aTrack.pos.xyz.y
512 ! y0(3) = aTrack.pos.xyz.z
513  y0(1:3) = atrack%pos%xyz%r(1:3)
514 
515 ! y0(4)= aTrack%vec%w%x
516 ! y0(5)= aTrack%vec%w%y
517 ! y0(6)= aTrack%vec%w%z
518  y0(4:6) = atrack%vec%w%r(1:3)
519 
520 ! pc in GeV
521  geomcnst = max(sqrt(atrack%p%fm%p(4)**2- atrack%p%mass**2), 1.d-9)
522 
523  geomcnst = atrack%p%charge/geomcnst/3.3358
524 
525  x0 = 0.d0
526  call krungekutg(cgeomfunc, 6, x0, y0, m, cbetat, 6, y, wk, ier)
527  if(ier .ne. 0) then
528  write(0,*) ' ier=',ier
529  endif
530 
531 
532 ! newpos.x = y(1,m)
533 ! newpos.y = y(2,m)
534 ! newpos.z = y(3,m)
535 ! for ibm
536  newpos%r(1) = y(1,m)
537  newpos%r(2) = y(2,m)
538  newpos%r(3) = y(3,m)
539  newpos%sys='xyz'
540 ! newdir.r(1) = y(4,m)
541 ! newdir.r(2) = y(5,m)
542 ! newdir.r(3) = y(6,m)
543 ! for ibm
544  newdir%r(1) = y(4,m)
545  newdir%r(2) = y(5,m)
546  newdir%r(3) = y(6,m)
547  newmom(:) = y(4:6,m)*
548  * sqrt( dot_product(atrack%p%fm%p(1:3),atrack%p%fm%p(1:3)) )
549  newdir%sys='xyz'
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
cgeomfunc
subroutine cgeomfunc(L, y, f)
Definition: cbDefByRK.f:103
track
Definition: Ztrack.h:44
leng
integer leng
Definition: interface2.h:1
krungekutg
subroutine krungekutg(func, n, x0, y0, m, h, ny, y, wk, ierr)
Definition: cbDefByRK.f:554
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35
coord
Definition: Zcoord.h:43
Here is the call graph for this function:
Here is the caller graph for this function:

◆ cbdefuser()

subroutine cbdefuser ( type(track)  aTrack,
real*8  leng,
type(coord)  newpos,
type(coord)  newdir,
real(8), intent(out)  newmom 
)

Definition at line 641 of file cbDefByRK.f.

References d.

Referenced by cmagdef().

641 ! ************************************************
642 !
643 ! another method to get the new position and direction
644 ! cosines when magnetic field exists.
645 ! This is called when UseRungeKutta is 8.
646 
647  implicit none
648 
649 #include "Ztrack.h"
650 #include "Zcode.h"
651 
652  real*8 geomcnst ! Z/pc/3.3358. Z charge, pc GeV
653  common /zgeomcnst/ geomcnst
654 
655  type(track)::atrack ! input. a charged ptcl
656  real(8),intent(out):: newmom ! changed momentum
657 !
658 ! current position
659 ! aTrack.pos.xyz.x, aTrack.pos.xyz.y, aTrack.pos.xyz.z (m)
660 ! current direction cosines
661 ! aTrack.vec.w.x, aTrack.vec.w.y, aTrack.vec.w.z
662 ! current momentum, total energy (GeV)
663 ! aTrack.p.fm.p(i) i=1,4
664 ! mass patcl code subcode charge
665 ! aTrack.p.mass aTrack.p.code aTrack.p.subcode,aTrack.p.charge
666 ! if code = kgnuc (isotope), subcode has mass number.
667 !
668 
669 
670 
671 ! ptcl and magfiled coord is in Exyz
672 
673  real*8 leng ! input. length travelled in m
674 
675  type(coord)::newpos ! output. new position after moving leng (m)
676  !
677  ! newpos.x newpos.y, newpos.z must be given
678  ! newpos.sys='xyz' must be gieven
679  type(coord)::newdir ! // new direction cose. //
680  ! newdir.x, newdir.y, newdir.z must be given
681  ! newdir.sys='xyz' must be given
682 
683 ! pc in GeV
684  geomcnst = max(sqrt(atrack%p%fm%p(4)**2- atrack%p%mass**2), 1.d-9)
685  geomcnst = atrack%p%charge/geomcnst/3.3358
686 !c
687 !c You may see cgeomfunc how to get magnetic field etc
688 !c
689 ! Note: cgeomag routine below gives B in north, east, down
690 ! coordinate system, so that it must be converted into
691 ! E-xyz system by calling ctransMagTo like below.
692 !
693 ! call cgeomag(YearOfGeomag, llh, B, icon) ! B is ned
694 ! call ctransMagTo('xyz', llh, B, B) ! to xyx
695 !
696 
697  newpos%sys = 'xyz'
698  newdir%sys = 'xyz'
track
Definition: Ztrack.h:44
leng
integer leng
Definition: interface2.h:1
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
coord
Definition: Zcoord.h:43
Here is the caller graph for this function:

◆ cgeomfunc()

subroutine cgeomfunc ( real*8  L,
real*8, dimension(6)  y,
real*8, dimension(6)  f 
)

Definition at line 103 of file cbDefByRK.f.

References cgeomag(), and ctransmagto().

Referenced by cbdefbyrk(), and cbdefbyrk2().

103  implicit none
104 ! this is used both by cbDefByRK and cbDefByRK2
105 ! Zobs.h is For Zobsp.h
106 ! Zobsp.h is For YearOfGeomag
107 #include "Zcoord.h"
108 #include "Zmagfield.h"
109 #include "Zobs.h"
110 #include "Zobsp.h"
111 
112  real*8 geomcnst ! Z/pc/3.3358. Z charge, pc GeV
113  common /zgeomcnst/ geomcnst
114 
115  type(coord)::llh
116  type(magfield)::b
117  integer icon
118  real*8 l, y(6), f(6)
119 ! L=cbetat is not used explicitly
120 
121  f(1) = y(4)
122  f(2) = y(5)
123  f(3) = y(6)
124 ! compute B at y(1),y(2), y(3)
125 !
126 ! llh.x = y(1)
127 ! llh.y = y(2)
128 ! llh.z = y(3)
129 ! for ibm
130  llh%r(1) = y(1)
131  llh%r(2) = y(2)
132  llh%r(3) = y(3)
133 
134  llh%sys = 'xyz' ! converted to 'llh' in cgeomag
135  call cgeomag(yearofgeomag, llh, b, icon) ! B is ned
136  call ctransmagto('xyz', llh, b, b) ! to xyx
137  f(4)= geomcnst * (y(5)*b%z - y(6)*b%y)
138  f(5) =geomcnst * (y(6)*b%x - y(4)*b%z)
139  f(6) =geomcnst * (y(4)*b%y - y(5)*b%x)
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
y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
ctransmagto
subroutine ctransmagto(sys, pos, a, b)
Definition: ctransMagTo.f:11
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
coord
Definition: Zcoord.h:43
f
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krkfehl()

subroutine krkfehl (   lindex,
  neq,
  func,
  x,
  h,
dimension(neq)  y0,
dimension(neq)  yn,
dimension(neq)  ak1,
dimension(neq)  ak2,
dimension(neq)  ak3,
dimension(neq)  ak4,
dimension(neq)  ak5,
dimension(neq)  ak6,
dimension(neq)  err,
dimension(neq)  w2,
dimension(neq)  w3,
dimension(neq)  w4,
dimension(neq)  w5,
dimension(neq)  w6 
)

Definition at line 411 of file cbDefByRK.f.

References a, d, h, i, o, parameter(), x, and z.

Referenced by krkfhstep(), and krkflstep().

411  implicit real*8(a-h,o-z)
412  dimension y0(neq),yn(neq),ak1(neq),ak2(neq),ak3(neq),ak4(neq),
413  & ak5(neq),ak6(neq),err(neq),w2(neq),w3(neq),w4(neq),
414  & w5(neq),w6(neq)
415  parameter(one = 1, two = 2, thr = 3, twl = 12)
416  parameter(al2 = one / 4, al3 = thr / 8, al4 = twl / 13,
417  & al5 = one, al6 = one / 2)
418  parameter(b21 = one / 4, b31 = thr / 32, b32 = 9.0d+0 / 32,
419  & b41 = 1932.0d+0 / 2197, b42 = -7200.0d+0 / 2197,
420  & b43 = 7296.0d+0 / 2197, b51 = 439.0d+0 / 216,
421  & b52 = -8, b53 = 3680.0d+0 / 513,
422  & b54 = -845.0d+0 / 4104, b61 = -8.0d+0 / 27,
423  & b62 = 2, b63 = -3544.0d+0 / 2565,
424  & b64 = 1859.0d+0 / 4104, b65 = -11.0d+0 / 40)
425  parameter(ga1 = 16.0d+0 / 135, ga3 = 6656.0d+0 / 12825,
426  & ga4 = 28561.0d+0 / 56430,
427  & ga5 = -9.0d+0 / 50, ga6 = two / 55)
428  parameter(da1 = one / 360, da3 = -128.0d+0 / 4275,
429  & da4 = -2197.0d+0 / 75240,
430  & da5 = one / 50, da6 = two / 55)
431 
432  call func(x,y0,ak1)
433  do 10 i = 1,neq
434  w2(i) = y0(i) + h * b21 * ak1(i)
435  10 continue
436  call func(x + al2 * h,w2,ak2)
437  do 20 i = 1,neq
438  w3(i) = y0(i) + h * (b31 * ak1(i) + b32 * ak2(i))
439  20 continue
440  call func(x + al3 * h,w3,ak3)
441  do 30 i = 1,neq
442  w4(i) = y0(i) + h * (b41 * ak1(i) + b42 * ak2(i) + b43 * ak3(i))
443  30 continue
444 
445  call func(x + al4 * h,w4,ak4)
446  do 40 i = 1,neq
447  w5(i) = y0(i) + h * (b51 * ak1(i) + b52 * ak2(i)
448  & + b53 * ak3(i) + b54 * ak4(i))
449  40 continue
450 
451  call func(x + al5 * h,w5,ak5)
452  do 50 i = 1,neq
453  w6(i) = y0(i) + h * (b61 * ak1(i) + b62 * ak2(i) + b63 * ak3(i)
454  & + b64 * ak4(i) + b65 * ak5(i))
455  50 continue
456  call func(x + al6 * h,w6,ak6)
457  do 60 i = 1,neq
458  yn(i) = y0(i) + h * (ga1 * ak1(i) + ga3 * ak3(i) + ga4 * ak4(i)
459  & + ga5 * ak5(i) + ga6 * ak6(i))
460  60 continue
461  if (lindex .eq. 1) then
462  do 70 i = 1,neq
463  err(i) = h * (da1 * ak1(i) + da3 * ak3(i) + da4 * ak4(i)
464  & + da5 * ak5(i) + da6 * ak6(i))
465  70 continue
466  endif
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
z
nodes z
Definition: cstdatmos0.fNew2.h:93
i
nodes i
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
o
struct ob o[NpMax]
Definition: Zprivate.h:34
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krkfhstep()

subroutine krkfhstep (   neq,
external  func,
  x,
  h,
dimension(neq)  y0,
dimension(neq)  yn,
dimension(neq)  ak1,
dimension(neq)  ak2,
dimension(neq)  ak3,
dimension(neq)  ak4,
dimension(neq)  ak5,
dimension(neq)  ak6,
dimension(neq)  err,
dimension(neq,5)  work 
)

Definition at line 390 of file cbDefByRK.f.

References a, d, h, i, krkfehl(), o, x, x1(), and z.

Referenced by krungekutfs().

390  implicit real*8(a-h,o-z)
391  external func
392  dimension y0(neq),yn(neq),ak1(neq),ak2(neq),ak3(neq),ak4(neq),
393  & ak5(neq),ak6(neq),err(neq),work(neq,5)
394  x1 = x
395  h1 = 0.5d+0 * h
396  lindex = 0
397 
398  call krkfehl(lindex,neq,func,x1,h1,y0,yn,ak1,ak2,ak3,ak4,ak5,ak6,
399  & err,work(1,1),work(1,2),work(1,3),work(1,4),work(1,5))
400  x1 = x1 + h1
401  do 10 i = 1,neq
402  y0(i) = yn(i)
403  10 continue
404  call krkfehl(lindex,neq,func,x1,h1,y0,yn,ak1,ak2,ak3,ak4,ak5,ak6,
405  & err,work(1,1),work(1,2),work(1,3),work(1,4),work(1,5))
406  return
z
nodes z
Definition: cstdatmos0.fNew2.h:93
x1
block data include Zlatfit h c fitting region data x1(1)/0.03/
i
nodes i
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
krkfehl
subroutine krkfehl(lindex, neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, w2, w3, w4, w5, w6)
Definition: cbDefByRK.f:411
o
struct ob o[NpMax]
Definition: Zprivate.h:34
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krkflstep()

subroutine krkflstep (   neq,
external  func,
  x,
  h,
dimension(neq)  y0,
dimension(neq)  yn,
dimension(neq)  ak1,
dimension(neq)  ak2,
dimension(neq)  ak3,
dimension(neq)  ak4,
dimension(neq)  ak5,
dimension(neq)  ak6,
dimension(neq)  err,
dimension(neq,5)  work 
)

Definition at line 377 of file cbDefByRK.f.

References a, h, krkfehl(), o, x, and z.

Referenced by krungekutfs().

377  implicit real*8(a-h,o-z)
378  external func
379  dimension y0(neq),yn(neq),ak1(neq),ak2(neq),ak3(neq),ak4(neq),
380  & ak5(neq),ak6(neq),err(neq),work(neq,5)
381  lindex = 1
382  call krkfehl(lindex,neq,func,x,h,y0,yn,ak1,ak2,ak3,ak4,ak5,ak6,
383  & err,work(1,1),work(1,2),work(1,3),work(1,4),work(1,5))
384 
385  return
z
nodes z
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
krkfehl
subroutine krkfehl(lindex, neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, w2, w3, w4, w5, w6)
Definition: cbDefByRK.f:411
o
struct ob o[NpMax]
Definition: Zprivate.h:34
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krungekutfds()

subroutine krungekutfds (   neq,
external  func,
  x0,
  xe,
dimension(neq)  y0,
  init,
  relerr,
  abserr,
dimension(neq)  yn,
dimension(neq)  esterr,
  nde,
  ier,
  h,
dimension(neq)  yl0,
dimension(neq)  yln,
dimension(neq)  err,
dimension(neq,12)  work 
)

Definition at line 145 of file cbDefByRK.f.

References a, h, init(), krungekutfs(), o, and z.

Referenced by cbdefbyrk().

145 !
146 ! adapted from Maruzen Library.
147 !
148 !*********************************************************************
149 ! subroutine krungekutfds numerically integrates a system of neq *
150 ! first order ordinary differential equations of the form *
151 ! dy(i)/dx = f(x, y(1),..., y(neq)), *
152 ! by the runge-kutta-fehlberg (4,5) formula, *
153 ! subroutine rkfds can detect the stiffness of the differential *
154 ! system. *
155 ! *
156 ! parameters *
157 ! === input === *
158 ! (1) neq: number of equations to be integrated *
159 ! (2) func: subroutine func(x,y,f) to evaluate derivatives *
160 ! f(i)=dy(i)/dx *
161 ! (3) x0: initial value of independent variable *
162 ! (4) xe: output point at which the solution is desired *
163 ! (5) y0(i) (i=1,..,neq): initial value at x0 *
164 ! (6) init: indicator to initialize the code *
165 ! init=1..the code will be initialized(first call). *
166 ! init=2..the code will not be initialized(subsequent call).*
167 ! (7) relerr: relative local error tolerance *
168 ! (8) abserr: absolute local error tolerance *
169 ! === output === *
170 ! (9) yn(i) (i=1,..,neq): approximate solution at xe *
171 ! (10) esterr(i) (i=1,..,neq): estimate of the global error in *
172 ! the approximate solution yhn(i) *
173 ! (11) nde: number of derivative evaluations *
174 ! (12) ier: indicator for status of integration *
175 ! ier=0..integration reached xe. indicates successful *
176 ! return and is the normal mode for continuing the *
177 ! integration. *
178 ! ier=10000..integration was not completed because too many *
179 ! derivatives evaluations were needed. if the user wants*
180 ! to continue the integration, he just calls rkf again. *
181 ! ier=20000..integration was not completed because error *
182 ! tolerance was inappropriate(=it was zero). must use *
183 ! non-zero abserr to continue the integration. *
184 ! ier=30000..integration was not completed because requested*
185 ! accuracy could not be achieved using smallest stepsize.*
186 ! must increase the error tolerance to continue the *
187 ! integration. *
188 ! ier=40000..integration was not completed because the *
189 ! stiffness of the differential system was detected. *
190 ! must use the subroutine, krsnat or krsat, which *
191 ! is designed to integrate stiff differential systems. *
192 ! === others === *
193 ! (13) h: variable to hold information internal to rkfds which is *
194 ! necessary for subsequent calls. *
195 ! (14) yl0(), yln(): array (size=neq) to hold information internal*
196 ! to rkfds which is necessary for subsequent calls. *
197 ! (15) err(): array (size=neq) to be used inside rkfds *
198 ! (16) work(): two-dimentional array (size=(neq,12)) to be *
199 ! used inside rkfds *
200 ! copyright: m. sugihara, november 15, 1989, v. 1 *
201 !*********************************************************************
202  implicit real*8(a-h,o-z)
203  external func
204  dimension y0(neq),yn(neq),esterr(neq),work(neq,12),
205  & yl0(neq),yln(neq),err(neq)
206  call krungekutfs(neq,func,x0,xe,y0,init,relerr,abserr,
207 ! Next original one is wrong. //// corrected by K.K
208 ! & yn,esterr,nde,ier,h,yl0(neq),yln(neq),err(neq),
209  & yn,esterr,nde,ier,h, yl0, yln, err,
210  & work(1,1),work(1,2),work(1,3),work(1,4),work(1,5),
211  & work(1,6),work(1,7),work(1,8))
212  return
z
nodes z
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
init
subroutine init
Definition: Gencol.f:62
o
struct ob o[NpMax]
Definition: Zprivate.h:34
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
krungekutfs
subroutine krungekutfs(neq, func, x0, xe, yh0, init, relerr, abserr, yhn, esterr, nde, ier, h, yl0, yln, err, ak1, ak2, ak3, ak4, ak5, ak6, et, work)
Definition: cbDefByRK.f:218
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krungekutfs()

subroutine krungekutfs (   neq,
external  func,
  x0,
  xe,
dimension(neq)  yh0,
  init,
  relerr,
  abserr,
dimension(neq)  yhn,
dimension(neq)  esterr,
  nde,
  ier,
  h,
dimension(neq)  yl0,
dimension(neq)  yln,
dimension(neq)  err,
dimension(neq)  ak1,
dimension(neq)  ak2,
dimension(neq)  ak3,
dimension(neq)  ak4,
dimension(neq)  ak5,
dimension(neq)  ak6,
dimension(neq)  et,
dimension(neq,5)  work 
)

Definition at line 218 of file cbDefByRK.f.

References a, d, h, i, init(), krkfhstep(), krkflstep(), o, parameter(), x, and z.

Referenced by krungekutfds().

218  implicit real*8(a-h,o-z)
219  parameter(c1 = 0.055455d+0, c2 = -0.035493d+0,
220  & c3 = -0.036571d+0, c4 = 0.023107d+0,
221  & c5 = -9.515d-3, c6 = 3.017d-3 )
222  external func
223  dimension yh0(neq),yhn(neq),yl0(neq),yln(neq),esterr(neq),
224  & ak1(neq),ak2(neq),ak3(neq),ak4(neq),ak5(neq),ak6(neq),
225  & err(neq),et(neq),work(neq,5)
226  data itemax / 100000 /
227  data epsmin / 1.0d-15 /
228 !
229  if (init .eq. 1) then
230 ! -------------- initialization (first call)-------------------------
231  do 10 i = 1,neq
232  yl0(i) = yh0(i)
233  10 continue
234 ! -------- set initial step size --------
235 
236  call func(x0,yl0,ak1)
237  toldy = 10.0d+74
238  do 20 i = 1,neq
239  if (ak1(i) .ne. 0.0d+0) then
240  tol = relerr * dabs(yl0(i)) + abserr
241  toldy = min(toldy, tol / dabs(ak1(i)))
242  end if
243  20 continue
244  hmin = epsmin * (xe - x0)
245  if (toldy .eq. 10.0d+74) then
246  h = hmin
247  else
248  h = min((xe - x0) / 100.0d+0, toldy ** 0.2d+0)
249  end if
250 ! ----------------------------------------
251  init = 2
252  nde = 1
253  else
254  nde = 0
255  endif
256 ! -------------------------------------------------------------------
257  istiff = 0
258  icount = 0
259  ier = 0
260  x = x0
261 !
262 !********************** main iteration *********************************
263  do 30 iter = 1,itemax
264  call krkflstep(neq,func,x,h,yl0,yln,ak1,ak2,ak3,ak4,ak5,
265  & ak6,err,work)
266  nde = nde + 6
267  erret = 0.0d+0
268  do 40 i = 1,neq
269  et(i) = relerr * 0.5d+0 * (dabs(yl0(i)) + dabs(yln(i))) + abserr
270  if (et(i) .eq. 0.0d+0) then
271  go to 20000
272  else
273  erret = max(erret, dabs(err(i)) / et(i))
274  end if
275  40 continue
276  if (erret .ge. 1.0d+0) then
277 !--------------- unsuccessful step ------------------------------
278  if (erret .ge. 59049.0d+0) then
279  h = 0.1d+0 * h
280  else
281  h = 0.9d+0 * h / (erret ** 0.2d+0)
282  end if
283  if (h .le. epsmin) go to 30000
284  else if ((x + h) .lt. xe) then
285 !-------------- successful step (x + h < the end point)------------
286 ! ********* routine for detecting stiffness *********
287  icount = icount + 1
288  erretl = 0.0d+0
289  do 100 i = 1,neq
290  erri = h * (c1 * ak1(i) + c2 * ak2(i) + c3 * ak3(i)
291  & + c4 * ak4(i) + c5 * ak5(i) + c6 * ak6(i))
292  erretl = max(erretl, dabs(erri) / et(i))
293  100 continue
294  if (erretl .le. 1.0d+0) istiff = istiff + 1
295  if (icount .eq. 50) then
296  if (istiff .ge. 25) then
297 ! ------ stiffness was detected ------
298  go to 40000
299  else
300  icount = 0
301  istiff = 0
302  end if
303  end if
304 ! **************************************************
305 
306  call krkfhstep(neq,func,x,h,yh0,yhn,ak1,ak2,ak3,ak4,ak5,
307  & ak6,err,work)
308  nde = nde + 12
309  x = x + h
310  x0 = x
311 
312  do 50 i = 1,neq
313  yl0(i) = yln(i)
314  yh0(i) = yhn(i)
315  50 continue
316 
317 ! ------- choose next step -----------
318  if (erret .le. 1.889568d-4) then
319  h = 5.0d+0 * h
320  else
321  h = 0.9d+0 * h / (erret ** 0.2d+0)
322  end if
323 ! ------------------------------------
324  else
325 !-------------- successful step (x + h > = the end point) --------------
326  h = xe-x
327  call krkflstep(neq,func,x,h,yl0,yln,ak1,ak2,ak3,ak4,ak5,
328  & ak6,err,work)
329  call krkfhstep(neq,func,x,h,yh0,yhn,ak1,ak2,ak3,ak4,ak5,
330  & ak6,err,work)
331 
332  nde = nde + 18
333 ! -------- estimate global error ----------
334  do 60 i = 1,neq
335  esterr(i) = (yln(i) - yhn(i)) / 31.0d+0
336  60 continue
337 ! -----------------------------------------
338  x0 = xe
339  do 70 i = 1,neq
340  yl0(i) = yln(i)
341  yh0(i) = yhn(i)
342  70 continue
343  return
344  endif
345  30 continue
346 !*********************************************************************
347 !
348  ier = 10000
349  write( * ,10001) x
350 10001 format(' ','(subr.-rkfds) trouble(too many iterations))',
351  & 'at x = ',1pe15.7)
352  return
353 20000 continue
354  ier = 20000
355  write( * ,20001) x
356 20001 format(' ','(subr.-rkfds) trouble(inappropriate error tolerance)'
357  & ,' at x = ',1pe15.7)
358  return
359 30000 continue
360  ier = 30000
361  write( * ,30001) x
362 30001 format(' ','(subr.-rkfds) trouble(too small step size)',
363  & ' at x = ',1pe15.7)
364  return
365 40000 continue
366  ier = 40000
367  write( 0 , * ) 'the equations are stiff!'
368  write( 0 ,40001) x
369 40001 format(' ','stiffness was detected at x = ',1pe15.7,'.')
370  write( 0 , * ) 'to use the subroutine (krsnat or krsat)'
371  &,' is recommened.'
372  return
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
z
nodes z
Definition: cstdatmos0.fNew2.h:93
i
nodes i
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
init
subroutine init
Definition: Gencol.f:62
o
struct ob o[NpMax]
Definition: Zprivate.h:34
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
krkfhstep
subroutine krkfhstep(neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, work)
Definition: cbDefByRK.f:390
krkflstep
subroutine krkflstep(neq, func, x, h, y0, yn, ak1, ak2, ak3, ak4, ak5, ak6, err, work)
Definition: cbDefByRK.f:377
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the call graph for this function:
Here is the caller graph for this function:

◆ krungekutg()

subroutine krungekutg ( external  func,
integer  n,
real*8  x0,
real*8, dimension(n)  y0,
integer  m,
real*8  h,
integer  ny,
real*8, dimension(ny,m)  y,
real*8, dimension(n,4)  wk,
integer  ierr 
)

Definition at line 554 of file cbDefByRK.f.

References d0.

Referenced by cbdefbyrk2().

554  implicit none
555 !
556 ! solvoes n sets of dy/dx = f(x,y)
557 !
558  external func ! input. func is a external subroutine name
559  integer n ! input. number of defferential equations
560  real*8 x0 ! input. initial value.
561  real*8 y0(n) ! input. initial value.
562  integer m ! input. total number of points where the
563  ! solutions are to be obtained (Note that
564  ! 1st point is the initial value
565  real*8 h ! input. step size
566  integer ny ! input. see below. ny >= n.
567  real*8 y(ny,m) ! output. solutions at m-points.
568  real*8 wk(n,4) ! output. workin area.
569  integer ierr ! output. condition code.
570  ! 0 ,no error was detected.
571  ! 1 n<1, n> ny, m<2
572 ! the structure of func is the same as for krungekutfds.
573 !-----------------------------------------------------------------------
574 !
575  real*8 zero, half, two, three, const, six
576  integer k, i
577  real*8 t2, t3, x
578 
579  data
580  * zero/0.d0/,
581  * half/0.5d0/, two/2.0d0/, three/3.0d0/, six/6.0d0/
582  data const/ 0.29289321881345247560d0 /
583 
584 
585 
586 
587 
588  if(n .le. 0 .or. n .gt. ny .or. m .le. 1) then
589  ierr=1
590  else
591  ierr=0
592  do k = 1, n
593  wk(k,1)=zero
594  y(k,1)=y0(k)
595  wk(k,2)=y0(k)
596  enddo
597  do i = 1, m-1
598  x=x0+(i-1)*h
599  call func(x, wk(1,2), wk(1,4))
600  do k=1,n
601  t2=h*wk(k, 4)
602  t3=half*t2-wk(k,1)
603  wk(k,3)=wk(k,2)+t3
604  t3=wk(k,3)-wk(k,2)
605  wk(k,1)=wk(k,1)+three*t3-half*t2
606  enddo
607 
608  call func(x+half*h, wk(1,3), wk(1,4))
609  do k=1,n
610  t2=h*wk(k,4)
611  t3=const*(t2-wk(k,1))
612  wk(k,2)=wk(k,3)+t3
613  t3=wk(k,2)-wk(k,3)
614  wk(k,1)=wk(k,1)+three*t3-const*t2
615  enddo
616 
617  call func(x+half*h, wk(1,2), wk(1,4))
618 
619  do k = 1, n
620  t2=h*wk(k,4)
621  t3=(two-const)*(t2-wk(k,1))
622  wk(k,3)=wk(k,2)+t3
623  t3=wk(k,3)-wk(k,2)
624  wk(k,1)=wk(k,1)+three*t3-(two-const)*t2
625  enddo
626 
627  call func(x+h, wk(1,3), wk(1,4))
628  do k=1,n
629  t2=h*wk(k,4)
630  t3=(t2-two*wk(k,1))/six
631  wk(k,2)=wk(k,3)+t3
632  y(k,i+1)=wk(k,2)
633  t3=wk(k,2)-wk(k,3)
634  wk(k,1)=wk(k,1)+three*t3-half*t2
635  enddo
636  enddo
637  endif
i
nodes i
Definition: cstdatmos0.fNew2.h:93
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
d0
block data cblkEvhnp ! currently usable models data RegMdls ad *special data *Cekaon d0
Definition: cblkEvhnp.h:5
y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35
n
integer n
Definition: Zcinippxc.h:1
softenpik::half
integer, parameter half
Definition: csoftenPiK.f:108
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
Here is the caller graph for this function: