[comp.sources.unix] v12i073: StarChart program and Yale star data, Part03/07

rsalz@uunet.UU.NET (Rich Salz) (11/30/87)

Submitted-by: awpaeth@watcgl.waterloo.edu (Alan W. Paeth)
Posting-number: Volume 12, Issue 73
Archive-name: starchart/part03

# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by watcgl!awpaeth on Mon Oct  5 19:01:46 EDT 1987
# Contents:  planet.c planet.1 yaleformat yale.omit
 
echo x - planet.c
sed 's/^@//' > "planet.c" <<'@//E*O*F planet.c//'
/*  Planets
 *           A program to determine the position of
 *	     the Sun., Mer., Ven., Mars, Jup, Sat & Ura
 *
 *  reference Jean Meesus's book, " Astronomical Formulae for
 *  Calculators
 *			program by
 *
 *			F.T. Mendenhall
 *			ihnp4!inuxe!fred
 *
 *		Copyright 1985 by F.T. Mendenhall
 *		all rights reserved
 *
 * (This copyright prevents you from selling the program - the
 *  author grants anyone the right to copy and install the program
 *  on any machine it will run on)
 !
 ! Modified by Alan Paeth, awpaeth@watcgl, January 1987 for
 !   (1) non-interactive mode (uses current GMT)
 !   (2) output in simplified Yale format for use with star charter software
 !   (3) lines 365-380 rewritten as a simplier C expression (awpaeth, Apr '87)
 !
 ! Modified by Petri Launiainen, pl@intrin.FI, February 1987 for
 !   easier LOGFILE definition/modification in Makefile
 !
 ! Modified by Alan Paeth, September 1987, for command line flags.
 !    The program continues to provide reduced Yale output in the old style,
 !     but new versions of "starchart" accept either. The old format allows
 !     representation of magnitudes below -.99 (useful for planets), but no
 !     spectral or stellar identification, which is not useful for planets.
 !
 ! See "Time Notes" section below if you plan to rewrite the user interface.
 */

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

#include <sys/time.h>	/* for getting current GMT */

#define CURYEAR	1987	/* default year -- needs to be maintained */

#ifndef PLANETFILE
#define PLANETFILE "./planet.star"
#endif

/* crude magnitudes of planets (x100) for output filtering by brightness */

#define MAGSOL 0
#define MAGMER 150
#define MAGVEN 0
#define MAGMAR 100
#define MAGJUP 0
#define MAGSAT 100
#define MAGURA 590
#define MAGNEP 800

double pie, rad;
double htod(), atof(), kepler(), truean();
double longi(), lati(), poly(), aint(), range();
int cls(), trans(), speak();
FILE *logfile;
char *progname;

/*
 * Times Notes.
 *
 * This software would ideally use a generic UNIX system call to which would
 * provide day and date info as integers, plus current time west of GMT. These
 * values would become the default settings for each of the -m -d -y -t & -z.
 * Unfortunately, time on across many UNIX systems does not (to my knowledge)
 * exist as a uniform subroutine call.
 *
 * Instead, we presume that the very low level "gettimeofday" (GMT in seconds,
 * from 1 Jan 1970) exists on all UNIX installations, and use this subroutine
 * to use the present time of day when no command line parameters are present,
 * as defaults for the current time have not been filled in.
 *
 * When any flag-style command line parameters are present, hardcoded defaults
 * are substitude for -m, -d, -y -t switches, which may then be overridden.
 * The -z flag is filled in by a companion return value from "gettimeofday".
 * This approach makes for poor (on account of hard coding) defaults.
 *
 * A worthwhile change would be to simplify the command line control structure
 * and default mechanism by using a higher level system call, or (safer),
 * write the routine to convert Unix GMT into day, month, year, so that it
 * lives within this program. This would allow for more reasonable defaults,
 * and still merely one UNIX system call to get GMT time. Perhaps sources
 * from a UNIX system may be studied to do this correctly.
 *
 * Expertise in "time" on a specific UNIX system is *NOT* what is required --
 * what *IS* required is knowledge and access to enough independent UNIX
 * systems to insure (to a high probability) that the code remains portable.
 * In other words, only low-level system calls known to work across a large
 * range of systems should be employed.
 *
 */

double currenttime()
    {
#define GMT1970 2440587.5
#define SECSPERDAY (60.0*60.0*24.0)
    struct timeval tv;
    struct timezone tz;
    gettimeofday(&tv, &tz);
    return(GMT1970 + tv.tv_sec/SECSPERDAY);
    }

double commandtime(argc, argv)
    char **argv;
    {
    struct timeval tv;
    struct timezone tz;
    int mm,dd,yy,aa1,bb1;
    double jd,year,month,day,time,timez;
    int j;
    static char *usage =
"\nusage:\tplanet {current time used}\nor\tplanet juliandate\nor\tplanet [-t hh.mm -z zone -y year -m mon -d day ]";
/*
 * check command line style
 */
    if (argc == 1) return (currenttime());
 /*
  * no args - use current time
  */
    if  ((argv[1][0] != '-'))
	{
/*
 * one numeric arg - use it as the julian date
 */
	if (argc == 2) return(atof(argv[1]));
	die("no switches - %s", usage);
	}
    else
	{
/*
 * flags present, set up defaults
 */
    gettimeofday(&tv, &tz);
    time = 0.0;				/* local time (0.0-23.999) */
    yy = CURYEAR;			/* local year */
    mm = 0;				/* local month */
    dd = 0;				/* local day of month */
    timez = tz.tz_minuteswest/60.0;	/* local time hrs west of Greenwich */
/*
 * override defaults
 */
	for (j=1; j<argc; j++)
	    {
	    if (argv[j][0] != '-')
		die("unknown argument - %s", usage /*argv[j]*/ );
	    switch (argv[j][1])
		{
		case 't': time = htod(argv[++j]); break;
		case 'z': timez = htod(argv[++j]); break;
		case 'y': yy = atoi(argv[++j]); break;
		case 'm': mm = atoi(argv[++j]); break;
		case 'd': dd = atoi(argv[++j]); break;
		default:  die("unknown switch - %s", usage /* argv[j] */ );
		}
	    if (j == argc) die("trailing command line flag - %s", argv[j-1]);
	    }
	}
    if (!mm || !dd) die("both month and day flags must be present");
    day = ((float)dd + (time + timez)/24.0);
    if (mm > 2)
	{
	year = yy;
	month = mm;
	}
    else
	{
	year = yy - 1;
	month = mm + 12;
	}
    aa1 = (int) (year/100);
    bb1 = 2 - aa1 + (int)(aa1/4);
    jd = aint(365.25*year) + aint(30.6001*(month + 1));
    jd = jd  + day + 1720994.5;
    if((year + month/100) > 1582.10) jd = jd + bb1;
    return(jd);
    }

