/* "eph.c" star/planet celestial navigation program August 30, 1989 */

// This program is copyright (c) 2016, P. Lutus and is released
// under the GPL (http://www.gnu.org/licenses/gpl-3.0.en.html).

/* Most recent modification: 03.21.2016 fix segfaulting error */

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

#include <string.h>

#define char_bell 7
#define char_bs 8
#define char_fs 21
#define char_del 127
#define char_nl 10
#define char_cr 13
#define char_cz 26

#define cls() myputs("\033[H\033[J")
#define cleol() myputs("\033[K")
#define gotoxy(x,y) printf("%c[%d;%dH",27,y,x)
#define dsr() printf("\033[6n")
#define cpr(s,x,y) sscanfptr = sscanf(s,"\033[%d;%dR",y,x)
#define home() myputs("\033[H")
#define normal() myputs("\033[0m")
#define inverse() myputs("\033[7m")
#define toupper(c) ((c >= 'a') && (c <= 'z')?c-32:c)
#define tolower(c) ((c >= 'A') && (c <= 'Z')?c+32:c)
#define isnum(c) (((c >= '0') && (c <= '9')) || (c == '.'))
#define isalpha(c) (((c >= 'A') && (c <= 'Z')) || ((c >= 'a') && (c <= 'z')))
#define sgn(x) ((x > 0)?1:(x < 0)?-1:0)

#define PI 3.141592653589793
#define PI2 6.283185307179586
#define PI4 12.56637061435917
#define PD2 1.570796326794897
#define DTOR 1.74532925199433e-2
#define RTOD 5.729577951308232e1
#define ASTUNIT 9.2957130e7

const int buflen = 128;
char inbuf[buflen];
char temp[buflen];
char lastin[2][buflen];

char *starstrings[] = {
	"acamar            3.1    315.6232  -0.00942  -40.3847   0.004",
	"achernar          0.6    335.7575  -0.00917  -57.3374   0.005",
	"acrux             1.1    173.6225  -0.0138   -62.9899  -0.0055",
	"adhara            1.6    255.5374  -0.00983  -28.9463  -0.0014",
	"aldebaran         1.1    291.3032  -0.01425   16.4688   0.002",
	"alioth            1.7    166.7119  -0.0108    56.0672  -0.0053",
	"alkaid            1.9    153.3111  -0.0097    49.4103  -0.005",
	"alnair            2.2     28.2531  -0.0156   -47.0563   0.0048",
	"alnilam           1.8    276.1968  -0.0126    -1.2157   0.0006",
	"alphard           2.2    218.3453  -0.0122    -8.5739  -0.004",
	"alphecca          2.3    126.5376  -0.0105    26.7818  -0.0033",
	"alpheratz         2.2    358.1578  -0.0128    28.9814   0.0055",
	"altair            0.9     62.5546  -0.01217    8.8171   0.00267",
	"ankaa             2.4    353.6739  -0.0123   -42.4132   0.0053",
	"antares           1.2    112.9507  -0.0153   -26.3876  -0.00217",
	"arcturus          0.2    146.3103  -0.0113    19.2857  -0.00517",
	"atria             1.9    108.3565  -0.026    -68.9915  -0.00175",
	"avior             1.7    234.4704  -0.0051   -59.4471  -0.0032",
	"bellatrix         1.7    278.9822  -0.0133     6.3307   0.0008",
	"betelgeuse        0.6    271.4740  -0.0134     7.4026   0.0002",
	"canopus          -0.9    264.1210  -0.0056   -52.6869  -0.0006",
	"capella           0.2    281.1924  -0.0183    45.9774   0.0009",
	"deneb             1.3     49.8092  -0.0084    45.2110   0.0036",
	"denebola          2.2    182.9871  -0.0127    14.6819  -0.0055",
	"diphda            2.2    349.3257  -0.0125   -18.0951   0.0054",
	"dubhe             2.0    194.3736  -0.0153    61.8572  -0.0054",
	"elnath            1.8    278.7392  -0.0157    28.5900   0.0008",
	"eltanin           2.4     90.9631  -0.0058    51.4931  -0.0001",
	"enif              2.5     34.1954  -0.0122     9.7850   0.0045",
	"fomalhaut         1.3     15.8601  -0.0138   -29.7264   0.0053",
	"gacrux            1.6    172.4811  -0.0138   -57.0025  -0.0055",
	"gienah            2.8    176.3022  -0.0128   -17.4326  -0.0055",
	"hadar             0.9    149.3918  -0.0176   -60.2775  -0.0048",
	"hamal             2.2    328.4849  -0.014     23.3688   0.0047",
	"kausastralis      2.0     84.2844  -0.0164   -34.3935  -0.0006",
	"kochab            2.2    137.3182   0.0007    74.2374  -0.0041",
	"markab            2.6     14.0554  -0.0124    15.0996   0.0053",
	"menkar            2.8    314.6885  -0.013      4.0117   0.0039",
	"menkent           2.3    148.6199  -0.0146   -36.2728  -0.0048",
	"miaplacidus       1.8    221.7482  -0.0028   -69.6374  -0.0041",
	"mirfak            1.9    309.2719  -0.0178    49.7907   0.0035",
	"nunki             2.1     76.4899  -0.0154   -26.3201   0.0013",
	"peacock           2.1     53.9778  -0.0196   -56.7979   0.0033",
	"pollux            1.2    243.9731  -0.0152    28.0732  -0.0024",
	"procyon           0.5    245.4326  -0.013      5.2746  -0.0026",
	"rasalhague        2.1     96.4954  -0.0115    12.5754  -0.0007",
	"regulus           1.3    208.1699  -0.0133    12.0632  -0.00492",
	"rigel             0.3    281.6028  -0.012     -8.2256   0.0011",
	"rigilkentaurus    0.1    140.4333  -0.017    -60.7521  -0.004",
	"sabik             2.6    102.6886  -0.0143   -15.6999  -0.0012",
	"schedar           2.5    350.1526  -0.0142    56.4293   0.0054",
	"shaula            1.7     96.9326  -0.0169   -37.0886  -0.0008",
	"sirius           -1.6    258.9308  -0.0109   -16.6907  -0.0014",
	"spica             1.2    158.9617  -0.0131   -11.0578  -0.0051",
	"suhail            2.2    223.1813  -0.009    -43.3539  -0.004",
	"vega              0.1     80.9324  -0.0084    38.7667   0.001",
	"zubenelgenubi     2.9    137.5535  -0.0138   -15.9592  -0.0041",
	""
};

struct STARS {
	char name[128];
	double magn;
	double semidiam;
	double dist;
	double sha;
	double shacor;
	double decl;
	double declcor;
};

struct STARS startable[70];

struct COMPLX {
	double rad;
	double lat;
	double lng;
	double x;
	double y;
	double z;
};

struct COMPLX vals[8],loc1,loc2;

struct mytm {
	double tm_sec;
	double tm_min;
	double tm_hour;
	double tm_mday;
	double tm_mon;
	double tm_year;
	double tm_wday;
	double tm_yday;
	double tm_isdst;
};

struct mytm tvals[8];

struct COMPLX here;

