geofax Subroutine

  subroutine geofax(rho, theta, t, e, d, nth, ntgrid)
    use geometry, only: rpofrho, geom
    implicit none
    integer, intent(in) :: nth, ntgrid
    real, intent (in) :: rho
    real, dimension(-ntgrid:ntgrid), intent(in) :: theta
    real, intent (out) :: t, e, d
    real, dimension(-nth:nth) :: rgrid
    real :: a, b, aa, bb, rp, hmaxu, hmaxl, rmaxu, rmaxl, rmid, Rlo, Rhi, Rinnr, Routr
    integer :: i
    logical, save :: first=.true.
    if (first) then
       write(*,'("WARNING: GEOFAX USES UNINITIALISED VALUES -- RESULT MAY NOT BE SENSIBLE")')
       first = .false.
    end if
    rp = rpofrho(rho)
    do i = -nth, nth
       rgrid(i) = geom%rfun(rp, theta(i))
    end do

    ! do not treat up-down asymmetry completely for now
    hmaxu = -1.e6 ; hmaxl = 1.e6
    rmaxu = 0. ; rmaxl = 0.

    do i = 0, nth
       a = geom%Rpos(rgrid(i), theta(i)) ; b = geom%Zpos(rgrid(i), theta(i))

       aa = geom%Rpos(rgrid(-i), theta(-i)) ; bb = geom%Zpos(rgrid(-i), theta(-i))

       if (hmaxu < b ) then
          hmaxu = b ; rmaxu = a
       end if

       if (hmaxl > bb ) then
          hmaxl = bb ; rmaxl = aa
       end if

       Rlo = min(a, aa) ; Rhi = max(a, aa)
       !<DD>WARNING : Rinnr has not been defined on first loop through
       if (Rinnr > Rlo) Rinnr = Rlo
       !<DD>WARNING : Routr has not been defined on first loop through
       if (Routr < Rhi) Routr = Rhi
    end do
    a = (geom%Rpos(rgrid(0), theta(0)) - geom%Rpos(rgrid(nth), theta(nth)))
    e = (hmaxu - hmaxl) / a
    ! Ro=(Routr+Rinnr)/2. should = rcenter(rp)
    rmid = (Routr - Rinnr) / 2. ! midplane minor radius of this surface
    ! Using only the upper triangularity:
    ! wasn't properly normalized before 29 Feb 08, and don't use asin:
    !    t = asin(rcenter(rp)-rmax)
    t = (geom%rcenter(rp) - rmaxu) / rmid
    d = geom%rcenter(rp) - geom%rcenter(geom%psi_a)
  end subroutine geofax