How to use the treecorr.CelestialCoord function in TreeCorr

def test_angle():

    # Again, let's start with some answers we can get by inspection.
    eq1 = treecorr.CelestialCoord(0.,0.)  # point on the equator
    eq2 = treecorr.CelestialCoord(1.,0.)  # 1 radian along equator
    eq3 = treecorr.CelestialCoord(pi,0.) # antipode of eq1
    north_pole = treecorr.CelestialCoord(0.,pi/2.)  # north pole
    south_pole = treecorr.CelestialCoord(0.,-pi/2.) # south pole

    numpy.testing.assert_almost_equal(north_pole.angleBetween(eq1,eq2), -1.)
    numpy.testing.assert_almost_equal(north_pole.angleBetween(eq2,eq1), 1.)
    numpy.testing.assert_almost_equal(north_pole.angleBetween(eq2,eq3), 1.-pi)
    numpy.testing.assert_almost_equal(north_pole.angleBetween(eq3,eq2), pi-1.)
    numpy.testing.assert_almost_equal(south_pole.angleBetween(eq1,eq2), 1.)
    numpy.testing.assert_almost_equal(south_pole.angleBetween(eq2,eq1), -1.)
    numpy.testing.assert_almost_equal(south_pole.angleBetween(eq2,eq3), pi-1.)
    numpy.testing.assert_almost_equal(south_pole.angleBetween(eq3,eq2), 1.-pi)

    numpy.testing.assert_almost_equal(eq1.angleBetween(north_pole,eq2), pi/2.)
    numpy.testing.assert_almost_equal(eq2.angleBetween(north_pole,eq1), -pi/2.)

    numpy.testing.assert_almost_equal(north_pole.area(eq1,eq2), 1.)
    numpy.testing.assert_almost_equal(north_pole.area(eq2,eq1), 1.)
    numpy.testing.assert_almost_equal(south_pole.area(eq1,eq2), 1.)