struct mytm datime;
double gmtdiff_s;
double gmtdiff_h;
double now;
double mag_dev;
double gha_a;
double horiz = 0.0;
double mag_lim = 2.0;
double eye_ht = 6;
double eye_ht_corr;
double speed;
double bearing;
double timtbl[4];
int ptrtbl[4];
char splitstr[12][128];
char *getsptr;
int sscanfptr;
int useloctm = 0;
int comline = 0;
char tzlbl[8];

int maxstars;

void myputs(char *s)
{	
printf("%s",s);
/*	char *p;
	p = s-1;
	while(*(++p));
	write(1,s,p-s); */
}

double myacos(double x)
{
	if(fabs(x) < 1.0)
		x = acos(x);
	return(x);
}

double myatan2(double a, double b)
{
	if((a != 0.0) || (b != 0.0))
		a = atan2(a,b);
	return(a);
}

char *mygets(char *s,int len)
{
	char *p = NULL;
	if(!comline)
		p = fgets(s,len,stdin);
	return(p);
}

double mdy_sect(struct mytm p)
{
	double r;
	if(p.tm_year > 1900)
		p.tm_year -= 1900;
	p.tm_mon++;
	if(p.tm_mon < 4)
	{
		p.tm_mon += 12;
		p.tm_year -= 1;
	}
	r = p.tm_mday + floor(p.tm_mon * 30.6001) + (p.tm_year * 365.25) - 63;
	r = floor(r);
	r = (r * 24) + p.tm_hour;
	r = (r * 60) + p.tm_min;
	r = (r * 60) + p.tm_sec;
	return(r);
}

struct mytm sect_mdy(double a)
{
	struct mytm r;
	r.tm_sec = fmod(a,60.0);
	a = floor(a/60.0);
	r.tm_min = fmod(a,60.0);
	a = floor(a/60.0);
	r.tm_hour = fmod(a,24.0);
	a = floor(a/24.0);
	a += 63;
	r.tm_year = floor((a -122.1) / 365.25);
	r.tm_mon = floor((a - floor(365.25 * r.tm_year)) / 30.6001);
	r.tm_mday = a - floor(365.25 * r.tm_year);
	r.tm_mday -= floor(30.6001 * r.tm_mon);
	if(r.tm_mon < 14)
		r.tm_mon -= 1;
	else
		r.tm_mon -= 13;
	if (r.tm_mon <= 2)
		r.tm_year += 1;
	r.tm_wday = a / 7.0;
	r.tm_wday = floor(7.0 * (r.tm_wday - floor(r.tm_wday)) + .5);
	r.tm_year += 1970;
	return(r);
}

double fpavals[] = { /* first point of aries, 1980 - 1999 */
	98.8256,
	99.5713,
	99.3317,
	99.0926,
	98.8540,
	99.6017,
	99.3641,
	99.1268,
	98.8897,
	99.6382,
	99.4008,
	99.1631,
	98.9250,
	99.6719,
	99.4326,
	99.1929,
	98.9529,
	99.6982,
	99.4579,
	99.2177
};

double meananom[] = { /* solar mean anomaly, 1980 - 1999 */
	-3.7737,
	-3.0452,
	-3.3020,
	-3.5583,
	-3.8140,
	-3.0836,
	-3.3383,
	-3.5927,
	-3.8470,
	-3.1157,
	-3.3702,
	-3.6251,
	-3.8804,
	-3.1507,
	-3.4071,
	-3.6640,
	-3.9212,
	-3.1930,
	-3.4505,
	-3.7078
};

double solperi[] = { /* sun's perihelion point, 1980 - 1999 */
	77.4006,
	77.3835,
	77.3663,
	77.3491,
	77.3320,
	77.3148,
	77.2976,
	77.2805,
	77.2633,
	77.2461,
	77.2289,
	77.2118,
	77.1946,
	77.1774,
	77.1603,
	77.1431,
	77.1260,
	77.1087,
	77.0916,
	77.0744
};

/* solar position in GHA form, so GHA is subtracted
before return to main program, where it is added again (sigh) */

struct COMPLX getsolpos(struct mytm datq,struct COMPLX r)
{
	struct mytm rt;
	double yd,dt,tm,m,ghaa;
	dt = mdy_sect(datq);
	rt = datq;
	rt.tm_mon = 1;
	rt.tm_mday = 1;
	rt.tm_hour = 0;
	rt.tm_min = 0;
	rt.tm_sec = 0;
	yd = mdy_sect(rt); /* get year day */
	yd = dt - yd;
	yd = floor(yd / 86400.0) + 1;
	tm = datq.tm_hour + (datq.tm_min / 60.0) + (datq.tm_sec / 3600.0);
	ghaa = fpavals[((int)datq.tm_year - 80)]
		+ (.985647 * (yd + (tm / 24.0))) + (15 * tm);
	ghaa *= DTOR;
	ghaa = fmod(ghaa+PI4,PI2);
	m = (.985647 * (yd + (tm / 24.0))) + meananom[((int)datq.tm_year - 80)];
	r.lng = m + (1.9160 * sin(m * DTOR)) + (0.02 * sin(m * 2 * DTOR));
	r.lng -= solperi[((int)datq.tm_year - 80)];
	r.lng *= DTOR;
	r.lng = fmod(r.lng+PI4,PI2);
	r.lat = asin(0.3978 * sin(r.lng));
	dt = atan(0.9175 * tan(r.lng));
	if(dt < 0.0) dt += PI;
	if(r.lng >= PI) dt += PI;
	r.lng = (m + (tm * 15)) * DTOR;
	r.lng -= dt;
	r.lng -= ((solperi[((int)datq.tm_year - 80)] + 180.0) * DTOR);
	r.lng -= ghaa;
	r.lng = fmod(r.lng+PI4,PI2);
	r.lng = PI2 - r.lng;
	return(r);
}

/* import planet ephemeris */

struct HMS
{
	unsigned hr,mn,sc;
};

double hms_radians(struct HMS hms)
{
	double digit;
	digit = hms.hr + (hms.mn / 60.0) + (hms.sc / 3600.0);
	digit *= DTOR;
	return(digit);
}

char *planet_name[] ={
	"sun",
	"mercury",
	"venus",
	"mars",
	"jupiter",
	"saturn",
	"uranus",
	"neptune",
	"pluto"
};

double orbital_period[] = {
	
	365.2596, /* Earth */
	87.9693, /* Mercury */
	224.7009, /* Venus */
	686.9796, /* Mars */
	4332.1248, /* Jupiter */
	10825.863, /* Saturn */
	30676.15, /* Uranus */
	59911.13, /* Neptune */
	90824.2 /* Pluto */
};

double eccentricity[] = {
	
	.016755, /* Earth */
	.205636, /* Mercury */
	.006759, /* Venus */
	.093348, /* Mars */
	.048061, /* Jupiter */
	.050822, /* Saturn */
	.047552, /* Uranus */
	.006608, /* Neptune */
	.253364 /* Pluto */
};

struct HMS ascending_node[] = {
	
	0, 0, 0, /* Earth */
		48, 9, 4, /* Mercury */
		76, 32, 42,/* Venus */
		49, 26, 35, /* Mars */
		100, 20, 20, /* Jupiter */
		113, 32, 20, /* Saturn */
		73, 59, 17, /* Uranus */
		131, 37, 34, /* Neptune */
		110, 12, 58   /* Pluto */
};

