import numpy as np
from datetime import datetime, timezone

def csheet_struc(R,theta,phi,XJSO,YJSO,ctime):
    """
    !Version 2004 08 15
	INPUT:  R,theta,phi  ->in System III
             XJSO,YJSO    ->in JSO
             ctime        ->UNIX cline time
    OUTPUT: ZNS3         ->height of the current sheet in System III
    """
    XH = 47
    dLamVIP4 = 5.92068
    dThetaVIP4 = 0.16615534478986

    RHO = R * np.sin(theta)

    stheta,sphi,phase = JSun(ctime)

    RLT = np.arctan2(YJSO,XJSO) # zero at noon

    dlamprime = (waveDelay(RLT,RHO)
                 + bendBack(RLT,RHO) + dLamVIP4)
    rho1 = np.sqrt((XH*np.tanh(np.abs(XJSO/XH)))**2 + YJSO**2)
    term2 = rho1*np.tan(dThetaVIP4)
    ZNS3 = (term2*np.cos(phi-dlamprime)
            + XJSO*(1.0 - np.tanh(np.abs(XH/XJSO)))*np.tan(stheta))

    return ZNS3


def waveDelay(RLT, RHO):

    B0 = 3.79406064748764e-002
    C1 = 1.34045565128326 #!C1=Domegaj*rho2/V0 with Domegaj=36.26125*pi/180
    rho2 = 83.3625411987305
    B1 = 0.748919927647621; Phi1 = -1.04848085782547
    B2 = 0.614373903493367; Phi2 = -1.95757783830202
    B3 = 0.541435668088830; Phi3 = -1.25514203423947
    B4 = 7.02600475148873e-002; Phi4 = -1.25542899066018
    B5 = 26.1889910697937

    rhoor2 = RHO/rho2
    rhooB5 = RHO/B5

    waveDelay = (B0 - C1*np.log(np.cosh(rhoor2))
                        *(1 + np.tanh(rhooB5)
                        *(B1*np.cos(RLT - Phi1)
                          + B2*np.cos(2*RLT - Phi3)
                          + B3*np.cos(3*RLT - Phi3)
                          + B4*np.cos(4*RLT - Phi4))))

    return waveDelay

def bendBack(RLT, RHO):

    rho0 = 33.0
    A1  = 1.062528141129624e-004 #!divided by rho1 (1.6*rho0=52.8)
    A2  = 5.950493765539412e-005 #!divided by rho1 (1.6*rho0=52.8)
    A3  = 6.788463593277820e-005 #!divided by rho1 (1.6*rho0=52.8)
    A4  = 3.009364679056859e-005 #!divided by rho1 (1.6*rho0=52.8)
    Phi1 = 1.79696364159234
    Phi2 = 0.41408648226213
    Phi3 =-0.22365844912454
    Phi4 =-0.42846193975833


    Cb1n   = -0.18591516533271       #!Cb1=Cb1*rho0
    Cb1logc= -8.534743746501581e-005 #!Cb1n*dlog(dcosh(dabs(1/rho0)))
    ddK2   = -0.01893939393939       #!dK2=-1/(1.6*rho0)

    ddK1 = A1*np.cos(  RLT - Phi1)
    ddK3 = A2*np.cos(2*RLT - Phi2)
    ddK5 = A3*np.cos(3*RLT - Phi3)
    ddK7 = A4*np.cos(4*RLT - Phi4)

    T3 = (np.exp(ddK2*RHO)/(ddK2**2)*(ddK2*RHO-1) -
          np.exp(ddK2    )/(ddK2**2)*(ddK2    -1))
    T4 = ddK1/2 + ddK3/2 + ddK5/2 + ddK7/2
    rhoSq = RHO**2

    bendBack = (Cb1n*np.log(np.cosh(np.abs(RHO/rho0)))
                - Cb1logc + T4*rhoSq + T3*(-ddK1-ddK3-ddK5-ddK7) - T4)

    return bendBack


