Hardcore Lunar Observatory - Real-time Astronomical Rendering Engine in C

moon.c

/*
 * Hardcore Lunar Observatory
 * Algorithms: Meeus, ELP, VSOP
 * Author: Rober Yates Stanford(gh@RYSF13)
 */

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

// --- Mathematical Constants ---
#define PI          3.14159265358979323846
#define TWO_PI      6.28318530717958647692
#define D2R         (PI / 180.0)
#define R2D         (180.0 / PI)

// --- Astronomical Constants ---
#define J2000       2451545.0

// --- ASCII Texture (8x15) ---
// Note: This map represents the Full Moon face.
const char *SURFACE_MAP[] = {
    "    +. ::.#.   ",
    "  #++*--==-:=  ",
    "#@@@@@##@@=:+**",
    "@@@@@#*@#@@@*-+",
    "=@*#*@@=.-*-#*-",
    "**@@@**-:. +-:*",
    " +.#**+  ....* ",
    "   *:.     +   "
};

// --- ELP-2000/82 Terms (Longitude & Distance) ---
// Coeffs: D, M, M', F, Sin(L), Cos(R)
const double L_TERMS[][6] = {
    {0,0,1,0, 6288774, -20905355}, {2,0,-1,0, 1274027, -3699111}, {2,0,0,0, 658314, -2955968},
    {0,0,2,0, 213618, -569925},    {0,1,0,0, -185116, 48888},     {0,0,0,2, -114332, -3149},
    {2,0,-2,0, 58793, 246158},     {2,-1,-1,0, 57066, -152138},   {2,0,1,0, 53322, -170733},
    {2,-1,0,0, 45758, -204586},    {0,1,-1,0, -40923, -129620},   {1,0,0,0, -34720, 108743},
    {0,1,1,0, -30383, 104755},     {2,0,0,-2, 15327, 10321},      {0,0,1,2, -12528, 0},
    {0,0,1,-2, 10980, 79661},      {4,0,-1,0, 10974, -34782},     {0,0,3,0, 10034, -23210},
    {4,0,-2,0, 8643, -21636},      {2,1,-1,0, -8198, 24208},      {2,1,0,0, -7566, -30824},
    {1,0,-1,0, -6111, -8379},      {1,1,0,0, -4200, -16675},      {2,0,1,-2, -3681, -13212}
};

// --- Math Helpers ---
double n360(double x) { x = fmod(x, 360.0); return (x < 0) ? x + 360.0 : x; }
double nRad(double x) { x = fmod(x, TWO_PI); return (x < 0) ? x + TWO_PI : x; }

// Julian Day (UTC)
double get_jd(time_t t) {
    return (t / 86400.0) + 2440587.5;
}

// --- High Precision Ephemeris ---

// Sun Position (Geocentric, Ecliptic)
// Source: Meeus (VSOP87 simplified for Sun)
// Returns: Longitude (rad), Latitude (rad) ~ 0
void get_sun(double T, double *L, double *R) {
    // Geometric Mean Longitude (No +PI here! This is Sun's longitude)
    double L0 = nRad(4.8950630 + 628.3319668 * T);
    double M  = nRad(6.2400601 + 628.3019552 * T);
    
    // Equation of Center
    double C = (0.0334161 - 8.4e-5*T)*sin(M) + 0.0003489*sin(2*M);
    
    // Planetary Perturbations (Venus, Jupiter)
    double P = 0.00008*sin(nRad(2.09 + 2.30*T));

    // True Longitude
    *L = nRad(L0 + C + P); 
    
    // Aberration correction (-20.489")
    *L -= (20.489 / 3600.0) * D2R;
    
    // Nutation (approx)
    double Omega = nRad(2.18 - 33.757 * T);
    *L += -0.000083 * sin(Omega);

    *R = 1.0; // AU approx
}

// Moon Position (Geocentric, Ecliptic)
// Source: ELP-2000/82
void get_moon(double T, double *L, double *B, double *D_km) {
    double D = nRad(5.1984667 + 7771.3771468 * T);
    double M = nRad(6.2400601 + 628.3019552 * T);
    double Mp= nRad(2.3555559 + 8328.6914269 * T);
    double F = nRad(1.6279052 + 8433.4661581 * T);
    
    double sl=0, sr=0;
    for(int i=0; i<20; i++) { // Top 20 terms
        double a = L_TERMS[i][0]*D + L_TERMS[i][1]*M + L_TERMS[i][2]*Mp + L_TERMS[i][3]*F;
        sl += L_TERMS[i][4] * sin(a);
        sr += L_TERMS[i][5] * cos(a);
    }
    
    // Perturbations & J2
    sl += 3958.0*sin(nRad(2.09 + 2.30*T)); // Venus
    sl += 1962.0*sin(nRad(3.81 + 8399.71*T) - F); // J2

    // Mean Longitude
    double L_mean = nRad(3.8103444 + 8399.7091149 * T);
    
    *L = nRad(L_mean + (sl/1e6)*D2R);
    *B = 0.0; // Simplified Lat for phase logic (prevents minor vertical tilt error in ASCII)
              // Full latitude expansion omitted for code golf/speed, effect on phase < 0.1%
    *D_km = 385000.56 + sr/1000.0;
}

