sublunar.almanac/lib/sweph/pmom.c

/**
* @preserve Copyright (c) 2018 T. F. Raaion, www.sublunar.space
* License MIT: http://www.opensource.org/licenses/MIT
*/

#define ONLINE 1
#define OFFLINE 2

#ifndef USECASE
#error "USECASE needs to be set"
#endif

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include "swephexp.h"

#if USECASE == OFFLINE
#include <emscripten.h>
#endif

/**
 * converts Julian Day to UNIX epoch
 */
long jul_to_epoch(double jd)
{
  struct tm   ts;
  int32       year, month, day, hour, minute, second;
  double      hours;
  time_t      epoch;

  // JD to year, month, day, hours fraction
  swe_revjul(jd, SE_GREG_CAL, &year, &month, &day, &hours);

  // convert hours fraction to full hour, minute, second
  hour = floor(hours);
  minute = floor((hours - hour) * 60.0);
  second = floor((((hours - hour) * 60.0) - minute) * 60.0);

  // fill tm struct and covert to time_t
  ts.tm_year = year - 1900;
  ts.tm_mon = month - 1;
  ts.tm_mday = day;
  ts.tm_hour = hour;
  ts.tm_min = minute;
  ts.tm_sec = second;
  ts.tm_isdst = -1;
  epoch = timegm(&ts);

  return (long) epoch;
}

/**
 * converts year, month, day, hour, minute, second to UNIX epoch
 */
long time_to_epoch(int32 year, int32 month, int32 day, int32 hour, int32 minute, int32 second)
{
  struct tm   ts;
  time_t      epoch;

  // fill tm struct and convert
  ts.tm_year = year - 1900;
  ts.tm_mon = month - 1;
  ts.tm_mday = day;
  ts.tm_hour = hour;
  ts.tm_min = minute;
  ts.tm_sec = second;
  ts.tm_isdst = -1;
  epoch = timegm(&ts);

  return (long) epoch;
}

/**
 * converts UNIX epoch to time array
 */
int * epoch_to_time(long ts, int *time)
{
    time_t        timestamp;
    struct tm     *tmp;

    // convert timestamp to tm struct
    timestamp = ts;
    tmp = gmtime(&timestamp);

    // read the data
    time[0] = tmp->tm_year+1900;
    time[1] = tmp->tm_mon+1;
    time[2] = tmp->tm_mday;
    time[3] = tmp->tm_hour;
    time[4] = tmp->tm_min;
    time[5] = tmp->tm_sec;

    return 0;
}

/**
 * converts UNIX epoch to JD
 */
double epoch_to_jul(long ts)
{
    int32         time[6];
    double        hours;

    epoch_to_time(ts, time);
    hours = (double) time[3] + (double) time[4] / 60 + (double) time[5] / 3600;

    // return JD
    return swe_julday(time[0],time[1],time[2],hours,SE_GREG_CAL);
}

/**
 * calculates planetary ephemeris data for a moment and loads it into arrays
 */
