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

Go to the source code of this file.

Functions/Subroutines

subroutine kpolintplogxyfe (xa, xstep, ya, ystep, nt, m, logxy, x, y, error)
 
subroutine kpolintplogfe (xa, xstep, ya, ystep, nt, m, x, y, error)
 
subroutine kpolintpfe (xa, xstep, ya, ystep, nt, m, x, y, error)
 
subroutine kpolintp (xa, xstep, ya, ystep, n, x, y, error)
 
subroutine kpolintpeqs (x0, dx, ya, ystep, n, x, y, error)
 

Function/Subroutine Documentation

◆ kpolintp()

subroutine kpolintp ( real*8, dimension(xstep, n)  xa,
integer  xstep,
real*8, dimension(ystep, n)  ya,
integer  ystep,
integer  n,
real*8  x,
real*8  y,
real*8  error 
)

Definition at line 163 of file kpolintp.f.

References parameter().

Referenced by cl2tintp(), kpolintp2(), kpolintpfe(), kpolintplogfe(), and kpolintplogxyfe().

163 !
164 ! integer n. input. number of points.
165 ! real*8 xa(xstep, n). input.
166 ! real*8 ya(ystep, n). input. function values at xa.
167 ! real*8 x. input.
168 ! real*8 y. output. interepolated functon value at x.
169 ! real*8 error. output. estiamted rough error.
170 !
171  implicit none
172  integer n, xstep, ystep
173  real*8 xa(xstep, n), ya(ystep, n), x, y, error
174 
175  integer i, maxm, j
176  parameter(maxm = 10)
177  real*8 c(maxm), d(maxm), diff, difft
178 
179  integer ns, m
180  real*8 h0, hp, w, den
181 
182  if(n .gt. maxm) then
183  write(*, *) ' kpolintp: use lesser number of points'
184  stop
185  endif
186  ns = 1
187  diff = abs(x - xa(1, 1))
188  do i = 1, n
189  difft= abs(x - xa(1, i))
190  if(difft .le. diff) then
191  ns = i
192  diff = difft
193  endif
194  c(i) = ya(1, i)
195  d(i) = ya(1, i)
196  enddo
197  y = ya(1, ns)
198  ns = ns-1
199  do m = 1, n-1
200  do i=1, n-m
201  h0 = xa(1,i) -x
202  hp = xa(1, i+m) - x
203  w = c(i+1) - d(i)
204  den = h0- hp
205  if(den .ne. 0.) then
206  den = w/den
207  d(i) = hp*den
208  c(i) = h0*den
209  else
210  write(0,*) ' error in kpolintp; den= ',den
211  write(0,*) 'x=',x
212  write(0,'(10G12.4)' ) ' xa=', (xa(1,j), j=1, n)
213  write(0,'(10G12.4)' ) ' ya=', (ya(1,j), j=1, n)
214  stop
215  endif
216  enddo
217  if(2*ns .le. n-m) then
218  error = c(ns+1)
219  else
220  error = d(ns)
221  ns = ns-1
222  endif
223  y = y + error
224  enddo
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
h0
! to be included just before the execution code ! density as a function of height real fd0 real h0
Definition: Zstdatmosf.h:3
i
nodes i
Definition: cstdatmos0.fNew2.h:93
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
nrl_atmos::den
real(4), dimension(:), allocatable, save den
Definition: cNRLAtmos.f:28
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
n
integer n
Definition: Zcinippxc.h:1
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
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ kpolintpeqs()

subroutine kpolintpeqs ( real*8  x0,
real*8  dx,
real*8, dimension(ystep, n)  ya,
integer  ystep,
integer  n,
real*8  x,
real*8  y,
real*8  error 
)

Definition at line 228 of file kpolintp.f.

References parameter().

Referenced by kpolintp2().