main(argc, argv)
    int argc;
    char **argv;
    {
	double jd, T0;
	double l,a,e,i,w1,w2,M,M0,M1,M2,M4,M5,M6,M7,M8;
	double a0,a1,a2,a3,RA,DEC;
	double ECC,nu,r,u,ll,b,lonpert,radpert;
	double esun,Lsun,Cen,Stheta,Snu,Sr;
	double N,D,epli,thapp,omeg;
	double nu2,P,Q,S,V,W,ze,l1pert,epert,w1pert,apert;
	double psi,H,G,eta,th;

#define WRITEMODE "w"
#define OPENFAIL 0
	progname = argv[0];
	if ((logfile = fopen(PLANETFILE, WRITEMODE)) == OPENFAIL)
	    die(stderr, "fail on opening log file %s\n", PLANETFILE);
	cls();
	pie = 3.1415926535898;
	rad = pie/180.0;

/*  calculate time T0 from 1900 January 0.5 ET */
	jd = commandtime(argc, argv);
	printf("Julian date: %f\n", jd);
	T0 = (jd - 2415020.0)/36525.0;

/* Calculate orbital elements for Mercury */
	a0=178.179078;
	a1=149474.07078;
	a2=0.0003011;
	a3=0.0;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 0.3870986;
	a0 = 0.20561421;
	a1 = 0.00002046;
	a2 = -0.000000030;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 7.002881;
	a1 = 0.0018608;
	a2 = -0.0000183;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 28.753753;
	a1 = 0.3702806;
	a2 = 0.0001208;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 47.145944;
	a1 = 1.1852083;
	a2= 0.0001739;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
	M1 = 102.27938 + 149472.51529*T0 + 0.000007*T0*T0;
	M1 = range(M1);
	
/* solving Kepler find the eccentric anomaly ECC    */

	ECC =kepler(e,M1);
	nu = truean(e,ECC);
	r = a*(1.0 - e * cos(ECC*pie/180.0));
	u = l + nu - M1 - w2;
	ll = longi(w2,i,u);
	b = lati(u,i);
	
/* Now make corrections due to perturbations */
	M2 = 212.60322 + 58517.80387*T0 +0.001286*T0*T0;
	M2 = range(M2);
	M4 = 319.51913 + 19139.85475*T0 + 0.000181*T0*T0;
	M4 = range(M4);
	M5 = 225.32833 + 3034.69202*T0 - 0.000722*T0*T0;
	M5 = range(M5);
	M6 = 175.46622 +1221.55147*T0 - 0.000502*T0*T0;
	M6 = range(M6);
	lonpert = 0.00204*cos((5*M2-2*M1+12.220)*rad)
		 +0.00103*cos((2*M2-M1-160.692)*rad)
		 +0.00091*cos((2*M5-M1-37.003)*rad)
		 +0.00078*cos((5*M2-3*M1+10.137)*rad);
		
	radpert = 0.000007525*cos((2*M5-M1+53.013)*rad)
		 +0.000006802*cos((5*M2-3*M1-259.918)*rad)
		 +0.000005457*cos((2*M2-2*M1-71.188)*rad)
		 +0.000003569*cos((5*M2-M1-77.75)*rad);
		
		r = r + radpert;
		ll = ll + lonpert;
		
/* Now find the Longi and radius vector of the sun */
		
		M= 358.47583 + 35999.04975*T0 - 0.000150*T0*T0
		   -0.0000033*T0*T0*T0;
		M = range(M);
		esun = 0.01675104 - 0.0000418*T0 -0.000000126*T0*T0;
		Lsun = 279.69668 + 36000.76892*T0 + 0.0003025*T0*T0;
		Lsun = range(Lsun);
		Cen= (1.919460-0.004789*T0-0.000014*T0*T0)*sin(M*rad)
		     +(0.020094 - 0.000100*T0)*sin(2*M*rad)
		     +0.000293*sin(3*M*rad);
		
		Stheta= Lsun + Cen;
		Stheta = range(Stheta);
		Snu = M + Cen;
		Sr = 1.0000002*(1.0 - esun*esun)/(1.0 + esun*cos(Snu*rad));

		omeg= 259.18 - 1934.142*T0;
		thapp = Stheta - 0.00569 - 0.00479* sin(omeg*rad);
		epli = 23.452294 - 0.0130125*T0
		       -0.00000164*T0*T0 + 0.000000503*T0*T0*T0
			+0.00256*cos(omeg*rad);
		printf("\nOBJECT     R.A.           DEC    \t DISTANCE\n");
		N = cos(epli*rad)*sin(thapp*rad);
		D = cos(thapp*rad);
		RA = atan2(N,D)/rad;
		DEC = asin(sin(epli*rad)*sin(thapp*rad))/rad;
		speak(RA,DEC,Sr, MAGSOL, "PS", "Sol");
/* tansformation of coordinates on Mercury and output */
		trans(r,b,ll,Stheta,Sr,epli, MAGMER, "PM", "Mercury");
		
/* Now start on Venus */
	a0 = 342.767053;
	a1 = 58519.21191;
	a2 = 0.0003097;
	a3 = 0.0;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 0.7233316;
	a0 = 0.00682069;
	a1 = -0.00004774;
	a2 = 0.000000091;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 3.393631;
	a1 = 0.0010058;
	a2 = -0.0000010;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 54.384186;
	a1 = 0.5081861;
	a2 = -0.0013864;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 75.779647;
	a1 = 0.8998500;
	a2 = 0.0004100;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
/* Venus has a long period pert that needs to be added before Kelper */
	lonpert = 0.00077 *sin((237.24 + 150.27*T0)*rad);
	l = l + lonpert;
	M0 = M2 + lonpert;
	ECC = kepler(e,M0);
	nu = truean(e,ECC);
	r = a*(1.0 - e * cos(ECC*rad));
	u = l + nu - M0- w2;
	ll = longi(w2,i,u);
	b = lati(u,i);
	
/* now Venus perturbations */
	
	lonpert = 0.00313*cos((2*M-2*M2 -148.225)*rad)
		 +0.00198*cos((3*M-3*M2 +2.565)*rad)
		 +0.00136*cos((M-M2-119.107)*rad)
		 +0.00096*cos((3*M-2*M2-135.912)*rad)
		 +0.00082*cos((M5-M2-208.087)*rad);
	ll = ll + lonpert;
	
	radpert = 0.000022501 * cos((2*M-2*M2-58.208)*rad)
		 +0.000019045 * cos((3*M-3*M2+92.577)*rad)
		 +0.000006887 * cos((M5-M2-118.090)*rad)
		 +0.000005172 * cos((M-M2-29.110)*rad)
		 +0.000003620 * cos((5*M-4*M2-104.208)*rad)
		 +0.000003283 * cos((4*M-4*M2+63.513)*rad)
		 +0.000003074 * cos((2*M5-2*M2-55.167)*rad);
	r = r + radpert;
	trans(r,b,ll,Stheta,Sr,epli, MAGVEN, "PV", "Venus");
	
/* Now  start the planet Mars */
	a0 = 293.737334;
	a1 = 19141.69551;
	a2 = 0.0003107;
	a3 = 0.0;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 1.5236883;
	a0 = 0.09331290;
	a1 = 0.000092064;
	a2 = -0.000000077;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 1.850333;
	a1 = -0.0006750;
	a2 = 0.0000126;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 285.431761;
	a1 = 1.0697667;
	a2 = 0.0001313;
	a3 = 0.00000414;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 48.786442;
	a1 = 0.7709917;
	a2 = -0.0000014;
	a3 = -0.00000533;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);

/* Mars has a long period perturbation */
	lonpert = -0.01133*sin((3*M5-8*M4 +4*M)*rad)
		  -0.00933*cos((3*M5-8*M4 +4*M)*rad);
	l = l + lonpert;
	M0 = M4 + lonpert;
	ECC = kepler(e,M0);
	nu = truean(e,ECC);
	r = a*(1.0 - e * cos(ECC*rad));
	u = l + nu - M0- w2;
	ll = longi(w2,i,u);
	b = lati(u,i);
	
