mirror of https://github.com/XEphem/XEphem.git
401 lines
9.0 KiB
C
401 lines
9.0 KiB
C
/* Fetch the desired region of xe2-formatted catalog stars.
|
|
* This format supports proper motion. The format is:
|
|
*
|
|
* Note 1: all multibyte quantities are big-endian, aka, network byte order
|
|
* Note 2: header consists of ASCII version number followed by newline.
|
|
* Byte Bits Description
|
|
* 0 2..0 Name code:
|
|
* 0: PPM
|
|
* 1: SAO
|
|
* 2: HD
|
|
* 3: Tycho 2
|
|
* 4: HIP
|
|
* 5: GSC22
|
|
* 6..7: reserved
|
|
* 0 5..3 Type code:
|
|
* 0: star
|
|
* 1: star-like
|
|
* 2: double
|
|
* 3..7: reserved
|
|
* 0 7..6 Reserved
|
|
* 1..4 Name value is just the number except for Tycho 2:
|
|
* TYC1 = (B[3]<<6)|(B[4]>>2)
|
|
* TYC2 = (B[1]<<8)|(B[2])
|
|
* TYC3 = (B[4]&3)
|
|
* 5..7 RA, J2000, rads, 0 .. 2*PI mapped to 0 .. 1<<24
|
|
* 8..10 Dec, J2000, rads, -PI/2 .. PI/2 mapped to 0 .. (1<<24)-1
|
|
* 11 Magnitude, -2..17 mapped to 0..255
|
|
* 12..13 Spectral class, ASCII characters, left justified, blank filled
|
|
* 14..15 Proper motion in RA on sky, 2*(mas/yr) + 32768
|
|
* 16..17 Proper motion in Dec, 2*(mas/yr) + 32768
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <string.h>
|
|
|
|
#include "xephem.h"
|
|
|
|
|
|
#define XE2PKTSZ 18 /* packet size - DO NOT CHANGE */
|
|
|
|
typedef unsigned char UC;
|
|
typedef unsigned long UL;
|
|
typedef UC PKT[XE2PKTSZ];
|
|
|
|
|
|
/* an array of ObjF which can be grown efficiently in mults of NINC */
|
|
typedef struct {
|
|
ObjF *mem; /* malloced array */
|
|
int max; /* max number of cells in mem[] or -1 to just count */
|
|
int used; /* number actually in use */
|
|
} ObjFArray;
|
|
|
|
#define NINC 32 /* grow mem mem this many at a time */
|
|
|
|
static FILE *xe2open (char *file, char msg[]);
|
|
static void unpackObj (PKT pkt, Now *np, Obj *op);
|
|
static int addOneObjF (ObjFArray *ap, ObjF *fop);
|
|
static void unpack (PKT pkt, double *ra, double *dec, double *pma,
|
|
double *pmd, double *mag, char name[MAXNM], char *spect, int *type);
|
|
|
|
/* return 0 if file looks like a xe2, else -1
|
|
*/
|
|
int
|
|
xe2chkfile (char *file, char *msg)
|
|
{
|
|
FILE *fp;
|
|
|
|
fp = xe2open (file, msg);
|
|
if (!fp)
|
|
return (-1);
|
|
fclose (fp);
|
|
return (0);
|
|
}
|
|
|
|
/* malloc a collection of XE2 stars around the given location at *opp.
|
|
* if trouble fill msg[] and return -1, else return count.
|
|
* N.B. memory is only actually malloced at *opp if we return > 0.
|
|
* N.B. we very much rely on file being sorted by increasing dec.
|
|
* TODO: proper motion. hard because sort is only valid at one epoch.
|
|
*/
|
|
int
|
|
xe2fetch (char *file, Now *np, double ra, double dec, double fov, double mag,
|
|
ObjF **opp, char *msg)
|
|
{
|
|
double rov;
|
|
double cdec;
|
|
double sdec;
|
|
double crov;
|
|
ObjFArray ofa;
|
|
FILE *fp;
|
|
long l0, l, u, m;
|
|
double ldec, udec;
|
|
Obj o;
|
|
PKT pkt;
|
|
|
|
/* initial setup */
|
|
rov = fov/2;
|
|
cdec = cos(dec);
|
|
sdec = sin(dec);
|
|
crov = cos(rov);
|
|
|
|
/* bug if not given a place to save new stars */
|
|
if (!opp) {
|
|
printf ("xe2fetch with opp == NULL\n");
|
|
abort();
|
|
}
|
|
|
|
/* open the XE2 file */
|
|
fp = xe2open (file, msg);
|
|
if (!fp)
|
|
return (-1);
|
|
l0 = ftell(fp);
|
|
|
|
/* binary search for star nearest lower dec */
|
|
ldec = dec - rov;
|
|
if (ldec < -PI/2)
|
|
ldec = -PI/2;
|
|
l = 0;
|
|
fseek (fp, 0, SEEK_END);
|
|
u = (ftell(fp)-l0) / XE2PKTSZ - 1;
|
|
while (l <= u) {
|
|
m = (l+u)/2;
|
|
if (fseek (fp, m*XE2PKTSZ+l0, SEEK_SET) < 0) {
|
|
(void) sprintf (msg, "%s:\nseek error: %s", file, syserrstr());
|
|
fclose (fp);
|
|
return (-1);
|
|
}
|
|
if (fread (pkt, XE2PKTSZ, 1, fp) != 1) {
|
|
(void) sprintf (msg, "%s:\nread error: %s", file, syserrstr());
|
|
fclose (fp);
|
|
return (-1);
|
|
}
|
|
unpackObj (pkt, np, &o);
|
|
if (ldec < o.f_dec)
|
|
u = m-1;
|
|
else
|
|
l = m+1;
|
|
}
|
|
|
|
/* add each entry from m up to upper dec */
|
|
udec = dec + rov;
|
|
if (udec > PI/2)
|
|
udec = PI/2;
|
|
ofa.mem = NULL;
|
|
ofa.max = 0;
|
|
ofa.used = 0;
|
|
while (o.f_dec <= udec) {
|
|
if (get_mag(&o) <= mag && sin(o.f_dec)*sdec +
|
|
cos(o.f_dec)*cdec*cos(ra-o.f_RA) >= crov) {
|
|
if (addOneObjF (&ofa, (ObjF *)&o) < 0) {
|
|
(void) sprintf (msg, "%s:\nno more memory", file);
|
|
fclose (fp);
|
|
if (ofa.mem)
|
|
free ((void *)ofa.mem);
|
|
return (-1);
|
|
}
|
|
}
|
|
|
|
if (fread (pkt, XE2PKTSZ, 1, fp) != 1) {
|
|
if (feof(fp))
|
|
break;
|
|
(void) sprintf (msg, "%s:\nread error: %s", file, syserrstr());
|
|
fclose (fp);
|
|
if (ofa.mem)
|
|
free ((void *)ofa.mem);
|
|
return (-1);
|
|
}
|
|
|
|
unpackObj (pkt, np, &o);
|
|
}
|
|
|
|
fclose (fp);
|
|
if (ofa.mem)
|
|
*opp = ofa.mem;
|
|
return (ofa.used);
|
|
}
|
|
|
|
/* open the xe2 file.
|
|
* return the FILE * positioned at first record if ok,
|
|
* else NULL with a reason in msg[].
|
|
*/
|
|
static FILE *
|
|
xe2open (file, msg)
|
|
char *file;
|
|
char msg[];
|
|
{
|
|
FILE *fp = fopenh (file, "rb");
|
|
char header[32];
|
|
|
|
if (!fp) {
|
|
(void) sprintf (msg, "%s:\n%s", file, syserrstr());
|
|
return (NULL);
|
|
}
|
|
if (fgets (header, sizeof(header), fp) == NULL) {
|
|
sprintf (msg, "%s:\nno header: %s", file, syserrstr());
|
|
fclose (fp);
|
|
return(NULL);
|
|
}
|
|
if (strncmp (header, "XE2.", 4)) {
|
|
sprintf (msg, "%s:\nnot an XE2 file", file);
|
|
fclose (fp);
|
|
return(NULL);
|
|
}
|
|
|
|
return (fp);
|
|
}
|
|
|
|
/* crack open pkt and fill in op */
|
|
static void
|
|
unpackObj (pkt, np, op)
|
|
PKT pkt;
|
|
Now *np;
|
|
Obj *op;
|
|
{
|
|
#define MASRAD(mas) degrad((mas)/(3600.*1000.))
|
|
double ra;
|
|
double dec;
|
|
double pma;
|
|
double pmd;
|
|
double mag;
|
|
char name[MAXNM];
|
|
char spect[32];
|
|
int type;
|
|
|
|
unpack (pkt, &ra, &dec, &pma, &pmd, &mag, name, spect, &type);
|
|
|
|
zero_mem ((void *)op, sizeof (*op));
|
|
|
|
(void) strncpy (op->o_name, name, MAXNM-1);
|
|
op->o_type = FIXED;
|
|
op->f_class = type == 2 ? 'D' : type == 1 ? 'T' : 'S';
|
|
op->f_spect[0] = spect[0] == ' ' ? '\0' : spect[0];
|
|
op->f_spect[1] = spect[1] == ' ' ? '\0' : spect[1];
|
|
op->f_pmdec = (float)(MASRAD(pmd)/365.24);
|
|
op->f_dec = (float)(dec);
|
|
op->f_pmRA = (float)(MASRAD(pma)/cos(op->f_dec)/365.24);
|
|
op->f_RA = (float)(ra);
|
|
op->f_epoch = (float)J2000;
|
|
set_fmag (op, mag);
|
|
#undef MASRAD
|
|
}
|
|
|
|
/* add one ObjF entry to ap[], growing if necessary.
|
|
* return 0 if ok else return -1
|
|
*/
|
|
static int
|
|
addOneObjF (ap, fop)
|
|
ObjFArray *ap;
|
|
ObjF *fop;
|
|
{
|
|
ObjF *newf;
|
|
|
|
if (ap->used >= ap->max) {
|
|
/* add room for NINC more */
|
|
char *newmem = ap->mem ? realloc ((void *)ap->mem,
|
|
(ap->max+NINC)*sizeof(ObjF))
|
|
: malloc (NINC*sizeof(ObjF));
|
|
if (!newmem)
|
|
return (-1);
|
|
ap->mem = (ObjF *)newmem;
|
|
ap->max += NINC;
|
|
}
|
|
|
|
newf = &ap->mem[ap->used++];
|
|
(void) memcpy ((void *)newf, (void *)fop, sizeof(ObjF));
|
|
|
|
return (0);
|
|
}
|
|
|
|
/* unpack the raw PKT into its basic parts.
|
|
*/
|
|
static void
|
|
unpack (pkt, ra, dec, pma, pmd, mag, name, spect, type)
|
|
PKT pkt;
|
|
double *ra;
|
|
double *dec;
|
|
double *pma;
|
|
double *pmd;
|
|
double *mag;
|
|
char name[MAXNM];
|
|
char *spect;
|
|
int *type;
|
|
{
|
|
UL t;
|
|
|
|
/* type code */
|
|
*type = (pkt[0]>>3) & 7;
|
|
|
|
/* RA, rads */
|
|
t = ((UL)pkt[5] << 16) | ((UL)pkt[6] << 8) | (UL)pkt[7];
|
|
*ra = 2*PI*t/(1L<<24);
|
|
|
|
/* Dec, rads */
|
|
t = ((UL)pkt[8] << 16) | ((UL)pkt[9] << 8) | (UL)pkt[10];
|
|
*dec = PI*t/((1L<<24)-1) - PI/2;
|
|
|
|
/* mag */
|
|
t = (UL)pkt[11];
|
|
*mag = (19./255.)*t - 2;
|
|
|
|
/* spect */
|
|
spect[0] = pkt[12];
|
|
spect[1] = pkt[13];
|
|
|
|
/* pm ra, mas/yr */
|
|
t = ((UL)pkt[14] << 8) | (UL)pkt[15];
|
|
*pma = (t-32768.)/2.;
|
|
|
|
/* pm dec, mas/yr */
|
|
t = ((UL)pkt[16] << 8) | (UL)pkt[17];
|
|
*pmd = (t-32768.)/2.;
|
|
|
|
/* name */
|
|
if ((pkt[0]&7) == 3) {
|
|
/* tycho format */
|
|
UL tyc1 = ((UL)pkt[3] << 6) | ((UL)pkt[4] >> 2);
|
|
UL tyc2 = ((UL)pkt[1] << 8) | ((UL)pkt[2]);
|
|
UL tyc3 = ((UL)pkt[4] & 3);
|
|
sprintf (name, "Tyc %ld-%ld-%ld", tyc1, tyc2, tyc3); /* MAXNM! */
|
|
} else {
|
|
/* just a number */
|
|
t = ((UL)pkt[1] << 24) | ((UL)pkt[2] << 16) | ((UL)pkt[3] << 8)
|
|
| (UL)pkt[4];
|
|
switch (pkt[0]&7) {
|
|
case 0: sprintf (name, "PPM %ld", t); break;
|
|
case 1: sprintf (name, "SAO %ld", t); break;
|
|
case 2: sprintf (name, "HD %ld", t); break;
|
|
case 4: sprintf (name, "HIP %ld", t); break;
|
|
case 5: sprintf (name, "GSC2201 %c %ld", *dec<0?'S':'N', t); break;
|
|
default: sprintf (name, "<Bogus>"); break;
|
|
}
|
|
}
|
|
}
|
|
|
|
#ifdef MAIN_TEST
|
|
/* dump entire named .xe2 in .edb format to stdout.
|
|
*/
|
|
|
|
#include <errno.h>
|
|
|
|
int
|
|
main (int ac, char *av[])
|
|
{
|
|
ObjF *ofp = NULL;
|
|
Now now, *np = &now;
|
|
char *file;
|
|
char msg[1024];
|
|
int i, n;
|
|
|
|
if (ac != 2) {
|
|
fprintf (stderr, "usage: %s x.xe2 > x.edb\n", av[0]);
|
|
abort();
|
|
}
|
|
file = av[1];
|
|
if (xe2chkfile (file, msg) < 0) {
|
|
fprintf (stderr, "%s\n", msg);
|
|
abort();
|
|
}
|
|
|
|
zero_mem (np, sizeof(*np));
|
|
np->n_mjd = J2000 + 3*365; /* ~ 1/1/2003 */
|
|
n = xe2fetch (file, np, 0.0, PI/2, PI, 100.0, &ofp, msg);
|
|
if (n < 0) {
|
|
fprintf (stderr, "%s\n", msg);
|
|
abort();
|
|
}
|
|
for (i = 0; i < n; i++) {
|
|
db_write_line ((Obj *)&ofp[i], msg);
|
|
printf ("%s\n", msg);
|
|
}
|
|
free (ofp);
|
|
ofp = NULL;
|
|
n = xe2fetch (file, np, 0.0, -PI/2, PI, 100.0, &ofp, msg);
|
|
if (n < 0) {
|
|
fprintf (stderr, "%s\n", msg);
|
|
abort();
|
|
}
|
|
for (i = 0; i < n; i++) {
|
|
db_write_line ((Obj *)&ofp[i], msg);
|
|
printf ("%s\n", msg);
|
|
}
|
|
return (0);
|
|
}
|
|
|
|
char *
|
|
syserrstr ()
|
|
{
|
|
return (strerror(errno));
|
|
}
|
|
|
|
FILE *
|
|
fopenh (char *name, char *how)
|
|
{
|
|
return (fopen (name, how));
|
|
}
|
|
|
|
#endif
|
|
|