struct HMS perihelion[] = {
	
	102, 42, 22, /* Earth */
		77, 13, 19, /* Mercury */
		131, 29, 24, /* Venus */
		335, 43, 16, /* Mars */
		15, 23, 56, /* Jupiter */
		93, 29, 13, /* Saturn */
		177, 7, 23, /* Uranus */
		354, 40, 12, /* Neptune */
		224, 15, 0   /* Pluto */
};

struct HMS inclination[] = {
	
	0, 0, 0, /* Earth */
		7, 0, 17, /* Mercury */
		3, 23, 40, /* Venus */
		1, 50, 59, /* Mars */
		1, 18, 19, /* Jupiter */
		2, 29, 8, /* Saturn */
		0, 46, 27, /* Uranus */
		1, 46, 15, /* Neptune */
		17, 7, 56   /* Pluto */
};

struct HMS epoch_longitude[] = {
	
	35, 32, 32, /* Earth */
		242, 6, 54, /* Mercury */
		298, 45, 24, /* Venus */
		329, 44, 5, /* Mars */
		293, 28, 58, /* Jupiter */
		224, 21, 44, /* Saturn */
		248, 5, 3, /* Uranus */
		271, 32, 56, /* Neptune */
		216, 55, 28   /* Pluto */
};

double distance[] = {
	
	.999994, /* Earth */
	.387098, /* Mercury */
	.723330, /* Venus */
	1.523712, /* Mars */
	5.20270, /* Jupiter */
	9.57542, /* Saturn */
	19.2998, /* Uranus */
	30.2813, /* Neptune */
	39.7138   /* Pluto */
};

double albedo[] = {
	
	.39, /* Earth */
	.06, /* Mercury */
	.72, /* Venus */
	.16, /* Mars */
	.70, /* Jupiter */
	.75, /* Saturn */
	.90, /* Uranus */
	.82, /* Neptune */
	.14   /* Pluto */
};

double radius[] = { /* miles */
	
	432560.0, /* Sun */
	1515.0, /* Mercury */
	3760.0, /* Venus */
	2108.4, /* Mars */
	44362.0, /* Jupiter */
	37280.0, /* Saturn */
	15800.0, /* Uranus */
	15100.0, /* Neptune */
	930.0   /* Pluto */
};

double epoch = 30981; /* oct 27 1984 */

struct COMPLX pol_rect(struct COMPLX pol)
{
	struct COMPLX r;
	double cp2;
	cp2 = cos(pol.lat);
	r.x = pol.rad * cos(pol.lng) * cp2;
	r.y = pol.rad * sin(pol.lng) * cp2;
	r.z = pol.rad * sin(pol.lat);
	return(r);
}

struct COMPLX rect_pol(struct COMPLX r)
{
	struct COMPLX pol;
	pol.lng = myatan2(r.y,r.x);
	pol.rad = hypot(r.y,r.x);
	if(pol.lng < 0) pol.lng = PI2 + pol.lng;
	pol.lat = myatan2(r.z,pol.rad);
	pol.rad = hypot(r.z,pol.rad);
	return(pol);
}

struct COMPLX cd_position(double cd,int n)
{
	double a,ec,e,oe,qe;
	int i;
	struct COMPLX pol;
	struct COMPLX r;
	a = (cd - epoch)/orbital_period[n];
	a -= floor(a);
	a *= PI2;
	a += hms_radians(epoch_longitude[n]);
	a -= hms_radians(perihelion[n]);
	a = (a < 0)?PI2 + a:a;
	ec = eccentricity[n];
	e = a;
	i = 0;
	do
	{
		oe = e;
		e = a + ec * sin(e); /* compute eccentric anomaly the normal way */
		if (i > 20)
		{
			qe = a + ec * sin(e); /* probably oscillating now ... */
			e = (e + qe) * .5;    /* take average of 2 extremes */
		}
	}
	while((fabs(oe - e) > .00001) && (i++ < 100)); /* eventually give up */
	a = sqrt((1 + ec)/(1 - ec));
	a = 2 * atan(a * tan(e * .5));
	pol.rad = distance[n] * (1 - ec * cos(e));
	a += hms_radians(perihelion[n]);
	pol.lng = a;
	e = a - hms_radians(ascending_node[n]);
	ec = hms_radians(inclination[n]);
	pol.lat = sin(e);
	pol.lat *= ec;
	r = pol_rect(pol);
	return(r);
}

struct COMPLX right_ascension(struct COMPLX arg)
{
	double r,ang;
	r = hypot(arg.z,arg.y);
	ang = myatan2(arg.z,arg.y);
	ang += .40927970959;
	arg.z = r * sin(ang);
	arg.y = r * cos(ang);
	return(arg);
}

void planet_table(double cd,int p,struct mytm datim)
{
	int n,a,b;
	double helio_dist,helio_angle,earth_dist,magnitude;
	struct COMPLX cart,pol;
	static struct COMPLX earth;
	static double oldcd = 0.0;
	cd /= 86400.0;
	if(cd != oldcd)
	{
		earth = cd_position(cd,0);
		earth.x = -earth.x;
		earth.y = -earth.y;
		earth.z = -earth.z;
		cart = right_ascension(earth);
		pol = rect_pol(cart);
		earth_dist = pol.rad;
		if((datim.tm_year >= 80) && (datim.tm_year <= 99))
			pol = getsolpos(datim,pol);
		startable[60].sha = PI2 - pol.lng;
		startable[60].decl = pol.lat;
		startable[60].magn = -26.8;
		startable[60].dist = earth_dist * ASTUNIT;
		startable[60].shacor = 0;
		startable[60].declcor = 0;
		startable[60].semidiam = asin(radius[0]/startable[60].dist);
		oldcd = cd;
	}
	a = 1;
	b = 8;
	if(p >= 0)
	{
		a = p-60;
		b = p-60;
	}
	for (n = a;n <= b;n++)
	{
		if(n)
		{
			cart = cd_position(cd,n);
			pol = rect_pol(cart);
			helio_dist = pol.rad;
			helio_angle = pol.lng;
			cart.x += earth.x;
			cart.y += earth.y;
			cart.z += earth.z;
			cart = right_ascension(cart);
			pol = rect_pol(cart);
			earth_dist = pol.rad;
			helio_angle -= pol.lng;
			startable[n+60].sha = PI2 - pol.lng;
			startable[n+60].decl = pol.lat;
			startable[n+60].dist = earth_dist * ASTUNIT;
			startable[n+60].shacor = 0;
			startable[n+60].declcor = 0;
			startable[n+60].semidiam = asin(radius[n]/startable[n+60].dist);
			magnitude = 2.8048e-1; /* relative sun brightness */
			magnitude *= (radius[n] * radius[n]);
			magnitude /= (helio_dist * helio_dist);
			magnitude *= albedo[n];
			magnitude /= (earth_dist * earth_dist);
			magnitude *= 1.7e-5;
			magnitude *= ((1.0 + cos(helio_angle)) / 2.0);
			magnitude = -2.5 * (log10(magnitude));
			startable[n+60].magn = magnitude;
		}
	} /* for */
}

/* end planet epehemeris import */

