PLplot
5.10.0
|
00001 // Copyright (C) 2009-2014 Alan W. Irwin 00002 // 00003 // This file is part of PLplot. 00004 // PLplot is free software; you can redistribute it and/or modify 00005 // it under the terms of the GNU Library General Public License as published 00006 // by the Free Software Foundation; either version 2 of the License, or 00007 // (at your option) any later version. 00008 // 00009 // PLplot is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 // GNU Library General Public License for more details. 00013 // 00014 // You should have received a copy of the GNU Library General Public License 00015 // along with PLplot; if not, write to the Free Software 00016 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00017 // 00018 // 00019 00020 // Program for generating data structure used for containing tai-utc 00021 // conversion information (linear transforms and leap seconds). 00022 // 00023 // The program assumes that argv[1] will be the input file, and 00024 // argv[2] the output file. This works cross-platform without 00025 // worrying about shell redirects of stdin and stdout that are 00026 // not accessible on Windows, apparently. 00027 00028 #include <stdio.h> 00029 #include <stdlib.h> 00030 #include <string.h> 00031 #include <math.h> 00032 00033 00034 //-------------------------------------------------------------------------- 00035 // Function-like macro definitions 00036 //-------------------------------------------------------------------------- 00037 00038 #define MemError1( a ) do { fprintf( stderr, "MEMORY ERROR %d\n" a "\n", __LINE__ ); exit( __LINE__ ); } while ( 0 ) 00039 00040 const char header[] = "" \ 00041 "/*\n" \ 00042 " This file is part of PLplot.\n" \ 00043 " \n" \ 00044 " PLplot is free software; you can redistribute it and/or modify\n" \ 00045 " it under the terms of the GNU Library General Public License as published\n" \ 00046 " by the Free Software Foundation; either version 2 of the License, or\n" \ 00047 " (at your option) any later version.\n" \ 00048 " \n" \ 00049 " PLplot is distributed in the hope that it will be useful,\n" \ 00050 " but WITHOUT ANY WARRANTY; without even the implied warranty of\n" \ 00051 " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n" \ 00052 " GNU Library General Public License for more details.\n" \ 00053 " \n" \ 00054 " You should have received a copy of the GNU Library General Public License\n" \ 00055 " along with PLplot; if not, write to the Free Software\n" \ 00056 " Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA\n" \ 00057 " \n" \ 00058 " \n" \ 00059 " This header file contains the table containing the linear transforms \n" \ 00060 " for converting between TAI and UTC.\n" \ 00061 " It is an automatically generated file, so please do\n" \ 00062 " not edit it directly. Make any changes to tai-utc.dat then use\n" \ 00063 " tai-utc-gen to recreate this header file.\n" \ 00064 " \n" \ 00065 " tai-utc.dat contains four essential fields to represent the following\n" \ 00066 " formula for the linear transformation between TAI and UTC: \n" \ 00067 " TAI-UTC (seconds) = offset1 + (MJD-offset2)*slope\n" \ 00068 " There are four essential fields per line in tai-utc.dat to represent\n" \ 00069 " this formula. They are the Julian date (UTC) where the linear\n" \ 00070 " transformation implied by the line is first applied;\n" \ 00071 " offset1 (seconds); offset2 (days), and slope (secs/day).\n" \ 00072 " \n" \ 00073 "*/"; 00074 00075 int main( int argc, char *argv[] ) 00076 { 00077 FILE *fr, *fw; 00078 char readbuffer[256]; 00079 int *MJDstart = NULL; 00080 double *offset1 = NULL; 00081 int *offset2 = NULL; 00082 double *slope = NULL; 00083 double sec, *leap_sec = NULL; 00084 int jd; 00085 int i = 0; 00086 int number_of_lines = 0; 00087 00088 if ( ( argc < 2 ) || ( fr = fopen( argv[1], "r" ) ) == NULL ) 00089 { 00090 fprintf( stderr, "Cannot open first file as readable\n" ); 00091 exit( 1 ); 00092 } 00093 00094 if ( ( argc < 3 ) || ( fw = fopen( argv[2], "w" ) ) == NULL ) 00095 { 00096 fprintf( stderr, "Cannot open second file as writable\n" ); 00097 exit( 1 ); 00098 } 00099 00100 // 00101 // Work out how many lines we have all up 00102 // 00103 00104 while ( ( fgets( readbuffer, 255, fr ) != NULL ) ) 00105 { 00106 ++number_of_lines; 00107 } 00108 00109 // 00110 // Allocate memory to the arrays which will hold the data 00111 // 00112 00113 if ( ( MJDstart = (int *) calloc( (size_t) number_of_lines, (size_t) sizeof ( int ) ) ) == NULL ) 00114 MemError1( "Allocating memory to the MJDstart table" ); 00115 00116 if ( ( offset1 = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL ) 00117 MemError1( "Allocating memory to the offset1 table" ); 00118 00119 if ( ( offset2 = (int *) calloc( (size_t) number_of_lines, (size_t) sizeof ( int ) ) ) == NULL ) 00120 MemError1( "Allocating memory to the offset2 table" ); 00121 00122 if ( ( slope = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL ) 00123 MemError1( "Allocating memory to the slope table" ); 00124 00125 if ( ( leap_sec = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL ) 00126 MemError1( "Allocating memory to the leap_sec table" ); 00127 00128 rewind( fr ); // Go back to the start of the file 00129 00130 // 00131 // Read in line by line, and copy the numbers into our arrays 00132 // 00133 00134 while ( ( fgets( readbuffer, 255, fr ) != NULL ) ) 00135 { 00136 sscanf( readbuffer, "%*s %*s %*s %*s %d.5 %*s %lf %*s %*s %*s %*s %d.) X %lf S", (int *) &jd, (double *) &offset1[i], (int *) &offset2[i], (double *) &slope[i] ); 00137 // Should be exact since all jd's in the file are integer+0.5 00138 MJDstart[i] = jd - 2400000; 00139 i++; 00140 } 00141 00142 fclose( fr ); 00143 00144 // 00145 // Write the data out to file ready to be included in our source 00146 // 00147 00148 00149 fprintf( fw, "%s\n", header ); 00150 00151 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" ); 00152 00153 fprintf( fw, "const int number_of_entries_in_tai_utc_table=%d;\n\n", number_of_lines ); 00154 00155 fprintf( fw, "const TAI_UTC TAI_UTC_lookup_table[%d] = {\n", number_of_lines ); 00156 for ( i = 0; i < number_of_lines; i++ ) 00157 { 00158 sec = offset1[i] + (double) ( MJDstart[i] - offset2[i] ) * slope[i]; 00159 if ( i == 0 ) 00160 leap_sec[i] = 0.; 00161 else 00162 // sec is TAI-UTC in seconds calculated from UTC transformation 00163 // (equation 1 in README.tai-utc). This calculation must be correct 00164 // for start of epoch range. However, near end of epoch range where 00165 // ambiguities in UTC occur, must use equivalent TAI transformation 00166 // (equation 2 from same source) to calculate the UTC discontinuity 00167 // unambiguously. 00168 leap_sec[i] = sec - ( offset1[i - 1] + (double) ( MJDstart[i] + sec / 86400. - offset2[i - 1] ) * slope[i - 1] ) / ( 1. + slope[i - 1] / 86400. ); 00169 if ( fabs( leap_sec[i] ) < 1.e-14 ) 00170 leap_sec[i] = 0.; 00171 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] ); 00172 } 00173 fprintf( fw, "};\n" ); 00174 00175 fclose( fw ); 00176 free( MJDstart ); 00177 free( offset1 ); 00178 free( offset2 ); 00179 free( slope ); 00180 free( leap_sec ); 00181 00182 return ( 0 ); 00183 } 00184