mirror of https://github.com/XEphem/XEphem.git
502 lines
11 KiB
C
502 lines
11 KiB
C
/* misc handy functions.
|
|
* every system has such, no?
|
|
* 4/20/98 now_lst() always just returns apparent time
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#include "astro.h"
|
|
|
|
/* 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 < NBCM; i++)
|
|
if (binary_class_map[i].classcode == op->f_class)
|
|
return (binary_class_map[i].desc);
|
|
}
|
|
return ("Binary system");
|
|
case PLANET: {
|
|
static char nsstr[16];
|
|
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)
|
|
{
|
|
return (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);
|
|
}
|
|
|