/*  * This file is part of the UCB release of Plan 9. It is subject to the license * terms in the LICENSE file found in the top-level directory of this * distribution and at http://akaros.cs.berkeley.edu/files/Plan9License. No * part of the UCB release of Plan 9, including this file, may be copied, * modified, propagated, or distributed except according to the terms contained * in the LICENSE file. */ #include "astro.h" char*	satlst[] ={	0,}; struct{	double	time;	double	tilt;	double	pnni;	double	psi;	double	ppi;	double	d1pp;	double	peri;	double	d1per;	double	e0;	double	deo;	double	rdp;	double	st;	double	ct;	double	rot;	double	semi;} satl; voidsatels(void){	double ifa[10], t, t1, t2, tinc;	char **satp;	int flag, f, i, n; 	satp = satlst; loop:	if(*satp == 0)		return;	f = open(*satp, 0);	if(f < 0) {		fprint(2, "cannot open %s\n", *satp);		satp += 2;		goto loop;	}	satp++;	rline(f);	tinc = atof(skip(6));	rline(f);	rline(f);	for(i=0; i<9; i++)		ifa[i] = atof(skip(i));	n = ifa[0];	i = ifa[1];	t = dmo[i-1] + ifa[2] - 1.;	if(n%4 == 0 && i > 2)		t += 1.;	for(i=1970; i<n; i++) {		t += 365.;		if(i%4 == 0)			t += 1.;	}	t = (t * 24. + ifa[3]) * 60. + ifa[4];	satl.time = t * 60.;	satl.tilt = ifa[5] * radian;	satl.pnni = ifa[6] * radian;	satl.psi = ifa[7];	satl.ppi = ifa[8] * radian;	rline(f);	for(i=0; i<5; i++)		ifa[i] = atof(skip(i));	satl.d1pp = ifa[0] * radian;	satl.peri = ifa[1];	satl.d1per = ifa[2];	satl.e0 = ifa[3];	satl.deo = 0;	satl.rdp = ifa[4]; 	satl.st = sin(satl.tilt);	satl.ct = cos(satl.tilt);	satl.rot = pipi / (1440. + satl.psi);	satl.semi = satl.rdp * (1. + satl.e0); 	n = PER*288.; /* 5 min steps */	t = day;	for(i=0; i<n; i++) {		if(sunel((t-day)/deld) > 0.)			goto out;		satel(t);		if(el > 0) {			t1 = t;			flag = 0;			do {				if(el > 30.)					flag++;				t -= tinc/(24.*3600.);				satel(t);			} while(el > 0.);			t2 = (t - day)/deld;			t = t1;			do {				t += tinc/(24.*3600.);				satel(t);				if(el > 30.)					flag++;			} while(el > 0.);			if(flag)				if((*satp)[0] == '-')					event("%s pass at ", (*satp)+1, "",						t2, SIGNIF+PTIME+DARK); else					event("%s pass at ", *satp, "",						t2, PTIME+DARK);		}	out:		t += 5./(24.*60.);	}	close(f);	satp++;	goto loop;} voidsatel(double time){	int i;	double amean, an, coc, csl, d, de, enom, eo;	double pp, q, rdp, slong, ssl, t, tp; 	i = 500;	el = -1;	time = (time-25567.5) * 86400;	t = (time - satl.time)/60;	if(t < 0)		return; /* too early for satelites */	an = floor(t/satl.peri);	while(an*satl.peri + an*an*satl.d1per/2. <= t) {		an += 1;		if(--i == 0)			return;	}	while((tp = an*satl.peri + an*an*satl.d1per/2.) > t) {		an -= 1;		if(--i == 0)			return;	}	amean = (t-tp)/(satl.peri+(an+.5)*satl.d1per);	pp = satl.ppi+(an+amean)*satl.d1pp;	amean *= pipi;	eo = satl.e0+satl.deo*an;	rdp = satl.semi/(1+eo);	enom = amean+eo*sin(amean);	do {		de = (amean-enom+eo*sin(enom))/(1.0-eo*cos(enom));		enom += de;		if(--i == 0)			return;	} while(fabs(de) >= 1.0e-6);	q = 3963.35*erad/(rdp*(1-eo*cos(enom))/(1-eo));	d = pp + 2*atan2(sqrt((1+eo)/(1-eo))*sin(enom/2),cos(enom/2));	slong = satl.pnni + t*satl.rot -		atan2(satl.ct*sin(d), cos(d));	ssl = satl.st*sin(d);	csl = pyth(ssl);	if(vis(time, atan2(ssl,csl), slong, q)) {		coc = ssl*sin(glat) + csl*cos(glat)*cos(wlong-slong);		el = atan2(coc-q, pyth(coc));		el /= radian;	}} intvis(double t, double slat, double slong, double q){	double t0, t1, t2, d; 	d = t/86400 - .005375 + 3653;	t0 = 6.238030674 + .01720196977*d;	t2 = t0 + .0167253303*sin(t0);	do {		t1 = (t0 - t2 + .0167259152*sin(t2)) /			(1 - .0167259152*cos(t2));		t2 = t2 + t1;	} while(fabs(t1) >= 1.e-4);	t0 = 2*atan2(1.01686816*sin(t2/2), cos(t2/2));	t0 = 4.926234925 + 8.214985538e-7*d + t0;	t1 = 0.91744599 * sin(t0);	t0 = atan2(t1, cos(t0));	if(t0 < -pi/2)		t0 = t0 + pipi;	d = 4.88097876 + 6.30038809*d - t0;	t0 = 0.43366079 * t1;	t1 = pyth(t0);	t2 = t1*cos(slat)*cos(d-slong) - t0*sin(slat);	if(t2 > 0.46949322e-2) {		if(0.46949322e-2*t2 + 0.999988979*pyth(t2) < q)			return 0;	}	t2 = t1*cos(glat)*cos(d-wlong) - t0*sin(glat);	if(t2 < .1)		return 0;	return 1;}