VMEC 8.52
3D Equilibrium solver with nested flux surfaces.
Loading...
Searching...
No Matches
spline_cubic_int.f
Go to the documentation of this file.
1
2 subroutine spline_cubic_int(x,y,xx,yy,n,iflag)
3 USE stel_kinds
4 implicit none
5!----------------------------------------------------------------------
6! iflag = 0 normal return
7! iflag =-1 x-request outside of bounds
8! iflag =-2 xx arrays with two equal entries or not in correct order
9!----------------------------------------------------------------------
10!-----------------------------------------------
11! D u m m y A r g u m e n t s
12!-----------------------------------------------
13 integer, intent(in) :: n
14 integer, intent(inout) :: iflag
15 real(rprec), intent(in) :: x
16 real(rprec), dimension(n), intent(in) :: xx,yy
17 real(rprec), intent(out) :: y
18 real(rprec), dimension(n) :: y2, dxx
19!-----------------------------------------------
20! L o c a l V a r i a b l e s
21!-----------------------------------------------
22 real(rprec) :: yp1, ypn
23 real(rprec) :: c
24!-----------------------------------------------
25
26 iflag = 0 !initialization
27 if(x < xx(1) .or. x > xx(n)) then
28 iflag = -1
29 y=0.d+0
30 return
31 endif
32 dxx(1:n-1)=xx(2:n)-xx(1:n-1)
33 if(minval(dxx(1:n-1)) <= 0.d+0) then
34 iflag=-2
35 return
36 endif
37
38! fix boundary derivatives by quadratic fit
39! left part
40 c=((yy(3)-yy(1))/(xx(3)-xx(1))-(yy(2)-yy(1))/(xx(2)-xx(1)))
41 > /(xx(3)-xx(2))
42 yp1=(yy(2)-yy(1))/(xx(2)-xx(1))-c*(xx(2)-xx(1))
43! right part
44 c=((yy(n-2)-yy(n))/(xx(n-2)-xx(n))
45 > -(yy(n-1)-yy(n))/(xx(n-1)-xx(n)))
46 > /(xx(n-2)-xx(n-1))
47 ypn=(yy(n-1)-yy(n))/(xx(n-1)-xx(n))-c*(xx(n-1)-xx(n))
48
49 call spline_int(xx,yy,n,yp1,ypn,y2)
50 call splint_int(xx,yy,y2,n,x,y)
51
52 return
53
54 contains
55
56 subroutine spline_int(x,y,n,yp1,ypn,y2)
57! taken from numerical recipes f77 and recoded.
58! Given the arrays x(1:n) and y(1:n) containing the tabulated function
59! with x(1) < x(2) <...< x(n) and given values yp1 and ypn for the first
60! derivative of the interpolating function at points q and n, respectively,
61! this routine returns an array y2(1:n) of length n which contains the
62! second derivatives of the interpolating function at the tabulated points x(i).
63! If yp1 and/or ypn are equatl to 1ed+30 or larger, the routine is signaled
64! to set the correspoinding boundary condition for a natural spline with zero
65! derivative on that boundary.
66! nmax is the largest anticipated value of n.
67 integer, intent(in) :: n
68 integer, parameter :: nmax=500
69 real(rprec), intent(in) :: yp1, ypn
70 real(rprec), dimension(n), intent(in) :: x, y
71 real(rprec), dimension(n), intent(out) :: y2
72 integer :: i,k
73 real(rprec) :: p, qn, sig, un
74 real(rprec), dimension(n) :: u
75
76 if(yp1 > .99d+30) then
77 y2(1)=0.d+0
78 u(1) =0.d+0
79 else
80 y2(1)=-0.5d+0
81 u(1)=(3.d+0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
82 endif
83
84 do i=2,n-1
85 sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
86 p=sig*y2(i-1)+2.d+0
87 y2(i)=(sig-1.d+0)/p
88 u(i)=(6.d+0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
89 > /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
90 enddo
91 if (ypn > .99d+30)then
92 qn=0.d+0
93 un=0.d+0
94 else
95 qn=0.5d+0
96 un=(3.d+0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
97 endif
98 y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d+0)
99 do k=n-1,1,-1
100 y2(k)=y2(k)*y2(k+1)+u(k)
101 enddo
102 return
103 end subroutine spline_int
104
105 subroutine splint_int(xa,ya,y2a,n,x,y)
106! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a function
107! with the xa(i)'s in order), and given the array y2a(1:n), which is the
108! output from spline above, and given a value of x, this routine returns
109! a cubic-spline interpolated value y.
110 implicit none
111 integer, intent(in) :: n
112 real(rprec), intent(in) :: x
113 real(rprec), intent(out) :: y
114 real(rprec), dimension(n), intent(in) :: xa, ya, y2a
115!- local -------------------
116 real(rprec) :: h, c, dx
117 integer :: k, khi, klo
118 real(rprec), dimension(n) :: dxa,dxa3, dya, dy2a, my2a, mya
119
120 klo=1
121 khi=n
122 do
123 if (khi-klo <= 1) exit !inverted num.rec. condition for endless do-loop exit
124 k=(khi+klo)/2
125 if(xa(k) > x) then
126 khi=k
127 else
128 klo=k
129 endif
130 enddo
131
132 dxa(1:n-1)=xa(2:n)-xa(1:n-1)
133 dxa3=dxa**3
134 dya(1:n-1)=ya(2:n)-ya(1:n-1)
135 mya(1:n-1)=ya(2:n)+ya(1:n-1)
136 dy2a(1:n-1)=y2a(2:n)-y2a(1:n-1)
137 my2a(1:n-1)=y2a(2:n)+y2a(1:n-1)
138! integral up to specific interval [klo:khi] in which x is.
139 c = .5d+0*dot_product(dxa(1:klo-1),mya(1:klo-1))
140 > -dot_product(dxa3(1:klo-1),my2a(1:klo-1))/24.0d+0
141 h=xa(khi)-xa(klo)
142 if(h ==0.d+0)
143 > stop "spline_cubic: bad xa input! xa(i) have to be distinct!"
144 dx=x-xa(klo)
145 y = ya(klo)*dx +.5d+0*dya(klo)*dx**2/h +
146 > y2a(klo)*dx**2*(2*dx-3*h)/12.d+0 +
147 > dy2a(klo)*dx**2*(dx**2-2*h**2)/(h*24.d+0)
148 y= y+c
149 return
150 end subroutine splint_int
151
152 end subroutine spline_cubic_int
153
stel_kinds
Definition stel_kinds.f90:2
spline_int
subroutine spline_int(x, y, n, yp1, ypn, y2)
Definition spline_cubic_int.f:57
spline_cubic_int
subroutine spline_cubic_int(x, y, xx, yy, n, iflag)
Definition spline_cubic_int.f:3