1 #include <u.h>
2 #include <libc.h>
3 #include <bio.h>
4 #include "sky.h"
5
6 double PI_180 = 0.0174532925199432957692369;
7 double TWOPI = 6.2831853071795864769252867665590057683943387987502;
8 double LN2 = 0.69314718055994530941723212145817656807550013436025;
9 static double angledangle=(180./PI)*MILLIARCSEC;
10
11 int
rint(char * p,int n)12 rint(char *p, int n)
13 {
14 int i=0;
15
16 while(*p==' ' && n)
17 p++, --n;
18 while(n--)
19 i=i*10+*p++-'0';
20 return i;
21 }
22
23 DAngle
dangle(Angle angle)24 dangle(Angle angle)
25 {
26 return angle*angledangle;
27 }
28
29 Angle
angle(DAngle dangle)30 angle(DAngle dangle)
31 {
32 return dangle/angledangle;
33 }
34
35 double
rfloat(char * p,int n)36 rfloat(char *p, int n)
37 {
38 double i, d=0;
39
40 while(*p==' ' && n)
41 p++, --n;
42 if(*p == '+')
43 return rfloat(p+1, n-1);
44 if(*p == '-')
45 return -rfloat(p+1, n-1);
46 while(*p == ' ' && n)
47 p++, --n;
48 if(n == 0)
49 return 0.0;
50 while(n-- && *p!='.')
51 d = d*10+*p++-'0';
52 if(n <= 0)
53 return d;
54 p++;
55 i = 1;
56 while(n--)
57 d+=(*p++-'0')/(i*=10.);
58 return d;
59 }
60
61 int
sign(int c)62 sign(int c)
63 {
64 if(c=='-')
65 return -1;
66 return 1;
67 }
68
69 char*
hms(Angle a)70 hms(Angle a)
71 {
72 static char buf[20];
73 double x;
74 int h, m, s, ts;
75
76 x=DEG(a)/15;
77 x += 0.5/36000.; /* round up half of 0.1 sec */
78 h = floor(x);
79 x -= h;
80 x *= 60;
81 m = floor(x);
82 x -= m;
83 x *= 60;
84 s = floor(x);
85 x -= s;
86 ts = 10*x;
87 sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts);
88 return buf;
89 }
90
91 char*
dms(Angle a)92 dms(Angle a)
93 {
94 static char buf[20];
95 double x;
96 int sign, d, m, s, ts;
97
98 x = DEG(a);
99 sign='+';
100 if(a<0){
101 sign='-';
102 x=-x;
103 }
104 x += 0.5/36000.; /* round up half of 0.1 arcsecond */
105 d = floor(x);
106 x -= d;
107 x *= 60;
108 m = floor(x);
109 x -= m;
110 x *= 60;
111 s = floor(x);
112 x -= s;
113 ts = floor(10*x);
114 sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts);
115 return buf;
116 }
117
118 char*
ms(Angle a)119 ms(Angle a)
120 {
121 static char buf[20];
122 double x;
123 int d, m, s, ts;
124
125 x = DEG(a);
126 x += 0.5/36000.; /* round up half of 0.1 arcsecond */
127 d = floor(x);
128 x -= d;
129 x *= 60;
130 m = floor(x);
131 x -= m;
132 x *= 60;
133 s = floor(x);
134 x -= s;
135 ts = floor(10*x);
136 if(d != 0)
137 sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);
138 else
139 sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
140 return buf;
141 }
142
143 char*
hm(Angle a)144 hm(Angle a)
145 {
146 static char buf[20];
147 double x;
148 int h, m, n;
149
150 x = DEG(a)/15;
151 x += 0.5/600.; /* round up half of tenth of minute */
152 h = floor(x);
153 x -= h;
154 x *= 60;
155 m = floor(x);
156 x -= m;
157 x *= 10;
158 n = floor(x);
159 sprint(buf, "%dh%.2d.%1dm", h, m, n);
160 return buf;
161 }
162
163 char*
hm5(Angle a)164 hm5(Angle a)
165 {
166 static char buf[20];
167 double x;
168 int h, m;
169
170 x = DEG(a)/15;
171 x += 2.5/60.; /* round up 2.5m */
172 h = floor(x);
173 x -= h;
174 x *= 60;
175 m = floor(x);
176 m -= m % 5;
177 sprint(buf, "%dh%.2dm", h, m);
178 return buf;
179 }
180
181 char*
dm(Angle a)182 dm(Angle a)
183 {
184 static char buf[20];
185 double x;
186 int sign, d, m, n;
187
188 x = DEG(a);
189 sign='+';
190 if(a<0){
191 sign='-';
192 x=-x;
193 }
194 x += 0.5/600.; /* round up half of tenth of arcminute */
195 d = floor(x);
196 x -= d;
197 x *= 60;
198 m = floor(x);
199 x -= m;
200 x *= 10;
201 n = floor(x);
202 sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);
203 return buf;
204 }
205
206 char*
deg(Angle a)207 deg(Angle a)
208 {
209 static char buf[20];
210 double x;
211 int sign, d;
212
213 x = DEG(a);
214 sign='+';
215 if(a<0){
216 sign='-';
217 x=-x;
218 }
219 x += 0.5; /* round up half degree */
220 d = floor(x);
221 sprint(buf, "%c%d°", sign, d);
222 return buf;
223 }
224
225 char*
getword(char * ou,char * in)226 getword(char *ou, char *in)
227 {
228 int c;
229
230 for(;;) {
231 c = *in++;
232 if(c == ' ' || c == '\t')
233 continue;
234 if(c == 0)
235 return 0;
236 break;
237 }
238
239 if(c == '\'')
240 for(;;) {
241 if(c >= 'A' && c <= 'Z')
242 c += 'a' - 'A';
243 *ou++ = c;
244 c = *in++;
245 if(c == 0)
246 return 0;
247 if(c == '\'') {
248 *ou = 0;
249 return in-1;
250 }
251 }
252 for(;;) {
253 if(c >= 'A' && c <= 'Z')
254 c += 'a' - 'A';
255 *ou++ = c;
256 c = *in++;
257 if(c == ' ' || c == '\t' || c == 0) {
258 *ou = 0;
259 return in-1;
260 }
261 }
262 }
263
264 /*
265 * Read formatted angle. Must contain no embedded blanks
266 */
267 Angle
getra(char * p)268 getra(char *p)
269 {
270 Rune r;
271 char *q;
272 Angle f, d;
273 int neg;
274
275 neg = 0;
276 d = 0;
277 while(*p == ' ')
278 p++;
279 for(;;) {
280 if(*p == ' ' || *p=='\0')
281 goto Return;
282 if(*p == '-') {
283 neg = 1;
284 p++;
285 }
286 if(*p == '+') {
287 neg = 0;
288 p++;
289 }
290 q = p;
291 f = strtod(p, &q);
292 if(q > p) {
293 p = q;
294 }
295 p += chartorune(&r, p);
296 switch(r) {
297 default:
298 Return:
299 if(neg)
300 d = -d;
301 return RAD(d);
302 case 'h':
303 d += f*15;
304 break;
305 case 'm':
306 d += f/4;
307 break;
308 case 's':
309 d += f/240;
310 break;
311 case L'°':
312 d += f;
313 break;
314 case '\'':
315 d += f/60;
316 break;
317 case '\"':
318 d += f/3600;
319 break;
320 }
321 }
322 }
323
324 double
xsqrt(double a)325 xsqrt(double a)
326 {
327
328 if(a < 0)
329 return 0;
330 return sqrt(a);
331 }
332
333 Angle
dist(Angle ra1,Angle dec1,Angle ra2,Angle dec2)334 dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
335 {
336 double a;
337
338 a = sin(dec1) * sin(dec2) +
339 cos(dec1) * cos(dec2) *
340 cos(ra1 - ra2);
341 a = atan2(xsqrt(1 - a*a), a);
342 if(a < 0)
343 a = -a;
344 return a;
345 }
346
347 int
dogamma(Pix c)348 dogamma(Pix c)
349 {
350 float f;
351
352 f = c - gam.min;
353 if(f < 1)
354 f = 1;
355
356 if(gam.absgamma == 1)
357 c = f * gam.mult2;
358 else
359 c = exp(log(f*gam.mult1) * gam.absgamma) * 255;
360 if(c > 255)
361 c = 255;
362 if(gam.neg)
363 c = 255-c;
364 return c;
365 }
366