tutCalendar2_cal.c - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tutCalendar2_cal.c (18607B)
---
1 /*
2 The CalCalcs routines, a set of C-language routines to perform
3 calendar calculations.
4
5 Version 1.01, released 9 January 2010
6
7 Copyright (C) 2010 David W. Pierce, dpierce@ucsd.edu
8
9 This program is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation, either version 3 of the License, or
12 (at your option) any later version.
13
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with this program. If not, see <http://www.gnu.org/licenses/>.
21 */
22
23 #include <stdio.h>
24 #include <unistd.h>
25 #include <stdlib.h>
26 #include <string.h>
27
28 #include "udunits2.h"
29 #include "calcalcs.h"
30
31 #include "utCalendar2_cal.h"
32
33 static int have_initted=0;
34 static calcalcs_cal *cal_std=NULL;
35
36 /* Following controls the rounding precision of the routines. I.e., if we end up with a value
37 such as 59.99999999 seconds, we are going to round it to 60 seconds. The value is given
38 in seconds, so, for example, 1.e-3 means round up to 1 second if the value is 0.999 seconds or greater,
39 and 1.e-6 means round up to 1 second if the value is 0.999999 seconds or greater.
40 */
41 static const double sec_rounding_value = 1.e-8;
42
43
44 /* Internal to this file only */
45 static void initialize( ut_system *units_system );
46 static void get_origin( ut_unit *dataunits, int *y0, int *mon0, int *d0, int *h0, int *min0, double *s0 );
47 static cv_converter *get_day_to_user_converter( ut_unit *uu, int y0, int mon0, int d0, int h0, int min0, double s0 );
48 static cv_converter *get_user_to_day_converter( ut_unit *uu, int y0, int mon0, int d0, int h0, int min0, double s0 );
49 static calcalcs_cal *getcal( const char *name );
50 static void unknown_cal_emit_warning( const char *calendar_name );
51
52 /* Stores previuosly initialized calendars and their names */
53 static const int maxcals_known=100;
54 static int ncals_known=0;
55 static calcalcs_cal **known_cal; /* ptr to array of calcals_cal ptrs */
56 static char **known_cal_name;
57
58 /* Stores previously emitted "unknown calendar" warnings */
59 #define UTC2_MAX_UNKCAL_WARNS 1000
60 static char *unknown_cal_emitted_warning_for[ UTC2_MAX_UNKCAL_WARNS ];
61 static int n_unkcal=0;
62
63 /*========================================================================================
64 * Turns the passed value into a Y/M/D date
65 */
66 int utCalendar2_cal( double val, ut_unit *dataunits, int *year, int *month, int *day, int *hour,
67 int *minute, double *second, const char *calendar_name )
68 {
69
70 int jdnew, ndinc, ierr, iorig, iround;
71 double fdays, extra_seconds, tot_extra_seconds;
72 int ndays;
73 calcalcs_cal *cal2use;
74
75 /* Following vars are saved between invocations and reused if the
76 * passed units are the same as last time.
77 */
78 static ut_unit *prev_units=NULL;
79 static cv_converter *conv_user_units_to_days=NULL;
80 static int y0, mon0, d0, h0, min0, jday0;
81 static double s0, extra_seconds0;
82 static char *prev_calendar=NULL;
83
84 if( dataunits == NULL ) {
85 fprintf( stderr, "Error, utCalendar2 passed a NULL units\n" );
86 return( UT_ENOINIT );
87 }
88
89 if( have_initted == 0 )
90 initialize( ut_get_system(dataunits) );
91
92 /* Get the calendar we will be using, based on the passed name
93 */
94 cal2use = getcal( calendar_name );
95 if( cal2use == NULL ) {
96 unknown_cal_emit_warning( calendar_name );
97 cal2use = getcal( "Standard" );
98 }
99
100 /* See if we are being passed the same units and calendar as last time. If so,
101 * we can optimize by not recomputing all this junk
102 */
103 if( (prev_units != NULL) && (prev_calendar != NULL)
104 && (strcmp(prev_calendar,cal2use->name)==0)
105 && (ut_compare( prev_units, dataunits ) == 0)) {
106 /* Units and calendar are same as used last call */
107 }
108 else
109 {
110 /* Units/calendar combo are different from previously saved units, must redo calcs */
111
112 if( prev_units != NULL )
113 ut_free( prev_units );
114
115 if( prev_calendar != NULL )
116 free( prev_calendar );
117
118 /* Get origin day of the data units */
119 get_origin( dataunits, &y0, &mon0, &d0, &h0, &min0, &s0 ); /* Note: static vars */
120
121 /* Number of seconds into the specified origin day */
122 extra_seconds0 = h0*3600.0 + min0*60.0 + s0; /* Note: static vars */
123
124 /* Convert the origin day to Julian Day number in the specified calendar */
125 if( (ierr = ccs_date2jday( cal2use, y0, mon0, d0, &jday0 )) != 0 ) {
126 fprintf( stderr, "Error in utCalendar2: %s\n", ccs_err_str(ierr) );
127 return( UT_EINVALID );
128 }
129
130 /* Get converter from user-specified units to "days" */
131 if( conv_user_units_to_days != NULL )
132 cv_free( conv_user_units_to_days );
133 conv_user_units_to_days = get_user_to_day_converter( dataunits, y0, mon0, d0, h0, min0, s0 );
134
135 /* Save these units so we can reuse our time-consuming
136 * calculations next time if they are the same units
137 */
138 prev_units = ut_clone( dataunits );
139 if( ut_compare( prev_units, dataunits ) != 0 ) {
140 fprintf( stderr, "error, internal error in udunits2 library found in routine utCalendar2: a clone of the user's units does not equal the original units!\n" );
141 exit(-1);
142 }
143
144 prev_calendar = (char *)malloc( sizeof(char) * (strlen(cal2use->name)+1 ));
145 strcpy( prev_calendar, cal2use->name );
146 }
147
148 /* Convert user value of offset to floating point days */
149 fdays = cv_convert_double( conv_user_units_to_days, val );
150
151 /* Get integer number of days and seconds past that */
152 ndays = fdays;
153 extra_seconds = (fdays - ndays)*86400.0;
154
155 /* Get new Julian day */
156 jdnew = jday0 + ndays;
157
158 /* Handle the sub-day part */
159 tot_extra_seconds = extra_seconds0 + extra_seconds;
160 ndinc = tot_extra_seconds / 86400.0;
161 jdnew += ndinc;
162 tot_extra_seconds -= ndinc*86400.0;
163 if( tot_extra_seconds < 0.0 ) {
164 tot_extra_seconds += 86400.0;
165 jdnew--;
166 }
167
168 /* Convert to a date */
169 if( (ierr = ccs_jday2date( cal2use, jdnew, year, month, day )) != 0 ) {
170 fprintf( stderr, "Error in utCalendar2: %s\n", ccs_err_str(ierr) );
171 return( UT_EINVALID );
172 }
173
174 *hour = tot_extra_seconds / 3600.0;
175 tot_extra_seconds -= *hour * 3600.0;
176 *minute = tot_extra_seconds / 60.0;
177 tot_extra_seconds -= *minute * 60.0;
178 *second = tot_extra_seconds;
179
180 /* Handle the rouding issues */
181 iorig = *second; /* Integer conversion */
182 iround = *second + sec_rounding_value;
183 if( iround > iorig ) {
184 /* printf( "rounding alg invoked, orig date: %04d-%02d-%02d %02d:%02d:%.20lf\n", *year, *month, *day, *hour, *minute, *second ); */
185 *second = (double)iround;
186 if( *second >= 60.0 ) {
187 *second -= 60.0;
188 *minute += 1.0;
189 if( *minute >= 60.0 ) {
190 *minute -= 60.0;
191 *hour += 1.0;
192 if( *hour >= 24.0 ) {
193 *hour -= 24.0;
194 if( (ierr = ccs_jday2date( cal2use, jdnew+1, year, month, day )) != 0 ) {
195 fprintf( stderr, "Error in utCalendar2: %s\n", ccs_err_str(ierr) );
196 return( UT_EINVALID );
197 }
198 }
199 }
200 }
201 /* printf( "after rounding alg here is new date: %04d-%02d-%02d %02d:%02d:%02.20lf\n", *year, *month, *day, *hour, *minute, *second ); */
202 }
203
204 return(0);
205 }
206
207 /*========================================================================================
208 * Turn the passed Y/M/D date into a value in the user's units
209 */
210 int utInvCalendar2_cal( int year, int month, int day, int hour, int minute, double second,
211 ut_unit *user_unit, double *value, const char *calendar_name )
212 {
213 int jday, ierr, diff_in_days;
214 double fdiff_in_days, val_days, val_partdays, fdiff_in_partdays, fpartday;
215 calcalcs_cal *cal2use;
216
217 /* Following vars are static and retained between invocations for efficiency */
218 static ut_unit *prev_units=NULL;
219 static int y0, mon0, d0, h0, min0, jday0;
220 static double s0, fpartday0;
221 static cv_converter *conv_days_to_user_units=NULL;
222 static char *prev_calendar=NULL;
223
224 if( have_initted == 0 )
225 initialize( ut_get_system(user_unit) );
226
227 /* Get the calendar we will be using, based on the passed name
228 */
229 cal2use = getcal( calendar_name );
230 if( cal2use == NULL ) {
231 unknown_cal_emit_warning( calendar_name );
232 cal2use = getcal( "Standard" );
233 }
234
235 if( (prev_units != NULL) && (prev_calendar != NULL)
236 && (strcmp(prev_calendar,cal2use->name)==0)
237 && (ut_compare( prev_units, user_unit ) == 0)) {
238 /* Units are same as used last call */
239 }
240 else
241 {
242 if( prev_units != NULL )
243 ut_free( prev_units );
244
245 if( prev_calendar != NULL )
246 free( prev_calendar );
247
248 if( conv_days_to_user_units != NULL )
249 cv_free( conv_days_to_user_units );
250
251 /* Get origin day of the data units */
252 get_origin( user_unit, &y0, &mon0, &d0, &h0, &min0, &s0 ); /* Note: static vars */
253
254 /* Convert the origin day to Julian Day number in the specified calendar */
255 if( (ierr = ccs_date2jday( cal2use, y0, mon0, d0, &jday0 )) != 0 ) {
256 fprintf( stderr, "Error in utCalendar2: %s\n", ccs_err_str(ierr) );
257 return( UT_EINVALID );
258 }
259
260 /* Get the origin's HMS in fractional (floating point) part of a Julian day */
261 fpartday0 = (double)h0/24.0 + (double)min0/1440.0 + s0/86400.0;
262
263 /* Get converter for turning days into user's units */
264 conv_days_to_user_units = get_day_to_user_converter( user_unit, y0, mon0, d0, h0, min0, s0 );
265
266 /* Save these units so we can reuse our time-consuming
267 * calculations next time if they are the same units
268 */
269 prev_units = ut_clone( user_unit );
270 if( ut_compare( prev_units, user_unit ) != 0 ) {
271 fprintf( stderr, "error, internal error in udunits2 library found in routine utInvCalendar2: a clone of the user's units does not equal the original units!\n" );
272 exit(-1);
273 }
274
275 prev_calendar = (char *)malloc( sizeof(char) * (strlen(cal2use->name)+1 ));
276 strcpy( prev_calendar, cal2use->name );
277 }
278
279 /* Turn passed date into a Julian day */
280 if( (ierr = ccs_date2jday( cal2use, year, month, day, &jday )) != 0 ) {
281 fprintf( stderr, "Error in utInvCalendar2: %s\n", ccs_err_str(ierr) );
282 return( UT_EINVALID );
283 }
284
285 /* jday and jday0 can be very large and nearly equal, so we difference
286 * them first to keep the precision high
287 */
288 diff_in_days = jday - jday0;
289 fdiff_in_days = (double)diff_in_days;
290
291 /* Get the fractional (floating point) part of a Julian day difference
292 */
293 fpartday = (double)hour/24.0 + (double)minute/1440.0 + second/86400.0;
294 fdiff_in_partdays = fpartday - fpartday0;
295
296 /* Convert days and partial days to user's units */
297 val_days = cv_convert_double( conv_days_to_user_units, fdiff_in_days );
298 val_partdays = cv_convert_double( conv_days_to_user_units, fdiff_in_partdays );
299
300 /* Hopefully this will minimize the roundoff errors */
301 *value = val_days + val_partdays;
302
303 return(0);
304 }
305
306 /*==============================================================================================
307 * Get a converter that turns the user's units into days
308 */
309 static cv_converter *get_user_to_day_converter( ut_unit *uu, int y0, int mon0, int d0, int h0, int min0, double s0 )
310 {
311 char daystr[1024];
312 ut_unit *udu_days;
313 cv_converter *conv_user_units_to_days;
314
315 sprintf( daystr, "days since %04d-%02d-%02d %02d:%02d:%f",
316 y0, mon0, d0, h0, min0, s0 );
317
318 udu_days = ut_parse( ut_get_system(uu), daystr, UT_ASCII );
319 if( udu_days == NULL ) {
320 fprintf( stderr, "internal error in utCalendar2/conv_to_days: failed to parse following string: \"%s\"\n",
321 daystr );
322 exit(-1);
323 }
324 conv_user_units_to_days = ut_get_converter( uu, udu_days );
325 if( conv_user_units_to_days == NULL ) {
326 fprintf( stderr, "internal error in utCalendar2/conv_to_days: cannot convert from \"%s\" to user units\n",
327 daystr );
328 exit(-1);
329 }
330
331 ut_free( udu_days );
332 return( conv_user_units_to_days );
333 }
334
335 /*==============================================================================================
336 * Get a converter that turns days into the user's units
337 */
338 static cv_converter *get_day_to_user_converter( ut_unit *uu, int y0, int mon0, int d0, int h0, int min0, double s0 )
339 {
340 char daystr[1024];
341 ut_unit *udu_days;
342 cv_converter *conv_days_to_user_units;
343
344 sprintf( daystr, "days since %04d-%02d-%02d %02d:%02d:%f",
345 y0, mon0, d0, h0, min0, s0 );
346
347 udu_days = ut_parse( ut_get_system(uu), daystr, UT_ASCII );
348 if( udu_days == NULL ) {
349 fprintf( stderr, "internal error in utCalendar2/conv_to_user_units: failed to parse following string: \"%s\"\n",
350 daystr );
351 exit(-1);
352 }
353 conv_days_to_user_units = ut_get_converter( udu_days, uu );
354 if( conv_days_to_user_units == NULL ) {
355 fprintf( stderr, "internal error in utCalendar2/conv_to_user_units: cannot convert from user units to \"%s\"\n",
356 daystr );
357 exit(-1);
358 }
359
360 ut_free( udu_days );
361 return( conv_days_to_user_units );
362 }
363
364 /*==========================================================================================
365 * The user specified some origin to the time units. For example, if the units string
366 * were "days since 2005-10-15", then the origin date is 2005-10-15. This routine
367 * deduces the specified origin date from the passed units structure
368 */
369 static void get_origin( ut_unit *dataunits, int *y0, int *mon0, int *d0, int *h0, int *min0, double *s0 )
370 {
371 double s0lib, rez, tval, tval_conv, resolution;
372 cv_converter *conv_user_date_to_ref_date;
373 ut_unit *tmpu;
374 int y0lib, mon0lib, d0lib, h0lib, min0lib;
375 char ustr[1024];
376
377 static ut_unit *udu_ref_date=NULL; /* Saved between invocations */
378
379 if( udu_ref_date == NULL ) {
380
381 /* Make a timestamp units that refers to the udunits2 library's intrinsic
382 * time origin. The first thing we do is parse a timestampe unit
383 * (any old timestamp unit) and immediately discard the result, to account
384 * for a bug in the udunits-2 library that fails to set the library's
385 * time origin unless this step is done first. Without this, a call to
386 * ut_decode_time with tval==0 returns year=-4713 rather than 2001.
387 */
388 if( (tmpu = ut_parse( ut_get_system(dataunits), "days since 2010-01-09", UT_ASCII )) == NULL ) {
389 fprintf( stderr, "Internal error in routnie utCalendar2/get_origin: failed to parse temp timestamp string\n" );
390 exit(-1);
391 }
392 ut_free( tmpu );
393
394 tval = 0.0;
395 ut_decode_time( tval, &y0lib, &mon0lib, &d0lib, &h0lib, &min0lib, &s0lib, &rez );
396 sprintf( ustr, "seconds since %04d-%02d-%02d %02d:%02d:%f",
397 y0lib, mon0lib, d0lib, h0lib, min0lib, s0lib );
398 udu_ref_date = ut_parse( ut_get_system(dataunits), ustr, UT_ASCII );
399 if( udu_ref_date == NULL ) {
400 fprintf( stderr, "internal error in routine utCalendar2/get_origin: could not parse origin string \"%s\"\n",
401 ustr );
402 exit(-1);
403 }
404 }
405
406 /* Get converter from passed units to the library's intrinsic time units */
407 conv_user_date_to_ref_date = ut_get_converter( dataunits, udu_ref_date );
408
409 /* convert tval=0 to ref date */
410 tval = 0.0;
411 tval_conv = cv_convert_double( conv_user_date_to_ref_date, tval );
412
413 /* Now decode the converted value */
414 ut_decode_time( tval_conv, y0, mon0, d0, h0, min0, s0, &resolution );
415
416 cv_free( conv_user_date_to_ref_date );
417 }
418
419 /*========================================================================================*/
420 static void initialize( ut_system *units_system )
421 {
422 int i;
423
424 /* Make space for our saved calendars */
425 known_cal = (calcalcs_cal **)malloc( sizeof( calcalcs_cal * ) * maxcals_known );
426 if( known_cal == NULL ) {
427 fprintf( stderr, "Error in utCalendar2 routines, could not allocate internal storage\n" );
428 exit(-1);
429 }
430 for( i=0; i<maxcals_known; i++ )
431 known_cal[i] = NULL;
432 known_cal_name = (char **)malloc( sizeof( char * ) * maxcals_known );
433 if( known_cal_name == NULL ) {
434 fprintf( stderr, "Error in utCalendar2 routines, could not allocate internal storage for calendar names\n" );
435 exit(-1);
436 }
437 for( i=0; i<maxcals_known; i++ )
438 known_cal_name[i] = NULL;
439
440 /* Make a standard calendar */
441 cal_std = ccs_init_calendar( "Standard" );
442
443 have_initted = 1; /* a global */
444 }
445
446 /*========================================================================================
447 * Returns NULL if the passed calendar name is both not found and not creatable
448 */
449 static calcalcs_cal *getcal( const char *name )
450 {
451 int i, new_index;
452 calcalcs_cal *new_cal;
453
454 if( cal_std == NULL ) {
455 fprintf( stderr, "Coding error in utCalendar2_cal routines, cal_std is null!\n" );
456 exit(-1);
457 }
458
459 if( strcasecmp( name, "standard" ) == 0 )
460 return( cal_std );
461
462 /* See if it is one of the previously known calendars */
463 for( i=0; i<ncals_known; i++ ) {
464 if( strcmp( known_cal_name[i], name ) == 0 )
465 return( known_cal[i] );
466 }
467
468 /* If we get here, the cal is not known, so create it if possible */
469 new_cal = ccs_init_calendar( name );
470 if( new_cal == NULL )
471 return( NULL ); /* unknown name */
472
473 /* We now know a new calendar */
474 new_index = ncals_known;
475 ncals_known++;
476
477 /* new_index is where we will be storing the new calendar */
478 if( ncals_known > maxcals_known ) {
479 ncals_known = maxcals_known;
480 new_index = strlen(name); /* some arbitrary value */
481 if( new_index >= maxcals_known )
482 new_index = 10;
483 }
484
485 /* If there was already a calendar stored in this slot
486 * (because we might be reusing slots) then clear it out
487 */
488 if( known_cal[new_index] != NULL )
489 ccs_free_calendar( known_cal[new_index] );
490
491 if( known_cal_name[new_index] != NULL )
492 free( known_cal_name[new_index] );
493
494 known_cal[new_index] = new_cal;
495 known_cal_name[new_index] = (char *)malloc( sizeof(char) * (strlen(name)+1 ));
496 strcpy( known_cal_name[new_index], name );
497
498 return( new_cal );
499 }
500
501 /*=============================================================================================
502 */
503 static void unknown_cal_emit_warning( const char *calendar_name )
504 {
505 int i;
506
507 for( i=0; i<n_unkcal; i++ ) {
508 if( strcmp( calendar_name, unknown_cal_emitted_warning_for[i] ) == 0 )
509 /* Already emitted a warning for this calendar */
510 return;
511 }
512
513 if( n_unkcal == UTC2_MAX_UNKCAL_WARNS )
514 /* Too many warnings already given, give up */
515 return;
516
517 /* OK, emit the warning now */
518 fprintf( stderr, "Error in utCalendar2_cal/utInvCalendar2_cal: unknown calendar \"%s\". Using Standard calendar instead\n",
519 calendar_name );
520
521 /* Save the fact that we have warned for this calendar */
522 unknown_cal_emitted_warning_for[ n_unkcal ] = (char *)malloc( sizeof(char) * (strlen(calendar_name)+1 ));
523 if( unknown_cal_emitted_warning_for[ n_unkcal ] == NULL ) {
524 fprintf( stderr, "Error in utCalendar_cal: could not allocate internal memory to store string \"%s\"\n",
525 calendar_name );
526 return;
527 }
528
529 strcpy( unknown_cal_emitted_warning_for[ n_unkcal ], calendar_name );
530 n_unkcal++;
531 }