V3FIT
order.f
1  SUBROUTINE order(rval, zval, xaxis, yaxis, inside)
2 C-----------------------------------------------
3 C M o d u l e s
4 C-----------------------------------------------
5  USE vname0
6  USE vname1
7  IMPLICIT NONE
8 C-----------------------------------------------
9 C D u m m y A r g u m e n t s
10 C-----------------------------------------------
11  INTEGER inside
12  REAL(rprec) xaxis, yaxis
13  REAL(rprec), DIMENSION(*) :: rval, zval
14 C-----------------------------------------------
15 C L o c a l V a r i a b l e s
16 C-----------------------------------------------
17  INTEGER :: i, ip1, i1, j, next
18  REAL(rprec), DIMENSION(nu) :: tempr, tempz
19  REAL(rprec) ::
20  1 newdist,olddist,shortest,saver,savez,residue,x,y,dx,dy
21 C-----------------------------------------------
22 c**********************************************************************
23 c Program ORDER : orders points on a magnetic surface at a
24 c fixed toroidal plane and assigns right-handed circulation
25 c around flux surface. XAXIS, YAXIS: Polar-TYPE axis (must lie
26 c inside curve to check sign of rotation)
27 c
28 c Modified 8/13/97, J. Breslau. No longer quits IF the axis is not
29 c contained by the points. instead, returns inside=0 so inguess can
30 c try again.
31 c**********************************************************************
32  olddist = 1.e20_dp
33  DO i = 1, ntheta - 1
34  ip1 = i + 1
35  i1 = i
36  shortest = 1.e20_dp
37  15 CONTINUE
38  DO j = ip1, ntheta
39  IF (i1 > 1) olddist = (rval(i1-1)-rval(j))**2 + (zval(i1-1)-
40  1 zval(j))**2
41  newdist = (rval(i1)-rval(j))**2 + (zval(i1)-zval(j))**2
42  IF (newdist.le.olddist .and. newdist<shortest) THEN
43  next = j
44  shortest = newdist
45  ENDIF
46  END DO
47 c**********************************************************************
48 c Swap nearest point (next) with current point (ip1)
49 c**********************************************************************
50  IF (shortest .ge. 1.e10_dp) THEN
51  saver = rval(i-1)
52  rval(i-1) = rval(i)
53  rval(i) = saver
54  savez = zval(i-1)
55  zval(i-1) = zval(i)
56  zval(i) = savez
57  i1 = i1 - 1
58  ip1 = ip1 - 1
59  GOTO 15
60  ENDIF
61  saver = rval(ip1)
62  rval(ip1) = rval(next)
63  rval(next) = saver
64  savez = zval(ip1)
65  zval(ip1) = zval(next)
66  zval(next) = savez
67  END DO
68 c**********************************************************************
69 c Check that xaxis,yaxis is inside surface and
70 c ascertain that the angle rotates counterclockwise
71 c using Cauchy theorem in "complex"-plane
72 c**********************************************************************
73  inside = 1
74  residue = zero
75  DO i = 1, ntheta - 1
76  x = 0.5_dp*(rval(i)+rval(i+1)) - xaxis
77  y = 0.5_dp*(zval(i)+zval(i+1)) - yaxis
78  dx = rval(i+1) - rval(i)
79  dy = zval(i+1) - zval(i)
80  residue = residue + (x*dy - y*dx)/(x**2 + y**2 + 1.e-10_dp)
81  END DO
82  x = 0.5_dp*(rval(1)+rval(ntheta)) - xaxis
83  y = 0.5_dp*(zval(1)+zval(ntheta)) - yaxis
84  dx = rval(1) - rval(ntheta)
85  dy = zval(1) - zval(ntheta)
86  residue = residue + (x*dy - y*dx)/(x**2 + y**2 + 1.e-10_dp)
87 
88  IF (residue < (-0.9_dp*twopi)) THEN
89  DO i = 2, ntheta
90  j = ntheta - i + 2
91  tempr(i) = rval(j)
92  tempz(i) = zval(j)
93  END DO
94  rval(2:ntheta) = tempr(2:ntheta)
95  zval(2:ntheta) = tempz(2:ntheta)
96  ELSE IF (abs(residue) < 0.9_dp*twopi) THEN
97  print *, ' mag. axis not enclosed by bndry; trying again'
98  WRITE (3, *) ' mag. axis not enclosed by bndry; trying again'
99  inside = 0
100 c STOP
101  ENDIF
102 
103  END SUBROUTINE order