/* Now Mars Perturbations */	
	lonpert = 0.00705*cos((M5-M4-48.958)*rad)
		 +0.00607*cos((2*M5-M4-188.350)*rad)
		 +0.00445*cos((2*M5-2*M4-191.897)*rad)
		 +0.00388*cos((M-2*M4+20.495)*rad)
		 +0.00238*cos((M-M4+35.097)*rad)
		 +0.00204*cos((2*M-3*M4+158.638)*rad)
		 +0.00177*cos((3*M4-M2-57.602)*rad)
		 +0.00136*cos((2*M-4*M4+154.093)*rad)
		 +0.00104*cos((M5+17.618)*rad);
	ll = ll + lonpert;
	
	radpert = 0.000053227*cos((M5-M4+41.1306)*rad)
		 +0.000050989*cos((2*M5-2*M4-101.9847)*rad)
		 +0.000038278*cos((2*M5-M4-98.3292)*rad)
		 +0.000015996*cos((M-M4-55.555)*rad)
		 +0.000014764*cos((2*M-3*M4+68.622)*rad)
		 +0.000008966*cos((M5-2*M4+43.615)*rad);
/*
 * broken into two pieces for wimpy C compilers
 */
	radpert+= 0.000007914*cos((3*M5-2*M4-139.737)*rad)
		 +0.000007004*cos((2*M5-3*M4-102.888)*rad)
		 +0.000006620*cos((M-2*M4+113.202)*rad)
		 +0.000004930*cos((3*M5-3*M4-76.243)*rad)
		 +0.000004693*cos((3*M-5*M4+190.603)*rad)
		 +0.000004571*cos((2*M-4*M4+244.702)*rad)
		 +0.000004409*cos((3*M5-M4-115.828)*rad);
	r = r + radpert;
	trans(r,b,ll,Stheta,Sr,epli, MAGMAR, "PM", "Mars");
	
/* Now start Jupiter */
	a0 = 238.049257;
	a1 = 3036.301986;
	a2 = 0.0003347;
	a3 = -0.00000165;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 5.202561;
	a0 = 0.04833475;
	a1 = 0.000164180;
	a2 = -0.0000004676;
	a3 = -0.0000000017;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 1.308736;
	a1 = -0.0056961;
	a2 = 0.0000039;
	a3 = 0.0;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 273.277558;
	a1 = 0.5994317;
	a2 = 0.00070405;
	a3 = 0.00000508;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 99.443414;
	a1 = 1.0105300;
	a2 = 0.00035222;
	a3 = -0.00000851;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
/*
	Now start perturbation calclations */
	
	nu2 = T0/5.0 +0.1;
	P = 237.47555 +3034.9061*T0;
	Q = 265.91650 + 1222.1139*T0;
	S = 243.51721 + 428.4677*T0;
	V = 5.0*Q -2.0*P;
	W = 2.0*P - 6.0*Q +3.0*S;
	ze = Q - P;
	psi = S - Q;
	
	P = range(P)*rad;
	Q = range(Q)*rad;
	S = range(S)*rad;
	V = range(V)*rad;
	W = range(W)*rad;
	ze = range(ze)*rad;
	psi = range(psi)*rad;
	
l1pert = (0.331364 - 0.010281*nu2 - 0.004692*nu2*nu2)*sin(V)
	+(0.003228 - 0.064436*nu2 + 0.002075*nu2*nu2)*cos(V)
	-(0.003083 + 0.000275*nu2 - 0.000489*nu2*nu2)*sin(2*V)
	+0.002472*sin(W)
	+0.013619*sin(ze)
	+0.018472*sin(2*ze)
	+0.006717*sin(3*ze)
	+0.002775*sin(4*ze)
	+(0.007275 - 0.001253*nu2)*sin(ze)*sin(Q)
	+0.006417*sin(2*ze)*sin(Q)
	+0.002439*sin(3*ze)*sin(Q);
	
l1pert = l1pert	-(0.033839 + 0.001125*nu2)*cos(ze)*sin(Q)
	-0.003767*cos(2*ze)*sin(Q)
	-(0.035681 + 0.001208*nu2)*sin(ze)*cos(Q)
	-0.004261*sin(2*ze)*cos(Q)
	+0.002178*cos(Q)
	+(-0.006333 + 0.001161*nu2)*cos(ze)*cos(Q)
	-0.006675*cos(2*ze)*cos(Q)
	-0.002664*cos(3*ze)*cos(Q)
	-0.002572*sin(ze)*sin(2*Q)
	-0.003567*sin(2*ze)*sin(2*Q)
	+0.002094*cos(ze)*cos(2*Q)
	+0.003342*cos(2*ze)*cos(2*Q);
	
epert =  (.0003606 + .0000130*nu2 - .0000043*nu2*nu2)*sin(V)
	+(.0001289 - .0000580*nu2)*cos(V)
	-.0006764*sin(ze)*sin(Q)
	-.0001110*sin(2*ze)*sin(Q)
	-.0000224*sin(3*ze)*sin(Q)
	-.0000204*sin(Q)
	+(.0001284 + .0000116*nu2)*cos(ze)*sin(Q)
	+.0000188*cos(2*ze)*sin(Q)
	+(.0001460 + .0000130*nu2)*sin(ze)*cos(Q)
	+.0000224*sin(2*ze)*cos(Q)
	-.0000817*cos(Q);
	
epert = epert	+.0006074*cos(ze)*cos(Q)
	+.0000992*cos(2*ze)*cos(Q)
	+.0000508*cos(3*ze)*cos(Q)
	+.0000230*cos(4*ze)*cos(Q)
	+.0000108*cos(5*ze)*cos(Q)
	-(.0000956 + .0000073*nu2)*sin(ze)*sin(2*Q)
	+.0000448*sin(2*ze)*sin(2*Q)
	+.0000137*sin(3*ze)*sin(2*Q)
	+(-.0000997 + .0000108*nu2)*cos(ze)*sin(2*Q)
	+.0000480*cos(2*ze)*sin(2*Q);

epert = epert	+.0000148*cos(3*ze)*sin(2*Q)
	+(-.0000956 +.0000099*nu2)*sin(ze)*cos(2*Q)
	+.0000490*sin(2*ze)*cos(2*Q)
	+.0000158*sin(3*ze)*cos(2*Q)
	+.0000179*cos(2*Q)
	+(.0001024 + .0000075*nu2)*cos(ze)*cos(2*Q)
	-.0000437*cos(2*ze)*cos(2*Q)
	-.0000132*cos(3*ze)*cos(2*Q);
	
w1pert = (0.007192 - 0.003147*nu2)*sin(V)
	+(-0.020428 - 0.000675*nu2 + 0.000197*nu2*nu2)*cos(V)
	+(0.007269 + 0.000672*nu2)*sin(ze)*sin(Q)
	-0.004344*sin(Q)
	+0.034036*cos(ze)*sin(Q)
	+0.005614*cos(2*ze)*sin(Q)
	+0.002964*cos(3*ze)*sin(Q)
	+0.037761*sin(ze)*cos(Q);
	
w1pert = w1pert
	+0.006158*sin(2*ze)*cos(Q)
	-0.006603*cos(ze)*cos(Q)
	-0.005356*sin(ze)*sin(2*Q)
	+0.002722*sin(2*ze)*sin(2*Q)
	+0.004483*cos(ze)*sin(2*Q)
	-0.002642*cos(2*ze)*sin(2*Q)
	+0.004403*sin(ze)*cos(2*Q)
	-0.002536*sin(2*ze)*cos(2*Q)
	+0.005547*cos(ze)*cos(2*Q)
	-0.002689*cos(2*ze)*cos(2*Q);
	
