xref: /plan9/sys/src/cmd/scat/util.c (revision 219b2ee8daee37f4aad58d63f21287faa8e4ffdc)
13e12c5d1SDavid du Colombier #include <u.h>
23e12c5d1SDavid du Colombier #include <libc.h>
3*219b2ee8SDavid du Colombier #include <bio.h>
43e12c5d1SDavid du Colombier #include "sky.h"
53e12c5d1SDavid du Colombier 
63e12c5d1SDavid du Colombier double	PI_180	= 0.0174532925199432957692369;
73e12c5d1SDavid du Colombier double	TWOPI	= 6.2831853071795864769252867665590057683943387987502;
8*219b2ee8SDavid du Colombier double	LN2	= 0.69314718055994530941723212145817656807550013436025;
9*219b2ee8SDavid du Colombier static double angledangle=(180./PI)*60.*60.*1000.;
103e12c5d1SDavid du Colombier 
113e12c5d1SDavid du Colombier int
123e12c5d1SDavid du Colombier rint(char *p, int n)
133e12c5d1SDavid du Colombier {
143e12c5d1SDavid du Colombier 	int i=0;
153e12c5d1SDavid du Colombier 
163e12c5d1SDavid du Colombier 	while(*p==' ' && n)
173e12c5d1SDavid du Colombier 		p++, --n;
183e12c5d1SDavid du Colombier 	while(n--)
193e12c5d1SDavid du Colombier 		i=i*10+*p++-'0';
203e12c5d1SDavid du Colombier 	return i;
213e12c5d1SDavid du Colombier }
223e12c5d1SDavid du Colombier 
233e12c5d1SDavid du Colombier DAngle
243e12c5d1SDavid du Colombier dangle(Angle angle)
253e12c5d1SDavid du Colombier {
263e12c5d1SDavid du Colombier 	return angle*angledangle;
273e12c5d1SDavid du Colombier }
283e12c5d1SDavid du Colombier 
293e12c5d1SDavid du Colombier Angle
303e12c5d1SDavid du Colombier angle(DAngle dangle)
313e12c5d1SDavid du Colombier {
323e12c5d1SDavid du Colombier 	return dangle/angledangle;
333e12c5d1SDavid du Colombier }
343e12c5d1SDavid du Colombier 
353e12c5d1SDavid du Colombier double
363e12c5d1SDavid du Colombier rfloat(char *p, int n)
373e12c5d1SDavid du Colombier {
383e12c5d1SDavid du Colombier 	double i, d=0;
393e12c5d1SDavid du Colombier 
403e12c5d1SDavid du Colombier 	while(*p==' ' && n)
413e12c5d1SDavid du Colombier 		p++, --n;
423e12c5d1SDavid du Colombier 	if(*p == '+')
433e12c5d1SDavid du Colombier 		return rfloat(p+1, n-1);
443e12c5d1SDavid du Colombier 	if(*p == '-')
453e12c5d1SDavid du Colombier 		return -rfloat(p+1, n-1);
463e12c5d1SDavid du Colombier 	while(*p == ' ' && n)
473e12c5d1SDavid du Colombier 		p++, --n;
483e12c5d1SDavid du Colombier 	if(n == 0)
493e12c5d1SDavid du Colombier 		return 0.0;
503e12c5d1SDavid du Colombier 	while(n-- && *p!='.')
513e12c5d1SDavid du Colombier 		d = d*10+*p++-'0';
523e12c5d1SDavid du Colombier 	if(n <= 0)
533e12c5d1SDavid du Colombier 		return d;
543e12c5d1SDavid du Colombier 	p++;
553e12c5d1SDavid du Colombier 	i = 1;
563e12c5d1SDavid du Colombier 	while(n--)
573e12c5d1SDavid du Colombier 		d+=(*p++-'0')/(i*=10.);
583e12c5d1SDavid du Colombier 	return d;
593e12c5d1SDavid du Colombier }
603e12c5d1SDavid du Colombier 
613e12c5d1SDavid du Colombier int
623e12c5d1SDavid du Colombier sign(int c)
633e12c5d1SDavid du Colombier {
643e12c5d1SDavid du Colombier 	if(c=='-')
653e12c5d1SDavid du Colombier 		return -1;
663e12c5d1SDavid du Colombier 	return 1;
673e12c5d1SDavid du Colombier }
683e12c5d1SDavid du Colombier 
693e12c5d1SDavid du Colombier char*
703e12c5d1SDavid du Colombier hms(Angle a)
713e12c5d1SDavid du Colombier {
723e12c5d1SDavid du Colombier 	static char buf[20];
73*219b2ee8SDavid du Colombier 	double x;
743e12c5d1SDavid du Colombier 	int h, m, s, ts;
753e12c5d1SDavid du Colombier 
76*219b2ee8SDavid du Colombier 	x=DEG(a)/15;
77*219b2ee8SDavid du Colombier 	x += 0.5/36000.;	/* round up half of 0.1 sec */
783e12c5d1SDavid du Colombier 	h = floor(x);
793e12c5d1SDavid du Colombier 	x -= h;
803e12c5d1SDavid du Colombier 	x *= 60;
813e12c5d1SDavid du Colombier 	m = floor(x);
823e12c5d1SDavid du Colombier 	x -= m;
833e12c5d1SDavid du Colombier 	x *= 60;
843e12c5d1SDavid du Colombier 	s = floor(x);
853e12c5d1SDavid du Colombier 	x -= s;
863e12c5d1SDavid du Colombier 	ts = 10*x;
873e12c5d1SDavid du Colombier 	sprint(buf, "%dh%.2dm%.2d.%ds", h, m, s, ts);
883e12c5d1SDavid du Colombier 	return buf;
893e12c5d1SDavid du Colombier }
903e12c5d1SDavid du Colombier 
913e12c5d1SDavid du Colombier char*
923e12c5d1SDavid du Colombier dms(Angle a)
933e12c5d1SDavid du Colombier {
943e12c5d1SDavid du Colombier 	static char buf[20];
95*219b2ee8SDavid du Colombier 	double x;
963e12c5d1SDavid du Colombier 	int sign, d, m, s, ts;
973e12c5d1SDavid du Colombier 
98*219b2ee8SDavid du Colombier 	x = DEG(a);
993e12c5d1SDavid du Colombier 	sign='+';
1003e12c5d1SDavid du Colombier 	if(a<0){
1013e12c5d1SDavid du Colombier 		sign='-';
1023e12c5d1SDavid du Colombier 		x=-x;
1033e12c5d1SDavid du Colombier 	}
104*219b2ee8SDavid du Colombier 	x += 0.5/36000.;	/* round up half of 0.1 arcsecond */
1053e12c5d1SDavid du Colombier 	d = floor(x);
1063e12c5d1SDavid du Colombier 	x -= d;
1073e12c5d1SDavid du Colombier 	x *= 60;
1083e12c5d1SDavid du Colombier 	m = floor(x);
1093e12c5d1SDavid du Colombier 	x -= m;
1103e12c5d1SDavid du Colombier 	x *= 60;
1113e12c5d1SDavid du Colombier 	s = floor(x);
1123e12c5d1SDavid du Colombier 	x -= s;
1133e12c5d1SDavid du Colombier 	ts = floor(10*x);
1143e12c5d1SDavid du Colombier 	sprint(buf, "%c%d°%.2d'%.2d.%d\"", sign, d, m, s, ts);
1153e12c5d1SDavid du Colombier 	return buf;
1163e12c5d1SDavid du Colombier }
1173e12c5d1SDavid du Colombier 
1183e12c5d1SDavid du Colombier char*
119*219b2ee8SDavid du Colombier ms(Angle a)
120*219b2ee8SDavid du Colombier {
121*219b2ee8SDavid du Colombier 	static char buf[20];
122*219b2ee8SDavid du Colombier 	double x;
123*219b2ee8SDavid du Colombier 	int d, m, s, ts;
124*219b2ee8SDavid du Colombier 
125*219b2ee8SDavid du Colombier 	x = DEG(a);
126*219b2ee8SDavid du Colombier 	x += 0.5/36000.;	/* round up half of 0.1 arcsecond */
127*219b2ee8SDavid du Colombier 	d = floor(x);
128*219b2ee8SDavid du Colombier 	x -= d;
129*219b2ee8SDavid du Colombier 	x *= 60;
130*219b2ee8SDavid du Colombier 	m = floor(x);
131*219b2ee8SDavid du Colombier 	x -= m;
132*219b2ee8SDavid du Colombier 	x *= 60;
133*219b2ee8SDavid du Colombier 	s = floor(x);
134*219b2ee8SDavid du Colombier 	x -= s;
135*219b2ee8SDavid du Colombier 	ts = floor(10*x);
136*219b2ee8SDavid du Colombier 	if(d != 0)
137*219b2ee8SDavid du Colombier 		sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);
138*219b2ee8SDavid du Colombier 	else
139*219b2ee8SDavid du Colombier 		sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
140*219b2ee8SDavid du Colombier 	return buf;
141*219b2ee8SDavid du Colombier }
142*219b2ee8SDavid du Colombier 
143*219b2ee8SDavid du Colombier char*
1443e12c5d1SDavid du Colombier hm(Angle a)
1453e12c5d1SDavid du Colombier {
1463e12c5d1SDavid du Colombier 	static char buf[20];
147*219b2ee8SDavid du Colombier 	double x;
148*219b2ee8SDavid du Colombier 	int h, m, n;
1493e12c5d1SDavid du Colombier 
150*219b2ee8SDavid du Colombier 	x = DEG(a)/15;
151*219b2ee8SDavid du Colombier 	x += 0.5/600.;	/* round up half of tenth of minute */
1523e12c5d1SDavid du Colombier 	h = floor(x);
1533e12c5d1SDavid du Colombier 	x -= h;
1543e12c5d1SDavid du Colombier 	x *= 60;
1553e12c5d1SDavid du Colombier 	m = floor(x);
156*219b2ee8SDavid du Colombier 	x -= m;
157*219b2ee8SDavid du Colombier 	x *= 10;
158*219b2ee8SDavid du Colombier 	n = floor(x);
159*219b2ee8SDavid du Colombier 	sprint(buf, "%dh%.2d.%1dm", h, m, n);
1603e12c5d1SDavid du Colombier 	return buf;
1613e12c5d1SDavid du Colombier }
1623e12c5d1SDavid du Colombier 
1633e12c5d1SDavid du Colombier char*
1643e12c5d1SDavid du Colombier dm(Angle a)
1653e12c5d1SDavid du Colombier {
1663e12c5d1SDavid du Colombier 	static char buf[20];
167*219b2ee8SDavid du Colombier 	double x;
1683e12c5d1SDavid du Colombier 	int sign, d, m, n;
1693e12c5d1SDavid du Colombier 
170*219b2ee8SDavid du Colombier 	x = DEG(a);
1713e12c5d1SDavid du Colombier 	sign='+';
1723e12c5d1SDavid du Colombier 	if(a<0){
1733e12c5d1SDavid du Colombier 		sign='-';
1743e12c5d1SDavid du Colombier 		x=-x;
1753e12c5d1SDavid du Colombier 	}
176*219b2ee8SDavid du Colombier 	x += 0.5/600.;	/* round up half of tenth of arcminute */
1773e12c5d1SDavid du Colombier 	d = floor(x);
1783e12c5d1SDavid du Colombier 	x -= d;
1793e12c5d1SDavid du Colombier 	x *= 60;
1803e12c5d1SDavid du Colombier 	m = floor(x);
1813e12c5d1SDavid du Colombier 	x -= m;
1823e12c5d1SDavid du Colombier 	x *= 10;
1833e12c5d1SDavid du Colombier 	n = floor(x);
1843e12c5d1SDavid du Colombier 	sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);
1853e12c5d1SDavid du Colombier 	return buf;
1863e12c5d1SDavid du Colombier }
187*219b2ee8SDavid du Colombier 
188*219b2ee8SDavid du Colombier char*
189*219b2ee8SDavid du Colombier getword(char *ou, char *in)
190*219b2ee8SDavid du Colombier {
191*219b2ee8SDavid du Colombier 	int c;
192*219b2ee8SDavid du Colombier 
193*219b2ee8SDavid du Colombier 	for(;;) {
194*219b2ee8SDavid du Colombier 		c = *in++;
195*219b2ee8SDavid du Colombier 		if(c == ' ' || c == '\t')
196*219b2ee8SDavid du Colombier 			continue;
197*219b2ee8SDavid du Colombier 		if(c == 0)
198*219b2ee8SDavid du Colombier 			return 0;
199*219b2ee8SDavid du Colombier 		break;
200*219b2ee8SDavid du Colombier 	}
201*219b2ee8SDavid du Colombier 
202*219b2ee8SDavid du Colombier 	if(c == '\'')
203*219b2ee8SDavid du Colombier 		for(;;) {
204*219b2ee8SDavid du Colombier 			if(c >= 'A' && c <= 'Z')
205*219b2ee8SDavid du Colombier 				c += 'a' - 'A';
206*219b2ee8SDavid du Colombier 			*ou++ = c;
207*219b2ee8SDavid du Colombier 			c = *in++;
208*219b2ee8SDavid du Colombier 			if(c == 0)
209*219b2ee8SDavid du Colombier 				return 0;
210*219b2ee8SDavid du Colombier 			if(c == '\'') {
211*219b2ee8SDavid du Colombier 				*ou = 0;
212*219b2ee8SDavid du Colombier 				return in-1;
213*219b2ee8SDavid du Colombier 			}
214*219b2ee8SDavid du Colombier 		}
215*219b2ee8SDavid du Colombier 	for(;;) {
216*219b2ee8SDavid du Colombier 		if(c >= 'A' && c <= 'Z')
217*219b2ee8SDavid du Colombier 			c += 'a' - 'A';
218*219b2ee8SDavid du Colombier 		*ou++ = c;
219*219b2ee8SDavid du Colombier 		c = *in++;
220*219b2ee8SDavid du Colombier 		if(c == ' ' || c == '\t' || c == 0) {
221*219b2ee8SDavid du Colombier 			*ou = 0;
222*219b2ee8SDavid du Colombier 			return in-1;
223*219b2ee8SDavid du Colombier 		}
224*219b2ee8SDavid du Colombier 	}
225*219b2ee8SDavid du Colombier }
226*219b2ee8SDavid du Colombier 
227*219b2ee8SDavid du Colombier /*
228*219b2ee8SDavid du Colombier  * Read formatted angle.  Must contain no embedded blanks
229*219b2ee8SDavid du Colombier  */
230*219b2ee8SDavid du Colombier Angle
231*219b2ee8SDavid du Colombier getra(char *p)
232*219b2ee8SDavid du Colombier {
233*219b2ee8SDavid du Colombier 	Rune r;
234*219b2ee8SDavid du Colombier 	char *q;
235*219b2ee8SDavid du Colombier 	Angle f, d;
236*219b2ee8SDavid du Colombier 	int neg;
237*219b2ee8SDavid du Colombier 
238*219b2ee8SDavid du Colombier 	neg = 0;
239*219b2ee8SDavid du Colombier 	d = 0;
240*219b2ee8SDavid du Colombier 	while(*p == ' ')
241*219b2ee8SDavid du Colombier 		p++;
242*219b2ee8SDavid du Colombier 	for(;;) {
243*219b2ee8SDavid du Colombier 		if(*p == ' ')
244*219b2ee8SDavid du Colombier 			goto Return;
245*219b2ee8SDavid du Colombier 		if(*p == '-') {
246*219b2ee8SDavid du Colombier 			neg = 1;
247*219b2ee8SDavid du Colombier 			p++;
248*219b2ee8SDavid du Colombier 		}
249*219b2ee8SDavid du Colombier 		q = p;
250*219b2ee8SDavid du Colombier 		f = strtod(p, &q);
251*219b2ee8SDavid du Colombier 		if(q > p) {
252*219b2ee8SDavid du Colombier 			p = q;
253*219b2ee8SDavid du Colombier 		}
254*219b2ee8SDavid du Colombier 		p += chartorune(&r, p);
255*219b2ee8SDavid du Colombier 		switch(r) {
256*219b2ee8SDavid du Colombier 		default:
257*219b2ee8SDavid du Colombier 		Return:
258*219b2ee8SDavid du Colombier 			if(neg)
259*219b2ee8SDavid du Colombier 				d = -d;
260*219b2ee8SDavid du Colombier 			return RAD(d);
261*219b2ee8SDavid du Colombier 		case 'h':
262*219b2ee8SDavid du Colombier 			d += f*15;
263*219b2ee8SDavid du Colombier 			break;
264*219b2ee8SDavid du Colombier 		case 'm':
265*219b2ee8SDavid du Colombier 			d += f/4;
266*219b2ee8SDavid du Colombier 			break;
267*219b2ee8SDavid du Colombier 		case 's':
268*219b2ee8SDavid du Colombier 			d += f/240;
269*219b2ee8SDavid du Colombier 			break;
270*219b2ee8SDavid du Colombier 		case L'°':
271*219b2ee8SDavid du Colombier 			d += f;
272*219b2ee8SDavid du Colombier 			break;
273*219b2ee8SDavid du Colombier 		case '\'':
274*219b2ee8SDavid du Colombier 			d += f/60;
275*219b2ee8SDavid du Colombier 			break;
276*219b2ee8SDavid du Colombier 		case '\"':
277*219b2ee8SDavid du Colombier 			d += f/3600;
278*219b2ee8SDavid du Colombier 			break;
279*219b2ee8SDavid du Colombier 		}
280*219b2ee8SDavid du Colombier 	}
281*219b2ee8SDavid du Colombier 	return 0;
282*219b2ee8SDavid du Colombier }
283*219b2ee8SDavid du Colombier 
284*219b2ee8SDavid du Colombier double
285*219b2ee8SDavid du Colombier xsqrt(double a)
286*219b2ee8SDavid du Colombier {
287*219b2ee8SDavid du Colombier 
288*219b2ee8SDavid du Colombier 	if(a < 0)
289*219b2ee8SDavid du Colombier 		return 0;
290*219b2ee8SDavid du Colombier 	return sqrt(a);
291*219b2ee8SDavid du Colombier }
292*219b2ee8SDavid du Colombier 
293*219b2ee8SDavid du Colombier Angle
294*219b2ee8SDavid du Colombier dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
295*219b2ee8SDavid du Colombier {
296*219b2ee8SDavid du Colombier 	double a;
297*219b2ee8SDavid du Colombier 
298*219b2ee8SDavid du Colombier 	a = sin(dec1) * sin(dec2) +
299*219b2ee8SDavid du Colombier 		cos(dec1) * cos(dec2) *
300*219b2ee8SDavid du Colombier 		cos(ra1 - ra2);
301*219b2ee8SDavid du Colombier 	a = atan2(xsqrt(1 - a*a), a);
302*219b2ee8SDavid du Colombier 	if(a < 0)
303*219b2ee8SDavid du Colombier 		a = -a;
304*219b2ee8SDavid du Colombier 	return a;
305*219b2ee8SDavid du Colombier }
306*219b2ee8SDavid du Colombier 
307*219b2ee8SDavid du Colombier int
308*219b2ee8SDavid du Colombier dogamma(Pix c)
309*219b2ee8SDavid du Colombier {
310*219b2ee8SDavid du Colombier 	float f;
311*219b2ee8SDavid du Colombier 
312*219b2ee8SDavid du Colombier 	f = c - gam.min;
313*219b2ee8SDavid du Colombier 	if(f < 1)
314*219b2ee8SDavid du Colombier 		f = 1;
315*219b2ee8SDavid du Colombier 
316*219b2ee8SDavid du Colombier 	if(gam.absgamma == 1)
317*219b2ee8SDavid du Colombier 		c = f * gam.mult2;
318*219b2ee8SDavid du Colombier 	else
319*219b2ee8SDavid du Colombier 		c = exp(log(f*gam.mult1) * gam.absgamma) * 255;
320*219b2ee8SDavid du Colombier 	if(c > 255)
321*219b2ee8SDavid du Colombier 		c = 255;
322*219b2ee8SDavid du Colombier 	if(gam.neg)
323*219b2ee8SDavid du Colombier 		c = 255-c;
324*219b2ee8SDavid du Colombier 	return c;
325*219b2ee8SDavid du Colombier }
326