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