#!/usr/bin/env python3

import ephem
from datetime import datetime
import xml.dom.minidom
import json
import math
import sys
from pprint import pprint
import argparse


def nearest_time(targettime, t1, t2):
    """Given a target ephem.time and two other datetimes,
       return the time closer to the target.
    """
    d1 = abs(targettime - t1)
    d2 = abs(targettime - t2)
    if d1 <= d2:
        return t1
    return t2


def find_rise_set(observer, obj, d=None):
    """Given an object (like Sun or Moon), find its rising and setting time
       closest to the given date d, either preceding or following it,
       for the observer's location.
       If date isn't specified, use the observer's date.
    """
    if d:
        observer.date = d
    prevrise = observer.previous_rising(obj)
    nextrise = observer.next_rising(obj)
    prevset = observer.previous_setting(obj)
    nextset = observer.next_setting(obj)

    riseset_ret = {}

    # Rise time
    observer.date = nearest_time(observer.date, prevrise, nextrise)
    obj.compute(observer)
    riseset_ret['rise'] = { 'az':   obj.az / ephem.degree,
                            'time': observer.date
                          }

    # Set time
    observer.date = nearest_time(observer.date, prevset, nextset)
    obj.compute(observer)
    riseset_ret['set'] = { 'az':   obj.az / ephem.degree,
                           'time': observer.date
                         }

    return riseset_ret


def find_azimuths(observer):
    riseset = {}

    # Find sunrise and sunset:
    riseset['sun'] = find_rise_set(observer, ephem.Sun())

    # Now find the full moon closest to the date,
    # which may be the next full moon or the previous one.
    lastfull = ephem.previous_full_moon(observer.date)
    nextfull = ephem.next_full_moon(observer.date)
    now = ephem.now()
    if now - lastfull > nextfull - now:
        observer.date = nextfull
    else:
        observer.date = lastfull

    riseset['full moon'] = find_rise_set(observer, ephem.Moon())

    return riseset


def bearing_to(wp1, wp2):
    # https://www.movable-type.co.uk/scripts/latlong.html
    # Don't trust any code you find for this: test it extensively;
    # most posted bearing finding code is bogus.
    # print("bearing from waypoints:", wp1, wp2)
    lat1, lon1 = math.radians(wp1[1]), math.radians(wp1[2])
    lat2, lon2 = math.radians(wp2[1]), math.radians(wp2[2])
    y = math.sin(lon2 - lon1) * math.cos(lat2)
    x = math.cos(lat1) * math.sin(lat2) - \
        math.sin(lat1) * math.cos(lat2) * math.cos(lon2-lon1)

    return math.degrees(math.atan2(y, x)) % 360


def find_alignments(observer, waypoints, year, allpoints=False):
    """Find all the alignments with solstice/equinox sun/moon rise/set.
       Returns a dict: { 'vernal equinox': { 'moon': { 'rise': 94.17... } } }
       of azimuth angles in decimal degrees
    """
    azimuths = {}

    # start_date = ephem.Date('%d/1/1' % year)
    start_date = ephem.Date((year, 1, 1))
    # date= ephem.date((-59000,1,1))

    observer.date = ephem.next_equinox(start_date)
    azimuths['vernal equinox'] = find_azimuths(observer)

    observer.date = ephem.next_solstice(observer.date)
    azimuths['summer solstice'] = find_azimuths(observer)

    observer.date = ephem.next_equinox(observer.date)
    azimuths['autumnal equinox'] = find_azimuths(observer)

    observer.date = ephem.next_solstice(observer.date)
    azimuths['winter solstice'] = find_azimuths(observer)

    # pprint(azimuths)

    # How many degrees is close enough?
    DEGREESLOP = 2.

    # If allpoints is set, check angles among all pairs of points.
    # Otherwise, only check angles from observer to other points.
    if allpoints:
        print("Looking for alignments among all points")
        observer_points = waypoints
    else:
        observer_points = [ [ observer.name,
                              observer.lat / ephem.degree,
                              observer.lon / ephem.degree,
                              observer.elevation ] ]

    # Now go through all the angles between waypoints and see if
    # any of them correspond to any of the astronomical angles.
    matches = []
    for wp1 in observer_points:
        print("\nChecking observer", wp1)
        for wp2 in waypoints:
            if wp1 == wp2:
                continue
            angle = bearing_to(wp1, wp2)
            print("  ... vs", wp2, angle)

            # Does that angle match any of our astronomical ones?
            for season in azimuths:  # vernal equinox, etc.
                for body in azimuths[season]:  # sun, full moon
                    for event in azimuths[season][body]:  # rise, set
                        event_az = azimuths[season][body][event]['az']
                        if abs(event_az - angle) < DEGREESLOP:
                            matches.append({
                                'observer': wp1[0],
                                'target': wp2[0],
                                'event': '%s %s%s' % (season, body, event),
                                'azimuth': event_az,
                                'slop': event_az - angle,
                                'time': azimuths[season][body][event]['time'],
                                'latitude': wp2[1],
                                'longitude': wp2[2],
                            })

    return matches


