astrognome/swe/src/swedate.c

/*********************************************************
  $Header: /home/dieter/sweph/RCS/swedate.c,v 1.75 2009/04/08 07:17:29 dieter Exp $
  version 15-feb-89 16:30
  
  swe_date_conversion()
  swe_revjul()
  swe_julday()

************************************************************/
/* Copyright (C) 1997 - 2008 Astrodienst AG, Switzerland.  All rights reserved.
  
  License conditions
  ------------------

  This file is part of Swiss Ephemeris.
  
  Swiss Ephemeris is distributed with NO WARRANTY OF ANY KIND.  No author
  or distributor accepts any responsibility for the consequences of using it,
  or for whether it serves any particular purpose or works at all, unless he
  or she says so in writing.  

  Swiss Ephemeris is made available by its authors under a dual licensing
  system. The software developer, who uses any part of Swiss Ephemeris
  in his or her software, must choose between one of the two license models,
  which are
  a) GNU public license version 2 or later
  b) Swiss Ephemeris Professional License
  
  The choice must be made before the software developer distributes software
  containing parts of Swiss Ephemeris to others, and before any public
  service using the developed software is activated.

  If the developer choses the GNU GPL software license, he or she must fulfill
  the conditions of that license, which includes the obligation to place his
  or her whole software project under the GNU GPL or a compatible license.
  See http://www.gnu.org/licenses/old-licenses/gpl-2.0.html

  If the developer choses the Swiss Ephemeris Professional license,
  he must follow the instructions as found in http://www.astro.com/swisseph/ 
  and purchase the Swiss Ephemeris Professional Edition from Astrodienst
  and sign the corresponding license contract.

  The License grants you the right to use, copy, modify and redistribute
  Swiss Ephemeris, but only under certain conditions described in the License.
  Among other things, the License requires that the copyright notices and
  this notice be preserved on all copies.

  Authors of the Swiss Ephemeris: Dieter Koch and Alois Treindl

  The authors of Swiss Ephemeris have no control or influence over any of
  the derived works, i.e. over software or services created by other
  programmers which use Swiss Ephemeris functions.

  The names of the authors or of the copyright holder (Astrodienst) must not
  be used for promoting any software, product or service which uses or contains
  the Swiss Ephemeris. This copyright notice is the ONLY place where the
  names of the authors can legally appear, except in cases where they have
  given special permission in writing.

  The trademarks 'Swiss Ephemeris' and 'Swiss Ephemeris inside' may be used
  for promoting such software, products or services.
*/

/*
  swe_date_conversion():
  This function converts some date+time input {d,m,y,uttime}
  into the Julian day number tjd.
  The function checks that the input is a legal combination
  of dates; for illegal dates like 32 January 1993 it returns ERR
  but still converts the date correctly, i.e. like 1 Feb 1993.
  The function is usually used to convert user input of birth data
  into the Julian day number. Illegal dates should be notified to the user.

  Be aware that we always use astronomical year numbering for the years
  before Christ, not the historical year numbering.
  Astronomical years are done with negative numbers, historical
  years with indicators BC or BCE (before common era).
  Year 0 (astronomical)  	= 1 BC historical.
  year -1 (astronomical) 	= 2 BC
  etc.
  Many users of Astro programs do not know about this difference.

  Return: OK or ERR (for illegal date)
*********************************************************/

# include "swephexp.h"
# include "sweph.h"

static AS_BOOL init_leapseconds_done = FALSE;


int FAR PASCAL_CONV swe_date_conversion(int y,
		     int m,
		     int d,		/* day, month, year */
		     double uttime, 	/* UT in hours (decimal) */
		     char c,		/* calendar g[regorian]|j[ulian] */
		     double *tjd)
{
  int rday, rmon, ryear;
  double rut, jd;
  int gregflag = SE_JUL_CAL;
  if (c == 'g')
    gregflag = SE_GREG_CAL;
  rut = uttime;		/* hours UT */
  jd = swe_julday(y, m, d, rut, gregflag);
  swe_revjul(jd, gregflag, &ryear, &rmon, &rday, &rut);
  *tjd = jd;
  if (rmon == m && rday == d && ryear == y) {
    return OK;
  } else {
    return ERR;
  }
}	/* end date_conversion */

