<filename>code/rt/rt_utils.py
import numpy


def get_rt_consts():
    """
    Get the constants used for RKCK!
    """
    return {
                    "Rav": 6370.,
                    "a": 6378.137,
                    "f": 1./298.257223563,
                    "ai": [0., 1./5., 3./10., 3./5., 1., 7./8.],
                    "ci": [37./378., 0., 250./621., 125./594., 0., 512./1771.],
                    "dci": [(37./378.)-(2825./27648.), 0.-0., (250./621.)-(18575./48384.), (125./594.)-(13525./55296.),0.-(277./14336.), (512./1771.)-(1./4.)],
                    "b2": [(1./5.), 0., 0., 0., 0.],
                    "b3": [3./40., 9./40., 0., 0., 0.],
                    "b4": [3./10., -9./10., 6./5., 0., 0.],
                    "b5": [-11./54., 5./2., -70./27., 35./27., 0.],
                    "b6": [1631./55296., 175./512., 575./13824., 44275./110592., 253./4096.],

                    "Notes" : "Generic constants used in the ray tracing code"
                }

def calc_pos(lati, longi, alti, azim, dist, elev):
    """
    Calculate new position for a given path distance,
    elevation angle and azimuth.
    Input and output positions are in degrees
    """
    # use temp variable to store latitude and longitude
    tlat = lati
    tlon = longi

    # Convert from geodetic to geocentric and find Earth radius
    glat, glon, Re = calc_gd2gc(tlat, tlon, gd2gc=True)
    # Adjusts azimuth and elevation for the oblateness of the Earth
    gaz, gel = calc_azel(tlat, tlon, azim, elev)
    # Pre-calculate sin and cos of lat and lon
    coslat = numpy.cos(numpy.deg2rad(glat))
    sinlat = numpy.sin(numpy.deg2rad(glat))
    coslon = numpy.cos(numpy.deg2rad(glon))
    sinlon = numpy.sin(numpy.deg2rad(glon))

    # Convert from glabal spherical to global cartesian
    rx = (Re + alti) * coslat * coslon
    ry = (Re + alti) * coslat * sinlon
    rz = (Re + alti) * sinlat

    # Convert from local spherical to local cartesian
    sx = -dist * numpy.cos(numpy.deg2rad(gel)) * numpy.cos(numpy.deg2rad(gaz))
    sy = dist * numpy.cos(numpy.deg2rad(gel)) * numpy.sin(numpy.deg2rad(gaz))
    sz = dist * numpy.sin(numpy.deg2rad(gel))

    # Convert from local cartesian to global cartesian
    tx = sinlat * sx + coslat * sz
    ty = sy
    tz = -coslat * sx + sinlat * sz
    sx = coslon * tx - sinlon * ty
    sy = sinlon * tx + coslon * ty
    sz = tz

    # Add vectors in global cartesian system
    tx = rx + sx
    ty = ry + sy
    tz = rz + sz

    # Convert from global cartesian to global spherical
    rho = numpy.sqrt( tx**2. + ty**2. + tz**2. )
    glat = 90. - numpy.rad2deg(numpy.arccos(tz/rho))
    glon = numpy.rad2deg(numpy.arctan2(ty, tx))

    latiout, longiout, Re = calc_gd2gc(glat, glon, gd2gc=False)

    return (latiout, longiout)


def calc_gd2gc(in_lat, in_lon, gd2gc=True):
    """
    in_lat, in_lon
    Convert between geocentric coordinates
    and geodetic (World Geodetic System 1984 (WGS84))
    gd2gc: True, geodetic to geocentric
    gd2gc: False, geocentric to geodetic
    """
    cnst_dict = get_rt_consts()
    # semi-minor axis (polar radius)
    b = cnst_dict["a"]*(1. - cnst_dict["f"])

    # first eccentricity squared
    e2 = cnst_dict["a"]**2./b**2. - 1.


    if gd2gc:
        glat = numpy.rad2deg(numpy.arctan( b**2./cnst_dict["a"]**2. * numpy.tan(numpy.deg2rad(in_lat)) ))
        glon = in_lon
        if (glon > 180.):
            glon = glon - 360.
        # calculate Earth radius at point (uses geocentric latitude)
        rho = cnst_dict["a"] / numpy.sqrt( 1. + e2*numpy.sin(numpy.deg2rad(glat))**2. )
        return (glat, glon, rho)
    else:
        gdlat = numpy.rad2deg(numpy.arctan( cnst_dict["a"]**2./b**2. * numpy.tan(numpy.deg2rad(in_lat)) ))
        gdlon = in_lon
        # calculate Earth radius at point (uses geocentric latitude)
        rho = cnst_dict["a"] / numpy.sqrt( 1. + e2*numpy.sin(numpy.deg2rad(in_lat))**2. )
        return (gdlat, gdlon, rho)

def calc_azel(lati, longi, azim, elev):
    """
    Calculate azimuth and elevation for oblate Earth.
    Input and output positions are in degrees
    """

    # use temp variable to store latitude and longitude
    tlat = lati
    tlon = longi

    # Convert from geodetic to geocentric and find Earth radius
    glat, glon, Re = calc_gd2gc(tlat, tlon, gd2gc=True)
    _del = tlat - glat

    # Ray k-vector
    kxg = numpy.cos(numpy.deg2rad(elev)) * numpy.sin(numpy.deg2rad(azim))
    kyg = numpy.cos(numpy.deg2rad(elev)) * numpy.cos(numpy.deg2rad(azim))
    kzg = numpy.sin(numpy.deg2rad(elev))

    # Correction to the k-vector due to oblateness
    kxr = kxg
    kyr = kyg * numpy.cos(numpy.deg2rad(_del)) + kzg * numpy.sin(numpy.deg2rad(_del))
    kzr = -kyg * numpy.sin(numpy.deg2rad(_del)) + kzg * numpy.cos(numpy.deg2rad(_del))

    # Finally compute corrected elevation and azimuth
    gaz = numpy.rad2deg(numpy.arctan2(kxr,kyr))
    gel = numpy.rad2deg(numpy.arctan(kzr / numpy.sqrt(kxr**2. + kyr**2.) ))

    return (gaz,gel)
