1 cd5bae78 2004-04-21 devnull #include "astro.h"
  2 cd5bae78 2004-04-21 devnull 
  3 cd5bae78 2004-04-21 devnull char*	satlst[] =
  4 cd5bae78 2004-04-21 devnull {
  5 cbeb0b26 2006-04-01 devnull 	0
  6 cd5bae78 2004-04-21 devnull };
  7 cd5bae78 2004-04-21 devnull 
  8 cd5bae78 2004-04-21 devnull struct
  9 cd5bae78 2004-04-21 devnull {
 10 cd5bae78 2004-04-21 devnull 	double	time;
 11 cd5bae78 2004-04-21 devnull 	double	tilt;
 12 cd5bae78 2004-04-21 devnull 	double	pnni;
 13 cd5bae78 2004-04-21 devnull 	double	psi;
 14 cd5bae78 2004-04-21 devnull 	double	ppi;
 15 cd5bae78 2004-04-21 devnull 	double	d1pp;
 16 cd5bae78 2004-04-21 devnull 	double	peri;
 17 cd5bae78 2004-04-21 devnull 	double	d1per;
 18 cd5bae78 2004-04-21 devnull 	double	e0;
 19 cd5bae78 2004-04-21 devnull 	double	deo;
 20 cd5bae78 2004-04-21 devnull 	double	rdp;
 21 cd5bae78 2004-04-21 devnull 	double	st;
 22 cd5bae78 2004-04-21 devnull 	double	ct;
 23 cd5bae78 2004-04-21 devnull 	double	rot;
 24 cd5bae78 2004-04-21 devnull 	double	semi;
 25 cd5bae78 2004-04-21 devnull } satl;
 26 cd5bae78 2004-04-21 devnull 
 27 cd5bae78 2004-04-21 devnull void
 28 cd5bae78 2004-04-21 devnull satels(void)
 29 cd5bae78 2004-04-21 devnull {
 30 cd5bae78 2004-04-21 devnull 	double ifa[10], t, t1, t2, tinc;
 31 cd5bae78 2004-04-21 devnull 	char **satp;
 32 cd5bae78 2004-04-21 devnull 	int flag, f, i, n;
 33 cd5bae78 2004-04-21 devnull 
 34 cd5bae78 2004-04-21 devnull 	satp = satlst;
 35 cd5bae78 2004-04-21 devnull 
 36 cd5bae78 2004-04-21 devnull loop:
 37 cd5bae78 2004-04-21 devnull 	if(*satp == 0)
 38 cd5bae78 2004-04-21 devnull 		return;
 39 cd5bae78 2004-04-21 devnull 	f = open(*satp, 0);
 40 cd5bae78 2004-04-21 devnull 	if(f < 0) {
 41 cd5bae78 2004-04-21 devnull 		fprint(2, "cannot open %s\n", *satp);
 42 cd5bae78 2004-04-21 devnull 		satp += 2;
 43 cd5bae78 2004-04-21 devnull 		goto loop;
 44 cd5bae78 2004-04-21 devnull 	}
 45 cd5bae78 2004-04-21 devnull 	satp++;
 46 cd5bae78 2004-04-21 devnull 	rline(f);
 47 cd5bae78 2004-04-21 devnull 	tinc = atof(skip(6));
 48 cd5bae78 2004-04-21 devnull 	rline(f);
 49 cd5bae78 2004-04-21 devnull 	rline(f);
 50 cd5bae78 2004-04-21 devnull 	for(i=0; i<9; i++)
 51 cd5bae78 2004-04-21 devnull 		ifa[i] = atof(skip(i));
 52 cd5bae78 2004-04-21 devnull 	n = ifa[0];
 53 cd5bae78 2004-04-21 devnull 	i = ifa[1];
 54 cd5bae78 2004-04-21 devnull 	t = dmo[i-1] + ifa[2] - 1.;
 55 cd5bae78 2004-04-21 devnull 	if(n%4 == 0 && i > 2)
 56 cd5bae78 2004-04-21 devnull 		t += 1.;
 57 cd5bae78 2004-04-21 devnull 	for(i=1970; i<n; i++) {
 58 cd5bae78 2004-04-21 devnull 		t += 365.;
 59 cd5bae78 2004-04-21 devnull 		if(i%4 == 0)
 60 cd5bae78 2004-04-21 devnull 			t += 1.;
 61 cd5bae78 2004-04-21 devnull 	}
 62 cd5bae78 2004-04-21 devnull 	t = (t * 24. + ifa[3]) * 60. + ifa[4];
 63 cd5bae78 2004-04-21 devnull 	satl.time = t * 60.;
 64 cd5bae78 2004-04-21 devnull 	satl.tilt = ifa[5] * radian;
 65 cd5bae78 2004-04-21 devnull 	satl.pnni = ifa[6] * radian;
 66 cd5bae78 2004-04-21 devnull 	satl.psi = ifa[7];
 67 cd5bae78 2004-04-21 devnull 	satl.ppi = ifa[8] * radian;
 68 cd5bae78 2004-04-21 devnull 	rline(f);
 69 cd5bae78 2004-04-21 devnull 	for(i=0; i<5; i++)
 70 cd5bae78 2004-04-21 devnull 		ifa[i] = atof(skip(i));
 71 cd5bae78 2004-04-21 devnull 	satl.d1pp = ifa[0] * radian;
 72 cd5bae78 2004-04-21 devnull 	satl.peri = ifa[1];
 73 cd5bae78 2004-04-21 devnull 	satl.d1per = ifa[2];
 74 cd5bae78 2004-04-21 devnull 	satl.e0 = ifa[3];
 75 cd5bae78 2004-04-21 devnull 	satl.deo = 0;
 76 cd5bae78 2004-04-21 devnull 	satl.rdp = ifa[4];
 77 cd5bae78 2004-04-21 devnull 
 78 cd5bae78 2004-04-21 devnull 	satl.st = sin(satl.tilt);
 79 cd5bae78 2004-04-21 devnull 	satl.ct = cos(satl.tilt);
 80 cd5bae78 2004-04-21 devnull 	satl.rot = pipi / (1440. + satl.psi);
 81 cd5bae78 2004-04-21 devnull 	satl.semi = satl.rdp * (1. + satl.e0);
 82 cd5bae78 2004-04-21 devnull 
 83 cd5bae78 2004-04-21 devnull 	n = PER*288.; /* 5 min steps */
 84 cd5bae78 2004-04-21 devnull 	t = day;
 85 cd5bae78 2004-04-21 devnull 	for(i=0; i<n; i++) {
 86 cd5bae78 2004-04-21 devnull 		if(sunel((t-day)/deld) > 0.)
 87 cd5bae78 2004-04-21 devnull 			goto out;
 88 cd5bae78 2004-04-21 devnull 		satel(t);
 89 cd5bae78 2004-04-21 devnull 		if(el > 0) {
 90 cd5bae78 2004-04-21 devnull 			t1 = t;
 91 cd5bae78 2004-04-21 devnull 			flag = 0;
 92 cd5bae78 2004-04-21 devnull 			do {
 93 cd5bae78 2004-04-21 devnull 				if(el > 30.)
 94 cd5bae78 2004-04-21 devnull 					flag++;
 95 cd5bae78 2004-04-21 devnull 				t -= tinc/(24.*3600.);
 96 cd5bae78 2004-04-21 devnull 				satel(t);
 97 cd5bae78 2004-04-21 devnull 			} while(el > 0.);
 98 cd5bae78 2004-04-21 devnull 			t2 = (t - day)/deld;
 99 cd5bae78 2004-04-21 devnull 			t = t1;
100 cd5bae78 2004-04-21 devnull 			do {
101 cd5bae78 2004-04-21 devnull 				t += tinc/(24.*3600.);
102 cd5bae78 2004-04-21 devnull 				satel(t);
103 cd5bae78 2004-04-21 devnull 				if(el > 30.)
104 cd5bae78 2004-04-21 devnull 					flag++;
105 cd5bae78 2004-04-21 devnull 			} while(el > 0.);
106 cd5bae78 2004-04-21 devnull 			if(flag)
107 cd5bae78 2004-04-21 devnull 				if((*satp)[0] == '-')
108 cd5bae78 2004-04-21 devnull 					event("%s pass at ", (*satp)+1, "",
109 cd5bae78 2004-04-21 devnull 						t2, SIGNIF+PTIME+DARK); else
110 cd5bae78 2004-04-21 devnull 					event("%s pass at ", *satp, "",
111 cd5bae78 2004-04-21 devnull 						t2, PTIME+DARK);
112 cd5bae78 2004-04-21 devnull 		}
113 cd5bae78 2004-04-21 devnull 	out:
114 cd5bae78 2004-04-21 devnull 		t += 5./(24.*60.);
115 cd5bae78 2004-04-21 devnull 	}
116 cd5bae78 2004-04-21 devnull 	close(f);
117 cd5bae78 2004-04-21 devnull 	satp++;
118 cd5bae78 2004-04-21 devnull 	goto loop;
119 cd5bae78 2004-04-21 devnull }
120 cd5bae78 2004-04-21 devnull 
121 cd5bae78 2004-04-21 devnull void
122 cd5bae78 2004-04-21 devnull satel(double time)
123 cd5bae78 2004-04-21 devnull {
124 cd5bae78 2004-04-21 devnull 	int i;
125 cd5bae78 2004-04-21 devnull 	double amean, an, coc, csl, d, de, enom, eo;
126 cd5bae78 2004-04-21 devnull 	double pp, q, rdp, slong, ssl, t, tp;
127 cd5bae78 2004-04-21 devnull 
128 cd5bae78 2004-04-21 devnull 	i = 500;
129 cd5bae78 2004-04-21 devnull 	el = -1;
130 cd5bae78 2004-04-21 devnull 	time = (time-25567.5) * 86400;
131 cd5bae78 2004-04-21 devnull 	t = (time - satl.time)/60;
132 cd5bae78 2004-04-21 devnull 	if(t < 0)
133 cd5bae78 2004-04-21 devnull 		return; /* too early for satelites */
134 cd5bae78 2004-04-21 devnull 	an = floor(t/satl.peri);
135 cd5bae78 2004-04-21 devnull 	while(an*satl.peri + an*an*satl.d1per/2. <= t) {
136 cd5bae78 2004-04-21 devnull 		an += 1;
137 cd5bae78 2004-04-21 devnull 		if(--i == 0)
138 cd5bae78 2004-04-21 devnull 			return;
139 cd5bae78 2004-04-21 devnull 	}
140 cd5bae78 2004-04-21 devnull 	while((tp = an*satl.peri + an*an*satl.d1per/2.) > t) {
141 cd5bae78 2004-04-21 devnull 		an -= 1;
142 cd5bae78 2004-04-21 devnull 		if(--i == 0)
143 cd5bae78 2004-04-21 devnull 			return;
144 cd5bae78 2004-04-21 devnull 	}
145 cd5bae78 2004-04-21 devnull 	amean = (t-tp)/(satl.peri+(an+.5)*satl.d1per);
146 cd5bae78 2004-04-21 devnull 	pp = satl.ppi+(an+amean)*satl.d1pp;
147 cd5bae78 2004-04-21 devnull 	amean *= pipi;
148 cd5bae78 2004-04-21 devnull 	eo = satl.e0+satl.deo*an;
149 cd5bae78 2004-04-21 devnull 	rdp = satl.semi/(1+eo);
150 cd5bae78 2004-04-21 devnull 	enom = amean+eo*sin(amean);
151 cd5bae78 2004-04-21 devnull 	do {
152 cd5bae78 2004-04-21 devnull 		de = (amean-enom+eo*sin(enom))/(1.0-eo*cos(enom));
153 cd5bae78 2004-04-21 devnull 		enom += de;
154 cd5bae78 2004-04-21 devnull 		if(--i == 0)
155 cd5bae78 2004-04-21 devnull 			return;
156 cd5bae78 2004-04-21 devnull 	} while(fabs(de) >= 1.0e-6);
157 cd5bae78 2004-04-21 devnull 	q = 3963.35*erad/(rdp*(1-eo*cos(enom))/(1-eo));
158 cd5bae78 2004-04-21 devnull 	d = pp + 2*atan2(sqrt((1+eo)/(1-eo))*sin(enom/2),cos(enom/2));
159 cd5bae78 2004-04-21 devnull 	slong = satl.pnni + t*satl.rot -
160 cd5bae78 2004-04-21 devnull 		atan2(satl.ct*sin(d), cos(d));
161 cd5bae78 2004-04-21 devnull 	ssl = satl.st*sin(d);
162 cd5bae78 2004-04-21 devnull 	csl = pyth(ssl);
163 cd5bae78 2004-04-21 devnull 	if(vis(time, atan2(ssl,csl), slong, q)) {
164 cd5bae78 2004-04-21 devnull 		coc = ssl*sin(glat) + csl*cos(glat)*cos(wlong-slong);
165 cd5bae78 2004-04-21 devnull 		el = atan2(coc-q, pyth(coc));
166 cd5bae78 2004-04-21 devnull 		el /= radian;
167 cd5bae78 2004-04-21 devnull 	}
168 cd5bae78 2004-04-21 devnull }
169 cd5bae78 2004-04-21 devnull 
170 cd5bae78 2004-04-21 devnull int
171 cd5bae78 2004-04-21 devnull vis(double t, double slat, double slong, double q)
172 cd5bae78 2004-04-21 devnull {
173 cd5bae78 2004-04-21 devnull 	double t0, t1, t2, d;
174 cd5bae78 2004-04-21 devnull 
175 cd5bae78 2004-04-21 devnull 	d = t/86400 - .005375 + 3653;
176 cd5bae78 2004-04-21 devnull 	t0 = 6.238030674 + .01720196977*d;
177 cd5bae78 2004-04-21 devnull 	t2 = t0 + .0167253303*sin(t0);
178 cd5bae78 2004-04-21 devnull 	do {
179 cd5bae78 2004-04-21 devnull 		t1 = (t0 - t2 + .0167259152*sin(t2)) /
180 cd5bae78 2004-04-21 devnull 			(1 - .0167259152*cos(t2));
181 cd5bae78 2004-04-21 devnull 		t2 = t2 + t1;
182 cd5bae78 2004-04-21 devnull 	} while(fabs(t1) >= 1.e-4);
183 cd5bae78 2004-04-21 devnull 	t0 = 2*atan2(1.01686816*sin(t2/2), cos(t2/2));
184 cd5bae78 2004-04-21 devnull 	t0 = 4.926234925 + 8.214985538e-7*d + t0;
185 cd5bae78 2004-04-21 devnull 	t1 = 0.91744599 * sin(t0);
186 cd5bae78 2004-04-21 devnull 	t0 = atan2(t1, cos(t0));
187 cd5bae78 2004-04-21 devnull 	if(t0 < -pi/2)
188 cd5bae78 2004-04-21 devnull 		t0 = t0 + pipi;
189 cd5bae78 2004-04-21 devnull 	d = 4.88097876 + 6.30038809*d - t0;
190 cd5bae78 2004-04-21 devnull 	t0 = 0.43366079 * t1;
191 cd5bae78 2004-04-21 devnull 	t1 = pyth(t0);
192 cd5bae78 2004-04-21 devnull 	t2 = t1*cos(slat)*cos(d-slong) - t0*sin(slat);
193 cd5bae78 2004-04-21 devnull 	if(t2 > 0.46949322e-2) {
194 cd5bae78 2004-04-21 devnull 		if(0.46949322e-2*t2 + 0.999988979*pyth(t2) < q)
195 cd5bae78 2004-04-21 devnull 			return 0;
196 cd5bae78 2004-04-21 devnull 	}
197 cd5bae78 2004-04-21 devnull 	t2 = t1*cos(glat)*cos(d-wlong) - t0*sin(glat);
198 cd5bae78 2004-04-21 devnull 	if(t2 < .1)
199 cd5bae78 2004-04-21 devnull 		return 0;
200 cd5bae78 2004-04-21 devnull 	return 1;
201 cd5bae78 2004-04-21 devnull }