/*************** swe_julday ********************************************
 * This function returns the absolute Julian day number (JD)
 * for a given calendar date.
 * The arguments are a calendar date: day, month, year as integers,
 * hour as double with decimal fraction.
 * If gregflag = SE_GREG_CAL (1), Gregorian calendar is assumed,
 * if gregflag = SE_JUL_CAL (0),Julian calendar is assumed.
  
 The Julian day number is a system of numbering all days continously
 within the time range of known human history. It should be familiar
 to every astrological or astronomical programmer. The time variable
 in astronomical theories is usually expressed in Julian days or
 Julian centuries (36525 days per century) relative to some start day;
 the start day is called 'the epoch'.
 The Julian day number is a double representing the number of
 days since JD = 0.0 on 1 Jan -4712, 12:00 noon (in the Julian calendar).
 
 Midnight has always a JD with fraction .5, because traditionally
 the astronomical day started at noon. This was practical because
 then there was no change of date during a night at the telescope.
 From this comes also the fact the noon ephemerides were printed
 before midnight ephemerides were introduced early in the 20th century.
 
 NOTE: The Julian day number must not be confused with the Julian 
 calendar system.

 Be aware the we always use astronomical year numbering for the years
 before Christ, not the historical year numbering.
 Astronomical years are done with negative numbers, historical
 years with indicators BC or BCE (before common era).
 Year 0 (astronomical)  	= 1 BC
 year -1 (astronomical) 	= 2 BC
 etc.

 Original author: Marc Pottenger, Los Angeles.
 with bug fix for year < -4711   15-aug-88 by Alois Treindl
 (The parameter sequence m,d,y still indicates the US origin,
  be careful because the similar function date_conversion() uses
  other parameter sequence and also Astrodienst relative juldate.)
 
 References: Oliver Montenbruck, Grundlagen der Ephemeridenrechnung,
             Verlag Sterne und Weltraum (1987), p.49 ff
 
 related functions: swe_revjul() reverse Julian day number: compute the
  			       calendar date from a given JD
	            date_conversion() includes test for legal date values
		    and notifies errors like 32 January.
 ****************************************************************/

double FAR PASCAL_CONV swe_julday(int year, int month, int day, double hour, int gregflag) 
{
  double jd;
  double u,u0,u1,u2;
  u = year;
  if (month < 3) u -=1;
  u0 = u + 4712.0;
  u1 = month + 1.0;
  if (u1 < 4) u1 += 12.0;
  jd = floor(u0*365.25)
     + floor(30.6*u1+0.000001)
     + day + hour/24.0 - 63.5;
  if (gregflag == SE_GREG_CAL) {
    u2 = floor(fabs(u) / 100) - floor(fabs(u) / 400);
    if (u < 0.0) u2 = -u2;
    jd = jd - u2 + 2;            
    if ((u < 0.0) && (u/100 == floor(u/100)) && (u/400 != floor(u/400)))
      jd -=1;
  }
  return jd;
}