lonpert = l1pert -(w1pert/e);
l = l + l1pert;
M0 = M5 + lonpert;
e = e + epert;
w1 = w1 + w1pert;

apert =  -.000263*cos(V)
	+.000205*cos(ze)
	+.000693*cos(2*ze)
	+.000312*cos(3*ze)
	+.000147*cos(4*ze)
	+.000299*sin(ze)*sin(Q)
	+.000181*cos(2*ze)*sin(Q)
	+.000204*sin(2*ze)*cos(Q)
	+.000111*sin(3*ze)*cos(Q)
	-.000337*cos(ze)*cos(Q)
	-.000111*cos(2*ze)*cos(Q);
	
a = a + apert;
ECC = kepler(e,M0);
nu = truean(e,ECC);
r = a*(1.0 - e * cos(ECC*rad));
u = l + nu - M0 - w2;
ll = longi(w2,i,u);
b = lati(u,i);

trans(r,b,ll,Stheta,Sr,epli, MAGJUP, "PJ", "Jupiter");

/* Now start Saturn */

	a0 = 266.564377;
	a1 = 1223.509884;
	a2 = 0.0003245;
	a3 = -0.0000058;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 9.554747;
	a0 = 0.05589232;
	a1 = -0.00034550;
	a2 = -0.000000728;
	a3 = 0.00000000074;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 2.492519;
	a1 = -0.0039189;
	a2 = -0.00001549;
	a3 = 0.00000004;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 338.307800;
	a1 = 1.0852207;
	a2 = 0.00097854;
	a3 = 0.00000992;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 112.790414;
	a1 = 0.8731951;
	a2 = -0.00015218;
	a3 = -0.00000531;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
	
/* Now Saturn's perturbations */
	
l1pert = (-0.814181 + 0.018150*nu2 + 0.016714*nu2*nu2)*sin(V)
	 +(-0.010497 + 0.160906*nu2 - 0.004100*nu2*nu2)*cos(V)
	 +0.007581*sin(2*V)
	 -0.007986*sin(W)
	 -0.148811*sin(ze)
	 -0.040786*sin(2*ze)
	 -0.015208*sin(3*ze)
	 -0.006339*sin(4*ze)
	 -0.006244*sin(Q);
l1pert = l1pert
	+(0.008931 + 0.002728*nu2)*sin(ze)*sin(Q)
	-0.016500*sin(2*ze)*sin(Q)
	-0.005775*sin(3*ze)*sin(Q)
	+(0.081344 + 0.003206*nu2)*cos(ze)*sin(Q)
	+0.015019*cos(2*ze)*sin(Q)
	+(0.085581 + 0.002494*nu2)*sin(ze)*cos(Q)
	+(0.025328 - 0.003117*nu2)*cos(ze)*cos(Q);
l1pert = l1pert
	+0.014394*cos(2*ze)*cos(Q)
	+0.006319*cos(3*ze)*cos(Q)
	+0.006369*sin(ze)*sin(2*Q)
	+0.009156*sin(2*ze)*sin(2*Q)
	+0.007525*sin(3*psi)*sin(2*Q)
	-0.005236*cos(ze)*cos(2*Q)
	-0.007736*cos(2*ze)*cos(2*Q)
	-0.007528*cos(3*psi)*cos(2*Q);
	
epert = (-.0007927 + .0002548*nu2 +.0000091*nu2*nu2)*sin(V)
	+(.0013381 + .0001226*nu2 -.0000253*nu2*nu2)*cos(V)
	+(.0000248 - .0000121*nu2)*sin(2*V)
	-(.0000305 + .0000091*nu2)*cos(2*V)
	+.0000412*sin(2*ze)
	+.0012415*sin(Q)
	+(.0000390 -.0000617*nu2)*sin(ze)*sin(Q)
	+(.0000165 - .0000204*nu2)*sin(2*ze)*sin(Q)
	+.0026599*cos(ze)*sin(Q)
	-.0004687*cos(2*ze)*sin(Q);
epert = epert
	-.0001870*cos(3*ze)*sin(Q)
	-.0000821*cos(4*ze)*sin(Q)
	-.0000377*cos(5*ze)*sin(Q)
	+.0000497*cos(2*psi)*sin(Q)
	+(.0000163 - .0000611*nu2)*cos(Q)
	-.0012696*sin(ze)*cos(Q)
	-.0004200*sin(2*ze)*cos(Q)
	-.0001503*sin(3*ze)*cos(Q)
	-.0000619*sin(4*ze)*cos(Q)
	-.0000268*sin(5*ze)*cos(Q);
epert = epert
	-(.0000282 + .0001306*nu2)*cos(ze)*cos(Q)
	+(-.0000086 + .0000230*nu2)*cos(2*ze)*cos(Q)
	+.0000461*sin(2*psi)*cos(Q)
	-.0000350*sin(2*Q)
	+(.0002211 - .0000286*nu2)*sin(ze)*sin(2*Q)
	-.0002208*sin(2*ze)*sin(2*Q)
	-.0000568*sin(3*ze)*sin(2*Q)
	-.0000346*sin(4*ze)*sin(2*Q)
	-(.0002780 + .0000222*nu2)*cos(ze)*sin(2*Q)
	+(.0002022 + .0000263*nu2)*cos(2*ze)*sin(2*Q);
epert = epert
	+.0000248*cos(3*ze)*sin(2*Q)
	+.0000242*sin(3*psi)*sin(2*Q)
	+.0000467*cos(3*psi)*sin(2*Q)
	-.0000490*cos(2*Q)
	-(.0002842 + .0000279*nu2)*sin(ze)*cos(2*Q)
	+(.0000128 + .0000226*nu2)*sin(2*ze)*cos(2*Q)
	+.0000224*sin(3*ze)*cos(2*Q)
	+(-.0001594 + .0000282*nu2)*cos(ze)*cos(2*Q)
	+(.0002162 - .0000207*nu2)*cos(2*ze)*cos(2*Q)
	+.0000561*cos(3*ze)*cos(2*Q);
epert = epert
	+.0000343*cos(4*ze)*cos(2*Q)
	+.0000469*sin(3*psi)*cos(2*Q)
	-.0000242*cos(3*psi)*cos(2*Q)
	-.0000205*sin(ze)*sin(3*Q)
	+.0000262*sin(3*ze)*sin(3*Q)
	+.0000208*cos(ze)*cos(3*Q)
	-.0000271*cos(3*ze)*cos(3*Q)
	-.0000382*cos(3*ze)*sin(4*Q)
	-.0000376*sin(3*ze)*cos(4*Q);

w1pert = (0.077108 + 0.007186*nu2 - 0.001533*nu2*nu2)*sin(V)
	+(0.045803 - 0.014766*nu2 - 0.000536*nu2*nu2)*cos(V)
	-0.007075*sin(ze)
	-0.075825*sin(ze)*sin(Q)
	-0.024839*sin(2*ze)*sin(Q)
	-0.008631*sin(3*ze)*sin(Q)
	-0.072586*cos(Q)
	-0.150383*cos(ze)*cos(Q)
	+0.026897*cos(2*ze)*cos(Q)
	+0.010053*cos(3*ze)*cos(Q);
