home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / x / volume19 / xephem / part11 / cal_mjd.c next >
Encoding:
C/C++ Source or Header  |  1993-05-15  |  3.6 KB  |  188 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* given a date in months, mn, days, dy, years, yr,
  6.  * return the modified Julian date (number of days elapsed since 1900 jan 0.5),
  7.  * *mjd.
  8.  */
  9. void
  10. cal_mjd (mn, dy, yr, mjd)
  11. int mn, yr;
  12. double dy;
  13. double *mjd;
  14. {
  15.     static double last_mjd, last_dy;
  16.     static int last_mn, last_yr;
  17.     int b, d, m, y;
  18.     long c;
  19.  
  20.     if (mn == last_mn && yr == last_yr && dy == last_dy) {
  21.         *mjd = last_mjd;
  22.         return;
  23.     }
  24.  
  25.     m = mn;
  26.     y = (yr < 0) ? yr + 1 : yr;
  27.     if (mn < 3) {
  28.         m += 12;
  29.         y -= 1;
  30.     }
  31.  
  32.     if (yr < 1582 || yr == 1582 && (mn < 10 || mn == 10 && dy < 15)) 
  33.         b = 0;
  34.     else {
  35.         int a;
  36.         a = y/100;
  37.         b = 2 - a + a/4;
  38.     }
  39.  
  40.     if (y < 0)
  41.         c = (long)((365.25*y) - 0.75) - 694025L;
  42.     else
  43.         c = (long)(365.25*y) - 694025L;
  44.  
  45.     d = 30.6001*(m+1);
  46.  
  47.     *mjd = b + c + d + dy - 0.5;
  48.  
  49.     last_mn = mn;
  50.     last_dy = dy;
  51.     last_yr = yr;
  52.     last_mjd = *mjd;
  53. }
  54.  
  55. /* given the modified Julian date (number of days elapsed since 1900 jan 0.5,),
  56.  * mjd, return the calendar date in months, *mn, days, *dy, and years, *yr.
  57.  */
  58. void
  59. mjd_cal (mjd, mn, dy, yr)
  60. double mjd;
  61. int *mn, *yr;
  62. double *dy;
  63. {
  64.     static double last_mjd, last_dy;
  65.     static int last_mn, last_yr;
  66.     double d, f;
  67.     double i, a, b, ce, g;
  68.  
  69.     if (mjd == last_mjd) {
  70.         *mn = last_mn;
  71.         *yr = last_yr;
  72.         *dy = last_dy;
  73.         return;
  74.     }
  75.  
  76.     d = mjd + 0.5;
  77.     i = floor(d);
  78.     f = d-i;
  79.     if (f == 1) {
  80.         f = 0;
  81.         i += 1;
  82.     }
  83.  
  84.     if (i > -115860.0) {
  85.         a = floor((i/36524.25)+.9983573)+14;
  86.         i += 1 + a - floor(a/4.0);
  87.     }
  88.  
  89.     b = floor((i/365.25)+.802601);
  90.     ce = i - floor((365.25*b)+.750001)+416;
  91.     g = floor(ce/30.6001);
  92.     *mn = g - 1;
  93.     *dy = ce - floor(30.6001*g)+f;
  94.     *yr = b + 1899;
  95.  
  96.     if (g > 13.5)
  97.         *mn = g - 13;
  98.     if (*mn < 2.5)
  99.         *yr = b + 1900;
  100.     if (*yr < 1)
  101.         *yr -= 1;
  102.  
  103.     last_mn = *mn;
  104.     last_dy = *dy;
  105.     last_yr = *yr;
  106.     last_mjd = mjd;
  107. }
  108.  
  109. /* given an mjd, set *dow to 0..6 according to which day of the week it falls
  110.  * on (0=sunday) or set it to -1 if can't figure it out.
  111.  */
  112. void
  113. mjd_dow (mjd, dow)
  114. double mjd;
  115. int *dow;
  116. {
  117.     /* cal_mjd() uses Gregorian dates on or after Oct 15, 1582.
  118.      * (Pope Gregory XIII dropped 10 days, Oct 5..14, and improved the leap-
  119.      * year algorithm). however, Great Britian and the colonies did not
  120.      * adopt it until Sept 14, 1752 (they dropped 11 days, Sept 3-13,
  121.      * due to additional accumulated error). leap years before 1752 thus
  122.      * can not easily be accounted for from the cal_mjd() number...
  123.      */
  124.     if (mjd < -53798.5) {
  125.         /* pre sept 14, 1752 too hard to correct */
  126.         *dow = -1;
  127.         return;
  128.     }
  129.     *dow = ((long)floor(mjd-.5) + 1) % 7;/* 1/1/1900 (mjd 0.5) is a Monday*/
  130.     if (*dow < 0)
  131.         *dow += 7;
  132. }
  133.  
  134. /* given a mjd, return the the number of days in the month.  */
  135. void
  136. mjd_dpm (mjd, ndays)
  137. double mjd;
  138. int *ndays;
  139. {
  140.     static short dpm[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};
  141.     int m, y;
  142.     double d;
  143.  
  144.     mjd_cal (mjd, &m, &d, &y);
  145.     *ndays = (m==2 && ((y%4==0 && y%100!=0)||y%400==0)) ? 29 : dpm[m-1];
  146. }
  147.  
  148. /* given a mjd, return the year as a double. */
  149. void
  150. mjd_year (mjd, yr)
  151. double mjd;
  152. double *yr;
  153. {
  154.     static double last_mjd, last_yr;
  155.     int m, y;
  156.     double d;
  157.     double e0, e1;    /* mjd of start of this year, start of next year */
  158.  
  159.     if (mjd == last_mjd) {
  160.         *yr = last_yr;
  161.         return;
  162.     }
  163.  
  164.     mjd_cal (mjd, &m, &d, &y);
  165.     if (y == -1) y = -2;
  166.     cal_mjd (1, 1.0, y, &e0);
  167.     cal_mjd (1, 1.0, y+1, &e1);
  168.     *yr = y + (mjd - e0)/(e1 - e0);
  169.  
  170.     last_mjd = mjd;
  171.     last_yr = *yr;
  172. }
  173.  
  174. /* given a decimal year, return mjd */
  175. void
  176. year_mjd (y, mjd)
  177. double y;
  178. double *mjd;
  179. {
  180.     double e0, e1;    /* mjd of start of this year, start of next year */
  181.     int yf = floor (y);
  182.     if (yf == -1) yf = -2;
  183.  
  184.     cal_mjd (1, 1.0, yf, &e0);
  185.     cal_mjd (1, 1.0, yf+1, &e1);
  186.     *mjd = e0 + (y - yf)*(e1-e0);
  187. }
  188.