void localt(int n)
{
	if(n)
		useloctm ^= 1;
	if(useloctm)
		strcpy(tzlbl,*tzname);
	else
		strcpy(tzlbl,"GMT");
}

void locase(char *s)
{
	while(*s = tolower(*s))
		s++;
}

/* #ifndef MSDOS

double fmod(a,b)
double a,b;
{
	a /= b;
	a -= floor(a);
	a *= b;
	return(a);
}

#endif */

struct COMPLX rect_to_pol(struct COMPLX p)
{
	p.lng = myatan2(p.x,p.z);
	p.rad = hypot(p.x,p.z);
	p.lat = myatan2(p.y,p.rad);
	return(p);
}

struct COMPLX comp_pr(struct COMPLX b,struct COMPLX a)
{
	double dlo,cdlo,sal,cal;
	dlo = a.lng - b.lng;
	cdlo = cos(dlo);
	sal = sin(a.lat);
	cal = cos(a.lat);
	a.rad = myacos((sal * sin(b.lat)) + (cal * cos(b.lat) * cdlo));
	a.rad *= RTOD * 60;
	a.lat = myatan2(sin(dlo),(cal * tan(b.lat) - sal * cdlo));
	if (a.lat < 0) a.lat = PI2 + a.lat;
	return(a);
}

struct COMPLX comp_pa(struct COMPLX b,struct COMPLX a)
{
	a = comp_pr(b,a);
	a.rad = 90.0 - (a.rad / 60.0);
	return(a);
}

struct COMPLX comp_rp(struct COMPLX b,struct COMPLX a)
{
	struct COMPLX c;
	double r,sr,ang;
	r = b.rad / (RTOD * 60.0);
	sr = sin(r);
	c.x = -sin(b.lat) * sr; /* create rect. of dist. & hdg. */
	c.y = cos(b.lat) * sr;
	c.z = cos(r);
	r = hypot(c.z,c.y);
	ang = myatan2(c.y,c.z) + a.lat; /* rotate y,z to orig. lat. */
	c.z = r * cos(ang);
	c.y = r * sin(ang);
	c = rect_to_pol(c);
	c.lng += a.lng; /* add long. */
	return(c);
}

void ztest(double *x)
{
	if(*x == 0.0)
		*x = 1.0e-15;
}

void conv_circles(double lat,double lon,double zen,double *a, double *b,double *c,double *d, double *k)
{
	*a = cos(lon) * cos(lat);
	*b = -sin(lon) * cos(lat);
	*c = sin(lat);
	*d = cos(zen);
	*k  = *d * (*a * *a +  *b * *b +  *c * *c);
}

double rad_dist(double a,double b,double c,double d)
{
	return(((a-b)*(a-b)) + ((c-d)*(c-d)));
}

struct COMPLX do_circles(struct COMPLX aa,struct COMPLX bb)
{
	struct COMPLX rr;
	double a1,b1,c1,d1,k1;
	double a2,b2,c2,d2,k2;
	double r,s,t,l,m,n,p,u,v,w,f,g,h,x,y,z,dd,ee;
	aa.rad = PD2 - aa.rad;
	bb.rad = PD2 - bb.rad;
	if(aa.lat == bb.lat)
		aa.lat += 1e-10;
	conv_circles(aa.lat,aa.lng,aa.rad,&a1,&b1,&c1,&d1,&k1);
	conv_circles(bb.lat,bb.lng,bb.rad,&a2,&b2,&c2,&d2,&k2);
	ztest(&b1);
	ztest(&b2);
	ztest(&c1);
	ztest(&c2);
	r = a1*c2/c1-a2;
	s = b1*c2/c1-b2;
	t = k1*c2/c1-k2;
	m = t/r;
	n = -t/s;
	p = (b1*t/s-a1*t/r)/c1;
	l = sqrt(m*m + n*n + p*p);
	ztest(&l);
	m /= l;
	n /= l;
	p /= l;
	u = a2*p/c2-m;
	v = b2*p/c2-n;
	w = k2*p/c2;
	ztest(&r);
	ztest(&v);
	f = (t*v/s-w)/(r*v/s-u);
	g = (t-r*f)/s;
	h = (k1-a1*f-b1*g)/c1;
	dd = sqrt(f*f + g*g + h*h);
	if((dd*dd) <= 1.0) /* if circles intersect */
	{
		rr.z = 0;
		ee = sqrt(1.0-dd*dd);
		/* first result */
		x = f + ee*m;
		y = g + ee*n;
		z = h + ee*p;
		rr.lat = asin(z);
		rr.lng = myatan2(-y,x);
		/* second result */
		x = f - ee*m;
		y = g - ee*n;
		z = h - ee*p;
		rr.x = asin(z);
		rr.y = myatan2(-y,x);
	}
	else
		rr.z = -1;
	return(rr);
}

void disp_deg_min(double v,int f)
{
	int q_d,q_m,vs;
	double q_m2;
	vs = (v < 0)?-1:1;
	v *= RTOD;
	v = fabs(v) + .00008;
	q_d = (int) v;
	q_m = (int) ((v - q_d) * 6000.0);
	q_m2 = q_m * 1e-2;
	if(f)
		q_d *= vs;
	printf("%3d deg %5.02lf min",q_d,q_m2);
}

void disp_hms(double v)
{
	int q_d,q_m,q_s;
	if(useloctm)
	{
		v -= gmtdiff_h;
		if (v < 0)
			v += 24;
		if(v > 24)
			v -= 24;
	}
	v = fabs(v) + .000138;
	q_d = (int) v;
	v = (v - q_d) * 60.0;
	q_m = (int) v;
	v = (v - q_m) * 60.0;
	q_s = (int) v;
	printf("%02d:%02d:%02d %s",q_d,q_m,q_s,tzlbl);
}

void disp_lat_lng(struct COMPLX q)
{
	int lat_s,lng_s;
	int lat_d,lng_d;
	int lat_m,lng_m;
	double lat_m2,lng_m2;
	q.lat *= RTOD;
	q.lng *= RTOD;
	if(q.lng > 180.0)
		q.lng -= 360.0;
	lat_s = (q.lat >= 0);
	lng_s = (q.lng >= 0);
	q.lat = fabs(q.lat) + .000008;
	q.lng = fabs(q.lng) + .000008;
	lat_d = (int) q.lat;
	lat_m = (int) ((q.lat - lat_d) * 6000.0);
	lat_m2 = lat_m * 1e-2;
	lng_d = (int) q.lng;
	lng_m = (int) ((q.lng - lng_d) * 6000.0);
	lng_m2 = lng_m * 1e-2;
	printf("lat %3d deg %5.02lf min %c lng %3d deg %5.02lf min %c"
	,lat_d,lat_m2,(lat_s)?'N':'S',lng_d,lng_m2,(lng_s)?'W':'E');
}

void print_datime(struct mytm t)
{
	double dt;
	if(useloctm)
	{
		dt = mdy_sect(t);
		t = sect_mdy(dt-gmtdiff_s);
	}
	printf("%02.0lf:%02.0lf:%02.0lf %s %02.0lf/%02.0lf/%04.0lf"
		,t.tm_hour,t.tm_min,t.tm_sec,tzlbl
		,t.tm_mon,t.tm_mday,t.tm_year);
}

