source: XIOS/trunk/extern/src_netcdf4/nctime.c @ 409

Last change on this file since 409 was 409, checked in by ymipsl, 11 years ago

Add improved nectdf internal library src

YM

  • Property svn:eol-style set to native
File size: 29.7 KB
Line 
1/*********************************************************************
2 *   Copyright 2008, University Corporation for Atmospheric Research
3 *   See netcdf/COPYRIGHT file for copying and redistribution conditions.
4 *   $Id: nctime.c,v 1.9 2010/05/05 22:15:39 dmh Exp $
5 *********************************************************************/
6
7/*
8 * This code was extracted with permission from the CDMS time
9 * conversion and arithmetic routines developed by Bob Drach, Lawrence
10 * Livermore National Laboratory as part of the cdtime library.  Russ
11 * Rew of the UCAR Unidata Program made changes and additions to
12 * support the "-t" option of the netCDF ncdump utility, including a
13 * 366-day climate calendar.
14 *
15 * For the complete time conversion and climate calendar facilities of
16 * the CDMS library, get the original sources from LLNL.
17 */
18
19#include <stdlib.h>
20#include <stdio.h>
21#include <ctype.h>
22#include <math.h>
23#include <string.h>
24#include <stdarg.h>
25#include <assert.h>
26#include "nctime.h"
27
28static int cuErrOpts;                        /* Error options */
29static int cuErrorOccurred = 0;              /* True iff cdError was called */
30
31#define CU_FATAL 1                           /* Exit immediately on fatal error */
32#define CU_VERBOSE 2                         /* Report errors */
33
34#define CD_DEFAULT_BASEYEAR "1979"           /* Default base year for relative time (no 'since' clause) */
35#define VALCMP(a,b) ((a)<(b)?-1:(b)<(a)?1:0)
36
37/* forward declarations */
38static void cdComp2Rel(cdCalenType timetype, cdCompTime comptime, char* relunits, double* reltime);
39static void cdRel2CompMixed(double reltime, cdUnitTime unit, cdCompTime basetime, cdCompTime *comptime);
40static void cdRel2Comp(cdCalenType timetype, char* relunits, double reltime, cdCompTime* comptime);
41
42/* Trim trailing whitespace, up to n characters. */
43/* If no whitespace up to the last character, set */
44/* the last character to null, else set the first */
45/* whitespace character to null. */
46static void
47cdTrim(char* s, int n)
48{
49        char* c;
50
51        if(s==NULL)
52                return;
53        for(c=s; *c && c<s+n-1 && !isspace(*c); c++);
54        *c='\0';
55        return;
56}
57
58static
59void cdError(char *fmt, ...){
60        va_list args;
61       
62        cuErrorOccurred = 1;
63        if(cuErrOpts & CU_VERBOSE){
64                va_start(args,fmt);
65                fprintf(stderr, "CDMS error: ");
66                vfprintf(stderr, fmt, args);
67                fprintf(stderr, "\n");
68                va_end(args);
69        }
70        if(cuErrOpts & CU_FATAL)
71                exit(1);
72        return;
73}
74
75#define ISLEAP(year,timeType)   ((timeType & Cd366) || (((timeType) & CdHasLeap) && (!((year) % 4) && (((timeType) & CdJulianType) || (((year) % 100) || !((year) % 400))))))
76static int mon_day_cnt[12] = {31,28,31,30,31,30,31,31,30,31,30,31};
77static int days_sum[12] = {0,31,59,90,120,151,181,212,243,273,304,334};
78
79/* Compute month and day from year and day-of-year.
80 *
81 *      Input:
82 *              doy          (int)  (day-of-year)
83 *              date->year   (long)  (year since 0 BC)
84 *              date->timeType (CdTimetype) (time type)
85 *              date->baseYear   base year for relative times
86 *      Output:
87 *              date->month  (short)  (month in year)
88 *              date->day    (short)  (day in month)
89 *
90 *
91 * Derived from NRL NEONS V3.6.
92 */
93
94static void
95CdMonthDay(int *doy, CdTime *date)
96{
97        int i;                          /* month counter */
98        int idoy;                       /* day of year counter */
99        long year;
100
101        if ((idoy = *doy) < 1) {
102                date->month = 0;
103                date->day   = 0;
104                return;
105        }
106
107        if(!(date->timeType & CdChronCal))   /* Ignore year for Clim calendar */
108                year = 0;
109        else if(!(date->timeType & CdBase1970)) /* year is offset from base for relative time */
110                year = date->baseYear + date->year;
111        else
112                year = date->year;
113
114        if (ISLEAP(year,date->timeType)) {
115                mon_day_cnt[1] = 29;
116        } else {
117                mon_day_cnt[1] = 28;
118        }
119        date->month     = 0;
120        for (i = 0; i < 12; i++) {
121                (date->month)++;
122                date->day       = idoy;
123                if ((idoy -= ((date->timeType & Cd365) ? (mon_day_cnt[date->month-1]) : 30)) <= 0) {
124                        return;
125                }
126        }
127        return;
128}
129
130/* Compute day-of-year from year, month and day
131 *
132 *      Input:
133 *              date->year  (long)  (year since 0 BC)
134 *              date->month (short)  (month in year)
135 *              date->day   (short)  (day in month)
136 *              date->baseYear   base year for relative times
137 *      Output: doy         (int)  (day-of-year)
138 *
139 * Derived from NRL NEONS V3.6
140 */
141
142static void
143CdDayOfYear(CdTime *date, int *doy)
144{
145        int leap_add = 0;               /* add 1 day if leap year */
146        int month;                      /* month */
147        long year;
148
149        month   = date->month;
150        if (month < 1 || month > 12) {
151                cdError( "Day-of-year error; month: %d\n", month);
152                month = 1;     
153        }
154
155        if(!(date->timeType & CdChronCal))   /* Ignore year for Clim calendar */
156                year = 0;
157        else if(!(date->timeType & CdBase1970)) /* year is offset from base for relative time */
158                year = date->baseYear + date->year;
159        else
160                year = date->year;
161
162        if (ISLEAP(year,date->timeType) && month > 2) leap_add = 1;
163        if( ((date->timeType) & Cd365) || ((date->timeType) & Cd366) ) {
164            *doy = days_sum[month-1] + date->day + leap_add ;
165        } else {                /* date->timeType & Cd360 */
166            *doy = 30*(month-1) + date->day + leap_add ;
167        }
168        return;
169}
170
171/* Convert epochal time (hours since 00 jan 1, 1970)
172 *   to human time (structured)
173 *
174 * Input:
175 *   etime = epochal time representation
176 *   timeType = time type (e.g., CdChron, CdClim, etc.) as defined in cdms.h
177 *   baseYear = base real, used for relative time types only
178 *
179 * Output: htime = human (structured) time representation
180 *
181 * Derived from NRL Neons V3.6
182 */
183void
184Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime)
185{
186        long    ytemp;                  /* temporary year holder */
187        int     yr_day_cnt;             /* count of days in year */
188        int     doy;                    /* day of year */
189        int     daysInLeapYear;              /* number of days in a leap year */
190        int     daysInYear;                  /* days in non-leap year */
191        extern void CdMonthDay(int *doy, CdTime *date);
192
193        doy     = (long) floor(etime / 24.) + 1;
194        htime->hour     = etime - (double) (doy - 1) * 24.;
195
196                                             /* Correct for goofy floor func on J90 */
197        if(htime->hour >= 24.){
198                doy += 1;
199                htime->hour -= 24.;
200        }
201
202        htime->baseYear = (timeType & CdBase1970) ? 1970 : baseYear;
203        if(!(timeType & CdChronCal)) htime->baseYear = 0; /* Set base year to 0 for Clim */
204        if(timeType & Cd366) {
205            daysInLeapYear = 366;
206            daysInYear = 366;
207        } else {
208            daysInLeapYear = (timeType & Cd365) ? 366 : 360;
209            daysInYear = (timeType & Cd365) ? 365 : 360;
210        }
211
212        if (doy > 0) {
213                for (ytemp = htime->baseYear; ; ytemp++) {
214                        yr_day_cnt = ISLEAP(ytemp,timeType) ? daysInLeapYear : daysInYear;
215                        if (doy <= yr_day_cnt) break;
216                        doy -= yr_day_cnt;
217                }
218        } else {
219                for (ytemp = htime->baseYear-1; ; ytemp--) {
220                        yr_day_cnt = ISLEAP(ytemp,timeType) ? daysInLeapYear : daysInYear;
221                        doy += yr_day_cnt;
222                        if (doy > 0) break;
223                }
224        }
225        htime->year = (timeType & CdBase1970) ? ytemp : (ytemp - htime->baseYear);
226        if(!(timeType & CdChronCal)) htime->year = 0; /* Set year to 0 for Clim */
227        htime->timeType = timeType;
228        CdMonthDay(&doy,htime);
229
230        return;
231}
232
233/* Add 'nDel' times 'delTime' to epochal time 'begEtm',
234 * return the result in epochal time 'endEtm'.
235 */
236static void
237CdAddDelTime(double begEtm, long nDel, CdDeltaTime delTime, CdTimeType timeType,
238             long baseYear, double *endEtm)
239{
240        double delHours;
241        long delMonths, delYears;
242        CdTime bhtime, ehtime;
243
244        extern void Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime);
245        extern void Cdh2e(CdTime *htime, double *etime);
246
247        switch(delTime.units){
248          case CdYear:
249                delMonths = 12;
250                break;
251          case CdSeason:
252                delMonths = 3;
253                break;
254          case CdMonth:
255                delMonths = 1;
256                break;
257          case CdWeek:
258                delHours = 168.0;
259                break;
260          case CdDay:
261                delHours = 24.0;
262                break;
263          case CdHour:
264                delHours = 1.0;
265                break;
266          case CdMinute:
267                delHours = 1./60.;
268                break;
269          case CdSecond:
270                delHours = 1./3600.;
271                break;
272          default:
273                cdError("Invalid delta time units: %d\n",delTime.units);
274                return;
275        }
276
277        switch(delTime.units){
278          case CdYear: case CdSeason: case CdMonth:
279                Cde2h(begEtm,timeType,baseYear,&bhtime);
280                delMonths = delMonths * nDel * delTime.count + bhtime.month - 1;
281                delYears = (delMonths >= 0 ? (delMonths/12) : (delMonths+1)/12 - 1);
282                ehtime.year = bhtime.year + delYears;
283                ehtime.month = delMonths - (12 * delYears) + 1;
284                ehtime.day = 1;
285                ehtime.hour = 0.0;
286                ehtime.timeType = timeType;
287                ehtime.baseYear = !(timeType & CdChronCal) ? 0 :
288                        (timeType & CdBase1970) ? 1970 : baseYear; /* base year is 0 for Clim, */
289                                                                   /* 1970 for Chron, */
290                                                                   /* or input base year for Rel */
291                Cdh2e(&ehtime,endEtm);
292                break;
293          case CdWeek: case CdDay: case CdHour: case CdMinute: case CdSecond:
294                delHours *= (nDel * delTime.count);
295                *endEtm = begEtm + delHours;
296                break;
297          default: break;
298        }
299        return;
300}
301
302/* Parse relative units, returning the unit and base component time. */
303/* Function returns 1 if error, 0 on success */
304int
305cdParseRelunits(cdCalenType timetype, char* relunits, cdUnitTime* unit, cdCompTime* base_comptime)
306{
307        char charunits[CD_MAX_RELUNITS];
308        char basetime_1[CD_MAX_CHARTIME];
309        char basetime_2[CD_MAX_CHARTIME];
310        char basetime[CD_MAX_CHARTIME];
311        int nconv1, nconv2, nconv;
312
313                                             /* Parse the relunits */
314        /* Allow ISO-8601 "T" date-time separator as well as blank separator */
315        nconv1 = sscanf(relunits,"%s since %[^T]T%s",charunits,basetime_1,basetime_2);
316        if(nconv1==EOF || nconv1==0){
317                cdError("Error on relative units conversion, string = %s\n",relunits);
318                return 1;
319        }
320        nconv2 = sscanf(relunits,"%s since %s %s",charunits,basetime_1,basetime_2);
321        if(nconv2==EOF || nconv2==0){
322                cdError("Error on relative units conversion, string = %s\n",relunits);
323                return 1;
324        }
325        if(nconv1 < nconv2) {
326            nconv = nconv2;
327        } else {
328            nconv = sscanf(relunits,"%s since %[^T]T%s",charunits,basetime_1,basetime_2);
329        }
330                                             /* Get the units */
331        cdTrim(charunits,CD_MAX_RELUNITS);
332        if(!strncmp(charunits,"sec",3) || !strcmp(charunits,"s")){
333                *unit = cdSecond;
334        }
335        else if(!strncmp(charunits,"min",3) || !strcmp(charunits,"mn")){
336                *unit = cdMinute;
337        }
338        else if(!strncmp(charunits,"hour",4) || !strcmp(charunits,"hr")){
339                *unit = cdHour;
340        }
341        else if(!strncmp(charunits,"day",3) || !strcmp(charunits,"dy")){
342                *unit = cdDay;
343        }
344        else if(!strncmp(charunits,"week",4) || !strcmp(charunits,"wk")){
345                *unit = cdWeek;
346        }
347        else if(!strncmp(charunits,"month",5) || !strcmp(charunits,"mo")){
348                *unit = cdMonth;
349        }
350        else if(!strncmp(charunits,"season",6)){
351                *unit = cdSeason;
352        }
353        else if(!strncmp(charunits,"year",4) || !strcmp(charunits,"yr")){
354                if(!(timetype & cdStandardCal)){
355                        cdError("Error on relative units conversion: climatological units cannot be 'years'.\n");
356                        return 1;
357                }
358                *unit = cdYear;
359        }
360        else {
361                cdError("Error on relative units conversion: invalid units = %s\n",charunits);
362                return 1;
363        }
364
365                                             /* Build the basetime, if any (default is 1979), */
366                                             /* or month 1 for climatological time. */
367        if(nconv == 1){
368                if(timetype & cdStandardCal)
369                        strcpy(basetime,CD_DEFAULT_BASEYEAR);
370                else
371                        strcpy(basetime,"1");
372        }
373                                             /* Convert the basetime to component, then epochal (hours since 1970) */
374        else{
375                if(nconv == 2){
376                        cdTrim(basetime_1,CD_MAX_CHARTIME);
377                        strcpy(basetime,basetime_1);
378                }
379                else{
380                        cdTrim(basetime_1,CD_MAX_CHARTIME);
381                        cdTrim(basetime_2,CD_MAX_CHARTIME);
382                        sprintf(basetime,"%s %s",basetime_1,basetime_2);
383                }
384        }
385
386        cdChar2Comp(timetype, basetime, base_comptime);
387
388        return 0;
389}
390
391/* ca - cb in Gregorian calendar */
392/* Result is in hours. */
393static double
394cdDiffGregorian(cdCompTime ca, cdCompTime cb){
395
396        double rela, relb;
397
398        cdComp2Rel(cdStandard, ca, "hours", &rela);
399        cdComp2Rel(cdStandard, cb, "hours", &relb);
400        return (rela - relb);
401}
402
403/* Return -1, 0, 1 as ca is less than, equal to, */
404/* or greater than cb, respectively. */
405static int 
406cdCompCompare(cdCompTime ca, cdCompTime cb){
407
408        int test;
409
410        if ((test = VALCMP(ca.year, cb.year)))
411                return test;
412        else if ((test = VALCMP(ca.month, cb.month)))
413                return test;
414        else if ((test = VALCMP(ca.day, cb.day)))
415                return test;
416        else
417                return (VALCMP(ca.hour, cb.hour));
418}
419
420/* ca - cb in Julian calendar.  Result is in hours. */
421static double
422cdDiffJulian(cdCompTime ca, cdCompTime cb){
423
424        double rela, relb;
425
426        cdComp2Rel(cdJulian, ca, "hours", &rela);
427        cdComp2Rel(cdJulian, cb, "hours", &relb);
428        return (rela - relb);
429}
430
431/* ca - cb in mixed Julian/Gregorian calendar. */
432/* Result is in hours. */
433static double
434cdDiffMixed(cdCompTime ca, cdCompTime cb){
435
436        static cdCompTime ZA = {1582, 10, 5, 0.0};
437        static cdCompTime ZB = {1582, 10, 15, 0.0};
438        double result;
439
440        if (cdCompCompare(cb, ZB) == -1){
441                if (cdCompCompare(ca, ZB) == -1) {
442                        result = cdDiffJulian(ca, cb);
443                }
444                else {
445                        result = cdDiffGregorian(ca, ZB) + cdDiffJulian(ZA, cb);
446                }
447        }
448        else {
449                if (cdCompCompare(ca, ZB) == -1){
450                        result = cdDiffJulian(ca, ZA) + cdDiffGregorian(ZB, cb);
451                }
452                else {
453                        result = cdDiffGregorian(ca, cb);
454                }
455        }
456        return result;
457}
458
459/* Divide ('endEtm' - 'begEtm') by 'delTime',
460 * return the integer portion of the result in 'nDel'.
461 */
462static void
463CdDivDelTime(double begEtm, double endEtm, CdDeltaTime delTime, CdTimeType timeType,
464             long baseYear, long *nDel)
465{
466        double delHours, frange;
467        long delMonths, range;
468        CdTime bhtime, ehtime;
469        int hoursInYear;
470       
471        extern void Cde2h(double etime, CdTimeType timeType, long baseYear, CdTime *htime);
472
473        switch(delTime.units){
474          case CdYear:
475                delMonths = 12;
476                break;
477          case CdSeason:
478                delMonths = 3;
479                break;
480          case CdMonth:
481                delMonths = 1;
482                break;
483          case CdWeek:
484                delHours = 168.0;
485                break;
486          case CdDay:
487                delHours = 24.0;
488                break;
489          case CdHour:
490                delHours = 1.0;
491                break;
492          case CdMinute:
493                delHours = 1./60.;
494                break;
495          case CdSecond:
496                delHours = 1./3600.;
497                break;
498          default:
499                cdError("Invalid delta time units: %d\n",delTime.units);
500                return;
501        }
502
503        switch(delTime.units){
504          case CdYear: case CdSeason: case CdMonth:
505                delMonths *= delTime.count;
506                Cde2h(begEtm,timeType,baseYear,&bhtime);
507                Cde2h(endEtm,timeType,baseYear,&ehtime);
508                if(timeType & CdChronCal){   /* Chron and Rel time */
509                        range = 12*(ehtime.year - bhtime.year)
510                                + (ehtime.month - bhtime.month);
511                }
512                else{                        /* Clim time, ignore year */
513                        range = (ehtime.month - bhtime.month);
514                        if(range < 0) range += 12;
515                }
516                *nDel = abs(range)/delMonths;
517                break;
518          case CdWeek: case CdDay: case CdHour: case CdMinute: case CdSecond:
519                delHours *= (double)delTime.count;
520                if(timeType & CdChronCal){   /* Chron and Rel time */
521                        frange = fabs(endEtm - begEtm);
522                }
523                else{                        /* Clim time, ignore year, but */
524                                             /* wraparound relative to hours-in-year*/
525                        frange = endEtm - begEtm;
526                        if(timeType & Cd366) {
527                            hoursInYear = 8784;
528                        } else {
529                            hoursInYear = (timeType & Cd365) ? 8760. : 8640.;
530                        }
531                                             /* Normalize frange to interval [0,hoursInYear) */
532                        if(frange < 0.0 || frange >= hoursInYear)
533                                frange -= hoursInYear * floor(frange/hoursInYear);
534                }
535                *nDel = (frange + 1.e-10*delHours)/delHours;
536                break;
537            default: break;
538        }
539        return;
540}
541
542/* Value is in hours. Translate to units. */
543static double
544cdFromHours(double value, cdUnitTime unit){
545        double result;
546
547        switch(unit){
548        case cdSecond:
549                result = value * 3600.0;
550                break;
551        case cdMinute:
552                result = value * 60.0;
553                break;
554        case cdHour:
555                result = value;
556                break;
557        case cdDay:
558                result = value/24.0;
559                break;
560        case cdWeek:
561                result = value/168.0;
562                break;
563        case cdMonth:
564        case cdSeason:
565        case cdYear:
566        case cdFraction:
567        default:
568                cdError("Error on conversion from hours to vague unit");
569                result = 0;
570                break;
571        }
572        return result;
573}
574                                             /* Map to old timetypes */
575static int
576cdToOldTimetype(cdCalenType newtype, CdTimeType* oldtype)
577{
578        switch(newtype){
579          case cdStandard:
580                *oldtype = CdChron;
581                break;
582          case cdJulian:
583                *oldtype = CdJulianCal;
584                break;
585          case cdNoLeap:
586                *oldtype = CdChronNoLeap;
587                break;
588          case cd360:
589                *oldtype = CdChron360;
590                break;
591          case cd366:
592                *oldtype = CdChron366;
593                break;
594          case cdClim:
595                *oldtype = CdClim;
596                break;
597          case cdClimLeap:
598                *oldtype = CdClimLeap;
599                break;
600          case cdClim360:
601                *oldtype = CdClim360;
602                break;
603          default:
604                cdError("Error on relative units conversion, invalid timetype = %d",newtype);
605                return 1;
606        }
607        return 0;
608}
609
610/* Convert human time to epochal time (hours since 00 jan 1, 1970)
611 *
612 * Input: htime = human time representation
613 *
614 * Output: etime = epochal time representation
615 *
616 * Derived from NRL Neons V3.6
617 */
618void
619Cdh2e(CdTime *htime, double *etime)
620{
621        long    ytemp, year;                    /* temporary year holder */
622        int     day_cnt;                /* count of days */
623        int     doy;                    /* day of year */
624        long    baseYear;                    /* base year for epochal time */
625        int     daysInLeapYear;              /* number of days in a leap year */
626        int     daysInYear;                  /* days in non-leap year */
627        extern void CdDayOfYear(CdTime *date, int *doy);
628
629        CdDayOfYear(htime,&doy);
630       
631        day_cnt = 0;
632
633        baseYear = ((htime->timeType) & CdBase1970) ? 1970 : htime->baseYear;
634        year = ((htime->timeType) & CdBase1970) ? htime->year : (htime->year + htime->baseYear);
635        if(!((htime->timeType) & CdChronCal)) baseYear = year = 0;      /* set year and baseYear to 0 for Clim */
636        if((htime->timeType) & Cd366) {
637            daysInLeapYear = 366;
638            daysInYear = 366;
639        } else {
640            daysInLeapYear = ((htime->timeType) & Cd365) ? 366 : 360;
641            daysInYear = ((htime->timeType) & Cd365) ? 365 : 360;
642        }
643       
644        if (year > baseYear) {
645                for (ytemp = year - 1; ytemp >= baseYear; ytemp--) {
646                        day_cnt += ISLEAP(ytemp,htime->timeType) ? daysInLeapYear : daysInYear;
647                }
648        } else if (year < baseYear) {
649                for (ytemp = year; ytemp < baseYear; ytemp++) {
650                        day_cnt -= ISLEAP(ytemp,htime->timeType) ? daysInLeapYear : daysInYear;
651                }
652        }       
653        *etime  = (double) (day_cnt + doy - 1) * 24. + htime->hour;
654        return;
655}
656
657/* Validate the component time, return 0 if valid, 1 if not */
658static int
659cdValidateTime(cdCalenType timetype, cdCompTime comptime)
660{
661        if(comptime.month<1 || comptime.month>12){
662                cdError("Error on time conversion: invalid month = %hd\n",comptime.month);
663                return 1;
664        }
665        if(comptime.day<1 || comptime.day>31){
666                cdError("Error on time conversion: invalid day = %hd\n",comptime.day);
667                return 1;
668        }
669        if(comptime.hour<0.0 || comptime.hour>24.0){
670                cdError("Error on time conversion: invalid hour = %lf\n",comptime.hour);
671                return 1;
672        }
673        return 0;
674}
675
676void
677cdChar2Comp(cdCalenType timetype, char* chartime, cdCompTime* comptime)
678{
679        double sec;
680        int ihr, imin, nconv;
681        long year;
682        short day;
683        short month;
684
685        comptime->year = CD_NULL_YEAR;
686        comptime->month = CD_NULL_MONTH;
687        comptime->day = CD_NULL_DAY;
688        comptime->hour = CD_NULL_HOUR;
689       
690        if(timetype & cdStandardCal){
691                nconv = sscanf(chartime,"%ld-%hd-%hd %d:%d:%lf",&year,&month,&day,&ihr,&imin,&sec);
692                if(nconv==EOF || nconv==0){
693                        cdError("Error on character time conversion, string = %s\n",chartime);
694                        return;
695                }
696                if(nconv >= 1){
697                        comptime->year = year;
698                }
699                if(nconv >= 2){
700                        comptime->month = month;
701                }
702                if(nconv >= 3){
703                        comptime->day = day;
704                }
705                if(nconv >= 4){
706                        if(ihr<0 || ihr>23){
707                                cdError("Error on character time conversion: invalid hour = %d\n",ihr);
708                                return;
709                        }
710                        comptime->hour = (double)ihr;
711                }
712                if(nconv >= 5){
713                        if(imin<0 || imin>59){
714                                cdError("Error on character time conversion: invalid minute = %d\n",imin);
715                                return;
716                        }
717                        comptime->hour += (double)imin/60.;
718                }
719                if(nconv >= 6){
720                        if(sec<0.0 || sec>60.0){
721                                cdError("Error on character time conversion: invalid second = %lf\n",sec);
722                                return;
723                        }
724                        comptime->hour += sec/3600.;
725                }
726        }
727        else{                                /* Climatological */
728                nconv = sscanf(chartime,"%hd-%hd %d:%d:%lf",&month,&day,&ihr,&imin,&sec);
729                if(nconv==EOF || nconv==0){
730                        cdError("Error on character time conversion, string = %s",chartime);
731                        return;
732                }
733                if(nconv >= 1){
734                        comptime->month = month;
735                }
736                if(nconv >= 2){
737                        comptime->day = day;
738                }
739                if(nconv >= 3){
740                        if(ihr<0 || ihr>23){
741                                cdError("Error on character time conversion: invalid hour = %d\n",ihr);
742                                return;
743                        }
744                        comptime->hour = (double)ihr;
745                }
746                if(nconv >= 4){
747                        if(imin<0 || imin>59){
748                                cdError("Error on character time conversion: invalid minute = %d\n",imin);
749                                return;
750                        }
751                        comptime->hour += (double)imin/60.;
752                }
753                if(nconv >= 5){
754                        if(sec<0.0 || sec>60.0){
755                                cdError("Error on character time conversion: invalid second = %lf\n",sec);
756                                return;
757                        }
758                        comptime->hour += sec/3600.;
759                }
760        }
761        (void)cdValidateTime(timetype,*comptime);
762        return;
763}
764
765/* Convert ct to relunits (unit, basetime) */
766/* in the mixed Julian/Gregorian calendar. */
767/* unit is anything but year, season, month. unit and basetime are */
768/* from the parsed relunits. Return result in reltime. */
769static void
770cdComp2RelMixed(cdCompTime ct, cdUnitTime unit, cdCompTime basetime, double *reltime){
771
772        double hourdiff;
773
774        hourdiff = cdDiffMixed(ct, basetime);
775        *reltime = cdFromHours(hourdiff, unit);
776        return;
777}
778
779static void
780cdComp2Rel(cdCalenType timetype, cdCompTime comptime, char* relunits, double* reltime)
781{
782        cdCompTime base_comptime;
783        CdDeltaTime deltime;
784        CdTime humantime;
785        CdTimeType old_timetype;
786        cdUnitTime unit;
787        double base_etm, etm, delta;
788        long ndel, hoursInYear;
789       
790                                             /* Parse the relunits */
791        if(cdParseRelunits(timetype, relunits, &unit, &base_comptime))
792                return;
793
794                                             /* Handle mixed Julian/Gregorian calendar */
795        if (timetype == cdMixed){
796                switch(unit){
797                case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
798                        cdComp2RelMixed(comptime, unit, base_comptime, reltime);
799                        return;
800                case cdYear: case cdSeason: case cdMonth:
801                        timetype = cdStandard;
802                        break;
803                case cdFraction:
804                        cdError("invalid unit in conversion");
805                        break;
806                default: break;
807                }
808        }
809       
810                                             /* Convert basetime to epochal */
811        humantime.year = base_comptime.year;
812        humantime.month = base_comptime.month;
813        humantime.day = base_comptime.day;
814        humantime.hour = base_comptime.hour;
815        humantime.baseYear = 1970;
816                                             /* Map to old-style timetype */
817        if(cdToOldTimetype(timetype,&old_timetype))
818                return;
819        humantime.timeType = old_timetype;
820        Cdh2e(&humantime,&base_etm);
821
822                                             /* Map end time to epochal */
823        humantime.year = comptime.year;
824        humantime.month = comptime.month;
825        humantime.day = comptime.day;
826        humantime.hour = comptime.hour;
827        Cdh2e(&humantime,&etm);
828                                             /* Calculate relative time value for months or hours */
829        deltime.count = 1;
830        deltime.units = (CdTimeUnit)unit;
831        switch(unit){
832          case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
833                delta = etm - base_etm;
834                if(!(timetype & cdStandardCal)){        /* Climatological time */
835                        hoursInYear = (timetype & cd365Days) ? 8760. : (timetype & cdHasLeap) ? 8784. : 8640.;
836                                             /* Normalize delta to interval [0,hoursInYear) */
837                        if(delta < 0.0 || delta >= hoursInYear)
838                                delta -= hoursInYear * floor(delta/hoursInYear);
839                }
840                break;
841          case cdYear: case cdSeason: case cdMonth:
842                CdDivDelTime(base_etm, etm, deltime, old_timetype, 1970, &ndel);
843                break;
844          case cdFraction:
845                cdError("invalid unit in conversion");
846                break;
847          default: break;
848        }
849
850                                             /* Convert to output units */
851        switch(unit){
852          case cdSecond:
853                *reltime = 3600.0 * delta;
854                break;
855          case cdMinute:
856                *reltime = 60.0 * delta;
857                break;
858          case cdHour:
859                *reltime = delta;
860                break;
861          case cdDay:
862                *reltime = delta/24.0;
863                break;
864          case cdWeek:
865                *reltime = delta/168.0;
866                break;
867          case cdMonth: case cdSeason: case cdYear: /* Already in correct units */
868                if(timetype & cdStandardCal)
869                        *reltime = (base_etm <= etm) ? (double)ndel : (double)(-ndel);
870                else                         /* Climatological time is already normalized*/
871                        *reltime = (double)ndel;
872                break;
873          default:
874                cdError("invalid unit in conversion");
875                break;
876        }
877
878        return;
879}
880
881/* Add (value,unit) to comptime. */
882/* value is in hours. */
883/* calendar is anything but cdMixed. */
884static void
885cdCompAdd(cdCompTime comptime, double value, cdCalenType calendar, cdCompTime *result){
886
887        double reltime;
888
889        cdComp2Rel(calendar, comptime, "hours", &reltime);
890        reltime += value;
891        cdRel2Comp(calendar, "hours", reltime, result);
892        return;
893}
894
895/* Add value in hours to ct, in the mixed Julian/Gregorian
896 * calendar. */
897static void
898cdCompAddMixed(cdCompTime ct, double value, cdCompTime *result){
899
900        static cdCompTime ZA = {1582, 10, 5, 0.0};
901        static cdCompTime ZB = {1582, 10, 15, 0.0};
902        double xj, xg;
903
904        if (cdCompCompare(ct, ZB) == -1){
905                xj = cdDiffJulian(ZA, ct);
906                if (value <= xj){
907                        cdCompAdd(ct, value, cdJulian, result);
908                }
909                else {
910                        cdCompAdd(ZB, value-xj, cdStandard, result);
911                }
912        }
913        else {
914                xg = cdDiffGregorian(ZB, ct);
915                if (value > xg){
916                        cdCompAdd(ct, value, cdStandard, result);
917                }
918                else {
919                        cdCompAdd(ZA, value-xg, cdJulian, result);
920                }
921        }
922        return;
923}
924
925/* Return value expressed in hours. */
926static double
927cdToHours(double value, cdUnitTime unit){
928
929        double result = 0;
930
931        switch(unit){
932        case cdSecond:
933                result = value/3600.0;
934                break;
935        case cdMinute:
936                result = value/60.0;
937                break;
938        case cdHour:
939                result = value;
940                break;
941        case cdDay:
942                result = 24.0 * value;
943                break;
944        case cdWeek:
945                result = 168.0 * value;
946                break;
947        default:
948                cdError("invalid unit in conversion");
949                break;
950
951        }
952        return result;
953}
954
955/* Convert relative time (reltime, unit, basetime) to comptime in the
956 * mixed Julian/Gregorian calendar. unit is anything but year, season,
957 * month. unit and basetime are from the parsed relunits. Return
958 * result in comptime. */
959static void
960cdRel2CompMixed(double reltime, cdUnitTime unit, cdCompTime basetime, cdCompTime *comptime){
961
962        reltime = cdToHours(reltime, unit);
963        cdCompAddMixed(basetime, reltime, comptime);
964        return;
965}
966
967
968static void
969cdRel2Comp(cdCalenType timetype, char* relunits, double reltime, cdCompTime* comptime)
970{
971        CdDeltaTime deltime;
972        CdTime humantime;
973        CdTimeType old_timetype;
974        cdCompTime base_comptime;
975        cdUnitTime unit, baseunits;
976        double base_etm, result_etm;
977        double delta;
978        long idelta;
979
980                                             /* Parse the relunits */
981        if(cdParseRelunits(timetype, relunits, &unit, &base_comptime))
982                return;
983
984        if (timetype == cdMixed){
985                switch(unit){
986                case cdWeek: case cdDay: case cdHour: case cdMinute: case cdSecond:
987                        cdRel2CompMixed(reltime, unit, base_comptime, comptime);
988                        return;
989                case cdYear: case cdSeason: case cdMonth:
990                        timetype = cdStandard;
991                        break;
992                case cdFraction:
993                        cdError("invalid unit in conversion");
994                        break;
995                default: break;
996                }
997        }
998
999        baseunits =cdBadUnit;
1000        switch(unit){
1001          case cdSecond:
1002                delta = reltime/3600.0;
1003                baseunits = cdHour;
1004                break;
1005          case cdMinute:
1006                delta = reltime/60.0;
1007                baseunits = cdHour;
1008                break;
1009          case cdHour:
1010                delta = reltime;
1011                baseunits = cdHour;
1012                break;
1013          case cdDay:
1014                delta = 24.0 * reltime;
1015                baseunits = cdHour;
1016                break;
1017          case cdWeek:
1018                delta = 168.0 * reltime;
1019                baseunits = cdHour;
1020                break;
1021          case cdMonth:
1022                idelta = (long)(reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1023                baseunits = cdMonth;
1024                break;
1025          case cdSeason:
1026                idelta = (long)(3.0 * reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1027                baseunits = cdMonth;
1028                break;
1029          case cdYear:
1030                idelta = (long)(12 * reltime + (reltime<0 ? -1.e-10 : 1.e-10));
1031                baseunits = cdMonth;
1032                break;
1033          default:
1034                cdError("invalid unit in conversion");
1035                break;
1036        }
1037
1038        deltime.count = 1;
1039        deltime.units = (CdTimeUnit)baseunits;
1040
1041        humantime.year = base_comptime.year;
1042        humantime.month = base_comptime.month;
1043        humantime.day = base_comptime.day;
1044        humantime.hour = base_comptime.hour;
1045        humantime.baseYear = 1970;
1046                                             /* Map to old-style timetype */
1047        if(cdToOldTimetype(timetype,&old_timetype))
1048                return;
1049        humantime.timeType = old_timetype;
1050
1051        Cdh2e(&humantime,&base_etm);
1052                                             /* If months, seasons, or years, */
1053        if(baseunits == cdMonth){
1054
1055                                             /* Calculate new epochal time from integer months. */
1056                                             /* Convert back to human, then comptime. */
1057                                             /* For zero reltime, just return the basetime*/
1058                if(reltime != 0.0){
1059                        CdAddDelTime(base_etm,idelta,deltime,old_timetype,1970,&result_etm);
1060                        Cde2h(result_etm, old_timetype, 1970, &humantime);
1061                }
1062        }
1063                                             /* Calculate new epochal time. */
1064                                             /* Convert back to human, then comptime. */
1065        else if(baseunits == cdHour){
1066                Cde2h(base_etm+delta, old_timetype, 1970, &humantime);
1067               
1068        }
1069        comptime->year = humantime.year;
1070        comptime->month = humantime.month;
1071        comptime->day = humantime.day;
1072        comptime->hour = humantime.hour;
1073
1074        return;
1075}
1076
1077/* rkr: output as ISO 8601 strings */
1078static void
1079cdComp2Iso(cdCalenType timetype, int separator, cdCompTime comptime, char* time)
1080{
1081        double dtmp, sec;
1082        int ihr, imin, isec;
1083        int nskip;
1084
1085        if(cdValidateTime(timetype,comptime))
1086                return;
1087       
1088        ihr = (int)comptime.hour;
1089        dtmp = 60.0 * (comptime.hour - (double)ihr);
1090        imin = (int)dtmp;
1091        sec = 60.0 * (dtmp - (double)imin);
1092        isec = sec;
1093
1094        if(sec == isec)
1095            if(isec == 0)
1096                if(imin == 0)
1097                    if(ihr == 0)
1098                        nskip = 4;
1099                    else
1100                        nskip = 3;
1101                else
1102                    nskip = 2;
1103            else
1104                nskip = 1;
1105        else
1106            nskip = 0;
1107
1108        if(timetype & cdStandardCal){
1109            switch (nskip) {
1110            case 0:             /* sec != 0 && (int)sec != sec */
1111                sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d:%lf",
1112                        comptime.year,comptime.month,comptime.day,separator,ihr,imin,sec);
1113                break;
1114            case 1:
1115                sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d:%2.2d",
1116                        comptime.year,comptime.month,comptime.day,separator,ihr,imin,isec);
1117                break;
1118            case 2:
1119                sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d:%2.2d",
1120                        comptime.year,comptime.month,comptime.day,separator,ihr,imin);
1121                break;
1122            case 3:
1123                sprintf(time,"%4.4ld-%2.2hd-%2.2hd%c%2.2d",
1124                        comptime.year,comptime.month,comptime.day,separator,ihr);
1125                break;
1126            case 4:
1127                sprintf(time,"%4.4ld-%2.2hd-%2.2hd",
1128                        comptime.year,comptime.month,comptime.day);
1129                break;
1130            }
1131        }
1132        else {                               /* Climatological */
1133            switch (nskip) {
1134            case 0:             /* sec != 0 && (int)sec != sec */
1135                sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d:%lf",
1136                        comptime.month,comptime.day,separator,ihr,imin,sec);
1137                break;
1138            case 1:
1139                sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d:%2.2d",
1140                        comptime.month,comptime.day,separator,ihr,imin,isec);
1141                break;
1142            case 2:
1143                sprintf(time,"%2.2hd-%2.2hd%c%2.2d:%2.2d",
1144                        comptime.month,comptime.day,separator,ihr,imin);
1145                break;
1146            case 3:
1147                sprintf(time,"%2.2hd-%2.2hd%c%2.2d",
1148                        comptime.month,comptime.day,separator,ihr);
1149                break;
1150            case 4:
1151                sprintf(time,"%2.2hd-%2.2hd",
1152                        comptime.month,comptime.day);
1153                break;
1154            }
1155        }
1156        return;
1157}
1158
1159/* rkr: added for output closer to ISO 8601 */
1160void
1161cdRel2Iso(cdCalenType timetype, char* relunits, int separator, double reltime, char* chartime)
1162{
1163        cdCompTime comptime;
1164
1165        cdRel2Comp(timetype, relunits, reltime, &comptime);
1166        cdComp2Iso(timetype, separator, comptime, chartime);
1167
1168        return;
1169}
Note: See TracBrowser for help on using the repository browser.