V3FIT
getangle.f
1  SUBROUTINE getangle(rval, zval, angle, rcenter, zcenter)
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  REAL(rprec), DIMENSION(ntheta, nphi) :: rval, zval, angle
12  REAL(rprec), DIMENSION(ntheta) :: rcenter, zcenter
13 C-----------------------------------------------
14 C L o c a l V a r i a b l e s
15 C-----------------------------------------------
16  INTEGER :: i, j, iterate
17  REAL(rprec), DIMENSION(nv) ::
18  1 rcos, rsin, zcos, zsin, phiangle
19  REAL(rprec) :: xc, yc, dnum, denom, delangle
20 C-----------------------------------------------
21 c**********************************************************************
22 c Compute angle offset consistent with constraint Z1n = Z1,-n
23 c Note: This is done iteratively, since elongation is unknown
24 c**********************************************************************
25  DO i = 1, nphi
26  DO j = 1, ntheta
27  angle(j,i) = twopi*(j - 1)/real(ntheta,rprec)
28  END DO
29  END DO
30  DO iterate = 1, 5
31  DO i = 1, nphi
32  rcos(i) = zero
33  rsin(i) = zero
34  zcos(i) = zero
35  zsin(i) = zero
36  DO j = 1, ntheta
37  xc = rval(j,i) - rcenter(i)
38  yc = zval(j,i) - zcenter(i)
39  rcos(i) = rcos(i) + cos(angle(j,i))*xc
40  rsin(i) = rsin(i) + sin(angle(j,i))*xc
41  zcos(i) = zcos(i) + cos(angle(j,i))*yc
42  zsin(i) = zsin(i) + sin(angle(j,i))*yc
43  END DO
44  END DO
45 c**********************************************************************
46 c Compute new angles starting from offset phiangle(i)
47 c**********************************************************************
48  dnum = zero
49  denom = zero
50  dnum = sum(zsin(:nphi))
51  denom = sum(rcos(:nphi))
52  IF (denom .ne. zero) THEN
53  elongate = dnum/denom
54  ELSE
55  elongate = 1.e10
56  END IF
57 
58  delangle = zero
59  DO i = 1, nphi
60  phiangle(i) = atan2(elongate*zcos(i)-rsin(i),
61  1 elongate*zsin(i)+rcos(i))
62  delangle = max(delangle,abs(phiangle(i)))
63  angle(:,i) = angle(:,i) + phiangle(i)
64  END DO
65  IF (delangle < 0.02_dp) EXIT
66  END DO
67  WRITE (*, 1010) elongate, raxis(1), zaxis(1)
68  WRITE (3, 1010) elongate, raxis(1), zaxis(1)
69  WRITE (*, 1020) ntheta, nphi
70  WRITE (3, 1020) ntheta, nphi
71  WRITE (*, 1030) mpol1, ntor
72  WRITE (3, 1030) mpol1, ntor
73  1010 FORMAT(' Average elongation = ',1pe10.4,/,' Raxis = ',1pe12.4,
74  1 ' Zaxis = ',1pe12.4)
75  1020 FORMAT(' Number of Theta Points Matched = ',i5,/,
76  1 ' Number of Phi Planes = ',i5)
77  1030 FORMAT(' Max Poloidal Mode Number = ',i5,/,
78  1 ' Max Toroidal Mode Number = ',i5)
79 
80  END SUBROUTINE getangle