void split(char *s)
{
	int i;
	for(i = 0; i < 12;i++)
		splitstr[i][0] = 0;
	sscanfptr = sscanf(s,
		"%s %s %s %s %s %s %s %s %s %s %s %s"
		,splitstr[0],splitstr[1],splitstr[2],splitstr[3]
			,splitstr[4],splitstr[5],splitstr[6],splitstr[7]
		,splitstr[8],splitstr[9],splitstr[10],splitstr[11]);
}

int isdelim(char *s,int c)
{
	int result = 0;
	do
		result = (*s == c);
	while((*(++s)) && (!result));
	return(result);
}

int llscan(int i,double *q)
{
	double x;
	int c;
	*q *= RTOD;
	if(*(splitstr+i)[0])
	{
		sscanf(*(splitstr+i),"%lf",q);
		i++;
	}
	if(*(splitstr+i)[0])
	{
		sscanf(*(splitstr+i),"%lf",&x);
		*q += (x/60.0);
		i++;
	}
	if(*(splitstr+i)[0])
	{
		c = toupper(*(splitstr+i)[0]);
		if ((c == 'S') || (c == 'N') || (c == 'E') || (c == 'W'))
		{
			if ((c == 'S') || (c == 'E'))
				*q = -*q;
			i++;
		}
	}
	*q *= DTOR;
	return(i);
}

int limbscan(int i,int *q)
{
	int c;
	c = toupper(*(splitstr+i)[0]);
	if((c == 'L') || (c == 'U'))
	{
		*q = (c == 'U');
		i++;
	}
	else *q = 0;
	return(i);
}

int getdate(int i,struct mytm *t)
{
	struct mytm q;
	q = *t;
	if(isdelim(splitstr[i],'/'))
	{
		sscanfptr = sscanf(splitstr[i],"%lf/%lf/%lf"
		,&q.tm_mon,&q.tm_mday,&q.tm_year);
		i++;
	}
	*t = q;
	return(i);
}

int gettime(int i,struct mytm *t)
{
	struct mytm q;
	q = *t;
	if(isdelim(splitstr[i],':'))
	{
		sscanfptr = sscanf(splitstr[i],"%lf:%lf:%lf"
		,&q.tm_hour,&q.tm_min,&q.tm_sec);
		i++;
	}
	*t = q;
	return(i);
}

struct mytm getdatime(struct mytm *t)
{
	int i = 0;
	unsigned long tm;
	double dtm;
	if(!comline)
		printf("Enter %s Time & date as [hh:mm:ss] [mm/dd/yy] (N=now):",tzlbl);
	getsptr = mygets(inbuf,buflen);
	if(inbuf[0])
	{
		if(toupper(inbuf[0]) == 'N')
		{
			tm = time(NULL);
			dtm = tm;
			*t = sect_mdy(dtm);
		}
		else
		{
			if(useloctm)
			{
				dtm = mdy_sect(*t);
				*t = sect_mdy(dtm-gmtdiff_s);
			}
			split(inbuf);
			i = gettime(i,t);
			i = getdate(i,t);
			if(useloctm)
			{
				dtm = mdy_sect(*t);
				*t = sect_mdy(dtm+gmtdiff_s);
			}
		}
	}
}

struct COMPLX getlatlng(struct COMPLX *x,int prompt)
{
	int i = 0;
	double a;
	if(prompt)
	{
		if(!comline)
			myputs(
		"Enter estimated position as (lat) dd mm.mm [N/S] (lng) ddd mm.mm [E/W]\n:");
	}
	getsptr = mygets(inbuf,buflen);
	split(inbuf);
	a = x->lat;
	i = llscan(i,&a);
	x->lat = a;
	a = x->lng;
	i = llscan(i,&a);
	x->lng = a;
}

double getghaa(struct mytm datq)
{
	struct mytm r;
	double cf = 99.2124; /* an average aries GHA in case no year value */
	double gha,th,yd,dt;
	dt = mdy_sect(datq);
	r = datq;
	r.tm_mon = 1;
	r.tm_mday = 1;
	r.tm_hour = 0;
	r.tm_min = 0;
	r.tm_sec = 0;
	yd = mdy_sect(r); /* get year day */
	yd = dt - yd;
	yd = floor(yd / 86400.0) + 1;
	if ((datq.tm_year >= 80) && (datq.tm_year <= 99))
		cf = fpavals[((int)datq.tm_year - 80)];
	th = datq.tm_hour + (datq.tm_min / 60.0) + (datq.tm_sec / 3600.0);
	gha = cf + (.985647 * (yd + (th/24.0))) + (15 * th);
	gha = fmod(gha,360.0);
	gha *= DTOR;
	return(gha);
}

double sight_corr(double x,int i,int uprlimb) /* eye height, parallax, limb, refraction */
{
	double a,r;
	r = -eye_ht_corr; /* sextant eye height correction */
	if(startable[i].dist)
	{
		a = cos(x) * 3963; /* small angles mean large parallax */
		a = asin(a/startable[i].dist); /* parallax angle */
		r += a;
	}
	if(uprlimb)
		r -= startable[i].semidiam; /* for upper limb */
	else
		r += startable[i].semidiam; /* for lower limb */
	a = x + 5.23598775598e-2; /* perform standard refraction corr. */
	a = atan(a * 12.0 * RTOD);
	r += (2.821615624e-4 * tan(x - a));
	return(r);
}

void liststars(int f,int p)
{
	int i,j = 0,ds;
	double cv,gha,dec,ghaa;
	if(!comline)
	{
		cls();
		printf(
			"Objects brighter than selected magnitude limit (now %.01lf):\n\n"
		,mag_lim);
	}
	cv = (now / 86400.0) - 29220.0;
	cv /= 365.25;
	ghaa = (f)?0.0:gha_a;
	i = (p)?60:0;
	while(startable[i].name[0])
	{
		if(p)
			planet_table(now,i,datime);
		if(startable[i].magn <= mag_lim)
		{
			gha = ghaa + (startable[i].sha + (startable[i].shacor * cv));
			gha = fmod(gha,PI2);
			dec = startable[i].decl + (startable[i].declcor * cv);
			printf("%-16.16s mag %+5.01lf  dec "
			,startable[i].name,startable[i].magn);
			ds = (dec < 0);
			dec = fabs(dec);
			disp_deg_min(dec,0);
			printf(" %c",(ds)?'S':'N');
			printf("  %s ",(f)?"sha":"gha");
			disp_deg_min(gha,0);
			myputs("\n");
			if((!comline) && (j) && (!(j % 18)))
			{
				myputs("(press return):");
				getsptr = mygets(inbuf,buflen);
			}
			j++;
		}
		i++;
	}
	if(!comline)
		myputs("(press return):");
	getsptr = mygets(inbuf,buflen);
}