/*** swe_revjul ******************************************************
  swe_revjul() is the inverse function to swe_julday(), see the description
  there.
  Arguments are julian day number, calendar flag (0=julian, 1=gregorian)
  return values are the calendar day, month, year and the hour of
  the day with decimal fraction (0 .. 23.999999).

  Be aware the we use astronomical year numbering for the years
  before Christ, not the historical year numbering.
  Astronomical years are done with negative numbers, historical
  years with indicators BC or BCE (before common era).
  Year  0 (astronomical)  	= 1 BC historical year
  year -1 (astronomical) 	= 2 BC historical year
  year -234 (astronomical) 	= 235 BC historical year
  etc.

  Original author Mark Pottenger, Los Angeles.
  with bug fix for year < -4711 16-aug-88 Alois Treindl
*************************************************************************/
void FAR PASCAL_CONV swe_revjul (double jd, int gregflag,
	     int *jyear, int *jmon, int *jday, double *jut)
{
  double u0,u1,u2,u3,u4;
  u0 = jd + 32082.5;
  if (gregflag == SE_GREG_CAL) {
    u1 = u0 + floor (u0/36525.0) - floor (u0/146100.0) - 38.0;
    if (jd >= 1830691.5) u1 +=1;
    u0 = u0 + floor (u1/36525.0) - floor (u1/146100.0) - 38.0;
  }
  u2 = floor (u0 + 123.0);
  u3 = floor ( (u2 - 122.2) / 365.25);
  u4 = floor ( (u2 - floor (365.25 * u3) ) / 30.6001);
  *jmon = (int) (u4 - 1.0);
  if (*jmon > 12) *jmon -= 12;
  *jday = (int) (u2 - floor (365.25 * u3) - floor (30.6001 * u4));
  *jyear = (int) (u3 + floor ( (u4 - 2.0) / 12.0) - 4800);
  *jut = (jd - floor (jd + 0.5) + 0.5) * 24.0;
}

/* transform local time to UTC or UTC to local time
 *
 * input 
 *   iyear ... dsec     date and time
 *   d_timezone		timezone offset
 * output
 *   iyear_out ... dsec_out
 * 
 * For time zones east of Greenwich, d_timezone is positive.
 * For time zones west of Greenwich, d_timezone is negative.
 * 
 * For conversion from local time to utc, use +d_timezone.
 * For conversion from utc to local time, use -d_timezone.
 */
void FAR PASCAL_CONV swe_utc_time_zone(
        int32 iyear, int32 imonth, int32 iday,
        int32 ihour, int32 imin, double dsec,
        double d_timezone,
        int32 *iyear_out, int32 *imonth_out, int32 *iday_out,
        int32 *ihour_out, int32 *imin_out, double *dsec_out
        )
{
  double tjd, d;
  AS_BOOL have_leapsec = FALSE;
  double dhour;
  if (dsec >= 60.0) {
    have_leapsec = TRUE;
    dsec -= 1.0;
  }
  dhour = ((double) ihour) + ((double) imin) / 60.0 + dsec / 3600.0;
  tjd = swe_julday(iyear, imonth, iday, 0, SE_GREG_CAL);
  dhour -= d_timezone;
  if (dhour < 0.0) {
    tjd -= 1.0;
    dhour += 24.0;
  }
  if (dhour >= 24.0) {
    tjd += 1.0;
    dhour -= 24.0;
  }
  swe_revjul(tjd + 0.001, SE_GREG_CAL, iyear_out, imonth_out, iday_out, &d);
  *ihour_out = (int) dhour;
  d = (dhour - (double) *ihour_out) * 60;
  *imin_out = (int) d;
  *dsec_out = (d - (double) *imin_out) * 60;
  if (have_leapsec)
    *dsec_out += 1.0;
}

/*
 * functions for the handling of UTC
 */

/* Leap seconds were inserted at the end of the following days:*/
#define NLEAP_SECONDS 24
#define NLEAP_SECONDS_SPACE 100
static int leap_seconds[NLEAP_SECONDS_SPACE] = {
19720630,
19721231,
19731231,
19741231,
19751231,
19761231,
19771231,
19781231,
19791231,
19810630,
19820630,
19830630,
19850630,
19871231,
19891231,
19901231,
19920630,
19930630,
19940630,
19951231,
19970630,
19981231,
20051231,
20081231,
0  /* keep this 0 as end mark */
};
#define J1972 2441317.5
#define NLEAP_INIT 10

/* Read additional leap second dates from external file, if given.
 */
