-
Notifications
You must be signed in to change notification settings - Fork 0
/
contour.f90
128 lines (108 loc) · 4.61 KB
/
contour.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
module contour_mod
contains
subroutine contour(x_lcfs,z_lcfs,np_lcfs,x_axis,z_axis,psival,x_contour,z_contour)
!given a value of the poloidal flux, psival, this subroutine find the magnetic surface corresponding to this poloidal flux
use constants,only:p_
use constants,only:zero,one,two,twopi
implicit none
integer,intent(in):: np_lcfs
real(p_),intent(in):: x_lcfs(np_lcfs),z_lcfs(np_lcfs),x_axis,z_axis,psival
real(p_),intent(out):: x_contour(np_lcfs),z_contour(np_lcfs)
real(p_),parameter:: xacc=1.0d-6 !tolerance used in bi-section root-finder
real(p_):: x1,x2,z1,z2
real(p_):: slope(np_lcfs),slope2(np_lcfs)
integer:: i
real(p_),parameter:: large_number=1d30
!do i=1,np_lcfs-1
do i=1,np_lcfs
if(x_lcfs(i)-x_axis .eq. 0._p_) then !since I use compiler option which catches all erroneous arithmetic operation, I need to avoid dividing by zero
slope(i)= large_number
else
slope(i)= (z_lcfs(i)-z_axis)/(x_lcfs(i)-x_axis) !the slope for function Z=Z(X)
endif
if(z_lcfs(i)-z_axis .eq. 0._p_) then
slope2(i)=large_number
else
slope2(i)=(x_lcfs(i)-x_axis)/(z_lcfs(i)-z_axis) !the slope for function X=X(Z)
endif
!write(*,*) i,slope(i),slope2(i)
enddo
!!$ write(*,*) maxval(slope),minval(slope)
!!$ write(*,*) maxloc(slope),minloc(slope)
!!$ write(*,*) x_axis,x_lcfs( maxloc(slope)),x_lcfs(minloc(slope))
do i=1,np_lcfs-1 !exclude i=np_lcfs because it is identical to i=1
if(abs(slope(i)).le.1.0_p_) then !use Z=Z(X) function, the reason that I switch between using function X=X(Z) and Z=Z(X) is to aviod large slope.
x1=x_axis
x2=x_lcfs(i) !+0.01 !shift left a little to gurrantee that the range is enough for a root to lie in
x_contour(i)=rtbis(one_dim_psi_func,x1,x2,xacc,x_axis,z_axis,slope(i),psival)
z_contour(i)=zfunc(x_axis,z_axis,slope(i),x_contour(i))
else !switch to using X=X(Z) function
z1=z_axis
z2=z_lcfs(i)
z_contour(i)=rtbis(one_dim_psi_func2,z1,z2,xacc,x_axis,z_axis,slope2(i),psival)
x_contour(i)=xfunc(x_axis,z_axis,slope2(i),z_contour(i))
endif
enddo
x_contour(np_lcfs)=x_contour(1) !i=1 and i=np_lcfs are identical
z_contour(np_lcfs)=z_contour(1) !i=1 and i=np_lcfs are identical
end subroutine contour
function one_dim_psi_func(x_axis,z_axis,slope,psival,x)
!poloidal flux as a function of x on a straight line with slope "slope" in poloidal plane
use constants,only:p_
use magnetic_field, only : psi_func
implicit none
real(p_):: one_dim_psi_func,x,x_axis,z_axis,slope,psival
one_dim_psi_func=psi_func(x,zfunc(x_axis,z_axis,slope,x))-psival
end function one_dim_psi_func
function zfunc(x_axis,z_axis,slope,x) !straight line Z=Z(x) with slope "slope" in poloidal plane starting from the location of magnetic axis
use constants,only:p_
implicit none
real(p_):: zfunc,x,x_axis,z_axis,slope
zfunc=z_axis+slope*(x-x_axis)
end function zfunc
function one_dim_psi_func2(x_axis,z_axis,slope,psival,z) result(fun_val)
!poloidal flux as a function of z on a straight line with slope "slope" in poloidal plane
use constants,only:p_
use magnetic_field, only : psi_func
implicit none
real(p_):: fun_val,z
real(p_):: x_axis,z_axis,slope,psival
fun_val=psi_func(xfunc(x_axis,z_axis,slope,z),z)-psival
end function one_dim_psi_func2
function xfunc(x_axis,z_axis,slope,z) !straight line X=X(Z) with slope "slope" in poloidal plane starting from the location of magnetic axis
use constants,only:p_
implicit none
real(p_):: xfunc, z
real(p_)::x_axis,z_axis,slope
xfunc=x_axis+slope*(z-z_axis)
end function xfunc
FUNCTION rtbis(func,x1,x2,xacc,xmaxis,zmaxis,slope,psival)
!find a root of func by using the bisection method
use constants,only: p_
implicit none
INTEGER, parameter :: JMAX=40
REAL(p_) rtbis, x1, x2, func, xacc, xmaxis, zmaxis, slope, psival
external :: func
INTEGER j
REAL(p_) dx,f,fmid,xmid
fmid=func(xmaxis,zmaxis,slope,psival,x2)
f= func(xmaxis,zmaxis,slope,psival,x1)
! write(*,*) 'f1=', f, 'f2=',fmid
if(f*fmid.ge.0.) stop 'root must be bracketed in rtbis'
if(f.lt.0.)then
rtbis=x1
dx=x2-x1
else
rtbis=x2
dx=x1-x2
endif
do j=1,JMAX
dx=dx*.5
xmid=rtbis+dx
fmid=func(xmaxis,zmaxis,slope,psival,xmid)
if(fmid.le.0.)rtbis=xmid
if(abs(dx).lt.xacc .or. fmid.eq.0.) return
enddo
stop 'too many bisections in rtbis'
end FUNCTION rtbis
end module contour_mod