int * planetary_moment(double lat, double lon, long unixtime, double *lunar, double *planets) {

    // double lunar[0] = lunar day
    // double lunar[1] = lunar phase angle
    // double lunar[2] = lunar phase
    // double planet[0] = ecliptic longitude sun
    // double planet[1] = ecliptic longitude moon
    // double planet[2] = ecliptic longitude mercury
    // double planet[3] = ecliptic longitude venus
    // double planet[4] = ecliptic longitude mars
    // double planet[5] = ecliptic longitude jupiter
    // double planet[6] = ecliptic longitude saturn

    double        latitude = 0.0;
    double        longitude = 0.0;
    long          epoch = 0;
    struct tm     *tmp;
    time_t        timestamp = 0;
    int32         year, month, day, hour, minute, second, iflag, iflgret;
    int           pl, ipl, rsmi, return_code, dl, lhour, nm;
    double        hours, tjd_ts, te, lastnew;
    double        x2[6], datm[2], geopos[3], moon[20], llp[20], lunarday[31], data[13], datb[10];
    char          serr[AS_MAXCH], starname[255];

    datm[0] = 1013.25; // atmospheric pressure
    datm[1] = 15; // atmospheric temperature;

    // exact new moon occurances taken from Steve Mosiher's website (moshier.org)
    // new moon moments from 1.1.1970 until 31.12.2100
    // replaced a formula for speed reasons
    double newmoons[2474] = {2415021.07774, 2415050.55740, 2415079.97593, 2415109.35449, 2415138.72449, 2415168.11786, 2415197.56062, 2415227.07140, 2415256.66155, 2415286.33129, 2415316.06058, 2415345.80353, 2415375.50085, 2415405.10821, 2415434.61475, 2415464.03679, 2415493.40102, 2415522.73438, 2415552.06442, 2415581.42396, 2415610.85238, 2415640.38793, 2415670.04948, 2415699.81533, 2415729.62018, 2415759.38510, 2415789.05661, 2415818.61816, 2415848.07638, 2415877.44794, 2415906.75758, 2415936.04107, 2415965.34522, 2415994.72177, 2416024.21468, 2416053.84272, 2416083.58638, 2416113.39218, 2416143.19346, 2416172.93024, 2416202.55976, 2416232.06343, 2416261.45109, 2416290.75750, 2416320.03204, 2416349.32707, 2416378.68798, 2416408.14604, 2416437.71530, 2416467.39295, 2416497.15739, 2416526.96159, 2416556.73555, 2416586.41167, 2416615.95688, 2416645.38211, 2416674.72728, 2416704.04034, 2416733.36298, 2416762.72560, 2416792.15042, 2416821.65722, 2416851.26199, 2416880.96246, 2416910.72175, 2416940.47452, 2416970.15949, 2416999.74751, 2417029.24295, 2417058.66842, 2417088.05091, 2417117.41625, 2417146.79006, 2417176.19928, 2417205.66922, 2417235.21484, 2417264.83141, 2417294.49436, 2417324.17114, 2417353.83384, 2417383.46221, 2417413.04100, 2417442.56066, 2417472.02320, 2417501.44624, 2417530.85857, 2417560.28773, 2417589.74789, 2417619.23816, 2417648.75333, 2417678.29571, 2417707.87445, 2417737.49296, 2417767.13683, 2417796.77527, 2417826.37758, 2417855.93089, 2417885.44356, 2417914.93218, 2417944.40512, 2417973.85863, 2418003.28949, 2418032.70982, 2418062.14798, 2418091.63512, 2418121.18852, 2418150.80329, 2418180.45739, 2418210.12440, 2418239.78220, 2418269.41172, 2418298.99289, 2418328.50815, 2418357.95279, 2418387.34109, 2418416.70227, 2418446.07085, 2418475.47801, 2418504.94768, 2418534.49623, 2418564.13090, 2418593.84254, 2418623.59595, 2418653.33241, 2418682.99374, 2418712.55062, 2418742.00847, 2418771.39228, 2418800.73110, 2418830.05311, 2418859.38895, 2418888.77559, 2418918.25390, 2418947.85546, 2418977.58056, 2419007.38242, 2419037.18135, 2419066.90588, 2419096.52147, 2419126.02619, 2419155.43402, 2419184.76677, 2419214.05534, 2419243.34176, 2419272.67651, 2419302.10926, 2419331.67317, 2419361.36764, 2419391.15295, 2419420.96509, 2419450.73883, 2419480.42262, 2419509.98619, 2419539.42612, 2419568.76638, 2419598.05082, 2419627.33169, 2419656.65866, 2419686.06989, 2419715.58669, 2419745.21301, 2419774.93634, 2419804.72339, 2419834.51552, 2419864.24170, 2419893.85018, 2419923.33116, 2419952.71255, 2419982.04024, 2420011.35972, 2420040.70612, 2420070.10344, 2420099.57032, 2420129.12411, 2420158.77358, 2420188.50140, 2420218.25629, 2420247.97336, 2420277.60746, 2420307.14796, 2420336.60993, 2420366.01833, 2420395.39825, 2420424.77336, 2420454.16809, 2420483.60766, 2420513.11242, 2420542.68825, 2420572.32107, 2420601.98321, 2420631.64656, 2420661.28971, 2420690.89629, 2420720.45296, 2420749.95323, 2420779.40428, 2420808.82806, 2420838.25257, 2420867.69829, 2420897.17053, 2420926.66506, 2420956.18130, 2420985.72845, 2421015.31756, 2421044.94675, 2421074.59394, 2421104.22541, 2421133.81538, 2421163.35910, 2421192.86835, 2421222.35497, 2421251.81942, 2421281.25628, 2421310.67021, 2421340.08430, 2421369.53251, 2421399.04315, 2421428.62497, 2421458.26455, 2421487.93575, 2421517.61184, 2421547.26979, 2421576.88707, 2421606.44147, 2421635.91991, 2421665.32807, 2421694.69054, 2421724.04227, 2421753.41859, 2421782.84870, 2421812.35387, 2421841.94710, 2421871.62865, 2421901.37607, 2421931.13835, 2421960.85012, 2421990.46317, 2422019.96634, 2422049.37837, 2422078.72945, 2422108.04995, 2422137.36993, 2422166.72355, 2422196.15080, 2422225.69027, 2422255.36071, 2422285.13866, 2422314.95506, 2422344.72704, 2422374.39917, 2422403.95545, 2422433.40494, 2422462.76745, 2422492.07034, 2422521.35070, 2422550.65550, 2422580.03599, 2422609.53519, 2422639.17025, 2422668.91936, 2422698.72695, 2422728.52564, 2422758.25640, 2422787.87858, 2422817.37598, 2422846.76019, 2422876.06683, 2422905.34546, 2422934.64796, 2422964.01833, 2422993.48527, 2

    latitude = lat;
    longitude = lon;
    epoch = unixtime;
    timestamp = unixtime;

    // take timestamp, convert it to year, month, day, decimal hour
    tmp = gmtime(&timestamp);
    year = tmp->tm_year+1900;
    month = tmp->tm_mon+1;
    day = tmp->tm_mday;
    hour = tmp->tm_hour;
    minute = tmp->tm_min;
    second = tmp->tm_sec;
    hours = (double) hour + (double) minute / 60 + (double) second / 3600;

    // set ephemeris path and flags
    swe_set_ephe_path("lib");
    iflag = SEFLG_SPEED;
    ipl = SE_MOON;
    rsmi = SE_CALC_RISE;

    // set julian day and calculate ephemeris time
    tjd_ts = epoch_to_jul(epoch);
    te = tjd_ts + swe_deltat(tjd_ts);

    // calculate geocentric ecliptic longitudes for planets
    for (pl = SE_SUN; pl <= SE_SATURN; pl++) {
      iflgret = swe_calc(te, pl, iflag, x2, serr);
      planets[pl] = x2[0];
      planets[pl+7] = x2[3];
    }

    // calculate geocentric ecliptic longitude for mean node
    iflgret = swe_calc(te, SE_MEAN_NODE, iflag, x2, serr);
    planets[17] = x2[0];

    // calculate house cusps, ascendant and mc
    swe_houses_ex(tjd_ts, 0, latitude, longitude, 'P', data, datb);
    planets[14] = datb[0]; // asc
    planets[15] = datb[1]; // mc
    planets[16] = datb[2]; // armc

    // calculate moon phase
    return_code = swe_pheno(te, SE_MOON, iflag, &moon, serr);
    if (return_code == ERR) printf("%s\n", serr);
    lunar[1] = moon[0];
    lunar[2] = moon[1];

    // calculate last New Moon from array data
    for (nm = 0; nm <= 2474; nm++) {
      if (newmoons[nm] > tjd_ts ) {
        lastnew = newmoons[nm-1];
        break;
      }
    }

    // calculate lunar day counting moonrises from last new moon
    while ( tjd_ts - lastnew > 29.530588853 ) lastnew += 29.530588853;
    lunarday[0] = lastnew;
    for (dl = 0; lunarday[0]+dl < lunarday[0]+29.530588853 && tjd_ts > lunarday[0]+dl; dl++) {
      swe_rise_trans(lunarday[0]+dl, ipl, starname, iflag, rsmi, geopos, datm[0], datm[1], &lunarday[dl+1], serr);
      lhour = dl+1;
    }
    lunar[0] = (double) lhour;

    return 0;
}

/**
 * calculates sunrise and sunset times to determine planetary day and hour
 */
int * solar_moment(double lat, double lon, long unixtime, long *solar) {

    long          epoch = 0;
    struct tm     *tmp;
    time_t        timestamp = 0;
    int32         year, month, day, hour, minute, second, iflag, pday = 0;
    int           ipl, rsmi, hourlength[3];
    double        hours, tjd, td;
    double        datm[2], geopos[3], lastset, lastrise, nextset, nextrise, onextrise, onextset;
    char          serr[AS_MAXCH], starname[255], pds[80];

    datm[0] = 1013.25; // atmospheric pressure
    datm[1] = 15; // atmospheric temperature;

    // take timestamp, convert it to year, month, day, decimal hour
    epoch = unixtime;
    timestamp = unixtime;
    tmp = gmtime(&timestamp);
    year = tmp->tm_year+1900;
    month = tmp->tm_mon+1;
    day = tmp->tm_mday;
    hour = tmp->tm_hour;
    minute = tmp->tm_min;
    second = tmp->tm_sec;
    hours = (double) hour + (double) minute / 60 + (double) second / 3600;
    pday = tmp->tm_wday;

    // prepare variables for rise and set calculations
    geopos[0] = lon;
    geopos[1] = lat;
    geopos[2] = 0;
    swe_set_topo(geopos[0], geopos[1], geopos[2]);
    tjd = swe_julday(year,month,day,0.0,SE_GREG_CAL);

    if ( lon > 70.0 ) tjd = tjd - 1;

    // calculate sunrise and sunset of the day and the days before and after tomorrow
    ipl = SE_SUN;

    // next rise and set
    swe_rise_trans(tjd, ipl, starname, iflag, SE_CALC_RISE, geopos, datm[0], datm[1], &nextrise, serr);
    swe_rise_trans(tjd, ipl, starname, iflag, SE_CALC_SET, geopos, datm[0], datm[1], &nextset, serr);

    // previous rise and set
    swe_rise_trans(tjd-1, ipl, starname, iflag, SE_CALC_RISE, geopos, datm[0], datm[1], &lastrise, serr);
    swe_rise_trans(tjd-1, ipl, starname, iflag, SE_CALC_SET, geopos, datm[0], datm[1], &lastset, serr);

    // next rise after next rise
    swe_rise_trans(tjd+1, ipl, starname, iflag, SE_CALC_RISE, geopos, datm[0], datm[1], &onextrise, serr);
    swe_rise_trans(tjd+1, ipl, starname, iflag, SE_CALC_SET, geopos, datm[0], datm[1], &onextset, serr);


    if ( nextset < nextrise ) {
        lastset = nextset;
        nextset = onextset;
    }

    hourlength[0] = (jul_to_epoch(nextrise) - jul_to_epoch(lastset)) / 12; // length of night hours yesterday
    hourlength[1] = (jul_to_epoch(nextset) - jul_to_epoch(nextrise)) / 12; // length of day hours today
    hourlength[2] = (jul_to_epoch(onextrise) - jul_to_epoch(nextset)) / 12; // length of night hours today

    // calculate planetary hour
    solar[0] = jul_to_epoch(lastset); // set yesterday
    solar[1] = jul_to_epoch(nextrise); // rise today
    solar[2] = jul_to_epoch(nextset); // set today
    solar[3] = jul_to_epoch(onextrise); // rise tomorrow
    solar[4] = hourlength[0]; // night hour length yesterday
    solar[5] = hourlength[1]; // day hour length today
    solar[6] = hourlength[2]; // night hour length today
    solar[7] = pday; // planetary day

    return 0;
}

/**
 * calculates JSON Ephemeris Object
 */
const char * pmom(int days, long epoch, double latitude, double longitude, int buflen) // days, timestamp, latitude, longitude
{
    long          queryday;
    struct tm     *tmp;
    time_t        timestamp = 0;
    int32         year, month, day, hour, minute, second, j, h, m, interval, start, hm, time[6];
    double        hours;
    long          solar[9] = {0};
    double        lunar[3] = {0.0}, planet[18] = {0.0};
    char*         Buffer = malloc(buflen);
    int           length = 0;

    timestamp = epoch;

    // take timestamp, convert it to year, month, day, decimal hour
    tmp = gmtime(&timestamp);
    year = tmp->tm_year+1900;
    month = tmp->tm_mon+1;
    day = tmp->tm_mday;
    hour = tmp->tm_hour;
    minute = tmp->tm_min;
    second = tmp->tm_sec;
    hours = (double) hour + (double) minute / 60 + (double) second / 3600;

    length += snprintf(Buffer+length, buflen-length, "{ ");

    // calculate timestamp moment
    solar_moment(latitude, longitude, epoch, solar);
    planetary_moment(latitude, longitude, epoch, lunar, planet);
    if ( epoch < solar[2]) {
      start = solar[1];
      interval = solar[5];
      m = 0;
    }
    if ( solar[2] < epoch ) {
      start = solar[2];
      interval = solar[6];
      m = 12;
    }
    for (hm = 1; hm <= 12; hm++) {
      if (epoch < (start + hm * interval) ) {
        hm = hm + m;
        break;
      }
    }
    length += snprintf(Buffer+length, buflen-length, "\"moment\":[{ \"ts\": %d, \"gst\": %f, \"lunar\": { \"day\": %d, \"angle\": %f, \"phase\": %f }, \"planetary\": { \"day\": { \"no\": %d, \"start\": %d, \"end\": %d }, \"night\": { \"start\": %d, \"end\": %d }, \"hour\": { \"no\": %d, \"start\": %d, \"end\": %d, \"length\": { \"day\": %d, \"night\": %d } } }, \"ephemeris\": { \"sun\": { \"deg\": %f, \"speed\": %f }, \"moon\": { \"deg\": %f, \"speed\": %f }, \"mercury\": { \"deg\": %f, \"speed\": %f }, \"venus\": { \"deg\": %f, \"speed\": %f }, \"mars\": { \"deg\": %f, \"speed\": %f }, \"jupiter\": { \"deg\": %f, \"speed\": %f }, \"saturn\": { \"deg\": %f, \"speed\": %f }, \"node\": { \"deg\": %f }, \"asc\": { \"deg\": %f }, \"mc\": { \"deg\": %f } } }], ", epoch, (planet[16] / 15.0), (int) lunar[0], lunar[1], lunar[2], solar[7], solar[1], solar[2]-1, solar[2], solar[3]-1, hm, (start + hm * interval), (start + hm * interval)+interval-1, solar[5], solar[6], planet[0], planet[7], planet[1], planet[8], planet[2], planet[9], planet[3], planet[10], planet[4], planet[11], planet[5], planet[12], planet[6], planet[13], planet[16], planet[14], planet[15] );

    if ( days > 0 ) {

      // start of Almanac
      queryday = time_to_epoch(year, month, day, 0, 0, 0);
      j = 0;
      length += snprintf(Buffer+length, buflen-length, "\"query\":[");

      // calculate all hours for all days
      while ( j <= days ) {
        solar_moment(latitude, longitude, queryday, solar);
        h = interval = start = m = 0;
        while ( h < 24 ) {
          if ( h < 12 ) {
            interval = solar[5];
            start = solar[1];
            m = h;
          } else {
            interval = solar[6];
            start = solar[2];
            m = h - 12;
          }
          hm = start + m * interval;

          planetary_moment(latitude, longitude, hm, lunar, planet);

          length += snprintf(Buffer+length, buflen-length, "{ \"ts\": %d, \"gst\": %f, \"lunar\": { \"day\": %d, \"angle\": %f, \"phase\": %f }, \"planetary\": { \"day\": { \"no\": %d, \"start\": %d, \"end\": %d }, \"night\": { \"start\": %d, \"end\": %d }, \"hour\": { \"no\": %d, \"start\": %d, \"end\": %d, \"length\": { \"day\": %d, \"night\": %d } } }, \"ephemeris\": { \"sun\": { \"deg\": %f, \"speed\": %f }, \"moon\": { \"deg\": %f, \"speed\": %f }, \"mercury\": { \"deg\": %f, \"speed\": %f }, \"venus\": { \"deg\": %f, \"speed\": %f }, \"mars\": { \"deg\": %f, \"speed\": %f }, \"jupiter\": { \"deg\": %f, \"speed\": %f }, \"saturn\": { \"deg\": %f, \"speed\": %f }, \"node\": { \"deg\": %f }, \"asc\": { \"deg\": %f }, \"mc\": { \"deg\": %f } } }", hm, (planet[16] / 15.0), (int) lunar[0], lunar[1], lunar[2], solar[7], solar[1], solar[2]-1, solar[2], solar[3]-1, h, hm, hm+interval-1, solar[5], solar[6], planet[0], planet[7], planet[1], planet[8], planet[2], planet[9], planet[3], planet[10], planet[4], planet[11], planet[5], planet[12], planet[6], planet[13], planet[16], planet[14], planet[15] );

          //length += snprintf(Buffer+length, buflen-length, "%d setyes %d risetod %d settod %d risetom %d nhly %d dhlt %d nhlt %d pday %d\n", hm, solar[0], solar[1], solar[2], solar[3], solar[4], solar[5], solar[6], solar[7]);

          if(j == days && h == 23) length += snprintf(Buffer+length, buflen-length, "");
          else length += snprintf(Buffer+length, buflen-length, ", ");

          h++;
        }
        queryday += 86400;
        j++;
  #if USECASE == OFFLINE
        EM_ASM_({
            NProgress.set($0/$1);
        }, j, days*2);
  #endif
      };
    }
    length += snprintf(Buffer+length, buflen-length, "]}");
    return Buffer;
}

#if USECASE == ONLINE
int main(int argc,char* argv[]) // days, timestamp, latitude, longitude
{
    double        latitude = 0.0;
    double        longitude = 0.0;
    long          epoch = 0;
    time_t        timestamp = 0;
    int32         days = 1, buflen;

    // read and assign commandline arguments
    if ( argc > 1 ) {
      days = atof(argv[1]);
    }
    if ( argc > 3 ) {
      latitude = atof(argv[3]);
      longitude = atof(argv[4]);
    }
    if ( argc > 2 ) {
      sscanf(argv[2], "%lu", &epoch);
      timestamp = epoch;
    }
    else {
      time(&timestamp);
      epoch = (int) timestamp;
    }

    timestamp = epoch;
    buflen = (1 + days) * 100000;
    printf(pmom( days, epoch, latitude, longitude, buflen));

}
#endif

#if USECASE == OFFLINE
EMSCRIPTEN_KEEPALIVE
const char * get(int days, long epoch, double latitude, double longitude) // days, timestamp, latitude, longitude
{
    time_t        timestamp = 0;
    int32         buflen;

    timestamp = epoch;
    buflen = (1 + days) * 100000;
    return pmom( days, epoch, latitude, longitude, buflen);
}
#endif