V3FIT
descur.f
1  PROGRAM descur
2 c THIS IS PROGRAM DESCUR - WHICH USES A STEEPEST DESCENT
3 c ALGORITHM TO FIND A LEAST SQUARES APPROXIMATION TO AN
4 c ARBITRARY 3-D SPACE CURVE. ANGLE CONSTRAINTS, BASED ON A
5 c MINIMUM SPECTRAL WIDTH CRITERION, ARE APPLIED.
6 c THE CONSTRAINT IS SATISFIED BY TANGENTIAL VARIATIONS
7 c ALONG THE CURVE.
8 c
9 c THE MAIN PROGRAM SETS UP THE INITIAL POINTS TO BE FIT AND
10 c THEN CALLS THE SUBROUTINE SCRUNCH, WHICH DOES THE ACTUAL
11 c CURVE-FITTING.
12 c
13 c***********************************************************************
14 c REVISION 1: January 26, 1989
15 c Compute Fit to individual toroidal planes separately
16 c and THEN Fourier decompose in phi, using optimized
17 c angle representation
18 c
19 c REVISION 2: July, 1995
20 c Generalized for up-down asymmetric CASE
21 c
22 c REVISION 3: July, 1997
23 c Converted to Fortran-90
24 c Eliminated Lagrange multiplier constraint
25 c by using an explicit representation
26 c
27 c
28 c PARAMETER VALUE INTERPRETATIONS:
29 c
30 c MPOL: ACTUAL NO. OF POLOIDAL MODES USED IN FIT FOR R and Z
31 c MRHO: ACTUAL NO. OF POLOIDAL MODES IN FIT FOR QUASI-POLAR RHO
32 c NTHETA: ACTUAL NO. OF THETA (POLOIDAL) POINTS USED IN FIT
33 c NPHI: ACTUAL NO. OF TOROIDAL PLANES IN FIT
34 c
35 c***********************************************************************
36 C-----------------------------------------------
37 C M o d u l e s
38 C-----------------------------------------------
39  USE vname0
40  USE vname1
41  IMPLICIT NONE
42 C-----------------------------------------------
43 C L o c a l P a r a m e t e r s
44 C-----------------------------------------------
45  INTEGER, PARAMETER :: nphi20 = 1 + nv/2
46  INTEGER, PARAMETER :: mntot = mu*(2*nphi20 + 1)
47  INTEGER :: IUNIT=33
48 C-----------------------------------------------
49 C L o c a l V a r i a b l e s
50 C-----------------------------------------------
51  INTEGER :: niter, nstep, itype, ncount,
52  1 i, k, j, m, n, mexp = 4, ns, iread, numargs, ier_flag
53  REAL(rprec), DIMENSION(nu*nv) :: rin, zin
54  REAL(rprec), DIMENSION(mntot) :: rbc, rbs, zbc, zbs
55  REAL(rprec), DIMENSION(0:nv) :: rmnaxis, zmnaxis
56  REAL(rprec), DIMENSION(0:mu-1,-nphi20:nphi20,2) ::
57  1 rbdy_3d, zbdy_3d
58  REAL(rprec), DIMENSION(4) :: rbdy, zbdy, rdshape, zdshape,
59  1 rbean, zbean, rbelt, zbelt, rellip, zellip, rsqr, zsqr
60  REAL(rprec), DIMENSION(0:11) :: rd3dc, zd3dc, rd3ds, zd3ds
61  REAL(rprec), DIMENSION(0:3,0:2) :: rheliac, zheliac
62  REAL(rprec) :: ftol, zeta, texp, dels, A2, VP,
63  1 dumr, dumphi, dumz, arg, PSI, EPS, KAPPA, BDIA, R1, A1, QPSI
64 ! 2 , XI, DELTA_R, rmin, rmax, denom
65  CHARACTER :: datafile*80
66  CHARACTER*120 :: command_arg(10)
67  CHARACTER*1 :: ch_yn
68 C-----------------------------------------------
69 
70  data rdshape/3.0, 0.991, 0.136, 0./
71  data zdshape/0.0, 1.409, -0.118, 0./
72  data rbean/3.0, 1.042, 0.502, -0.0389/
73  data zbean/0.0, 1.339, 0.296, -0.0176/
74  data rbelt/3.0, 0.453, 0., 0./
75  data zbelt/0.0, 0.60, 0., 0.196/
76  data rellip/3.0e0, 1.0e0, 0., 0./
77  data zellip/0.0, 3.0e0, 0., 0./
78  DATA rd3dc/1.600721e+00, 5.675630e-01, 8.285081e-02, 8.140693e-03,
79  1 -5.756955e-03, 6.116248e-03, -9.446896e-04, 2.829636e-04,
80  2 9.626338e-04, 1.882352e-04, 2.479652e-04, 2.384830e-04 /
81  DATA zd3ds/0.000000e+00, 1.022053e+00,-6.887837e-02,-1.375168e-02,
82  1 1.524003e-02,-3.278997e-03, -3.202863e-03, 2.162696e-03,
83  2 -2.481337e-04,-8.761237e-04, 2.479652e-04, 2.384830e-04 /
84  DATA rd3ds/0.000000e+00, 5.456112e-02, 4.568940e-02,-1.259973e-02,
85  1 1.985606e-03, 1.203369e-03,-5.408820e-04, 6.187615e-04,
86  2 1.184803e-03, 1.187193e-04, 5.690797e-04, 7.653277e-05 /
87  DATA zd3dc/-1.274712e-01,5.456112e-02, 7.422938e-04,-1.394692e-02,
88  1 -9.163748e-04, 4.643816e-03,-7.520882e-04,-9.033089e-04,
89  2 1.593907e-03, 2.228495e-04,-5.690797e-04,-7.653277e-05 /
90  data rheliac/4.115, 0.4753, 0., 0., 0.3225, -0.06208, -0.1358,
91  1 -0.04146, 0.0136, 0.0806, -0.0205, 0.0445/
92  data zheliac/0., -0.5045, 0., 0., 0.337, 0.0637, 0.1325, -0.04094
93  1 , 0.01152, 0.06186, -0.03819, - 0.02366/
94  data rsqr/ 3.0, 0.4268, 0., 0.07322/,
95  1 zsqr/ 0.0, 0.4268, 0., -0.07322/
96 
97 c**********************************************************************
98 c CONTROL DATA THAT CAN BE SPECIFIED BY USER
99 c**********************************************************************
100  data niter/1500/
101  data nstep/100/
102 
103  data ftol/1.e-5_dp/
104 
105  data datafile/'none'/
106 
107  twopi = 8*atan(one)
108  hb_parameter = 1
109  mexp = 4
110 !**********************************************************************
111 ! CREATE ARRAYS OF CURVE POINTS, RIN(THETA,PHI) , ZIN(THETA,PHI),
112 ! TO BE FIT TO AN OPTIMIZED FOURIER SERIES BY CALLING "SCRUNCH".
113 ! NOTE: N IS IN UNITS NORMED TO NFP (NUMBER OF FIELD PERIODS) IN 20 LOOP
114 !**********************************************************************
115 !
116 ! Read in data from file or command line i numargs > 0
117 !
118  CALL getcarg(1, command_arg(1), numargs)
119  IF (numargs .GT. 0) THEN
120 !
121 ! READ INPUT DATA FILE ON COMMAND LINE, RATHER THAN INTERACTIVE MODE
122  CALL read_indata(command_arg(1), iunit, ier_flag,
123  1 nphi20, rbdy_3d, zbdy_3d)
124  IF (ier_flag .ne. 0) stop 'ERROR READING INPUT FILE'
125  GO TO 1006
126  ENDIF
127 
128 !
129 ! INTERACTIVE MODE
130 !
131  WRITE (6, '(a)') ' Enter spectral convergence parameter'
132  WRITE (6, '(a)') ' =0 for polar, =1 for equal arclength'
133  WRITE (6, '(a)') ' >1 for smaller spectral width'
134  WRITE (6, '(a)', advance='no')
135  1 ' =4 corresponds to STELLOPT choice): '
136  READ (*,*) texp
137  mexp = texp
138 
139  hb_parameter = 1
140 
141  WRITE (6, '(a,/,a)', advance='no')
142  1 ' Use (default) VMEC-compatible compression (V)',
143  2 ' or Advanced Hirshman-Breslau compression (A): '
144  READ (*,'(a)') ch_yn
145 
146  IF (ch_yn == 'A' .or. ch_yn == 'a') hb_parameter = 0
147  nfp = 1
148 
149  WRITE (6, '(5(a,/))')
150  1 ' Select source of curve data to be fit: ',
151  2 ' 0 : Points from file',
152  3 ' 1 : Fourier coeffs from file (wout file allowed)',
153  4 ' 2 : Solove''ev Equilibrium',
154  5 ' 3 : Assorted 2-D shapes'
155 
156  READ (*,*) itype
157 
158  IF (itype .eq. 3) THEN
159  WRITE (*, *)
160  1 ' Shape to fit: 1=ellipse; 2=D; 3=bean; 4=belt; 5=square;',
161  2 ' 6=D3D-asym; 7=heliac)'
162  READ (*, *) itype
163  IF (itype<1 .or. itype>7) itype = 1
164  IF (itype .eq. 7) THEN
165  WRITE (*, *)
166  1 ' Enter toroidal cross section to plot in degrees: '
167  READ (*, *) zeta
168  zeta = zeta*twopi/360
169  ENDIF
170  GOTO 1001
171 
172  ELSE IF (itype .eq. 0) THEN
173  WRITE(*, '(a)', advance='no')
174  1 ' Enter file name containing R, PHI, Z data to fit:'
175  READ(*, *) datafile
176  IF (datafile .eq. 'none') stop 'Must enter datafile name!'
177  iread = 0
178  GOTO 1002
179  ELSE IF (itype .eq. 1) THEN
180  WRITE(6, '(a)', advance='no')
181  1 ' Enter full file name containing Rmn, Zmn',
182  2 ' coefficients to convert: '
183  READ(*, *) datafile
184  IF (datafile .eq. 'none') stop 'Must enter datafile name!'
185  GOTO 1003
186  ELSE IF (itype .eq. 2) THEN
187  WRITE(6, '(a,/,a)')
188  1 ' Normalized Soloveev boundary: ',
189 ! 2 ' Psi = .5*(1 + xi*r**2)*z**2 + (Psi/delta_r**4)*(r**2 - 1)**2',
190  2 ' Psi = eps^-2[(r^2 + bdia^2)(z/kappa)^2 + (r^2 - 1)^2/4] '
191  WRITE(6, '(a)')
192 ! 1 ' Psi (>0), Xi, Delta-r (<1): (blank or tab delimited)'
193  1 ' Psi (0<Psi<1), eps (<0.5), kappa, bdia: (blank delimited)'
194 ! READ (*,*) PSI, XI, DELTA_R
195  READ (*,*) psi, eps, kappa, bdia
196  GOTO 1004
197  ELSE
198  stop 'Must enter a number between 0 and 2!'
199  END IF
200 
201 
202  1001 CONTINUE
203 
204  ntheta = nu
205  nphi = 1
206  nphi2 = 1+nphi/2
207  mpol = mu
208  mpol1 = mpol-1
209 
210  IF (itype .LT. 6) THEN !!bdy coefficients: m<4
211  SELECT CASE (itype)
212  CASE (1)
213  rbdy = rellip
214  zbdy = zellip
215  CASE (2)
216  rbdy = rdshape
217  zbdy = zdshape
218  CASE (3)
219  rbdy = rbean
220  zbdy = zbean
221  CASE (4)
222  rbdy = rbelt
223  zbdy = zbelt
224  CASE (5)
225  rbdy = rsqr
226  zbdy = zsqr
227 
228  END SELECT
229  rin(:ntheta) = rbdy(1)
230  zin(:ntheta) = zbdy(1)
231  DO m = 2, 4
232  DO j = 1, ntheta
233  arg = (m-1)*twopi*(j - 1)/real(ntheta,rprec)
234  rin(j) = rin(j) + rbdy(m)*cos(arg)
235  zin(j) = zin(j) + zbdy(m)*sin(arg)
236  END DO
237  END DO
238 
239  ELSE IF (itype .EQ. 6) THEN
240  rin(:ntheta) = zero; zin(:ntheta) = zero
241 
242  DO m = 0, 11
243  DO j = 1, ntheta
244  arg = twopi*(j-1)/real(ntheta,rprec)
245  rin(j)=rin(j) + rd3dc(m)*cos(m*arg) + rd3ds(m)*sin(m*arg)
246  zin(j)=zin(j) + zd3ds(m)*sin(m*arg) + zd3dc(m)*cos(m*arg)
247  END DO
248  END DO
249 
250  ELSE IF (itype .EQ. 7) THEN
251  rin(:ntheta) = zero; zin(:ntheta) = zero
252 
253  DO m = 0, 2
254  DO n = 0, 3
255  DO j = 1, ntheta
256  arg = twopi*(j-1)/real(ntheta,rprec)
257  rin(j) = rin(j) + rheliac(n,m)*cos(m*arg - n*zeta)
258  zin(j) = zin(j) + zheliac(n,m)*sin(m*arg - n*zeta)
259  END DO
260  END DO
261  END DO
262 
263  ENDIF
264  3852 FORMAT(4e18.11)
265 
266  GOTO 4000
267 
268  1002 CONTINUE
269 
270  iread = iread+1
271  OPEN(unit=20,file=trim(datafile),status='old', iostat=itype)
272  IF (itype .ne. 0) stop 'error reading datafile'
273  READ (20, *, iostat=itype) ntheta, nphi, nfp
274  IF (itype .ne. 0) THEN
275  print *,'Error reading first line of ',trim(datafile)
276  stop
277  END IF
278  IF (ntheta .gt. nu) stop 'ntheta > nu'
279  IF (nphi .gt. nv) stop 'nphi > nv'
280  mpol = mu
281  mpol1 = mpol-1
282  nphi2 = 1+nphi/2
283  i = 0
284  DO k = 1, nphi
285  DO j = 1, ntheta
286  i = i + 1
287  IF (iread .eq. 1) THEN
288  READ (20, *, iostat=itype) rin(i), dumphi, zin(i)
289  ELSE
290  READ (20, *, iostat=itype) rin(i), zin(i)
291  END IF
292  IF (itype .ne. 0) THEN
293  IF (iread .eq. 1) THEN
294  CLOSE (20)
295  GOTO 1002
296  ELSE
297  print *,' Error reading datafile: iostat = ',itype
298  stop
299  END IF
300  END IF
301  END DO
302  END DO
303 
304  CLOSE(20)
305 
306  IF (nphi .eq. 1
307  1 .and. (all(zin .ge. 0.) .or. all(zin .le. 0.))) THEN
308  ncount = ntheta
309  DO i = ntheta,1,-1
310  IF (zin(i) .eq. 0.) CONTINUE
311  ncount = ncount+1
312  rin(ncount) = rin(i)
313  zin(ncount) =-zin(i)
314  END DO
315  ntheta = ncount
316  IF (ntheta .gt. nu) stop 'ntheta > nu'
317  END IF
318 
319  GOTO 4000
320 
321  1003 CONTINUE
322 
323 !
324 ! try to read data from a wout file
325 !
326  CALL read_wout_curve (datafile, i, nphi20, rbdy_3d, zbdy_3d)
327 
328  IF (i .ne. 0) THEN
329 
330  OPEN(unit=20,file=datafile,iostat=itype)
331  IF (itype .ne. 0) stop 'error reading datafile'
332  READ (20, *) mpol, nphi2, nfp
333  IF (mpol .gt. mu) stop 'mpol-input > mu'
334 
335  DO WHILE (itype .eq. 0)
336  READ(20, *, iostat=itype, END=200) n, m, dumr, dumz
337  IF (itype .ne. 0) cycle
338  IF (nfp .ne. 0) n = n/nfp
339  rbdy_3d(m,n,1) = dumr
340  zbdy_3d(m,n,1) = dumz
341  END DO
342 
343  200 CONTINUE
344  CLOSE(20)
345 
346  END IF
347 
348  1006 CONTINUE
349  ntheta = nu
350  mpol1 = mpol-1
351  nphi = 2*(nphi2+1)
352  IF (nphi2 .le. 0) nphi = 1
353  IF (nphi .gt. nv) stop 'nphi > nv'
354 
355  rin = 0
356  zin = 0
357 
358  DO m = 0, mpol1
359  DO n = -nphi2, nphi2
360  i = 0
361  DO k = 1, nphi
362  zeta = twopi*(k-1)/real(nphi,rprec)
363  DO j = 1, ntheta
364  i = i + 1
365  arg = twopi*(j - 1)/real(ntheta,rprec)
366  rin(i) = rin(i) + rbdy_3d(m,n,1)*cos(m*arg - n*zeta)
367  1 + rbdy_3d(m,n,2)*sin(m*arg - n*zeta)
368  zin(i) = zin(i) + zbdy_3d(m,n,1)*sin(m*arg - n*zeta)
369  1 + zbdy_3d(m,n,2)*cos(m*arg - n*zeta)
370  END DO
371  END DO
372  END DO
373  END DO
374 
375  IF (hb_parameter == 0) THEN !!extra resolution, NOT VMEC compatible
376  mpol = mpol+2
377  ELSE
378  mpol = mpol+1
379  END IF
380 
381  mpol1 = mpol-1
382 
383  GOTO 4000
384 
385  1004 CONTINUE
386 
387  mpol = mu
388  mpol1 = mpol-1
389  ntheta = nu
390  nphi = 1
391  nphi2 = 1+nphi/2
392 ! IF (DELTA_R .LE. ZERO) STOP 'Delta_r must be > 0'
393 ! IF (DELTA_R .GT. one) STOP 'Delta_r must be < 1'
394  IF (psi .LE. zero) stop 'Psi must be > 0'
395  IF (psi .GT. one) stop 'Psi must be <= 1'
396  IF (eps .GT. 0.5) stop 'Eps must be <= 0.5'
397 
398  a1 = eps*sqrt(abs(psi))
399  DO j = 1, ntheta
400  arg = twopi*(j - 1)/real(ntheta,rprec)
401  rin(j) = sqrt(1 + 2*a1*cos(arg))
402  zin(j) = a1*kappa*sin(arg)/sqrt(rin(j)**2 + bdia**2)
403  END DO
404 
405 !COMPUTE q(PSI)/q0
406  ns = 11
407  dels = one/(ns-1)
408  DO i = 1, ns
409  a1 = (i-1)*dels
410  a2 = eps*sqrt(a1)
411  qpsi = 0
412  vp = 0
413  DO j = 1, ntheta
414  arg = twopi*(j-1)/ntheta
415  r1 = sqrt(1 + 2*a2*cos(arg))
416  qpsi = qpsi + one/r1**3
417  vp = vp + one/r1
418  END DO
419 
420  qpsi = qpsi/ntheta
421  vp = (vp/ntheta) * (kappa*eps**2)/2
422  WRITE (6, '(a,1p,3e12.3)')' PSI,Q/Q0, VP: ', a1, qpsi, vp
423  WRITE (33,'(a,1p,3e12.3)')' PSI,Q/Q0, VP: ', a1, qpsi, vp
424  END DO
425 
426 
427 
428 ! rmin = SQRT(one - delta_r**2)
429 ! rmax = SQRT(one + delta_r**2)
430 ! zin = zero
431 ! DO j = 1, ntheta/2
432 ! rin(j) = rmin + REAL(j - 1)*(rmax - rmin)/REAL(ntheta/2)
433 ! rin(j+ntheta/2) =
434 ! 1 rmax - REAL(j - 1)*(rmax - rmin)/REAL(ntheta/2)
435 ! END DO
436 ! DO j = 1, ntheta
437 ! denom = one + XI*rin(j)**2
438 ! IF (denom .le. zero) CYCLE !!avoid separatrix...
439 ! zin(j) = 2*PSI*(one - (rin(j)**2 - one)**2/delta_r**4)/denom
440 ! IF (zin(j) .lt. zero) zin(j) = zero
441 ! zin(j) = SQRT(zin(j))
442 ! IF (j .gt. ntheta/2) zin(j) = -zin(j)
443 ! END DO
444 
445 !
446 ! Attempt to redistribute points more evenly, based on arc-length...
447 !
448 ! DO j = 1, ntheta/2
449 ! length(j) = (rin(j+1) - rin(j))**2 + (zin(j+1) - zin(j))**2
450 ! length(j) = SQRT(length(j)) + 1.E-10_dp
451 ! END DO
452 
453 ! denom = SUM(one/length(:ntheta/2))
454 ! denom = (rmax - rmin)/denom
455 
456 ! DO j = 1, ntheta/2
457 ! rin(j+1) = rin(j) + denom/length(j)
458 ! END DO
459 
460 ! DO j = 2, ntheta/2
461 ! rin(ntheta - j + 2) = rin(j)
462 ! END DO
463 
464 ! DO j = 1,ntheta
465 ! denom = one + XI*rin(j)**2
466 ! IF (denom .le. zero) CYCLE !!avoid separatrix...
467 ! zin(j) = 2*PSI*(one - (rin(j)**2 - one)**2/delta_r**4)/denom
468 ! IF (zin(j) .lt. 0.) zin(j) = 0.
469 ! zin(j) = SQRT(zin(j))
470 ! IF (j .gt. ntheta/2) zin(j) = -zin(j)
471 ! END DO
472 
473 
474 
475  4000 CONTINUE
476 c**********************************************************************
477 c PERFORM OPTIMIZING CURVE-FIT
478 c**********************************************************************
479  CALL scrunch (rin, zin, rbc, zbs, rbs, zbc, rmnaxis,
480  1 zmnaxis, ftol, niter, nstep, mexp, mntot)
481 
482 c**********************************************************************
483 c PRINT OPTIMIZED TRANSFORM COEFFICIENTS TO "OUTCURVE" FILE
484 c AND STORE RIN, ZIN, RBC, ZBS DATA IN PLOTOUT FILE
485 c**********************************************************************
486  CALL printit(rin,zin,rbc,zbs,rbs,zbc,rmnaxis,zmnaxis)
487  CLOSE(10)
488 
489 #ifndef WIN32
490 c**********************************************************************
491 c PLOT FIT TO DATA POINTS
492 c**********************************************************************
493  WRITE (6, '(/,a)', advance='no')
494  1 ' Do you wish to plot this data (Y/N)? '
495  READ (*,*) ch_yn
496 
497  IF (ch_yn == 'Y' .or. ch_yn == 'y') CALL system ("xdes_plot")
498 
499 #endif
500 
501  END PROGRAM descur
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17