void findstars(int p)
{
	int i,j = 0;
	struct COMPLX x;
	double cv;
	if(!comline)
	{
		cls();
		printf(
			"Objects brighter than selected magnitude limit (now %.01lf)\n"
		,mag_lim);
		printf(
			"and above selected horizon limit angle (now %.01lf deg):\n\n"
		,horiz * RTOD);
	}
	cv = (now / 86400.0) - 29220.0; /* data are epoch 1980 */
	cv /= 365.25;
	i = (p)?60:0;
	while(startable[i].name[0])
	{
		if(p)
			planet_table(now,i,datime);
		x.lng = gha_a + (startable[i].sha + (startable[i].shacor * cv));
		x.lng = fmod(x.lng,PI2);
		x.lat = startable[i].decl + (startable[i].declcor * cv);
		x = comp_pa(x,here);
		x.lat += mag_dev;
		x.rad *= DTOR;
		x.rad -= sight_corr(x.rad,i,0); /* in reverse, assume lower limb */
		if((x.rad >= horiz) && (startable[i].magn < mag_lim))
		{
			printf("%-16.16s mag %+5.01lf  alt "
			,startable[i].name,startable[i].magn);
			disp_deg_min(x.rad,1);
			myputs("  az ");
			disp_deg_min(x.lat,0);
			myputs("\n");
			if((!comline) && (j) && (!(j % 18)))
			{
				myputs("(press return):");
				getsptr = mygets(inbuf,buflen);
			}
			j++;
		}
		i++;
	}
	if(!comline)
		myputs("(press return):");
	getsptr = mygets(inbuf,buflen);
}

int srchname(char *str)
{
	char *a,*b;
	int i = 0,found = -1,p = 0;
	while((found < 0) && ((startable[i].name[0]) || (!p)))
	{
		if(startable[i].name[0])
		{
			a = str;
			b = startable[i].name;
			while((*a) && (*b) && (*a == *b))
			{
				a++;
				b++;
			}
			if(*a == 0) found = i;
			i++;
		}
		else if(!p)
		{
			p++;
			i = 60;
		}
	}
	return(found);
}

void disp_fix(struct COMPLX c,int f)
{
	if(f)
		myputs("best result --> ");
	else
		myputs("                ");
	disp_lat_lng(c);
	myputs("\n");
}

void starfix()
{
	int i = 0,j,k,uprlimb;
	double dd,cv,ghaa,datemp;
	struct COMPLX c,r;
	r.lat = 0.0;
	if(!comline)
	{
		cls();
		myputs("Compute position by two sextant sightings.\n\n");
		myputs(
			"1. Greatest accuracy is achieved when two sightings take place at the same"
		);
		putc(char_nl,stdout);
		myputs(
			"time (within minutes), or the vessel is not in motion. If this is not  the"
		);
		putc(char_nl,stdout);
		myputs(
			"case, accurate  results  require  the entry of a  good  position  estimate,"
		);
		putc(char_nl,stdout);
		myputs(
			"magnetic  deviation,  vessel  speed and bearing  (main menu options B - E)."
		);
		myputs("\n\n");
		myputs(
			"2. Greatest accuracy is obtained when the sun or stars are used for sights,"
		);
		putc(char_nl,stdout);
		myputs(
			"their positions are calculated with greater precision than the planets.\n\n"
		);
	}
	while(i < 2)
	{
		if(!comline)
		{
			printf(
			"Enter name of object %d, sextant altitude,\n",i+1);
			printf(
			"lower or upper limb (for sun or planets), time %s, date\n",tzlbl);
			myputs("as name (alt) dd mm.mm [L/U] (time) hh:mm:ss (date) mm/dd/yy\n(type \"=\" to use last entry):");
		}
		getsptr = mygets(inbuf,buflen);
		if(inbuf[0])
		{
			tvals[i] = datime;
			if(useloctm)
			{
				datemp = mdy_sect(tvals[i]);
				tvals[i] = sect_mdy(datemp-gmtdiff_s);
			}
			if((inbuf[0] == '=') && (inbuf[1] == 0) && (lastin[i][0]))
				strcpy(inbuf,lastin[i]);
			strcpy(lastin[i],inbuf);
			split(inbuf);
			getsptr = strcpy(temp,splitstr[0]);
			j = 1;
			j = llscan(j,&vals[i].rad);
			j = limbscan(j,&uprlimb);
			j = gettime(j,&tvals[i]);
			j = getdate(j,&tvals[i]);
			if(useloctm)
			{
				datemp = mdy_sect(tvals[i]);
				tvals[i] = sect_mdy(datemp+gmtdiff_s);
			}
			locase(temp);
			datemp = mdy_sect(tvals[i]);
			timtbl[i] = datemp;
			if((j = srchname(temp)) < 0) /* not found */
			{
				j = 58+i;
				getsptr = strcpy(startable[j].name,temp);
				startable[j].semidiam = 0;
				if(!comline)
				{
					printf(
					"\n\"%s\" not found in internal almanac table. Please enter\n",temp);
					myputs(
					"this object's declination and GHA for the time entered above\n");
					myputs(
					"as (decl) ddd mm.mm N/S (GHA) ddd mm.mm:");
				}
				vals[i].lat = 0;
				vals[i].lng = 0;
				getlatlng(&vals[i],0);
			}
			else
			{
				ghaa = getghaa(tvals[i]);
				if(j >= 60)
					planet_table(datemp,j,tvals[i]);
				cv = (datemp / 86400.0) - 29220.0; /* data are epoch 1980 */
				cv /= 365.25;
				vals[i].lng = ghaa + (startable[j].sha + (startable[j].shacor * cv));
				vals[i].lng = fmod(vals[i].lng,PI2);
				vals[i].lat = startable[j].decl + (startable[j].declcor * cv);
			}
			ptrtbl[i] = j;
			vals[i].rad += sight_corr(vals[i].rad,j,uprlimb);
			i++;
		}
		else i = 8; /* quit */
	}
	if(i == 2)
	{
		if(!comline)
			cls();
		myputs("\nSelected objects, positions, and sextant altitudes:\n");
		for(j = 0;j < i;j++)
		{
			k = ptrtbl[j];
			printf("%d: %-16.16s "
			,j+1,startable[k].name);
			disp_lat_lng(vals[j]);
			myputs("\nCorr. Sext. Alt. ");
			disp_deg_min(vals[j].rad,0);
			myputs(" at ");
			print_datime(tvals[j]);
			myputs("\n");
		}
		if(speed != 0.0)
		{
			if(now != timtbl[0])
			{
				dd = (timtbl[0] - now) / 3600.0; /* hours between pos and sight 1 */
				c.rad = speed * dd; /* shift position of object along course */
				c.lat = bearing - mag_dev;
				loc1 = comp_rp(c,here);
			}
			else loc1 = here;
			if(now != timtbl[1])
			{
				dd = (timtbl[1] - now) / 3600.0; /* hours between pos and sight 2 */
				c.rad = speed * dd; /* shift position of object along course */
				c.lat = bearing - mag_dev;
				loc2 = comp_rp(c,here);
			}
			else loc2 = here;
			c = comp_pr(loc1,vals[0]); /* radius for orig. pos. */
			dd = c.rad;
			c = comp_pr(loc2,vals[0]); /* estimated radius for new pos. */
			if(dd != 0)
			{
				dd = c.rad / dd; /* ratio of old and new radii */
				vals[0].rad = PD2 - ((PD2 - vals[0].rad) * dd);
			}
		}
		else loc2 = here;
		c = do_circles(vals[1],vals[0]);
		if(c.z == 0)
		{
			if(rad_dist(c.lat,here.lat,c.lng,here.lng) < rad_dist(c.x,here.lat,c.y,here.lng))
			{
				r = c;
				j = 0;
			}
			else
			{
				r.lat = c.x;
				r.lng = c.y;
				j = 1;
			}
			myputs(
			"\nPossible positions at time of sighting 2:\n\n");
			disp_fix(c,j == 0);
			c.lat = c.x;
			c.lng = c.y;
			disp_fix(c,j == 1);
		}
		else
			myputs(
		"\nError: position circles do not intersect.\n");
		if(!comline)
		{
			if(r.lat)
			{
				myputs("\nWant to replace estimated position with this result (Y/N):");
				getsptr = mygets(inbuf,buflen);
				if(toupper(inbuf[0]) == 'Y')
					here = r;
			}
			else
			{
				myputs("\n(return):");
				getsptr = mygets(inbuf,buflen);
			}
		}
	}
}

