V3FIT
spline_akima_int.f
1  SUBROUTINE spline_akima_int(x,y,xx,yy,npts,iflag)
2 ! returns the integral from xx(1) to x of the function
3 ! given by the Akima spline.
4 ! 2012-07-27: Typo correction, J Geiger
5  USE stel_kinds
6  implicit none
7 !-----------------------------------------------
8 ! D u m m y A r g u m e n t s
9 !-----------------------------------------------
10  integer, intent(in) :: npts ! number of active points
11  integer, intent(inout) :: iflag
12  real(rprec), intent(in) :: x
13  real(rprec), intent(out) :: y
14  real(rprec), dimension(npts), intent(in) :: xx, yy
15 !-----------------------------------------------
16 ! L o c a l V a r i a b l e s
17 !-----------------------------------------------
18  real(rprec), dimension(-1:size(xx)+2) :: a, b, c, d
19  real(rprec), dimension(-1:size(xx)+2) :: xloc, yloc,dxloc
20  real(rprec), dimension(-1:size(xx)+2) :: m, t, dm, p, q
21  real(rprec), dimension(-1:size(xx)+2) :: int_f
22  integer :: i, ix, iv, ivm1
23  real(rprec) :: cl,bl,cr,br,dx
24  real(rprec), parameter :: onethird = 1._dp/3._dp
25 !-----------------------------------------------
26 
27  iflag = 0 !initialization
28  ix = size(xx)
29  if(npts > ix)
30  > stop 'spline_akima_int: more active points
31  > requested than available'
32  if(ix /= size(yy)) stop 'size mismatch of xx and yy!'
33  iv=npts
34  ivm1 = iv-1
35 ! initialize local variables
36  a = 0._dp ; b = 0._dp ; c = 0._dp ; d = 0._dp
37  xloc = 0._dp ; yloc = 0._dp
38  m = 0._dp ; t = 0._dp ; dm = 0._dp
39  p = 0._dp ; q = 0._dp
40 
41  xloc(1:iv)=xx
42  xloc(-1)= 2*xloc(1)-xloc(3)
43  xloc( 0)= xloc(1)+xloc(2)-xloc(3)
44  xloc(iv+2)= 2*xloc(iv)-xloc(iv-2)
45  xloc(iv+1)= xloc(iv)+xloc(iv-1)-xloc(iv-2)
46  dxloc(-1:iv+1)=xloc(0:iv+2)-xloc(-1:iv+1)
47  yloc(1:iv)=yy
48 ! calculate linear derivatives as far as existent
49  m(1:iv-1) = (yloc(2:iv)-yloc(1:iv-1))/
50  > dxloc(1:iv-1)
51 ! values for i=0, -1 and iv, iv+1 by quadratic extrapolation:
52  cl = (m(2)-m(1))/(xloc(3)-xloc(1))
53  bl = m(1) - cl*(xloc(2)-xloc(1))
54  cr = (m(iv-2)-m(iv-1))/(xloc(iv)-xloc(iv-2))
55  br = m(iv-2) - cr*(xloc(iv-1)-xloc(iv-2))
56  yloc( 0)=yloc(1)+bl*(xloc( 0)-xloc(1))+
57  > cl*(xloc( 0)-xloc(1))**2
58  yloc(-1)=yloc(1)+bl*(xloc(-1)-xloc(1))+
59  > cl*(xloc(-1)-xloc(1))**2
60  yloc(iv+1)=yloc(iv)+br*(xloc(iv+1)-xloc(iv))+
61  > cr*(xloc(iv+1)-xloc(iv))**2
62  yloc(iv+2)=yloc(iv)+br*(xloc(iv+2)-xloc(iv))+
63  > cr*(xloc(iv+2)-xloc(iv))**2
64 ! rest of linear derivatives
65  m(-1) = (yloc(0)-yloc(-1))/(xloc(0)-xloc(-1))
66  m( 0) = (yloc(1)-yloc( 0))/(xloc(1)-xloc( 0))
67  m(iv ) = (yloc(iv+1)-yloc(iv ))/(xloc(iv+1)-xloc(iv ))
68  m(iv+1) = (yloc(iv+2)-yloc(iv+1))/(xloc(iv+2)-xloc(iv+1))
69 ! calculate weights for derivatives
70  dm(-1:iv)= abs(m(0:iv+1)-m(-1:iv))
71  where (dm /= 0._dp) !exclude division by zero
72  p(1:iv) = dm(1:iv)/(dm(1:iv)+dm(-1:iv-2))
73  end where
74  where (dm /= 0._dp) !exclude division by zero
75  q(1:iv) = dm(-1:iv-2)/(dm(1:iv)+dm(-1:iv-2))
76  end where
77  t(1:iv) = p(1:iv)*m(0:iv-1)+q(1:iv)*m(1:iv)
78  where ( p(1:iv)+q(1:iv) < tiny(1._dp)) ! in case of two zeros give equal weight
79  t(1:iv) = 0.5_dp*m(0:iv-1)+0.5_dp*m(1:iv)
80  end where
81 ! fix coefficients
82  a = yloc
83  b = t
84  c(1:iv-1) = (3*m(1:iv-1)-t(2:iv)-2*t(1:iv-1))/
85  > (dxloc(1:iv-1))
86  d(1:iv-1) = (t(2:iv)+t(1:iv-1)-2*m(1:iv-1))/
87  > (dxloc(1:iv-1))**2
88 
89 ! calculation
90  int_f(1:iv-1) = a(1:iv-1)*dxloc(1:iv-1)+
91  > 0.5_dp*b(1:iv-1)*dxloc(1:iv-1)**2 +
92  > onethird*c(1:iv-1)*dxloc(1:iv-1)**3 +
93  > 0.25_dp*d(1:iv-1)*dxloc(1:iv-1)**4
94  if(x >= xloc(iv)) then
95  y = sum(int_f(1:iv-1))
96  iflag = 0
97  if(x>xloc(iv)) iflag=-1
98  return
99  endif
100  if(x <= xloc(1)) then
101  y = 0._dp
102  iflag = 0
103  if(x<xloc(1)) iflag=-1
104  return
105  endif
106  y=0._dp
107  do i=1,iv-1
108  if(x >= xloc(i+1))then
109  y=y + int_f(i)
110  else
111  dx = x-xloc(i)
112  y=y+dx*(a(i)+dx*(.5_dp*b(i)+dx*(
113  > onethird*c(i)+.25_dp*d(i)*dx)))
114  return
115  endif
116  enddo
117 
118  END SUBROUTINE spline_akima_int
119