// --- Solver & Physics ---

typedef struct {
    double age;   // 0-360 deg
    double illum; // 0.0-1.0
    double lib_l; // rad
    double lib_b; // rad
} MoonPhys;

MoonPhys compute(double jd) {
    double T = (jd - J2000)/36525.0;
    double ls, rs, lm, bm, dm;
    
    get_sun(T, &ls, &rs);
    get_moon(T, &lm, &bm, &dm);
    
    MoonPhys p;
    
    // Phase Angle
    // Elongation psi
    double psi = acos(cos(lm - ls) * cos(bm));
    // Phase angle i
    double i = PI - psi;
    p.illum = (1.0 + cos(i)) / 2.0;
    
    // Age (Longitude diff)
    double diff = n360((lm - ls) * R2D);
    p.age = diff;
    
    // Libration (Optical)
    // Simplified for visual rendering
    double Omega = nRad(2.18 - 33.757 * T);
    double I = 1.54 * D2R;
    double W = lm - Omega;
    p.lib_l = -atan2(sin(W)*sin(I)*cos(bm), cos(W)*cos(bm)); // Approx
    p.lib_b = asin(sin(W)*sin(I));
    
    return p;
}

// Newton Solver
double solve_event(double jd, double target_deg) {
    double t = jd;
    for(int k=0; k<6; k++) {
        MoonPhys p = compute(t);
        double err = n360(p.age - target_deg + 180.0) - 180.0;
        if(fabs(err) < 1e-5) break;
        t -= err / 12.19075;
    }
    return t;
}

void fmt_dt(double dt, char *b) {
    double s = fabs(dt)*86400;
    int d=s/86400, h=(s-d*86400)/3600, m=(s-d*86400-h*3600)/60, sc=s-d*86400-h*3600-m*60;
    sprintf(b, "%2d %02d:%02d:%02d", d, h, m, sc);
}

// Shader: The heart of the visual
char render(int r, int c, MoonPhys p) {
    // Sphere Coordinates
    double y = (r-3.5)/4.2;
    double x = (c-7.0)/7.8; // Aspect ratio correction for text
    double z2 = 1.0 - x*x - y*y;
    if(z2 < 0) return ' ';
    double z = sqrt(z2);

    // Texture Rotation (Libration)
    double x_t = x*cos(p.lib_l) - z*sin(p.lib_l);
    // Simple 2D shift approx for the static ASCII map

    
    // Lighting
    // Sun Vector determination based on Age
    double age_rad = p.age * D2R;
    double sun_x = -sin(age_rad);
    
    double s_ang = -age_rad;
    double sx = -sin(s_ang); // Flip sign to match visual expectation
    double sz = -cos(s_ang);
    
    double dot = x*sx + z*sz;
    
    // Shadow Application
    if (dot < -0.1) return ' '; // Shadow
    
    return SURFACE_MAP[r][c];
}

const char* name_phase(double d) {
    if(d>355||d<5) return "New";
    if(d<85) return "Waxing Crescent"; if(d<95) return "First Quarter";
    if(d<175) return "Waxing Gibbous"; if(d<185) return "Full";
    if(d<265) return "Waning Gibbous"; if(d<275) return "Last Quarter";
    return "Waning Crescent";
}

int main() {
    time_t now = time(NULL);
    double jd = get_jd(now);
    
    MoonPhys p = compute(jd);
    int pct = (int)(p.illum * 100.0 + 0.5);

    // Event Solver
    double q1 = solve_event(jd, 90.0);
    double q2 = solve_event(jd-20, 90.0);
    double q = (fabs(q1-jd) < fabs(jd-q2)) ? q1 : q2;
    char qs = (q > jd) ? '+' : '-';

    double f1 = solve_event(jd, 180.0);
    double f2 = solve_event(jd-20, 180.0);
    double f = (fabs(f1-jd) < fabs(jd-f2)) ? f1 : f2;
    char fs = (f > jd) ? '+' : '-';

    char bq[32], bf[32];
    fmt_dt(q-jd, bq); fmt_dt(f-jd, bf);

    // Draw
    char grid[8][16];
    for(int r=0; r<8; r++) {
        for(int c=0; c<15; c++) grid[r][c] = render(r, c, p);
        grid[r][15] = 0;
    }

    char st[64];
    sprintf(st, "The Moon is %s(%d%%)", name_phase(p.age), pct);

    // Strict format output
    printf("%-15s  %s\n", grid[0], st);
    printf("%-15s    \n",     grid[1]);
    printf("%-15s  First Quater %c\n", grid[2], qs);
    printf("%-15s  %s\n",    grid[3], bq);
    printf("%-15s  Full Moon %c\n",    grid[4], fs);
    printf("%-15s  %s\n",    grid[5], bf);
    printf("%-15s    \n",     grid[6]);
    printf("%-15s    \n",     grid[7]);

    return 0;
}