geo_global Subroutine

  pure subroutine geo_global(self, r, theta, Rpos, Zpos)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real, intent (out) :: Rpos, Zpos
    real :: rho0, thAdj, rCyl, Rc, Zc
    real :: Cm, Cn, psiN, xAng, yAng, rAng

    rho0 = r/self%surf%aSurf

    thAdj = theta + self%thetaShift

    Cm=(self%surf%delm**2-1)/(self%surf%delm**2+1)
    Cn=(self%surf%deln**2-1)/(self%surf%deln**3+1)

    psiN=rho0**2+Cm*rho0**2+Cn*rho0**3

    if (abs(Cn*cos(3*(thAdj+self%surf%thn)))<10e-6) then
      rCyl=self%surf%aSurf*sqrt(psiN/(1+Cm*cos(2*(thAdj+self%surf%thm))))
    else
      xAng=2*(1+Cm*cos(2*(thAdj+self%surf%thm)))**3-27*Cn**2*psiN*cos(3*(thAdj+self%surf%thn))**2
      yAng=3*Cn*sqrt(3*psiN)*cos(3*(thAdj+self%surf%thn))*sqrt(xAng+2*(1+Cm*cos(2*(thAdj+self%surf%thm)))**3)

      rAng=(1/3.0)*atan2(yAng,xAng)

      rCyl=self%surf%aSurf*(1+Cm*cos(2*(thAdj+self%surf%thm)))/(3*Cn*cos(3*(thAdj+self%surf%thn))) &
           *(cos(rAng)+sqrt(3.0)*sin(rAng)-1)
    end if

    Rc=self%surf%sHorz-(self%surf%sHorz-self%surf%Rmaj)*rho0**2
    Zc=self%surf%sVert-(self%surf%sVert)*rho0**2

    Rpos=Rc+rCyl*cos(thAdj)
    Zpos=Zc+rCyl*sin(thAdj)

  end subroutine geo_global