static int init_leapsec(void)
{
  FILE *fp;
  int ndat, ndat_last;
  int tabsiz = 0;
  int i;
  char s[AS_MAXCH];
  char *sp;
  if (!init_leapseconds_done) {
    init_leapseconds_done = TRUE;
    tabsiz = NLEAP_SECONDS;
    ndat_last = leap_seconds[NLEAP_SECONDS - 1];
    /* no error message if file is missing */
    if ((fp = swi_fopen(-1, "seleapsec.txt", swed.ephepath, NULL)) == NULL)
      return NLEAP_SECONDS; 
    while(fgets(s, AS_MAXCH, fp) != NULL) {
      sp = s;
      while (*sp == ' ' || *sp == '\t') sp++;
        sp++;
      if (*sp == '#' || *sp == '\n')
        continue;
      ndat = atoi(s);
      if (ndat <= ndat_last)
        continue;
      /* table space is limited. no error msg, if exceeded */
      if (tabsiz >= NLEAP_SECONDS_SPACE)
        return tabsiz;
      leap_seconds[tabsiz] = ndat;
      tabsiz++;
    }
    if (tabsiz > NLEAP_SECONDS) leap_seconds[tabsiz] = 0; /* end mark */
    fclose(fp);
    return tabsiz;
  }
  /* find table size */
  tabsiz = 0;
  for (i = 0; i < NLEAP_SECONDS_SPACE; i++) {
    if (leap_seconds[i] == 0) 
      break;
    else
      tabsiz++;
  }
  return tabsiz;
}

/*
 * Input:  Clock time UTC, year, month, day, hour, minute, second (decimal).
 *         gregflag  Calendar flag
 *         serr      error string
 * Output: An array of doubles:
 *         dret[0] = Julian day number TT (ET)
 *         dret[1] = Julian day number UT1
 *
 * Function returns OK or Error.
 *
 * - Before 1972, swe_utc_to_jd() treats its input time as UT1.
 *   Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
 *   UTC second was regularly changed, so that UTC remained very close to UT1.
 * - From 1972 on, input time is treated as UTC.
 * - If delta_t - nleap - 32.184 > 1, the input time is treated as UT1.
 *   Note: Like this we avoid errors greater than 1 second in case that
 *   the leap seconds table (or the Swiss Ephemeris version) is not updated
 *   for a long time.
*/
int32 FAR PASCAL_CONV swe_utc_to_jd(int32 iyear, int32 imonth, int32 iday, int32 ihour, int32 imin, double dsec, int32 gregflag, double *dret, char *serr)
{
  double tjd_ut1, tjd_et, tjd_et_1972, dhour, d;
  int iyear2, imonth2, iday2;
  int i, j, ndat, nleap, tabsiz_nleap;
  /* 
   * error handling: invalid iyear etc. 
   */
  tjd_ut1 = swe_julday(iyear, imonth, iday, 0, gregflag);
  swe_revjul(tjd_ut1, gregflag, &iyear2, &imonth2, &iday2, &d);
  if (iyear != iyear2 || imonth != imonth2 || iday != iday2) {
    if (serr != NULL)
      sprintf(serr, "invalid date: year = %d, month = %d, day = %d", iyear, imonth, iday);
    return ERR;
  }
  if (ihour < 0 || ihour > 23 
   || imin < 0 || imin > 59 
   || dsec < 0 || dsec >= 61
   || (dsec >= 60 && (imin < 59 || ihour < 23 || tjd_ut1 < J1972))) {
    if (serr != NULL)
      sprintf(serr, "invalid time: %d:%d:%.2f", ihour, imin, dsec);
    return ERR;
  }
  dhour = (double) ihour + ((double) imin) / 60.0 + dsec / 3600.0;
  /* 
   * before 1972, we treat input date as UT1 
   */
  if (tjd_ut1 < J1972) {
    dret[1] = swe_julday(iyear, imonth, iday, dhour, gregflag);
    dret[0] = dret[1] + swe_deltat(dret[1]);
    return OK;
  }
  /* 
   * if gregflag = Julian calendar, convert to gregorian calendar 
   */
  if (gregflag == SE_JUL_CAL) {
    gregflag = SE_GREG_CAL;
    swe_revjul(tjd_ut1, gregflag, &iyear, &imonth, &iday, &d);
  }
  /* 
   * number of leap seconds since 1972: 
   */
  tabsiz_nleap = init_leapsec();
  nleap = NLEAP_INIT; /* initial difference between UTC and TAI in 1972 */
  ndat = iyear * 10000 + imonth * 100 + iday;
  for (i = 0; i < tabsiz_nleap; i++) {
    if (ndat <= leap_seconds[i])
      break;
    nleap++;
  }
  /*
   * For input dates > today:
   * If leap seconds table is not up to date, we'd better interpret the
   * input time as UT1, not as UTC. How do we find out? 
   * Check, if delta_t - nleap - 32.184 > 0.9
   */
  d = swe_deltat(tjd_ut1) * 86400.0;
  if (d - (double) nleap - 32.184 >= 1.0) {
    dret[1] = tjd_ut1 + dhour / 24.0;
    dret[0] = dret[1] + swe_deltat(dret[1]);
    return OK;
  }
  /* 
   * if input second is 60: is it a valid leap second ? 
   */
  if (dsec >= 60) {
    j = 0;
    for (i = 0; i < tabsiz_nleap; i++) {
      if (ndat == leap_seconds[i]) {
	j = 1;
	break;
      }
    }
    if (j != 1) {
      if (serr != NULL)
	sprintf(serr, "invalid time (no leap second!): %d:%d:%.2f", ihour, imin, dsec);
      return ERR;
    }
  }
  /* 
   * convert UTC to ET and UT1 
   */
  /* the number of days between input date and 1 jan 1972: */
  d = tjd_ut1 - J1972;
  /* SI time since 1972, ignoring leap seconds: */
  d += (double) ihour / 24.0 + (double) imin / 1440.0 + dsec / 86400.0; 
  /* ET (TT) */
  tjd_et_1972 = J1972 + (32.184 + NLEAP_INIT) / 86400.0;
  tjd_et = tjd_et_1972 + d + ((double) (nleap - NLEAP_INIT)) / 86400.0;
  d = swe_deltat(tjd_et);
  tjd_ut1 = tjd_et - swe_deltat(tjd_et - d);
  dret[0] = tjd_et;
  dret[1] = tjd_ut1;
  return OK;
}

