V3FIT
cntourp.f
1  SUBROUTINE cntourp(x,nw,y,nh,cspln,xc,yc,gsq,ipts,dx,dy,ntmx,psivl
2  $ ,nxd)
3 c...........................................................................
4 c 02/04/97 (yr)
5 c add a new output array gsq to cntour and rename it to cntourp
6 c the original cntour is still avaible in the lib
7 c 06/14/96
8 c cntour generates a contour of ordered points (xc(i),yc(i)) for
9 c i = 1,ipts along a contour psivl of a FUNCTION psi(x,y). The function
10 c values psi(x,y) are defined from the spline coefficients in cspln.
11 c The contour psivl must fully encircle a given poINT (xaxd,yaxd).
12 c
13 c The contour is determined by defining a series of rays from
14 c (xaxd,yaxd) INTersecting with the contour.
15 c Each ray is defined by an equation y = a*x + b or x = a*y + b,
16 c depending on whether the angle is in the left or right quadrants
17 c or in the top or bottom quadrants respectively.
18 c For a given ray, the INTersection is first bounded by a coarse sliding
19 c INTerval search along the ray to find two points on the ray bounding
20 c the INTersection. The search is limited to a rectangle specified by
21 c xmin,xmax,ymin,ymax. A sliding INTerval search, rather than a
22 c binary search, is USEd for the coarse grid so that non-monotonic
23 c psi can be handled. Newtons Method is THEN USEd to converge to the
24 c INTersection poINT and the coordinates are stored as (xc(i),yc(i)).
25 c This task is performed by the SUBROUTINE 'rayatt'.
26 c
27 c The ray angular spacing is controlled dynamically and so the total
28 c number of contour points (xc(i),yc(i)) is determined by the routine.
29 c Normally, the increment is defined as the minimum of arcl/rad and
30 c 2pi*dang*(rad0/rad), WHERE dang is an angle fraction of 2pi and
31 c arcl is an arclength spacing between successive points, both input
32 c through common/cntd/, and rad and rad0 are the computed radii at the
33 c respective and first (i.e. theta = 0) contour points respectively.
34 c After each poINT on the contour is found, the relative change in the
35 c poloidal field from the previous poINT is checked to ensure that it
36 c is less than bpers. If the change in poloidal field exceeds bpers,
37 c the angular increment is successively reduced, up to ihlfmx times,
38 c and the poINT is recalculated.
39 c If the number of contour points exceeds the maximum ALLowed DIMENSION
40 c ntmx, the mapping for this contour is restarted with bpers increased
41 c and appropriate warnings are PRINTed. The mapping can restart
42 c up to jcmax tries, after which, the routine aborts.
43 c
44 c Input parameters :
45 c x,y : Vectors containing the rectangular grid poINT values on
46 c which the contoured FUNCTION is defined.
47 c nw,nh : Respective DIMENSIONs of the vectors x and y.
48 c cspln : Spline coefficients for the contoured FUNCTION.
49 c dx,dy : Increment (in metres) for the coarse grid search of
50 c the INTersection between the contour and rays.
51 c dx and dy should normally be chosen so that, over the
52 c distance dx or dy, psi is SINgle valued. This can
53 c usually be ensured by setting dx and dy to the x and y
54 c grid increments.
55 c ntmx : Maximum DIMENSION for the vector containing the contour
56 c points.
57 c psivl : Value of the contour to be mapped.
58 c
59 c Input paramaters passed through common/cntd/ :
60 c xaxd,yaxd : Interior poINT for the contour. The rays INTersecting
61 c the contour emmanate from the poINT (xaxd,yaxd).
62 c xmin,xmax : The search for INTersection points is limited
63 c horizontally to (xmin,xmax).
64 c ymin,ymax : The search for INTersection points is limited vertically
65 c to (ymin,ymax).
66 c dang : The maximum angular increment ALLowed for the generated
67 c contour points is limited to
68 c 2pi*dang*(rad(theta=0)/rad(theta)).
69 c In general, dang should be set according to shape of the
70 c equilibrium. If dang is too small, mANY unnecessary
71 c points will be generated. If dang is too large for
72 c bperr (defined below), the routine wastes computation
73 c time in cutting the step angle DOwn to SIZE.
74 c dang = 0.03 near psilim and dang = 0.10 near psimax
75 c is usually satisfactory.
76 c arcl : The maximum arclength ALLowed between generated contour
77 c points is limited to arcl. For highly elongated plasmas
78 c control with dang is difficult to set for ALL contours.
79 c For such CASEs, arcl should be USEd instead to set the
80 c angular increment. arcl can be set to a large number to
81 c overide this option so that dang is always USEd.
82 c bpers : Relative change in poloidal b field between (xc(i),yc(i))
83 c and (xc(i+1),yc(i+1)). If dang yields a relative change
84 c in poloidal b less than bperr, THEN the the computed
85 c poINT is retained. Otherwise, dang is successively
86 c reduced by stfrac until this condition is met. bpers
87 c is started at bperr (set in chali). After each failure
88 c (i.e. IF ipts gets too large), THEN the calculation is
89 c restarted with larger bpers, up to a maximum of jcmax
90 c tries.
91 c
92 c Output parameters :
93 c xc,yc : Contour points generated at INTersections of the
94 c psivl contour and the rays emmanating from (xaxd,yaxd).
95 c gsq : square of grad(psi)--changed to bpsquared rlm 2/11/97
96 c ipts : Number of contour points generated. ipts must be less
97 c than or equal to ntmx.
98 c
99 c Output parameters passed through common/cntd/ :
100 c xemin,xemax : Minimum and maximum x values of the contour.
101 c yemin,yemax : Minimum and maximum y values of the contour.
102 c xymin,xymax : x values at y=yemin and y=yemax respectively.
103 c yxmin,yxmax : y values at x=xemin and x=xemax respectively.
104 c
105 c jcmax is the maximum number of restarts ALLowed.
106 c ihlfmx is the maximum number of increment halvings ALLowed for
107 c each contour poINT calculation.
108 c rndoff0 is a roundoff error for checking for division by zero.
109 c stfrac is the fractional reduction in step SIZE imposed when the
110 c test for the poloidal field fails.
111 c
112  USE precision
113  IMPLICIT NONE
114  INTEGER nw,nh,nxd,ntmx,ipts
115  REAL(rprec) x(nw),y(nh),cspln(2,nxd,nh),xc(ntmx)
116  $ ,yc(ntmx),gsq(ntmx), dx,dy
117  REAL(rprec) psivl
118  REAL(rprec) pds(6)
119  REAL(rprec) xaxd,yaxd,psiaxd
120  REAL(rprec) xemin,xemax,yemin,yemax,yxmin,yxmax,xymin,xymax
121  REAL(rprec) dang,arcl,bperr,xmin,xmax,ymin,ymax
122  common/cntd/xaxd,yaxd,
123  $ xemin,xemax,yemin,yemax,yxmin,yxmax,xymin,xymax,
124  $ dang,arcl,bperr,xmin,xmax,ymin,ymax
125  INTEGER jcmax,ihlfmx
126  INTEGER i, j
127  REAL(rprec) stfrac,delta
128  parameter(jcmax=5, ihlfmx=4)
129  parameter(stfrac=0.5,delta=0.1)
130 c........................................................................
131 c local variables:
132  REAL(rprec) pi,twopi
133  REAL(rprec) thet,dthet,theta0,dthet0,dthet1,dthet2,xn,yn,rad0,rad
134  REAL(rprec) bp1,bp2,bpers,bpdiff
135  INTEGER IFail
136  INTEGER ntmx1,ier,jcount,ihalf,ipt1,ik
137  REAL(rprec) xcfn,ycfn,deli2,deln2,delfn2,delti2,deltn2
138 c.........................................................................
139 c Initialization
140 c.........................................................................
141  pi = acos(-1._dbl)
142  twopi = 2.00*pi
143  xemin = +1.0e+10
144  xemax = -1.0e+10
145  yemin = +1.0e+10
146  yemax = -1.0e+10
147  xymin = 0.0
148  xymax = 0.0
149  yxmin = 0.0
150  yxmax = 0.0
151  ntmx1 = ntmx - 1
152  dthet0 = dang*twopi
153 
154 ! MRC: Add hack to debug efit gfile issues. This dumps out the flux surface
155 ! contours to a file.
156  OPEN(unit=1, file="dump.txt", action='WRITE')
157  do i = 1, nxd
158  do j = 1, nxd
159  CALL bcsplint(x,y,cspln,nw,nh,nxd,x(i),y(j),pds,1,ier)
160  WRITE (1,*) x(i),y(j),pds(1)
161  enddo
162  enddo
163  CLOSE(unit=1)
164 
165  CALL bcsplint(x,y,cspln,nw,nh,nxd,xaxd,yaxd,pds,1,ier)
166  psiaxd = pds(1)
167  IF(psiaxd .gt. psivl) THEN
168  WRITE(6,*)'errorin cntourp-cnt1'
169  stop
170  ENDIF
171 c.........................................................................
172 c find (x,y) at thet=0 and initialize some quantities
173 c.........................................................................
174  thet=0.
175  CALL rayatt(x,nw,y,nh,cspln,nxd,dx,dy,xaxd,yaxd,psiaxd,
176  $ xmin,xmax,ymin,ymax,psivl,thet,xn,yn,pds,ifail)
177  IF (ifail .gt. 0)stop
178  bp1 = sqrt(pds(2)*pds(2) + pds(3)*pds(3))/xn
179  xemin = min(xn,xemin)
180  xemax = max(xn,xemax)
181  yemin = min(yn,yemin)
182  yemax = max(yn,yemax)
183  IF(xn .eq. xemax) yxmax = yn
184  IF(xn .eq. xemin) yxmin = yn
185  IF(yn .eq. yemax) xymax = xn
186  IF(yn .eq. yemin) xymin = xn
187  xc(1)=xn
188  yc(1)=yn
189  gsq(1)=(pds(2)**2+pds(3)**2)/xn**2
190  rad0=sqrt((xn-xaxd)**2+(yn-yaxd)**2)
191 c...........................................................................
192 c main loop
193 c...........................................................................
194 c Set up the loop to increment bpers up to jcmax tries, in CASE that
195 c bpers is too samll and there are too mANY contour points found.
196  jcount = 0
197  5 jcount = jcount + 1 ! restart poINT
198  bpers = bperr*jcount
199  ipts = 1
200  ihalf = 0
201  dthet = min(dthet0,arcl/rad0)
202  thet = dthet
203  theta0 = thet
204  10 CONTINUE ! starting poINT of main loop
205  CALL rayatt(x,nw,y,nh,cspln,nxd,dx,dy,xaxd,yaxd,psiaxd,
206  $ xmin,xmax,ymin,ymax,psivl,thet,xn,yn,pds,ifail)
207  IF (ifail .gt. 0)stop
208 c
209 c Check for sufficient accuracy in the poINT spacing for theta.
210 c The accuracy test is based on a relative error in the poloidal
211 c magnetic field of bpers. If the spacing is too large, set theta
212 c back to its previous value theta0, decrease dtheta, and recalculate
213 c the point. The spacing can be repeatedly reduced ihlfmx times.
214 c Excessive accumulation near the x points is avoided by keeping
215 c ihlfmx small; IF ihlfmx is exceeded, the mapping simply CONTINUEs
216 c to the next point.
217 c
218  bp2 = sqrt(pds(2)*pds(2) + pds(3)*pds(3))/xn
219  bpdiff = abs(bp2-bp1)/max(bp2,bp1)
220  IF(bpdiff .ge. bpers) THEN
221  ihalf = ihalf + 1
222  IF(ihalf .le. ihlfmx) THEN
223  thet = theta0
224  dthet = stfrac*dthet
225  thet = thet + dthet
226  go to 10
227  ENDIF
228  ENDIF
229 c
230 c Increment ipts and check IF it has exceeded the DIMENSIONs.
231 c Fill in the new poINT (xc,zc) and CONTINUE.
232 c
233  ipts = ipts + 1
234  IF(ipts .gt. ntmx1)then
235  IF (jcount .ge. jcmax)then
236  WRITE(6,*)'error in cntourp--cnt2' ! quit
237  stop
238  ENDIF
239  go to 5 ! restart
240  ENDIF
241  bp1 = bp2
242  xemin = min(xn,xemin)
243  xemax = max(xn,xemax)
244  yemin = min(yn,yemin)
245  yemax = max(yn,yemax)
246  IF(xn .eq. xemax) yxmax = yn
247  IF(xn .eq. xemin) yxmin = yn
248  IF(yn .eq. yemax) xymax = xn
249  IF(yn .eq. yemin) xymin = xn
250  xc(ipts) = xn
251  yc(ipts) = yn
252  gsq(ipts)= (pds(2)**2+pds(3)**2)/xn**2
253  ihalf=0
254 c
255 c Determine the increment of thet for the next point
256 c The angle increment is limited to dang*2pi*(rad0/rad) and to an
257 c arclength of arcl metres.
258  rad = sqrt((xn-xaxd)**2 + (yn-yaxd)**2)
259  dthet1 = dthet0*(rad0/rad)
260  dthet2 = arcl/rad
261  dthet = min(dthet1,dthet2)
262  theta0 = thet
263  thet = thet + dthet
264  IF(thet .lt. twopi) go to 10
265 c END of main loop
266 c.........................................................................
267 c Close the contour at theta = twopi.
268 c Eliminate the last calculated poINT IF it is too CLOSE to the
269 c closure point.
270 c.........................................................................
271  IF(ipts .gt. 1) THEN
272  ipt1 = ipts - 1
273  xcfn = xc(1)
274  ycfn = yc(1)
275  deli2 = (xc( 2 )-xc( 1 ))*(xc( 2 )-xc( 1 ))
276  $ + (yc( 2 )-yc( 1 ))*(yc( 2 )-yc( 1 ))
277  deln2 = (xc(ipts)-xc(ipt1))*(xc(ipts)-xc(ipt1))
278  $ + (yc(ipts)-yc(ipt1))*(yc(ipts)-yc(ipt1))
279  delfn2 = (xcfn -xc(ipts))*(xcfn -xc(ipts))
280  $ + (ycfn -yc(ipts))*(ycfn -yc(ipts))
281  delti2 = deli2*delta*delta
282  deltn2 = deln2*delta*delta
283  IF((delfn2 .ge. delti2) .and. (delfn2 .ge. deltn2))
284  $ ipts = ipts + 1
285  ENDIF
286  thet = twopi
287  xc(ipts) = xc(1)
288  yc(ipts) = yc(1)
289  gsq(ipts) = gsq(1)
290 c..........................................................................
291 c Write error messages IF bpers needed to be increased.
292 c..........................................................................
293  IF(jcount .gt. 1) THEN
294  WRITE(3,2010) psivl,jcount,jcmax,bperr,bpers
295  WRITE(3,2020) xaxd,yaxd,psiaxd,psivl,ipts
296  WRITE(3,2021) ipts,(xc(ik),ik=1,ipts)
297  WRITE(3,2022) ipts,(yc(ik),ik=1,ipts)
298  ENDIF
299  2010 FORMAT(1x,'needed to increment bpers for psivl =',e12.5,/,1x
300  $ ,'jcount,jcmax,bperr,bpers ='
301  $ ,1x,i3,1x,i3,2x,e12.5,2x,e12.5)
302  2020 FORMAT(1x,/,'xaxd/yaxd/psiaxd',/,1x,3e16.8,/,
303  $ 'psivl',1x,e16.8,/,'ipts',i7,/)
304  2021 FORMAT(/,1x,'(xc(i),i=1,',i5,')',/,(10(1x,e11.4)))
305  2022 FORMAT(/,1x,'(yc(i),i=1,',i5,')',/,(10(1x,e11.4)))
306 c..........................................................................
307  END SUBROUTINE cntourp
308 
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17