w1pert = w1pert
	-(0.013597 +0.001719*nu2)*sin(ze)*sin(2*Q)
	+(-0.007742 + 0.001517*nu2)*cos(ze)*sin(2*Q)
	+(0.013586 - 0.001375*nu2)*cos(2*ze)*sin(2*Q)
	+(-0.013667 + 0.001239*nu2)*sin(ze)*cos(2*Q)
	+0.011981*sin(2*ze)*cos(2*Q)
	+(0.014861 + 0.001136*nu2)*cos(ze)*cos(2*Q)
	-(0.013064 + 0.001628*nu2)*cos(2*ze)*cos(2*Q);
	
lonpert = l1pert -(w1pert/e);
l = l + l1pert;
M0 = M6 + lonpert;
e = e + epert;
w1 = w1 + w1pert;

apert = .000572*sin(V) -.001590*sin(2*ze)*cos(Q)
	+.002933*cos(V) -.000647*sin(3*ze)*cos(Q)
	+.033629*cos(ze) -.000344*sin(4*ze)*cos(Q)
	-.003081*cos(2*ze) +.002885*cos(ze)*cos(Q)
	-.001423*cos(3*ze) +(.002172 + .000102*nu2)*cos(2*ze)*cos(Q)
	-.000671*cos(4*ze) +.000296*cos(3*ze)*cos(Q)
	-.000320*cos(5*ze) -.000267*sin(2*ze)*sin(2*Q);
apert = apert
	+.001098*sin(Q) -.000778*cos(ze)*sin(2*Q)
	-.002812*sin(ze)*sin(Q) +.000495*cos(2*ze)*sin(2*Q)
	+.000688*sin(2*ze)*sin(Q) +.000250*cos(3*ze)*sin(2*Q);
apert = apert
	-.000393*sin(3*ze)*sin(Q)
	-.000228*sin(4*ze)*sin(Q)
	+.002138*cos(ze)*sin(Q)
	-.000999*cos(2*ze)*sin(Q)
	-.000642*cos(3*ze)*sin(Q)
	-.000325*cos(4*ze)*sin(Q)
	-.000890*cos(Q)
	+.002206*sin(ze)*cos(Q);
apert = apert
	-.000856*sin(ze)*cos(2*Q)
	+.000441*sin(2*ze)*cos(2*Q)
	+.000296*cos(2*ze)*cos(2*Q)
	+.000211*cos(3*ze)*cos(2*Q)
	-.000427*sin(ze)*sin(3*Q)
	+.000398*sin(3*ze)*sin(3*Q)
	+.000344*cos(ze)*cos(3*Q)
	-.000427*cos(3*ze)*cos(3*Q);
	
a = a + apert;
ECC = kepler(e,M0);
nu = truean(e,ECC);
r = a*(1.0 - e * cos(ECC*rad));
u = l + nu - M0 - w2;
ll = longi(w2,i,u);
b = lati(u,i);

b = b
	+0.000747*cos(ze)*sin(Q)
	+0.001069*cos(ze)*cos(Q)
	+0.002108*sin(2*ze)*sin(2*Q)
	+0.001261*cos(2*ze)*sin(2*Q)
	+0.001236*sin(2*ze)*cos(2*Q)
	-0.002075*cos(2*ze)*cos(2*Q);

trans(r,b,ll,Stheta,Sr,epli, MAGSAT, "PS", "Saturn");

/* Now Start on Uranus */
	a0 = 244.197470;
	a1 = 429.863546;
	a2 = 0.0003160;
	a3 = -0.00000060;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 19.21814;
	a0 = 0.0463444;
	a1 = -0.00002658;
	a2 = 0.000000077;
	a3 = 0.0;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 0.772464;
	a1 = 0.0006253;
	a2 = 0.0000395;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 98.071581;
	a1 = 0.9857650;
	a2 = -0.0010745;
	a3 = -0.00000061;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 73.477111;
	a1 = 0.4986678;
	a2 = 0.0013117;
	a3 = 0.0;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
	M7 = 72.64878 + 428.37911*T0 + 0.000079*T0*T0;
	M7 = range(M7);
/* now perturbation corrections for Uranus */
	G = 83.76922 + 218.4901*T0;
	S = S/rad;
	P = P/rad;
	Q = Q/rad;
	H = 2.0*G - S;
	ze = S - P;
	eta = S - Q;
	th = G - S;
	S = S*rad;
	G = range(G)*rad;
	P = P*rad;
	Q = Q*rad;
	H = range(H)*rad;
	ze = range(ze)*rad;
	eta = range(eta)*rad;
	th = range(th)*rad;
	
l1pert = (0.864319 - 0.001583*nu2)*sin(H)
	+(0.082222 - 0.006833*nu2)*cos(H)
	+0.036017*sin(2*H)
	-0.003019*cos(2*H)
	+0.008122*sin(W);
	
epert = (-.0003349 + .0000163*nu2)*sin(H)
	+.0020981*cos(H)
	+.0001311*cos(H);
	
w1pert = 0.120303*sin(H)
	+(0.019472 - 0.000947*nu2)*cos(H)
	+0.006197*sin(2*H);
	
lonpert = l1pert -(w1pert/e);
l = l + l1pert;
M0 = M7 + lonpert;
e = e + epert;
w1 = w1 + w1pert;

a = a - 0.003825*cos(H);
ECC = kepler(e,M0);
nu = truean(e,ECC);
r = a*(1.0 - e * cos(ECC*rad));
u = l + nu - M0 - w2;
ll = longi(w2,i,u);
b = lati(u,i);

ll = ll
	+(0.010122 - 0.000988*nu2)*sin(S+eta)
	+(-0.038581 + 0.002031*nu2 - 0.001910*nu2*nu2)*cos(S+eta)
	+(0.034964 - 0.001038*nu2 + 0.000868*nu2*nu2)*cos(2*S+eta)
	+0.005594*sin(S +3*th);
ll = ll
	-0.014808*sin(ze)
	-0.005794*sin(eta)
	+0.002347*cos(eta)
	+0.009872*sin(th)
	+0.008803*sin(2*th)
	-0.004308*sin(3*th);
	
b = b
	+(0.000458*sin(eta) - 0.000642*cos(eta) - 0.000517*cos(4*th))
	*sin(S)
	-(0.000347*sin(eta) + 0.000853*cos(eta) + 0.000517*sin(4*eta))
	*cos(S)
	+0.000403*(cos(2*th)*sin(2*S) + sin(2*th)*cos(2*S));
	
r = r
	-.025948
	+.004985*cos(ze)
	-.001230*cos(S)
	+.003354*cos(eta)
	+(.005795*cos(S) - .001165*sin(S) + .001388*cos(2*S))*sin(eta)
	+(.001351*cos(S) + .005702*sin(S) + .001388*sin(2*S))*cos(eta)
	+.000904*cos(2*th)
	+.000894*(cos(th) - cos(3*th));
trans(r,b,ll,Stheta,Sr,epli, MAGURA, "PU", "Uranus");

/* Now start Neptune */
	a0 = 84.457994;
	a1 = 219.885914;
	a2 = 0.0003205;
	a3 = -0.00000060;
	l = poly(a0,a1,a2,a3,T0);
	l = range(l);
	a = 30.10957;
	a0 = 0.00899704;
	a1 = 0.000006330;
	a2 = -0.000000002;
	a3 = 0.0;
	e = poly(a0,a1,a2,a3,T0);
	a0 = 1.779242;
	a1 = -0.0095436;
	a2 = -0.0000091;
	i = poly(a0,a1,a2,a3,T0);
	a0 = 276.045975;
	a1 = 0.3256394;
	a2 = 0.00014095;
	a3 = 0.000004113;
	w1 = poly(a0,a1,a2,a3,T0);
	w1 = range(w1);
	a0 = 130.681389;
	a1 = 1.0989350;
	a2 = 0.00024987;
	a3 = -0.000004718;
	w2 = poly(a0,a1,a2,a3,T0);
	w2 = range(w2);
	M8 = 37.73063 + 218.46134*T0 -0.000070*T0*T0;
	
/* now perturbation corrections for neptune */
	G = G/rad;
	P = P/rad;
	Q = Q/rad;
	S = S/rad;
	ze = G - P;
	eta = G - Q;
	th = G - S;
	G = G*rad;
	P = P*rad;
	Q = Q*rad;
	S = S*rad;
	ze = range(ze)*rad;
	eta = range(eta)*rad;
	th = range(th)*rad;
	
l1pert = (-0.589833 + 0.001089*nu2)*sin(H)
	+(-0.056094 + 0.004658*nu2)*cos(H)
	-0.024286*sin(2*H);
	
epert = .0004389*sin(H)
	+.0004262*cos(H)
	+.0001129*sin(2*H)
	+.0001089*cos(2*H);
	
w1pert = 0.024039*sin(H)
	-0.025303*cos(H)
	+0.006206*sin(2*H)
	-0.005992*cos(2*H);
	
	
lonpert = l1pert -(w1pert/e);
l = l + l1pert;
M0 = M8 + lonpert;
e = e + epert;
w1 = w1 + w1pert;

a = a - 0.000817*sin(H)
	+0.008189*cos(H)
	+0.000781*cos(2*H);
	
ECC = kepler(e,M0);
nu = truean(e,ECC);
r = a*(1.0 - e * cos(ECC*rad));
u = l + nu - M0 - w2;
ll = longi(w2,i,u);
b = lati(u,i);

ll = ll
	-0.009556*sin(ze)
	-0.005178*sin(eta)
	+0.002572*sin(2*th)
	-0.002972*cos(2*th)*sin(G)
	-0.002833*sin(2*th)*cos(G);
	
b = b
	+0.000336*cos(2*th)*sin(G)
	+0.000364*sin(2*th)*cos(G);
	
r = r
	-.040596
	+.004992*cos(ze)
	+.002744*cos(eta)
	+.002044*cos(th)
	+.001051*cos(2*th);
trans(r,b,ll,Stheta,Sr,epli, MAGNEP, "PN", "Neptune");

putchar('\n');
close(logfile);
exit(0);
} /* end of program main */


int cls()
    {
#define CLEARCNT 1
    int i;
    for (i=0; i < CLEARCNT; i++) putchar('\n');
    }

double poly(a0,a1,a2,a3,t)
    double a0,a1,a2,a3,t;
    {
    return(a0+a1*t+a2*t*t+a3*t*t*t);
    }
double aint(z)
	double z;
{
	int trunk;
	
	trunk = (int)z;
	z = (double) trunk;
	return(z);
}
double kepler(e,M)
	double e,M;
{
	double corr,e0,E0,E1;
	e0 = e /rad;
	corr = 1;
	E0=M;
	while (corr > 0.000001)
	{
		corr = (M + e0* sin(E0*rad) -E0)/(1-e*cos(E0*rad));
		E1 = E0 + corr;
		if (corr < 0)
		{
			corr = -1.0* corr ;
		}
		E0 = E1;
	}
	
	return(E1);
}

double truean(e,E)
	double e,E;
{
	double nu;

	nu = 2.0 * atan(sqrt((1+e)/(1-e))*tan(E*rad/2));
	nu = nu/rad;
	if (nu < 0.0)
	{
		nu = nu + 360.0;
	}
	if (fabs(nu-E) > 90.0)
	{
		if (nu > 180.0)
		{
			nu = nu - 180.0;
		}
		else
		{
			nu = nu + 180.0;
		}
	}
return(nu);
}

double longi(w2,i,u)
	double w2,i,u;
{
	double x,y,l;
	
	y=cos(i*rad)*sin(u*rad);
	x=cos(u*rad);
	l = atan2(y,x);
	l = l/rad;
	if (l < 0.0)
	{
		l = l + 360.0;
	}
	return(l+ w2);
}

double lati(u,i)
	double u,i;
{
	double b;
	b = asin(sin(u*rad)*sin(i*rad))/rad;
	return(b);
}

int speak(ra,dec,dis, mag, sym, str)
	double ra,dec,dis;
	char *sym, *str;
{
	int h,m,s,deg;
	
	if ( ra < 0)
	{
		ra = ra + 360.0;
	}
	ra = ra /15.0; /* convert from degs to hours */
	h = (int) aint(ra);
	m = (int)((ra-h)*60);
	s = (int) (((ra-h)*60-m)*60);
	printf("%s\t%2dh %2dm %2ds ",str, h,m,s);
	fprintf(logfile, "%02d%02d%02d", h, m, s);
	deg = (int)aint(dec);
	m   =(int)((dec - deg)*60);
	s   = (int) (((dec - deg)*60 - m)*60);
	if ( m < 0)
	{
		m = m * -1;
	}
	if (s < 0)
	{
		s = s * -1;
	}
	printf("   %+3ddeg %2dm %2ds ",deg,m,s);
	printf("\t  %.3f\n",dis);
	fprintf(logfile, "%+03d%02d+%03d%s%s\n", deg, m, mag, sym, str);
	return(0);
}

double range(val)
	double val;
{
	while (val < 0.0)
	{
		val = val + 360.0;
	}
	if (val > 360.0)
	{
		val = val - (aint(val/360.0)*360.0);
	}
	return(val);
}
int trans(r,b,ll,Stheta,Sr,epli, mag, sym, str)
	double r,b,ll,Stheta,Sr,epli;
	char *sym, *str;
{
	double N,D,lambda,delta,beta,RA,DEC;
	
		N = r * cos(b*rad) * sin((ll - Stheta)*rad);
		D = r * cos(b*rad) * cos((ll - Stheta)*rad) + Sr;
		lambda = atan2(N,D)/rad;
		if (lambda < 0.0)
		{
			lambda = lambda + 360.0;
		}
		lambda = range(lambda+Stheta);
		delta = sqrt(N*N + D*D + (r*sin(b*rad))*(r*sin(b*rad)));
		beta = asin(r*sin(b*rad)/delta)/rad;
		N = sin(lambda*rad)*cos(epli*rad)
		    - tan(beta*rad)*sin(epli*rad);
		D = cos(lambda*rad);
		RA = atan2(N,D)/rad;
		DEC = asin(sin(beta*rad)*cos(epli*rad)
		     +cos(beta*rad)*sin(epli*rad)*sin(lambda*rad))/rad;
		speak(RA,DEC,delta, mag, sym, str);
		return(0);
}

double htod(s)
    char *s;
    {
/*
 * htod(str) reads floating point strings of the form {+-}hh.{mm} thus allowing
 * for fractional values in base 60 (ie, degrees/minutes).
 */
    double x, sign;
    int full, frac;
    char *t;
    t = s-1;
    while(*++t)
	{
	if ( ( (*t<'0') || (*t>'9') ) && (*t!='.') && (*t!='+') && (*t!='-'))
	    die("non-digit in dd.mm style numeric argument");
	}
    if (s == NULL) return 0.0;
    full = frac = 0;
    sign = 1.0;
    if (*s == '-')
	{
	sign = -1.0;
	s++;
	}
    else if (*s == '+') s++;
    while (*s && *s != '.') full = 10 * full + *s++ - '0';
    if (*s++ == '.')
	{
	if (*s) frac = 10 * (*s++ - '0');
	if (*s) frac += *s++ - '0';
	if (frac > 59) frac = 59;
	}
    x = (double) full + ((double) frac) / 60.0;
    return sign * x;
    }

die(a,b)
    char *a, *b;
    {
    fprintf(stderr,"%s: ", progname);
    fprintf(stderr,a,b);
    fprintf(stderr,"\n");
    exit(1);
    }
@//E*O*F planet.c//
chmod u=rwx,g=rwx,o=rwx planet.c
 