def read_waypoint_file_CSV(filename):
    """Read a CSV waypoint file. Ignore tracks.
       Return a list of [name, lat, lon, ele] floats for waypoints.
    """
    import csv

    points = []
    observer = None
    with open(filename) as csvfp:
        reader = csv.DictReader(csvfp)
        for row in reader:
            # Each row is an OrderedDict
            try:
                if 'elevation' in row:
                    ele = float(row['elevation'])
                else:
                    ele = 0
                point = [ row['name'],
                          float(row['latitude']), float(row['longitude']),
                          ele ]
            except ValueError:
                print("ValueError on", row)
            if row['name'].lower() == "observer":
                observer = point
            else:
                points.append(point)

    return observer, points


def read_waypoint_file_GPX(filename):
    """Read a GPX waypoint file. Ignore tracks.
       Return a list of [name, lat, lon, ele] floats for waypoints.
    """

    dom = xml.dom.minidom.parse(filename)
    first_segment_name = None
    observer = None

    # Handle waypoints
    waypts = dom.getElementsByTagName("wpt")
    if not waypts:
        return []

    waypoints = []
    pointno = 0
    for pt in waypts:
        lat = float(pt.getAttribute("lat"))
        lon = float(pt.getAttribute("lon"))
        try:
            ele = float(get_DOM_text(pt, "ele"))
        except:
            ele = 500    # meters

        # Some of these names have newlines in them.
        # Just keep the part before the newline.
        name = get_DOM_text(pt, "name").split('\n')[0]

        # print("  { 'name': '%s', 'lat': %f, 'lon': %f }," % (name, lat, lon))
        if not name:
            pointno += 1
            name = "Point %d" % pointno

        if name.lower() == "observer":
            observer = [ 'Observer', lat, lon, ele ]
        else:
            waypoints.append([ name, lat, lon, ele ])

    # pprint(waypoints)
    return observer, waypoints


def get_DOM_text(node, childname=None):
    """Get the text out of a DOM node.
       Or, if childname is specified, get the text out of a child
       node with node name childname.
    """
    if childname:
        nodes = node.getElementsByTagName(childname)
        # print "node has", len(nodes), childname, "children"
        if not nodes:
            return None
        node = nodes[0]
    if not node:
        return None
    n = node.childNodes
    if len(n) >= 1 and n[0].nodeType == n[0].TEXT_NODE:
        return n[0].data
    return None


def find_observer_point(obsname, waypoints):
    """Find the waypoint with the name obsname."""
    for wp in waypoints:
        if wp[0] == obsname:
            return wp
    return None


def save_alignments_as_JSON(observer, alignments, waypoints, filename):
    """Given a list of alignments,
       save a JSON file that can be plotted in various ways.
    """
    print("Saving as JSON")
    out = []
    for a in alignments:
        if observer:
            obsname = observer.name
            obslat = math.degrees(observer.lat)
            obslon = math.degrees(observer.lon)
        else:
            obs = find_observer_point(a['observer'], waypoints)
            if not obs:
                print("Can't find observer for", a)
                continue
            obsname = obs[0]
            obslat = obs[1]
            obslon = obs[2]

        if (a['latitude'] == obslat and
            a['longitude'] == obslon):
            continue

        out.append({ 'observer_name': obsname,
                     'event': a['event'],
                     'observer_lat': math.degrees(obslat),
                     'observer_lon': math.degrees(obslon),
                     'target_name': a['target'],
                     'target_lat': a['latitude'],
                     'target_lon': a['longitude'],
                    })

    with open(filename, 'w') as outfp:
        outfp.write(json.dumps(out))


