#!/usr/bin/env python3
# ---------------------------------------------------------------------------- #
## \file sundial3D.py
## \author Sebastien Beaugrand
## \sa http://beaugrand.chez.com/
## \copyright CeCILL 2.1 Free Software license
# ---------------------------------------------------------------------------- #
import sys
import argparse
import math as m
from abc import ABC, abstractmethod
from datetime import datetime
from threading import Thread
from pymeeus.Angle import *
from pymeeus.Epoch import Epoch
from pymeeus.Coordinates import true_obliquity, nutation_longitude, equatorial2horizontal
from pymeeus.Sun import Sun
from solid import *

MARGIN = 0.1
RSTYLE = 5

parser = argparse.ArgumentParser(
    epilog='''
example: pip3 install pymeeus
         pip3 install solidpython
         ./sundial3D.py -s sphere     --draft -o ~/sundial3Ds.scad
         ./sundial3D.py -s polyhedron --draft -o ~/sundial3Dp.scad
         ./sundial3D.py -s wall       --draft -o ~/sundial3Dw.scad
         openscad ~/sundial3Dw.scad &
''',
    formatter_class=argparse.RawDescriptionHelpFormatter)
parser.add_argument('--lat',
                    type=float,
                    default=48.44728,
                    help='default: 48.44728 (deg)')
parser.add_argument('--lon',
                    type=float,
                    default=1.48749,
                    help='default: 1.48749 (deg)')
parser.add_argument('-s',
                    '--shape',
                    default='wall',
                    choices=['sphere', 'polyhedron', 'wall'],
                    help='default: wall')
parser.add_argument('-R',
                    '--R',
                    type=int,
                    default=100,
                    help='R sphere, default: 100 (mm)')
parser.add_argument('-r',
                    '--r',
                    type=int,
                    default=50,
                    help='r sphere, default: 50 (mm)')
parser.add_argument('--draft', action='store_true')
parser.add_argument('-o', '--output')
parser.add_argument('--debug', action='store_true')
parser.add_argument('-y', '--year', type=int, default=datetime.today().year)
parser.add_argument('-m', '--month', type=int)
parser.add_argument('-d', '--day', type=int)
parser.add_argument('-H', '--hour', type=int)
args = parser.parse_args()


# ---------------------------------------------------------------------------- #
## \fn debugf
# ---------------------------------------------------------------------------- #
def debugf(dic):
    if not args.debug:
        return
    print('debug: ', end='', file=sys.stderr)
    print(', '.join('{} = {:.7f}'.format(k, v) for k, v in dic.items()),
          file=sys.stderr)


# ---------------------------------------------------------------------------- #
## \fn azi_ele
# ---------------------------------------------------------------------------- #
def azi_ele(year, month, day, hour):
    day += hour / 24
    epoch = Epoch(year, month, day)

    ra, dec, _ = Sun.apparent_rightascension_declination_coarse(epoch)
    debugf({'ra': ra._deg, 'dec': dec._deg})

    theta0 = epoch.apparent_sidereal_time(true_obliquity(year, month, day),
                                          nutation_longitude(year, month, day))
    theta0 *= 24 * 15
    azi, ele = equatorial2horizontal(
        Angle(Angle.reduce_deg(theta0 + args.lon - ra)), dec, Angle(args.lat))
    debugf({'theta0': theta0, 'azi': azi._deg, 'ele': ele._deg})

    if ele < 0:
        return None

    azi = m.radians(azi._deg)
    ele = m.radians(ele._deg)
    return (azi, ele)


