import sys
import time
import math
from math import pi, sin, cos, tan, atan2, asin, acos

import jdcal
from skyfield.api import Star, Topos, load

twopi = 2 * pi


def getCurrentJd():
    utc = time.gmtime()
    jd = float(sum(jdcal.gcal2jd(utc[0], utc[1], utc[2]))) + (utc[5] + 60. * (utc[4] + 24. * utc[3])) / 86400. - 0.5
    return jd


def distsph(long1, lat1, long2, lat2):
    cs = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos (long1 - long2)
    return acos(cs)
    

def distsphh(long1, lat1, long2, lat2):
    long1, lat1 = hr2rad(long1), deg2rad(lat1)
    long2, lat2 = hr2rad(long2), deg2rad(lat2)
    return rad2deg(distsph(long1, lat1, long2, lat2))
    

# elementary conversion
def hr2rad(a):
    return twopi * a / 24.


def deg2rad(a):
    return twopi * a / 360.


def rad2deg(a):
    return 360. * a / twopi


def rad2hr(a):
    return 24. * a / twopi


def regul(theta):
    # fold theta between [0; 2pi[
    return theta - twopi * math.floor(theta / twopi)


def regul0(theta):
    theta = regul(theta)
    if theta<pi:
        return theta
    else:
        return theta - twopi


def regulh(hr):
    return hr - 24. * math.floor(hr / 24.)


def sregul(theta):
    theta = regul(theta)
    if theta > pi:
        return theta-twopi
    else:
        return theta


# time misc
jd2000 = 2451545.


def getGMST0(jd):
    mjd = jd - jd2000
    return regulh(18.697374558 + 24.06570982441908 * mjd)


getGMST = getGMST0


def getMeanSiderealTimeFull(ts, yyyy, mm, dd, hr, min, sec, long):
    t = ts.utc(yyy, mm, dd, hr, min, sec)
    return hr2rad(t.gmst) + long


def getMeanSiderealTimeNow(ts, long):
    t = ts.now()
    return hr2rad(t.gmst) + long


def getMeanSiderealTime(t):
    return t.gmst

def getApparentSiderealTime(t):
    return t.gast

def getSiderealTime(t):
    return getApparentSiderealTime(t)


def getLocalSiderealTime(t, long):
    return t.gmst + rad2hr(long)


# conversion between coordinate systems
def radec2hadec(long, lat, t, ra, dec):
    # wiki "Hour Angle"
    return regul(hr2rad(getSiderealTime(t)) + long - ra), dec


def hadec2radec(long, lat, t, ha, dec):
    # wiki "Hour Angle"
    return regul(hr2rad(getSiderealTime(t)) + long - ha), dec


def hadec2azelev(long, lat, ha, dec):
    # from wiki "coordonnees celeste"
    az = atan2(sin(ha), (cos(ha) * sin(lat) - tan(dec) * cos(lat))) + pi
    el = asin(sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(ha))
    return regul(az), regul(el)


def azelev2hadec(long, lat, az, el):
    # from wiki "coordonnees celeste"
    az = az - pi
    ha = atan2(sin(az), cos(az) * sin(lat) + tan(el) * cos(lat))
    dec = asin(sin(lat) * sin(el) - cos(lat) * cos(el) * cos(az))
    return sregul(ha), sregul(dec)


def radec2azelev(long, lat, t, ra, dec):
    ha, dec = radec2hadec(long, lat, t, ra, dec)
    az, elev = hadec2azelev(long, lat, ha, dec)
    return regul(az), regul(elev)


def azelev2radec(long, lat, t, az, el):
    ha, dec = azelev2hadec(long, lat, az, el)
    return hadec2radec(long, lat, t, ha, dec)


################################################################
# whole file conversion utilities

def convert_file_azelev2hadec(file, long, lat):
    # print("# convert")
    f = open(file, 'r')
    for line in f:
        if line[0] != '#':
            a = line.split()
            t = float(a[0])
            az = deg2rad(float(a[1]))
            el = deg2rad(float(a[2]))
            ha, dec = azelev2hadec(long, lat, az, el)
            print(t, rad2hr(ha), rad2deg(dec))
        # else: print(line.rstrip('\n'))


def convert_file_hadec2azelev(file, long, lat):
    # print("# convert")
    f = open(file, 'r')
    for line in f:
        if line[0] != '#':
            a = line.split()
            t = float(a[0])
            a1 = hr2rad(float(a[1]))
            a2 = deg2rad(float(a[2]))
            b1, b2 = hadec2azelev(long, lat, a1, a2)
            print(t, rad2deg(b1), rad2deg(b2))
        # else: print(line.rstrip('\n'))


def convert_file_radec2azelev(file, long, lat, ts):
    # print("# convert")
    f = open(file, 'r')
    for line in f:
        if line[0] != '#':
            a = line.split()
            t = float(a[0])
            a1 = hr2rad(float(a[1]))
            a2 = deg2rad(float(a[2]))
            b1, b2 = radec2azelev(long, lat, ts.tt_jd(t), a1, a2)
            print(t, rad2deg(b1), rad2deg(b2))
        # else: print(line.rstrip('\n'))


def convert_file_azelev2radec(file, long, lat, ts):
    # print("# convert")
    f = open(file, 'r')
    for line in f:
        if line[0] != '#':
            a = line.split()
            t = float(a[0])
            az = deg2rad(float(a[1]))
            el = deg2rad(float(a[2]))
            ra, dec = azelev2radec(long, lat, ts.tt_jd(t), az, el)
            print(t, rad2hr(ra), rad2deg(dec))
        # else: print(line.rstrip('\n'))


# simple direct conversions
# in natural units

def convert_azelev2radec(long, lat, t, az, elev):
    az = deg2rad(az)
    elev = deg2rad(elev)
    b1, b2 = azelev2radec(long, lat, t, az, elev)
    b1 = rad2hr(b1); b2 = rad2deg(b2)
    print(t, b1, "hr  ", b2, "°")
    return b1, b2


def convert_azelev2hadec(long, lat, t, az, elev):
    az = deg2rad(az)
    elev = deg2rad(elev)
    b1, b2 = azelev2hadec(long, lat, az, elev)
    b1 = rad2hr(b1); b2 = rad2deg(b2)
    print(t, b1, "hr  ", b2, "°")
    return b1, b2



def convert_hadec2radec(long, lat, t, ha, dec):
    ha = hr2rad(ha)
    dec = deg2rad(dec)
    b1, b2 = hadec2radec(long, lat, t, ha, dec)
    b1 = rad2hr(b1); b2 = rad2deg(b2)
    print(t, b1,  "hr  ", b2, "°")
    return b1, b2


def convert_hadec2azelev(long, lat, t, ha, dec):
    ha = hr2rad(ha)
    dec = deg2rad(dec)
    b1, b2 = hadec2azelev(long, lat, ha, dec)
    b1 = rad2deg(b1); b2 = rad2deg(b2)
    print(t, b1, "deg  ", b2, "°")
    return b1, b2