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

Go to the source code of this file.

Functions/Subroutines

subroutine kpolintpslogxyfe (xa, xstep, ya, ystep, nt, m, logxy, x, y, error)
 
subroutine kpolintpslogfe (xa, xstep, ya, ystep, nt, m, x, y, error)
 
subroutine kpolintpsfe (xa, xstep, ya, ystep, nt, m, x, y, error)
 
subroutine kpolintps (xa, xstep, ya, ystep, n, x, y, error)
 
subroutine kpolintpseqs (x0, dx, ya, ystep, n, x, y, error)
 

Function/Subroutine Documentation

◆ kpolintps()

subroutine kpolintps ( real(4), dimension(xstep, n)  xa,
integer  xstep,
real(4), dimension(ystep, n)  ya,
integer  ystep,
integer  n,
real(4)  x,
real(4)  y,
real(4)  error 
)

Definition at line 164 of file kpolintpS.f.

References parameter().

Referenced by kpolintp2s(), kpolintpsfe(), kpolintpslogfe(), and kpolintpslogxyfe().

164 !
165 ! integer n. input. number of points.
166 ! real(4):: xa(xstep, n). input.
167 ! real(4):: ya(ystep, n). input. function values at xa.
168 ! real(4):: x. input.
169 ! real(4):: y. output. interepolated functon value at x.
170 ! real(4):: error. output. estiamted rough error.
171 !
172  implicit none
173  integer n, xstep, ystep
174  real(4):: xa(xstep, n), ya(ystep, n), x, y, error
175 
176  integer i, maxm, j
177  parameter(maxm = 10)
178  real(4):: c(maxm), d(maxm), diff, difft
179 
180  integer ns, m
181  real(4):: h0, hp, w, den
182 
183  if(n .gt. maxm) then
184  write(*, *) ' kpolintpS: use lesser number of points'
185  stop
186  endif
187  ns = 1
188  diff = abs(x - xa(1, 1))
189  do i = 1, n
190  difft= abs(x - xa(1, i))
191  if(difft .le. diff) then
192  ns = i
193  diff = difft
194  endif
195  c(i) = ya(1, i)
196  d(i) = ya(1, i)
197  enddo
198  y = ya(1, ns)
199  ns = ns-1
200  do m = 1, n-1
201  do i=1, n-m
202  h0 = xa(1,i) -x
203  hp = xa(1, i+m) - x
204  w = c(i+1) - d(i)
205  den = h0- hp
206  if(den .ne. 0.) then
207  den = w/den
208  d(i) = hp*den
209  c(i) = h0*den
210  else
211  write(0,*) ' error in kpolintpS'
212  write(0,*) 'x=',x
213  write(0,'(10G12.4)' ) ' xa=', (xa(j,1), j=1, n)
214  write(0,'(10G12.4)' ) ' ya=', (ya(j,1), j=1, n)
215  stop
216  endif
217  enddo
218  if(2*ns .le. n-m) then
219  error = c(ns+1)
220  else
221  error = d(ns)
222  ns = ns-1
223  endif
224  y = y + error
225  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:

◆ kpolintpseqs()

subroutine kpolintpseqs ( real(4)  x0,
real(4)  dx,
real(4), dimension(ystep, n)  ya,
integer  ystep,
integer  n,
real(4)  x,
real(4)  y,
real(4)  error 
)

Definition at line 229 of file kpolintpS.f.

References parameter().

Referenced by kpolintp2s().

229 !
230 ! integer n. input. number of points.
231 ! real(4):: x0. input. x0, x0+dx, x0+2dx, ...x0+(n-1)dx
232 ! are given data points
233 ! real(4):: ya(ystep,n) input. function values at x0,..
234 ! real(4):: x. input.
235 ! real(4):: y. output. interepolated functon value at x.
236 ! real(4):: error. output. estiamted rough error.
237 !
238  implicit none
239  integer n, ystep
240  real(4):: x0, ya(ystep, n), x, y, error, dx
241 
242  integer i, maxm
243  parameter(maxm = 10)
244  real(4):: c(maxm), d(maxm), diff, difft
245 
246  integer ns, m
247  real(4):: h0, hp, w, den
248  integer p, q
249  real(4):: xa
250  xa(p, q) = (q-1)*dx + x0
251 
252  if(n .gt. maxm) then
253  write(*, *)
254  * ' kpolintpSeqs: use lesser number of points'
255  stop
256  endif
257  ns = 1
258  diff = abs(x - xa(1, 1))
259  do i = 1, n
260  difft= abs(x - xa(1, i))
261  if(difft .le. diff) then
262  ns = i
263  diff = difft
264  endif
265  c(i) = ya(1, i)
266  d(i) = ya(1, i)
267  enddo
268  y = ya(1, ns)
269  ns = ns-1
270  do m = 1, n-1
271  do i=1, n-m
272  h0 = xa(1,i) -x
273  hp = xa(1, i+m) - x
274  w = c(i+1) - d(i)
275  den = h0- hp
276  if(den .eq. 0.) then
277  write(0,*) ' error in kpolintpSeqs'
278  stop
279  endif
280  den = w/den
281  d(i) = hp*den
282  c(i) = h0*den
283  enddo
284  if(2*ns .le. n-m) then
285  error = c(ns+1)
286  else
287  error = d(ns)
288  ns = ns-1
289  endif
290  y = y + error
291  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:

◆ kpolintpsfe()

subroutine kpolintpsfe ( real(4), dimension(xstep, nt)  xa,
integer  xstep,
real(4), dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
real(4)  x,
real(4)  y,
real(4)  error 
)

Definition at line 135 of file kpolintpS.f.

References kpolintps(), and kwhereis().

Referenced by kpolintpslogfe(), and kpolintpslogxyfe().

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

◆ kpolintpslogfe()

subroutine kpolintpslogfe ( real(4), dimension(xstep, nt)  xa,
integer  xstep,
real(4), dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
real(4)  x,
real(4)  y,
real(4)  error 
)

Definition at line 92 of file kpolintpS.f.

References kpolintps(), kpolintpsfe(), and kwhereis().

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

◆ kpolintpslogxyfe()

subroutine kpolintpslogxyfe ( real(4), dimension(xstep, nt)  xa,
integer  xstep,
real(4), dimension(ystep, nt)  ya,
integer  ystep,
integer  nt,
integer  m,
integer  logxy,
real(4)  x,
real(4)  y,
real(4)  error 
)

Definition at line 23 of file kpolintpS.f.

References kpolintps(), kpolintpsfe(), and kwhereis().

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