/trunk/lib/qsastime/deltaT-gen.c

# · C · 170 lines · 95 code · 29 blank · 46 comment · 20 complexity · b254d5740eba0699db9427e0c67d8245 MD5 · raw file

  1. // $Id: deltaT-gen.c 11971 2011-10-14 09:28:58Z andrewross $
  2. //
  3. // Copyright (C) 2009 Alan W. Irwin
  4. //
  5. // This file is part of PLplot.
  6. // PLplot is free software; you can redistribute it and/or modify
  7. // it under the terms of the GNU Library General Public License as published
  8. // by the Free Software Foundation; either version 2 of the License, or
  9. // (at your option) any later version.
  10. //
  11. // PLplot is distributed in the hope that it will be useful,
  12. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14. // GNU Library General Public License for more details.
  15. //
  16. // You should have received a copy of the GNU Library General Public License
  17. // along with PLplot; if not, write to the Free Software
  18. // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19. //
  20. //
  21. // Program for generating spline representation (xspline, yspline,
  22. // and y2spline arrays) header from deltaT.dat.
  23. //
  24. // The program assumes that argv[1] will be the input file, and
  25. // argv[2] the output file. This works cross-platform without
  26. // worrying about shell redirects of stdin and stdout that are
  27. // not accessible on Windows, apparently.
  28. #include <stdio.h>
  29. #include <stdlib.h>
  30. #include <string.h>
  31. #include <math.h>
  32. #include "dspline.h"
  33. //--------------------------------------------------------------------------
  34. // Function-like macro definitions
  35. //--------------------------------------------------------------------------
  36. #define MemError1( a ) do { fprintf( stderr, "MEMORY ERROR %d\n" a "\n", __LINE__ ); exit( __LINE__ ); } while ( 0 )
  37. const char header[] = "" \
  38. "/*\n" \
  39. " This file is part of PLplot.\n" \
  40. " \n" \
  41. " PLplot is free software; you can redistribute it and/or modify\n" \
  42. " it under the terms of the GNU Library General Public License as published\n" \
  43. " by the Free Software Foundation; either version 2 of the License, or\n" \
  44. " (at your option) any later version.\n" \
  45. " \n" \
  46. " PLplot is distributed in the hope that it will be useful,\n" \
  47. " but WITHOUT ANY WARRANTY; without even the implied warranty of\n" \
  48. " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n" \
  49. " GNU Library General Public License for more details.\n" \
  50. " \n" \
  51. " You should have received a copy of the GNU Library General Public License\n" \
  52. " along with PLplot; if not, write to the Free Software\n" \
  53. " Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA\n" \
  54. " \n" \
  55. " \n" \
  56. " This header file contains spline data (xspline, yspline, and y2spline)\n" \
  57. " for converting between UT1 and ephemeris time.\n" \
  58. " It is an automatically generated file, so please do\n" \
  59. " not edit it directly. Make any changes to deltaT.dat then use\n" \
  60. " deltaT-gen to recreate this header file.\n" \
  61. " \n" \
  62. "*/";
  63. int main( int argc, char *argv[] )
  64. {
  65. FILE *fr, *fw;
  66. char readbuffer[256];
  67. double *xspline = NULL;
  68. double *yspline = NULL;
  69. double *y2spline = NULL;
  70. int i = 0;
  71. int number_of_lines = 0;
  72. if ( ( argc < 2 ) || ( fr = fopen( argv[1], "r" ) ) == NULL )
  73. {
  74. fprintf( stderr, "Cannot open first file as readable\n" );
  75. exit( 1 );
  76. }
  77. if ( ( argc < 3 ) || ( fw = fopen( argv[2], "w" ) ) == NULL )
  78. {
  79. fprintf( stderr, "Cannot open second file as writable\n" );
  80. exit( 1 );
  81. }
  82. //
  83. // Work out how many lines we have all up
  84. //
  85. while ( ( fgets( readbuffer, 255, fr ) != NULL ) )
  86. {
  87. ++number_of_lines;
  88. }
  89. //
  90. // Allocate memory to the arrays which will hold the data
  91. //
  92. if ( ( xspline = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL )
  93. MemError1( "Allocating memory to the xspline table" );
  94. if ( ( yspline = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL )
  95. MemError1( "Allocating memory to the yspline table" );
  96. if ( ( y2spline = (double *) calloc( (size_t) number_of_lines, (size_t) sizeof ( double ) ) ) == NULL )
  97. MemError1( "Allocating memory to the y2spline table" );
  98. rewind( fr ); // Go back to the start of the file
  99. //
  100. // Read in line by line, and copy the numbers into our arrays
  101. //
  102. while ( ( fgets( readbuffer, 255, fr ) != NULL ) )
  103. {
  104. sscanf( readbuffer, "%lf %lf", (double *) &xspline[i], (double *) &yspline[i] );
  105. i++;
  106. }
  107. fclose( fr );
  108. // Calculate spline representation using second derivative condition
  109. // on end points that is consistent with overall average parabolic shape
  110. // of delta T curve (Morrison and Stephenson, 2004) with second
  111. // derivative = 6.4e-3 secs/year/year.
  112. dspline( xspline, yspline, number_of_lines, 2, 6.4e-3, 2, 6.4e-3, y2spline );
  113. //
  114. // Write the data out to file ready to be included in our source
  115. //
  116. fprintf( fw, "%s\n", header );
  117. fprintf( fw, "const int number_of_entries_in_spline_tables=%d;\n\n", number_of_lines );
  118. fprintf( fw, "const double xspline[%d] = {\n", number_of_lines );
  119. for ( i = 0; i < number_of_lines; i++ )
  120. {
  121. fprintf( fw, "%10.0f,\n", xspline[i] );
  122. }
  123. fprintf( fw, "};\n" );
  124. fprintf( fw, "const double yspline[%d] = {\n", number_of_lines );
  125. for ( i = 0; i < number_of_lines; i++ )
  126. {
  127. fprintf( fw, "%10.0f,\n", yspline[i] );
  128. }
  129. fprintf( fw, "};\n" );
  130. fprintf( fw, "const double y2spline[%d] = {\n", number_of_lines );
  131. for ( i = 0; i < number_of_lines; i++ )
  132. {
  133. fprintf( fw, "%25.15e,\n", y2spline[i] );
  134. }
  135. fprintf( fw, "};\n" );
  136. fclose( fw );
  137. free( xspline );
  138. free( yspline );
  139. free( y2spline );
  140. return ( 0 );
  141. }