228 !
229 ! integer n. input. number of points.
230 ! real*8 x0. input. x0, x0+dx, x0+2dx, ...x0+(n-1)dx
231 ! are given data points
232 ! real*8 ya(ystep,n) input. function values at x0,..
233 ! real*8 x. input.
234 ! real*8 y. output. interepolated functon value at x.
235 ! real*8 error. output. estiamted rough error.
236 !
237  implicit none
238  integer n, ystep
239  real*8 x0, ya(ystep, n), x, y, error, dx
240 
241  integer i, maxm
242  parameter(maxm = 10)
243  real*8 c(maxm), d(maxm), diff, difft
244 
245  integer ns, m
246  real*8 h0, hp, w, den
247  integer p, q
248  real*8 xa
249  xa(p, q) = (q-1)*dx + x0
250 
251  if(n .gt. maxm) then
252  write(*, *)
253  * ' kpolintpeqs: use lesser number of points'
254  stop
255  endif
256  ns = 1
257  diff = abs(x - xa(1, 1))
258  do i = 1, n
259  difft= abs(x - xa(1, i))
260  if(difft .le. diff) then
261  ns = i
262  diff = difft
263  endif
264  c(i) = ya(1, i)
265  d(i) = ya(1, i)
266  enddo
267  y = ya(1, ns)
268  ns = ns-1
269  do m = 1, n-1
270  do i=1, n-m
271  h0 = xa(1,i) -x
272  hp = xa(1, i+m) - x
273  w = c(i+1) - d(i)
274  den = h0- hp
275  if(den .eq. 0.) then
276  write(0,*) ' error in kpolintpeqs'
277  stop
278  endif
279  den = w/den
280  d(i) = hp*den
281  c(i) = h0*den
282  enddo
283  if(2*ns .le. n-m) then
284  error = c(ns+1)
285  else
286  error = d(ns)
287  ns = ns-1
288  endif
289  y = y + error
290  enddo
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
h0
! to be included just before the execution code ! density as a function of height real fd0 real h0
Definition: Zstdatmosf.h:3
i
nodes i
Definition: cstdatmos0.fNew2.h:93
dx
integer npitbl real *nx dx real dx
Definition: Zcinippxc.h:10
p
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
nrl_atmos::den
real(4), dimension(:), allocatable, save den
Definition: cNRLAtmos.f:28
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
n
integer n
Definition: Zcinippxc.h:1
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
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function:

◆ kpolintpfe()

subroutine kpolintpfe ( real*8, dimension(xstep, nt)  xa,
integer  xstep,
real*8, dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
real*8  x,
real*8  y,
real*8  error 
)

Definition at line 134 of file kpolintp.f.

References kdwhereis(), and kpolintp().

Referenced by ckmpelaxs(), ckmptotxs(), ckpnelaxs(), ckppelaxs(), cnpelaxs(), cnptotxs(), cpbarpelaxs(), cpimpelaxs(), cpimptotxs(), cpippelaxs(), cpipptotxs(), kpolintplogfe(), kpolintplogxyfe(), and proctimegettf().

134  implicit none
135 ! This is a front end for kpolintp for which we must give
136 ! some few to several points around x. This manages such
137 ! business automatically.
138 !
139  integer xstep ! input. see below
140  integer nt ! input. total number of points
141  integer m ! input. the number of points to be used
142  ! for interpolation. must be <=10.
143  real*8 xa(xstep, nt) ! input. values of x-coordinate at xa(1, i)
144  ! (i=1, nt) are valid x data.
145  integer ystep ! input. see below
146  real*8 ya(ystep, nt) ! input. values of y-coordinate at ya(1, i)
147  ! (i=1, nt) are valid y data.
148  real*8 x ! input. x-value where an interpolated y
149  ! value is wanted
150  real*8 y ! output. see above
151  real*8 error ! output. estimated error
152 
153 
154  integer loc, k
155 ! find location of x in xa
156  call kdwhereis(x, nt, xa, xstep, loc)
157  k = min(max(loc - (m-1)/2,1), nt+1-m) ! max of m points from k
158  call kpolintp(xa(1, k), xstep, ya(1, k), ystep, m, x, y, error)
kpolintp
subroutine kpolintp(xa, xstep, ya, ystep, n, x, y, error)
Definition: kpolintp.f:163
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
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
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:

◆ kpolintplogfe()

subroutine kpolintplogfe ( real*8, dimension(xstep, nt)  xa,
integer  xstep,
real*8, dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
real*8  x,
real*8  y,
real*8  error 
)

Definition at line 91 of file kpolintp.f.

References kdwhereis(), kpolintp(), and kpolintpfe().

91  implicit none
92 ! This is a front end for kpolintp for which we must give
93 ! some few to several points around x. This manages such
94 ! business automatically. This version takes log of y before
95 ! kpolintp is called.
96 !
97  integer xstep ! input. see below
98  integer nt ! input. total number of points
99  integer m ! input. the number of points to be used
100  ! for interpolation. must be <=10.
101  real*8 xa(xstep, nt) ! input. values of x-coordinate at xa(1, i)
102  ! (i=1, nt) are valid x data.
103  integer ystep ! input. see below
104  real*8 ya(ystep, nt) ! input. values of y-coordinate at ya(1, i)
105  ! (i=1, nt) are valid y data.
106  real*8 x ! input. x-value where an interpolated y
107  ! value is wanted
108  real*8 y ! output. see above
109  real*8 error ! output. estimated error
110 ! -----------------------------------------------
111 
112  real*8 logy(10) ! working array.
113  real*8 yy
114 
115  integer loc, k, i
116 ! find location of x in xa
117  call kdwhereis(x, nt, xa, xstep, loc)
118  k = min(max(loc - (m-1)/2,1), nt+1-m) ! max of m points from k
119 
120  do i = k, m+k-1
121  if(ya(1,i) .gt. 0.) then
122  logy(i-k+1) = log(ya(1,i))
123  else
124 ! if some of y is <= 0; we don't use log
125  call kpolintpfe(xa, xstep, ya, ystep, nt, m,
126  * x, y, error)
127  return ! **************
128  endif
129  enddo
130  call kpolintp(xa(1, k), xstep, logy, 1, m, x, yy, error)
131  y = exp(yy)
i
nodes i
Definition: cstdatmos0.fNew2.h:93
kpolintp
subroutine kpolintp(xa, xstep, ya, ystep, n, x, y, error)
Definition: kpolintp.f:163
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
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
kpolintpfe
subroutine kpolintpfe(xa, xstep, ya, ystep, nt, m, x, y, error)
Definition: kpolintp.f:134
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:

