From: <ai...@us...> - 2009-05-07 18:17:00
|
Revision: 9936 http://plplot.svn.sourceforge.net/plplot/?rev=9936&view=rev Author: airwin Date: 2009-05-07 18:16:53 +0000 (Thu, 07 May 2009) Log Message: ----------- Store TAI-UTC data in just one consolidated lookup table containing revised TAI_UTC structs which have members for time_secs_tai and time_secs_utc rather than just time_secs alone. The result builds and gives the same results for my (not yet committed) extended python example 29. Modified Paths: -------------- trunk/lib/qsastime/README.tai-utc trunk/lib/qsastime/qsastime.c trunk/lib/qsastime/tai-utc-gen.c Modified: trunk/lib/qsastime/README.tai-utc =================================================================== --- trunk/lib/qsastime/README.tai-utc 2009-05-07 12:43:52 UTC (rev 9935) +++ trunk/lib/qsastime/README.tai-utc 2009-05-07 18:16:53 UTC (rev 9936) @@ -23,10 +23,10 @@ The tai-utc-gen application reads in these (UTC) epoch, offset1, offset2, and slope data from tai-utc.dat, calculates the discontinuous change in UTC at the unique TAI epoch of the UTC discontinuity using equation 2 below, and -stores the results in the header file tai-utc.h with two tables which have -identical offset1, offset2, slope, and discontinuous change data, but which -have epochs of the discontinuity expressed in MJD(TAI) for -TAI_TO_UTC_lookup_table and MJD(UTC) for UTC_TO_TAI_lookup_table. +stores the results in the header file tai-utc.h as the TAI_UTC_lookup_table +with data for MJD(TAI) (int base_day, double time_sec_tai), MJD(UTC) (int +base_day, double time_sec_utc = 0.), the discontinuous change, offset1, +offset2, and slope stored for each epoch. The qsastime library uses equation 1 to calculate MJD(TAI) as a function of MJD(UTC), and Modified: trunk/lib/qsastime/qsastime.c =================================================================== --- trunk/lib/qsastime/qsastime.c 2009-05-07 12:43:52 UTC (rev 9935) +++ trunk/lib/qsastime/qsastime.c 2009-05-07 18:16:53 UTC (rev 9936) @@ -63,7 +63,8 @@ static const int MonthStartDOY_L[] = {0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335}; /* Static function declarations. */ -static int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2); +static int geMJDtime_TAI(const MJDtime *number1, const TAI_UTC *number2); +static int geMJDtime_UTC(const MJDtime *number1, const TAI_UTC *number2); static double leap_second_TAI(const MJDtime *MJD_TAI, int *inleap, int *index); /* End of static function declarations. */ @@ -913,7 +914,7 @@ return posn; } -int geMJDtimeTAI_UTC(const MJDtime *number1, const TAI_UTC *number2) { +int geMJDtime_TAI(const MJDtime *number1, const TAI_UTC *number2) { /* Returns true if number1 >= number2. */ /* N.B. both number1 and number2 must be normalized. */ if(number1->base_day > number2->base_day) { @@ -921,10 +922,22 @@ } else if (number1->base_day < number2->base_day) { return 0; } else { - return (number1->time_sec >= number2->time_sec); + return (number1->time_sec >= number2->time_sec_tai); } } +int geMJDtime_UTC(const MJDtime *number1, const TAI_UTC *number2) { + /* Returns true if number1 >= number2. */ + /* N.B. both number1 and number2 must be normalized. */ + if(number1->base_day > number2->base_day) { + return 1; + } else if (number1->base_day < number2->base_day) { + return 0; + } else { + return (number1->time_sec >= number2->time_sec_utc); + } +} + double leap_second_TAI(const MJDtime *MJD_TAI, int *inleap, int *index) { /* Logic assumes input MJD_TAI is in TAI */ /* *inleap lets the calling routine know whether MJD_TAI corresponds @@ -933,17 +946,17 @@ MJDtime MJD_value, *MJD = &MJD_value; double leap; int debug=0; - /* N.B. geMJDtimeTAI_UTC only works for normalized values. */ + /* N.B. geMJDtime_TAI only works for normalized values. */ *MJD = *MJD_TAI; normalize_MJD(MJD); - /* Search for index such that TAI_TO_UTC_lookup_table[*index] <= MJD(TAI) < TAI_TO_UTC_lookup_table[*index+1] */ - bhunt_search(MJD, TAI_TO_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), index, (int (*)(const void *, const void *)) geMJDtimeTAI_UTC); + /* Search for index such that TAI_UTC_lookup_table[*index] <= MJD(TAI) < TAI_UTC_lookup_table[*index+1] */ + bhunt_search(MJD, TAI_UTC_lookup_table, number_of_entries_in_tai_utc_table, sizeof(TAI_UTC), index, (int (*)(const void *, const void *)) geMJDtime_TAI); if(debug == 2) fprintf(stderr, "*index = %d\n", *index); if(*index == -1) { /* MJD is less than first table entry. */ /* Debug: check that condition is met */ - if(debug && geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[*index+1])) { + if(debug && geMJDtime_TAI(MJD, &TAI_UTC_lookup_table[*index+1])) { fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index); exit(EXIT_FAILURE); } @@ -954,11 +967,11 @@ /* Calculate this offset strictly from offset1. The slope term doesn't enter because offset2 is the same as the UTC of the first epoch of the table. */ - return -TAI_TO_UTC_lookup_table[*index+1].offset1; + return -TAI_UTC_lookup_table[*index+1].offset1; } else if(*index == number_of_entries_in_tai_utc_table-1) { /* MJD is greater than or equal to last table entry. */ /* Debug: check that condition is met */ - if(debug && !geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[*index])) { + if(debug && !geMJDtime_TAI(MJD, &TAI_UTC_lookup_table[*index])) { fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index); exit(EXIT_FAILURE); } @@ -967,26 +980,26 @@ /* Use final offset for MJD values after last table entry. The slope term doesn't enter because modern values of the slope are zero.*/ - return -TAI_TO_UTC_lookup_table[*index].offset1; + return -TAI_UTC_lookup_table[*index].offset1; } else if(*index >= 0 && *index < number_of_entries_in_tai_utc_table) { /* table[*index] <= MJD < table[*index+1]. */ /* Debug: check that condition is met */ - if(debug && !(geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[*index]) && !geMJDtimeTAI_UTC(MJD, &TAI_TO_UTC_lookup_table[*index+1]))) { + if(debug && !(geMJDtime_TAI(MJD, &TAI_UTC_lookup_table[*index]) && !geMJDtime_TAI(MJD, &TAI_UTC_lookup_table[*index+1]))) { fprintf(stderr, "MJD = {%d, %f}\n", MJD->base_day, MJD->time_sec); fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad condition for *index = %d\n", *index); exit(EXIT_FAILURE); } - leap = -(TAI_TO_UTC_lookup_table[*index].offset1 + ((MJD->base_day-TAI_TO_UTC_lookup_table[*index].offset2) + MJD->time_sec/SecInDay)*TAI_TO_UTC_lookup_table[*index].slope)/(1. + TAI_TO_UTC_lookup_table[*index].slope/SecInDay); + leap = -(TAI_UTC_lookup_table[*index].offset1 + ((MJD->base_day-TAI_UTC_lookup_table[*index].offset2) + MJD->time_sec/SecInDay)*TAI_UTC_lookup_table[*index].slope)/(1. + TAI_UTC_lookup_table[*index].slope/SecInDay); /* Convert MJD(TAI) to normalized MJD(UTC). */ MJD->time_sec += leap; normalize_MJD(MJD); /* If MJD(UT) is in the next interval of the corresponding - UTC_TO_TAI_lookup_table, then we are right in the middle of a + TAI_UTC_lookup_table, then we are right in the middle of a leap interval (recently a second but for earlier epochs it could be less) insertion. Note this logic even works when leap intervals are taken away from UTC (i.e., leap is positive) since in that case the UTC *index always corresponds to the TAI *index. */ - *inleap = geMJDtimeTAI_UTC(MJD, &UTC_TO_TAI_lookup_table[*index+1]); + *inleap = geMJDtime_UTC(MJD, &TAI_UTC_lookup_table[*index+1]); return leap; } else { fprintf(stderr, "libqsastime (leap_second_TAI) logic ERROR: bad *index = %d\n", *index); Modified: trunk/lib/qsastime/tai-utc-gen.c =================================================================== --- trunk/lib/qsastime/tai-utc-gen.c 2009-05-07 12:43:52 UTC (rev 9935) +++ trunk/lib/qsastime/tai-utc-gen.c 2009-05-07 18:16:53 UTC (rev 9936) @@ -148,11 +148,11 @@ fprintf(fw, "%s\n",header); + fprintf(fw, "typedef struct {\n\tint base_day;\n\tdouble time_sec_tai;\n\tdouble time_sec_utc;\n\tdouble size_prev_leap_sec;\n\tdouble offset1;\n\tint offset2;\n\tdouble slope;\n} TAI_UTC;\n\n"); + fprintf(fw, "const int number_of_entries_in_tai_utc_table=%d;\n\n",number_of_lines); - fprintf(fw, "typedef struct {\n\tint base_day;\n\tdouble time_sec;\n\tdouble size_prev_leap_sec;\n\tdouble offset1;\n\tint offset2;\n\tdouble slope;\n} TAI_UTC;\n\n"); - - fprintf(fw, "const TAI_UTC TAI_TO_UTC_lookup_table[%d] = {\n",number_of_lines); + fprintf(fw, "const TAI_UTC TAI_UTC_lookup_table[%d] = {\n",number_of_lines); for (i=0;i<number_of_lines;i++) { sec = offset1[i] + (double)(MJDstart[i]-offset2[i])*slope[i]; if(i==0) @@ -166,16 +166,10 @@ unambiguously. */ leap_sec[i] = sec - (offset1[i-1] + (double)(MJDstart[i]+sec/86400.-offset2[i-1])*slope[i-1])/(1. + slope[i-1]/86400.); if(fabs(leap_sec[i]) < 1.e-14) leap_sec[i] = 0.; - fprintf(fw,"{%d, %15.8f, %20.14f, %15.8f, %d, %15.8f},\n", MJDstart[i], sec, leap_sec[i], offset1[i], offset2[i], slope[i]); + fprintf(fw,"{%d, %15.8f, 0., %20.14f, %15.8f, %d, %15.8f},\n", MJDstart[i], sec, leap_sec[i], offset1[i], offset2[i], slope[i]); } fprintf(fw,"};\n"); - fprintf(fw, "const TAI_UTC UTC_TO_TAI_lookup_table[%d] = {\n",number_of_lines); - for (i=0;i<number_of_lines;i++) { - fprintf(fw,"{%d, %15.8f, %20.14f, %15.8f, %d, %15.8f},\n", MJDstart[i], 0., leap_sec[i], offset1[i], offset2[i], slope[i]); - } - fprintf(fw,"};\n"); - fclose(fw); free(MJDstart); free(offset1); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |