# -*- coding: utf-8 -*-
"""
Created on Mon Dec 04 17:06:16 2017

@author: N.Sirdeshmukh
"""

    # Equation 17, coordinates are already translated to origin of triangle
    
    Azprime = atan2(xOrigin, yOrigin)
    
    # Equation 18
    
    rho = sqrt((pow(xOrigin,2) + pow(yOrigin,2)))
    
    # Adjust Azprime to fall within 0 to 120 degrees
    
    Azprime_adjust_multiples = 0;
    while Azprime < 0.0:
        Azprime += DEG120
        Azprime_adjust_multiples-=1
    while (Azprime > DEG120):
        Azprime -= DEG120
        Azprime_adjust_multiples+=1


    #Equation 19
    
    AG = (pow(RPRIME,2) * pow(TANLOWERG,2)) / (2 * ((1/(tan(Azprime))) + COTTHETA))

    # 4 iterations, Azprime (plane) converges to Az (sphere)
    
    for i in range(4):
    
        H = acos((sin(Azprime) * SINUPPERG * COSLOWERG) - (cos(Azprime) * COSUPPERG))

        FAZ = (AG - G - H - Azprime + DEG180)
    
        FPRIMEAZ = (((cos(Azprime) * SINUPPERG * COSLOWERG) + (sin(Azprime) * COSUPPERG)) / (sin(H))) - 1
    
        DeltaAzprime = -FAZ/(FPRIMEAZ)
        
        Azprime = Azprime + DeltaAzprime


    Az = Azprime
 

    # Equations 9-11, 23 to obtain z 
    
    q = atan((TANLOWERG)/(cos(Az) + (sin(Az)*COTTHETA)))
    
    
    # eq 10
    
    dprime = ((RPRIME * TANLOWERG) / (cos(Az) + (sin(Az) * COTTHETA)))

    # eq 11
    
    f = dprime / (2.0 * RPRIME * sin(q / 2.0))
    
    #eq 23, obtain z- does NOT give the same Z value as the z 
    # in the ISEA forward projection... slightly off !
    
    z = 2 * asin((rho)/(2*RPRIME*f))
 
    
    # Add back 120 degree adjustments to Az

    Az += DEG120 * Azprime_adjust_multiples
     

    # Adjust Az to be clockwise from north (needed for final calculation)
    if (face >=1 and face<=5) or (face>=11 and face<=15):
        if (Az <0):
            Az = (M_PI - (Az * (-1))) + M_PI 
    else:
        if (Az <0):
            Az = M_PI - (Az * (-1))
        else:
            Az = Az + M_PI
    
    # R is the radius of the Earth ellipsoid GRS80

    z = z * R
    
    # triangle center of face on which point lies 
    center = icostriangles[face]

    # Convert to final latitude & longitude
    # using geographicLib Python library


    geod = Geodesic(R, flattening)
    d = geod.Direct(center.lat * RAD2DEG, center.lon * RAD2DEG, Az * RAD2DEG, z)
    
    
    return d['lat2'], d['lon2']