◆ kpolintplogxyfe()

subroutine kpolintplogxyfe ( real*8, dimension(xstep, nt)  xa,
integer  xstep,
real*8, dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
integer  logxy,
real*8  x,
real*8  y,
real*8  error 
)

Definition at line 22 of file kpolintp.f.

References kdwhereis(), kpolintp(), and kpolintpfe().

Referenced by cgpxs1(), ckmpelaxs(), ckmptotxs(), ckpnelaxs(), ckpntotxs(), ckppelaxs(), ckpptotxs(), cnpelaxs(), cpbarpelaxs(), cpbarptotxs(), cpimpelaxs(), cpimptotxs(), cpippelaxs(), cpipptotxs(), cppelaxs(), and cpptotxs().

22  implicit none
23 ! This is a front end for kpolintp for which we must give
24 ! some few to several points around x. This manages such
25 ! business automatically. This version takes log of x and or y before
26 ! kpolintp is called.
27 !
28  integer xstep ! input. see below
29  integer nt ! input. total number of points
30  integer m ! input. the number of points to be used
31  ! for interpolation. must be <=10.
32  real*8 xa(xstep, nt) ! input. values of x-coordinate at xa(1, i)
33  ! (i=1, nt) are valid x data.
34  integer ystep ! input. see below
35  real*8 ya(ystep, nt) ! input. values of y-coordinate at ya(1, i)
36  ! (i=1, nt) are valid y data.
37  integer logxy ! input. 1-->log(x) 2-->log(y) 3-->log(x),log(y)
38  real*8 x ! input. x-value where an interpolated y
39  ! value is wanted
40  real*8 y ! output. see above
41 
42  real*8 error ! output. estimated error
43 ! -----------------------------------------------
44 
45  real*8 logx(10), logy(10) ! working array.
46  real*8 xx, yy
47 
48  integer loc, k, i
49  logical kbitest
50 ! find location of x in xa
51  call kdwhereis(x, nt, xa, xstep, loc)
52  k = min(max(loc - (m-1)/2,1), nt+1-m) ! max of m points from k
53 
54  do i = k, m+k-1
55  if( kbitest(logxy, 1)) then
56  if( xa(1, i) .gt. 0.) then
57  logx(i-k+1) = log(xa(1,i))
58  xx = log(x)
59  else
60 ! if some of x is <= 0; we don't use log
61  call kpolintpfe(xa, xstep, ya, ystep, nt, m,
62  * x, y, error)
63  return ! **************
64  endif
65  else
66  logx(i-k+1) = xa(1,i)
67  xx = x
68  endif
69  if( kbitest(logxy, 2)) then
70  if( ya(1, i) .gt. 0.) then
71  logy(i-k+1) = log(ya(1,i))
72  else
73 ! if some of y is <= 0; we don't use log
74  call kpolintpfe(xa, xstep, ya, ystep, nt, m,
75  * x, y, error)
76  return ! **************
77  endif
78  else
79  logy(i-k+1) = ya(1,i)
80  endif
81  enddo
82  call kpolintp(logx, 1, logy, 1, m, xx, yy, error)
83  if( kbitest(logxy, 2) ) then
84  y = exp(yy)
85  else
86  y = yy
87  endif
i
nodes i
Definition: cstdatmos0.fNew2.h:93
kpolintp
subroutine kpolintp(xa, xstep, ya, ystep, n, x, y, error)
Definition: kpolintp.f:163
kdwhereis
subroutine kdwhereis(x, in, a, step, loc)
Definition: kdwhereis.f:27
y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
kbitest
logical function kbitest(i, bit)
Definition: kmanbit.f:13
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
kpolintpfe
subroutine kpolintpfe(xa, xstep, ya, ystep, nt, m, x, y, error)
Definition: kpolintp.f:134
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: