home *** CD-ROM | disk | FTP | other *** search
/ Collection of Education / collectionofeducationcarat1997.iso / SCIENCE / AA_51.ZIP / KFILES.C < prev    next >
C/C++ Source or Header  |  1993-02-13  |  7KB  |  323 lines
  1.  
  2. /* Disc file input routines to read initialization parameters
  3.  * or file containing orbital elements.
  4.  */
  5. #include "kep.h"
  6.  
  7. extern char *intfmt, *strfmt;/* see dms.c */
  8.  
  9. static char starnam[80] = {'s','t','a','r','.','c','a','t','\0'};
  10. static char orbnam[80] = {'o','r','b','i','t','.','c','a','t','\0'};
  11. static int linenum = 1;
  12.  
  13. /* Read initialization file aa.ini
  14.  * and adjust topocentric coordinates of observer.
  15.  *
  16.  * The following items will be read in automatically from the disc file
  17.  * named aa.ini, if one is provided.  The file contains one ASCII
  18.  * string number per line so is easily edited.
  19.  *
  20.  * Terrestrial geocentric latitude and longitude of observer
  21.  * in degrees East of Greenwich, North of equator
  22.  * (kfiles.c converts from geodetic latitude input.)
  23.  */
  24. double tlong = -71.13;    /* Cambridge, Massachusetts */
  25. double tlat = 42.38; /* geocentric */
  26. double glat =42.17; /* geodetic */
  27. /* Parameters for calculation of azimuth and elevation
  28.  */
  29. double attemp = 20.0;    /* atmospheric temperature, degrees Centigrade */
  30. double atpress = 1013.0; /* atmospheric pressure, millibars */
  31. /* Distance from observer to center of earth, in earth radii
  32.  */
  33. double trho = 0.9985;
  34. static double flat = 298.257222;
  35. static double aearth = 6378137.;
  36. static double height = 0.0;
  37. extern double tlong, tlat, glat, trho, attemp, atpress, dtgiven;
  38.  
  39. int kinit()
  40. {
  41. double a, b, fl, co, si, u;
  42. FILE *f, *fopen();
  43. char s[84];
  44.  
  45. printf( "\n\tSteve Moshier's Ephemeris Program v5.1\n\n" );
  46. #if BandS
  47. printf( "Planetary perturbations are from Bretagnon and Simon.\n" );
  48. #else
  49. printf( "Planetary perturbations are from Meeus.\n" );
  50. #endif
  51. printf( "Extended Chapront Lunar theory.\n" );
  52.  
  53. f = fopen( "aa.ini", "r" );
  54. if( f )
  55.     {
  56.     fgets( s, 80, f );
  57.     sscanf( s, "%lf", &tlong );
  58.     fgets( s, 80, f );
  59.     sscanf( s, "%lf", &glat );
  60.     fgets( s, 80, f );
  61.     sscanf( s, "%lf", &height );
  62.     u = glat * DTR;
  63.  
  64. /* Reduction from geodetic latitude to geocentric latitude
  65.  * AA page K5
  66.  */
  67.     co = cos(u);
  68.     si = sin(u);
  69.     fl = 1.0 - 1.0/flat;
  70.     fl = fl*fl;
  71.     si = si*si;
  72.     u = 1.0/sqrt( co*co + fl*si );
  73.     a = aearth*u + height;
  74.     b = aearth*fl*u  +  height;
  75.     trho = sqrt( a*a*co*co + b*b*si );
  76.     tlat = RTD * acos( a*co/trho );
  77.     if( glat < 0.0 )
  78.         tlat = -tlat;
  79.     trho /= aearth;
  80. /* Reduction from geodetic latitude to geocentric latitude
  81.  * AA page K5
  82.  */
  83. /*
  84.     tlat = glat
  85.         - 0.19242861 * sin(2.0*u)
  86.         + 0.00032314 * sin(4.0*u)
  87.         - 0.00000072 * sin(6.0*u);
  88.  
  89.     trho =    0.998327073
  90.         + 0.001676438 * cos(2.0*u)
  91.         - 0.000003519 * cos(4.0*u)
  92.         + 0.000000008 * cos(6.0*u);
  93.     trho += height/6378160.;
  94. */
  95.     printf( "Terrestrial east longitude %.4lf deg\n", tlong );
  96.     printf( "geocentric latitude %.4lf deg\n", tlat );
  97.     printf( "Earth radius %.5lf\n", trho );
  98.  
  99.     fgets( s, 80, f );
  100.     sscanf( s, "%lf", &attemp );
  101.     printf( "temperature %.1lf C\n", attemp );
  102.     fgets( s, 80, f );
  103.     sscanf( s, "%lf", &atpress );
  104.     printf( "pressure %.0lf mb\n", atpress );
  105.     fgets( s, 80, f );
  106.     sscanf( s, "%d", &jdflag );
  107.     switch( jdflag )
  108.         {
  109.         case 0: printf("TDT and UT assumed equal.\n");
  110.             break;
  111.         case 1: printf("Input time is TDT.\n" );
  112.             break;
  113.         case 2: printf("Input time is UT.\n" );
  114.             break;
  115.         default: printf("Illegal jdflag\n" );
  116.         exit(0);
  117.         }
  118.     fgets( s, 80, f );
  119.     sscanf( s, "%lf", &dtgiven );
  120.     if( dtgiven != 0.0 )
  121.         printf( "Using deltaT = %.2lfs.\n", dtgiven );
  122.     fclose(f);
  123.     }
  124. return(0);
  125. }
  126.  
  127.  
  128.  
  129. /* Program to read in a file containing orbital parameters
  130.  */
  131. extern struct orbit earth;
  132.  
  133. int getorbit(el)
  134. struct orbit *el;
  135. {
  136. FILE *f, *fincat();
  137. char s1[100], s2[100], *u, *v;
  138. int i;
  139.  
  140.  
  141. getnum( "Name of orbit catalogue file: ", orbnam, strfmt );
  142. f = fincat( orbnam, 2, s1, s2 );
  143. if( f == 0 )
  144.     return(-1); /* failure flag */
  145.  
  146.  
  147. printf( "%s\n", s1 );
  148. printf( "%s\n", s2 );
  149.  
  150. /* Read in ASCII floating point numbers
  151.  */
  152. sscanf( s1, "%lf %lf %lf %lf %lf %lf",
  153.     &el->epoch, &el->i, &el->W, &el->w, &el->a, &el->dm );
  154.  
  155. sscanf( s2, "%lf %lf %lf %lf %lf %15s", &el->ecc, &el->M,
  156.     &el->equinox, &el->mag, &el->sdiam, &el->obname[0] );
  157.  
  158. el->obname[15] = '\0';
  159.  
  160. /* Clear out the rest of the data structure
  161.  */
  162. el->oelmnt = 0;
  163. el->celmnt = 0;
  164. el->L = 0.0;
  165. el->r = 0.0;
  166. el->plat = 0.0;
  167. if( strcmp( &el->obname[0], "Earth" )  )
  168.     {
  169.     return(0);
  170.     }
  171. else
  172.     {
  173.     u = (char *)&earth;
  174.     v = (char *)el;
  175.     for( i=0; i<sizeof(struct orbit); i++ )
  176.         *u++ = *v++;
  177.     printf( "Read in earth orbit\n" );
  178.     return(1);
  179.     }
  180. }
  181.  
  182.  
  183.  
  184. int getstar(el)
  185. struct star *el;
  186. {
  187. int sign;
  188. char s[100];
  189. double rh, rm, rs, dd, dm, ds, x, z;
  190. FILE *f, *fincat();
  191. char *p;
  192. int i;
  193.  
  194. getnum( "Name of star catalogue file: ", starnam, strfmt );
  195. f = fincat( starnam, 1, s );
  196. if( f == 0 )
  197.     return(-1); /* failure flag */
  198.  
  199.  
  200. printf( "%s\n", s );
  201. /* Read in the ASCII string data and name of the object
  202.  */
  203. sscanf( s, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %s",
  204.     &el->epoch, &rh, &rm, &rs, &dd, &dm, &ds,
  205.     &el->mura, &el->mudec, &el->v, &el->px, &x, &el->obname[0] );
  206.  
  207. x = el->epoch;
  208. if( x == 2000.0 )
  209.     x = J2000;
  210. else if( x == 1950.0 )
  211.     x = B1950;
  212. else if( x == 1900.0 )
  213.     x = J1900;
  214. else
  215.     x = J2000  +  365.25 * (x - 2000.0);
  216. el->epoch = x;
  217.  
  218. /* read the right ascension */
  219. el->ra = 2.0 * PI * (3600.0*rh + 60.0*rm + rs)/86400.0;
  220.  
  221. /* read the declination */
  222. sign = 1;
  223. if( (dd < 0.0) || (dm < 0.0) || (ds < 0.0) )
  224.     sign = -1;
  225. z = (3600.0*fabs(dd) + 60.0*fabs(dm) + fabs(ds))/RTS;
  226. if( dd == 0.0 )
  227.     {
  228. /* Scan the text for possible minus sign in front of declination 0 */
  229.     p = s;
  230. /* skip over 4 fields */
  231.     for( i=0; i<4; i++ )
  232.         {
  233.         while( *p++ == ' ' )
  234.             ;
  235.         while( *p++ != ' ' )
  236.             ;
  237.         }
  238.     while( *p++ == ' ' )
  239.         ;
  240.     --p;
  241.     if( *p == '-' )
  242.         sign = -1;
  243.     }
  244. if( sign < 0 )
  245.     z = -z;
  246. el->dec = z;
  247.  
  248. #if DEBUG
  249. printf( "%.2lf\n", el->epoch );
  250. printf( "%.0lf %.0lf %.3lf\n", rh, rm, rs );
  251. printf( "%.8lf\n", el->ra );
  252. printf( "%.0lf %.0lf %.3lf\n", dd, dm, ds );
  253. printf( "%.8lf\n", el->dec );
  254. printf( "d %.3lf mua %.3lf mud %.3lf v %.3lf\n",
  255.     el->px, el->mura, el->mudec, el->v );
  256. #endif
  257.  
  258. el->mura *= 15.0/RTS;    /* s/century -> "/century -> rad/century */
  259. el->mudec /= RTS;
  260. z = el->px;
  261. if( z < 1.0 )
  262.     {
  263.     if( z <= 0.0 )
  264.         el->px = 0.0;
  265.     else
  266.         el->px = STR * z;  /* assume px in arc seconds */
  267.     }
  268. else
  269.     {
  270.     el->px = 1.0/(RTS * z);    /* parsecs -> radians */
  271.     }
  272. return(0);
  273. }
  274.  
  275.  
  276.  
  277.  
  278. /* Open catalogue and find line number
  279.  */
  280. FILE *fincat( name, n, str1, str2 )
  281. char *name;
  282. int n;   /* number of lines per catalogue entry */
  283. char *str1, *str2;
  284. {
  285. int i;
  286. FILE *f, *fopen();
  287.  
  288. f = fopen( name, "r" );
  289. if( f == 0 )
  290.     {
  291.     printf( "Can't find file %s\n", name );
  292.     return(0); /* failure flag */
  293.     }
  294.  
  295. getnum( "Line number", &linenum, intfmt );
  296. if( linenum <= 0 )
  297.     goto failure;
  298.  
  299. for( i=0; i<linenum; i++ )
  300.     {
  301.  
  302.     fgets( str1, 98, f );
  303.     if( *str1 == '-' )
  304.         goto endf;
  305.     if( n > 1 )
  306.         {
  307.         fgets( str2, 98, f );
  308.         if( *str2 == '-' )
  309.             goto endf;
  310.         }
  311.     }
  312. fclose(f);
  313. return( f );
  314.  
  315.  
  316. endf:
  317. printf( "End of file reached.\n" );
  318. failure:
  319. fclose(f);
  320. return(0);
  321. }
  322.  
  323.