source: Sophya/trunk/SophyaExt/XephemAstroLib/misc.c@ 3891

Last change on this file since 3891 was 3654, checked in by cmv, 16 years ago

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