-
Notifications
You must be signed in to change notification settings - Fork 6
/
regcoil_init_Fourier_modes_mod.f90
129 lines (102 loc) · 3.19 KB
/
regcoil_init_Fourier_modes_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
module regcoil_init_Fourier_modes_mod
implicit none
contains
subroutine regcoil_init_Fourier_modes(mpol, ntor, mnmax, xm, xn, include_00)
implicit none
integer :: mpol, ntor, mnmax
integer, dimension(:), allocatable :: xm, xn
logical, intent(in) :: include_00
integer :: jn, jm, index, iflag
integer, dimension(:), allocatable :: xm_temp, xn_temp
! xm is nonnegative.
! xn can be negative, zero, or positive.
! When xm is 0, xn must be 0 or positive.
mnmax = mpol*(ntor*2+1) + ntor
if (include_00) mnmax = mnmax + 1
if (allocated(xm)) deallocate(xm)
allocate(xm(mnmax),stat=iflag)
if (iflag .ne. 0) stop 'init_Fourier Allocation error!'
if (allocated(xn)) deallocate(xn)
allocate(xn(mnmax),stat=iflag)
if (iflag .ne. 0) stop 'init_Fourier Allocation error!'
xm = 0
xn = 0
! Handle the xm=0 modes:
index = 0
if (include_00) index = 1
do jn = 1, ntor
index = index + 1
xn(index)=jn
end do
! Handle the xm>0 modes:
do jm = 1,mpol
do jn = -ntor, ntor
index = index + 1
xn(index) = jn
xm(index) = jm
end do
end do
if (index .ne. mnmax) then
print *,"Error! index=",index," but mnmax=",mnmax
stop
end if
end subroutine regcoil_init_Fourier_modes
subroutine regcoil_init_Fourier_modes_sensitivity &
(mpol,ntor,mnmax,nomega,xm,xn,omega,sensitivity_symmetry_option)
implicit none
integer :: mpol, ntor, mnmax, iomega, i, nomega
integer, dimension(:), allocatable :: xm, xn, omega
integer :: minSymmetry, maxSymmetry, sensitivity_symmetry_option
integer :: jn, jm, iflag
! xm is nonnegative.
! xn can be negative, zero, or positive.
! When xm is 0, xn must be non-negative
mnmax = mpol*(ntor*2+1) + ntor+1
! nomega is the length of the number of fourier coefficients
select case (sensitivity_symmetry_option)
case (1) ! stellarator symmetric
nomega = mnmax*2
minSymmetry = 1
maxSymmetry = 2
case (2) ! even in theta and zeta
nomega = mnmax*2
minSymmetry = 3
maxSymmetry = 4
case (3) ! no symmetry
nomega = mnmax*4
minSymmetry = 1
maxSymmetry = 4
end select
allocate(xm(nomega),stat=iflag)
if (iflag .ne. 0) stop 'Allocation error!'
allocate(xn(nomega),stat=iflag)
if (iflag .ne. 0) stop 'Allocation error!'
allocate(omega(nomega),stat=iflag)
if (iflag .ne. 0) stop 'Allocation error!'
! Handle the xm=0 modes:
xm=0
iomega = 0
do jn=1,ntor+1
do i=minSymmetry,maxSymmetry
iomega = iomega + 1
omega(iomega) = i
xn(iomega)=jn-1
enddo
end do
! Handle the xm>0 modes:
do jm = 1,mpol
do jn = -ntor, ntor
do i=minSymmetry,maxSymmetry
iomega = iomega + 1
xn(iomega) = jn
xm(iomega) = jm
omega(iomega) = i
enddo
end do
end do
if (iomega .ne. nomega) then
print *,"Error! iomega=",iomega," but nomega=",nomega
stop
end if
end subroutine regcoil_init_Fourier_modes_sensitivity
end module regcoil_init_Fourier_modes_mod