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

Go to the source code of this file.

Functions/Subroutines

subroutine kscsplcoef (x, y, n, coef, nc)
 

Function/Subroutine Documentation

◆ kscsplcoef()

subroutine kscsplcoef ( real(4), dimension(n)  x,
real(4), dimension(n)  y,
integer  n,
real(4), dimension(nc, 3)  coef,
integer  nc 
)

Definition at line 2 of file kScsplCoef.f.

References cerrormsg().

Referenced by cnrldatamanip().

2 ! compute the coefficients of the cubic spline
3  implicit none
4  integer n ! input. size of x, y
5  real(4):: x(n) ! input.
6  real(4):: y(n) ! input.
7  integer nc ! input. size of 1st dim. of coef. >= n-1
8  real(4):: coef(nc, 3) ! output. obtained coeff.
9 
10  character*100 msg
11 
12  real(4):: ec(4), d(2)
13 
14 
15  integer i, ii, ks, ke, ip, i1, i2, k, j
16  real(4):: a2, a1, x1, x2, x3, h, hy, dy1, dy2, y1, y2, hh, h1, h2
17  real(4):: pinv
18 !
19  if (n .lt. 2 .or. nc .lt. n-1 ) then
20  write(0,*) ' n=', n, ' nc=',nc
21  write(0,*) .lt..or..lt." if (n 2 nc n-1 ) wrong"
22  call cerrormsg('input is wrong for kScslpCoef',0)
23  endif
24 
25  do i = 1, n-1
26  if ( x(i) .eq. x(i+1)) then
27  call cerrormsg(
28  * 'x must be ascending: kScsplCoef',1)
29  write(msg,*) ' i=',i, ' x(i)=',x(i), ' x(i+1)=',x(i+1)
30  call cerrormsg(msg, 0)
31  endif
32  enddo
33 
34 
35  ii = 2
36  ks = 1
37  ke = min(4, n)
38  ip = 1
39  do i = 1, 2
40  i1 = 2 * i - 1
41  i2 = 2 * i
42  if (i .ne. 1) then
43  ks = max(1, n-3)
44  ke = n
45  ip = n
46  endif
47  a2 = 0.
48  d(i) = 0.
49  do k = ks, ke
50  if (ip .ne. k) then
51  a1 = 1.0
52  do j = ks, ke
53  if (j .ne. ip .and. j .ne. k) then
54  x1 = x(ip) - x(j)
55  x2 = x(k) - x(j)
56  a1 = a1 * x1 / x2
57  endif
58  enddo
59  x3 = x(k) - x(ip)
60  d(i) = d(i) + a1 * y(k) / x3
61  a2 = a2 - 1.0 / x3
62  endif
63  enddo
64  d(i) = d(i) + y(ip) * a2
65  if (i .eq. 2) ii = n
66  h = x(ii) - x(ii-1)
67  ec(i1) = 1.0
68  hy = y(ii) - y(ii-1)
69  ec(i2) = 6.0 * (hy / h - d(i)) / h
70  if (i .eq. 2) ec(i2) = - ec(i2)
71  enddo
72 
73  if (n .ne. 2) then
74  h1 = x(2) - x(1)
75  y1 = y(2) - y(1)
76  do i = 2, n-1
77  h2 = x(i+1) - x(i)
78  y2 = y(i+1) - y(i)
79  hh = h1 + h2
80  coef(i,1) = h2 / hh
81  coef(i,2) = 1.0 - coef(i,1)
82  coef(i,3) = 6.0 * (y2 / h2 - y1 / h1) / hh
83  h1 = h2
84  y1 = y2
85  enddo
86  endif
87 
88  coef(1,1) = - ec(1) * 0.5
89  coef(1,2) = ec(2) * 0.5
90  if (n .ne. 2) then
91  do k = 2, n-1
92  pinv = 2.0 + coef(k,2) * coef(k-1,1)
93  coef(k,1) = - coef(k,1) / pinv
94  coef(k,2) =
95  * (coef(k,3) - coef(k,2) * coef(k-1,2)) / pinv
96  enddo
97  endif
98  dy1 = (ec(4) - ec(3) * coef(n-1,2)) /
99  * (2.0 + ec(3) * coef(n-1,1))
100  do i = 1, n-1
101  k = n - i
102  dy2 = coef(k,1) * dy1 + coef(k,2)
103  h = x(k+1) - x(k)
104  coef(k,3) = (dy1 - dy2) / (6.0 * h)
105  coef(k,2) = 0.5 * dy2
106  coef(k,1) = (y(k+1) - y(k)) / h -
107  * (coef(k,2) + coef(k,3) * h) * h
108  dy1 = dy2
109  enddo
cerrormsg
subroutine cerrormsg(msg, needrtn)
Definition: cerrorMsg.f:4
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
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
ke
integer maxn LabEquivE real * ke(maxn) integer indx(maxn) integer nevent integer outzero
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
nc
real cut integer nc
Definition: Zprivate.h:1
x2
block data include Zlatfit h c fitting region data x2(1)/0.5/data x1(2)/0.3/
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
Here is the call graph for this function:
Here is the caller graph for this function: