try:
    from . import convlib
    from . import matlabport
    from . import skyplot
except ImportError:
    import convlib
    import matlabport
    import skyplot
import numpy as np
import math
import matplotlib.pyplot as plt
import pathlib

SCRIPTDIR = pathlib.Path(__file__).parent

LAT, LON = math.radians(48.775556), math.radians(9.182778)
# LAT, LON = math.radians(52.52), math.radians(13.405)
# LAT, LON = 0, 0


elevation_values = []
azimuth_values = []
hour_values = []

for hour in range(0, 2400):
    hour /= 100
    hour -= 2

    t, _, _ = matlabport.julianjh(2025, 7, 19, hour)

    GAST, _ = matlabport.jul2gast(hour, t)

    GAST_rad = math.radians(GAST * 15)


    aapp,dapp,r = matlabport.solarpos(t)

    aapp = math.radians(aapp)
    dapp = math.radians(dapp)

    epsilon, delta_epsilon, delta_psi = matlabport.nutwink(t)
    z, Theta, zeta_0 = matlabport.prezwink(t)
    epsilon = math.radians(epsilon)
    delta_epsilon = math.radians(delta_epsilon)
    delta_psi = math.radians(delta_psi)
    z = math.radians(z)
    Theta = math.radians(Theta)
    zeta_0 = math.radians(zeta_0)

    r_i0 = np.array([math.cos(dapp) * math.cos(aapp), math.cos(dapp) * math.sin(aapp), math.sin(dapp)])


    N = convlib.get_N(epsilon, delta_epsilon, delta_psi)
    P = convlib.get_P(z, Theta, zeta_0)
    r_g = convlib.get_r_g(LAT, LON, GAST_rad, N, P, r_i0)

    if False:
        print(hour)
        print(t)
        print(GAST)
        print(aapp)
        print(dapp)
        print(N)
        print(P)
        print(r_i0)
        print(r_g)
    #print(r_g)

    x, y, z = r_g

    elev, azim = convlib.elev_azim(x, y, z)

    elevation_values.append(elev)
    azimuth_values.append(azim)
    hour_values.append(hour)


fig, ax = plt.subplots()
for elev, azim in zip(elevation_values, azimuth_values):
    ax.plot(elevation_values, azimuth_values)
ax.axhline(0)
ax.axhline(-6)

fig = skyplot.skyplot(plt, elevation_values, azimuth_values)

fig.savefig(SCRIPTDIR / ".." / "output" / "fig5.svg")

if __name__ == "__main__":
    plt.show()