program miller ! calculate some quantities relevant to Miller geometry ! see Miller et.al., Phys. Plasmas 5, 973 (1998) implicit none real, allocatable, dimension(:) :: theta real :: R0, r, a_norm, kappa, tri, eps, shift_prime, s_kappa, s_delta, q real :: Vstar, Vstar0,Vstar1, Vstar2, dVstardr, dr, tri2, tri0, eps2, eps0 real :: pi, delta_theta, sum, top, g, denom, denom_approx, maxerr, g_v2 integer :: i integer :: n = 1000 namelist /miller_geo/ R0, r, a_norm, kappa, tri, shift_prime, s_kappa, & s_delta, q ! Initializations: pi = abs(acos(-1.0)) dr=1.e-5 open(unit=20,file="miller.in",status="old") read(20,miller_geo) close(unit=20) eps=r/R0 Vstar1=Vstar(eps,tri) tri2=tri+s_delta*dr/r eps2=(r+dr)/(R0+dr*shift_prime) tri0=tri-s_delta*dr/r eps0=(r-dr)/(R0-dr*shift_prime) Vstar2=Vstar(eps2,tri2) Vstar0=Vstar(eps0,tri0) dVstardr=(Vstar2-Vstar0)/2/dr g_v2= (1+eps/2*shift_prime+s_kappa/2 + r*dVstardr/Vstar1/2) write(*,*) "Volume = 2 pi R0 pi r^2 kappa * g_v" write(*,*) "g_v = ", Vstar1 write(*,*) " " write(*,*) "s_hat_mhd = s_hat/g_v2" write(*,*) "dV/dr = 2 (V/r) * g_v2" write(*,*) "g_v2 = ", g_v2 write(*,*) " " ! Calculate g (related to dPsi/dr): ! Using Eq. 16 and 37 of Miller, we can derive ! ! dpsi/dr = f*kappa*r*g/q/R0, ! ! where f=R*BT is a flux function, and g is calculated below: ! allocate(theta(n)) delta_theta=2*pi/n do i=1, n theta(i)=(i-1)*delta_theta enddo sum=0.0 do i=1, n ! "top" is the denominator of Eq. 37: top=cos(tri*sin(theta(i)))+shift_prime*cos(theta(i)) & +(s_kappa - s_delta*cos(theta(i)) + (1+s_kappa)*tri*cos(theta(i))) & *sin(theta(i))*sin(theta(i)+tri*sin(theta(i))) sum=sum+top/(1+eps*cos(theta(i)+tri*sin(theta(i)))) enddo g=sum*delta_theta/2/pi write(*,*) "dpsi/dr = f*kappa*r*g/q/R0, where f=R_0*B_T0, and" write(*,*) "g = ", g write(*,*) " " write(*,*) "d(beta)/drho = - alpha * ", & g**2/Vstar1**1.5/g_v2*sqrt(kappa)/q**2*a_norm/R0 write(*,*) "d(beta)/drho = - alpha * ", & g**2/Vstar1**1.5/g_v2*sqrt(kappa)/q**2*a_norm/R0 write(*,*) "(8*pi/B_unit^2)dp/drho = - alpha * ", & a_norm/(q**2*R0*kappa**1.5*Vstar1**1.5*g_v2) ! debugging: ! sum=0.0 ! do i=1, n ! sum=sum+cos(theta(i))**2*sin(theta(i))**2 ! enddo ! sum=sum*delta_theta/2/pi ! write(*,*) " =", sum ! debugging: ! sum=0.0 ! maxerr=0.0 ! do i=1, n ! denom=1+eps*cos(theta(i)) ! denom_approx = 1/(1-eps*cos(theta(i))) ! maxerr=max(maxerr,abs(denom-denom_approx)) ! sum=sum+(1-s_delta*cos(theta(i))*sin(theta(i))**2) & ! /denom ! enddo ! sum=sum*delta_theta/2/pi ! write(*,*) "max err in 1-eps*cos(theta) = ", maxerr ! write(*,*) " eps=", eps ! write(*,*) " s_delta=", s_delta ! write(*,*) "eps*s_delta/8 ~ ", sum-1 end program miller real function Vstar(eps,tri) ! ! Calculate the normalized volume of a flux surface. The physical volume ! will be: ! ! Volume = 2*pi*R0*pi*r**2*kappa * Vstar ! implicit none real, allocatable, dimension(:) :: theta real :: tri, eps real :: pi, delta_theta, sum integer :: i integer :: n = 1000 pi = abs(acos(-1.0)) allocate(theta(n)) delta_theta=2*pi/n do i=1, n theta(i)=(i-1)*delta_theta enddo sum=0.0 do i=1, n sum=sum+delta_theta*sin(theta(i))*sin(theta(i)+tri*sin(theta(i))) & *(1+tri*cos(theta(i)))*(1+eps*cos(theta(i)+tri*sin(theta(i)))) enddo Vstar=2*sum/2/pi return end function Vstar