numpy.testing.assert_almost_equal(eq1.distanceTo(eq2), 1.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq1), 1.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(eq3), pi)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq3), pi-1.)

    numpy.testing.assert_almost_equal(north_pole.distanceTo(south_pole), pi)

    numpy.testing.assert_almost_equal(eq1.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq3.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(south_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(south_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq3.distanceTo(south_pole), pi/2.)

    c1 = treecorr.CelestialCoord(0.234, 0.342)  # Some random point
    c2 = treecorr.CelestialCoord(0.234, -1.093) # Same meridian
    c3 = treecorr.CelestialCoord(pi + 0.234, -0.342) # Antipode
    c4 = treecorr.CelestialCoord(pi + 0.234, 0.832) # Different point on opposide meridian

    numpy.testing.assert_almost_equal(c1.distanceTo(c1), 0.)
    numpy.testing.assert_almost_equal(c1.distanceTo(c2), 1.435)
    numpy.testing.assert_almost_equal(c1.distanceTo(c3), pi)
    numpy.testing.assert_almost_equal(c1.distanceTo(c4), pi-1.174)

    # Now some that require spherical trig calculations.
    # Importantly, this uses the more straightforward spherical trig formula, the cosine rule.
    # The CelestialCoord class uses a different formula that is more stable for very small
    # distances, which are typical in the correlation function calculation.
    c5 = treecorr.CelestialCoord(1.832, -0.723)  # Some other random point
    # The standard formula is:
    # cos(d) = sin(dec1) sin(dec2) + cos(dec1) cos(dec2) cos(delta ra)
    d = arccos(sin(0.342) * sin(-0.723) + cos(0.342) * cos(-0.723) * cos(1.832 - 0.234))

def test_distance():

    # First, let's test some distances that are easy to figure out
    # without any spherical trig.
    eq1 = treecorr.CelestialCoord(0.,0.)  # point on the equator
    eq2 = treecorr.CelestialCoord(1.,0.)  # 1 radian along equator
    eq3 = treecorr.CelestialCoord(pi,0.) # antipode of eq1
    north_pole = treecorr.CelestialCoord(0.,pi/2.)  # north pole
    south_pole = treecorr.CelestialCoord(0.,-pi/2.) # south pole

    numpy.testing.assert_almost_equal(eq1.distanceTo(eq2), 1.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq1), 1.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(eq3), pi)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq3), pi-1.)

    numpy.testing.assert_almost_equal(north_pole.distanceTo(south_pole), pi)

    numpy.testing.assert_almost_equal(eq1.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq3.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(south_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(south_pole), pi/2.)

print('south.galactic = ',south.galactic())
    el,b = south.galactic()
    numpy.testing.assert_almost_equal(b, -pi/2., decimal=3)

    # The anti-center is at 5h:42.6m, 28.92d
    anticenter = treecorr.CelestialCoord( (5.+45.6/60.) * treecorr.hours,
                                          28.94 * treecorr.degrees )
    print('anticenter.galactic = ',anticenter.galactic())
    el,b = anticenter.galactic()
    numpy.testing.assert_almost_equal(el, pi, decimal=3)
    numpy.testing.assert_almost_equal(b, 0., decimal=3)

    # Now test with some random values using this NASA website to get the correct answers:
    # http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm
    # While I'm at it, they also give me eclipticc coords, so test those as well.
    coord1 = treecorr.CelestialCoord(13.234 * treecorr.degrees, -73.438 * treecorr.degrees)
    el, b = coord1.galactic()
    print('coord1 = ',coord1)
    print('el, b = ',el/treecorr.degrees, b/treecorr.degrees)
    numpy.testing.assert_almost_equal(el, (302.78430-360.) * treecorr.degrees, decimal=5)
    numpy.testing.assert_almost_equal(b, -43.68987 * treecorr.degrees, decimal=5)
    lam, beta = coord1.ecliptic()
    print('lam, beta = ',lam,beta)
    print('lam, beta = ',lam/treecorr.degrees, beta/treecorr.degrees)
    numpy.testing.assert_almost_equal(lam, (310.80568-360.) * treecorr.degrees, decimal=5)
    numpy.testing.assert_almost_equal(beta, -64.87384 * treecorr.degrees, decimal=5)

    coord2_1950 = treecorr.CelestialCoord(122.93 * treecorr.degrees, 16.01 * treecorr.degrees)
    coord2_2000 = coord2_1950.precess(1950., 2000.)
    print('coord2_1950 = ',coord2_1950)
    print('coord2_2000 = ',coord2_2000)
    numpy.testing.assert_almost_equal(coord2_2000.ra, 123.63697 * treecorr.degrees, decimal=5)

def test_projection():

    # Test that a small triangle has the correct properties for each kind of projection
    center = treecorr.CelestialCoord(0.234,0.342)
    cA = treecorr.CelestialCoord(-0.193,0.882)
    cB = treecorr.CelestialCoord(-0.193 + 1.7e-6,0.882 + 1.2e-6)
    cC = treecorr.CelestialCoord(-0.193 - 2.4e-6,0.882 + 3.1e-6)

    a = cB.distanceTo(cC)
    b = cC.distanceTo(cA)
    c = cA.distanceTo(cB)
    A = cA.angleBetween(cB,cC)
    B = cB.angleBetween(cC,cA)
    C = cC.angleBetween(cA,cB)
    E = cA.area(cB,cC)

    #
    # The lambert is supposed to preserve area
    #

    pA = center.project(cA, projection='lambert')
    pB = center.project(cB, projection='lambert')
    pC = center.project(cC, projection='lambert')

def test_distance():

    # First, let's test some distances that are easy to figure out
    # without any spherical trig.
    eq1 = treecorr.CelestialCoord(0.,0.)  # point on the equator
    eq2 = treecorr.CelestialCoord(1.,0.)  # 1 radian along equator
    eq3 = treecorr.CelestialCoord(pi,0.) # antipode of eq1
    north_pole = treecorr.CelestialCoord(0.,pi/2.)  # north pole
    south_pole = treecorr.CelestialCoord(0.,-pi/2.) # south pole

    numpy.testing.assert_almost_equal(eq1.distanceTo(eq2), 1.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq1), 1.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(eq3), pi)
    numpy.testing.assert_almost_equal(eq2.distanceTo(eq3), pi-1.)

    numpy.testing.assert_almost_equal(north_pole.distanceTo(south_pole), pi)

    numpy.testing.assert_almost_equal(eq1.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq3.distanceTo(north_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq1.distanceTo(south_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq2.distanceTo(south_pole), pi/2.)
    numpy.testing.assert_almost_equal(eq3.distanceTo(south_pole), pi/2.)

    c1 = treecorr.CelestialCoord(0.234, 0.342)  # Some random point

center = treecorr.CelestialCoord( (17.+45.6/60.) * treecorr.hours,
                                      -28.94 * treecorr.degrees )
    print('center.galactic = ',center.galactic())
    el,b = center.galactic()
    numpy.testing.assert_almost_equal(el, 0., decimal=3)
    numpy.testing.assert_almost_equal(b, 0., decimal=3)

    # The north pole is at 12h:51.4m, 27.13d
    north = treecorr.CelestialCoord( (12.+51.4/60.) * treecorr.hours,
                                     27.13 * treecorr.degrees )
    print('north.galactic = ',north.galactic())
    el,b = north.galactic()
    numpy.testing.assert_almost_equal(b, pi/2., decimal=3)

    # The south pole is at 0h:51.4m, -27.13d
    south = treecorr.CelestialCoord( (0.+51.4/60.) * treecorr.hours,
                                     -27.13 * treecorr.degrees )
    print('south.galactic = ',south.galactic())
    el,b = south.galactic()
    numpy.testing.assert_almost_equal(b, -pi/2., decimal=3)

    # The anti-center is at 5h:42.6m, 28.92d
    anticenter = treecorr.CelestialCoord( (5.+45.6/60.) * treecorr.hours,
                                          28.94 * treecorr.degrees )
    print('anticenter.galactic = ',anticenter.galactic())
    el,b = anticenter.galactic()
    numpy.testing.assert_almost_equal(el, pi, decimal=3)
    numpy.testing.assert_almost_equal(b, 0., decimal=3)

    # Now test with some random values using this NASA website to get the correct answers:
    # http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm
    # While I'm at it, they also give me eclipticc coords, so test those as well.

numpy.testing.assert_almost_equal(c1.distanceTo(c2), 1.435)
    numpy.testing.assert_almost_equal(c1.distanceTo(c3), pi)
    numpy.testing.assert_almost_equal(c1.distanceTo(c4), pi-1.174)

    # Now some that require spherical trig calculations.
    # Importantly, this uses the more straightforward spherical trig formula, the cosine rule.
    # The CelestialCoord class uses a different formula that is more stable for very small
    # distances, which are typical in the correlation function calculation.
    c5 = treecorr.CelestialCoord(1.832, -0.723)  # Some other random point
    # The standard formula is:
    # cos(d) = sin(dec1) sin(dec2) + cos(dec1) cos(dec2) cos(delta ra)
    d = arccos(sin(0.342) * sin(-0.723) + cos(0.342) * cos(-0.723) * cos(1.832 - 0.234))
    numpy.testing.assert_almost_equal(c1.distanceTo(c5), d)

    # Tiny displacements should have dsq = (dra^2 cos^2 dec) + (ddec^2)
    c6 = treecorr.CelestialCoord(0.234 + 1.7e-9, 0.342)
    c7 = treecorr.CelestialCoord(0.234, 0.342 + 1.9e-9)
    c8 = treecorr.CelestialCoord(0.234 + 2.3e-9, 0.342 + 1.2e-9)

    # Note that the standard formula gets thsse wrong.  d comes back as 0.
    d = arccos(sin(0.342) * sin(0.342) + cos(0.342) * cos(0.342) * cos(1.7e-9))
    print('d(c6) = ',1.7e-9 * cos(0.342), c1.distanceTo(c6), d)
    d = arccos(sin(0.342) * sin(0.342+1.9e-9) + cos(0.342) * cos(0.342+1.9e-9) * cos(0.))
    print('d(c7) = ',1.9e-9, c1.distanceTo(c7), d)
    d = arccos(sin(0.342) * sin(0.342) + cos(0.342) * cos(0.342) * cos(1.2e-9))
    true_d = sqrt( (2.3e-9 * cos(0.342))**2 + 1.2e-9**2)
    print('d(c7) = ',true_d, c1.distanceTo(c8), d)
    numpy.testing.assert_almost_equal(c1.distanceTo(c6)/(1.7e-9 * cos(0.342)), 1.0)
    numpy.testing.assert_almost_equal(c1.distanceTo(c7)/1.9e-9, 1.0)
    numpy.testing.assert_almost_equal(c1.distanceTo(c8)/true_d, 1.0)

# The north pole is at 12h:51.4m, 27.13d
    north = treecorr.CelestialCoord( (12.+51.4/60.) * treecorr.hours,
                                     27.13 * treecorr.degrees )
    print('north.galactic = ',north.galactic())
    el,b = north.galactic()
    numpy.testing.assert_almost_equal(b, pi/2., decimal=3)

    # The south pole is at 0h:51.4m, -27.13d
    south = treecorr.CelestialCoord( (0.+51.4/60.) * treecorr.hours,
                                     -27.13 * treecorr.degrees )
    print('south.galactic = ',south.galactic())
    el,b = south.galactic()
    numpy.testing.assert_almost_equal(b, -pi/2., decimal=3)

    # The anti-center is at 5h:42.6m, 28.92d
    anticenter = treecorr.CelestialCoord( (5.+45.6/60.) * treecorr.hours,
                                          28.94 * treecorr.degrees )
    print('anticenter.galactic = ',anticenter.galactic())
    el,b = anticenter.galactic()
    numpy.testing.assert_almost_equal(el, pi, decimal=3)
    numpy.testing.assert_almost_equal(b, 0., decimal=3)

    # Now test with some random values using this NASA website to get the correct answers:
    # http://lambda.gsfc.nasa.gov/toolbox/tb_coordconv.cfm
    # While I'm at it, they also give me eclipticc coords, so test those as well.
    coord1 = treecorr.CelestialCoord(13.234 * treecorr.degrees, -73.438 * treecorr.degrees)
    el, b = coord1.galactic()
    print('coord1 = ',coord1)
    print('el, b = ',el/treecorr.degrees, b/treecorr.degrees)
    numpy.testing.assert_almost_equal(el, (302.78430-360.) * treecorr.degrees, decimal=5)
    numpy.testing.assert_almost_equal(b, -43.68987 * treecorr.degrees, decimal=5)
    lam, beta = coord1.ecliptic()