def save_alignments_as_GPX(observer, alignments, waypoints, filename):
    """Given a list of alignments,
       save them as a GPX file with tracks between observer and each target,
       and a waypoint for each target.
    """
    with open(filename, 'w') as outfp:
        print('''<?xml version='1.0' encoding='UTF-8' standalone='yes' ?>
<gpx version="1.1" creator="SkyAlignments~" xmlns="http://www.topografix.com/GPX/1/1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
<trk>''', file=outfp)

        for a in alignments:
            if observer:
                obsname = observer.name
                obslat = math.degrees(observer.lat)
                obslon = math.degrees(observer.lon)
            else:
                obs = find_observer_point(a['observer'], waypoints)
                if not obs:
                    print("Can't find observer for", a)
                    continue
                obsname = obs[0]
                obslat = obs[1]
                obslon = obs[2]

            if (a['latitude'] == obslat and
                a['longitude'] == obslon):
                continue
            print('''  <trkseg>
    <name>%s, %s</name>
    <trkpt lat="%f" lon="%f">
      <time>%s</time>
      <name>%s</name>
    </trkpt>
    <trkpt lat="%f" lon="%f">
      <time>%s</time>
      <name>%s, %s</name>
    </trkpt>
  </trkseg>''' % (a['target'], a['event'],
                  obslat, obslon,
                 str(a['time']), obsname,
                   a['latitude'], a['longitude'], str(a['time']),
                   a['target'], a['event']),
                  file=outfp)

        print("</trk>", file=outfp)

        # Apparently naming trksegs doesn't do anything, at least in pytopo.
        # Set waypoints too.
        for a in alignments:
            if (a['latitude'] == obslat and
                a['longitude'] == obslon):
                continue
            print('''  <wpt lat="%f" lon="%f">
    <time>%s</time>
    <name>%s, %s</name>
  </wpt>''' % (a['latitude'], a['longitude'], str(a['time']),
               a['target'], a['event']),
                  file=outfp)

        print("</gpx>", file=outfp)
        print("Saved alignments to", filename)

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=""
        """Find alignments between latitude/longitude coordinate pairs
and the sun. moon, and other objects on special dates such as
solstices and equinoxes.

Observer location may be specified either with -o lat,lon,ele or by
naming one of the GPX waypoints 'Observer'; otherwise the first
waypoint in the first file will be taken as the observer location.

When specifying location on the command line, latitude and longitude
are in decimal degrees. Elevation is optional; it will be assumed to be
meters unless followed by the letter f,
e.g. -o 34.8086585,-103.2011914,1650f""",
                                     formatter_class=argparse.RawTextHelpFormatter)
    # Specify an observer:
    parser.add_argument('-o', '--observer', action="store", dest="observer",
                        help='Observer location (lat,lon[,ele])')

    # Specify an observer name:
    parser.add_argument('-n', '--observername', action="store",
                        dest="observername", help='Observer name')

    # Different year from now:
    parser.add_argument('-y', '--year', action="store", type=int,
                        dest="year", help='Year')

    # Output JSON file:
    parser.add_argument('-j', '--writejson', action="store",
                        dest="writejson", help='Save JSON to output file')

    # Output GPX file:
    parser.add_argument('-g', '--writegpx', action="store",
                        dest="writegpx", help='Save GPX to output file')

    # Don't use an observer, check angles between all pairs of points:
    parser.add_argument('-a', "--all", dest="allpoints", default=False,
                        action="store_true",
                        help="Don't use an observer, check angles "
                              "between all pairs of points")

    parser.add_argument('waypointfiles', nargs='+',
                        help='GPX or CSV files containing waypoints')
    args = parser.parse_args(sys.argv[1:])

    if args.observer:
        floats = args.observer.split(',')
        lat = float(floats[0].strip())
        lon = float(floats[1].strip())
        if len(floats) > 2:
            if floats[2].endswith('f'):    # ends with f, convert feet to meters
                ele = float(floats[2][:-1].strip()) * 0.3048
            elif floats[2].endswith('m'):  # ends with m, already meters
                ele = float(floats[2][:-1].strip())
            else:                          # assume meters
                ele = float(floats[2].strip())
        else:
            ele = 0.
        observer_point = [ 'Observer', lat, lon, ele ]
    else:
        observer_point = None

    if args.year:
        year = args.year
    else:
        year = datetime.now().year


    waypoints = []
    for filename in args.waypointfiles:
        if filename.lower().endswith('gpx'):
            obs, wp = read_waypoint_file_GPX(filename)
        elif filename.lower().endswith('csv'):
            obs, wp = read_waypoint_file_CSV(filename)
        else:
            print("Unknown file type %s: skipping" % filename)
            continue
        if wp:
            waypoints += wp
        else:
            print("No waypoints in", filename)
        if obs:
            observer_point = obs

    if not waypoints:
        parser.print_help()
        sys.exit(1)

    if not observer_point:
        print("Using first waypoint for observer:", waypoints[0])
        observer_point = waypoints[0]
        waypoints = waypoints[1:]

    # find_alignments will need an initial observer to calculate
    # things like the dates of equinoxes and solstices, and the
    # azimuths. So set it even if allpoints is set.
    observer = ephem.Observer()
    # Observer will take degrees as a string, but if you pass it floats
    # it expects radians, though that's undocumented.
    observer.lat = observer_point[1] * ephem.degree
    observer.lon = observer_point[2] * ephem.degree
    if len(observer_point) > 3:
        observer.elevation = observer_point[3]
    else:
        observer.elevation = 500.0  # meters
    if args.observername:
        observer.name = args.observername
    else:
        observer.name = "%s %f, %f, %dm" % (observer_point[0],
                                            observer.lat / ephem.degree,
                                            observer.lon / ephem.degree,
                                            observer.elevation)
    # print(observer)

    alignments = find_alignments(observer, waypoints,
                                 year=year, allpoints=args.allpoints)

    # Now we no longer need the ephem.observer, and if allpoints is set,
    # don't want to pass observer to save_alignments*.
    if args.allpoints:
        observer = None

    if alignments:
        # pprint(alignments)
        cur_observer = None
        for a in alignments:
            if a['observer'] != cur_observer:
                cur_observer = a['observer']
                print("\nFrom %s:" % cur_observer)

            print("%s, %s at %s (az %d +/- %.2f)"
                  % ( a['event'], a['target'],
                      a['time'],
                      a['azimuth'], a['slop']))

        if args.writejson:
            save_alignments_as_JSON(observer, alignments, waypoints,
                                    args.writejson)

        if args.writegpx:
            save_alignments_as_GPX(observer, alignments, waypoints,
                                   args.writegpx)

    else:
        print("Couldn't find any alignments with %s" % observer.name)

