Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

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