# ---------------------------------------------------------------------------- #
## \class Triangle
# ---------------------------------------------------------------------------- #
class Triangle:
    def __init__(self, a, b, c):
        self.a = a
        self.b = b
        self.c = c
        self.v0 = (b[0] - a[0], b[1] - a[1], b[2] - a[2])
        self.v1 = (c[0] - a[0], c[1] - a[1], c[2] - a[2])
        self.n1 = self.v0[1] * self.v1[2] - self.v0[2] * self.v1[1]
        self.n2 = self.v0[2] * self.v1[0] - self.v0[0] * self.v1[2]
        self.n3 = self.v0[0] * self.v1[1] - self.v0[1] * self.v1[0]
        n = (self.n1 * self.n1 + self.n2 * self.n2 + self.n3 * self.n3)**0.5
        self.n1 /= n
        self.n2 /= n
        self.n3 /= n
        self.d = -self.n1 * a[0] - self.n2 * a[1] - self.n3 * a[2]

    def contains(self, v2):
        v0 = self.v0
        v1 = self.v1
        # Dot products
        dot00 = v0[0] * v0[0] + v0[1] * v0[1] + v0[2] * v0[2]
        dot01 = v0[0] * v1[0] + v0[1] * v1[1] + v0[2] * v1[2]
        dot02 = v0[0] * v2[0] + v0[1] * v2[1] + v0[2] * v2[2]
        dot11 = v1[0] * v1[0] + v1[1] * v1[1] + v1[2] * v1[2]
        dot12 = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]
        # Barycentric coordinates
        invDenom = 1 / (dot00 * dot11 - dot01 * dot01)
        u = (dot11 * dot02 - dot01 * dot12) * invDenom
        v = (dot00 * dot12 - dot01 * dot02) * invDenom

        return u >= 0 and v >= 0 and u + v < 1

    def intersect(self, u, v, w):
        t = -self.d / (self.n1 * u + self.n2 * v + self.n3 * w)
        x = u * t
        y = v * t
        z = w * t
        v2 = (x - self.a[0], y - self.a[1], z - self.a[2])
        if self.contains(v2):
            return (x, y, z, self.n1, self.n2, self.n3)
        else:
            return None


# ---------------------------------------------------------------------------- #
## \class Quadrilateral
# ---------------------------------------------------------------------------- #
class Quadrilateral:
    def __init__(self, a, b, c, d):
        self.a = a
        self.b = b
        self.c = c
        self.d = d
        self.t = list()
        self.t.append(Triangle(a, b, c))
        self.t.append(Triangle(a, c, d))

    def contains(self, p):
        for t in self.t:
            v2 = (p[0] - t.a[0], p[1] - t.a[1], p[2] - t.a[2])
            if t.contains(v2):
                return (p[0], p[1], p[2], t.n1, t.n2, t.n3)
        return None

    def intersect(self, u, v, w):
        for t in self.t:
            p = t.intersect(u, v, w)
            if p is not None:
                return p
        return None

    def polyhedron(self):
        tr = self.t[0]
        points = [
            self.a, self.b, self.c, self.d,
            [self.a[0] - tr.n1, self.a[1] - tr.n2, self.a[2] - tr.n3],
            [self.b[0] - tr.n1, self.b[1] - tr.n2, self.b[2] - tr.n3],
            [self.c[0] - tr.n1, self.c[1] - tr.n2, self.c[2] - tr.n3],
            [self.d[0] - tr.n1, self.d[1] - tr.n2, self.d[2] - tr.n3]
        ]
        faces = [[0, 1, 2, 3], [4, 5, 1, 0], [7, 6, 5, 4], [5, 6, 2, 1],
                 [6, 7, 3, 2], [7, 4, 0, 3]]
        return polyhedron(points=points, faces=faces)


# ---------------------------------------------------------------------------- #
## \class Shape
# ---------------------------------------------------------------------------- #
class Shape(ABC):
    @abstractmethod
    def intersect(self, year, month, day, hour):
        pass

    def styledroit(self):
        return translate([0, 0,
                          -args.r - MARGIN])(cylinder(h=args.r + MARGIN * 2,
                                                      r1=RSTYLE,
                                                      r2=0,
                                                      segments=180))