double rshour(double l,double d,double subh)
{
	double r;
	r = sin(subh) - (sin(l) * sin(d));
	r /= cos(l) * cos(d);
	if(fabs(r) < 1.0)
		r = (acos(r) * RTOD) / 15.0; /* expressed as delta hours */
	else r *= 1000.0; /* +- error value */
	return(r);
}

void rise_set(int p)
{
	struct COMPLX q;
	struct mytm t,u;
	int i,j = 0;
	double cv,ghaa,datemp,datempb,tval,setpt,rsa;
	if(!comline)
	{
		cls();
		printf(
			"Objects brighter than selected magnitude limit (now %.01lf):\n\n"
		,mag_lim);
	}
	i = (p)?60:0;
	t = datime;
	t.tm_hour = 0;
	t.tm_min = 0;
	t.tm_sec = 0;
	u = t;
	ghaa = getghaa(t);
	datemp = mdy_sect(t);
	cv = (datemp / 86400.0) - 29220.0; /* data are epoch 1980 */
	cv /= 365.25;
	while(startable[i].name[0])
	{
		if(p)
			planet_table(datemp,i,t);
		if(startable[i].magn <= mag_lim)
		{
			setpt = (i == 60)?-.8333:-.55; /* 50 min below for sun, 33 min others */
			setpt *= DTOR;
			q.lng = ghaa + (startable[i].sha + (startable[i].shacor * cv));
			q.lng = fmod(q.lng,PI2);
			q.lat = startable[i].decl + (startable[i].declcor * cv);
			tval = ((here.lng - q.lng) / PI2) * 24.0;
			if (tval < 0) tval += 24.0;
			if(p)
			{
				u.tm_hour = tval;
				datempb = mdy_sect(u);
				planet_table(datempb,i,u);
				q.lng = ghaa + startable[i].sha;
				q.lng = fmod(q.lng,PI2);
				q.lat = startable[i].decl;
				tval = ((here.lng - q.lng) / PI2) * 24.0;
				if (tval < 0) tval += 24.0;
			}
			tval *= .997269567;
			tval = fmod(tval,24.0);
			printf("%-16.16s mer. pass. ",startable[i].name);
			disp_hms(tval);
			rsa = rshour(here.lat,q.lat,setpt);
			if(rsa >= 1000.0)
				myputs(", always below horizon this date & loc.");
			else if(rsa <= -1000.0)
				myputs(", always above horizon this date & loc.");
			else
			{
				myputs(", rise ");
				cv = tval - rsa;
				if(cv < 0.0)
					cv += 24.0;
				disp_hms(cv);
				myputs(", set ");
				cv = tval + rsa;
				if(cv >= 24.0)
					cv -= 24.0;
				disp_hms(cv);
			}
			myputs("\n");
			if((!comline) && (j) && (!(j % 18)))
			{
				myputs("(press return):");
				getsptr = mygets(inbuf,buflen);
			}
			j++;
		}
		i++;
	}
	if(!comline)
		myputs("(press return):");
	getsptr = mygets(inbuf,buflen);
}

void noonsight()
{
	struct COMPLX sp,ss;
	struct mytm st;
	double qt,ghaa,ssign,ssdeg;
	int j = 0,uprlimb;
	st = datime;
	ss.rad = 0;
	if(!comline)
	{
		cls();
		myputs(
		"Compute position by noon sun sight.\n\n");
		myputs(
		"Enter sextant altitude at meridian passage, whether sun was north or\n");
		printf(
		"south of vessel, whether upper or lower limb, %s time and date\n",tzlbl);
		myputs(
		"of meridian passage as\n");
		myputs(
		"(alt) dd mm.mm [N/S] (limb) [L/U] (time) hh:mm:ss (date) mm/dd/yy\n\n:");
	}
	getsptr = mygets(inbuf,buflen);
	split(inbuf);
	j = llscan(j,&ss.rad);
	j = limbscan(j,&uprlimb);
	if(useloctm)
	{
		qt = mdy_sect(st);
		st = sect_mdy(qt-gmtdiff_s);
	}
	j = gettime(j,&st);
	j = getdate(j,&st);
	if(useloctm)
	{
		qt = mdy_sect(st);
		st = sect_mdy(qt+gmtdiff_s);
	}
	qt = mdy_sect(st);
	ghaa = getghaa(st);
	planet_table(qt,60,st);
	sp.lng = startable[60].sha + ghaa;
	sp.lng = fmod(sp.lng,PI2);
	sp.lat = startable[60].decl;
	ssign = sgn(ss.rad);
	ss.rad = fabs(ss.rad);
	ss.rad += sight_corr(ss.rad,60,uprlimb);
	ssdeg = ss.rad;
	ss.rad *= RTOD;
	ss.rad = (90.0 - ss.rad) * 60.0;
	ss.rad *= ssign;
	ss.lat = PI; /* south, .rad has -sign for north */
	sp = comp_rp(ss,sp);
	myputs("Corr. sext. alt. ");
	disp_deg_min(ssdeg,0);
	printf("\nsun %s of vessel, %s limb"
		,(sgn(ss.rad) > 0)?"north":"south"
		,(uprlimb)?"upper":"lower");
		myputs(" at ");
		print_datime(st);
		myputs("\n");
		myputs("vessel location: ");
		disp_lat_lng(sp);
		if(!comline)
		{
			myputs("\n\nWant to replace estimated position with this result (Y/N):");
			getsptr = mygets(inbuf,buflen);
			if(toupper(inbuf[0]) == 'Y')
				here = sp;
	}
}

void putdat()
{
	FILE *fp;
	fp = fopen("eph.dat","w");
	if(fp != NULL)
	{
		fprintf(fp,"%lf %lf\n",here.lat * RTOD,here.lng * RTOD);
		fprintf(fp,"%lf\n",mag_dev * RTOD);
		fprintf(fp,"%lf\n",speed);
		fprintf(fp,"%lf\n",bearing * RTOD);
		fprintf(fp,"%lf\n",eye_ht);
		fprintf(fp,"%lf\n",mag_lim);
		fprintf(fp,"%lf\n",horiz * RTOD);
		fprintf(fp,"%d\n",useloctm);
		fclose(fp);
	}
	else myputs("unable to save program data.\n");
}

