home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / unix / volume12 / starcharts / part03 < prev    next >
Encoding:
Internet Message Format  |  1987-11-29  |  41.8 KB
  1. Subject:  v12i073:  StarChart program and Yale star data, Part03/07
  2. Newsgroups: comp.sources.unix
  3. Sender: sources
  4. Approved: rs@uunet.UU.NET
  5.  
  6. Submitted-by: awpaeth@watcgl.waterloo.edu (Alan W. Paeth)
  7. Posting-number: Volume 12, Issue 73
  8. Archive-name: starchart/part03
  9.  
  10. # This is a shell archive.  Remove anything before this line,
  11. # then unpack it by saving it in a file and typing "sh file".
  12. #
  13. # Wrapped by watcgl!awpaeth on Mon Oct  5 19:01:46 EDT 1987
  14. # Contents:  planet.c planet.1 yaleformat yale.omit
  15.  
  16. echo x - planet.c
  17. sed 's/^@//' > "planet.c" <<'@//E*O*F planet.c//'
  18. /*  Planets
  19.  *           A program to determine the position of
  20.  *         the Sun., Mer., Ven., Mars, Jup, Sat & Ura
  21.  *
  22.  *  reference Jean Meesus's book, " Astronomical Formulae for
  23.  *  Calculators
  24.  *            program by
  25.  *
  26.  *            F.T. Mendenhall
  27.  *            ihnp4!inuxe!fred
  28.  *
  29.  *        Copyright 1985 by F.T. Mendenhall
  30.  *        all rights reserved
  31.  *
  32.  * (This copyright prevents you from selling the program - the
  33.  *  author grants anyone the right to copy and install the program
  34.  *  on any machine it will run on)
  35.  !
  36.  ! Modified by Alan Paeth, awpaeth@watcgl, January 1987 for
  37.  !   (1) non-interactive mode (uses current GMT)
  38.  !   (2) output in simplified Yale format for use with star charter software
  39.  !   (3) lines 365-380 rewritten as a simplier C expression (awpaeth, Apr '87)
  40.  !
  41.  ! Modified by Petri Launiainen, pl@intrin.FI, February 1987 for
  42.  !   easier LOGFILE definition/modification in Makefile
  43.  !
  44.  ! Modified by Alan Paeth, September 1987, for command line flags.
  45.  !    The program continues to provide reduced Yale output in the old style,
  46.  !     but new versions of "starchart" accept either. The old format allows
  47.  !     representation of magnitudes below -.99 (useful for planets), but no
  48.  !     spectral or stellar identification, which is not useful for planets.
  49.  !
  50.  ! See "Time Notes" section below if you plan to rewrite the user interface.
  51.  */
  52.  
  53. #include <stdio.h>
  54. #include <math.h>
  55.  
  56. #include <sys/time.h>    /* for getting current GMT */
  57.  
  58. #define CURYEAR    1987    /* default year -- needs to be maintained */
  59.  
  60. #ifndef PLANETFILE
  61. #define PLANETFILE "./planet.star"
  62. #endif
  63.  
  64. /* crude magnitudes of planets (x100) for output filtering by brightness */
  65.  
  66. #define MAGSOL 0
  67. #define MAGMER 150
  68. #define MAGVEN 0
  69. #define MAGMAR 100
  70. #define MAGJUP 0
  71. #define MAGSAT 100
  72. #define MAGURA 590
  73. #define MAGNEP 800
  74.  
  75. double pie, rad;
  76. double htod(), atof(), kepler(), truean();
  77. double longi(), lati(), poly(), aint(), range();
  78. int cls(), trans(), speak();
  79. FILE *logfile;
  80. char *progname;
  81.  
  82. /*
  83.  * Times Notes.
  84.  *
  85.  * This software would ideally use a generic UNIX system call to which would
  86.  * provide day and date info as integers, plus current time west of GMT. These
  87.  * values would become the default settings for each of the -m -d -y -t & -z.
  88.  * Unfortunately, time on across many UNIX systems does not (to my knowledge)
  89.  * exist as a uniform subroutine call.
  90.  *
  91.  * Instead, we presume that the very low level "gettimeofday" (GMT in seconds,
  92.  * from 1 Jan 1970) exists on all UNIX installations, and use this subroutine
  93.  * to use the present time of day when no command line parameters are present,
  94.  * as defaults for the current time have not been filled in.
  95.  *
  96.  * When any flag-style command line parameters are present, hardcoded defaults
  97.  * are substitude for -m, -d, -y -t switches, which may then be overridden.
  98.  * The -z flag is filled in by a companion return value from "gettimeofday".
  99.  * This approach makes for poor (on account of hard coding) defaults.
  100.  *
  101.  * A worthwhile change would be to simplify the command line control structure
  102.  * and default mechanism by using a higher level system call, or (safer),
  103.  * write the routine to convert Unix GMT into day, month, year, so that it
  104.  * lives within this program. This would allow for more reasonable defaults,
  105.  * and still merely one UNIX system call to get GMT time. Perhaps sources
  106.  * from a UNIX system may be studied to do this correctly.
  107.  *
  108.  * Expertise in "time" on a specific UNIX system is *NOT* what is required --
  109.  * what *IS* required is knowledge and access to enough independent UNIX
  110.  * systems to insure (to a high probability) that the code remains portable.
  111.  * In other words, only low-level system calls known to work across a large
  112.  * range of systems should be employed.
  113.  *
  114.  */
  115.  
  116. double currenttime()
  117.     {
  118. #define GMT1970 2440587.5
  119. #define SECSPERDAY (60.0*60.0*24.0)
  120.     struct timeval tv;
  121.     struct timezone tz;
  122.     gettimeofday(&tv, &tz);
  123.     return(GMT1970 + tv.tv_sec/SECSPERDAY);
  124.     }
  125.  
  126. double commandtime(argc, argv)
  127.     char **argv;
  128.     {
  129.     struct timeval tv;
  130.     struct timezone tz;
  131.     int mm,dd,yy,aa1,bb1;
  132.     double jd,year,month,day,time,timez;
  133.     int j;
  134.     static char *usage =
  135. "\nusage:\tplanet {current time used}\nor\tplanet juliandate\nor\tplanet [-t hh.mm -z zone -y year -m mon -d day ]";
  136. /*
  137.  * check command line style
  138.  */
  139.     if (argc == 1) return (currenttime());
  140.  /*
  141.   * no args - use current time
  142.   */
  143.     if  ((argv[1][0] != '-'))
  144.     {
  145. /*
  146.  * one numeric arg - use it as the julian date
  147.  */
  148.     if (argc == 2) return(atof(argv[1]));
  149.     die("no switches - %s", usage);
  150.     }
  151.     else
  152.     {
  153. /*
  154.  * flags present, set up defaults
  155.  */
  156.     gettimeofday(&tv, &tz);
  157.     time = 0.0;                /* local time (0.0-23.999) */
  158.     yy = CURYEAR;            /* local year */
  159.     mm = 0;                /* local month */
  160.     dd = 0;                /* local day of month */
  161.     timez = tz.tz_minuteswest/60.0;    /* local time hrs west of Greenwich */
  162. /*
  163.  * override defaults
  164.  */
  165.     for (j=1; j<argc; j++)
  166.         {
  167.         if (argv[j][0] != '-')
  168.         die("unknown argument - %s", usage /*argv[j]*/ );
  169.         switch (argv[j][1])
  170.         {
  171.         case 't': time = htod(argv[++j]); break;
  172.         case 'z': timez = htod(argv[++j]); break;
  173.         case 'y': yy = atoi(argv[++j]); break;
  174.         case 'm': mm = atoi(argv[++j]); break;
  175.         case 'd': dd = atoi(argv[++j]); break;
  176.         default:  die("unknown switch - %s", usage /* argv[j] */ );
  177.         }
  178.         if (j == argc) die("trailing command line flag - %s", argv[j-1]);
  179.         }
  180.     }
  181.     if (!mm || !dd) die("both month and day flags must be present");
  182.     day = ((float)dd + (time + timez)/24.0);
  183.     if (mm > 2)
  184.     {
  185.     year = yy;
  186.     month = mm;
  187.     }
  188.     else
  189.     {
  190.     year = yy - 1;
  191.     month = mm + 12;
  192.     }
  193.     aa1 = (int) (year/100);
  194.     bb1 = 2 - aa1 + (int)(aa1/4);
  195.     jd = aint(365.25*year) + aint(30.6001*(month + 1));
  196.     jd = jd  + day + 1720994.5;
  197.     if((year + month/100) > 1582.10) jd = jd + bb1;
  198.     return(jd);
  199.     }
  200.  
  201. main(argc, argv)
  202.     int argc;
  203.     char **argv;
  204.     {
  205.     double jd, T0;
  206.     double l,a,e,i,w1,w2,M,M0,M1,M2,M4,M5,M6,M7,M8;
  207.     double a0,a1,a2,a3,RA,DEC;
  208.     double ECC,nu,r,u,ll,b,lonpert,radpert;
  209.     double esun,Lsun,Cen,Stheta,Snu,Sr;
  210.     double N,D,epli,thapp,omeg;
  211.     double nu2,P,Q,S,V,W,ze,l1pert,epert,w1pert,apert;
  212.     double psi,H,G,eta,th;
  213.  
  214. #define WRITEMODE "w"
  215. #define OPENFAIL 0
  216.     progname = argv[0];
  217.     if ((logfile = fopen(PLANETFILE, WRITEMODE)) == OPENFAIL)
  218.         die(stderr, "fail on opening log file %s\n", PLANETFILE);
  219.     cls();
  220.     pie = 3.1415926535898;
  221.     rad = pie/180.0;
  222.  
  223. /*  calculate time T0 from 1900 January 0.5 ET */
  224.     jd = commandtime(argc, argv);
  225.     printf("Julian date: %f\n", jd);
  226.     T0 = (jd - 2415020.0)/36525.0;
  227.  
  228. /* Calculate orbital elements for Mercury */
  229.     a0=178.179078;
  230.     a1=149474.07078;
  231.     a2=0.0003011;
  232.     a3=0.0;
  233.     l = poly(a0,a1,a2,a3,T0);
  234.     l = range(l);
  235.     a = 0.3870986;
  236.     a0 = 0.20561421;
  237.     a1 = 0.00002046;
  238.     a2 = -0.000000030;
  239.     e = poly(a0,a1,a2,a3,T0);
  240.     a0 = 7.002881;
  241.     a1 = 0.0018608;
  242.     a2 = -0.0000183;
  243.     i = poly(a0,a1,a2,a3,T0);
  244.     a0 = 28.753753;
  245.     a1 = 0.3702806;
  246.     a2 = 0.0001208;
  247.     w1 = poly(a0,a1,a2,a3,T0);
  248.     w1 = range(w1);
  249.     a0 = 47.145944;
  250.     a1 = 1.1852083;
  251.     a2= 0.0001739;
  252.     w2 = poly(a0,a1,a2,a3,T0);
  253.     w2 = range(w2);
  254.     M1 = 102.27938 + 149472.51529*T0 + 0.000007*T0*T0;
  255.     M1 = range(M1);
  256.     
  257. /* solving Kepler find the eccentric anomaly ECC    */
  258.  
  259.     ECC =kepler(e,M1);
  260.     nu = truean(e,ECC);
  261.     r = a*(1.0 - e * cos(ECC*pie/180.0));
  262.     u = l + nu - M1 - w2;
  263.     ll = longi(w2,i,u);
  264.     b = lati(u,i);
  265.     
  266. /* Now make corrections due to perturbations */
  267.     M2 = 212.60322 + 58517.80387*T0 +0.001286*T0*T0;
  268.     M2 = range(M2);
  269.     M4 = 319.51913 + 19139.85475*T0 + 0.000181*T0*T0;
  270.     M4 = range(M4);
  271.     M5 = 225.32833 + 3034.69202*T0 - 0.000722*T0*T0;
  272.     M5 = range(M5);
  273.     M6 = 175.46622 +1221.55147*T0 - 0.000502*T0*T0;
  274.     M6 = range(M6);
  275.     lonpert = 0.00204*cos((5*M2-2*M1+12.220)*rad)
  276.          +0.00103*cos((2*M2-M1-160.692)*rad)
  277.          +0.00091*cos((2*M5-M1-37.003)*rad)
  278.          +0.00078*cos((5*M2-3*M1+10.137)*rad);
  279.         
  280.     radpert = 0.000007525*cos((2*M5-M1+53.013)*rad)
  281.          +0.000006802*cos((5*M2-3*M1-259.918)*rad)
  282.          +0.000005457*cos((2*M2-2*M1-71.188)*rad)
  283.          +0.000003569*cos((5*M2-M1-77.75)*rad);
  284.         
  285.         r = r + radpert;
  286.         ll = ll + lonpert;
  287.         
  288. /* Now find the Longi and radius vector of the sun */
  289.         
  290.         M= 358.47583 + 35999.04975*T0 - 0.000150*T0*T0
  291.            -0.0000033*T0*T0*T0;
  292.         M = range(M);
  293.         esun = 0.01675104 - 0.0000418*T0 -0.000000126*T0*T0;
  294.         Lsun = 279.69668 + 36000.76892*T0 + 0.0003025*T0*T0;
  295.         Lsun = range(Lsun);
  296.         Cen= (1.919460-0.004789*T0-0.000014*T0*T0)*sin(M*rad)
  297.              +(0.020094 - 0.000100*T0)*sin(2*M*rad)
  298.              +0.000293*sin(3*M*rad);
  299.         
  300.         Stheta= Lsun + Cen;
  301.         Stheta = range(Stheta);
  302.         Snu = M + Cen;
  303.         Sr = 1.0000002*(1.0 - esun*esun)/(1.0 + esun*cos(Snu*rad));
  304.  
  305.         omeg= 259.18 - 1934.142*T0;
  306.         thapp = Stheta - 0.00569 - 0.00479* sin(omeg*rad);
  307.         epli = 23.452294 - 0.0130125*T0
  308.                -0.00000164*T0*T0 + 0.000000503*T0*T0*T0
  309.             +0.00256*cos(omeg*rad);
  310.         printf("\nOBJECT     R.A.           DEC    \t DISTANCE\n");
  311.         N = cos(epli*rad)*sin(thapp*rad);
  312.         D = cos(thapp*rad);
  313.         RA = atan2(N,D)/rad;
  314.         DEC = asin(sin(epli*rad)*sin(thapp*rad))/rad;
  315.         speak(RA,DEC,Sr, MAGSOL, "PS", "Sol");
  316. /* tansformation of coordinates on Mercury and output */
  317.         trans(r,b,ll,Stheta,Sr,epli, MAGMER, "PM", "Mercury");
  318.         
  319. /* Now start on Venus */
  320.     a0 = 342.767053;
  321.     a1 = 58519.21191;
  322.     a2 = 0.0003097;
  323.     a3 = 0.0;
  324.     l = poly(a0,a1,a2,a3,T0);
  325.     l = range(l);
  326.     a = 0.7233316;
  327.     a0 = 0.00682069;
  328.     a1 = -0.00004774;
  329.     a2 = 0.000000091;
  330.     e = poly(a0,a1,a2,a3,T0);
  331.     a0 = 3.393631;
  332.     a1 = 0.0010058;
  333.     a2 = -0.0000010;
  334.     i = poly(a0,a1,a2,a3,T0);
  335.     a0 = 54.384186;
  336.     a1 = 0.5081861;
  337.     a2 = -0.0013864;
  338.     w1 = poly(a0,a1,a2,a3,T0);
  339.     w1 = range(w1);
  340.     a0 = 75.779647;
  341.     a1 = 0.8998500;
  342.     a2 = 0.0004100;
  343.     w2 = poly(a0,a1,a2,a3,T0);
  344.     w2 = range(w2);
  345. /* Venus has a long period pert that needs to be added before Kelper */
  346.     lonpert = 0.00077 *sin((237.24 + 150.27*T0)*rad);
  347.     l = l + lonpert;
  348.     M0 = M2 + lonpert;
  349.     ECC = kepler(e,M0);
  350.     nu = truean(e,ECC);
  351.     r = a*(1.0 - e * cos(ECC*rad));
  352.     u = l + nu - M0- w2;
  353.     ll = longi(w2,i,u);
  354.     b = lati(u,i);
  355.     
  356. /* now Venus perturbations */
  357.     
  358.     lonpert = 0.00313*cos((2*M-2*M2 -148.225)*rad)
  359.          +0.00198*cos((3*M-3*M2 +2.565)*rad)
  360.          +0.00136*cos((M-M2-119.107)*rad)
  361.          +0.00096*cos((3*M-2*M2-135.912)*rad)
  362.          +0.00082*cos((M5-M2-208.087)*rad);
  363.     ll = ll + lonpert;
  364.     
  365.     radpert = 0.000022501 * cos((2*M-2*M2-58.208)*rad)
  366.          +0.000019045 * cos((3*M-3*M2+92.577)*rad)
  367.          +0.000006887 * cos((M5-M2-118.090)*rad)
  368.          +0.000005172 * cos((M-M2-29.110)*rad)
  369.          +0.000003620 * cos((5*M-4*M2-104.208)*rad)
  370.          +0.000003283 * cos((4*M-4*M2+63.513)*rad)
  371.          +0.000003074 * cos((2*M5-2*M2-55.167)*rad);
  372.     r = r + radpert;
  373.     trans(r,b,ll,Stheta,Sr,epli, MAGVEN, "PV", "Venus");
  374.     
  375. /* Now  start the planet Mars */
  376.     a0 = 293.737334;
  377.     a1 = 19141.69551;
  378.     a2 = 0.0003107;
  379.     a3 = 0.0;
  380.     l = poly(a0,a1,a2,a3,T0);
  381.     l = range(l);
  382.     a = 1.5236883;
  383.     a0 = 0.09331290;
  384.     a1 = 0.000092064;
  385.     a2 = -0.000000077;
  386.     e = poly(a0,a1,a2,a3,T0);
  387.     a0 = 1.850333;
  388.     a1 = -0.0006750;
  389.     a2 = 0.0000126;
  390.     i = poly(a0,a1,a2,a3,T0);
  391.     a0 = 285.431761;
  392.     a1 = 1.0697667;
  393.     a2 = 0.0001313;
  394.     a3 = 0.00000414;
  395.     w1 = poly(a0,a1,a2,a3,T0);
  396.     w1 = range(w1);
  397.     a0 = 48.786442;
  398.     a1 = 0.7709917;
  399.     a2 = -0.0000014;
  400.     a3 = -0.00000533;
  401.     w2 = poly(a0,a1,a2,a3,T0);
  402.     w2 = range(w2);
  403.  
  404. /* Mars has a long period perturbation */
  405.     lonpert = -0.01133*sin((3*M5-8*M4 +4*M)*rad)
  406.           -0.00933*cos((3*M5-8*M4 +4*M)*rad);
  407.     l = l + lonpert;
  408.     M0 = M4 + lonpert;
  409.     ECC = kepler(e,M0);
  410.     nu = truean(e,ECC);
  411.     r = a*(1.0 - e * cos(ECC*rad));
  412.     u = l + nu - M0- w2;
  413.     ll = longi(w2,i,u);
  414.     b = lati(u,i);
  415.     
  416. /* Now Mars Perturbations */    
  417.     lonpert = 0.00705*cos((M5-M4-48.958)*rad)
  418.          +0.00607*cos((2*M5-M4-188.350)*rad)
  419.          +0.00445*cos((2*M5-2*M4-191.897)*rad)
  420.          +0.00388*cos((M-2*M4+20.495)*rad)
  421.          +0.00238*cos((M-M4+35.097)*rad)
  422.          +0.00204*cos((2*M-3*M4+158.638)*rad)
  423.          +0.00177*cos((3*M4-M2-57.602)*rad)
  424.          +0.00136*cos((2*M-4*M4+154.093)*rad)
  425.          +0.00104*cos((M5+17.618)*rad);
  426.     ll = ll + lonpert;
  427.     
  428.     radpert = 0.000053227*cos((M5-M4+41.1306)*rad)
  429.          +0.000050989*cos((2*M5-2*M4-101.9847)*rad)
  430.          +0.000038278*cos((2*M5-M4-98.3292)*rad)
  431.          +0.000015996*cos((M-M4-55.555)*rad)
  432.          +0.000014764*cos((2*M-3*M4+68.622)*rad)
  433.          +0.000008966*cos((M5-2*M4+43.615)*rad);
  434. /*
  435.  * broken into two pieces for wimpy C compilers
  436.  */
  437.     radpert+= 0.000007914*cos((3*M5-2*M4-139.737)*rad)
  438.          +0.000007004*cos((2*M5-3*M4-102.888)*rad)
  439.          +0.000006620*cos((M-2*M4+113.202)*rad)
  440.          +0.000004930*cos((3*M5-3*M4-76.243)*rad)
  441.          +0.000004693*cos((3*M-5*M4+190.603)*rad)
  442.          +0.000004571*cos((2*M-4*M4+244.702)*rad)
  443.          +0.000004409*cos((3*M5-M4-115.828)*rad);
  444.     r = r + radpert;
  445.     trans(r,b,ll,Stheta,Sr,epli, MAGMAR, "PM", "Mars");
  446.     
  447. /* Now start Jupiter */
  448.     a0 = 238.049257;
  449.     a1 = 3036.301986;
  450.     a2 = 0.0003347;
  451.     a3 = -0.00000165;
  452.     l = poly(a0,a1,a2,a3,T0);
  453.     l = range(l);
  454.     a = 5.202561;
  455.     a0 = 0.04833475;
  456.     a1 = 0.000164180;
  457.     a2 = -0.0000004676;
  458.     a3 = -0.0000000017;
  459.     e = poly(a0,a1,a2,a3,T0);
  460.     a0 = 1.308736;
  461.     a1 = -0.0056961;
  462.     a2 = 0.0000039;
  463.     a3 = 0.0;
  464.     i = poly(a0,a1,a2,a3,T0);
  465.     a0 = 273.277558;
  466.     a1 = 0.5994317;
  467.     a2 = 0.00070405;
  468.     a3 = 0.00000508;
  469.     w1 = poly(a0,a1,a2,a3,T0);
  470.     w1 = range(w1);
  471.     a0 = 99.443414;
  472.     a1 = 1.0105300;
  473.     a2 = 0.00035222;
  474.     a3 = -0.00000851;
  475.     w2 = poly(a0,a1,a2,a3,T0);
  476.     w2 = range(w2);
  477. /*
  478.     Now start perturbation calclations */
  479.     
  480.     nu2 = T0/5.0 +0.1;
  481.     P = 237.47555 +3034.9061*T0;
  482.     Q = 265.91650 + 1222.1139*T0;
  483.     S = 243.51721 + 428.4677*T0;
  484.     V = 5.0*Q -2.0*P;
  485.     W = 2.0*P - 6.0*Q +3.0*S;
  486.     ze = Q - P;
  487.     psi = S - Q;
  488.     
  489.     P = range(P)*rad;
  490.     Q = range(Q)*rad;
  491.     S = range(S)*rad;
  492.     V = range(V)*rad;
  493.     W = range(W)*rad;
  494.     ze = range(ze)*rad;
  495.     psi = range(psi)*rad;
  496.     
  497. l1pert = (0.331364 - 0.010281*nu2 - 0.004692*nu2*nu2)*sin(V)
  498.     +(0.003228 - 0.064436*nu2 + 0.002075*nu2*nu2)*cos(V)
  499.     -(0.003083 + 0.000275*nu2 - 0.000489*nu2*nu2)*sin(2*V)
  500.     +0.002472*sin(W)
  501.     +0.013619*sin(ze)
  502.     +0.018472*sin(2*ze)
  503.     +0.006717*sin(3*ze)
  504.     +0.002775*sin(4*ze)
  505.     +(0.007275 - 0.001253*nu2)*sin(ze)*sin(Q)
  506.     +0.006417*sin(2*ze)*sin(Q)
  507.     +0.002439*sin(3*ze)*sin(Q);
  508.     
  509. l1pert = l1pert    -(0.033839 + 0.001125*nu2)*cos(ze)*sin(Q)
  510.     -0.003767*cos(2*ze)*sin(Q)
  511.     -(0.035681 + 0.001208*nu2)*sin(ze)*cos(Q)
  512.     -0.004261*sin(2*ze)*cos(Q)
  513.     +0.002178*cos(Q)
  514.     +(-0.006333 + 0.001161*nu2)*cos(ze)*cos(Q)
  515.     -0.006675*cos(2*ze)*cos(Q)
  516.     -0.002664*cos(3*ze)*cos(Q)
  517.     -0.002572*sin(ze)*sin(2*Q)
  518.     -0.003567*sin(2*ze)*sin(2*Q)
  519.     +0.002094*cos(ze)*cos(2*Q)
  520.     +0.003342*cos(2*ze)*cos(2*Q);
  521.     
  522. epert =  (.0003606 + .0000130*nu2 - .0000043*nu2*nu2)*sin(V)
  523.     +(.0001289 - .0000580*nu2)*cos(V)
  524.     -.0006764*sin(ze)*sin(Q)
  525.     -.0001110*sin(2*ze)*sin(Q)
  526.     -.0000224*sin(3*ze)*sin(Q)
  527.     -.0000204*sin(Q)
  528.     +(.0001284 + .0000116*nu2)*cos(ze)*sin(Q)
  529.     +.0000188*cos(2*ze)*sin(Q)
  530.     +(.0001460 + .0000130*nu2)*sin(ze)*cos(Q)
  531.     +.0000224*sin(2*ze)*cos(Q)
  532.     -.0000817*cos(Q);
  533.     
  534. epert = epert    +.0006074*cos(ze)*cos(Q)
  535.     +.0000992*cos(2*ze)*cos(Q)
  536.     +.0000508*cos(3*ze)*cos(Q)
  537.     +.0000230*cos(4*ze)*cos(Q)
  538.     +.0000108*cos(5*ze)*cos(Q)
  539.     -(.0000956 + .0000073*nu2)*sin(ze)*sin(2*Q)
  540.     +.0000448*sin(2*ze)*sin(2*Q)
  541.     +.0000137*sin(3*ze)*sin(2*Q)
  542.     +(-.0000997 + .0000108*nu2)*cos(ze)*sin(2*Q)
  543.     +.0000480*cos(2*ze)*sin(2*Q);
  544.  
  545. epert = epert    +.0000148*cos(3*ze)*sin(2*Q)
  546.     +(-.0000956 +.0000099*nu2)*sin(ze)*cos(2*Q)
  547.     +.0000490*sin(2*ze)*cos(2*Q)
  548.     +.0000158*sin(3*ze)*cos(2*Q)
  549.     +.0000179*cos(2*Q)
  550.     +(.0001024 + .0000075*nu2)*cos(ze)*cos(2*Q)
  551.     -.0000437*cos(2*ze)*cos(2*Q)
  552.     -.0000132*cos(3*ze)*cos(2*Q);
  553.     
  554. w1pert = (0.007192 - 0.003147*nu2)*sin(V)
  555.     +(-0.020428 - 0.000675*nu2 + 0.000197*nu2*nu2)*cos(V)
  556.     +(0.007269 + 0.000672*nu2)*sin(ze)*sin(Q)
  557.     -0.004344*sin(Q)
  558.     +0.034036*cos(ze)*sin(Q)
  559.     +0.005614*cos(2*ze)*sin(Q)
  560.     +0.002964*cos(3*ze)*sin(Q)
  561.     +0.037761*sin(ze)*cos(Q);
  562.     
  563. w1pert = w1pert
  564.     +0.006158*sin(2*ze)*cos(Q)
  565.     -0.006603*cos(ze)*cos(Q)
  566.     -0.005356*sin(ze)*sin(2*Q)
  567.     +0.002722*sin(2*ze)*sin(2*Q)
  568.     +0.004483*cos(ze)*sin(2*Q)
  569.     -0.002642*cos(2*ze)*sin(2*Q)
  570.     +0.004403*sin(ze)*cos(2*Q)
  571.     -0.002536*sin(2*ze)*cos(2*Q)
  572.     +0.005547*cos(ze)*cos(2*Q)
  573.     -0.002689*cos(2*ze)*cos(2*Q);
  574.     
  575. lonpert = l1pert -(w1pert/e);
  576. l = l + l1pert;
  577. M0 = M5 + lonpert;
  578. e = e + epert;
  579. w1 = w1 + w1pert;
  580.  
  581. apert =  -.000263*cos(V)
  582.     +.000205*cos(ze)
  583.     +.000693*cos(2*ze)
  584.     +.000312*cos(3*ze)
  585.     +.000147*cos(4*ze)
  586.     +.000299*sin(ze)*sin(Q)
  587.     +.000181*cos(2*ze)*sin(Q)
  588.     +.000204*sin(2*ze)*cos(Q)
  589.     +.000111*sin(3*ze)*cos(Q)
  590.     -.000337*cos(ze)*cos(Q)
  591.     -.000111*cos(2*ze)*cos(Q);
  592.     
  593. a = a + apert;
  594. ECC = kepler(e,M0);
  595. nu = truean(e,ECC);
  596. r = a*(1.0 - e * cos(ECC*rad));
  597. u = l + nu - M0 - w2;
  598. ll = longi(w2,i,u);
  599. b = lati(u,i);
  600.  
  601. trans(r,b,ll,Stheta,Sr,epli, MAGJUP, "PJ", "Jupiter");
  602.  
  603. /* Now start Saturn */
  604.  
  605.     a0 = 266.564377;
  606.     a1 = 1223.509884;
  607.     a2 = 0.0003245;
  608.     a3 = -0.0000058;
  609.     l = poly(a0,a1,a2,a3,T0);
  610.     l = range(l);
  611.     a = 9.554747;
  612.     a0 = 0.05589232;
  613.     a1 = -0.00034550;
  614.     a2 = -0.000000728;
  615.     a3 = 0.00000000074;
  616.     e = poly(a0,a1,a2,a3,T0);
  617.     a0 = 2.492519;
  618.     a1 = -0.0039189;
  619.     a2 = -0.00001549;
  620.     a3 = 0.00000004;
  621.     i = poly(a0,a1,a2,a3,T0);
  622.     a0 = 338.307800;
  623.     a1 = 1.0852207;
  624.     a2 = 0.00097854;
  625.     a3 = 0.00000992;
  626.     w1 = poly(a0,a1,a2,a3,T0);
  627.     w1 = range(w1);
  628.     a0 = 112.790414;
  629.     a1 = 0.8731951;
  630.     a2 = -0.00015218;
  631.     a3 = -0.00000531;
  632.     w2 = poly(a0,a1,a2,a3,T0);
  633.     w2 = range(w2);
  634.     
  635. /* Now Saturn's perturbations */
  636.     
  637. l1pert = (-0.814181 + 0.018150*nu2 + 0.016714*nu2*nu2)*sin(V)
  638.      +(-0.010497 + 0.160906*nu2 - 0.004100*nu2*nu2)*cos(V)
  639.      +0.007581*sin(2*V)
  640.      -0.007986*sin(W)
  641.      -0.148811*sin(ze)
  642.      -0.040786*sin(2*ze)
  643.      -0.015208*sin(3*ze)
  644.      -0.006339*sin(4*ze)
  645.      -0.006244*sin(Q);
  646. l1pert = l1pert
  647.     +(0.008931 + 0.002728*nu2)*sin(ze)*sin(Q)
  648.     -0.016500*sin(2*ze)*sin(Q)
  649.     -0.005775*sin(3*ze)*sin(Q)
  650.     +(0.081344 + 0.003206*nu2)*cos(ze)*sin(Q)
  651.     +0.015019*cos(2*ze)*sin(Q)
  652.     +(0.085581 + 0.002494*nu2)*sin(ze)*cos(Q)
  653.     +(0.025328 - 0.003117*nu2)*cos(ze)*cos(Q);
  654. l1pert = l1pert
  655.     +0.014394*cos(2*ze)*cos(Q)
  656.     +0.006319*cos(3*ze)*cos(Q)
  657.     +0.006369*sin(ze)*sin(2*Q)
  658.     +0.009156*sin(2*ze)*sin(2*Q)
  659.     +0.007525*sin(3*psi)*sin(2*Q)
  660.     -0.005236*cos(ze)*cos(2*Q)
  661.     -0.007736*cos(2*ze)*cos(2*Q)
  662.     -0.007528*cos(3*psi)*cos(2*Q);
  663.     
  664. epert = (-.0007927 + .0002548*nu2 +.0000091*nu2*nu2)*sin(V)
  665.     +(.0013381 + .0001226*nu2 -.0000253*nu2*nu2)*cos(V)
  666.     +(.0000248 - .0000121*nu2)*sin(2*V)
  667.     -(.0000305 + .0000091*nu2)*cos(2*V)
  668.     +.0000412*sin(2*ze)
  669.     +.0012415*sin(Q)
  670.     +(.0000390 -.0000617*nu2)*sin(ze)*sin(Q)
  671.     +(.0000165 - .0000204*nu2)*sin(2*ze)*sin(Q)
  672.     +.0026599*cos(ze)*sin(Q)
  673.     -.0004687*cos(2*ze)*sin(Q);
  674. epert = epert
  675.     -.0001870*cos(3*ze)*sin(Q)
  676.     -.0000821*cos(4*ze)*sin(Q)
  677.     -.0000377*cos(5*ze)*sin(Q)
  678.     +.0000497*cos(2*psi)*sin(Q)
  679.     +(.0000163 - .0000611*nu2)*cos(Q)
  680.     -.0012696*sin(ze)*cos(Q)
  681.     -.0004200*sin(2*ze)*cos(Q)
  682.     -.0001503*sin(3*ze)*cos(Q)
  683.     -.0000619*sin(4*ze)*cos(Q)
  684.     -.0000268*sin(5*ze)*cos(Q);
  685. epert = epert
  686.     -(.0000282 + .0001306*nu2)*cos(ze)*cos(Q)
  687.     +(-.0000086 + .0000230*nu2)*cos(2*ze)*cos(Q)
  688.     +.0000461*sin(2*psi)*cos(Q)
  689.     -.0000350*sin(2*Q)
  690.     +(.0002211 - .0000286*nu2)*sin(ze)*sin(2*Q)
  691.     -.0002208*sin(2*ze)*sin(2*Q)
  692.     -.0000568*sin(3*ze)*sin(2*Q)
  693.     -.0000346*sin(4*ze)*sin(2*Q)
  694.     -(.0002780 + .0000222*nu2)*cos(ze)*sin(2*Q)
  695.     +(.0002022 + .0000263*nu2)*cos(2*ze)*sin(2*Q);
  696. epert = epert
  697.     +.0000248*cos(3*ze)*sin(2*Q)
  698.     +.0000242*sin(3*psi)*sin(2*Q)
  699.     +.0000467*cos(3*psi)*sin(2*Q)
  700.     -.0000490*cos(2*Q)
  701.     -(.0002842 + .0000279*nu2)*sin(ze)*cos(2*Q)
  702.     +(.0000128 + .0000226*nu2)*sin(2*ze)*cos(2*Q)
  703.     +.0000224*sin(3*ze)*cos(2*Q)
  704.     +(-.0001594 + .0000282*nu2)*cos(ze)*cos(2*Q)
  705.     +(.0002162 - .0000207*nu2)*cos(2*ze)*cos(2*Q)
  706.     +.0000561*cos(3*ze)*cos(2*Q);
  707. epert = epert
  708.     +.0000343*cos(4*ze)*cos(2*Q)
  709.     +.0000469*sin(3*psi)*cos(2*Q)
  710.     -.0000242*cos(3*psi)*cos(2*Q)
  711.     -.0000205*sin(ze)*sin(3*Q)
  712.     +.0000262*sin(3*ze)*sin(3*Q)
  713.     +.0000208*cos(ze)*cos(3*Q)
  714.     -.0000271*cos(3*ze)*cos(3*Q)
  715.     -.0000382*cos(3*ze)*sin(4*Q)
  716.     -.0000376*sin(3*ze)*cos(4*Q);
  717.  
  718. w1pert = (0.077108 + 0.007186*nu2 - 0.001533*nu2*nu2)*sin(V)
  719.     +(0.045803 - 0.014766*nu2 - 0.000536*nu2*nu2)*cos(V)
  720.     -0.007075*sin(ze)
  721.     -0.075825*sin(ze)*sin(Q)
  722.     -0.024839*sin(2*ze)*sin(Q)
  723.     -0.008631*sin(3*ze)*sin(Q)
  724.     -0.072586*cos(Q)
  725.     -0.150383*cos(ze)*cos(Q)
  726.     +0.026897*cos(2*ze)*cos(Q)
  727.     +0.010053*cos(3*ze)*cos(Q);
  728. w1pert = w1pert
  729.     -(0.013597 +0.001719*nu2)*sin(ze)*sin(2*Q)
  730.     +(-0.007742 + 0.001517*nu2)*cos(ze)*sin(2*Q)
  731.     +(0.013586 - 0.001375*nu2)*cos(2*ze)*sin(2*Q)
  732.     +(-0.013667 + 0.001239*nu2)*sin(ze)*cos(2*Q)
  733.     +0.011981*sin(2*ze)*cos(2*Q)
  734.     +(0.014861 + 0.001136*nu2)*cos(ze)*cos(2*Q)
  735.     -(0.013064 + 0.001628*nu2)*cos(2*ze)*cos(2*Q);
  736.     
  737. lonpert = l1pert -(w1pert/e);
  738. l = l + l1pert;
  739. M0 = M6 + lonpert;
  740. e = e + epert;
  741. w1 = w1 + w1pert;
  742.  
  743. apert = .000572*sin(V) -.001590*sin(2*ze)*cos(Q)
  744.     +.002933*cos(V) -.000647*sin(3*ze)*cos(Q)
  745.     +.033629*cos(ze) -.000344*sin(4*ze)*cos(Q)
  746.     -.003081*cos(2*ze) +.002885*cos(ze)*cos(Q)
  747.     -.001423*cos(3*ze) +(.002172 + .000102*nu2)*cos(2*ze)*cos(Q)
  748.     -.000671*cos(4*ze) +.000296*cos(3*ze)*cos(Q)
  749.     -.000320*cos(5*ze) -.000267*sin(2*ze)*sin(2*Q);
  750. apert = apert
  751.     +.001098*sin(Q) -.000778*cos(ze)*sin(2*Q)
  752.     -.002812*sin(ze)*sin(Q) +.000495*cos(2*ze)*sin(2*Q)
  753.     +.000688*sin(2*ze)*sin(Q) +.000250*cos(3*ze)*sin(2*Q);
  754. apert = apert
  755.     -.000393*sin(3*ze)*sin(Q)
  756.     -.000228*sin(4*ze)*sin(Q)
  757.     +.002138*cos(ze)*sin(Q)
  758.     -.000999*cos(2*ze)*sin(Q)
  759.     -.000642*cos(3*ze)*sin(Q)
  760.     -.000325*cos(4*ze)*sin(Q)
  761.     -.000890*cos(Q)
  762.     +.002206*sin(ze)*cos(Q);
  763. apert = apert
  764.     -.000856*sin(ze)*cos(2*Q)
  765.     +.000441*sin(2*ze)*cos(2*Q)
  766.     +.000296*cos(2*ze)*cos(2*Q)
  767.     +.000211*cos(3*ze)*cos(2*Q)
  768.     -.000427*sin(ze)*sin(3*Q)
  769.     +.000398*sin(3*ze)*sin(3*Q)
  770.     +.000344*cos(ze)*cos(3*Q)
  771.     -.000427*cos(3*ze)*cos(3*Q);
  772.     
  773. a = a + apert;
  774. ECC = kepler(e,M0);
  775. nu = truean(e,ECC);
  776. r = a*(1.0 - e * cos(ECC*rad));
  777. u = l + nu - M0 - w2;
  778. ll = longi(w2,i,u);
  779. b = lati(u,i);
  780.  
  781. b = b
  782.     +0.000747*cos(ze)*sin(Q)
  783.     +0.001069*cos(ze)*cos(Q)
  784.     +0.002108*sin(2*ze)*sin(2*Q)
  785.     +0.001261*cos(2*ze)*sin(2*Q)
  786.     +0.001236*sin(2*ze)*cos(2*Q)
  787.     -0.002075*cos(2*ze)*cos(2*Q);
  788.  
  789. trans(r,b,ll,Stheta,Sr,epli, MAGSAT, "PS", "Saturn");
  790.  
  791. /* Now Start on Uranus */
  792.     a0 = 244.197470;
  793.     a1 = 429.863546;
  794.     a2 = 0.0003160;
  795.     a3 = -0.00000060;
  796.     l = poly(a0,a1,a2,a3,T0);
  797.     l = range(l);
  798.     a = 19.21814;
  799.     a0 = 0.0463444;
  800.     a1 = -0.00002658;
  801.     a2 = 0.000000077;
  802.     a3 = 0.0;
  803.     e = poly(a0,a1,a2,a3,T0);
  804.     a0 = 0.772464;
  805.     a1 = 0.0006253;
  806.     a2 = 0.0000395;
  807.     i = poly(a0,a1,a2,a3,T0);
  808.     a0 = 98.071581;
  809.     a1 = 0.9857650;
  810.     a2 = -0.0010745;
  811.     a3 = -0.00000061;
  812.     w1 = poly(a0,a1,a2,a3,T0);
  813.     w1 = range(w1);
  814.     a0 = 73.477111;
  815.     a1 = 0.4986678;
  816.     a2 = 0.0013117;
  817.     a3 = 0.0;
  818.     w2 = poly(a0,a1,a2,a3,T0);
  819.     w2 = range(w2);
  820.     M7 = 72.64878 + 428.37911*T0 + 0.000079*T0*T0;
  821.     M7 = range(M7);
  822. /* now perturbation corrections for Uranus */
  823.     G = 83.76922 + 218.4901*T0;
  824.     S = S/rad;
  825.     P = P/rad;
  826.     Q = Q/rad;
  827.     H = 2.0*G - S;
  828.     ze = S - P;
  829.     eta = S - Q;
  830.     th = G - S;
  831.     S = S*rad;
  832.     G = range(G)*rad;
  833.     P = P*rad;
  834.     Q = Q*rad;
  835.     H = range(H)*rad;
  836.     ze = range(ze)*rad;
  837.     eta = range(eta)*rad;
  838.     th = range(th)*rad;
  839.     
  840. l1pert = (0.864319 - 0.001583*nu2)*sin(H)
  841.     +(0.082222 - 0.006833*nu2)*cos(H)
  842.     +0.036017*sin(2*H)
  843.     -0.003019*cos(2*H)
  844.     +0.008122*sin(W);
  845.     
  846. epert = (-.0003349 + .0000163*nu2)*sin(H)
  847.     +.0020981*cos(H)
  848.     +.0001311*cos(H);
  849.     
  850. w1pert = 0.120303*sin(H)
  851.     +(0.019472 - 0.000947*nu2)*cos(H)
  852.     +0.006197*sin(2*H);
  853.     
  854. lonpert = l1pert -(w1pert/e);
  855. l = l + l1pert;
  856. M0 = M7 + lonpert;
  857. e = e + epert;
  858. w1 = w1 + w1pert;
  859.  
  860. a = a - 0.003825*cos(H);
  861. ECC = kepler(e,M0);
  862. nu = truean(e,ECC);
  863. r = a*(1.0 - e * cos(ECC*rad));
  864. u = l + nu - M0 - w2;
  865. ll = longi(w2,i,u);
  866. b = lati(u,i);
  867.  
  868. ll = ll
  869.     +(0.010122 - 0.000988*nu2)*sin(S+eta)
  870.     +(-0.038581 + 0.002031*nu2 - 0.001910*nu2*nu2)*cos(S+eta)
  871.     +(0.034964 - 0.001038*nu2 + 0.000868*nu2*nu2)*cos(2*S+eta)
  872.     +0.005594*sin(S +3*th);
  873. ll = ll
  874.     -0.014808*sin(ze)
  875.     -0.005794*sin(eta)
  876.     +0.002347*cos(eta)
  877.     +0.009872*sin(th)
  878.     +0.008803*sin(2*th)
  879.     -0.004308*sin(3*th);
  880.     
  881. b = b
  882.     +(0.000458*sin(eta) - 0.000642*cos(eta) - 0.000517*cos(4*th))
  883.     *sin(S)
  884.     -(0.000347*sin(eta) + 0.000853*cos(eta) + 0.000517*sin(4*eta))
  885.     *cos(S)
  886.     +0.000403*(cos(2*th)*sin(2*S) + sin(2*th)*cos(2*S));
  887.     
  888. r = r
  889.     -.025948
  890.     +.004985*cos(ze)
  891.     -.001230*cos(S)
  892.     +.003354*cos(eta)
  893.     +(.005795*cos(S) - .001165*sin(S) + .001388*cos(2*S))*sin(eta)
  894.     +(.001351*cos(S) + .005702*sin(S) + .001388*sin(2*S))*cos(eta)
  895.     +.000904*cos(2*th)
  896.     +.000894*(cos(th) - cos(3*th));
  897. trans(r,b,ll,Stheta,Sr,epli, MAGURA, "PU", "Uranus");
  898.  
  899. /* Now start Neptune */
  900.     a0 = 84.457994;
  901.     a1 = 219.885914;
  902.     a2 = 0.0003205;
  903.     a3 = -0.00000060;
  904.     l = poly(a0,a1,a2,a3,T0);
  905.     l = range(l);
  906.     a = 30.10957;
  907.     a0 = 0.00899704;
  908.     a1 = 0.000006330;
  909.     a2 = -0.000000002;
  910.     a3 = 0.0;
  911.     e = poly(a0,a1,a2,a3,T0);
  912.     a0 = 1.779242;
  913.     a1 = -0.0095436;
  914.     a2 = -0.0000091;
  915.     i = poly(a0,a1,a2,a3,T0);
  916.     a0 = 276.045975;
  917.     a1 = 0.3256394;
  918.     a2 = 0.00014095;
  919.     a3 = 0.000004113;
  920.     w1 = poly(a0,a1,a2,a3,T0);
  921.     w1 = range(w1);
  922.     a0 = 130.681389;
  923.     a1 = 1.0989350;
  924.     a2 = 0.00024987;
  925.     a3 = -0.000004718;
  926.     w2 = poly(a0,a1,a2,a3,T0);
  927.     w2 = range(w2);
  928.     M8 = 37.73063 + 218.46134*T0 -0.000070*T0*T0;
  929.     
  930. /* now perturbation corrections for neptune */
  931.     G = G/rad;
  932.     P = P/rad;
  933.     Q = Q/rad;
  934.     S = S/rad;
  935.     ze = G - P;
  936.     eta = G - Q;
  937.     th = G - S;
  938.     G = G*rad;
  939.     P = P*rad;
  940.     Q = Q*rad;
  941.     S = S*rad;
  942.     ze = range(ze)*rad;
  943.     eta = range(eta)*rad;
  944.     th = range(th)*rad;
  945.     
  946. l1pert = (-0.589833 + 0.001089*nu2)*sin(H)
  947.     +(-0.056094 + 0.004658*nu2)*cos(H)
  948.     -0.024286*sin(2*H);
  949.     
  950. epert = .0004389*sin(H)
  951.     +.0004262*cos(H)
  952.     +.0001129*sin(2*H)
  953.     +.0001089*cos(2*H);
  954.     
  955. w1pert = 0.024039*sin(H)
  956.     -0.025303*cos(H)
  957.     +0.006206*sin(2*H)
  958.     -0.005992*cos(2*H);
  959.     
  960.     
  961. lonpert = l1pert -(w1pert/e);
  962. l = l + l1pert;
  963. M0 = M8 + lonpert;
  964. e = e + epert;
  965. w1 = w1 + w1pert;
  966.  
  967. a = a - 0.000817*sin(H)
  968.     +0.008189*cos(H)
  969.     +0.000781*cos(2*H);
  970.     
  971. ECC = kepler(e,M0);
  972. nu = truean(e,ECC);
  973. r = a*(1.0 - e * cos(ECC*rad));
  974. u = l + nu - M0 - w2;
  975. ll = longi(w2,i,u);
  976. b = lati(u,i);
  977.  
  978. ll = ll
  979.     -0.009556*sin(ze)
  980.     -0.005178*sin(eta)
  981.     +0.002572*sin(2*th)
  982.     -0.002972*cos(2*th)*sin(G)
  983.     -0.002833*sin(2*th)*cos(G);
  984.     
  985. b = b
  986.     +0.000336*cos(2*th)*sin(G)
  987.     +0.000364*sin(2*th)*cos(G);
  988.     
  989. r = r
  990.     -.040596
  991.     +.004992*cos(ze)
  992.     +.002744*cos(eta)
  993.     +.002044*cos(th)
  994.     +.001051*cos(2*th);
  995. trans(r,b,ll,Stheta,Sr,epli, MAGNEP, "PN", "Neptune");
  996.  
  997. putchar('\n');
  998. close(logfile);
  999. exit(0);
  1000. } /* end of program main */
  1001.  
  1002.  
  1003. int cls()
  1004.     {
  1005. #define CLEARCNT 1
  1006.     int i;
  1007.     for (i=0; i < CLEARCNT; i++) putchar('\n');
  1008.     }
  1009.  
  1010. double poly(a0,a1,a2,a3,t)
  1011.     double a0,a1,a2,a3,t;
  1012.     {
  1013.     return(a0+a1*t+a2*t*t+a3*t*t*t);
  1014.     }
  1015. double aint(z)
  1016.     double z;
  1017. {
  1018.     int trunk;
  1019.     
  1020.     trunk = (int)z;
  1021.     z = (double) trunk;
  1022.     return(z);
  1023. }
  1024. double kepler(e,M)
  1025.     double e,M;
  1026. {
  1027.     double corr,e0,E0,E1;
  1028.     e0 = e /rad;
  1029.     corr = 1;
  1030.     E0=M;
  1031.     while (corr > 0.000001)
  1032.     {
  1033.         corr = (M + e0* sin(E0*rad) -E0)/(1-e*cos(E0*rad));
  1034.         E1 = E0 + corr;
  1035.         if (corr < 0)
  1036.         {
  1037.             corr = -1.0* corr ;
  1038.         }
  1039.         E0 = E1;
  1040.     }
  1041.     
  1042.     return(E1);
  1043. }
  1044.  
  1045. double truean(e,E)
  1046.     double e,E;
  1047. {
  1048.     double nu;
  1049.  
  1050.     nu = 2.0 * atan(sqrt((1+e)/(1-e))*tan(E*rad/2));
  1051.     nu = nu/rad;
  1052.     if (nu < 0.0)
  1053.     {
  1054.         nu = nu + 360.0;
  1055.     }
  1056.     if (fabs(nu-E) > 90.0)
  1057.     {
  1058.         if (nu > 180.0)
  1059.         {
  1060.             nu = nu - 180.0;
  1061.         }
  1062.         else
  1063.         {
  1064.             nu = nu + 180.0;
  1065.         }
  1066.     }
  1067. return(nu);
  1068. }
  1069.  
  1070. double longi(w2,i,u)
  1071.     double w2,i,u;
  1072. {
  1073.     double x,y,l;
  1074.     
  1075.     y=cos(i*rad)*sin(u*rad);
  1076.     x=cos(u*rad);
  1077.     l = atan2(y,x);
  1078.     l = l/rad;
  1079.     if (l < 0.0)
  1080.     {
  1081.         l = l + 360.0;
  1082.     }
  1083.     return(l+ w2);
  1084. }
  1085.  
  1086. double lati(u,i)
  1087.     double u,i;
  1088. {
  1089.     double b;
  1090.     b = asin(sin(u*rad)*sin(i*rad))/rad;
  1091.     return(b);
  1092. }
  1093.  
  1094. int speak(ra,dec,dis, mag, sym, str)
  1095.     double ra,dec,dis;
  1096.     char *sym, *str;
  1097. {
  1098.     int h,m,s,deg;
  1099.     
  1100.     if ( ra < 0)
  1101.     {
  1102.         ra = ra + 360.0;
  1103.     }
  1104.     ra = ra /15.0; /* convert from degs to hours */
  1105.     h = (int) aint(ra);
  1106.     m = (int)((ra-h)*60);
  1107.     s = (int) (((ra-h)*60-m)*60);
  1108.     printf("%s\t%2dh %2dm %2ds ",str, h,m,s);
  1109.     fprintf(logfile, "%02d%02d%02d", h, m, s);
  1110.     deg = (int)aint(dec);
  1111.     m   =(int)((dec - deg)*60);
  1112.     s   = (int) (((dec - deg)*60 - m)*60);
  1113.     if ( m < 0)
  1114.     {
  1115.         m = m * -1;
  1116.     }
  1117.     if (s < 0)
  1118.     {
  1119.         s = s * -1;
  1120.     }
  1121.     printf("   %+3ddeg %2dm %2ds ",deg,m,s);
  1122.     printf("\t  %.3f\n",dis);
  1123.     fprintf(logfile, "%+03d%02d+%03d%s%s\n", deg, m, mag, sym, str);
  1124.     return(0);
  1125. }
  1126.  
  1127. double range(val)
  1128.     double val;
  1129. {
  1130.     while (val < 0.0)
  1131.     {
  1132.         val = val + 360.0;
  1133.     }
  1134.     if (val > 360.0)
  1135.     {
  1136.         val = val - (aint(val/360.0)*360.0);
  1137.     }
  1138.     return(val);
  1139. }
  1140. int trans(r,b,ll,Stheta,Sr,epli, mag, sym, str)
  1141.     double r,b,ll,Stheta,Sr,epli;
  1142.     char *sym, *str;
  1143. {
  1144.     double N,D,lambda,delta,beta,RA,DEC;
  1145.     
  1146.         N = r * cos(b*rad) * sin((ll - Stheta)*rad);
  1147.         D = r * cos(b*rad) * cos((ll - Stheta)*rad) + Sr;
  1148.         lambda = atan2(N,D)/rad;
  1149.         if (lambda < 0.0)
  1150.         {
  1151.             lambda = lambda + 360.0;
  1152.         }
  1153.         lambda = range(lambda+Stheta);
  1154.         delta = sqrt(N*N + D*D + (r*sin(b*rad))*(r*sin(b*rad)));
  1155.         beta = asin(r*sin(b*rad)/delta)/rad;
  1156.         N = sin(lambda*rad)*cos(epli*rad)
  1157.             - tan(beta*rad)*sin(epli*rad);
  1158.         D = cos(lambda*rad);
  1159.         RA = atan2(N,D)/rad;
  1160.         DEC = asin(sin(beta*rad)*cos(epli*rad)
  1161.              +cos(beta*rad)*sin(epli*rad)*sin(lambda*rad))/rad;
  1162.         speak(RA,DEC,delta, mag, sym, str);
  1163.         return(0);
  1164. }
  1165.  
  1166. double htod(s)
  1167.     char *s;
  1168.     {
  1169. /*
  1170.  * htod(str) reads floating point strings of the form {+-}hh.{mm} thus allowing
  1171.  * for fractional values in base 60 (ie, degrees/minutes).
  1172.  */
  1173.     double x, sign;
  1174.     int full, frac;
  1175.     char *t;
  1176.     t = s-1;
  1177.     while(*++t)
  1178.     {
  1179.     if ( ( (*t<'0') || (*t>'9') ) && (*t!='.') && (*t!='+') && (*t!='-'))
  1180.         die("non-digit in dd.mm style numeric argument");
  1181.     }
  1182.     if (s == NULL) return 0.0;
  1183.     full = frac = 0;
  1184.     sign = 1.0;
  1185.     if (*s == '-')
  1186.     {
  1187.     sign = -1.0;
  1188.     s++;
  1189.     }
  1190.     else if (*s == '+') s++;
  1191.     while (*s && *s != '.') full = 10 * full + *s++ - '0';
  1192.     if (*s++ == '.')
  1193.     {
  1194.     if (*s) frac = 10 * (*s++ - '0');
  1195.     if (*s) frac += *s++ - '0';
  1196.     if (frac > 59) frac = 59;
  1197.     }
  1198.     x = (double) full + ((double) frac) / 60.0;
  1199.     return sign * x;
  1200.     }
  1201.  
  1202. die(a,b)
  1203.     char *a, *b;
  1204.     {
  1205.     fprintf(stderr,"%s: ", progname);
  1206.     fprintf(stderr,a,b);
  1207.     fprintf(stderr,"\n");
  1208.     exit(1);
  1209.     }
  1210. @//E*O*F planet.c//
  1211. chmod u=rwx,g=rwx,o=rwx planet.c
  1212.  
  1213. echo x - planet.1
  1214. sed 's/^@//' > "planet.1" <<'@//E*O*F planet.1//'
  1215. @.TH starchart LOCAL 6/2/87
  1216. @.ad b
  1217. @.SH NAME
  1218. planet, moonphase, epoch
  1219. @.br
  1220. \- print phase of the Moon and planetary information, precess data in .star file format.
  1221. @.SH SYNOPSIS
  1222. \fBplanet [juliandate]\fR
  1223. @.br
  1224. or
  1225. @.br
  1226. \fBplanet [-z timezone -m month -d day -y year -t hh.mm ]\fR
  1227. @.br
  1228. \fBmoonphase\fR
  1229. @.br
  1230. \fBepoch [ inputepoch ] [ outputepoch ]\fR
  1231. @.SH DESCRIPTION
  1232. These programs generate and maintain star chart data in the format used by
  1233. the \fBstarchart\fR software tools.
  1234. @.PP
  1235. Planetary information is generated by the program \fBplanet\fR, which calculates
  1236. locations a given Julian date, or for the current time (given no arguments).
  1237. In common use, an alternate syntax is provided to simplify the calculation of
  1238. the Julian date, in which integers for month, day, (mandatory), and year are
  1239. given. Time is specified in the form \fBhh.mm\fP or \fBhh\fP on a 24-hour clock
  1240. and defaults to 0.00 (midnight). The \fB-z\fP flag (same format as -t) gives
  1241. the number of hours the observer is west of GMT (see below).
  1242. @.PP
  1243. A summary is then printed on the standard output giving the right ascension,
  1244. declination and distance (in AU's) for the Sun and all planets except Pluto.
  1245. This data is additionally placed in the file \fBplanet.star\fR for subsequent
  1246. use by charting software.
  1247. @.PP
  1248. The program \fBmoonphase\fP gives the current lunar phase. At present it does
  1249. not provide coordinate data for subsequent plotting.
  1250. @.PP
  1251. The program \fBepoch\fP converts right-ascension and declination data within
  1252. \fByale.star\fR style datasets from epoch 1975.0 into the epoch 2000.0.
  1253. Either the first or both epoch dates may be overridden with command line
  1254. parameters (floating point values). The data is read and written from
  1255. the standard input and output. Because only the coordinate data is updated
  1256. within each line item, the program may be used on either versions of the
  1257. reduced Yale data.
  1258. @.SH BUGS
  1259. Daylight savings time may have to be reckoned for in \fBplanet\fP, as the
  1260. default \fB-z\fP value during periods of daylight time may or may not be
  1261. adjusted relative to GMT. The required code changes would make \fBplanet\fP
  1262. installation dependent.
  1263. @.SH AUTHORS
  1264. \fBplanet\fP - F.T. Mendenhall (ihnp4!inuxe!fred)
  1265. @.br
  1266. \fBmoonphase\fP - Keith E. Brandt
  1267. @.br
  1268. \fBepoch\fP - Alan Paeth, University of Waterloo (AWPaeth@watCGL)
  1269. @.br
  1270. \fByale.star\fP - Robert Tidd (inp@VIOLET.BERKELEY.EDU)
  1271. @.br
  1272. \fBmessier.star\fP - Bob Tidd and Alan Paeth
  1273. @.SH SOURCES
  1274. \fBplanet\fP - \fIAstronomical Formulae for Calculators\fP by Jean Meesus
  1275. @.br
  1276. \fBmoonphase\fP - \fIPractical Astronomy with Your Calculator\fP by Duffett-Smith.
  1277. @.br
  1278. \fBepoch\fP -  \fICelestial BASIC\fP by Eric Burgess (SYBEX 1982)
  1279. @//E*O*F planet.1//
  1280. chmod u=rwx,g=rwx,o=rwx planet.1
  1281.  
  1282. echo x - yaleformat
  1283. sed 's/^@//' > "yaleformat" <<'@//E*O*F yaleformat//'
  1284. Overview
  1285.  
  1286. The Yale Bright Star Catalog is an electronic catalog providing extensive
  1287. stellar information for the brightest stars. It is large (megabytes), and
  1288. includes extensive detail for each star (e.g. B-V color temperature, proper
  1289. motion, catalogue names), much of which is not needed by charting software.
  1290. It is normally distributed on magnetic tape, and has not (to my knowledge)
  1291. been posted to the net.
  1292.  
  1293. The highly reduced and hand-revised versions included with the "starchart"
  1294. software family have appeared in two forms. The first (~160K, posted February,
  1295. 1987) included merely location, magnitude and familiar names for a small number
  1296. of stars. This recent reduction (~190K, posted October (?), 1987) adds spectral
  1297. designation for potential color applications, double and variable star
  1298. flags, Bayer and Flamsteed indices and increases the number of named stars.
  1299.  
  1300. The current dataset was reduced directly from the original Yale data by
  1301. Robert Tidd, using custom software. His file format is given below [1]. The
  1302. raw output was hand-edited by Alan Paeth to remove a number of anomolies [2].
  1303. Users wishing to update or extend yale.star should look in part [3].
  1304.  
  1305. Disclaimer
  1306.  
  1307. The original Yale Bright Star Catalog may not be used for commercial purposes
  1308. without express previous written consent from the Department of Astronomy,
  1309. Yale University.
  1310.  
  1311. --
  1312. Part [1]
  1313.  
  1314. File format for yale.star (reduced Yale star catalog, version 2)
  1315.  
  1316. Name    Col     Len    Type     Desc
  1317. -------------------------------------------------
  1318. ra    1    6    i6    ra  2000 hhmmss
  1319. dec    7    5    +-i4    dec 2000 -mmss
  1320. mag    12    3    {-}i    magnitude*100 -dd or -ddd
  1321. type    15    2    c2    object type
  1322. type1    15    1    c1    object type 1    / star,ga,cl,neb,etc
  1323. type2    16    1    c1    object type 2    / dbl,open,cl,etc
  1324. spec2    17    2    c1d1    spectral 2 'G2'
  1325. letter    19    2    c2    greek Bayer letter(s), or Flamsteed number
  1326. const    21    3    c3    constellation (I.A.U. designation)
  1327. name    24    30    c*    name,descrip.
  1328. newline    54    1    c1    newline (max final loc)
  1329. $
  1330. --------
  1331. 'type1'/'type' 2 coding table
  1332. CG-Glob.Cluster        CO-Open  Cluster     GC-Galac.Cluster
  1333. GP-Sphere Galaxy    GS-Spiral Galaxy
  1334. ND-Difuse Nebula    NP-Planetary Nebula
  1335. P*-Planet
  1336. SS-star            SB-Binary Star        SD-Double Star
  1337. SV-Variable Star
  1338. VM-vector move        VS-vector draw (solid)    VD-vector draw (dotted)
  1339. VH-vector draw (hyphens = dashed)
  1340. I*-Invisible (for annotation)
  1341. --------
  1342. 'letter' coding table
  1343.    
  1344. # Flamsteed
  1345. a-alpha b-beta g-gamma d-delta e-epsilon z-zeta h-eta @-theta i-iota k-kappa
  1346. l-lambda m-mu n-nu E-xi o-omicron p-pi r-rho s-sigma t-tau u-upsilon 0-phi
  1347. x-chi %-psi w-omega $-????
  1348. --------
  1349.  
  1350. Example:
  1351.  
  1352. ---yale.star---
  1353. 141540+1911-11SSK2a BOOArcturus
  1354. 062357-5241-72SSF0a CAR
  1355. 051641+4600006SDG8a AURCapella
  1356. @.
  1357. @.
  1358. @.
  1359. ---messier.star---
  1360. 053430+2202840ND    TAUm1  Crab Nebula
  1361. 175629-1901600CO    SGRm23 
  1362. 202355+3832710CO    CYGm29 ,good at low pwr
  1363. 084026+1959370CO    CNCm44 Praesepe, Beehive Cl
  1364. @.
  1365. @.
  1366. @.
  1367.  
  1368. [2] Hand edits done on "raw" reduced yale.star
  1369.  
  1370. [a] Proper names added:
  1371.  
  1372. 062357-5241-72SSF0a CARCanopus
  1373. 050229+4105000SSK5z AURCapella
  1374. 195047+0852075SDA7a AQLAltair
  1375. 173336-3706162SSB1l SCOShaula
  1376. 125602+3819288SDB9a2CVNCor Caroli
  1377.  
  1378. [b] Magnitude and Spectral data changes (based on "Norton"):
  1379.  
  1380. 122636-6306158SDB1a1CRU
  1381. 122636-6306080SDB2a1CRUAcrux
  1382.  
  1383. [c] Deleted (0.00 magnitude - obviously wrong. The position is in the extreme
  1384.          S.W. of Vela, near non-descript NGC 2547)
  1385.  
  1386. 081334-5012000SSK6
  1387.  
  1388. [d] Name Changed
  1389.  
  1390. 052617+2836165SSB7b TAUHath (old)
  1391. 052617+2836165SSB7b TAUAlnath (new)
  1392.  
  1393. Note: beta Taurii == lambda Aurigae (above)
  1394.  
  1395. [e] Mag limit:
  1396.            _
  1397. 064509-1643146SDA1a CMASirius        (-1.46 not representable)
  1398. 064509-1643-99SDA1a CMASirius
  1399.  
  1400. Sirus has been reduced from mag -1.46 to -.99, owing to limitations in
  1401. the three character positions used to give magnitude. This may be rectified
  1402. by using the old style format (style is checked on a line and not a file basis),
  1403. but this would remove spectral information and identification.
  1404.  
  1405. [f] Known Problems and "Features"
  1406.  
  1407. For certain double stars (when the companion is rather bright), duplicate
  1408. entries exist in the dataset, Ra and Dec values matching exactly. The
  1409. plotting of the second visible double star in an exact overstrike position
  1410. can cause problems, particularly with PostScript, which will inadvertently
  1411. form an annulus symbols which resembles the variable star indicator.
  1412.  
  1413. The fix is to remove all duplicate entries from yale.star which are exact
  1414. matches in RA and DECL, omitting those stars of lower value, while insuring
  1415. that the parent star has a "SD" (star double) indicator. 
  1416.  
  1417. This has been done to roughly seventy stars in the database, notibly to
  1418. Orion's belt and his Trapezium, and to Alpha Centaurii. The omitted entries
  1419. are included in "yale.omit". These remain useful for checking the combined
  1420. magnitudes of bright double stars, or in searching for pairs with striking
  1421. difference in star spectral class. The best known example of the latter in the
  1422. Northern Hemisphere is beta Cygnii [Alberio], but as the locations of both
  1423. companions differ by 2" in RA, the parent star of spectral class K5 is
  1424. represented, together with its B8 compainion.
  1425.  
  1426. [g] Bayer and Flamsteed numbers
  1427.  
  1428. A two character field represents Bayer letters, which by convention are
  1429. assigned as Greek lower case letters from alpha to omega, followed by
  1430. assignments in single and double digit Roman type. The database contains
  1431. entries such as " A" and "A " and "a ", and it is not always clear that the
  1432. lower case letter in the first column is an exact "trigger" for Greek fonts.
  1433. (Because entries such as "@" indicate a theta, and "E" also occurs).
  1434. In practice, this is not a problem unless one desires an authoratative code
  1435. for a quite dim star (ie, beyond the 24th Bayer assignment). For virtally all
  1436. stars given labels in standard Atlases, the StarChart output provides an exact
  1437. match.
  1438.  
  1439.  
  1440. [3] User Extentions to yale.star
  1441.  
  1442. Barring the correction of obvious errors in yale.star, users may wish to
  1443. add additional stars or ephemeride data in .star format. In any case, data
  1444. should be for the epoch E2000.0. The program "epoch" may be used if this is
  1445. not the case. The data records may be in either the new or old reduced
  1446. file format convention (planet.star still produces data in the old format).
  1447. The old format omits provision for spectral class or secondary label
  1448. information, but allow stars of magnitude at or below 10.0, as four (and not
  1449. three) characters are allocated to magnitude. Because testing of format is
  1450. done on a line-by-line basis, entries may be mixed and matched. Thus, the
  1451. present yale.star could be extended with data for dim, nameless stars below
  1452. magnitude 10. Note that yale.star must be sorted in order of decreasing
  1453. magnitude. This is only important for stellar data records.
  1454.  
  1455. For other cases, data produced in .star format may be included for printing
  1456. under the "-f file.star" option in the starcharting software. The important
  1457. fields are for Ra, Decl, Magnitude, Type (code and subcode) and name, as each
  1458. directly controls the location and appearence of the output. Two addition
  1459. record types have been added to facilitate general plotting (such as lunar
  1460. limb profiles or the ecliptic). The "I" type is an invisible object of some
  1461. specified magnitude and location which serves as a means to annotate output.
  1462. The "V" type is a vector object with subcodes "M", "S", "D", and "H" for move,
  1463. solid-draw, dot-draw and hyphen(dash)-draw. Annotations are ignored for vector
  1464. records, because the task of disentangling the intermingling of text and vector
  1465. commands to a large class of output devices is formidable. Again, magnitude
  1466. comes into play in performing clipping of objects too dim for display. A sample
  1467. file of these record types appears in "ephem.star"
  1468. @//E*O*F yaleformat//
  1469. chmod u=rwx,g=rwx,o=rwx yaleformat
  1470.  
  1471. echo x - yale.omit
  1472. sed 's/^@//' > "yale.omit" <<'@//E*O*F yale.omit//'
  1473. /* double star duplicates (in ra and decl, <mag omitted) */
  1474.  
  1475. 003133-6258453SDA2b2TUC
  1476. 004953+2743629SDF265PSC
  1477. 013947-5612600SDK0P ERI
  1478. 015331+1917483SDA g1ARI
  1479. 020202+0246523SDA a PSCKaitain
  1480. 025913+2120555SDA2e ARI
  1481. 030052+5221674SDB9
  1482. 034835-3738542SDA0
  1483. 035417-0257633SDA132ERI
  1484. 044335-0848682SDF255ERI
  1485. 053201-0018686SDB2d ORIMintaka
  1486. 053508+0956556SDB0l ORI
  1487. 053845-0236662SDB3
  1488. 054046-0157421SDB3z ORI
  1489. 061137+4843682SDA041AUR
  1490. 062346+0436672SDF5 8MON
  1491. 064813+5542633SDF6
  1492. 072022-5219660SDF2
  1493. 072852-3151724SDB3
  1494. 073419-2328601SDF6
  1495. 073436+3153285SDA1a GEMCastor
  1496. 073522-7417726SDB9
  1497. 074529-1442680SDA0 2PUP
  1498. 082618-3904728SDA0
  1499. 082640+2432764SDF624CNC
  1500. 082647+2656632SDA402CNC
  1501. 083551+0637715SDG5
  1502. 085530-0758691SDA317HYA
  1503. 093046-3153721SDA0z1ANT
  1504. 094706-6504603SDF0u CAR
  1505. 101959+1951190SDG7g2LEO
  1506. 105537+2445630SDA154LEO
  1507. 111811+3132487SDG0E UMA
  1508. 113216-2916586SDF7
  1509. 124140-0127500SDF0g VIR
  1510. 132356+5456395SDA z UMAMizar
  1511. 143936-6050170SDK1a CENRigil Kentaurus
  1512. 144044+1625581SDA p2BOO
  1513. 144109+1343483SDA2z BOO
  1514. 144459+2705512SDA2e BOOPulcherrima
  1515. 150507-4703482SDB5p LUP
  1516. 152312+3017608SDG2h CRB
  1517. 153449+1032516SDF0d SER
  1518. 153840-0848650SDF7
  1519. 153923+3639600SDB6z1CRB
  1520. 154754-6527652SDA5
  1521. 155654-3358573SDA0E2LUP
  1522. 160422-1123477SDF5E SCO
  1523. 160805+1703652SDK2k HERMarsik
  1524. 160834-3906671SDA0
  1525. 161441+3352666SDG1s CRB
  1526. 162440-2942594SDG0
  1527. 162535-2327522SDB3r OPH
  1528. 163614+5256553SDB917DRA
  1529. 170520+5428583SDF6m DRA
  1530. 171521-2636533SDK136OPH
  1531. 171801-2418680SDF6o OPH
  1532. 172341+3708547SDB9r HER
  1533. 175905-3016700SDG8
  1534. 180305-0811604SDF2t OPH
  1535. 180750+2606590SDA3  HER
  1536. 183323-3844655SDB8k1CRA
  1537. 184421+3940602SDA4e1LYR
  1538. 184423+3936537SDA5e2LYR
  1539. 190625-3703501SDF8g CRA
  1540. 191205+4951675SDG4
  1541. 210404-0549731SDA312AQR
  1542. 214408+2844608SDF3m2CYG
  1543. 222850-0001442SDF2z2AQR
  1544.  
  1545. /* triple star duplicates */
  1546.  
  1547. 081213+1739620SDG2z2CNC
  1548. 081213+1739626SDG0z1CNC
  1549.  
  1550. 062849-0702522SDB3b MON
  1551. 062849-0702540SDB3b
  1552.  
  1553. /* quadruple star duplicates (the Trapezium) */
  1554.  
  1555. 053517-0523670SDB0@1ORI
  1556. 053517-0523673SDB0@1ORI
  1557. 053517-0523810SDB3@1ORI
  1558.  
  1559. /* exact duplicates (match in mag and stellar class) -- possible glitches? */
  1560.  
  1561. 030527+2515611SDB752ARI    
  1562. 054046-0157421SDB3
  1563. 130959+1731522SDF5a COM
  1564. 180650-4326577SDA5
  1565. 205139-6226659SDA2
  1566. 235930+3343658SDG0
  1567. @//E*O*F yale.omit//
  1568. chmod u=rw,g=r,o=r yale.omit
  1569.  
  1570. exit 0
  1571.  
  1572.