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(*,*) "<cos**2(theta)*sin**2(theta)> =", 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