#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
//#define TRUE 1
//#define FALSE 0
/* Date for which to calculate phase *//* Illuminated fraction *//* Age of moon in days *//* Distance in kilometres *//* Angular diameter in degrees *//* Distance to Sun *//* Sun's angular diameter *//* Astronomical constants */
#define epoch 2444238.5 /* 1980 January 0.0 */
/* Constants defining the Sun's apparent orbit */
#define elonge 278.833540 /* Ecliptic longitude of the Sun
at epoch 1980.0 */
#define elongp 282.596403 /* Ecliptic longitude of the Sun at
perigee */
#define eccent 0.016718 /* Eccentricity of Earth's orbit */
#define sunsmax 1.495985e8 /* Semi-major axis of Earth's orbit, km */
#define sunangsiz 0.533128 /* Sun's angular size, degrees, at
semi-major axis distance */
/* Elements of the Moon's orbit, epoch 1980.0 */
#define mmlong 64.975464 /* Moon's mean longitude at the epoch */
#define mmlongp 349.383063 /* Mean longitude of the perigee at the
epoch */
#define mlnode 151.950429 /* Mean longitude of the node at the
epoch */
#define minc 5.145396 /* Inclination of the Moon's orbit */
#define mecc 0.054900 /* Eccentricity of the Moon's orbit */
#define mangsiz 0.5181 /* Moon's angular size at distance a
from Earth */
#define msmax 384401.0 /* Semi-major axis of Moon's orbit in km */
#define mparallax 0.9507 /* Parallax at distance a from Earth */
#define synmonth 29.53058868 /* Synodic month (new Moon to new Moon) */
#define lunatbase 2423436.0 /* Base date for E. W. Brown's numbered
series of lunations (1923 January 16) */
/* Properties of the Earth */
#define earthrad 6378.16 /* Radius of Earth in kilometres */
#define PI 3.14159265358979323846 /* Assume not near black hole nor in
Tennessee */
/* Handy mathematical functions */
#define sgn(x) (((x) < 0) ? -1 : ((x) > 0 ? 1 : 0)) /* Extract sign */
//#define abs(x) ((x) < 0 ? (-(x)) : (x)) /* Absolute val */
#define fixangle(a) ((a) - 360.0 * (floor((a) / 360.0))) /* Fix angle */
#define torad(d) ((d) * (PI / 180.0)) /* Deg->Rad */
#define todeg(d) ((d) * (180.0 / PI)) /* Rad->Deg */
#define dsin(x) (sin(torad((x)))) /* Sin from deg */
#define dcos(x) (cos(torad((x)))) /* Cos from deg */
// int mm, dd;
"time=%d\n""jtime=%f\n",jd);
// mm = gm->tm_mon + 1;
// dd = gm->tm_mday;
// drawmoon(hDC, "colourmoon", 64, 32, 27, 28,
// tm.tmMaxCharWidth * mgeom.moonbitx, 10, p, mm, dd);
/* Update textual information for open window */
/* Local time */"Local time %02d-%02d-%4d %02d:%02d:%02d\n""Moon age %d day, %02d:%02d\n""Moon phase %.0f%%\n""Moon dist. %.0f km %.2f e.rd.\n""Sun dist. %.0f km %.2f au\n"/* JDATE -- Convert internal GMT date and time to Julian day and fraction. *//* Compute century *//* JTIME -- Convert internal GMT date and time to astronomical
Julian time (i.e. Julian date plus day fraction,
expressed as a double). */"jdate=%d\n"/* JYEAR -- Convert Julian date to year, month, day, which are
returned via integer pointers to integers. *//* Astronomical to civil *//* JHMS -- Convert Julian time to hour, minutes, and seconds. *//* Astronomical to civil *//* KEPLER -- Solve the equation of Kepler. */#define EPSILON 1E-6
/* PHASE -- Calculate phase of moon as a fraction:
The argument is the time for which the phase is requested,
expressed as a Julian date and fraction. Returns the terminator
phase angle as a percentage of a full circle (i.e., 0 to 1), and
stores into pointer arguments the illuminated fraction of the
Moon's disc, the Moon's age in days and fraction, the distance of
the Moon from the centre of the Earth, and the angular diameter
subtended by the Moon as seen by an observer at the centre of the
Earth.
*//* Date for which to calculate phase *//* Illuminated fraction *//* Age of moon in days *//* Distance in kilometres *//* Angular diameter in degrees *//* Distance to Sun *//* Sun's angular diameter *//* MN,*/ Ev, Ae, A3, MmP,
mEc, A4, lP, V, lPP, // NP, y, x, Lambdamoon, BetaM,
MoonAge, MoonPhase,
MoonDist, MoonDFrac, MoonAng, // MoonPar,
F, SunDist, SunAng;
/* Calculation of the Sun's position */
Day = pdate - epoch; /* Date within epoch */
N = fixangle((360 / 365.2422) * Day); /* Mean anomaly of the Sun *//* Convert from perigee
co-ordinates to epoch 1980.0 */
Ec = kepler(M, eccent); /* Solve equation of Kepler */
Ec = sqrt((1 + eccent) / (1 - eccent)) * tan(Ec / 2);
Ec = 2 * todeg(atan(Ec)); /* True anomaly *//* Sun's geocentric ecliptic
longitude */
/* Orbital distance factor */
F = ((1 + eccent * cos(torad(Ec))) / (1 - eccent * eccent));
SunDist = sunsmax / F; /* Distance to Sun in km */
SunAng = F * sunangsiz; /* Sun's angular size in degrees */
/* Calculation of the Moon's position */
/* Moon's mean longitude *//* Moon's mean anomaly *//* Moon's ascending node mean longitude */
// MN = fixangle(mlnode - 0.0529539 * Day);
/* Evection */
Ev = 1.2739 * sin(torad(2 * (ml - Lambdasun) - MM));
/* Annual equation */
Ae = 0.1858 * sin(torad(M));
/* Correction term */
A3 = 0.37 * sin(torad(M));
/* Corrected anomaly */
MmP = MM + Ev - Ae - A3;
/* Correction for the equation of the centre */
mEc = 6.2886 * sin(torad(MmP));
/* Another correction term */
A4 = 0.214 * sin(torad(2 * MmP));
/* Corrected longitude */
lP = ml + Ev + mEc - Ae + A4;
/* Variation */
V = 0.6583 * sin(torad(2 * (lP - Lambdasun)));
/* True longitude */
lPP = lP + V;
/* Corrected longitude of the node */
// NP = MN - 0.16 * sin(torad(M));
/* Y inclination coordinate */
// y = sin(torad(lPP - NP)) * cos(torad(minc));
/* X inclination coordinate */
// x = cos(torad(lPP - NP));
/* Ecliptic longitude */
// Lambdamoon = todeg(atan2(y, x));
// Lambdamoon += NP;
/* Ecliptic latitude */
// BetaM = todeg(asin(sin(torad(lPP - NP)) * sin(torad(minc))));
/* Calculation of the phase of the Moon */
/* Age of the Moon in degrees */
MoonAge = lPP - Lambdasun;
/* Phase of the Moon */
MoonPhase = (1 - cos(torad(MoonAge))) / 2;
/* Calculate distance of moon from the centre of the Earth */
MoonDist = (msmax * (1 - mecc * mecc)) /
(1 + mecc * cos(torad(MmP + mEc)));
/* Calculate Moon's angular diameter */
MoonDFrac = MoonDist / msmax;
MoonAng = mangsiz / MoonDFrac;
/* Calculate Moon's parallax */
// MoonPar = mparallax / MoonDFrac;
"%lf\n"