1 8a3b2ceb 2004-04-24 devnull #include <u.h>
2 8a3b2ceb 2004-04-24 devnull #include <libc.h>
3 8a3b2ceb 2004-04-24 devnull #include <bio.h>
4 8a3b2ceb 2004-04-24 devnull #include "sky.h"
6 8a3b2ceb 2004-04-24 devnull double PI_180 = 0.0174532925199432957692369;
7 8a3b2ceb 2004-04-24 devnull double TWOPI = 6.2831853071795864769252867665590057683943387987502;
8 8a3b2ceb 2004-04-24 devnull double LN2 = 0.69314718055994530941723212145817656807550013436025;
9 8a3b2ceb 2004-04-24 devnull static double angledangle=(180./PI)*MILLIARCSEC;
11 8a3b2ceb 2004-04-24 devnull #define rint scatrint
14 8a3b2ceb 2004-04-24 devnull rint(char *p, int n)
18 8a3b2ceb 2004-04-24 devnull while(*p==' ' && n)
19 8a3b2ceb 2004-04-24 devnull p++, --n;
20 8a3b2ceb 2004-04-24 devnull while(n--)
21 8a3b2ceb 2004-04-24 devnull i=i*10+*p++-'0';
22 8a3b2ceb 2004-04-24 devnull return i;
26 8a3b2ceb 2004-04-24 devnull dangle(Angle angle)
28 8a3b2ceb 2004-04-24 devnull return angle*angledangle;
32 8a3b2ceb 2004-04-24 devnull angle(DAngle dangle)
34 8a3b2ceb 2004-04-24 devnull return dangle/angledangle;
38 8a3b2ceb 2004-04-24 devnull rfloat(char *p, int n)
40 8a3b2ceb 2004-04-24 devnull double i, d=0;
42 8a3b2ceb 2004-04-24 devnull while(*p==' ' && n)
43 8a3b2ceb 2004-04-24 devnull p++, --n;
44 8a3b2ceb 2004-04-24 devnull if(*p == '+')
45 8a3b2ceb 2004-04-24 devnull return rfloat(p+1, n-1);
46 8a3b2ceb 2004-04-24 devnull if(*p == '-')
47 8a3b2ceb 2004-04-24 devnull return -rfloat(p+1, n-1);
48 8a3b2ceb 2004-04-24 devnull while(*p == ' ' && n)
49 8a3b2ceb 2004-04-24 devnull p++, --n;
50 8a3b2ceb 2004-04-24 devnull if(n == 0)
51 8a3b2ceb 2004-04-24 devnull return 0.0;
52 8a3b2ceb 2004-04-24 devnull while(n-- && *p!='.')
53 8a3b2ceb 2004-04-24 devnull d = d*10+*p++-'0';
54 8a3b2ceb 2004-04-24 devnull if(n <= 0)
55 8a3b2ceb 2004-04-24 devnull return d;
58 8a3b2ceb 2004-04-24 devnull while(n--)
59 8a3b2ceb 2004-04-24 devnull d+=(*p++-'0')/(i*=10.);
60 8a3b2ceb 2004-04-24 devnull return d;
64 8a3b2ceb 2004-04-24 devnull sign(int c)
66 8a3b2ceb 2004-04-24 devnull if(c=='-')
67 8a3b2ceb 2004-04-24 devnull return -1;
68 8a3b2ceb 2004-04-24 devnull return 1;
72 8a3b2ceb 2004-04-24 devnull hms(Angle a)
74 8a3b2ceb 2004-04-24 devnull static char buf[20];
75 8a3b2ceb 2004-04-24 devnull double x;
76 8a3b2ceb 2004-04-24 devnull int h, m, s, ts;
78 8a3b2ceb 2004-04-24 devnull x=DEG(a)/15;
79 8a3b2ceb 2004-04-24 devnull x += 0.5/36000.; /* round up half of 0.1 sec */
80 8a3b2ceb 2004-04-24 devnull h = floor(x);
83 8a3b2ceb 2004-04-24 devnull m = floor(x);
86 8a3b2ceb 2004-04-24 devnull s = floor(x);
88 8a3b2ceb 2004-04-24 devnull ts = 10*x;
89 8a3b2ceb 2004-04-24 devnull sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts);
90 8a3b2ceb 2004-04-24 devnull return buf;
94 8a3b2ceb 2004-04-24 devnull dms(Angle a)
96 8a3b2ceb 2004-04-24 devnull static char buf[20];
97 8a3b2ceb 2004-04-24 devnull double x;
98 8a3b2ceb 2004-04-24 devnull int sign, d, m, s, ts;
100 8a3b2ceb 2004-04-24 devnull x = DEG(a);
101 8a3b2ceb 2004-04-24 devnull sign='+';
102 8a3b2ceb 2004-04-24 devnull if(a<0){
103 8a3b2ceb 2004-04-24 devnull sign='-';
106 8a3b2ceb 2004-04-24 devnull x += 0.5/36000.; /* round up half of 0.1 arcsecond */
107 8a3b2ceb 2004-04-24 devnull d = floor(x);
109 8a3b2ceb 2004-04-24 devnull x *= 60;
110 8a3b2ceb 2004-04-24 devnull m = floor(x);
112 8a3b2ceb 2004-04-24 devnull x *= 60;
113 8a3b2ceb 2004-04-24 devnull s = floor(x);
115 8a3b2ceb 2004-04-24 devnull ts = floor(10*x);
116 8a3b2ceb 2004-04-24 devnull sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts);
117 8a3b2ceb 2004-04-24 devnull return buf;
121 8a3b2ceb 2004-04-24 devnull ms(Angle a)
123 8a3b2ceb 2004-04-24 devnull static char buf[20];
124 8a3b2ceb 2004-04-24 devnull double x;
125 8a3b2ceb 2004-04-24 devnull int d, m, s, ts;
127 8a3b2ceb 2004-04-24 devnull x = DEG(a);
128 8a3b2ceb 2004-04-24 devnull x += 0.5/36000.; /* round up half of 0.1 arcsecond */
129 8a3b2ceb 2004-04-24 devnull d = floor(x);
131 8a3b2ceb 2004-04-24 devnull x *= 60;
132 8a3b2ceb 2004-04-24 devnull m = floor(x);
134 8a3b2ceb 2004-04-24 devnull x *= 60;
135 8a3b2ceb 2004-04-24 devnull s = floor(x);
137 8a3b2ceb 2004-04-24 devnull ts = floor(10*x);
138 8a3b2ceb 2004-04-24 devnull if(d != 0)
139 8a3b2ceb 2004-04-24 devnull sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);
141 8a3b2ceb 2004-04-24 devnull sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
142 8a3b2ceb 2004-04-24 devnull return buf;
146 8a3b2ceb 2004-04-24 devnull hm(Angle a)
148 8a3b2ceb 2004-04-24 devnull static char buf[20];
149 8a3b2ceb 2004-04-24 devnull double x;
150 8a3b2ceb 2004-04-24 devnull int h, m, n;
152 8a3b2ceb 2004-04-24 devnull x = DEG(a)/15;
153 8a3b2ceb 2004-04-24 devnull x += 0.5/600.; /* round up half of tenth of minute */
154 8a3b2ceb 2004-04-24 devnull h = floor(x);
156 8a3b2ceb 2004-04-24 devnull x *= 60;
157 8a3b2ceb 2004-04-24 devnull m = floor(x);
159 8a3b2ceb 2004-04-24 devnull x *= 10;
160 8a3b2ceb 2004-04-24 devnull n = floor(x);
161 8a3b2ceb 2004-04-24 devnull sprint(buf, "%dh%.2d.%1dm", h, m, n);
162 8a3b2ceb 2004-04-24 devnull return buf;
166 8a3b2ceb 2004-04-24 devnull hm5(Angle a)
168 8a3b2ceb 2004-04-24 devnull static char buf[20];
169 8a3b2ceb 2004-04-24 devnull double x;
170 8a3b2ceb 2004-04-24 devnull int h, m;
172 8a3b2ceb 2004-04-24 devnull x = DEG(a)/15;
173 8a3b2ceb 2004-04-24 devnull x += 2.5/60.; /* round up 2.5m */
174 8a3b2ceb 2004-04-24 devnull h = floor(x);
176 8a3b2ceb 2004-04-24 devnull x *= 60;
177 8a3b2ceb 2004-04-24 devnull m = floor(x);
178 8a3b2ceb 2004-04-24 devnull m -= m % 5;
179 8a3b2ceb 2004-04-24 devnull sprint(buf, "%dh%.2dm", h, m);
180 8a3b2ceb 2004-04-24 devnull return buf;
184 8a3b2ceb 2004-04-24 devnull dm(Angle a)
186 8a3b2ceb 2004-04-24 devnull static char buf[20];
187 8a3b2ceb 2004-04-24 devnull double x;
188 8a3b2ceb 2004-04-24 devnull int sign, d, m, n;
190 8a3b2ceb 2004-04-24 devnull x = DEG(a);
191 8a3b2ceb 2004-04-24 devnull sign='+';
192 8a3b2ceb 2004-04-24 devnull if(a<0){
193 8a3b2ceb 2004-04-24 devnull sign='-';
196 8a3b2ceb 2004-04-24 devnull x += 0.5/600.; /* round up half of tenth of arcminute */
197 8a3b2ceb 2004-04-24 devnull d = floor(x);
199 8a3b2ceb 2004-04-24 devnull x *= 60;
200 8a3b2ceb 2004-04-24 devnull m = floor(x);
202 8a3b2ceb 2004-04-24 devnull x *= 10;
203 8a3b2ceb 2004-04-24 devnull n = floor(x);
204 8a3b2ceb 2004-04-24 devnull sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);
205 8a3b2ceb 2004-04-24 devnull return buf;
209 8a3b2ceb 2004-04-24 devnull deg(Angle a)
211 8a3b2ceb 2004-04-24 devnull static char buf[20];
212 8a3b2ceb 2004-04-24 devnull double x;
213 8a3b2ceb 2004-04-24 devnull int sign, d;
215 8a3b2ceb 2004-04-24 devnull x = DEG(a);
216 8a3b2ceb 2004-04-24 devnull sign='+';
217 8a3b2ceb 2004-04-24 devnull if(a<0){
218 8a3b2ceb 2004-04-24 devnull sign='-';
221 8a3b2ceb 2004-04-24 devnull x += 0.5; /* round up half degree */
222 8a3b2ceb 2004-04-24 devnull d = floor(x);
223 8a3b2ceb 2004-04-24 devnull sprint(buf, "%c%d°", sign, d);
224 8a3b2ceb 2004-04-24 devnull return buf;
228 8a3b2ceb 2004-04-24 devnull getword(char *ou, char *in)
232 8a3b2ceb 2004-04-24 devnull for(;;) {
233 8a3b2ceb 2004-04-24 devnull c = *in++;
234 8a3b2ceb 2004-04-24 devnull if(c == ' ' || c == '\t')
235 8a3b2ceb 2004-04-24 devnull continue;
236 8a3b2ceb 2004-04-24 devnull if(c == 0)
237 8a3b2ceb 2004-04-24 devnull return 0;
241 8a3b2ceb 2004-04-24 devnull if(c == '\'')
242 8a3b2ceb 2004-04-24 devnull for(;;) {
243 8a3b2ceb 2004-04-24 devnull if(c >= 'A' && c <= 'Z')
244 8a3b2ceb 2004-04-24 devnull c += 'a' - 'A';
245 8a3b2ceb 2004-04-24 devnull *ou++ = c;
246 8a3b2ceb 2004-04-24 devnull c = *in++;
247 8a3b2ceb 2004-04-24 devnull if(c == 0)
248 8a3b2ceb 2004-04-24 devnull return 0;
249 8a3b2ceb 2004-04-24 devnull if(c == '\'') {
250 8a3b2ceb 2004-04-24 devnull *ou = 0;
251 8a3b2ceb 2004-04-24 devnull return in-1;
254 8a3b2ceb 2004-04-24 devnull for(;;) {
255 8a3b2ceb 2004-04-24 devnull if(c >= 'A' && c <= 'Z')
256 8a3b2ceb 2004-04-24 devnull c += 'a' - 'A';
257 8a3b2ceb 2004-04-24 devnull *ou++ = c;
258 8a3b2ceb 2004-04-24 devnull c = *in++;
259 8a3b2ceb 2004-04-24 devnull if(c == ' ' || c == '\t' || c == 0) {
260 8a3b2ceb 2004-04-24 devnull *ou = 0;
261 8a3b2ceb 2004-04-24 devnull return in-1;
267 8a3b2ceb 2004-04-24 devnull * Read formatted angle. Must contain no embedded blanks
270 8a3b2ceb 2004-04-24 devnull getra(char *p)
273 8a3b2ceb 2004-04-24 devnull char *q;
274 8a3b2ceb 2004-04-24 devnull Angle f, d;
275 8a3b2ceb 2004-04-24 devnull int neg;
277 8a3b2ceb 2004-04-24 devnull neg = 0;
279 8a3b2ceb 2004-04-24 devnull while(*p == ' ')
281 8a3b2ceb 2004-04-24 devnull for(;;) {
282 8a3b2ceb 2004-04-24 devnull if(*p == ' ' || *p=='\0')
283 8a3b2ceb 2004-04-24 devnull goto Return;
284 8a3b2ceb 2004-04-24 devnull if(*p == '-') {
285 8a3b2ceb 2004-04-24 devnull neg = 1;
288 8a3b2ceb 2004-04-24 devnull if(*p == '+') {
289 8a3b2ceb 2004-04-24 devnull neg = 0;
293 8a3b2ceb 2004-04-24 devnull f = strtod(p, &q);
294 8a3b2ceb 2004-04-24 devnull if(q > p) {
297 8a3b2ceb 2004-04-24 devnull p += chartorune(&r, p);
298 8a3b2ceb 2004-04-24 devnull switch(r) {
299 8a3b2ceb 2004-04-24 devnull default:
303 8a3b2ceb 2004-04-24 devnull return RAD(d);
304 8a3b2ceb 2004-04-24 devnull case 'h':
305 8a3b2ceb 2004-04-24 devnull d += f*15;
307 8a3b2ceb 2004-04-24 devnull case 'm':
308 8a3b2ceb 2004-04-24 devnull d += f/4;
310 8a3b2ceb 2004-04-24 devnull case 's':
311 8a3b2ceb 2004-04-24 devnull d += f/240;
313 8a3b2ceb 2004-04-24 devnull case 0xB0: /* ° */
316 8a3b2ceb 2004-04-24 devnull case '\'':
317 8a3b2ceb 2004-04-24 devnull d += f/60;
319 8a3b2ceb 2004-04-24 devnull case '\"':
320 8a3b2ceb 2004-04-24 devnull d += f/3600;
324 8a3b2ceb 2004-04-24 devnull return 0;
328 8a3b2ceb 2004-04-24 devnull xsqrt(double a)
331 8a3b2ceb 2004-04-24 devnull if(a < 0)
332 8a3b2ceb 2004-04-24 devnull return 0;
333 8a3b2ceb 2004-04-24 devnull return sqrt(a);
337 8a3b2ceb 2004-04-24 devnull dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
339 8a3b2ceb 2004-04-24 devnull double a;
341 8a3b2ceb 2004-04-24 devnull a = sin(dec1) * sin(dec2) +
342 8a3b2ceb 2004-04-24 devnull cos(dec1) * cos(dec2) *
343 8a3b2ceb 2004-04-24 devnull cos(ra1 - ra2);
344 8a3b2ceb 2004-04-24 devnull a = atan2(xsqrt(1 - a*a), a);
345 8a3b2ceb 2004-04-24 devnull if(a < 0)
347 8a3b2ceb 2004-04-24 devnull return a;
351 8a3b2ceb 2004-04-24 devnull dogamma(Pix c)
353 8a3b2ceb 2004-04-24 devnull float f;
355 8a3b2ceb 2004-04-24 devnull f = c - gam.min;
356 8a3b2ceb 2004-04-24 devnull if(f < 1)
359 8a3b2ceb 2004-04-24 devnull if(gam.absgamma == 1)
360 8a3b2ceb 2004-04-24 devnull c = f * gam.mult2;
362 8a3b2ceb 2004-04-24 devnull c = exp(log(f*gam.mult1) * gam.absgamma) * 255;
363 8a3b2ceb 2004-04-24 devnull if(c > 255)
364 8a3b2ceb 2004-04-24 devnull c = 255;
365 8a3b2ceb 2004-04-24 devnull if(gam.neg)
366 8a3b2ceb 2004-04-24 devnull c = 255-c;
367 8a3b2ceb 2004-04-24 devnull return c;