void getdat()
{
	FILE *fp;
	fp = fopen("eph.dat","r");
	if(fp != NULL)
	{
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf %lf",&here.lat,&here.lng);
		here.lat *= DTOR;
		here.lng *= DTOR;
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &mag_dev);
		mag_dev *= DTOR;
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &speed);
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &bearing);
		bearing *= DTOR;
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &eye_ht);
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &mag_lim);
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%lf", &horiz);
		horiz *= DTOR;
		fgets(inbuf,120,fp);
		sscanf(inbuf,"%d",&useloctm);
                fclose(fp);
	}
	localt(0);
}

int menu()
{
	int c;
	cls();
	myputs("Celestial Navigation Program (c) P. Lutus 1989-2005\n\n");
	myputs("A. Enter time & date (now ");
	print_datime(datime);
	myputs(")\n");
	myputs("B. Enter est. position (now ");
	disp_lat_lng(here);
	myputs(")\n");
	printf("C. Enter local magnetic deviation (now %.01lf deg)\n"
	,mag_dev * RTOD);
	printf("D. Enter vessel speed (now %.01lf knots)\n",speed);
	myputs("E. Enter vessel magnetic bearing (now ");
	disp_deg_min(bearing,0);
	myputs(")\n");
	printf("F. Enter sextant eye height (now %.01lf feet)\n",eye_ht);
	printf("G. Enter star/planet list magnitude limit (now %.01lf)\n"
	,mag_lim);
	myputs("H. List navigation stars'   current SHA\n");
	myputs("I. List navigation planets' current SHA\n");
	myputs("J. List navigation stars'   current GHA\n");
	myputs("K. List navigation planets' current GHA\n");
	myputs("L. List navigation stars'   meridian, rise, set times\n");
	myputs("M. List navigation planets' meridian, rise, set times\n");
	printf("N. Enter object-locator horizon limit angle (now %.01lf deg)\n"
	,horiz * RTOD);
	myputs("O. Locate navigation stars at current location and time\n");
	myputs("P. Locate navigation planets at current location and time\n");
	myputs("Q. Compute position by two celestial sightings\n");
	myputs("R. Compute position by noon sun sight\n");
	printf("S. Set display and entry time zone (now %s)\n",tzlbl);
	myputs("T. Exit program\n");
	myputs("\nEnter a letter :");
	mygets(inbuf,buflen);
	c = inbuf[0];
	c = toupper(c);
	return(c);
}

void docom(int c)
{
	int i;
	if(!comline)
	{
		gotoxy(1,24);
		cleol();
	}
	switch (c) {
		case 'A':
			getdatime(&datime);
			now = mdy_sect(datime);
			gha_a = getghaa(datime);
			break;
		case 'B':
			getlatlng(&here,1);
			break;
		case 'C':
			if(!comline)
			{
				cls();
				myputs(
				" Magnetic deviation is listed on most navigation charts. Enter magnetic\n");
				myputs(
				" deviation at estimated position as dd.dd (- for easterly):");
		}
		getsptr = mygets(inbuf,buflen);
		mag_dev *= RTOD;
		sscanfptr = sscanf(inbuf,"%lf",&mag_dev);
		mag_dev *= DTOR;
		break;
		case 'D':
			if(!comline)
				myputs("Enter vessel speed (knots):");
			getsptr = mygets(inbuf,buflen);
			sscanfptr = sscanf(inbuf,"%lf",&speed);
			break;
		case 'E':
			if(!comline)
				myputs("Enter vessel magnetic bearing as ddd [mm.mm]:");
			getsptr = mygets(inbuf,buflen);
			split(inbuf);
			i = llscan(0,&bearing);
			break;
		case 'F':
			if(!comline)
				myputs("Enter sextant eye height (ft):");
			getsptr = mygets(inbuf,buflen);
			sscanfptr = sscanf(inbuf,"%lf",&eye_ht);
			eye_ht_corr = 2.821615624e-4 * sqrt(eye_ht); /* sextant eye height corr */
			break;
		case 'G':
			if(!comline)
				myputs("Enter navigation list magnitude limit:");
			getsptr = mygets(inbuf,buflen);
			sscanfptr = sscanf(inbuf,"%lf",&mag_lim);
			break;
		case 'H':
			liststars(1,0);
			break;
		case 'I':
			liststars(1,1);
			break;
		case 'J':
			liststars(0,0);
			break;
		case 'K':
			liststars(0,1);
			break;
		case 'L':
			rise_set(0);
			break;
		case 'M':
			rise_set(1);
			break;
		case 'N':
			if(!comline)
				myputs("Enter object-locator horizon limit angle:");
			getsptr = mygets(inbuf,buflen);
			horiz *= RTOD;
			sscanfptr = sscanf(inbuf,"%lf",&horiz);
			horiz *= DTOR;
			break;
		case 'O':
			findstars(0);
			break;
		case 'P':
			findstars(1);
			break;
		case 'Q':
			starfix();
			break;
		case 'R':
			noonsight();
			break;
		case 'S':
			localt(1);
			break;
		default:
			break;
	}
}

int init()
{
	int i = 0;
	unsigned long qt;
	getdat();
	while(*(starstrings+i)[0])
	{
		sscanfptr = sscanf(*(starstrings+i),"%s %lf %lf %lf %lf %lf"
			,startable[i].name,&startable[i].magn,&startable[i].sha
			,&startable[i].shacor,&startable[i].decl,&startable[i].declcor);
		startable[i].semidiam = 0.0;
		startable[i].dist = 0.0;
		startable[i].sha *= DTOR;
		startable[i].shacor *= DTOR;
		startable[i].decl *= DTOR;
		startable[i].declcor *= DTOR;
		i++;
	}
	for(i = 0;i < 9;i++)
		getsptr = strcpy(startable[i+60].name,*(planet_name+i));
	startable[69].name[0] = 0;
	qt = time(NULL);
	now = qt;
	datime = sect_mdy(now);
	gha_a = getghaa(datime);
	eye_ht_corr = 2.821615624e-4 * sqrt(eye_ht); /* sextant eye height corr */
	return(i);
}

int ctest(char **c)
{
	return((c[0][0] == '-') && (isalpha(c[0][1])));
}

void proc_coml(int argc,char **argv)
{
	int com;
	comline++;
	while(--argc > 0)
	{
		argv++;
		if(ctest(argv))
		{
			com = toupper(argv[0][1]);
			if(argc > 1)
			{
				inbuf[0] = 0;
				do
				{
					argv++;
					argc--;
					if(!ctest(argv))
					{
						strcat(inbuf,*argv);
						strcat(inbuf," ");
					}
				}
				while((argc > 0) && (!ctest(argv)));
				if(ctest(argv))
				{
					argv--;
					argc++;
				}
			}
		}
		docom(com);
	}
}

int main(int argc,char **argv)
{
	int c = 0;
	tzset();
	maxstars = init();
	gmtdiff_s = (double)timezone;
	gmtdiff_h = gmtdiff_s / 3600.0;
	lastin[0][0] = 0;
	lastin[1][0] = 0;
	if(argc > 1)
		proc_coml(argc,argv);
	else
	{
		while(c != 'T') /* as in "nearly forever" */
		{
			c = menu();
			docom(c);
		}
	}
	putdat();
	return 0;
}