1 /* misc handy functions.
2 * every system has such, no?
3 * 4/20/98 now_lst() always just returns apparent time
4 */
5
6 #include "astro.h"
7
8 #include <stdio.h>
9 #include <math.h>
10 #include <stdlib.h>
11 #include <string.h>
12
13 static union {
14 unsigned char bytes[sizeof(double)];
15 double value;
16 } _nan = {
17 #if (defined(__s390__) || defined(__s390x__) || defined(__zarch__))
18 {0x7f, 0xf8, 0, 0, 0, 0, 0, 0}
19 #else
20 {0, 0, 0, 0, 0, 0, 0xf8, 0x7f}
21 #endif
22 };
23
24 double ascii_strtod(const char *s00, char **se); /* for PyEphem */
25
26 /* zero from loc for len bytes */
27 void
zero_mem(void * loc,unsigned len)28 zero_mem (void *loc, unsigned len)
29 {
30 (void) memset (loc, 0, len);
31 }
32
33 /* given min and max and an approximate number of divisions desired,
34 * fill in ticks[] with nicely spaced values and return how many.
35 * N.B. return value, and hence number of entries to ticks[], might be as
36 * much as 2 more than numdiv.
37 */
38 int
tickmarks(double min,double max,int numdiv,double ticks[])39 tickmarks (double min, double max, int numdiv, double ticks[])
40 {
41 static int factor[] = { 1, 2, 5 };
42 double minscale;
43 double delta;
44 double lo;
45 double v;
46 int n;
47
48 minscale = fabs(max - min);
49 delta = minscale/numdiv;
50 for (n=0; n < (int)(sizeof(factor)/sizeof(factor[0])); n++) {
51 double scale;
52 double x = delta/factor[n];
53 if ((scale = (pow(10.0, ceil(log10(x)))*factor[n])) < minscale)
54 minscale = scale;
55 }
56 delta = minscale;
57
58 lo = floor(min/delta);
59 for (n = 0; (v = delta*(lo+n)) < max+delta; )
60 ticks[n++] = v;
61
62 return (n);
63 }
64
65 /* given an Obj *, return its type as a descriptive string.
66 * if it's of type fixed then return its class description.
67 * N.B. we return the address of static storage -- do not free or change.
68 */
69 char *
obj_description(Obj * op)70 obj_description (Obj *op)
71 {
72 typedef struct {
73 char classcode;
74 char *desc;
75 } CC;
76
77 #define NFCM ((int)(sizeof(fixed_class_map)/sizeof(fixed_class_map[0])))
78 static CC fixed_class_map[] = {
79 {'A', "Cluster of Galaxies"},
80 {'B', "Binary System"},
81 {'C', "Globular Cluster"},
82 {'D', "Double Star"},
83 {'F', "Diffuse Nebula"},
84 {'G', "Spiral Galaxy"},
85 {'H', "Spherical Galaxy"},
86 {'J', "Radio"},
87 {'K', "Dark Nebula"},
88 {'L', "Pulsar"},
89 {'M', "Multiple Star"},
90 {'N', "Bright Nebula"},
91 {'O', "Open Cluster"},
92 {'P', "Planetary Nebula"},
93 {'Q', "Quasar"},
94 {'R', "Supernova Remnant"},
95 {'S', "Star"},
96 {'T', "Star-like Object"},
97 {'U', "Cluster, with nebulosity"},
98 {'V', "Variable Star"},
99 {'Y', "Supernova"},
100 };
101
102 #define NBCM ((int)(sizeof(binary_class_map)/sizeof(binary_class_map[0])))
103 static CC binary_class_map[] = {
104 {'a', "Astrometric binary"},
105 {'c', "Cataclysmic variable"},
106 {'e', "Eclipsing binary"},
107 {'x', "High-mass X-ray binary"},
108 {'y', "Low-mass X-ray binary"},
109 {'o', "Occultation binary"},
110 {'s', "Spectroscopic binary"},
111 {'t', "1-line spectral binary"},
112 {'u', "2-line spectral binary"},
113 {'v', "Spectrum binary"},
114 {'b', "Visual binary"},
115 {'d', "Visual binary, apparent"},
116 {'q', "Visual binary, optical"},
117 {'r', "Visual binary, physical"},
118 {'p', "Exoplanet"},
119 };
120
121 switch (op->o_type) {
122 case FIXED:
123 if (op->f_class) {
124 int i;
125 for (i = 0; i < NFCM; i++)
126 if (fixed_class_map[i].classcode == op->f_class)
127 return (fixed_class_map[i].desc);
128 }
129 return ("Fixed");
130 case PARABOLIC:
131 return ("Solar - Parabolic");
132 case HYPERBOLIC:
133 return ("Solar - Hyperbolic");
134 case ELLIPTICAL:
135 return ("Solar - Elliptical");
136 case BINARYSTAR:
137 if (op->f_class) {
138 int i;
139 for (i = 0; i < NFCM; i++)
140 if (binary_class_map[i].classcode == op->f_class)
141 return (binary_class_map[i].desc);
142 }
143 return ("Binary system");
144 case PLANET: {
145 static char nsstr[MAXNM + 9];
146 static Obj *biop;
147
148 if (op->pl_code == SUN)
149 return ("Star");
150 if (op->pl_code == MOON)
151 return ("Moon of Earth");
152 if (op->pl_moon == X_PLANET)
153 return ("Planet");
154 if (!biop)
155 getBuiltInObjs (&biop);
156 sprintf (nsstr, "Moon of %s", biop[op->pl_code].o_name);
157 return (nsstr);
158 }
159 case EARTHSAT:
160 return ("Earth Sat");
161 default:
162 printf ("obj_description: unknown type: 0x%x\n", op->o_type);
163 abort();
164 return (NULL); /* for lint */
165 }
166 }
167
168 /* given a Now *, find the local apparent sidereal time, in hours.
169 */
170 void
now_lst(Now * np,double * lstp)171 now_lst (Now *np, double *lstp)
172 {
173 static double last_mjd = -23243, last_lng = 121212, last_lst;
174 double eps, lst, deps, dpsi;
175
176 if (last_mjd == mjd && last_lng == lng) {
177 *lstp = last_lst;
178 return;
179 }
180
181 utc_gst (mjd_day(mjd), mjd_hr(mjd), &lst);
182 lst += radhr(lng);
183
184 obliquity(mjd, &eps);
185 nutation(mjd, &deps, &dpsi);
186 lst += radhr(dpsi*cos(eps+deps));
187
188 range (&lst, 24.0);
189
190 last_mjd = mjd;
191 last_lng = lng;
192 *lstp = last_lst = lst;
193 }
194
195 /* convert ra to ha, in range 0..2*PI
196 * need dec too if not already apparent.
197 */
198 void
radec2ha(Now * np,double ra,double dec,double * hap)199 radec2ha (Now *np, double ra, double dec, double *hap)
200 {
201 double ha, lst;
202
203 if (epoch != EOD)
204 as_ap (np, epoch, &ra, &dec);
205 now_lst (np, &lst);
206 ha = hrrad(lst) - ra;
207 if (ha < 0)
208 ha += 2*PI;
209 *hap = ha;
210 }
211
212 /* find Greenwich Hour Angle of the given object at the given time, 0..2*PI.
213 */
214 void
gha(Now * np,Obj * op,double * ghap)215 gha (Now *np, Obj *op, double *ghap)
216 {
217 Now n = *np;
218 Obj o = *op;
219 double tmp;
220
221 n.n_epoch = EOD;
222 n.n_lng = 0.0;
223 n.n_lat = 0.0;
224 obj_cir (&n, &o);
225 now_lst (&n, &tmp);
226 tmp = hrrad(tmp) - o.s_ra;
227 if (tmp < 0)
228 tmp += 2*PI;
229 *ghap = tmp;
230 }
231
232 /* given a circle and a line segment, find a segment of the line inside the
233 * circle.
234 * return 0 and the segment end points if one exists, else -1.
235 * We use a parametric representation of the line:
236 * x = x1 + (x2-x1)*t and y = y1 + (y2-y1)*t, 0 < t < 1
237 * and a centered representation of the circle:
238 * (x - xc)**2 + (y - yc)**2 = r**2
239 * and solve for the t's that work, checking for usual conditions.
240 */
241 int
lc(int cx,int cy,int cw,int x1,int y1,int x2,int y2,int * sx1,int * sy1,int * sx2,int * sy2)242 lc (
243 int cx, int cy, int cw, /* circle bbox corner and width */
244 int x1, int y1, int x2, int y2, /* line segment endpoints */
245 int *sx1, int *sy1, int *sx2, int *sy2) /* segment inside the circle */
246 {
247 int dx = x2 - x1;
248 int dy = y2 - y1;
249 int r = cw/2;
250 int xc = cx + r;
251 int yc = cy + r;
252 int A = x1 - xc;
253 int B = y1 - yc;
254 double a = dx*dx + dy*dy; /* O(2 * 2**16 * 2**16) */
255 double b = 2*(dx*A + dy*B); /* O(4 * 2**16 * 2**16) */
256 double c = A*A + B*B - r*r; /* O(2 * 2**16 * 2**16) */
257 double d = b*b - 4*a*c; /* O(2**32 * 2**32) */
258 double sqrtd;
259 double t1, t2;
260
261 if (d <= 0)
262 return (-1); /* containing line is purely outside circle */
263
264 sqrtd = sqrt(d);
265 t1 = (-b - sqrtd)/(2.0*a);
266 t2 = (-b + sqrtd)/(2.0*a);
267
268 if (t1 >= 1.0 || t2 <= 0.0)
269 return (-1); /* segment is purely outside circle */
270
271 /* we know now that some part of the segment is inside,
272 * ie, t1 < 1 && t2 > 0
273 */
274
275 if (t1 <= 0.0) {
276 /* (x1,y1) is inside circle */
277 *sx1 = x1;
278 *sy1 = y1;
279 } else {
280 *sx1 = (int)(x1 + dx*t1);
281 *sy1 = (int)(y1 + dy*t1);
282 }
283
284 if (t2 >= 1.0) {
285 /* (x2,y2) is inside circle */
286 *sx2 = x2;
287 *sy2 = y2;
288 } else {
289 *sx2 = (int)(x1 + dx*t2);
290 *sy2 = (int)(y1 + dy*t2);
291 }
292
293 return (0);
294 }
295
296 /* compute visual magnitude using the H/G parameters used in the Astro Almanac.
297 * these are commonly used for asteroids.
298 */
299 void
hg_mag(double h,double g,double rp,double rho,double rsn,double * mp)300 hg_mag (
301 double h, double g,
302 double rp, /* sun-obj dist, AU */
303 double rho, /* earth-obj dist, AU */
304 double rsn, /* sun-earth dist, AU */
305 double *mp)
306 {
307 double psi_t, Psi_1, Psi_2, beta;
308 double c;
309 double tb2;
310
311 c = (rp*rp + rho*rho - rsn*rsn)/(2*rp*rho);
312 if (c <= -1)
313 beta = PI;
314 else if (c >= 1)
315 beta = 0;
316 else
317 beta = acos(c);;
318 tb2 = tan(beta/2.0);
319 /* psi_t = exp(log(tan(beta/2.0))*0.63); */
320 psi_t = pow (tb2, 0.63);
321 Psi_1 = exp(-3.33*psi_t);
322 /* psi_t = exp(log(tan(beta/2.0))*1.22); */
323 psi_t = pow (tb2, 1.22);
324 Psi_2 = exp(-1.87*psi_t);
325 *mp = h + 5.0*log10(rp*rho);
326 if (Psi_1 || Psi_2) *mp -= 2.5*log10((1-g)*Psi_1 + g*Psi_2);
327 }
328
329 /* given faintest desired mag, mag step magstp, image scale and object
330 * magnitude and size, return diameter to draw object, in pixels, or 0 if
331 * dimmer than fmag.
332 */
333 int
magdiam(int fmag,int magstp,double scale,double mag,double size)334 magdiam (
335 int fmag, /* faintest mag */
336 int magstp, /* mag range per dot size */
337 double scale, /* rads per pixel */
338 double mag, /* magnitude */
339 double size) /* rads, or 0 */
340 {
341 int diam, sized;
342
343 if (mag > fmag)
344 return (0);
345 diam = (int)((fmag - mag)/magstp + 1);
346 sized = (int)(size/scale + 0.5);
347 if (sized > diam)
348 diam = sized;
349
350 return (diam);
351 }
352
353 /* computer visual magnitude using the g/k parameters commonly used for comets.
354 */
355 void
gk_mag(double g,double k,double rp,double rho,double * mp)356 gk_mag (
357 double g, double k,
358 double rp, /* sun-obj dist, AU */
359 double rho, /* earth-obj dist, AU */
360 double *mp)
361 {
362 *mp = g + 5.0*log10(rho) + 2.5*k*log10(rp);
363 }
364
365 /* given a string convert to floating point and return it as a double.
366 * this is to isolate possible unportabilities associated with declaring atof().
367 * it's worth it because atof() is often some 50% faster than sscanf ("%lf");
368 */
369 double
atod(char * buf)370 atod (char *buf)
371 {
372 if (*buf == '\0') return _nan.value;
373 return (ascii_strtod(buf, NULL));
374 }
375
376 /* solve a spherical triangle:
377 * A
378 * / \
379 * / \
380 * c / \ b
381 * / \
382 * / \
383 * B ____________ C
384 * a
385 *
386 * given A, b, c find B and a in range 0..B..2PI and 0..a..PI, respectively..
387 * cap and Bp may be NULL if not interested in either one.
388 * N.B. we pass in cos(c) and sin(c) because in many problems one of the sides
389 * remains constant for many values of A and b.
390 */
391 void
solve_sphere(double A,double b,double cc,double sc,double * cap,double * Bp)392 solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp)
393 {
394 double cb = cos(b), sb = sin(b);
395 double sA, cA = cos(A);
396 double x, y;
397 double ca;
398 double B;
399
400 ca = cb*cc + sb*sc*cA;
401 if (ca > 1.0) ca = 1.0;
402 if (ca < -1.0) ca = -1.0;
403 if (cap)
404 *cap = ca;
405
406 if (!Bp)
407 return;
408
409 if (sc < 1e-7)
410 B = cc < 0 ? A : PI-A;
411 else {
412 sA = sin(A);
413 y = sA*sb*sc;
414 x = cb - ca*cc;
415 B = y ? (x ? atan2(y,x) : (y>0 ? PI/2 : -PI/2)) : (x>=0 ? 0 : PI);
416 }
417
418 *Bp = B;
419 range (Bp, 2*PI);
420 }
421
422 /* #define WANT_MATHERR if your system supports it. it gives SGI fits.
423 */
424 #undef WANT_MATHERR
425 #if defined(WANT_MATHERR)
426 /* attempt to do *something* reasonable when a math function blows.
427 */
428 matherr (xp)
429 struct exception *xp;
430 {
431 static char *names[8] = {
432 "acos", "asin", "atan2", "pow",
433 "exp", "log", "log10", "sqrt"
434 };
435 int i;
436
437 /* catch-all */
438 xp->retval = 0.0;
439
440 for (i = 0; i < sizeof(names)/sizeof(names[0]); i++)
441 if (strcmp (xp->name, names[i]) == 0)
442 switch (i) {
443 case 0: /* acos */
444 xp->retval = xp->arg1 >= 1.0 ? 0.0 : -PI;
445 break;
446 case 1: /* asin */
447 xp->retval = xp->arg1 >= 1.0 ? PI/2 : -PI/2;
448 break;
449 case 2: /* atan2 */
450 if (xp->arg1 == 0.0)
451 xp->retval = xp->arg2 < 0.0 ? PI : 0.0;
452 else if (xp->arg2 == 0.0)
453 xp->retval = xp->arg1 < 0.0 ? -PI/2 : PI/2;
454 else
455 xp->retval = 0.0;
456 break;
457 case 3: /* pow */
458 /* FALLTHRU */
459 case 4: /* exp */
460 xp->retval = xp->o_type == OVERFLOW ? 1e308 : 0.0;
461 break;
462 case 5: /* log */
463 /* FALLTHRU */
464 case 6: /* log10 */
465 xp->retval = xp->arg1 <= 0.0 ? -1e308 : 0;
466 break;
467 case 7: /* sqrt */
468 xp->retval = 0.0;
469 break;
470 }
471
472 return (1); /* suppress default error handling */
473 }
474 #endif
475
476 /* given the difference in two RA's, in rads, return their difference,
477 * accounting for wrap at 2*PI. caller need *not* first force it into the
478 * range 0..2*PI.
479 */
480 double
delra(double dra)481 delra (double dra)
482 {
483 double fdra = fmod(fabs(dra), 2*PI);
484
485 if (fdra > PI)
486 fdra = 2*PI - fdra;
487 return (fdra);
488 }
489
490 /* return 1 if object is considered to be "deep sky", else 0.
491 * The only things deep-sky are fixed objects other than stars.
492 */
493 int
is_deepsky(Obj * op)494 is_deepsky (Obj *op)
495 {
496 int deepsky = 0;
497
498 if (is_type(op, FIXEDM)) {
499 switch (op->f_class) {
500 case 'T':
501 case 'B':
502 case 'D':
503 case 'M':
504 case 'S':
505 case 'V':
506 break;
507 default:
508 deepsky = 1;
509 break;
510 }
511 }
512
513 return (deepsky);
514 }
515
516