#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#define PI 3.141592653589793
#define EARTH_MAJOR_AXIS 6378140.
#define EARTH_MINOR_AXIS 6356755.
#define EARTH_AXIS_RATIO (EARTH_MINOR_AXIS / EARTH_MAJOR_AXIS)

/* The following code reads in the MPC file ObsCodes.txt,  which gives
longitudes and parallax constants;  and outputs the same file with new
columns for latitude and altitude.  Be warned that the parallax
constants are usually given to five places,  or a precision of about
64 meters.  Therefore,  the latitude and longitude are good to a
similar precision.  The longitude is given to .0001 degree precision,
and is therefore good to about 11 meters.       */

void main( int argc, char **argv)
{
   double lon, x, y, lat;
   char buff[190];
   int code, in_code_section = 0;
   FILE *ifile = fopen( "stations.txt", "rb");

   if( !ifile)
      {
      printf( "stations.txt not opened\n");
      exit( 0);
      }
   while( fgets( buff, sizeof( buff), ifile))
      {
      if( !memcmp( buff, "000  ", 5))
         in_code_section = 1;
      if( !memcmp( buff, "</pre>", 6))
         in_code_section = 0;
      if( in_code_section && strlen( buff) > 30 && buff[3] == ' ')
         {
         memmove( buff + 45, buff + 30, strlen( buff + 29));
         memset( buff + 30, ' ', 15);
         if( sscanf( buff + 1, "%d %lf%lf%lf", &code, &lon, &x, &y) == 4)
            if( x)
               {
               char tbuff[37], alt_buff[7];

               if( buff[19] != ' ' && buff[28] != ' ')
                  {        /* make sure there's precision for altitude */
                  double radius = sqrt( x * x + y * y);
                  double ang = atan( y / x);
                  double r_earth = cos( ang) * cos( ang) +
                           sin( ang) * sin( ang) * EARTH_AXIS_RATIO;

                  sprintf( alt_buff, "%5.0lf",
                               (radius - r_earth) * EARTH_MAJOR_AXIS);
                  }
               else        /* put blanks in for altitude */
                  strcpy( alt_buff, "     ");

               lat = atan( y / (x * EARTH_AXIS_RATIO * EARTH_AXIS_RATIO));
               if( lon > 180.)
                  lon -= 360.;
               memcpy( tbuff, buff, 30);
                  tbuff[30] = '\0';
               sprintf( buff + 30, "%8.4lf %s", lat * 180. / PI,
                           alt_buff);
               buff[44] = ' ';
               }
         if( in_code_section)
            printf( "%s", buff);
         }
      }
}