xref: /plan9/sys/src/cmd/scat/util.c (revision b85a83648eec38fe82b6f00adfd7828ceec5ee8d)
13e12c5d1SDavid du Colombier #include <u.h>
23e12c5d1SDavid du Colombier #include <libc.h>
3219b2ee8SDavid 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;
8219b2ee8SDavid du Colombier double	LN2	= 0.69314718055994530941723212145817656807550013436025;
9*59cc4ca5SDavid du Colombier static double angledangle=(180./PI)*MILLIARCSEC;
103e12c5d1SDavid du Colombier 
113e12c5d1SDavid du Colombier int
rint(char * p,int n)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
dangle(Angle angle)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
angle(DAngle dangle)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
rfloat(char * p,int n)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
sign(int c)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*
hms(Angle a)703e12c5d1SDavid du Colombier hms(Angle a)
713e12c5d1SDavid du Colombier {
723e12c5d1SDavid du Colombier 	static char buf[20];
73219b2ee8SDavid du Colombier 	double x;
743e12c5d1SDavid du Colombier 	int h, m, s, ts;
753e12c5d1SDavid du Colombier 
76219b2ee8SDavid du Colombier 	x=DEG(a)/15;
77219b2ee8SDavid 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*
dms(Angle a)923e12c5d1SDavid du Colombier dms(Angle a)
933e12c5d1SDavid du Colombier {
943e12c5d1SDavid du Colombier 	static char buf[20];
95219b2ee8SDavid du Colombier 	double x;
963e12c5d1SDavid du Colombier 	int sign, d, m, s, ts;
973e12c5d1SDavid du Colombier 
98219b2ee8SDavid 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 	}
104219b2ee8SDavid 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*
ms(Angle a)119219b2ee8SDavid du Colombier ms(Angle a)
120219b2ee8SDavid du Colombier {
121219b2ee8SDavid du Colombier 	static char buf[20];
122219b2ee8SDavid du Colombier 	double x;
123219b2ee8SDavid du Colombier 	int d, m, s, ts;
124219b2ee8SDavid du Colombier 
125219b2ee8SDavid du Colombier 	x = DEG(a);
126219b2ee8SDavid du Colombier 	x += 0.5/36000.;	/* round up half of 0.1 arcsecond */
127219b2ee8SDavid du Colombier 	d = floor(x);
128219b2ee8SDavid du Colombier 	x -= d;
129219b2ee8SDavid du Colombier 	x *= 60;
130219b2ee8SDavid du Colombier 	m = floor(x);
131219b2ee8SDavid du Colombier 	x -= m;
132219b2ee8SDavid du Colombier 	x *= 60;
133219b2ee8SDavid du Colombier 	s = floor(x);
134219b2ee8SDavid du Colombier 	x -= s;
135219b2ee8SDavid du Colombier 	ts = floor(10*x);
136219b2ee8SDavid du Colombier 	if(d != 0)
137219b2ee8SDavid du Colombier 		sprint(buf, "%d°%.2d'%.2d.%d\"", d, m, s, ts);
138219b2ee8SDavid du Colombier 	else
139219b2ee8SDavid du Colombier 		sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
140219b2ee8SDavid du Colombier 	return buf;
141219b2ee8SDavid du Colombier }
142219b2ee8SDavid du Colombier 
143219b2ee8SDavid du Colombier char*
hm(Angle a)1443e12c5d1SDavid du Colombier hm(Angle a)
1453e12c5d1SDavid du Colombier {
1463e12c5d1SDavid du Colombier 	static char buf[20];
147219b2ee8SDavid du Colombier 	double x;
148219b2ee8SDavid du Colombier 	int h, m, n;
1493e12c5d1SDavid du Colombier 
150219b2ee8SDavid du Colombier 	x = DEG(a)/15;
151219b2ee8SDavid 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);
156219b2ee8SDavid du Colombier 	x -= m;
157219b2ee8SDavid du Colombier 	x *= 10;
158219b2ee8SDavid du Colombier 	n = floor(x);
159219b2ee8SDavid 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*
hm5(Angle a)164*59cc4ca5SDavid du Colombier hm5(Angle a)
165*59cc4ca5SDavid du Colombier {
166*59cc4ca5SDavid du Colombier 	static char buf[20];
167*59cc4ca5SDavid du Colombier 	double x;
168*59cc4ca5SDavid du Colombier 	int h, m;
169*59cc4ca5SDavid du Colombier 
170*59cc4ca5SDavid du Colombier 	x = DEG(a)/15;
171*59cc4ca5SDavid du Colombier 	x += 2.5/60.;	/* round up 2.5m */
172*59cc4ca5SDavid du Colombier 	h = floor(x);
173*59cc4ca5SDavid du Colombier 	x -= h;
174*59cc4ca5SDavid du Colombier 	x *= 60;
175*59cc4ca5SDavid du Colombier 	m = floor(x);
176*59cc4ca5SDavid du Colombier 	m -= m % 5;
177*59cc4ca5SDavid du Colombier 	sprint(buf, "%dh%.2dm", h, m);
178*59cc4ca5SDavid du Colombier 	return buf;
179*59cc4ca5SDavid du Colombier }
180*59cc4ca5SDavid du Colombier 
181*59cc4ca5SDavid du Colombier char*
dm(Angle a)1823e12c5d1SDavid du Colombier dm(Angle a)
1833e12c5d1SDavid du Colombier {
1843e12c5d1SDavid du Colombier 	static char buf[20];
185219b2ee8SDavid du Colombier 	double x;
1863e12c5d1SDavid du Colombier 	int sign, d, m, n;
1873e12c5d1SDavid du Colombier 
188219b2ee8SDavid du Colombier 	x = DEG(a);
1893e12c5d1SDavid du Colombier 	sign='+';
1903e12c5d1SDavid du Colombier 	if(a<0){
1913e12c5d1SDavid du Colombier 		sign='-';
1923e12c5d1SDavid du Colombier 		x=-x;
1933e12c5d1SDavid du Colombier 	}
194219b2ee8SDavid du Colombier 	x += 0.5/600.;	/* round up half of tenth of arcminute */
1953e12c5d1SDavid du Colombier 	d = floor(x);
1963e12c5d1SDavid du Colombier 	x -= d;
1973e12c5d1SDavid du Colombier 	x *= 60;
1983e12c5d1SDavid du Colombier 	m = floor(x);
1993e12c5d1SDavid du Colombier 	x -= m;
2003e12c5d1SDavid du Colombier 	x *= 10;
2013e12c5d1SDavid du Colombier 	n = floor(x);
2023e12c5d1SDavid du Colombier 	sprint(buf, "%c%d°%.2d.%.1d'", sign, d, m, n);
2033e12c5d1SDavid du Colombier 	return buf;
2043e12c5d1SDavid du Colombier }
205219b2ee8SDavid du Colombier 
206219b2ee8SDavid du Colombier char*
deg(Angle a)207*59cc4ca5SDavid du Colombier deg(Angle a)
208*59cc4ca5SDavid du Colombier {
209*59cc4ca5SDavid du Colombier 	static char buf[20];
210*59cc4ca5SDavid du Colombier 	double x;
211*59cc4ca5SDavid du Colombier 	int sign, d;
212*59cc4ca5SDavid du Colombier 
213*59cc4ca5SDavid du Colombier 	x = DEG(a);
214*59cc4ca5SDavid du Colombier 	sign='+';
215*59cc4ca5SDavid du Colombier 	if(a<0){
216*59cc4ca5SDavid du Colombier 		sign='-';
217*59cc4ca5SDavid du Colombier 		x=-x;
218*59cc4ca5SDavid du Colombier 	}
219*59cc4ca5SDavid du Colombier 	x += 0.5;	/* round up half degree */
220*59cc4ca5SDavid du Colombier 	d = floor(x);
221*59cc4ca5SDavid du Colombier 	sprint(buf, "%c%d°", sign, d);
222*59cc4ca5SDavid du Colombier 	return buf;
223*59cc4ca5SDavid du Colombier }
224*59cc4ca5SDavid du Colombier 
225*59cc4ca5SDavid du Colombier char*
getword(char * ou,char * in)226219b2ee8SDavid du Colombier getword(char *ou, char *in)
227219b2ee8SDavid du Colombier {
228219b2ee8SDavid du Colombier 	int c;
229219b2ee8SDavid du Colombier 
230219b2ee8SDavid du Colombier 	for(;;) {
231219b2ee8SDavid du Colombier 		c = *in++;
232219b2ee8SDavid du Colombier 		if(c == ' ' || c == '\t')
233219b2ee8SDavid du Colombier 			continue;
234219b2ee8SDavid du Colombier 		if(c == 0)
235219b2ee8SDavid du Colombier 			return 0;
236219b2ee8SDavid du Colombier 		break;
237219b2ee8SDavid du Colombier 	}
238219b2ee8SDavid du Colombier 
239219b2ee8SDavid du Colombier 	if(c == '\'')
240219b2ee8SDavid du Colombier 		for(;;) {
241219b2ee8SDavid du Colombier 			if(c >= 'A' && c <= 'Z')
242219b2ee8SDavid du Colombier 				c += 'a' - 'A';
243219b2ee8SDavid du Colombier 			*ou++ = c;
244219b2ee8SDavid du Colombier 			c = *in++;
245219b2ee8SDavid du Colombier 			if(c == 0)
246219b2ee8SDavid du Colombier 				return 0;
247219b2ee8SDavid du Colombier 			if(c == '\'') {
248219b2ee8SDavid du Colombier 				*ou = 0;
249219b2ee8SDavid du Colombier 				return in-1;
250219b2ee8SDavid du Colombier 			}
251219b2ee8SDavid du Colombier 		}
252219b2ee8SDavid du Colombier 	for(;;) {
253219b2ee8SDavid du Colombier 		if(c >= 'A' && c <= 'Z')
254219b2ee8SDavid du Colombier 			c += 'a' - 'A';
255219b2ee8SDavid du Colombier 		*ou++ = c;
256219b2ee8SDavid du Colombier 		c = *in++;
257219b2ee8SDavid du Colombier 		if(c == ' ' || c == '\t' || c == 0) {
258219b2ee8SDavid du Colombier 			*ou = 0;
259219b2ee8SDavid du Colombier 			return in-1;
260219b2ee8SDavid du Colombier 		}
261219b2ee8SDavid du Colombier 	}
262219b2ee8SDavid du Colombier }
263219b2ee8SDavid du Colombier 
264219b2ee8SDavid du Colombier /*
265219b2ee8SDavid du Colombier  * Read formatted angle.  Must contain no embedded blanks
266219b2ee8SDavid du Colombier  */
267219b2ee8SDavid du Colombier Angle
getra(char * p)268219b2ee8SDavid du Colombier getra(char *p)
269219b2ee8SDavid du Colombier {
270219b2ee8SDavid du Colombier 	Rune r;
271219b2ee8SDavid du Colombier 	char *q;
272219b2ee8SDavid du Colombier 	Angle f, d;
273219b2ee8SDavid du Colombier 	int neg;
274219b2ee8SDavid du Colombier 
275219b2ee8SDavid du Colombier 	neg = 0;
276219b2ee8SDavid du Colombier 	d = 0;
277219b2ee8SDavid du Colombier 	while(*p == ' ')
278219b2ee8SDavid du Colombier 		p++;
279219b2ee8SDavid du Colombier 	for(;;) {
280*59cc4ca5SDavid du Colombier 		if(*p == ' ' || *p=='\0')
281219b2ee8SDavid du Colombier 			goto Return;
282219b2ee8SDavid du Colombier 		if(*p == '-') {
283219b2ee8SDavid du Colombier 			neg = 1;
284219b2ee8SDavid du Colombier 			p++;
285219b2ee8SDavid du Colombier 		}
286*59cc4ca5SDavid du Colombier 		if(*p == '+') {
287*59cc4ca5SDavid du Colombier 			neg = 0;
288*59cc4ca5SDavid du Colombier 			p++;
289*59cc4ca5SDavid du Colombier 		}
290219b2ee8SDavid du Colombier 		q = p;
291219b2ee8SDavid du Colombier 		f = strtod(p, &q);
292219b2ee8SDavid du Colombier 		if(q > p) {
293219b2ee8SDavid du Colombier 			p = q;
294219b2ee8SDavid du Colombier 		}
295219b2ee8SDavid du Colombier 		p += chartorune(&r, p);
296219b2ee8SDavid du Colombier 		switch(r) {
297219b2ee8SDavid du Colombier 		default:
298219b2ee8SDavid du Colombier 		Return:
299219b2ee8SDavid du Colombier 			if(neg)
300219b2ee8SDavid du Colombier 				d = -d;
301219b2ee8SDavid du Colombier 			return RAD(d);
302219b2ee8SDavid du Colombier 		case 'h':
303219b2ee8SDavid du Colombier 			d += f*15;
304219b2ee8SDavid du Colombier 			break;
305219b2ee8SDavid du Colombier 		case 'm':
306219b2ee8SDavid du Colombier 			d += f/4;
307219b2ee8SDavid du Colombier 			break;
308219b2ee8SDavid du Colombier 		case 's':
309219b2ee8SDavid du Colombier 			d += f/240;
310219b2ee8SDavid du Colombier 			break;
311219b2ee8SDavid du Colombier 		case L'°':
312219b2ee8SDavid du Colombier 			d += f;
313219b2ee8SDavid du Colombier 			break;
314219b2ee8SDavid du Colombier 		case '\'':
315219b2ee8SDavid du Colombier 			d += f/60;
316219b2ee8SDavid du Colombier 			break;
317219b2ee8SDavid du Colombier 		case '\"':
318219b2ee8SDavid du Colombier 			d += f/3600;
319219b2ee8SDavid du Colombier 			break;
320219b2ee8SDavid du Colombier 		}
321219b2ee8SDavid du Colombier 	}
322219b2ee8SDavid du Colombier }
323219b2ee8SDavid du Colombier 
324219b2ee8SDavid du Colombier double
xsqrt(double a)325219b2ee8SDavid du Colombier xsqrt(double a)
326219b2ee8SDavid du Colombier {
327219b2ee8SDavid du Colombier 
328219b2ee8SDavid du Colombier 	if(a < 0)
329219b2ee8SDavid du Colombier 		return 0;
330219b2ee8SDavid du Colombier 	return sqrt(a);
331219b2ee8SDavid du Colombier }
332219b2ee8SDavid du Colombier 
333219b2ee8SDavid du Colombier Angle
dist(Angle ra1,Angle dec1,Angle ra2,Angle dec2)334219b2ee8SDavid du Colombier dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
335219b2ee8SDavid du Colombier {
336219b2ee8SDavid du Colombier 	double a;
337219b2ee8SDavid du Colombier 
338219b2ee8SDavid du Colombier 	a = sin(dec1) * sin(dec2) +
339219b2ee8SDavid du Colombier 		cos(dec1) * cos(dec2) *
340219b2ee8SDavid du Colombier 		cos(ra1 - ra2);
341219b2ee8SDavid du Colombier 	a = atan2(xsqrt(1 - a*a), a);
342219b2ee8SDavid du Colombier 	if(a < 0)
343219b2ee8SDavid du Colombier 		a = -a;
344219b2ee8SDavid du Colombier 	return a;
345219b2ee8SDavid du Colombier }
346219b2ee8SDavid du Colombier 
347219b2ee8SDavid du Colombier int
dogamma(Pix c)348219b2ee8SDavid du Colombier dogamma(Pix c)
349219b2ee8SDavid du Colombier {
350219b2ee8SDavid du Colombier 	float f;
351219b2ee8SDavid du Colombier 
352219b2ee8SDavid du Colombier 	f = c - gam.min;
353219b2ee8SDavid du Colombier 	if(f < 1)
354219b2ee8SDavid du Colombier 		f = 1;
355219b2ee8SDavid du Colombier 
356219b2ee8SDavid du Colombier 	if(gam.absgamma == 1)
357219b2ee8SDavid du Colombier 		c = f * gam.mult2;
358219b2ee8SDavid du Colombier 	else
359219b2ee8SDavid du Colombier 		c = exp(log(f*gam.mult1) * gam.absgamma) * 255;
360219b2ee8SDavid du Colombier 	if(c > 255)
361219b2ee8SDavid du Colombier 		c = 255;
362219b2ee8SDavid du Colombier 	if(gam.neg)
363219b2ee8SDavid du Colombier 		c = 255-c;
364219b2ee8SDavid du Colombier 	return c;
365219b2ee8SDavid du Colombier }
366