subroutine halo
      real px0(3),px1(3),px2(3),px3(3)
      real py0(3),py1(3),py2(3),py3(3)
      common /haloxy/px0,px1,px2,px3,py0,py1,py2,py3
      data px0/ 2.95e1, 1.93e2, 5.95e2 /
      data px1/ 4.19e0, -4.54e2, -1.92e3 /
      data px2/ -2.40e5, -1.59e6, -4.98e6 /
      data px3/ 3.53e6, 3.53e6, 3.52e6 /
      data py0/ 1.58e1, 5.27e1, 2.14e2 /
      data py1/ -1.06e3, -7.48e2, -6.97e2 /
      data py2/ -2.26e5, -6.21e5, -2.10e6 /
      data py3/ 7.06e6, 7.06e6, 7.06e6 /
      end

      real function halox(x,n)
      real x
      integer n
      real px0(3),px1(3),px2(3),px3(3)
      real py0(3),py1(3),py2(3),py3(3)
      common /haloxy/px0,px1,px2,px3,py0,py1,py2,py3
      halox = px0(n)+px1(n)*x+px2(n)*x**2
      end

      real function haloy(y,n)
      real y
      integer n
      real px0(3),px1(3),px2(3),px3(3)
      real py0(3),py1(3),py2(3),py3(3)
      common /haloxy/px0,px1,px2,px3,py0,py1,py2,py3
      haloy = py0(n)+py1(n)*y+py2(n)*y**2
      end