echo x - planet.1
sed 's/^@//' > "planet.1" <<'@//E*O*F planet.1//'
@.TH starchart LOCAL 6/2/87
@.ad b
@.SH NAME
planet, moonphase, epoch
@.br
\- print phase of the Moon and planetary information, precess data in .star file format.
@.SH SYNOPSIS
\fBplanet [juliandate]\fR
@.br
or
@.br
\fBplanet [-z timezone -m month -d day -y year -t hh.mm ]\fR
@.br
\fBmoonphase\fR
@.br
\fBepoch [ inputepoch ] [ outputepoch ]\fR
@.SH DESCRIPTION
These programs generate and maintain star chart data in the format used by
the \fBstarchart\fR software tools.
@.PP
Planetary information is generated by the program \fBplanet\fR, which calculates
locations a given Julian date, or for the current time (given no arguments).
In common use, an alternate syntax is provided to simplify the calculation of
the Julian date, in which integers for month, day, (mandatory), and year are
given. Time is specified in the form \fBhh.mm\fP or \fBhh\fP on a 24-hour clock
and defaults to 0.00 (midnight). The \fB-z\fP flag (same format as -t) gives
the number of hours the observer is west of GMT (see below).
@.PP
A summary is then printed on the standard output giving the right ascension,
declination and distance (in AU's) for the Sun and all planets except Pluto.
This data is additionally placed in the file \fBplanet.star\fR for subsequent
use by charting software.
@.PP
The program \fBmoonphase\fP gives the current lunar phase. At present it does
not provide coordinate data for subsequent plotting.
@.PP
The program \fBepoch\fP converts right-ascension and declination data within
\fByale.star\fR style datasets from epoch 1975.0 into the epoch 2000.0.
Either the first or both epoch dates may be overridden with command line
parameters (floating point values). The data is read and written from
the standard input and output. Because only the coordinate data is updated
within each line item, the program may be used on either versions of the
reduced Yale data.
@.SH BUGS
Daylight savings time may have to be reckoned for in \fBplanet\fP, as the
default \fB-z\fP value during periods of daylight time may or may not be
adjusted relative to GMT. The required code changes would make \fBplanet\fP
installation dependent.
@.SH AUTHORS
\fBplanet\fP - F.T. Mendenhall (ihnp4!inuxe!fred)
@.br
\fBmoonphase\fP - Keith E. Brandt
@.br
\fBepoch\fP - Alan Paeth, University of Waterloo (AWPaeth@watCGL)
@.br
\fByale.star\fP - Robert Tidd (inp@VIOLET.BERKELEY.EDU)
@.br
\fBmessier.star\fP - Bob Tidd and Alan Paeth
@.SH SOURCES
\fBplanet\fP - \fIAstronomical Formulae for Calculators\fP by Jean Meesus
@.br
\fBmoonphase\fP - \fIPractical Astronomy with Your Calculator\fP by Duffett-Smith.
@.br
\fBepoch\fP -  \fICelestial BASIC\fP by Eric Burgess (SYBEX 1982)
@//E*O*F planet.1//
chmod u=rwx,g=rwx,o=rwx planet.1
 
echo x - yaleformat
sed 's/^@//' > "yaleformat" <<'@//E*O*F yaleformat//'
Overview

The Yale Bright Star Catalog is an electronic catalog providing extensive
stellar information for the brightest stars. It is large (megabytes), and
includes extensive detail for each star (e.g. B-V color temperature, proper
motion, catalogue names), much of which is not needed by charting software.
It is normally distributed on magnetic tape, and has not (to my knowledge)
been posted to the net.

The highly reduced and hand-revised versions included with the "starchart"
software family have appeared in two forms. The first (~160K, posted February,
1987) included merely location, magnitude and familiar names for a small number
of stars. This recent reduction (~190K, posted October (?), 1987) adds spectral
designation for potential color applications, double and variable star
flags, Bayer and Flamsteed indices and increases the number of named stars.

The current dataset was reduced directly from the original Yale data by
Robert Tidd, using custom software. His file format is given below [1]. The
raw output was hand-edited by Alan Paeth to remove a number of anomolies [2].
Users wishing to update or extend yale.star should look in part [3].

Disclaimer

The original Yale Bright Star Catalog may not be used for commercial purposes
without express previous written consent from the Department of Astronomy,
Yale University.

--
Part [1]

File format for yale.star (reduced Yale star catalog, version 2)

Name    Col     Len    Type     Desc
-------------------------------------------------
ra	1	6	i6	ra  2000 hhmmss
dec	7	5	+-i4	dec 2000 -mmss
mag	12	3	{-}i	magnitude*100 -dd or -ddd
type	15	2	c2	object type
type1	15	1	c1	object type 1	/ star,ga,cl,neb,etc
type2	16	1	c1	object type 2	/ dbl,open,cl,etc
spec2	17	2	c1d1	spectral 2 'G2'
letter	19	2	c2	greek Bayer letter(s), or Flamsteed number
const	21	3	c3	constellation (I.A.U. designation)
name	24	30	c*	name,descrip.
newline	54	1	c1	newline (max final loc)
$
--------
'type1'/'type' 2 coding table
CG-Glob.Cluster		CO-Open  Cluster	 GC-Galac.Cluster
GP-Sphere Galaxy	GS-Spiral Galaxy
ND-Difuse Nebula	NP-Planetary Nebula
P*-Planet
SS-star			SB-Binary Star		SD-Double Star
SV-Variable Star
VM-vector move		VS-vector draw (solid)	VD-vector draw (dotted)
VH-vector draw (hyphens = dashed)
I*-Invisible (for annotation)
--------
'letter' coding table
   
# Flamsteed
a-alpha b-beta g-gamma d-delta e-epsilon z-zeta h-eta @-theta i-iota k-kappa
l-lambda m-mu n-nu E-xi o-omicron p-pi r-rho s-sigma t-tau u-upsilon 0-phi
x-chi %-psi w-omega $-????
--------

Example:

---yale.star---
141540+1911-11SSK2a BOOArcturus
062357-5241-72SSF0a CAR
051641+4600006SDG8a AURCapella
@.
@.
@.
---messier.star---
053430+2202840ND    TAUm1  Crab Nebula
175629-1901600CO    SGRm23 
202355+3832710CO    CYGm29 ,good at low pwr
084026+1959370CO    CNCm44 Praesepe, Beehive Cl
@.
@.
@.

[2] Hand edits done on "raw" reduced yale.star

[a] Proper names added:

062357-5241-72SSF0a CARCanopus
050229+4105000SSK5z AURCapella
195047+0852075SDA7a AQLAltair
173336-3706162SSB1l SCOShaula
125602+3819288SDB9a2CVNCor Caroli

[b] Magnitude and Spectral data changes (based on "Norton"):

122636-6306158SDB1a1CRU
122636-6306080SDB2a1CRUAcrux

[c] Deleted (0.00 magnitude - obviously wrong. The position is in the extreme
	     S.W. of Vela, near non-descript NGC 2547)

081334-5012000SSK6

[d] Name Changed

052617+2836165SSB7b TAUHath (old)
052617+2836165SSB7b TAUAlnath (new)

Note: beta Taurii == lambda Aurigae (above)

[e] Mag limit:
           _
064509-1643146SDA1a CMASirius		(-1.46 not representable)
064509-1643-99SDA1a CMASirius

Sirus has been reduced from mag -1.46 to -.99, owing to limitations in
the three character positions used to give magnitude. This may be rectified
by using the old style format (style is checked on a line and not a file basis),
but this would remove spectral information and identification.

[f] Known Problems and "Features"

For certain double stars (when the companion is rather bright), duplicate
entries exist in the dataset, Ra and Dec values matching exactly. The
plotting of the second visible double star in an exact overstrike position
can cause problems, particularly with PostScript, which will inadvertently
form an annulus symbols which resembles the variable star indicator.

The fix is to remove all duplicate entries from yale.star which are exact
matches in RA and DECL, omitting those stars of lower value, while insuring
that the parent star has a "SD" (star double) indicator. 

This has been done to roughly seventy stars in the database, notibly to
Orion's belt and his Trapezium, and to Alpha Centaurii. The omitted entries
are included in "yale.omit". These remain useful for checking the combined
magnitudes of bright double stars, or in searching for pairs with striking
difference in star spectral class. The best known example of the latter in the
Northern Hemisphere is beta Cygnii [Alberio], but as the locations of both
companions differ by 2" in RA, the parent star of spectral class K5 is
represented, together with its B8 compainion.

[g] Bayer and Flamsteed numbers

A two character field represents Bayer letters, which by convention are
assigned as Greek lower case letters from alpha to omega, followed by
assignments in single and double digit Roman type. The database contains
entries such as " A" and "A " and "a ", and it is not always clear that the
lower case letter in the first column is an exact "trigger" for Greek fonts.
(Because entries such as "@" indicate a theta, and "E" also occurs).
In practice, this is not a problem unless one desires an authoratative code
for a quite dim star (ie, beyond the 24th Bayer assignment). For virtally all
stars given labels in standard Atlases, the StarChart output provides an exact
match.


[3] User Extentions to yale.star

Barring the correction of obvious errors in yale.star, users may wish to
add additional stars or ephemeride data in .star format. In any case, data
should be for the epoch E2000.0. The program "epoch" may be used if this is
not the case. The data records may be in either the new or old reduced
file format convention (planet.star still produces data in the old format).
The old format omits provision for spectral class or secondary label
information, but allow stars of magnitude at or below 10.0, as four (and not
three) characters are allocated to magnitude. Because testing of format is
done on a line-by-line basis, entries may be mixed and matched. Thus, the
present yale.star could be extended with data for dim, nameless stars below
magnitude 10. Note that yale.star must be sorted in order of decreasing
magnitude. This is only important for stellar data records.

For other cases, data produced in .star format may be included for printing
under the "-f file.star" option in the starcharting software. The important
fields are for Ra, Decl, Magnitude, Type (code and subcode) and name, as each
directly controls the location and appearence of the output. Two addition
record types have been added to facilitate general plotting (such as lunar
limb profiles or the ecliptic). The "I" type is an invisible object of some
specified magnitude and location which serves as a means to annotate output.
The "V" type is a vector object with subcodes "M", "S", "D", and "H" for move,
solid-draw, dot-draw and hyphen(dash)-draw. Annotations are ignored for vector
records, because the task of disentangling the intermingling of text and vector
commands to a large class of output devices is formidable. Again, magnitude
comes into play in performing clipping of objects too dim for display. A sample
file of these record types appears in "ephem.star"
@//E*O*F yaleformat//
chmod u=rwx,g=rwx,o=rwx yaleformat
 
echo x - yale.omit
sed 's/^@//' > "yale.omit" <<'@//E*O*F yale.omit//'
/* double star duplicates (in ra and decl, <mag omitted) */

003133-6258453SDA2b2TUC
004953+2743629SDF265PSC
013947-5612600SDK0P ERI
015331+1917483SDA g1ARI
020202+0246523SDA a PSCKaitain
025913+2120555SDA2e ARI
030052+5221674SDB9
034835-3738542SDA0
035417-0257633SDA132ERI
044335-0848682SDF255ERI
053201-0018686SDB2d ORIMintaka
053508+0956556SDB0l ORI
053845-0236662SDB3
054046-0157421SDB3z ORI
061137+4843682SDA041AUR
062346+0436672SDF5 8MON
064813+5542633SDF6
072022-5219660SDF2
072852-3151724SDB3
073419-2328601SDF6
073436+3153285SDA1a GEMCastor
073522-7417726SDB9
074529-1442680SDA0 2PUP
082618-3904728SDA0
082640+2432764SDF624CNC
082647+2656632SDA402CNC
083551+0637715SDG5
085530-0758691SDA317HYA
093046-3153721SDA0z1ANT
094706-6504603SDF0u CAR
101959+1951190SDG7g2LEO
105537+2445630SDA154LEO
111811+3132487SDG0E UMA
113216-2916586SDF7
124140-0127500SDF0g VIR
132356+5456395SDA z UMAMizar
143936-6050170SDK1a CENRigil Kentaurus
144044+1625581SDA p2BOO
144109+1343483SDA2z BOO
144459+2705512SDA2e BOOPulcherrima
150507-4703482SDB5p LUP
152312+3017608SDG2h CRB
153449+1032516SDF0d SER
153840-0848650SDF7
153923+3639600SDB6z1CRB
154754-6527652SDA5
155654-3358573SDA0E2LUP
160422-1123477SDF5E SCO
160805+1703652SDK2k HERMarsik
160834-3906671SDA0
161441+3352666SDG1s CRB
162440-2942594SDG0
162535-2327522SDB3r OPH
163614+5256553SDB917DRA
170520+5428583SDF6m DRA
171521-2636533SDK136OPH
171801-2418680SDF6o OPH
172341+3708547SDB9r HER
175905-3016700SDG8
180305-0811604SDF2t OPH
180750+2606590SDA3  HER
183323-3844655SDB8k1CRA
184421+3940602SDA4e1LYR
184423+3936537SDA5e2LYR
190625-3703501SDF8g CRA
191205+4951675SDG4
210404-0549731SDA312AQR
214408+2844608SDF3m2CYG
222850-0001442SDF2z2AQR

/* triple star duplicates */

081213+1739620SDG2z2CNC
081213+1739626SDG0z1CNC

062849-0702522SDB3b MON
062849-0702540SDB3b

/* quadruple star duplicates (the Trapezium) */

053517-0523670SDB0@1ORI
053517-0523673SDB0@1ORI
053517-0523810SDB3@1ORI

/* exact duplicates (match in mag and stellar class) -- possible glitches? */

030527+2515611SDB752ARI	
054046-0157421SDB3
130959+1731522SDF5a COM
180650-4326577SDA5
205139-6226659SDA2
235930+3343658SDG0
@//E*O*F yale.omit//
chmod u=rw,g=r,o=r yale.omit
 
ey")at %