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"
  5 8a3b2ceb 2004-04-24 devnull 
  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;
 10 8a3b2ceb 2004-04-24 devnull 
 11 8a3b2ceb 2004-04-24 devnull #define rint scatrint
 12 8a3b2ceb 2004-04-24 devnull 
 13 8a3b2ceb 2004-04-24 devnull int
 14 8a3b2ceb 2004-04-24 devnull rint(char *p, int n)
 15 8a3b2ceb 2004-04-24 devnull {
 16 8a3b2ceb 2004-04-24 devnull 	int i=0;
 17 8a3b2ceb 2004-04-24 devnull 
 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;
 23 8a3b2ceb 2004-04-24 devnull }
 24 8a3b2ceb 2004-04-24 devnull 
 25 8a3b2ceb 2004-04-24 devnull DAngle
 26 8a3b2ceb 2004-04-24 devnull dangle(Angle angle)
 27 8a3b2ceb 2004-04-24 devnull {
 28 8a3b2ceb 2004-04-24 devnull 	return angle*angledangle;
 29 8a3b2ceb 2004-04-24 devnull }
 30 8a3b2ceb 2004-04-24 devnull 
 31 8a3b2ceb 2004-04-24 devnull Angle
 32 8a3b2ceb 2004-04-24 devnull angle(DAngle dangle)
 33 8a3b2ceb 2004-04-24 devnull {
 34 8a3b2ceb 2004-04-24 devnull 	return dangle/angledangle;
 35 8a3b2ceb 2004-04-24 devnull }
 36 8a3b2ceb 2004-04-24 devnull 
 37 8a3b2ceb 2004-04-24 devnull double
 38 8a3b2ceb 2004-04-24 devnull rfloat(char *p, int n)
 39 8a3b2ceb 2004-04-24 devnull {
 40 8a3b2ceb 2004-04-24 devnull 	double i, d=0;
 41 8a3b2ceb 2004-04-24 devnull 
 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;
 56 8a3b2ceb 2004-04-24 devnull 	p++;
 57 8a3b2ceb 2004-04-24 devnull 	i = 1;
 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;
 61 8a3b2ceb 2004-04-24 devnull }
 62 8a3b2ceb 2004-04-24 devnull 
 63 8a3b2ceb 2004-04-24 devnull int
 64 8a3b2ceb 2004-04-24 devnull sign(int c)
 65 8a3b2ceb 2004-04-24 devnull {
 66 8a3b2ceb 2004-04-24 devnull 	if(c=='-')
 67 8a3b2ceb 2004-04-24 devnull 		return -1;
 68 8a3b2ceb 2004-04-24 devnull 	return 1;
 69 8a3b2ceb 2004-04-24 devnull }
 70 8a3b2ceb 2004-04-24 devnull 
 71 8a3b2ceb 2004-04-24 devnull char*
 72 8a3b2ceb 2004-04-24 devnull hms(Angle a)
 73 8a3b2ceb 2004-04-24 devnull {
 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;
 77 8a3b2ceb 2004-04-24 devnull 
 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);
 81 8a3b2ceb 2004-04-24 devnull 	x -= h;
 82 8a3b2ceb 2004-04-24 devnull 	x *= 60;
 83 8a3b2ceb 2004-04-24 devnull 	m = floor(x);
 84 8a3b2ceb 2004-04-24 devnull 	x -= m;
 85 8a3b2ceb 2004-04-24 devnull 	x *= 60;
 86 8a3b2ceb 2004-04-24 devnull 	s = floor(x);
 87 8a3b2ceb 2004-04-24 devnull 	x -= s;
 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;
 91 8a3b2ceb 2004-04-24 devnull }
 92 8a3b2ceb 2004-04-24 devnull 
 93 8a3b2ceb 2004-04-24 devnull char*
 94 8a3b2ceb 2004-04-24 devnull dms(Angle a)
 95 8a3b2ceb 2004-04-24 devnull {
 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;
 99 8a3b2ceb 2004-04-24 devnull 
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='-';
104 8a3b2ceb 2004-04-24 devnull 		x=-x;
105 8a3b2ceb 2004-04-24 devnull 	}
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);
108 8a3b2ceb 2004-04-24 devnull 	x -= d;
109 8a3b2ceb 2004-04-24 devnull 	x *= 60;
110 8a3b2ceb 2004-04-24 devnull 	m = floor(x);
111 8a3b2ceb 2004-04-24 devnull 	x -= m;
112 8a3b2ceb 2004-04-24 devnull 	x *= 60;
113 8a3b2ceb 2004-04-24 devnull 	s = floor(x);
114 8a3b2ceb 2004-04-24 devnull 	x -= s;
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;
118 8a3b2ceb 2004-04-24 devnull }
119 8a3b2ceb 2004-04-24 devnull 
120 8a3b2ceb 2004-04-24 devnull char*
121 8a3b2ceb 2004-04-24 devnull ms(Angle a)
122 8a3b2ceb 2004-04-24 devnull {
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;
126 8a3b2ceb 2004-04-24 devnull 
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);
130 8a3b2ceb 2004-04-24 devnull 	x -= d;
131 8a3b2ceb 2004-04-24 devnull 	x *= 60;
132 8a3b2ceb 2004-04-24 devnull 	m = floor(x);
133 8a3b2ceb 2004-04-24 devnull 	x -= m;
134 8a3b2ceb 2004-04-24 devnull 	x *= 60;
135 8a3b2ceb 2004-04-24 devnull 	s = floor(x);
136 8a3b2ceb 2004-04-24 devnull 	x -= s;
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);
140 8a3b2ceb 2004-04-24 devnull 	else
141 8a3b2ceb 2004-04-24 devnull 		sprint(buf, "%.2d'%.2d.%d\"", m, s, ts);
142 8a3b2ceb 2004-04-24 devnull 	return buf;
143 8a3b2ceb 2004-04-24 devnull }
144 8a3b2ceb 2004-04-24 devnull 
145 8a3b2ceb 2004-04-24 devnull char*
146 8a3b2ceb 2004-04-24 devnull hm(Angle a)
147 8a3b2ceb 2004-04-24 devnull {
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;
151 8a3b2ceb 2004-04-24 devnull 
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);
155 8a3b2ceb 2004-04-24 devnull 	x -= h;
156 8a3b2ceb 2004-04-24 devnull 	x *= 60;
157 8a3b2ceb 2004-04-24 devnull 	m = floor(x);
158 8a3b2ceb 2004-04-24 devnull 	x -= m;
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;
163 8a3b2ceb 2004-04-24 devnull }
164 8a3b2ceb 2004-04-24 devnull 
165 8a3b2ceb 2004-04-24 devnull char*
166 8a3b2ceb 2004-04-24 devnull hm5(Angle a)
167 8a3b2ceb 2004-04-24 devnull {
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;
171 8a3b2ceb 2004-04-24 devnull 
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);
175 8a3b2ceb 2004-04-24 devnull 	x -= h;
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;
181 8a3b2ceb 2004-04-24 devnull }
182 8a3b2ceb 2004-04-24 devnull 
183 8a3b2ceb 2004-04-24 devnull char*
184 8a3b2ceb 2004-04-24 devnull dm(Angle a)
185 8a3b2ceb 2004-04-24 devnull {
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;
189 8a3b2ceb 2004-04-24 devnull 
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='-';
194 8a3b2ceb 2004-04-24 devnull 		x=-x;
195 8a3b2ceb 2004-04-24 devnull 	}
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);
198 8a3b2ceb 2004-04-24 devnull 	x -= d;
199 8a3b2ceb 2004-04-24 devnull 	x *= 60;
200 8a3b2ceb 2004-04-24 devnull 	m = floor(x);
201 8a3b2ceb 2004-04-24 devnull 	x -= m;
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;
206 8a3b2ceb 2004-04-24 devnull }
207 8a3b2ceb 2004-04-24 devnull 
208 8a3b2ceb 2004-04-24 devnull char*
209 8a3b2ceb 2004-04-24 devnull deg(Angle a)
210 8a3b2ceb 2004-04-24 devnull {
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;
214 8a3b2ceb 2004-04-24 devnull 
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='-';
219 8a3b2ceb 2004-04-24 devnull 		x=-x;
220 8a3b2ceb 2004-04-24 devnull 	}
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;
225 8a3b2ceb 2004-04-24 devnull }
226 8a3b2ceb 2004-04-24 devnull 
227 8a3b2ceb 2004-04-24 devnull char*
228 8a3b2ceb 2004-04-24 devnull getword(char *ou, char *in)
229 8a3b2ceb 2004-04-24 devnull {
230 8a3b2ceb 2004-04-24 devnull 	int c;
231 8a3b2ceb 2004-04-24 devnull 
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;
238 8a3b2ceb 2004-04-24 devnull 		break;
239 8a3b2ceb 2004-04-24 devnull 	}
240 8a3b2ceb 2004-04-24 devnull 
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;
252 8a3b2ceb 2004-04-24 devnull 			}
253 8a3b2ceb 2004-04-24 devnull 		}
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;
262 8a3b2ceb 2004-04-24 devnull 		}
263 8a3b2ceb 2004-04-24 devnull 	}
264 8a3b2ceb 2004-04-24 devnull }
265 8a3b2ceb 2004-04-24 devnull 
266 8a3b2ceb 2004-04-24 devnull /*
267 8a3b2ceb 2004-04-24 devnull  * Read formatted angle.  Must contain no embedded blanks
268 8a3b2ceb 2004-04-24 devnull  */
269 8a3b2ceb 2004-04-24 devnull Angle
270 8a3b2ceb 2004-04-24 devnull getra(char *p)
271 8a3b2ceb 2004-04-24 devnull {
272 8a3b2ceb 2004-04-24 devnull 	Rune r;
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;
276 8a3b2ceb 2004-04-24 devnull 
277 8a3b2ceb 2004-04-24 devnull 	neg = 0;
278 8a3b2ceb 2004-04-24 devnull 	d = 0;
279 8a3b2ceb 2004-04-24 devnull 	while(*p == ' ')
280 8a3b2ceb 2004-04-24 devnull 		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;
286 8a3b2ceb 2004-04-24 devnull 			p++;
287 8a3b2ceb 2004-04-24 devnull 		}
288 8a3b2ceb 2004-04-24 devnull 		if(*p == '+') {
289 8a3b2ceb 2004-04-24 devnull 			neg = 0;
290 8a3b2ceb 2004-04-24 devnull 			p++;
291 8a3b2ceb 2004-04-24 devnull 		}
292 8a3b2ceb 2004-04-24 devnull 		q = p;
293 8a3b2ceb 2004-04-24 devnull 		f = strtod(p, &q);
294 8a3b2ceb 2004-04-24 devnull 		if(q > p) {
295 8a3b2ceb 2004-04-24 devnull 			p = q;
296 8a3b2ceb 2004-04-24 devnull 		}
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:
300 8a3b2ceb 2004-04-24 devnull 		Return:
301 8a3b2ceb 2004-04-24 devnull 			if(neg)
302 8a3b2ceb 2004-04-24 devnull 				d = -d;
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;
306 8a3b2ceb 2004-04-24 devnull 			break;
307 8a3b2ceb 2004-04-24 devnull 		case 'm':
308 8a3b2ceb 2004-04-24 devnull 			d += f/4;
309 8a3b2ceb 2004-04-24 devnull 			break;
310 8a3b2ceb 2004-04-24 devnull 		case 's':
311 8a3b2ceb 2004-04-24 devnull 			d += f/240;
312 8a3b2ceb 2004-04-24 devnull 			break;
313 8a3b2ceb 2004-04-24 devnull 		case 0xB0:	/* ° */
314 8a3b2ceb 2004-04-24 devnull 			d += f;
315 8a3b2ceb 2004-04-24 devnull 			break;
316 8a3b2ceb 2004-04-24 devnull 		case '\'':
317 8a3b2ceb 2004-04-24 devnull 			d += f/60;
318 8a3b2ceb 2004-04-24 devnull 			break;
319 8a3b2ceb 2004-04-24 devnull 		case '\"':
320 8a3b2ceb 2004-04-24 devnull 			d += f/3600;
321 8a3b2ceb 2004-04-24 devnull 			break;
322 8a3b2ceb 2004-04-24 devnull 		}
323 8a3b2ceb 2004-04-24 devnull 	}
324 8a3b2ceb 2004-04-24 devnull 	return 0;
325 8a3b2ceb 2004-04-24 devnull }
326 8a3b2ceb 2004-04-24 devnull 
327 8a3b2ceb 2004-04-24 devnull double
328 8a3b2ceb 2004-04-24 devnull xsqrt(double a)
329 8a3b2ceb 2004-04-24 devnull {
330 8a3b2ceb 2004-04-24 devnull 
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);
334 8a3b2ceb 2004-04-24 devnull }
335 8a3b2ceb 2004-04-24 devnull 
336 8a3b2ceb 2004-04-24 devnull Angle
337 8a3b2ceb 2004-04-24 devnull dist(Angle ra1, Angle dec1, Angle ra2, Angle dec2)
338 8a3b2ceb 2004-04-24 devnull {
339 8a3b2ceb 2004-04-24 devnull 	double a;
340 8a3b2ceb 2004-04-24 devnull 
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)
346 8a3b2ceb 2004-04-24 devnull 		a = -a;
347 8a3b2ceb 2004-04-24 devnull 	return a;
348 8a3b2ceb 2004-04-24 devnull }
349 8a3b2ceb 2004-04-24 devnull 
350 8a3b2ceb 2004-04-24 devnull int
351 8a3b2ceb 2004-04-24 devnull dogamma(Pix c)
352 8a3b2ceb 2004-04-24 devnull {
353 8a3b2ceb 2004-04-24 devnull 	float f;
354 8a3b2ceb 2004-04-24 devnull 
355 8a3b2ceb 2004-04-24 devnull 	f = c - gam.min;
356 8a3b2ceb 2004-04-24 devnull 	if(f < 1)
357 8a3b2ceb 2004-04-24 devnull 		f = 1;
358 8a3b2ceb 2004-04-24 devnull 
359 8a3b2ceb 2004-04-24 devnull 	if(gam.absgamma == 1)
360 8a3b2ceb 2004-04-24 devnull 		c = f * gam.mult2;
361 8a3b2ceb 2004-04-24 devnull 	else
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;
368 8a3b2ceb 2004-04-24 devnull }