// Encoding: UTF-8
#include "calendar.h"
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <iomanip>

// 常量
#define day_secs 86400
#define delta_time 1e-11

/* calendar.h规定之量 */
Julian julian;
parameter p;
double epsilon;

FILE *fp;

pair<double, double> getSunPos(double t, double flag = false) {
    pair<double, double> sun = p.sun_longitude(t);

    double alpha, delta;                          // 太阳赤经赤纬
    double lambda = sun.first, beta = sun.second; // 太阳黄经黄纬
    // 黄赤交角epsilon
    epsilon = p.get_epsilon(t);

    delta =
        asin(sin(epsilon) * sin(lambda) * cos(beta) + cos(epsilon) * sin(beta));
    alpha = atan2(
        (cos(epsilon) * sin(lambda) * cos(beta) - sin(epsilon) * sin(beta)),
        (cos(lambda) * cos(beta)));
    // 时角
    double H = p.getHourAngle(t, LONGITUDE, LATITUDE, alpha);
    double A, h; //方位角和高度角
    h = asin(sin(LATITUDE) * sin(delta) + cos(LATITUDE) * cos(delta) * cos(H));
    A = atan2(-cos(delta) * sin(H),
              cos(LATITUDE) * sin(delta) - sin(LATITUDE) * cos(delta) * cos(H));
    if (flag) {
        fprintf(
            fp, "%s\t%s\t%s\t%s\t%s\t%s\n", radToDMS(0, 360, lambda).c_str(),
            radToDMS(-90, 90, beta).c_str(), radToDMS(0, 360, alpha).c_str(),
            radToDMS(-90, 90, delta).c_str(), radToDMS(0, 360, A).c_str(),
            radToDMS(-90, 90, h).c_str());
    }
    return make_pair(A, h + M_PI / 216);
}

void binarySearch(string message, time_t begin, Date date) {
    double t = julian.getJulianKiloYear(begin);
    double t_end = julian.getJulianKiloYear(begin + day_secs / 2);
    double middle;
    time_t res;
    pair<double, double> sunPos;
    while (t_end - t >= delta_time) {
        middle = (t + t_end) / 2;
        sunPos = getSunPos(middle, false);
        if (sunPos.second * getSunPos(t).second > 0) {
            t = middle;
        } else {
            t_end = middle;
        }
    }
    res = julian.kiloYearToTime(middle, date.tm_year + 1900);
    Date *tm = localtime(&res);
    printf("%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n", message.c_str(),
           tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday, tm->tm_hour,
           tm->tm_min, tm->tm_sec, radToDMS(0, 360, sunPos.first).c_str());
    fprintf(fp, "%s\t%4d-%02d-%02d %02d:%02d:%02d\tAzimuth\t%s\n",
            message.c_str(), tm->tm_year + 1900, tm->tm_mon + 1, tm->tm_mday,
            tm->tm_hour, tm->tm_min, tm->tm_sec,
            radToDMS(0, 360, sunPos.first).c_str());
}

void getSunTime(Date date) {
    time_t time = mktime(&date);
    binarySearch("Sunrize", time, date);
    binarySearch("Sunset", time + day_secs / 2, date);
}

void showSunPos(Date date) {
    time_t time = mktime(&date);
    double t;
    fprintf(fp, "\n\nSun Trace Table\n\n");
    fprintf(fp, "--------------------------------------------------------------"
                "--------------------------------------------------------------"
                "----------\n");
    fprintf(fp, "Time     |\tlongitude\t\tlatitude\t\tright "
                "ascension\tdeclination\tAzimuth\t\tAltitude\n");
    fprintf(fp, "--------------------------------------------------------------"
                "--------------------------------------------------------------"
                "----------\n");
    for (int i = 0; i < 24; i++, time += 3600, date.tm_hour++) {
        t = julian.getJulianKiloYear(time);
        fprintf(fp, "%02d:%02d:%02d |\t", date.tm_hour, date.tm_min,
                date.tm_sec);
        getSunPos(t, true);
    }
    fprintf(fp, "--------------------------------------------------------------"
                "--------------------------------------------------------------"
                "----------\n");
}

int main() {
    Date date;
    cout << "Input a day you want(in YYYY-MM-DD): ";
    scanf("%d-%d-%d", &date.tm_year, &date.tm_mon, &date.tm_mday);
    date.tm_year -= 1900;
    date.tm_mon -= 1;
    date.tm_hour = date.tm_min = date.tm_sec = 0;

    fp = fopen("sunPos.txt", "w");
    if (fp == NULL) {
        printf("File open failed\n");
        exit(1);
    }
    fprintf(fp, "Sun trace File (Date: %4d-%02d-%02d)\n\n\n",
            date.tm_year + 1900, date.tm_mon + 1, date.tm_mday);
    getSunTime(date);

    showSunPos(date);
    fclose(fp);

    return 0;
}