#include <stdio.h>
#include <string.h>
#include <stdbool.h>
#include <math.h>
#include <time.h>

#include "../include/sun.h"

double deg_to_rad(double angle) {
    return angle * (M_PI / 180);
}

double rad_to_deg(double angle) {
    return angle * (180 / M_PI);
}

double sun_calculation(struct tm curr_time, double longitude, double latitude) {
    // Calulate day of the year
    int day_in_year = floor(275 * (curr_time.tm_mon + 1) / 9) - 
                        ((floor(((curr_time.tm_mon + 1) + 9) / 12)) * 
                        (1 + floor(((curr_time.tm_year + 1900) - 4 * floor((curr_time.tm_year + 1900) / 4) + 2) / 3))) + 
                        curr_time.tm_mday - 30;

    // Longitude to hour, get approx time
    double long_hour = longitude / 15;
    double approx_time;

    if (curr_time.tm_hour < 12)
        approx_time = day_in_year + ((6 - long_hour) / 24);
    else
        approx_time = day_in_year + ((18 - long_hour) / 24);

    // Sun's mean anomaly
    double mean_anomaly = (0.9856 * approx_time) - 3.289;

    // Sun's true longitude
    double true_longitude = mean_anomaly + (1.916 * sin(deg_to_rad(mean_anomaly)) + 
                            (0.02 * sin(deg_to_rad(2 * mean_anomaly)) + 282.634));

    // Keep longitude between [0, 360)
    if (true_longitude >= 360)
        true_longitude -= 360;
    else if (true_longitude < 0)
        true_longitude += 360;
    
    // Sun's right ascension
    double right_ascension = rad_to_deg(atan(0.91764 * tan(deg_to_rad(true_longitude))));

    // Keep Right Ascension between [0, 360)
    if (right_ascension >= 360)
        right_ascension -= 360;
    else if (right_ascension < 0)
        right_ascension += 360;

    // Quadrant calibration
    int longitude_quadrant = (floor(true_longitude / 90)) * 90;
    int right_ascension_quadrant = (floor(right_ascension / 90)) * 90;

    right_ascension += (longitude_quadrant - right_ascension_quadrant);

    right_ascension /= 15; // convert to hours

    // Calculate Sun's declination
    double sin_declination = 0.39782 * sin(deg_to_rad(true_longitude));
    double cos_declination = cos(asin(sin_declination));

    // Calculate Sun's local hour angle
    double sun_hour;

    double cos_hour = (-0.01454 - (sin_declination * sin(deg_to_rad(latitude)))) / 
                (cos_declination * cos(deg_to_rad(latitude)));
    
    if (curr_time.tm_hour < 12)
        sun_hour = 360 - rad_to_deg(acos(cos_hour));
    else
        sun_hour = rad_to_deg(acos(cos_hour));

    sun_hour /= 15; // convert to hours

    // Local time
    double local_mean_time = sun_hour + right_ascension - (0.06571 * approx_time) - 6.622;

    // Adjust to UTC
    double sun_UTC = local_mean_time - long_hour;

    // Keep Sun UTC between [0, 24)
    if (sun_UTC >= 24)
        sun_UTC -= 24;
    else if (sun_UTC < 0)
        sun_UTC += 24;

    double sun_time = sun_UTC + LOCAL_OFFSET;
    
    if (sun_time < 0)
        return (24 + sun_time);
    else if (sun_time > 24)
        return (sun_time - 24);
    else    
        return sun_time;
}