# ---------------------------------------------------------------------------- #
## \class Sphere
# ---------------------------------------------------------------------------- #
class Sphere(Shape):
    def __init__(self):
        pass

    def intersect(self, year, month, day, hour):
        c = azi_ele(year, month, day, hour)
        if c is None:
            return None
        azi, ele = c

        theta = m.atan(args.r * m.tan(azi) / args.R)
        if azi > -m.pi / 2 and azi < m.pi / 2:
            theta += m.pi

        delta = m.atan(args.R * m.tan(ele) * m.sin(azi) / m.sin(theta) /
                       args.r)

        y = -args.R * m.cos(delta) * m.cos(theta)
        x = -args.R * m.cos(delta) * m.sin(theta)
        z = args.r * m.sin(delta)
        return (x, y, z, 0, 0, -1)

    def normal_vector(self, p):
        return (p[0], p[1], p[2], 0, 0, -1)

    def draw(self):
        s = self.styledroit()
        demisphere = s
        demisphere += resize(newsize=[args.R * 2, args.R * 2, args.r * 2])(
            sphere(r=args.r, segments=180))
        demisphere -= translate([-args.R, -args.R, 0])(cube(
            [args.R * 2, args.R * 2, args.r + MARGIN]))
        demisphere -= s
        return color('white', 0.2)(demisphere)


# ---------------------------------------------------------------------------- #
## \class Polyhedron
# ---------------------------------------------------------------------------- #
class Polyhedron(Shape):
    def __init__(self):
        self.quadrilaterals = list()
        self.quadrilaterals.append(
            Quadrilateral((-args.R, RSTYLE, -args.r), (-args.R, args.R, 0),
                          (args.R, args.R, 0), (args.R, RSTYLE, -args.r)))
        a = self.intersect(args.year, 6, 21, 8.5)
        b = self.intersect(args.year, 12, 21, 8.5)
        if a[0] is None or b[0] is None:
            print(a)
            print(b)
            exit(1)
        a = (a[0], RSTYLE, -args.r)
        b = (b[0], args.R, 0)
        c = (-b[0], b[1], b[2])
        d = (-a[0], a[1], a[2])
        e = (b[0], -b[1], b[2])
        f = (a[0], -a[1], a[2])
        self.quadrilaterals = list()
        self.quadrilaterals.append(Quadrilateral(a, b, c, d))
        self.quadrilaterals.append(Quadrilateral(f, e, b, a))
        if not args.draft:
            g = (-b[0], -b[1], b[2])
            h = (-a[0], -a[1], a[2])
            self.quadrilaterals.append(Quadrilateral(d, c, g, h))

    def normal_vector(self, point):
        for q in self.quadrilaterals:
            p = q.contains(point)
            if p is not None:
                return p
        return None

    def intersect(self, year, month, day, hour):
        c = azi_ele(year, month, day, hour)
        if c is None:
            return None
        azi, ele = c

        azi = m.pi * 3 / 2 - azi
        u = m.cos(ele) * m.cos(azi)
        v = m.cos(ele) * m.sin(azi)
        w = m.sin(ele)

        for q in self.quadrilaterals:
            p = q.intersect(u, v, w)
            if p is not None:
                return p
        return None

    def draw(self):
        o = color('white', 0.2)(self.styledroit())
        for q in self.quadrilaterals:
            o += color('white', 0.2)(q.polyhedron())
        return o


# ---------------------------------------------------------------------------- #
## \class Wall
# ---------------------------------------------------------------------------- #
class Wall(Polyhedron):
    def __init__(self):
        a = (-300, 65, -280)
        b = (-300, 65, 0)
        c = (300, 65, 0)
        d = (300, 65, -280)
        e = (-300, -180, 0)
        f = (-300, -180, -280)
        self.quadrilaterals = list()
        self.quadrilaterals.append(Quadrilateral(a, b, c, d))
        self.quadrilaterals.append(Quadrilateral(f, e, b, a))
        if not args.draft:
            g = (300, -180, 0)
            h = (300, -180, -280)
            self.quadrilaterals.append(Quadrilateral(d, c, g, h))


