VMEC 8.52
3D Equilibrium solver with nested flux surfaces.
Loading...
Searching...
No Matches
explore_spline_akima.f90
Go to the documentation of this file.
1program hello
2 implicit none
3 INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
4 integer, parameter :: dp = rprec
5
6 integer, parameter :: l = 91
7
8 real(rprec), dimension(1:l) :: am_aux_s = (/ 0.0000000e+00_dp, &
9 1.8181818e-02_dp, 3.6363636e-02_dp, 5.4545455e-02_dp, 7.2727273e-02_dp, 9.0909091e-02_dp, &
10 1.0909091e-01_dp, 1.2727273e-01_dp, 1.4545455e-01_dp, 1.6363636e-01_dp, 1.8181818e-01_dp, &
11 2.0000000e-01_dp, 2.1818182e-01_dp, 2.3636364e-01_dp, 2.5454545e-01_dp, 2.7272727e-01_dp, &
12 2.9090909e-01_dp, 3.0909091e-01_dp, 3.2727273e-01_dp, 3.4545455e-01_dp, 3.6363636e-01_dp, &
13 3.8181818e-01_dp, 4.0000000e-01_dp, 4.1818182e-01_dp, 4.3636364e-01_dp, 4.5454545e-01_dp, &
14 4.7272727e-01_dp, 4.9090909e-01_dp, 5.0909091e-01_dp, 5.2727273e-01_dp, 5.4545455e-01_dp, &
15 5.6363636e-01_dp, 5.8181818e-01_dp, 6.0000000e-01_dp, 6.1818182e-01_dp, 6.3636364e-01_dp, &
16 6.5454545e-01_dp, 6.7272727e-01_dp, 6.9090909e-01_dp, 7.0909091e-01_dp, 7.2727273e-01_dp, &
17 7.4545455e-01_dp, 7.6363636e-01_dp, 7.8181818e-01_dp, 8.0000000e-01_dp, 8.0444444e-01_dp, &
18 8.0888889e-01_dp, 8.1333333e-01_dp, 8.1777778e-01_dp, 8.2222222e-01_dp, 8.2666667e-01_dp, &
19 8.3111111e-01_dp, 8.3555556e-01_dp, 8.4000000e-01_dp, 8.4444444e-01_dp, 8.4888889e-01_dp, &
20 8.5333333e-01_dp, 8.5777778e-01_dp, 8.6222222e-01_dp, 8.6666667e-01_dp, 8.7111111e-01_dp, &
21 8.7555556e-01_dp, 8.8000000e-01_dp, 8.8444444e-01_dp, 8.8888889e-01_dp, 8.9333333e-01_dp, &
22 8.9777778e-01_dp, 9.0222222e-01_dp, 9.0666667e-01_dp, 9.1111111e-01_dp, 9.1555556e-01_dp, &
23 9.2000000e-01_dp, 9.2444444e-01_dp, 9.2888889e-01_dp, 9.3333333e-01_dp, 9.3777778e-01_dp, &
24 9.4222222e-01_dp, 9.4666667e-01_dp, 9.5111111e-01_dp, 9.5555556e-01_dp, 9.6000000e-01_dp, &
25 9.6444444e-01_dp, 9.6888889e-01_dp, 9.7333333e-01_dp, 9.7777778e-01_dp, 9.8222222e-01_dp, &
26 9.8666667e-01_dp, 9.9111111e-01_dp, 9.9555556e-01_dp, 1.0000000e+00_dp, 0.0_dp /)
27
28 real(rprec), dimension(1:l) :: am_aux_f = (/ 1.1343073e+05_dp, &
29 1.0511264e+05_dp, 9.7719901e+04_dp, 9.1142143e+04_dp, 8.5238372e+04_dp, 7.9893902e+04_dp, &
30 7.5021684e+04_dp, 7.0555497e+04_dp, 6.6444168e+04_dp, 6.2647239e+04_dp, 5.9132114e+04_dp, &
31 5.5871944e+04_dp, 5.2844250e+04_dp, 5.0029862e+04_dp, 4.7412214e+04_dp, 4.4976739e+04_dp, &
32 4.2710495e+04_dp, 4.0601792e+04_dp, 3.8639948e+04_dp, 3.6815020e+04_dp, 3.5117615e+04_dp, &
33 3.3538679e+04_dp, 3.2069304e+04_dp, 3.0700518e+04_dp, 2.9423071e+04_dp, 2.8227198e+04_dp, &
34 2.7102369e+04_dp, 2.6037002e+04_dp, 2.5018180e+04_dp, 2.4031951e+04_dp, 2.3067082e+04_dp, &
35 2.2117542e+04_dp, 2.1181454e+04_dp, 2.0260068e+04_dp, 1.9356991e+04_dp, 1.8477613e+04_dp, &
36 1.7628651e+04_dp, 1.6817841e+04_dp, 1.6053656e+04_dp, 1.5345090e+04_dp, 1.4699130e+04_dp, &
37 1.4113009e+04_dp, 1.3572603e+04_dp, 1.3055378e+04_dp, 1.2532538e+04_dp, 1.2400160e+04_dp, &
38 1.2264938e+04_dp, 1.2126344e+04_dp, 1.1983859e+04_dp, 1.1836908e+04_dp, 1.1684914e+04_dp, &
39 1.1527343e+04_dp, 1.1363692e+04_dp, 1.1193491e+04_dp, 1.1016295e+04_dp, 1.0831470e+04_dp, &
40 1.0638256e+04_dp, 1.0435957e+04_dp, 1.0224325e+04_dp, 1.0003778e+04_dp, 9.7747634e+03_dp, &
41 9.5363620e+03_dp, 9.2816431e+03_dp, 9.0022662e+03_dp, 8.6903284e+03_dp, 8.3409886e+03_dp, &
42 7.9527438e+03_dp, 7.5237757e+03_dp, 7.0538828e+03_dp, 6.5554498e+03_dp, 6.0461036e+03_dp, &
43 5.5435531e+03_dp, 5.0630772e+03_dp, 4.6140337e+03_dp, 4.2046228e+03_dp, 3.8422577e+03_dp, &
44 3.5274335e+03_dp, 3.2543658e+03_dp, 3.0169701e+03_dp, 2.8089885e+03_dp, 2.6245185e+03_dp, &
45 2.4577789e+03_dp, 2.3032608e+03_dp, 2.1572214e+03_dp, 2.0204623e+03_dp, 1.8947411e+03_dp, &
46 1.7810514e+03_dp, 1.6799853e+03_dp, 1.5917883e+03_dp, 1.5163600e+03_dp, 0.0_dp /)
47
48 integer :: i, j, iflag
49 integer, parameter :: ns = 25
50 real(rprec) :: hs, si
51 real(rprec) :: y
52
53 j = minloc(am_aux_s(2:), dim=1)
54 print *, "length : ", j
55
56 hs = 1.0_rprec/(ns-1)
57
58 DO i = 2, ns
59 si = hs * (i-1.5_rprec)
60
61 call spline_akima(si, y, am_aux_s, am_aux_f, j, iflag)
62
63 if (iflag .eq. 0) then
64 print *, "si = ", si, " => y = ", y
65 else
66 print *, "error: iflag = ", iflag
67 end if
68
69 end do
70
71end program hello
72
73SUBROUTINE spline_akima(x, y, xx, yy, npts, iflag)
74 implicit none
75
76 INTEGER, PARAMETER :: rprec = selected_real_kind(12,100)
77 INTEGER, PARAMETER :: dp = rprec
78
79!-----------------------------------------------
80! D u m m y A r g u m e n t s
81!-----------------------------------------------
82 integer, intent(in) :: npts !< number of active points
83 real(rprec), intent(in) :: x !< location at which to evaluate
84 real(rprec), intent(out) :: y !< evaluated function value
85 real(rprec), dimension(npts), intent(in) :: xx !< knots of data to interpolate
86 real(rprec), dimension(npts), intent(in) :: yy !< values of data to interpolate
87 integer, intent(inout) :: iflag !< error flag; 0 means all good
88!-----------------------------------------------
89! L o c a l V a r i a b l e s
90!-----------------------------------------------
91 real(rprec), dimension(-1:size(xx)+2) :: a !< constant coefficients
92 real(rprec), dimension(-1:size(xx)+2) :: b !< 1st-order coefficients
93 real(rprec), dimension(-1:size(xx)+2) :: c !< 2nd-order coefficients
94 real(rprec), dimension(-1:size(xx)+2) :: d !< 3rd-order coefficients
95
96 real(rprec), dimension(-1:size(xx)+2) :: xloc !< copy of xx, extended by 2 points at start and end
97 real(rprec), dimension(-1:size(xx)+2) :: yloc !< copy of yy, quadratically-extrapolated beyond edge
98
99 real(rprec), dimension(-1:size(xx)+2) :: m !< linear derivatives
100 real(rprec), dimension(-1:size(xx)+2) :: t !< intermediate quantity
101 real(rprec), dimension(-1:size(xx)+2) :: dm !< intermediate quantity
102 real(rprec), dimension(-1:size(xx)+2) :: p !< intermediate quantity
103 real(rprec), dimension(-1:size(xx)+2) :: q !< intermediate quantity
104
105 integer :: i !< iteration variable
106 integer :: ix !< size of xx, yy
107 integer :: iv !< number of points to actually use: 1:iv
108 integer :: ivm1 !< iv-1
109
110 real(rprec) :: cl !< intermediate quantity
111 real(rprec) :: bl !< intermediate quantity
112 real(rprec) :: cr !< intermediate quantity
113 real(rprec) :: br !< intermediate quantity
114
115 real(rprec) :: dx !< temporary variable
116
117!-----------------------------------------------
118
119 !WRITE(6,*) x, y
120 !WRITE(6,*)npts
121 !WRITE(6,*) xx
122 !WRITE(6,*) yy
123 !stop 'Test!'
124
125 iflag = 0 !initialization
126
127 ix = size(xx)
128! print *, size(xx)
129! print *, size(yy)
130 if(npts > ix) &
131 stop'spline_akima: more active points requested than available'
132 if(ix /= size(yy)) stop 'size mismatch of xx and yy!'
133
134 iv=npts
135
136 ivm1 = iv-1
137! initialize local variables
138 a = 0._dp ; b = 0._dp ; c = 0._dp ; d = 0._dp
139 xloc = 0._dp ; yloc = 0._dp
140 m = 0._dp ; t = 0._dp ; dm = 0._dp
141 p = 0._dp ; q = 0._dp
142
143 xloc(1:iv)=xx
144 xloc(-1)= 2*xloc(1)-xloc(3)
145 xloc( 0)= xloc(1)+xloc(2)-xloc(3)
146 xloc(iv+2)= 2*xloc(iv)-xloc(iv-2)
147 xloc(iv+1)= xloc(iv)+xloc(iv-1)-xloc(iv-2)
148
149 !print *, xloc(-1), xloc(0), xloc(iv+1), xloc(iv+2)
150
151 yloc(1:iv)=yy
152
153! calculate linear derivatives as far as existent
154 m(-1:iv+1) = (yloc(0:iv+2)-yloc(-1:iv+1))/ &
155 (xloc(0:iv+2)-xloc(-1:iv+1))
156! values for i=0, -1 and iv, iv+1 by quadratic extrapolation:
157 cl = (m(2)-m(1))/(xloc(3)-xloc(1))
158 bl = m(1) - cl*(xloc(2)-xloc(1))
159 cr = (m(iv-2)-m(iv-1))/(xloc(iv)-xloc(iv-2))
160 br = m(iv-2) - cr*(xloc(iv-1)-xloc(iv-2))
161
162 !print *, cl, bl, cr, br
163
164 yloc( 0)=yloc(1)+bl*(xloc( 0)-xloc(1))+ &
165 cl*(xloc( 0)-xloc(1))**2
166 yloc(-1)=yloc(1)+bl*(xloc(-1)-xloc(1))+ &
167 cl*(xloc(-1)-xloc(1))**2
168 yloc(iv+1)=yloc(iv)+br*(xloc(iv+1)-xloc(iv))+ &
169 cl*(xloc(iv+1)-xloc(iv))**2
170 yloc(iv+2)=yloc(iv)+br*(xloc(iv+2)-xloc(iv))+ &
171 cl*(xloc(iv+2)-xloc(iv))**2
172
173 !print *, yloc(-1), yloc(0), yloc(iv+1), yloc(iv+2)
174
175! rest of linear derivatives
176 m(-1) = (yloc(0)-yloc(-1))/(xloc(0)-xloc(-1))
177 m( 0) = (yloc(1)-yloc( 0))/(xloc(1)-xloc( 0))
178 m(iv ) = (yloc(iv+1)-yloc(iv ))/(xloc(iv+1)-xloc(iv ))
179 m(iv+1) = (yloc(iv+2)-yloc(iv+1))/(xloc(iv+2)-xloc(iv+1))
180
181 !print *, m(-1), m(0), m(iv), m(iv+1)
182
183! calculate weights for derivatives
184 dm(-1:iv)= abs(m(0:iv+1)-m(-1:iv))
185 where (dm /= 0._dp) !exclude division by zero
186 p(1:iv) = dm(1:iv)/(dm(1:iv)+dm(-1:iv-2))
187 end where
188 where (dm /= 0._dp) !exclude division by zero
189 q(1:iv) = dm(-1:iv-2)/(dm(1:iv)+dm(-1:iv-2))
190 end where
191
192 print *, " i p(i) q(i)"
193 do i=1, iv
194 print *, i, p(i), q(i)
195 end do
196
197 stop
198
199 t(1:iv) = p(1:iv)*m(0:iv-1)+q(1:iv)*m(1:iv)
200 where ( p(1:iv)+q(1:iv) < tiny(1._dp)) ! in case of two zeros give equal weight
201 t(1:iv) = 0.5_dp*m(0:iv-1)+0.5_dp*m(1:iv)
202 end where
203! fix coefficients
204 a = yloc
205 b = t
206 c(1:iv-1) = (3*m(1:iv-1)-t(2:iv)-2*t(1:iv-1))/ &
207 (xloc(2:iv)-xloc(1:iv-1))
208 d(1:iv-1) = (t(2:iv)+t(1:iv-1)-2*m(1:iv-1))/ &
209 (xloc(2:iv)-xloc(1:iv-1))**2
210
211! calculation
212 if(x<xloc(1) .or. x>xloc(iv)) then
213 y=0.0
214 iflag=-1
215 return
216 endif
217
218 if(x == xloc(iv)) then
219 y = yy(iv)
220 iflag = 0
221 return
222 endif
223
224 do i=1,iv-1
225 if(x >= xloc(i) .and. x < xloc(i+1))then
226 dx=x-xloc(i)
227 y=a(i)+dx*(b(i)+dx*(c(i)+d(i)*dx))
228 iflag = 0
229 return
230 endif
231 enddo
232
233 END SUBROUTINE spline_akima
spline_akima
subroutine spline_akima(x, y, xx, yy, npts, iflag)
Definition explore_spline_akima.f90:74
hello
program hello
Definition explore_spline_akima.f90:1