def JSun(ctime):
    """
    INPUT:  ctime of the data point 
    OUTPUTS: stheta, sphi, latitude and longitude  (in radians) of the Sun in system III (RH).  
    OUTPUTS: phase, Orbital phase angle of Jupiter (in radians) 
    The equations are written in etime the J200 time convention followed by PDS. 
    We first convert ctime to etime

    Last updated August 12, 2002.

    The program first calculates the direction of the sun in non-rotating  
    coordinates.from equations of the type:
    theta=a1*cos(omegay*t)+a2*sin(omegay*t)+a3*cos(2.*omegay*t)+
    a4*sin(2.*omegay*t)+a5
    fphi=b1*cos(omegay*t)+b2*sin(omegay*t)+b3*cos(2.*omegay*t)+
    b4*sin(2.*omegay*t)+b5
    Then we rotate into the System III coordinates

    ctime is in double precision. theta and phi are single precision 
    +variables.

    omega is jupiter's rotation rate
    omegay is jupiter's yearly orbital rate
    """
    radian = np.pi/180.; degree=180./np.pi
    yrjup = .1185652502e2*86400e0*.36525e3
    omega = 870.536/86400.0; omegay=2.*np.pi/yrjup
    year = 86400e0*.36525e3; D360=360.0
    etime1 = -8.25767955817479e8
    three = 3.123*radian
    tan3=0.054560676; sin3=0.054479647; cos3=0.99851488
    aa = [0.14347029, 3.1145815, -0.12025561, 0.093909436,
         -0.39321884e-5, 0.10194945e-3, -0.12799464]
    bb = [-4.5467523, 3.1848875, -0.16329986, -0.09776818,
           0.17556527e-3, -0.01978317, 44.55915]
    x = np.empty_like(bb)

    # First calculate the latitude and longitude in non-rotating Jupiter coordinates.

    # Calculate the best fit theta and fphi
    t = etime(ctime)-etime1
    x[0] = np.cos(omegay*t)
    x[1] = np.sin(omegay*t)
    x[2] = np.cos(2.*omegay*t)
    x[3] = np.sin(2.*omegay*t)
    x[4] = (t/year)**2
    x[5] = t/year
    x[6] = 1.0
    stheta = 0.
    # fphi is phi in Jupiter fixed (non-rotating) coordinate
    fphi = 0.0

    fphi += np.dot(bb,x)
    stheta += np.dot(aa,x)

    # Now rotate the longitude to Jupiter System III
    # First Add the rotation of Jupiter around the Sun.
    # fphi is the phi of the Sun as unspinning Jupiter goes around the sun
    fphi = np.mod(fphi + t/yrjup*360., D360)
    # Next add the rotation of Jupiter around its axis.
    sphi = np.mod(fphi - t*omega, D360)
    if (sphi < 0):
        sphi = sphi + 360.0
    sphi = sphi*radian
    stheta = stheta*radian
    # Now compute the orbital phase (called phi2 or phase here)
    # There are two solutions to the problem. Only one that is close to phi2b is correct.
    if (stheta >= three):
        stheta=(three)
    if (-stheta >= three):
        stheta=-(three)
    phi21=(sphi+np.pi)+np.arccos(np.tan(stheta)/tan3)
    phi22=(sphi+np.pi)-np.arccos(np.tan(stheta)/tan3)
    phi21=np.mod(phi21,2*np.pi)
    phi22=np.mod(phi22,2*np.pi)
    if (phi22 < 0):
        phi22 = phi22 + 2*np.pi
    dphi = fphi + 48.23012
    phi2b = sphi - dphi*radian
    phi2b = np.mod(phi2b, 2*np.pi)
    if (phi2b < 0):
        phi2b = phi2b + 2*np.pi
    phase = phi21
    dif2 = np.abs(phi22-phi2b)
    if dif2 > 350*radian:
        dif2 = 2*np.pi-dif2
    dif1 = np.abs(phi21-phi2b)
    if (dif1 > 350*radian):
        dif1 = 2*np.pi-dif1
    if (dif2 < dif1):
        phase = phi22
    # phase=phi2b
    phase = np.mod(phase, 2*np.pi)

    # Write(*,*)stheta*degree,phase*degree,phi2b*degree
    return stheta,sphi,phase

def etime(ctime):
    Tcor=0.
    if(ctime >= 189302400.000): Tcor = Tcor+10
    if(ctime >= 205027200.000): Tcor = Tcor+1
    if(ctime >= 220924800.000): Tcor = Tcor+1
    if(ctime >= 252460800.000): Tcor = Tcor+1
    if(ctime >= 283996800.000): Tcor = Tcor+1
    if(ctime >= 315532800.000): Tcor = Tcor+1
    if(ctime >= 347155200.000): Tcor = Tcor+1
    if(ctime >= 378691200.000): Tcor = Tcor+1
    if(ctime >= 410227200.000): Tcor = Tcor+1
    if(ctime >= 441763200.000): Tcor = Tcor+1
    if(ctime >= 489024000.000): Tcor = Tcor+1
    if(ctime >= 520560000.000): Tcor = Tcor+1
    if(ctime >= 552096000.000): Tcor = Tcor+1
    if(ctime >= 615254400.000): Tcor = Tcor+1
    if(ctime >= 694224000.000): Tcor = Tcor+1
    if(ctime >= 757382400.000): Tcor = Tcor+1
    if(ctime >= 788918400.000): Tcor = Tcor+1
    if(ctime >= 836179200.000): Tcor = Tcor+1
    if(ctime >= 867715200.000): Tcor = Tcor+1
    if(ctime >= 899251200.000): Tcor = Tcor+1
    if(ctime >= 946684800.000): Tcor = Tcor+1
    if(ctime >= 993945600.000): Tcor = Tcor+1
    if(ctime >= 1041379200.00): Tcor = Tcor+1
    etime = ctime + Tcor - .1072958367816e10
    return etime

def ctimer(iyr,imon,iday,ihr,imin,sec):
    dt = (iyr,imon,iday,ihr,imin,sec)
    ctimer = datetime(*dt, tzinfo=timezone.utc).timestamp()

    # Uncomment if leap seconds should be included
    #  from astropy.time import Time
    #  t = Time(datetime(*dt), scale='utc')
    #  ctimer = t.unix_tai

    # Add epoch differences between Unix time and 1966
    tepoch0 = datetime(1966, 1, 1, tzinfo=timezone.utc).timestamp()
    tepoch1 = datetime(1970, 1, 1, tzinfo=timezone.utc).timestamp()
    epoch_delta = tepoch0-tepoch1 # = 126230400.0
    ctimer += np.abs(epoch_delta)

    return ctimer 

def doy(iyr,imon,iday):
    return datetime(iyr,imon,iday).timetuple().tm_yday

if __name__ == "__main__":

    # Test the calculation for a single point
    R,theta,phi = 10, np.pi, 0
    XJSO, YJSO = 10, 0
    ctime = 1041379200.0
    ZNS3 = csheet_struc(R,theta,phi,XJSO,YJSO,ctime)
    print(f'ZNS3 = {ZNS3}')