/*
 * Input:  tjd_et   Julian day number, terrestrial time (ephemeris time).
 *         gregfalg Calendar flag
 * Output: UTC year, month, day, hour, minute, second (decimal).
 *
 * - Before 1 jan 1972 UTC, output UT1.
 *   Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
 *   UTC second was regularly changed, so that UTC remained very close to UT1.
 * - From 1972 on, output is UTC.
 * - If delta_t - nleap - 32.184 > 1, the output is UT1.
 *   Note: Like this we avoid errors greater than 1 second in case that
 *   the leap seconds table (or the Swiss Ephemeris version) has not been
 *   updated for a long time.
 */
void FAR PASCAL_CONV swe_jdet_to_utc(double tjd_et, int32 gregflag, int32 *iyear, int32 *imonth, int32 *iday, int32 *ihour, int32 *imin, double *dsec) 
{
  int i;
  int second_60 = 0;
  int iyear2, imonth2, iday2, nleap, ndat, tabsiz_nleap;
  double d, tjd, tjd_et_1972, tjd_ut, dret[10];
  /* 
   * if tjd_et is before 1 jan 1972 UTC, return UT1
   */
  tjd_et_1972 = J1972 + (32.184 + NLEAP_INIT) / 86400.0; 
  d = swe_deltat(tjd_et);
  tjd_ut = tjd_et - swe_deltat(tjd_et - d);
  if (tjd_et < tjd_et_1972) {
    swe_revjul(tjd_ut, gregflag, iyear, imonth, iday, &d);
    *ihour = (int32) d;
    d -= (double) *ihour;
    d *= 60;
    *imin = (int32) d;
    *dsec = (d - (double) *imin) * 60.0;
    return;
  }
  /* 
   * minimum number of leap seconds since 1972; we may be missing one leap
   * second
   */
  tabsiz_nleap = init_leapsec();
  swe_revjul(tjd_ut-1, SE_GREG_CAL, &iyear2, &imonth2, &iday2, &d);
  ndat = iyear2 * 10000 + imonth2 * 100 + iday2;
  nleap = 0; 
  for (i = 0; i < tabsiz_nleap; i++) {
    if (ndat <= leap_seconds[i])
      break;
    nleap++;
  }
  /* date of potentially missing leapsecond */
  if (nleap < tabsiz_nleap) {
    i = leap_seconds[nleap];
    iyear2 = i / 10000;
    imonth2 = (i % 10000) / 100;;
    iday2 = i % 100;
    tjd = swe_julday(iyear2, imonth2, iday2, 0, SE_GREG_CAL);
    swe_revjul(tjd+1, SE_GREG_CAL, &iyear2, &imonth2, &iday2, &d);
    swe_utc_to_jd(iyear2,imonth2,iday2, 0, 0, 0, SE_GREG_CAL, dret, NULL);
    d = tjd_et - dret[0];
    if (d >= 0) {
      nleap++;
    } else if (d < 0 && d > -1.0/86400.0) {
      second_60 = 1;
    }
  }
  /*
   * UTC, still unsure about one leap second
   */
  tjd = J1972 + (tjd_et - tjd_et_1972) - ((double) nleap + second_60) / 86400.0;
  swe_revjul(tjd, SE_GREG_CAL, iyear, imonth, iday, &d);
  *ihour = (int32) d;
  d -= (double) *ihour;
  d *= 60;
  *imin = (int32) d;
  *dsec = (d - (double) *imin) * 60.0 + second_60;
  /*
   * For input dates > today:
   * If leap seconds table is not up to date, we'd better interpret the
   * input time as UT1, not as UTC. How do we find out? 
   * Check, if delta_t - nleap - 32.184 > 0.9
   */
  d = swe_deltat(tjd_et);
  d = swe_deltat(tjd_et - d);
  if (d * 86400.0 - (double) (nleap + NLEAP_INIT) - 32.184 >= 1.0) {
    swe_revjul(tjd_et - d, SE_GREG_CAL, iyear, imonth, iday, &d);
    *ihour = (int32) d;
    d -= (double) *ihour;
    d *= 60;
    *imin = (int32) d;
    *dsec = (d - (double) *imin) * 60.0;
  }
  if (gregflag == SE_JUL_CAL) {
    tjd = swe_julday(*iyear, *imonth, *iday, 0, SE_GREG_CAL);
    swe_revjul(tjd, gregflag, iyear, imonth, iday, &d);
  }
}

/*
 * Input:  tjd_ut   Julian day number, universal time (UT1).
 *         gregfalg Calendar flag
 * Output: UTC year, month, day, hour, minute, second (decimal).
 *
 * - Before 1 jan 1972 UTC, output UT1.
 *   Note: UTC was introduced in 1961. From 1961 - 1971, the length of the
 *   UTC second was regularly changed, so that UTC remained very close to UT1.
 * - From 1972 on, output is UTC.
 * - If delta_t - nleap - 32.184 > 1, the output is UT1.
 *   Note: Like this we avoid errors greater than 1 second in case that
 *   the leap seconds table (or the Swiss Ephemeris version) has not been
 *   updated for a long time.
 */
void FAR PASCAL_CONV swe_jdut1_to_utc(double tjd_ut, int32 gregflag, int32 *iyear, int32 *imonth, int32 *iday, int32 *ihour, int32 *imin, double *dsec) 
{
  double tjd_et = tjd_ut + swe_deltat(tjd_ut);
  swe_jdet_to_utc(tjd_et, gregflag, iyear, imonth, iday, ihour, imin, dsec);
}