/* misc handy functions. * every system has such, no? * 4/20/98 now_lst() always just returns apparent time */ #include "astro.h" #include #include #include #include static union { unsigned char bytes[sizeof(double)]; double value; } _nan = { #if (defined(__s390__) || defined(__s390x__) || defined(__zarch__)) {0x7f, 0xf8, 0, 0, 0, 0, 0, 0} #else {0, 0, 0, 0, 0, 0, 0xf8, 0x7f} #endif }; double ascii_strtod(const char *s00, char **se); /* for PyEphem */ /* zero from loc for len bytes */ void zero_mem (void *loc, unsigned len) { (void) memset (loc, 0, len); } /* given min and max and an approximate number of divisions desired, * fill in ticks[] with nicely spaced values and return how many. * N.B. return value, and hence number of entries to ticks[], might be as * much as 2 more than numdiv. */ int tickmarks (double min, double max, int numdiv, double ticks[]) { static int factor[] = { 1, 2, 5 }; double minscale; double delta; double lo; double v; int n; minscale = fabs(max - min); delta = minscale/numdiv; for (n=0; n < (int)(sizeof(factor)/sizeof(factor[0])); n++) { double scale; double x = delta/factor[n]; if ((scale = (pow(10.0, ceil(log10(x)))*factor[n])) < minscale) minscale = scale; } delta = minscale; lo = floor(min/delta); for (n = 0; (v = delta*(lo+n)) < max+delta; ) ticks[n++] = v; return (n); } /* given an Obj *, return its type as a descriptive string. * if it's of type fixed then return its class description. * N.B. we return the address of static storage -- do not free or change. */ char * obj_description (Obj *op) { typedef struct { char classcode; char *desc; } CC; #define NFCM ((int)(sizeof(fixed_class_map)/sizeof(fixed_class_map[0]))) static CC fixed_class_map[] = { {'A', "Cluster of Galaxies"}, {'B', "Binary System"}, {'C', "Globular Cluster"}, {'D', "Double Star"}, {'F', "Diffuse Nebula"}, {'G', "Spiral Galaxy"}, {'H', "Spherical Galaxy"}, {'J', "Radio"}, {'K', "Dark Nebula"}, {'L', "Pulsar"}, {'M', "Multiple Star"}, {'N', "Bright Nebula"}, {'O', "Open Cluster"}, {'P', "Planetary Nebula"}, {'Q', "Quasar"}, {'R', "Supernova Remnant"}, {'S', "Star"}, {'T', "Star-like Object"}, {'U', "Cluster, with nebulosity"}, {'V', "Variable Star"}, {'Y', "Supernova"}, }; #define NBCM ((int)(sizeof(binary_class_map)/sizeof(binary_class_map[0]))) static CC binary_class_map[] = { {'a', "Astrometric binary"}, {'c', "Cataclysmic variable"}, {'e', "Eclipsing binary"}, {'x', "High-mass X-ray binary"}, {'y', "Low-mass X-ray binary"}, {'o', "Occultation binary"}, {'s', "Spectroscopic binary"}, {'t', "1-line spectral binary"}, {'u', "2-line spectral binary"}, {'v', "Spectrum binary"}, {'b', "Visual binary"}, {'d', "Visual binary, apparent"}, {'q', "Visual binary, optical"}, {'r', "Visual binary, physical"}, {'p', "Exoplanet"}, }; switch (op->o_type) { case FIXED: if (op->f_class) { int i; for (i = 0; i < NFCM; i++) if (fixed_class_map[i].classcode == op->f_class) return (fixed_class_map[i].desc); } return ("Fixed"); case PARABOLIC: return ("Solar - Parabolic"); case HYPERBOLIC: return ("Solar - Hyperbolic"); case ELLIPTICAL: return ("Solar - Elliptical"); case BINARYSTAR: if (op->f_class) { int i; for (i = 0; i < NFCM; i++) if (binary_class_map[i].classcode == op->f_class) return (binary_class_map[i].desc); } return ("Binary system"); case PLANET: { static char nsstr[MAXNM + 9]; static Obj *biop; if (op->pl_code == SUN) return ("Star"); if (op->pl_code == MOON) return ("Moon of Earth"); if (op->pl_moon == X_PLANET) return ("Planet"); if (!biop) getBuiltInObjs (&biop); sprintf (nsstr, "Moon of %s", biop[op->pl_code].o_name); return (nsstr); } case EARTHSAT: return ("Earth Sat"); default: printf ("obj_description: unknown type: 0x%x\n", op->o_type); abort(); return (NULL); /* for lint */ } } /* given a Now *, find the local apparent sidereal time, in hours. */ void now_lst (Now *np, double *lstp) { static double last_mjd = -23243, last_lng = 121212, last_lst; double eps, lst, deps, dpsi; if (last_mjd == mjd && last_lng == lng) { *lstp = last_lst; return; } utc_gst (mjd_day(mjd), mjd_hr(mjd), &lst); lst += radhr(lng); obliquity(mjd, &eps); nutation(mjd, &deps, &dpsi); lst += radhr(dpsi*cos(eps+deps)); range (&lst, 24.0); last_mjd = mjd; last_lng = lng; *lstp = last_lst = lst; } /* convert ra to ha, in range 0..2*PI * need dec too if not already apparent. */ void radec2ha (Now *np, double ra, double dec, double *hap) { double ha, lst; if (epoch != EOD) as_ap (np, epoch, &ra, &dec); now_lst (np, &lst); ha = hrrad(lst) - ra; if (ha < 0) ha += 2*PI; *hap = ha; } /* find Greenwich Hour Angle of the given object at the given time, 0..2*PI. */ void gha (Now *np, Obj *op, double *ghap) { Now n = *np; Obj o = *op; double tmp; n.n_epoch = EOD; n.n_lng = 0.0; n.n_lat = 0.0; obj_cir (&n, &o); now_lst (&n, &tmp); tmp = hrrad(tmp) - o.s_ra; if (tmp < 0) tmp += 2*PI; *ghap = tmp; } /* given a circle and a line segment, find a segment of the line inside the * circle. * return 0 and the segment end points if one exists, else -1. * We use a parametric representation of the line: * x = x1 + (x2-x1)*t and y = y1 + (y2-y1)*t, 0 < t < 1 * and a centered representation of the circle: * (x - xc)**2 + (y - yc)**2 = r**2 * and solve for the t's that work, checking for usual conditions. */ int lc ( int cx, int cy, int cw, /* circle bbox corner and width */ int x1, int y1, int x2, int y2, /* line segment endpoints */ int *sx1, int *sy1, int *sx2, int *sy2) /* segment inside the circle */ { int dx = x2 - x1; int dy = y2 - y1; int r = cw/2; int xc = cx + r; int yc = cy + r; int A = x1 - xc; int B = y1 - yc; double a = dx*dx + dy*dy; /* O(2 * 2**16 * 2**16) */ double b = 2*(dx*A + dy*B); /* O(4 * 2**16 * 2**16) */ double c = A*A + B*B - r*r; /* O(2 * 2**16 * 2**16) */ double d = b*b - 4*a*c; /* O(2**32 * 2**32) */ double sqrtd; double t1, t2; if (d <= 0) return (-1); /* containing line is purely outside circle */ sqrtd = sqrt(d); t1 = (-b - sqrtd)/(2.0*a); t2 = (-b + sqrtd)/(2.0*a); if (t1 >= 1.0 || t2 <= 0.0) return (-1); /* segment is purely outside circle */ /* we know now that some part of the segment is inside, * ie, t1 < 1 && t2 > 0 */ if (t1 <= 0.0) { /* (x1,y1) is inside circle */ *sx1 = x1; *sy1 = y1; } else { *sx1 = (int)(x1 + dx*t1); *sy1 = (int)(y1 + dy*t1); } if (t2 >= 1.0) { /* (x2,y2) is inside circle */ *sx2 = x2; *sy2 = y2; } else { *sx2 = (int)(x1 + dx*t2); *sy2 = (int)(y1 + dy*t2); } return (0); } /* compute visual magnitude using the H/G parameters used in the Astro Almanac. * these are commonly used for asteroids. */ void hg_mag ( double h, double g, double rp, /* sun-obj dist, AU */ double rho, /* earth-obj dist, AU */ double rsn, /* sun-earth dist, AU */ double *mp) { double psi_t, Psi_1, Psi_2, beta; double c; double tb2; c = (rp*rp + rho*rho - rsn*rsn)/(2*rp*rho); if (c <= -1) beta = PI; else if (c >= 1) beta = 0; else beta = acos(c);; tb2 = tan(beta/2.0); /* psi_t = exp(log(tan(beta/2.0))*0.63); */ psi_t = pow (tb2, 0.63); Psi_1 = exp(-3.33*psi_t); /* psi_t = exp(log(tan(beta/2.0))*1.22); */ psi_t = pow (tb2, 1.22); Psi_2 = exp(-1.87*psi_t); *mp = h + 5.0*log10(rp*rho); if (Psi_1 || Psi_2) *mp -= 2.5*log10((1-g)*Psi_1 + g*Psi_2); } /* given faintest desired mag, mag step magstp, image scale and object * magnitude and size, return diameter to draw object, in pixels, or 0 if * dimmer than fmag. */ int magdiam ( int fmag, /* faintest mag */ int magstp, /* mag range per dot size */ double scale, /* rads per pixel */ double mag, /* magnitude */ double size) /* rads, or 0 */ { int diam, sized; if (mag > fmag) return (0); diam = (int)((fmag - mag)/magstp + 1); sized = (int)(size/scale + 0.5); if (sized > diam) diam = sized; return (diam); } /* computer visual magnitude using the g/k parameters commonly used for comets. */ void gk_mag ( double g, double k, double rp, /* sun-obj dist, AU */ double rho, /* earth-obj dist, AU */ double *mp) { *mp = g + 5.0*log10(rho) + 2.5*k*log10(rp); } /* given a string convert to floating point and return it as a double. * this is to isolate possible unportabilities associated with declaring atof(). * it's worth it because atof() is often some 50% faster than sscanf ("%lf"); */ double atod (char *buf) { if (*buf == '\0') return _nan.value; return (ascii_strtod(buf, NULL)); } /* solve a spherical triangle: * A * / \ * / \ * c / \ b * / \ * / \ * B ____________ C * a * * given A, b, c find B and a in range 0..B..2PI and 0..a..PI, respectively.. * cap and Bp may be NULL if not interested in either one. * N.B. we pass in cos(c) and sin(c) because in many problems one of the sides * remains constant for many values of A and b. */ void solve_sphere (double A, double b, double cc, double sc, double *cap, double *Bp) { double cb = cos(b), sb = sin(b); double sA, cA = cos(A); double x, y; double ca; double B; ca = cb*cc + sb*sc*cA; if (ca > 1.0) ca = 1.0; if (ca < -1.0) ca = -1.0; if (cap) *cap = ca; if (!Bp) return; if (sc < 1e-7) B = cc < 0 ? A : PI-A; else { sA = sin(A); y = sA*sb*sc; x = cb - ca*cc; B = y ? (x ? atan2(y,x) : (y>0 ? PI/2 : -PI/2)) : (x>=0 ? 0 : PI); } *Bp = B; range (Bp, 2*PI); } /* #define WANT_MATHERR if your system supports it. it gives SGI fits. */ #undef WANT_MATHERR #if defined(WANT_MATHERR) /* attempt to do *something* reasonable when a math function blows. */ matherr (xp) struct exception *xp; { static char *names[8] = { "acos", "asin", "atan2", "pow", "exp", "log", "log10", "sqrt" }; int i; /* catch-all */ xp->retval = 0.0; for (i = 0; i < sizeof(names)/sizeof(names[0]); i++) if (strcmp (xp->name, names[i]) == 0) switch (i) { case 0: /* acos */ xp->retval = xp->arg1 >= 1.0 ? 0.0 : -PI; break; case 1: /* asin */ xp->retval = xp->arg1 >= 1.0 ? PI/2 : -PI/2; break; case 2: /* atan2 */ if (xp->arg1 == 0.0) xp->retval = xp->arg2 < 0.0 ? PI : 0.0; else if (xp->arg2 == 0.0) xp->retval = xp->arg1 < 0.0 ? -PI/2 : PI/2; else xp->retval = 0.0; break; case 3: /* pow */ /* FALLTHRU */ case 4: /* exp */ xp->retval = xp->o_type == OVERFLOW ? 1e308 : 0.0; break; case 5: /* log */ /* FALLTHRU */ case 6: /* log10 */ xp->retval = xp->arg1 <= 0.0 ? -1e308 : 0; break; case 7: /* sqrt */ xp->retval = 0.0; break; } return (1); /* suppress default error handling */ } #endif /* given the difference in two RA's, in rads, return their difference, * accounting for wrap at 2*PI. caller need *not* first force it into the * range 0..2*PI. */ double delra (double dra) { double fdra = fmod(fabs(dra), 2*PI); if (fdra > PI) fdra = 2*PI - fdra; return (fdra); } /* return 1 if object is considered to be "deep sky", else 0. * The only things deep-sky are fixed objects other than stars. */ int is_deepsky (Obj *op) { int deepsky = 0; if (is_type(op, FIXEDM)) { switch (op->f_class) { case 'T': case 'B': case 'D': case 'M': case 'S': case 'V': break; default: deepsky = 1; break; } } return (deepsky); }