1 /***************************************************************************
2                           calendar.cpp  -  basic calendar transformations
3                              -------------------
4     begin                : 2017
5     copyright            : (C) 2017 by G. Duvert
6     email                : gilles.duvert@free.fr
7  ***************************************************************************/
8 
9 /***************************************************************************
10  *                                                                         *
11  *   This program is free software; you can redistribute it and/or modify  *
12  *   it under the terms of the GNU General Public License as published by  *
13  *   the Free Software Foundation; either version 2 of the License, or     *
14  *   (at your option) any later version.                                   *
15  *                                                                         *
16  ***************************************************************************/
17 
18 /* The intent of this file is to unify all calendar-related pieces of code used in various
19  * places of GDL. */
20 
21 #include "includefirst.hpp"
22 #include "calendar.hpp"
23 #include "gdlexception.hpp"
24 using namespace std;
25 
Gregorian2Julian(struct tm * ts)26   double Gregorian2Julian(struct tm *ts)
27   {
28     DDouble jd;
29     DLong day=ts->tm_mday;
30     DLong mon=ts->tm_mon+1;
31     DLong year=ts->tm_year+1900;
32     DLong hour=ts->tm_hour;
33     DLong min=ts->tm_min;
34     DDouble sec=ts->tm_sec;
35     if (!dateToJD(jd, day, mon, year, hour, min, sec))
36         throw GDLException("Invalid Julian date input.");
37     return jd;
38   }
dateToJD(DDouble & jd,DLong day,DLong month,DLong year,DLong hour,DLong minute,DDouble second)39   bool dateToJD(DDouble &jd, DLong day, DLong month, DLong year, DLong hour, DLong minute, DDouble second) {
40     if (year < -4716 || year > 5000000 || year == 0) return false;
41 
42     // the following tests seem to be NOT active in IDL. We choose to mimic it.
43     //    if (month < 1 || month > 12) return false;
44     //    if (day < 0 || day > 31) return false;
45     //these one too...
46     // if (hour < 0 || hour > 24) return false;
47     // if (minute < 0 || minute > 60) return false;
48     // if (second < 0 || second > 60) return false;
49 
50     DDouble a, y, b;
51     DLong m;
52     y = (year > 0) ? year : year + 1; //we use here a calendar with no year 0 (not astronomical)
53     m = month;
54     b = 0.0;
55     if (month <= 2) {
56       y = y - 1.0;
57       m = m + 12;
58     }
59     if (y >= 0) {
60       if (year > 1582 || (year == 1582 && (month > 10 ||
61         (month == 10 && day > 14)))) {
62         a = floor(y / 100.0);
63         b = 2.0 - a + floor(a / 4.0);
64       } else if (year == 1582 && month == 10 && day >= 5 && day <= 14) {
65         jd = 2299161; //date does not move
66         return true;
67       }
68     }
69     jd = floor(365.25 * y) + floor(30.6001 * (m + 1)) + day + (hour * 1.0) / 24.0 + (minute * 1.0) / 1440.0 +
70       (second * 1.0) / 86400.0 + 1720994.50 + b;
71     return true;
72   }
73 // C code ****************************************************
j2ymdhms(DDouble jd,DLong & iMonth,DLong & iDay,DLong & iYear,DLong & iHour,DLong & iMinute,DDouble & Second,DLong & dow,DLong & icap)74     bool j2ymdhms(DDouble jd, DLong &iMonth, DLong &iDay , DLong &iYear ,
75                   DLong &iHour , DLong &iMinute, DDouble &Second, DLong &dow, DLong &icap)
76     {
77     DDouble JD,Z,F;
78     DLong A,B,C,D,E;
79     JD = jd + 0.5;
80     Z = floor(JD);
81     if (Z < -1095 || Z > 1827933925 ) return FALSE;
82     F = JD - Z;
83     // note that IDL dow is false before Sun dec 31 12:00:00 -4714, (type: a=[-2,-1]& PRINT, FORMAT='(C())',a)
84     // ...and ... we are not!
85     if ((DLong)Z > 0) dow = ((DLong)Z) % 7; else dow = ((DLong)Z+1099) % 7; //just translate axis...
86 
87     if (Z < 2299161) A = (DLong)Z;
88     else {
89       DDouble a;
90       a = (DLong) ((Z - 1867216.25) / 36524.25);
91       A = (DLong) (Z + 1 + a - (DLong)(a / 4));
92     }
93 
94     B = A + 1524;
95     C = (DLong) ((B - 122.1) / 365.25);
96     D = (DLong) (365.25 * C);
97     E = (DLong) ((B - D) / 30.6001);
98 
99     // month
100     iMonth = E < 14 ? E - 1 : E - 13;
101     iMonth--; //to get a zero-based index;
102     // iday
103     iDay=B - D - (DLong)(30.6001 * E);
104 
105     // year
106 //    iYear = iMonth > 2 ? C - 4716 : C - 4715;
107     iYear = iMonth > 1 ? C - 4716 : C - 4715; //with a zero-based index
108     if (iYear < 1) iYear--; //No Year Zero
109     // hours
110     iHour = (DLong) (F * 24);
111     { //this prevents interpreting 04:00:00 as 03:59:60 !
112       //this kind of rounding up is explained in IDL doc.
113       DDouble FF=F+6E-10;
114       DLong test= (DLong) (FF * 24);
115       if (test > iHour) {iHour=test;F=FF;}
116     }
117 
118     icap = (iHour > 11);
119 
120     F -= (DDouble)iHour / 24;
121     // minutes
122     iMinute = (DLong) (F * 1440);
123     { //this prevents interpreting 04:00:00 as 03:59:60 !
124       //this kind of rounding up is explained in IDL doc.
125       DDouble FF=F+6E-10;
126       DLong test= (DLong) (FF * 1440);
127       if (test > iMinute) {iMinute=test;F=FF;}
128     }
129     F -= (DDouble)iMinute / (DDouble)1440;
130     // seconds
131     Second = F * 86400;
132     return TRUE;
133   }
134 
135 
136