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