# ---------------------------------------------------------------------------- #
## \class Analemme
# ---------------------------------------------------------------------------- #
class Analemme(Thread):
    WIDTH = 0.5

    def __init__(self, shape, year, hour, minu=0):
        Thread.__init__(self)
        self.shape = shape
        self.year = year
        self.hour = hour + minu / 60
        self.minu = minu
        if minu == 0:
            self.depth = 2
        else:
            self.depth = 1
        self.objects = union()
        if type(self.shape) is Wall:
            self.min = (0, 0, 0)
            self.max = (0, 0, -args.r)
        else:
            self.min = (0, args.R, 0)
            self.max = (0, -args.R, 0)

    def number(self, n):
        return linear_extrude(height=self.depth)(text(str(n),
                                                      font='DejaVu Sans',
                                                      size=3))

    def draw_numbers(self):
        rx = 0
        ry = 0
        rz = 0
        p = self.shape.normal_vector(self.min)
        if p is not None:
            x, y, z, n1, n2, n3 = p
            if n3 == -1:
                pass
            elif n1 == 0:
                rx = 180 - m.degrees(m.atan2(n2, n3))
            elif n2 == 0:
                ry = 180 + m.degrees(m.atan2(n1, n3))
        self.objects += translate([self.min[0], self.min[1], self.min[2]])(
            rotate([rx, ry,
                    rz])(translate([0, -5, -self.depth
                                    ])(self.number(int(self.hour) + 2))))
        if self.max[2] > -self.depth:
            return
        rx = 0
        ry = 0
        rz = 0
        p = self.shape.normal_vector(self.min)
        if p is not None:
            x, y, z, n1, n2, n3 = p
            if n3 == -1:
                pass
            elif n1 == 0:
                rx = 180 - m.degrees(m.atan2(n2, n3))
            elif n2 == 0:
                ry = 180 + m.degrees(m.atan2(n1, n3))
        self.objects += translate([self.max[0], self.max[1], self.max[2]])(
            rotate([rx, ry,
                    rz])(translate([0, 5, -self.depth
                                    ])(self.number(int(self.hour) + 1))))

    def run(self):
        for dnum in range(1, 366):
            date = datetime.strptime('{:04d}{:03d}'.format(self.year, dnum),
                                     '%Y%j')
            if date.month > 6:
                col = 'red'
            else:
                col = 'blue'
            rx = 0
            ry = 0
            rz = 0
            p = self.shape.intersect(self.year, date.month, date.day,
                                     self.hour)
            if p is None:
                continue
            x, y, z, n1, n2, n3 = p
            if n3 == -1:
                pass
            elif n1 == 0:
                rx = -m.degrees(m.atan2(n2, n3))
            elif n2 == 0:
                ry = m.degrees(m.atan2(n1, n3))
            if type(self.shape) is Wall:
                if self.min[2] > z:
                    self.min = (x, y, z)
                if self.max[2] < z:
                    self.max = (x, y, z)
            else:
                if self.min[1] > y:
                    self.min = (x, y, z)
                if self.max[1] < y:
                    self.max = (x, y, z)
            self.objects += color(col)(translate([x, y, z])(rotate(
                [rx, ry,
                 rz])(cube([self.WIDTH, self.WIDTH, self.depth + MARGIN]))))
        if self.minu == 0:
            self.draw_numbers()


# ---------------------------------------------------------------------------- #
# main
# ---------------------------------------------------------------------------- #
if args.debug:
    azi_ele(args.year, args.month, args.day, args.hour)
    exit(0)

if args.shape == 'sphere':
    shape = Sphere()
elif args.shape == 'polyhedron':
    shape = Polyhedron()
elif args.shape == 'wall':
    shape = Wall()

objects = union()

threads = list()
for hour in range(5, 13 if args.draft else 20):
    t = Analemme(shape, args.year, hour)
    t.start()
    threads.append(t)
    if not args.draft and hour < 19:
        t = Analemme(shape, args.year, hour, minu=30)
        t.start()
        threads.append(t)
for t in threads:
    t.join()
    objects += t.objects

objects += shape.draw()
scad_render_to_file(objects, args.output)