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