Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 3338) +++ /trunk/Ohana/src/addstar/Makefile (revision 3339) @@ -19,13 +19,17 @@ ADDSTAR = \ -$(SRC)/addstar.$(ARCH).o $(SRC)/gcatalog.$(ARCH).o \ -$(SRC)/gregions.$(ARCH).o $(SRC)/gimages.$(ARCH).o \ -$(SRC)/wimage.$(ARCH).o $(SRC)/gstars.$(ARCH).o \ -$(SRC)/find_matches.$(ARCH).o $(SRC)/wcatalog.$(ARCH).o \ -$(SRC)/sort_lists.$(ARCH).o $(SRC)/ConfigInit.$(ARCH).o \ -$(SRC)/aregion.$(ARCH).o \ -$(SRC)/check_permissions.$(ARCH).o $(SRC)/calibrate.$(ARCH).o \ -$(SRC)/args.$(ARCH).o $(SRC)/parse_time.$(ARCH).o \ -$(SRC)/mkcatalog.$(ARCH).o +$(SRC)/addstar.$(ARCH).o $(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/gregion_image.$(ARCH).o $(SRC)/gregion_star.$(ARCH).o \ +$(SRC)/gimages.$(ARCH).o $(SRC)/wimage.$(ARCH).o \ +$(SRC)/gstars.$(ARCH).o $(SRC)/grefstars.$(ARCH).o \ +$(SRC)/find_matches.$(ARCH).o $(SRC)/find_matches_refstars.$(ARCH).o \ +$(SRC)/wcatalog.$(ARCH).o $(SRC)/conversions.$(ARCH).o \ +$(SRC)/sort_lists.$(ARCH).o $(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/aregion.$(ARCH).o $(SRC)/find_subset.$(ARCH).o \ +$(SRC)/load_subpix.$(ARCH).o $(SRC)/airmass.$(ARCH).o \ +$(SRC)/mkcatalog.$(ARCH).o $(SRC)/calibrate.$(ARCH).o \ +$(SRC)/args.$(ARCH).o $(SRC)/parse_time.$(ARCH).o \ +$(SRC)/match_refstars.$(ARCH).o \ +$(SRC)/check_permissions.$(ARCH).o OBJ = $(ADDSTAR) Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 3338) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 3339) @@ -24,4 +24,5 @@ double DEFAULT_RADIUS, NSIGMA, SNLIMIT; double ZeroPt; +double Latitude, SiderealTime; int VERBOSE; int SKIP_MISSED; @@ -40,4 +41,6 @@ char AirmassKeyword[64]; char CCDNumKeyword[64]; +char STKeyword[64]; +char SubpixDatafile[256]; int CalReference; @@ -47,11 +50,22 @@ int ONLY_IMAGES; int FORCE_READ; +int SKYPROBE; +int ADDREFS; +int REPLACE; +int ONLY_MATCH; +time_t TIMEREF; + +int main (int argc, char **argv); int Shutdown (); void TrapSignal (int sig); -int main (int argc, char **argv); int SetSignals (); +void help (); +int args (int argc, char **argv); +void SetProtect (int mode); + int gcatalog (Catalog *catalog); -GSCRegion *gregions (Image *image, int *Nregions); +void gregion_star (Stars *star, GSCRegion *region); +GSCRegion *gregion_image (Image *image, int *Nregions); Image *gimages (FILE *f, Image *image, int *Npimage); int edge_check (double *x1, double *y1, double *x2, double *y2); @@ -59,21 +73,49 @@ void wimage (FILE *f, int dbstate, Image *image, int Nstars); Stars *gstars (char *file, int *NSTARS, Image *image); +void find_matches_refstars (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog); void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage); int in_image (double r, double d, Image *image); void wcatalog (Catalog *catalog); +void ConfigInit (int *argc, char **argv); +void aregion (GSCRegion *region, FILE *f, double ra, double dec); +void wimage (FILE *f, int dbstate, Image *image, int Nstars); +double airmass (double ra, double dec, double st, double latitude); +void InitCalibration (); +void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y); +void FindCalibration (Image *image); +double get_subpix (double x, double y); +void load_subpix (); +double scat_subpix (double x, double y); +void mkcatalog (GSCRegion *region, Catalog *catalog); +Stars *grefstars (char *file, int *Nstars); +Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, int *NSTARS); +int match_refstars (Stars *stars, int Nstars); +int parse_time (Header *header); + +void check_permissions (char *basefile); +void make_backup (char *filename); +void uppercase (char *string); + void sort_lists (float *X, float *Y, int *S, int N); void fsort (float *X, int N); -void ConfigInit (int *argc, char **argv); +void fsort2 (float *X, float *Y, int N); + +char *sec_to_date (time_t second); +double sec_to_jd (time_t second); +int chk_time (char *line); +int dms_to_ddd (double *Value, char *string); int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string); -void aregion (GSCRegion *region, FILE *f, double ra, double dec); -void wimage (FILE *f, int dbstate, Image *image, int Nstars); -void check_permissions (char *basefile); -void make_backup (char *filename); -void InitCalibration (); -void SaveCalibration (float M, float dM, float Mr, float Mc, float A, int N, float ra, float dec, float x, float y); -void FindCalibration (Image *image); -void help (); -int args (int argc, char **argv); -int parse_time (Header *header); -void uppercase (char *string); -void mkcatalog (GSCRegion *region, Catalog *catalog); +int str_to_dtime (char *line, double *second); +int str_to_radec (double *ra, double *dec, char *str1, char *str2); +int str_to_time (char *line, time_t *second); +time_t jd_to_sec (double jd); +time_t date_to_sec (char *date); + +/** + there is an inconsistency to be resolved: fixed structures (like Image) + need a fixed bit-length time (e_time), but these functions all use the + UNIX time_t types, which may be 32 or 64 bits, depending on the machine. + This can be resolved by using time_t with these functions, but casting + between e_time and time_t when necessary (ie, cannot return data to an + e_time pointer from one of these functions) +**/ Index: /trunk/Ohana/src/addstar/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/addstar/src/ConfigInit.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/ConfigInit.c (revision 3339) @@ -41,4 +41,5 @@ ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + /* used by parse_time to find time-related keywords */ ScanConfig (config, "DATE-KEYWORD", "%s", 0, DateKeyword); ScanConfig (config, "DATE-MODE", "%s", 0, DateMode); @@ -51,4 +52,8 @@ ScanConfig (config, "CCDNUM-KEYWORD", "%s", 0, CCDNumKeyword); + ScanConfig (config, "ST-KEYWORD", "%s", 0, STKeyword); + ScanConfig (config, "OBSERVATORY-LATITUDE", "%lf", 0, &Latitude); + ScanConfig (config, "SUBPIX_DATAFILE", "%s", 0, SubpixDatafile); + if (!LoadPhotcodes (PhotCodeFile)) { fprintf (stderr, "error loading photcodes\n"); Index: /trunk/Ohana/src/addstar/src/RegionOps.c =================================================================== --- /trunk/Ohana/src/addstar/src/RegionOps.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/RegionOps.c (revision 3339) @@ -0,0 +1,242 @@ +# include "addstar.h" +# define NBANDS +double DecBands[] = { +-90.0, -82.5, -75.0, -67.5, -60.0, +-52.5, -45.0, -37.5, -30.0, -22.5, +-15.0, -07.5, +00.0, +07.5, +15.0, ++22.5, +30.0, +37.5, +45.0, +52.5, ++60.0, +67.5, +75.0, +82.5, +90.0 +}; + +char *DecNames[] = { +"S8230", "S7500", "S6730", "S6000", "S5230" +"S4500", "S3730", "S3000", "S2230", "S1500", +"S0730", "S0000", "N0000", "N0730", "N1500", +"N2230", "N3000", "N3730", "N4500", "N5230", +"N6000", "N6730", "N7500", "N8230", +}; + +GSCRegion *LoadRegions (int *nregions) { + + Header header; + Header theader; + Matrix matrix; + FTable table; + + char buffer[28800], temp[50], file[256]; + double RA0, RA1, DEC0, DEC1; + int i, N, NDec, Nx, Ny; + + FILE *f; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC region file %s\n", GSCFILE); + exit (1); + } + + /* init & load in table data */ + if (!fits_fread_header (f, &header)) escape (UNLOCK, "ERROR: can't read primary header"); + if (!fits_fread_matrix (f, &matrix, &header)) escape (UNLOCK, "ERROR: can't read primary matrix"); + if (!fits_fread_ftable (f, &table, "REGIONS")) escape (UNLOCK, "ERROR: can't read table"); + + fits_scan (table.header, "NAXIS1", "%d", 1, &Nx); + fits_scan (table.header, "NAXIS2", "%d", 1, &Ny); + + Nregions = 0; + ALLOCATE (regions, GSCRegion, Ny); + + /* pole region is artificial - not in table */ + sprintf (file, "%s/n8230/pole.cpt", CATDIR); + region[Nregions].DEC[0] = 86.25; + region[Nregions].DEC[1] = 93.75; + region[Nregions].RA[0] = 0.0; + region[Nregions].RA[1] = 360.0; + strcpy (region[Nregions].filename, file); + Nregions ++; + + for (i = 0; i < Ny; i++) { + strncpy (temp, &table.buffer[i*48], 48); + temp[49] = 0; + hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); + if (RA1 < RA0) RA1 += 360.0; + + if (DEC0 < DEC1) { + regions[Nregions].DEC[0] = DEC0; + regions[Nregions].DEC[1] = DEC1; + } else { + regions[Nregions].DEC[0] = DEC1; + regions[Nregions].DEC[1] = DEC0; + } + /* skip the pole entries */ + if (regions[Nregions].DEC[1] > 86.25) continue; + + regions[Nregions].RA[0] = RA0; + regions[Nregions].RA[1] = RA1; + + /** convert DEC0, DEC1 to NBigDec **/ + NDec = -1; + for (N = 0; N < NBANDS; N++) { + if (DEC0 >= DecBand[N]) { + NDec = N; + break; + } + } + if (NDec < 0) { + fprintf (stderr, "programming / table error\n"); + exit (0); + } + + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt", CATDIR, Dec2Sections[NBigDec], &temp[1]); + + strcpy (regions[Nregions].filename, file); + Nregions ++; + } + + sort_regions (regions, Nregions); + *nregions = Nregions; + return (regions); +} + +FindRegionByPoint (double ra, double dec, GSCRegion *gsc, int Ngsc) { + + Ns = FindRegionDecBandStart (gsc, Ngsc, dec); + + for (i = Ns; i < Ngsc; i++) { + if (gsc[i].RA[0] >= ra) continue; + if (gsc[i].RA[1] < ra) continue; + if (gsc[i].DEC[0] >= dec) continue; + if (gsc[i].DEC[1] < dec) continue; + return (value); + } + + error: + return (-1); +} + +int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi, N1, N2; + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[0] < dec) { + Nlo = N; + } else { + Nhi = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] < v <= v[Nhi] */ + + /* find the last entry < dec */ + for (N = N2; N >= N1; N--) { + if (region[N].DEC[0] < dec) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +int FindRegionDecBandStop (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi, N1, N2; + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[1] > dec) { + Nhi = N; + } else { + Nlo = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] <= v < v[Nhi] */ + + /* find the last entry < dec */ + for (N = N1; N < N2; N++) { + if (region[N].DEC[1] > dec) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +/* bracket search: + + Nlo T N Nhi + X X X X X X X V X X X X X X X X X X X X X X + Nhi + + Nlo N T Nhi + X X X X X X X X X X X X X X X V X X X X X X + Nlo + +*/ + +/**********/ +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) { + + int flag_d0, flag_d1, flag_h0, flag_h1; + double tmp; + + *d0 = *h0 = *d1 = *h1 = 0; + + flag_h0 = dparse (h0, 1, string); + flag_h1 = dparse (h1, 4, string); + flag_d0 = dparse (d0, 7, string); + flag_d1 = dparse (d1, 9, string); + *h0 *= flag_h0; + *h1 *= flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + dparse (&tmp, 2, string); + *h0 += tmp/60.0; + dparse (&tmp, 3, string); + *h0 += tmp/3600.0; + + dparse (&tmp, 5, string); + *h1 += tmp/60.0; + dparse (&tmp, 6, string); + *h1 += tmp/3600.0; + + dparse (&tmp, 8, string); + *d0 += tmp/60.0; + + dparse (&tmp, 10, string); + *d1 += tmp/60.0; + + *h0 *= 15*flag_h0; + *h1 *= 15*flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + return (TRUE); +} + + +/* + if (buffer[i*48 + 19] == ' ') continue; + strncpy (temp, &buffer[i*48 + 19], 20); temp[20] = 0; + hms_to_deg (&RA, &DEC, temp, ' ', 3); + + + dBigRA = 360.0 / (int)(0.5 + 48*cos(3.1415927*0.5*(BigDecBounds[NBigDec] + BigDecBounds[NBigDec + 1])/180.0)); + + NBigRA = ra / dBigRA; + + NBig = NBigRA; + for (i = 0; (i < 12) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + for (i = 13; (i < 24) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + + + fprintf (stderr, "%d %d %d %d %f %f -> %f %f %f\n", NBigDec, NBigRA, NBigRASections[NBigDec], NBig, ra, dec, BigDecBounds[NBigDec], BigDecBounds[NBigDec + 1], dBigRA); +*/ Index: /trunk/Ohana/src/addstar/src/SetSignals.c =================================================================== --- /trunk/Ohana/src/addstar/src/SetSignals.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/SetSignals.c (revision 3339) @@ -0,0 +1,118 @@ +# include "addstar.h" + +static int Protect = FALSE; +static int Trapped = FALSE; + +void TrapSignal (int sig) { + fprintf (stderr, "trapped signal %d, exiting gracefully\n", sig); + if (sig == 11) { + fprintf (stderr, "seg fault\n"); + exit (1); + } + if (Protect) { + Trapped = TRUE; + fprintf (stderr, "blocking until protected sections are clear\n"); + return; + } + Shutdown (); +} + +void SetProtect (int mode) { + Protect = mode; + if (Trapped && !Protect) Shutdown (); +} + +int SetSignals () { + + int i; + + /* disable almost all signal interrupts */ + for (i = 0; i < 36; i++) { + switch (i) { + /* can't redirect this signals */ + case SIGKILL: /* kill -9: cannot be caught or ignored */ + case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ + case SIGCHLD: /* child halted: ignore */ + case SIGPWR: /* power failure - why ignore this? */ + case SIGWINCH: /* window resized */ + case SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ + break; + + default: + signal (i, TrapSignal); + } + } + return (TRUE); +} +/* + + Signal Value Action Comment + ------------------------------------------------------------------------- + SIGHUP 1 A Hangup detected on controlling terminal + or death of controlling process + SIGINT 2 A Interrupt from keyboard + SIGQUIT 3 A Quit from keyboard + SIGILL 4 A Illegal Instruction + SIGABRT 6 C Abort signal from abort(3) + SIGFPE 8 C Floating point exception + SIGKILL 9 AEF Kill signal + SIGSEGV 11 C Invalid memory reference + SIGPIPE 13 A Broken pipe: write to pipe with no readers + SIGALRM 14 A Timer signal from alarm(2) + SIGTERM 15 A Termination signal + SIGUSR1 30,10,16 A User-defined signal 1 + SIGUSR2 31,12,17 A User-defined signal 2 + SIGCHLD 20,17,18 B Child stopped or terminated + SIGCONT 19,18,25 Continue if stopped + SIGSTOP 17,19,23 DEF Stop process + SIGTSTP 18,20,24 D Stop typed at tty + SIGTTIN 21,21,26 D tty input for background process + SIGTTOU 22,22,27 D tty output for background process + + Next various other signals. + + Signal Value Action Comment + --------------------------------------------------------------------- + SIGTRAP 5 CG Trace/breakpoint trap + SIGIOT 6 CG IOT trap. A synonym for SIGABRT + SIGEMT 7,-,7 G + SIGBUS 10,7,10 AG Bus error + SIGSYS 12,-,12 G Bad argument to routine (SVID) + SIGSTKFLT -,16,- AG Stack fault on coprocessor + SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) + SIGIO 23,29,22 AG I/O now possible (4.2 BSD) + SIGPOLL AG A synonym for SIGIO (System V) + SIGCLD -,-,18 G A synonym for SIGCHLD + SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) + SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) + SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) + SIGPROF 27,27,29 AG Profile alarm clock + SIGPWR 29,30,19 AG Power failure (System V) + SIGINFO 29,-,- G A synonym for SIGPWR + SIGLOST -,-,- AG File lock lost + SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) + SIGUNUSED -,31,- AG Unused signal + (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and + sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a + sparc.) + + The letters in the "Action" column have the following meanings: + + A Default action is to terminate the process. + + B Default action is to ignore the signal. + + C Default action is to dump core. + + D Default action is to stop the process. + + E Signal cannot be caught. + + F Signal cannot be ignored. + + G Not a POSIX.1 conformant signal. + +*/ Index: /trunk/Ohana/src/addstar/src/addstar-old.c =================================================================== --- /trunk/Ohana/src/addstar/src/addstar-old.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/addstar-old.c (revision 3339) @@ -0,0 +1,186 @@ +# include "addstar.h" + +int main (int argc, char **argv) { + + int i, mode, dbstate; + int Nstars, Nimage, Nregions, Nmissed; + Stars *stars; + Image image, *pimage; + GSCRegion *region; + Catalog catalog; + + SetSignals (); + args (argc, argv); + + /* load the input datafile */ + if (ADDREFS) { + /** fixed-format ASCII (R,D,M,dM) = (1,2,3,4) **/ + stars = grefstars (argv[1], &Nstars); + /** all-sky catalogs in idiosyncratic formats **/ + stars = grefcat (catstats, &Nstars); + } else { + stars = gstars (argv[1], &Nstars, &image); + } + if (Nstars == 0) { + if (VERBOSE) fprintf (stderr, "no stars in data file, skipping\n"); + fprintf (stderr, "SUCCESS\n"); + exit (0); + } + + if (ADDREFS) { + match_refstars (stars, Nstars); + goto escape; + } + + /* load images */ + pimage = NULL; + if (dbstate == LCK_EMPTY) { + Nimage = 0; + } else { + pimage = gimages (f, &image, &Nimage); + } + region = gregion_image (&image, &Nregions); + fprintf (stderr, "region %f %f %f %f\n", region[0].RA[0], region[0].RA[1], region[0].DEC[0], region[0].DEC[1]); + + for (i = 0; i < Nregions; i++) { + check_permissions (region[i].filename); + } + + if (CALIBRATE) { InitCalibration (); } + + for (i = 0; i < Nregions; i++) { + catalog.filename = region[i].filename; /* don't free region before catalog! */ + fprintf (stderr, "adding to %s\n", region[i].filename); + + switch (lock_catalog (&catalog, LCK_XCLD)) { + case 0: + /* this is a serious error because ImageCat was NOT locked */ + fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); + exit (1); + case 1: + gcatalog (&catalog); /* load from disk */ + break; + case 2: + mkcatalog (®ion[i], &catalog); /* fills in new header info */ + break; + } + find_matches (®ion[i], stars, Nstars, &catalog, &image, pimage, Nimage); + + /* protect wcatalog from interrupt signals */ + Protect = TRUE; + if (!DUMP_MATCHES && !ONLY_IMAGES) wcatalog (&catalog); + if (Trapped) Shutdown (); + Protect = FALSE; + unlock_catalog (&catalog); + } + + if (CALIBRATE) { FindCalibration (&image); } + if (DUMP_MATCHES) Shutdown (); + + for (Nmissed = i = 0; i < Nstars; i++) { + if (stars[i].found == -1) { + fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM); + Nmissed ++; + } + } + if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed); + + fprintf (stderr, "SUCCESS\n"); + exit (0); + +} + +int SetSignals () { + + int i; + + /* disable almost all signal interrupts */ + for (i = 0; i < 36; i++) { + switch (i) { + /* can't redirect this signals */ + case SIGKILL: /* kill -9: cannot be caught or ignored */ + case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ + case SIGCHLD: /* child halted: ignore */ + case SIGPWR: /* power failure - why ignore this? */ + case SIGWINCH: /* window resized */ + case SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ + break; + + default: + signal (i, TrapSignal); + } + } + return (TRUE); +} +/* + + Signal Value Action Comment + ------------------------------------------------------------------------- + SIGHUP 1 A Hangup detected on controlling terminal + or death of controlling process + SIGINT 2 A Interrupt from keyboard + SIGQUIT 3 A Quit from keyboard + SIGILL 4 A Illegal Instruction + SIGABRT 6 C Abort signal from abort(3) + SIGFPE 8 C Floating point exception + SIGKILL 9 AEF Kill signal + SIGSEGV 11 C Invalid memory reference + SIGPIPE 13 A Broken pipe: write to pipe with no readers + SIGALRM 14 A Timer signal from alarm(2) + SIGTERM 15 A Termination signal + SIGUSR1 30,10,16 A User-defined signal 1 + SIGUSR2 31,12,17 A User-defined signal 2 + SIGCHLD 20,17,18 B Child stopped or terminated + SIGCONT 19,18,25 Continue if stopped + SIGSTOP 17,19,23 DEF Stop process + SIGTSTP 18,20,24 D Stop typed at tty + SIGTTIN 21,21,26 D tty input for background process + SIGTTOU 22,22,27 D tty output for background process + + Next various other signals. + + Signal Value Action Comment + --------------------------------------------------------------------- + SIGTRAP 5 CG Trace/breakpoint trap + SIGIOT 6 CG IOT trap. A synonym for SIGABRT + SIGEMT 7,-,7 G + SIGBUS 10,7,10 AG Bus error + SIGSYS 12,-,12 G Bad argument to routine (SVID) + SIGSTKFLT -,16,- AG Stack fault on coprocessor + SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) + SIGIO 23,29,22 AG I/O now possible (4.2 BSD) + SIGPOLL AG A synonym for SIGIO (System V) + SIGCLD -,-,18 G A synonym for SIGCHLD + SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) + SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) + SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) + SIGPROF 27,27,29 AG Profile alarm clock + SIGPWR 29,30,19 AG Power failure (System V) + SIGINFO 29,-,- G A synonym for SIGPWR + SIGLOST -,-,- AG File lock lost + SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) + SIGUNUSED -,31,- AG Unused signal + (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and + sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a + sparc.) + + The letters in the "Action" column have the following meanings: + + A Default action is to terminate the process. + + B Default action is to ignore the signal. + + C Default action is to dump core. + + D Default action is to stop the process. + + E Signal cannot be caught. + + F Signal cannot be ignored. + + G Not a POSIX.1 conformant signal. + +*/ Index: /trunk/Ohana/src/addstar/src/addstar.c =================================================================== --- /trunk/Ohana/src/addstar/src/addstar.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/addstar.c (revision 3339) @@ -1,32 +1,3 @@ # include "addstar.h" - -/* these variables are needed by Shutdown */ -static FILE *f = (FILE *) NULL; -static int Protect = FALSE; -static int Trapped = FALSE; - -/* clean up open / locked ImageCat before shutting down */ -int Shutdown () { - int dbstate; - - Protect = TRUE; - fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); - fprintf (stderr, "ERROR: addstar halted\n"); - exit (1); -} - -void TrapSignal (int sig) { - fprintf (stderr, "trapped signal %d, exiting gracefully\n", sig); - if (sig == 11) { - fprintf (stderr, "seg fault\n"); - exit (1); - } - if (Protect) { - Trapped = TRUE; - fprintf (stderr, "blocking until protected sections are clear\n"); - return; - } - Shutdown (); -} int main (int argc, char **argv) { @@ -34,68 +5,49 @@ int i, mode, dbstate; int Nstars, Nimage, Nregions, Nmissed; - Stars *stars, *gstars(); - Image image, *pimage, *gimages(); - GSCRegion *region, *gregions(); + Stars *stars; + Image image, *image_db; + GSCRegion *region; Catalog catalog; - + SetSignals (); + ConfigInit (&argc, argv); args (argc, argv); - /* load the input datafile */ - stars = gstars (argv[1], &Nstars, &image); - if (Nstars == 0) { - if (VERBOSE) fprintf (stderr, "no stars in data file, skipping\n"); - fprintf (stderr, "SUCCESS\n"); - exit (0); + lock_image_db (); + + switch (mode) { + case IMAGE: + stars = gstars (argv[1], &Nstars, &image); + region = gregion_image (&image, &Nregions); + image_db = gimages (&image, &Nimage_db); + break; + case REFLIST: + stars = grefstars (argv[1], &Nstars); + region = gregion_stars (stars, &Nregion); + break; + case PATCH : + regions = gregion_patch (patch, &Nregions); + break; + case MATCH: + regions = gregion_match (&Nregions); + break; } - /* lock the image catalog */ - check_permissions (ImageCat); - f = fsetlockfile (ImageCat, 3600.0, LCK_HARD, &dbstate); - if ((f == (FILE *) NULL) && (dbstate != LCK_EMPTY)) { - fprintf (stderr, "ERROR: can't lock image catalog\n"); - exit (1); - } - fseek (f, 0, SEEK_SET); + for (i = 0; i < Nregions; i++) { + load_catalog (&catalog, ®ion[i]); - /* load images */ - if (dbstate == LCK_EMPTY) { - Nimage = 0; - } else { - pimage = gimages (f, &image, &Nimage); - } - region = gregions (&image, &Nregions); - fprintf (stderr, "region %f %f %f %f\n", region[0].RA[0], region[0].RA[1], region[0].DEC[0], region[0].DEC[1]); - - for (i = 0; i < Nregions; i++) { - check_permissions (region[i].filename); - } - - if (CALIBRATE) { InitCalibration (); } - - for (i = 0; i < Nregions; i++) { - catalog.filename = region[i].filename; /* don't free region before catalog! */ - fprintf (stderr, "adding to %s\n", region[i].filename); - - switch (lock_catalog (&catalog, LCK_XCLD)) { - case 0: - /* this is a serious error because ImageCat was NOT locked */ - fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); - exit (1); - case 1: - gcatalog (&catalog); /* load from disk */ + switch (mode) { + case IMAGE: + find_matches (®ion[i], stars, Nstars, &catalog, &image, pimage, Nimage); break; - case 2: - mkcatalog (®ion[i], &catalog); /* fills in new header info */ + case PATCH: + case MATCH: + stars = grefcat (region, &Nstars); + case REFLIST: + subset = find_subset (®ion, stars, Nstars, &subref, &Nsubset); + find_matches_refstars (®ion, subset, Nsubset, &catalog); break; } - find_matches (®ion[i], stars, Nstars, &catalog, &image, pimage, Nimage); - - /* protect wcatalog from interrupt signals */ - Protect = TRUE; - if (!DUMP_MATCHES && !ONLY_IMAGES) wcatalog (&catalog); - if (Trapped) Shutdown (); - Protect = FALSE; - unlock_catalog (&catalog); + wcatalog (&catalog); } @@ -103,120 +55,14 @@ if (DUMP_MATCHES) Shutdown (); - for (Nmissed = i = 0; i < Nstars; i++) { - if (stars[i].found == -1) { - fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM); - Nmissed ++; - } - } - if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed); - - /* protect wimage from interrupt signals */ - Protect = TRUE; - wimage (f, dbstate, &image, Nstars); - if (Trapped) Shutdown (); - Protect = FALSE; - - fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); - - mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; - chmod (ImageCat, mode); - - fprintf (stderr, "SUCCESS\n"); + if (mode == IMAGE) wimage (f, dbstate, &image, Nstars); + unlock_image_db (); exit (0); - } -int SetSignals () { +/* names: - int i; - - /* disable almost all signal interrupts */ - for (i = 0; i < 36; i++) { - switch (i) { - /* can't redirect this signals */ - case SIGKILL: /* kill -9: cannot be caught or ignored */ - case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ - /* ignore these signals */ - case SIGCHLD: /* child halted: ignore */ - case SIGPWR: /* power failure - why ignore this? */ - case SIGWINCH: /* window resized */ - case SIGCONT: /* continue - maintain this action */ - case SIGTSTP: /* stop signal sent from tty - why ignore? */ - case SIGURG: /* socket signal, ignore this */ - break; - - default: - signal (i, TrapSignal); - } - } - return (TRUE); -} -/* - - Signal Value Action Comment - ------------------------------------------------------------------------- - SIGHUP 1 A Hangup detected on controlling terminal - or death of controlling process - SIGINT 2 A Interrupt from keyboard - SIGQUIT 3 A Quit from keyboard - SIGILL 4 A Illegal Instruction - SIGABRT 6 C Abort signal from abort(3) - SIGFPE 8 C Floating point exception - SIGKILL 9 AEF Kill signal - SIGSEGV 11 C Invalid memory reference - SIGPIPE 13 A Broken pipe: write to pipe with no readers - SIGALRM 14 A Timer signal from alarm(2) - SIGTERM 15 A Termination signal - SIGUSR1 30,10,16 A User-defined signal 1 - SIGUSR2 31,12,17 A User-defined signal 2 - SIGCHLD 20,17,18 B Child stopped or terminated - SIGCONT 19,18,25 Continue if stopped - SIGSTOP 17,19,23 DEF Stop process - SIGTSTP 18,20,24 D Stop typed at tty - SIGTTIN 21,21,26 D tty input for background process - SIGTTOU 22,22,27 D tty output for background process - - Next various other signals. - - Signal Value Action Comment - --------------------------------------------------------------------- - SIGTRAP 5 CG Trace/breakpoint trap - SIGIOT 6 CG IOT trap. A synonym for SIGABRT - SIGEMT 7,-,7 G - SIGBUS 10,7,10 AG Bus error - SIGSYS 12,-,12 G Bad argument to routine (SVID) - SIGSTKFLT -,16,- AG Stack fault on coprocessor - SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) - SIGIO 23,29,22 AG I/O now possible (4.2 BSD) - SIGPOLL AG A synonym for SIGIO (System V) - SIGCLD -,-,18 G A synonym for SIGCHLD - SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) - SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) - SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) - SIGPROF 27,27,29 AG Profile alarm clock - SIGPWR 29,30,19 AG Power failure (System V) - SIGINFO 29,-,- G A synonym for SIGPWR - SIGLOST -,-,- AG File lock lost - SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) - SIGUNUSED -,31,- AG Unused signal - (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and - sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a - sparc.) - - The letters in the "Action" column have the following meanings: - - A Default action is to terminate the process. - - B Default action is to ignore the signal. - - C Default action is to dump core. - - D Default action is to stop the process. - - E Signal cannot be caught. - - F Signal cannot be ignored. - - G Not a POSIX.1 conformant signal. - + catalog - existing object db table + region - sky area which may or may contain data + patch - RA,DEC bounded portion of sky + */ Index: /trunk/Ohana/src/addstar/src/airmass.c =================================================================== --- /trunk/Ohana/src/addstar/src/airmass.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/airmass.c (revision 3339) @@ -0,0 +1,18 @@ +# include "addstar.h" + +double airmass (double ra, double dec, double st, double latitude) { + + double hour, cosz, secz; + + /*** make this optional? we may not have ST... ***/ + /* ra, dec, latitude in dec deg; st in dec hours */ + + /* hour : hour angle in degrees */ + hour = 15.0*st - ra; + cosz = sin (RAD_DEG*dec) * sin (RAD_DEG*latitude) + cos (RAD_DEG*dec) * cos (RAD_DEG*hour) * cos (RAD_DEG*latitude); + + secz = 1.0 / cosz; + + return (secz); + +} Index: /trunk/Ohana/src/addstar/src/aregion.c =================================================================== --- /trunk/Ohana/src/addstar/src/aregion.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/aregion.c (revision 3339) @@ -88,5 +88,5 @@ if (!done) { - fprintf (stderr, "error in search: %f %f\n", ra, dec); + fprintf (stderr, "ERROR in search: %f %f\n", ra, dec); exit (1); } Index: /trunk/Ohana/src/addstar/src/args.c =================================================================== --- /trunk/Ohana/src/addstar/src/args.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/args.c (revision 3339) @@ -15,5 +15,5 @@ int args (int argc, char **argv) { - int N; + int i, N; /* check for help request */ @@ -22,7 +22,4 @@ help (); } - - /* configuration info */ - ConfigInit (&argc, argv); /* check for command line options */ @@ -43,4 +40,9 @@ remove_argument (N, &argc, argv); } + ONLY_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-image"))) { + ONLY_IMAGES = TRUE; + remove_argument (N, &argc, argv); + } DUMP_MATCHES = FALSE; if ((N = get_argument (argc, argv, "-dump"))) { @@ -50,6 +52,6 @@ FORCE_READ = FALSE; if ((N = get_argument (argc, argv, "-force"))) { + FORCE_READ = TRUE; remove_argument (N, &argc, argv); - FORCE_READ = TRUE; } SKIP_MISSED = FALSE; @@ -58,8 +60,32 @@ remove_argument (N, &argc, argv); } - ONLY_IMAGES = FALSE; - if ((N = get_argument (argc, argv, "-image"))) { - ONLY_IMAGES = TRUE; + SKYPROBE = FALSE; + if ((N = get_argument (argc, argv, "-skyprobe"))) { + SKYPROBE = TRUE; remove_argument (N, &argc, argv); + } + ADDREFS = FALSE; + if ((N = get_argument (argc, argv, "-addrefs"))) { + ADDREFS = TRUE; + remove_argument (N, &argc, argv); + } + REPLACE = FALSE; + if ((i = get_argument (argc, argv, "-replace"))) { + REPLACE = TRUE; + remove_argument (i, &argc, argv); + } + TIMEREF = 0; + if ((i = get_argument (argc, argv, "-time"))) { + remove_argument (i, &argc, argv); + if (!str_to_time (argv[i], &TIMEREF)) { + fprintf (stderr, "syntax error in time\n"); + exit (1); + } + remove_argument (i, &argc, argv); + } + ONLY_MATCH = FALSE; + if ((i = get_argument (argc, argv, "-only-match"))) { + ONLY_MATCH = TRUE; + remove_argument (i, &argc, argv); } CALIBRATE = FALSE; @@ -85,4 +111,6 @@ } + if (SKYPROBE) load_subpix (); + if (argc != 2) { fprintf (stderr, "USAGE: addstar (filename) [-dump] [-image] [-v] [-p photcode] [-accept] [-cal code1 code2]\n"); @@ -92,2 +120,16 @@ } +/** addstar modes: + + addstar (image.smp) - add cmp/smp image data to db + addstar -ref (file.dat) (photcode) + addstar -cat (USNO/2MASS/GSC) -region (ra dec - ra dec) + + -replace : ref/cat - replace existing match (photcode/time) + -match : ref/cat - only add measures to existing averages + + ref types: + ASCII - RA,DEC,M,dM in a table + +**/ + Index: /trunk/Ohana/src/addstar/src/bracket.c =================================================================== --- /trunk/Ohana/src/addstar/src/bracket.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/bracket.c (revision 3339) @@ -0,0 +1,37 @@ +int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi, N1, N2; + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[0] < dec) { + Nlo = N; + } else { + Nhi = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] < v <= v[Nhi] */ + + /* find the first entry >= dec */ + for (N = N1; N < N2; N++) { + if (region[N].DEC[0] >= T) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +/* bracket search: + + Nlo T N Nhi + X X X X X X X V X X X X X X X X X X X X X X + Nhi + + Nlo N T Nhi + X X X X X X X X X X X X X X X V X X X X X X + Nlo + +*/ + Index: /trunk/Ohana/src/addstar/src/calibrate.c =================================================================== --- /trunk/Ohana/src/addstar/src/calibrate.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/calibrate.c (revision 3339) @@ -2,5 +2,5 @@ static int Ncal, NCAL; -static float *Mobs, *dMobs, *Mref, *Cref, *Airm; +static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm; static float *RA, *DEC, *X, *Y; int *Nstar; @@ -14,4 +14,5 @@ ALLOCATE (dMobs, float, NCAL); ALLOCATE (Mref, float, NCAL); + ALLOCATE (dMref, float, NCAL); ALLOCATE (Cref, float, NCAL); ALLOCATE (Airm, float, NCAL); @@ -23,9 +24,10 @@ } -void SaveCalibration (float M, float dM, float Mr, float Mc, float A, int N, float ra, float dec, float x, float y) { +void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y) { Mobs[Ncal] = M; dMobs[Ncal] = dM; Mref[Ncal] = Mr; + dMref[Ncal] = dMr; Cref[Ncal] = Mc; Airm[Ncal] = A; @@ -43,4 +45,5 @@ REALLOCATE (Mobs, float, NCAL); REALLOCATE (Mref, float, NCAL); + REALLOCATE (dMref, float, NCAL); REALLOCATE (Cref, float, NCAL); REALLOCATE (Airm, float, NCAL); @@ -56,6 +59,6 @@ int i, MaxN, *Nlist, Nkeep; - float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag; - float dAs, dMs, dYs; + float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag, *dDmag; + float dMo, dMr, Mw, Dmed, W1, W2, NSigma; /* reject multiple matched-stars */ @@ -80,4 +83,5 @@ /* accumulate delta mags */ ALLOCATE (Dmag, float, Ncal); + ALLOCATE (dDmag, float, Ncal); Nkeep = 0; Klam = thiscode[0].K; @@ -89,22 +93,81 @@ Mabs = 0.001*(Mobs[i] + Klam*(Airm[i] - 1000) + Clam + Xlam*(Mref[i] - Cref[i])) - ZeroPt; - /* subpix correction : make parameters global! */ - dAs = 0; - if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i]; - dYs = Y[i] - (int) (Y[i]) + 0.55; - dMs = dAs * sin (6.28*dYs); - /* Mabs -= dMs; */ + /* note: subpix correction is applied in gstars */ if (DUMP_MATCHES) - fprintf (stdout, "%d %6.3f %6.3f %6.3f %10.6f %10.6f %7.2f %7.2f\n", - i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i]); + fprintf (stdout, "%d %6.3f %6.3f %6.3f %10.6f %10.6f %7.2f %7.2f %7.2f %7.2f\n", + i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i], Airm[i], 0.001*dMobs[i]); /* skip stars brighter than 8.0 */ - /* if (Mabs < 8.0) continue; */ + if (SKYPROBE) { + if (Mabs > 9.0) continue; + if (Mabs < 5.0) continue; + } + dMr = MAX (0.005, 0.001*dMref[i]); + dMo = MAX (0.005, 0.001*dMobs[i]); + Dmag[Nkeep] = (Mabs - 0.001*Mref[i]); + dDmag[Nkeep] = (dMr*dMr + dMo*dMo); Nkeep ++; } + if (Nkeep < 5) { + fprintf (stderr, "too few stars\n"); + image[0].Mcal = 10000; + image[0].dMcal = 10000; + return; + } + fsort2 (Dmag, dDmag, Nkeep); + + /* take sort list of Dmag, find median */ + Dmed = Dmag[(int)(0.5*Nkeep)]; + + /* exclude points with abs(Dmag - Dmed) / dDmag > 2.5 */ + + /* accumulate delta mags (25% - 75% of clipped range) */ + W1 = 0.0; + W2 = 0.0; + M1 = 0.0; + M2 = 0.0; + N = 0.0; + for (i = 0; i < Nkeep; i++) { + NSigma = fabs (Dmag[i] - Dmed) / sqrt (dDmag[i]); + if (NSigma > 2.5) continue; + W1 += Dmag[i] / dDmag[i]; + W2 += 1 / dDmag[i]; + M1 += Dmag[i]; + M2 += SQ (Dmag[i]); + N += 1.0; + } + + if (N > 1) { + M1 = M1 / N; + M2 = sqrt (fabs(M2/N - M1*M1)); + Mw = W1 / W2; + fprintf (stderr, "N: %.0f, mean: %f, wt mean: %f, stdev: %f, precision: %f\n", N, M1, Mw, M2, M2 / sqrt (N)); + image[0].Mcal = 1000 * M1; + image[0].dMcal = 1000 * M2 / sqrt (N); + image[0].Mxxxx = N; + } else { + fprintf (stderr, "too few stars\n"); + image[0].Mcal = 10000; + image[0].dMcal = 10000; + image[0].Mxxxx = 0; + } +} + +# if (0) + +/* old correction - superceeded by subpix grid */ + + /* find lo (Dmed - 0.1) + hi (Dmed + 0.1) + lo = -1; hi = -1; + for (i = 0; (lo == -1) && (i < Nkeep); i++) { if (Dmed - Dmag[i] < 0.1) lo = i; } + for (i = Nkeep - 1; (hi == -1) && (i >= 0); i--) { if (Dmag[i] - Dmed < 0.1) hi = i; } + if (lo == -1) lo = 0; + if (hi == -1) hi = Nkeep - 1; + */ + if (Nkeep < 3) { fprintf (stderr, "too few stars\n"); @@ -113,5 +176,5 @@ return; } - fsort (Dmag, Nkeep); + fsort2 (Dmag, dDmag, Nkeep); /* accumulate delta mags (25% - 75% of list) */ @@ -139,5 +202,12 @@ image[0].dMcal = 10000; } -} - - + + + /* subpix correction : make parameters global! */ + dAs = 0; + if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i]; + dYs = Y[i] - (int) (Y[i]) + 0.55; + dMs = dAs * sin (6.28*dYs); + /* Mabs -= dMs; */ + +# endif Index: /trunk/Ohana/src/addstar/src/conversions.c =================================================================== --- /trunk/Ohana/src/addstar/src/conversions.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/conversions.c (revision 3339) @@ -0,0 +1,299 @@ +# include "addstar.h" +enum {TIME_NONE, TIME_DATE, TIME_HOURS, TIME_DAYS, TIME_JD, TIME_MJD, TIME_SEC}; + +/***** convert [-]00:00:00 to 0.0000 ****/ +int dms_to_ddd (double *Value, char *string) { + + int valid, neg, status; + double tmp, value; + char *p1, *p2, *px; + + valid = FALSE; + neg = FALSE; + stripwhite (string); + p1 = string; + px = string + strlen(string); + + if (string[0] == '-') { + valid = TRUE; + neg = TRUE; + p1 = &string[1]; + } + if (string[0] == '+') { + valid = TRUE; + neg = FALSE; + p1 = &string[1]; + } + if (isdigit(string[0])) { + valid = TRUE; + p1 = &string[0]; + } + if (!valid) { return (FALSE); } + + status = 1; + tmp = strtod (p1, &p2); + if (p2 == p1) return (FALSE); /* entry not a number: +fred */ + value = tmp; + if (p2 == px) goto escape; /* entry only number: +1.0 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:fred */ + status = 2; + value += tmp / 60.0; + if (p2 == px) goto escape; /* entry only number: +1:1 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */ + value += tmp / 3600.0; + + escape: + if (neg) { + value *= -1; + } + *Value = value; + + return (status); +} + +/**********/ +int str_to_radec (double *ra, double *dec, char *str1, char *str2) { + + double Ra, Dec; + + *ra = *dec = 0; + switch (dms_to_ddd (&Ra, str1)) { + case 0: + fprintf (stderr, "syntax error in RA\n"); + return (FALSE); + case 1: + break; + case 2: + Ra = Ra * 15; + break; + } + switch (dms_to_ddd (&Dec, str2)) { + case 0: + fprintf (stderr, "syntax error in DEC\n"); + return (FALSE); + case 1: + case 2: + break; + } + *ra = Ra; + *dec = Dec; + return (TRUE); +} + +# define FORMAT_DAYS 1 +# define FORMAT_HOURS 2 +# define FORMAT_MINUTES 3 +# define FORMAT_SECONDS 4 +# define FORMAT_JD 5 +# define FORMAT_DATE 6 + +/**********/ +int chk_time (char *line) { + + char *p1, *p2; + double tmp; + int mode; + + p1 = line; + tmp = strtod (p1, &p2); + mode = FORMAT_DATE; + if (p2 == p1 + strlen (p1) - 1) { + if (*p2 == 'd') { + mode = FORMAT_DAYS; + } + if (*p2 == 'h') { + mode = FORMAT_HOURS; + } + if (*p2 == 'm') { + mode = FORMAT_MINUTES; + } + if (*p2 == 's') { + mode = FORMAT_SECONDS; + } + if (*p2 == 'j') { + mode = FORMAT_JD; + } + } + return (mode); +} + +/**********/ +int str_to_time (char *line, time_t *second) { + + char *tmpline; + struct tm *gmt; + struct timeval now; + double jd; + time_t tsec; + + if (!strcasecmp (line, "NOW")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = now.tv_sec; + return (TRUE); + } + + if (!strncasecmp (line, "TODAY", 5)) { + gettimeofday (&now, (struct timezone *) NULL); + if (line[5]) { /* line has extra data (ie, hh:mm:ss) */ + tsec = now.tv_sec; + ALLOCATE (tmpline, char, 64); + gmt = gmtime (&tsec); + sprintf (tmpline, "%04d/%02d/%02d,%s", 1900 + gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, &line[6]); + *second = date_to_sec (tmpline); + free (tmpline); + return (TRUE); + } else { + *second = 86400 * ((int)(now.tv_sec / 86400)); + return (TRUE); + } + } + + switch (chk_time (line)) { + case 0: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + case FORMAT_JD: + jd = strtod (line, 0); + *second = jd_to_sec (jd); + return (TRUE); + case FORMAT_DATE: + *second = date_to_sec (line); + return (TRUE); + } + return (FALSE); +} + + +/**********/ +int str_to_dtime (char *line, double *second) { + + switch (chk_time (line)) { + case 0: + case FORMAT_JD: + case FORMAT_DATE: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + } + return (FALSE); +} + +/**********/ +double sec_to_jd (time_t second) { + + double jd; + + jd = second/86400.0 + 2440587.5; + return (jd); +} + +/**********/ +time_t jd_to_sec (double jd) { + + time_t second; + + second = (jd - 2440587.5)*86400; + return (second); +} + +/**********/ +char *sec_to_date (time_t second) { + + struct tm *gmt; + char *line; + + ALLOCATE (line, char, 64); + gmt = gmtime (&second); + sprintf (line, "%04d/%02d/%02d,%02d:%02d:%02d", + 1900 + gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, + gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec); + return (line); + +} + +/***** date in format yyyy/mm/dd,hh:mm:ss *****/ +time_t date_to_sec (char *date) { + + time_t second; + double tmp, jd; + struct tm now; + char *p1, *p2, *px; + + p1 = date; + px = date + strlen(date); + bzero (&now, sizeof(now)); + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_year = tmp; + if (now.tm_year > 1000) now.tm_year -= 1900; + if (now.tm_year < 50) now.tm_year += 100; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mon = tmp - 1; /* mon runs from 0 - 11 */ + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mday = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + p1 = p2 + 1; + now.tm_hour = tmp; + if (p2 == px) goto escape; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_min = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_sec = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + escape: + jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4) + + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12) + - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5; + + second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec; + + return (second); +} Index: /trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/find_matches.c (revision 3339) @@ -7,7 +7,8 @@ void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage) { + double secz; int i, j, k, n, m, N, first_j, found, found0, found1; double X, Y, RADIUS, RADIUS2; - float *X1, *Y1, *X2, *Y2, CalM, CalC, *Xs, *Ys; + float *X1, *Y1, *X2, *Y2, CalM, CalC, dCalM, *Xs, *Ys; float dX, dY, dR; int *N1, *N2, *next, *next_miss, last, last_miss; @@ -19,18 +20,20 @@ short int Mrel, *Mval; - /** allocate local arrays **/ + /* photcode data */ + CalM = CalC = dCalM = 0; + Nsecfilt = GetPhotcodeNsecfilt (); + Nsec = GetPhotcodeNsec (thiscode[0].code); + + /** allocate local arrays (stars) **/ ALLOCATE (X1, float, Nstars); ALLOCATE (Y1, float, Nstars); - ALLOCATE (N1, int, Nstars); - - /** allocate local arrays **/ - - Nsecfilt = GetPhotcodeNsecfilt (); - + ALLOCATE (N1, int, Nstars); + + /** allocate local arrays (catalog) **/ Nave = catalog[0].Naverage; NAVE = Nave + 1000; ALLOCATE (X2, float, NAVE); ALLOCATE (Y2, float, NAVE); - ALLOCATE (N2, int, NAVE); + ALLOCATE (N2, int, NAVE); ALLOCATE (Xs, float, NAVE); ALLOCATE (Ys, float, NAVE); @@ -51,7 +54,4 @@ REALLOCATE (catalog[0].missing, Missing, NMISS); - /* check on photcode */ - Nsec = GetPhotcodeNsec (thiscode[0].code); - /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ /* reference for coords is this image */ @@ -63,5 +63,5 @@ for (i = 0; i < Nstars; i++) { - fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); + fRD_to_XY (&X1[i], &Y1[i], stars[0][i].R, stars[0][i].D, &tcoords); N1[i] = i; } @@ -90,11 +90,11 @@ /* choose a radius for matches */ - /* 0.02 corrects cerror to arcsec from storage units */ if (DEFAULT_RADIUS == 0) { - RADIUS = NSIGMA * 0.02 * image[0].cerror; + RADIUS = NSIGMA * 0.02 * image[0].cerror; /* 0.02 corrects cerror to arcsec from storage units */ } else { RADIUS = DEFAULT_RADIUS; } RADIUS2 = RADIUS*RADIUS; + /* correct instrumental mags for exposure time */ if (image[0].exptime > 0) { @@ -108,5 +108,4 @@ dX = X1[i] - X2[j]; - /* fprintf (stderr, "%f %f %f %f %f %f %f %f %f\n", dX, X1[i], Y1[i], X2[j], Y2[j], stars[i].R, stars[i].D, catalog[0].average[j].R, catalog[0].average[j].D); */ if (dX <= -2*RADIUS) { @@ -127,16 +126,22 @@ if (dR < RADIUS2) { /* new measurement of this star */ Nmatch ++; - /** insert star in measurement list */ n = N2[j]; N = N1[i]; - m = catalog[0].average[n].offset; /* first measurement of this star */ - for (k = 0; k < catalog[0].average[n].Nm - 1; k++) - m = next[m]; + + /** insert star in measurement list */ + /* find last measurement of this star */ + m = catalog[0].average[n].offset; + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; + /* set up references */ next[Nmeas] = next[m]; next[m] = Nmeas; - if (next[Nmeas] == -1) { /* last just was moved */ - last = Nmeas; - } + /* last just was moved */ + if (next[Nmeas] == -1) last = Nmeas; + /* calculate accurate per-star airmass */ + secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude); + + /* find ref measurements & add to calibration table */ + /* is this search right? it does not use the next[] list */ if (CALIBRATE) { found = FALSE; @@ -145,13 +150,20 @@ m = catalog[0].average[n].offset; /* first measurement of this star */ for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) { - if (catalog[0].measure[m+k].source == CalReference) { found0 = TRUE; CalM = catalog[0].measure[m+k].M; } - if (catalog[0].measure[m+k].source == CalColor) { found1 = TRUE; CalC = catalog[0].measure[m+k].M; } + if (catalog[0].measure[m+k].source == CalReference) { + found0 = TRUE; + CalM = catalog[0].measure[m+k].M; + dCalM = catalog[0].measure[m+k].dM; + } + if (catalog[0].measure[m+k].source == CalColor) { + found1 = TRUE; + CalC = catalog[0].measure[m+k].M; + } if (found0 && found1) { found = TRUE; - SaveCalibration (MIN (NO_MAG, 1000*stars[N].M + MTIME), - MIN (255, stars[N].dM), - CalM, CalC, - image[0].secz, N, - stars[N].R, stars[N].D, Xs[n], Ys[n]); + SaveCalibration (MIN (NO_MAG, 1000*stars[0][N].M + MTIME), + MIN (255, stars[0][N].dM), + CalM, dCalM, CalC, + secz, N, + stars[0][N].R, stars[0][N].D, Xs[n], Ys[n]); } } @@ -159,21 +171,21 @@ /** add measurements for this star **/ - catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); - catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); - catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000*stars[N].M + MTIME); + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[0][N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[0][N].D); + catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000*stars[0][N].M + MTIME); catalog[0].measure[Nmeas].Mcal = image[0].Mcal; - catalog[0].measure[Nmeas].dM = MIN (255, stars[N].dM); /* error in input files stored in thousandths of mag */ - catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */ + catalog[0].measure[Nmeas].dM = MIN (255, stars[0][N].dM); /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate; /* trate is in 0.1 msec / row */ catalog[0].measure[Nmeas].averef = n; catalog[0].measure[Nmeas].source = image[0].source; /* photometry source */ - catalog[0].measure[Nmeas].dophot = stars[N].dophot; + catalog[0].measure[Nmeas].dophot = stars[0][N].dophot; catalog[0].measure[Nmeas].flags = 0; catalog[0].measure[Nmeas].dt = MTIME; - catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000*stars[N].Mgal + MTIME); - catalog[0].measure[Nmeas].airmass = image[0].secz; - catalog[0].measure[Nmeas].FWx = MIN (NO_MAG, 100*stars[N].fx); - catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[N].fy / stars[N].fx)); - catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[N].df); + catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000*stars[0][N].Mgal + MTIME); + catalog[0].measure[Nmeas].airmass = secz; + catalog[0].measure[Nmeas].FWx = MIN (NO_MAG, 100*stars[0][N].fx); + catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx)); + catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[0][N].df); /* refers to same number as first measurement */ @@ -185,13 +197,13 @@ /*** handle multiple stars */ /* this image star matches more than one catalog star */ - if (stars[N].found > -1) { - catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE; + if (stars[0][N].found > -1) { + catalog[0].measure[stars[0][N].found].flags |= BLEND_IMAGE; catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; } - if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ + if (stars[0][N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; } - if (stars[N].found == -1) { /* this image star matches only this catalog star */ - stars[N].found = Nmeas; /* save first match, in case coincidences are found */ + if (stars[0][N].found == -1) { /* this image star matches only this catalog star */ + stars[0][N].found = Nmeas; /* save first match, in case coincidences are found */ } /* this catalog star matches more than one image star */ @@ -291,7 +303,7 @@ for (i = 0; i < Nstars; i++) { N = N1[i]; - if ((stars[N].found < 0) && IN_CATALOG (stars[N].R, stars[N].D)) { - catalog[0].average[Nave].R = stars[N].R; - catalog[0].average[Nave].D = stars[N].D; + if ((stars[0][N].found < 0) && IN_CATALOG (stars[0][N].R, stars[0][N].D)) { + catalog[0].average[Nave].R = stars[0][N].R; + catalog[0].average[Nave].D = stars[0][N].D; catalog[0].average[Nave].M = NO_MAG; for (j = 0; j < Nsecfilt; j++) { @@ -300,7 +312,8 @@ } - Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K; - Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M; - if (*Mval == NO_MAG) *Mval = Mrel; + if (SKYPROBE) + secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude); + else + secz = image[0].secz; catalog[0].average[Nave].Nm = 1; @@ -315,20 +328,25 @@ catalog[0].measure[Nmeas].dR = 0.0; catalog[0].measure[Nmeas].dD = 0.0; - catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000.0*stars[N].M + MTIME); + catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000.0*stars[0][N].M + MTIME); catalog[0].measure[Nmeas].Mcal = image[0].Mcal; - catalog[0].measure[Nmeas].dM = MIN (255, stars[N].dM); - catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */ + catalog[0].measure[Nmeas].dM = MIN (255, stars[0][N].dM); + catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate; /* trate is in 0.1 msec / row */ catalog[0].measure[Nmeas].averef = Nave; catalog[0].measure[Nmeas].source = image[0].source; /* photometry source */ - catalog[0].measure[Nmeas].dophot = stars[N].dophot; + catalog[0].measure[Nmeas].dophot = stars[0][N].dophot; catalog[0].measure[Nmeas].flags = 0; catalog[0].measure[Nmeas].dt = MTIME; - catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000.0*stars[N].Mgal + MTIME); - catalog[0].measure[Nmeas].airmass = image[0].secz; - catalog[0].measure[Nmeas].FWx = MIN (0x7fff, 100*stars[N].fx); - catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[N].fy / stars[N].fx)); - catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[N].df); - stars[N].found = Nmeas; + catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000.0*stars[0][N].Mgal + MTIME); + catalog[0].measure[Nmeas].airmass = secz; + catalog[0].measure[Nmeas].FWx = MIN (0x7fff, 100*stars[0][N].fx); + catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx)); + catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[0][N].df); + + Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K; + Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M; + if (*Mval == NO_MAG) *Mval = Mrel; + + stars[0][N].found = Nmeas; next[last] = Nmeas; next[Nmeas] = -1; @@ -340,4 +358,5 @@ REALLOCATE (catalog[0].measure, Measure, NMEAS); } + /** now add references from all previous non-detection observations of this spot on the sky */ for (j = 0; (j < Nimage) && !SKIP_MISSED; j++) { @@ -408,6 +427,6 @@ /* note stars which have been found in this catalog */ for (i = 0; i < Nstars; i++) { - if (stars[i].found > -1) { - stars[i].found = -2; + if (stars[0][i].found > -1) { + stars[0][i].found = -2; } } Index: /trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 3339) @@ -0,0 +1,387 @@ +# include "addstar.h" + +# define IN_CATALOG(R,D) ( \ +((D) >= region[0].DEC[0]) && ((D) < region[0].DEC[1]) && \ +((R) >= region[0].RA[0]) && ((R) < region[0].RA[1])) + +void find_matches_refstars (GSCRegion *region, Stars **stars, int Nstars, Catalog *catalog) { + + int i, j, k, n, m, N, first_j; + double RADIUS, RADIUS2; + float *X1, *Y1, *X2, *Y2; + float dX, dY, dR; + int *N1, *N2, *next, *next_miss, last, last_miss; + int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch; + Measure *tmpmeasure; + Missing *tmpmissing; + Coords tcoords; + int Nsecfilt, Nsec; + + /* photcode data */ + Nsecfilt = GetPhotcodeNsecfilt (); + Nsec = GetPhotcodeNsec (thiscode[0].code); + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, float, Nstars); + ALLOCATE (Y1, float, Nstars); + ALLOCATE (N1, int, Nstars); + + /** allocate local arrays (catalog) **/ + Nave = catalog[0].Naverage; + NAVE = Nave + 1000; + ALLOCATE (X2, float, NAVE); + ALLOCATE (Y2, float, NAVE); + ALLOCATE (N2, int, NAVE); + ALLOCATE (catalog[0].found, int, NAVE); + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); + /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */ + + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + + Nmiss = catalog[0].Nmissing; + NMISS = Nmiss + 1000; + ALLOCATE (next_miss, int, NMISS); + REALLOCATE (catalog[0].missing, Missing, NMISS); + + /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ + /* reference for coords is catalog center */ + tcoords.crval1 = 0.5*(region[0].RA[0] + region[0].RA[1]); + tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]); + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + tcoords.Npolyterms = 0; + + for (i = 0; i < Nstars; i++) { + fRD_to_XY (&X1[i], &Y1[i], stars[i][0].R, stars[i][0].D, &tcoords); + N1[i] = i; + } + if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars); + + for (i = 0; i < Nave; i++) { + fRD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + catalog[0].found[N2[i]] = -1; + } + if (Nave > 1) sort_lists (X2, Y2, N2, Nave); + + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + for (i = 0; i < Nmiss - 1; i++) { + next_miss[i] = i+1; + } + next_miss[i] = -1; + last_miss = i; + + /* choose a radius for matches */ + if (DEFAULT_RADIUS == 0) { + RADIUS = 2.0; + } else { + RADIUS = DEFAULT_RADIUS; + } + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + + dX = X1[i] - X2[j]; + + if (dX <= -2*RADIUS) { + i++; + continue; + } + if (dX >= 2*RADIUS) { + j++; + continue; + } + + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*RADIUS) && (j < Nave); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < RADIUS2) { /* new measurement of this star */ + Nmatch ++; + n = N2[j]; + N = N1[i]; + + /** in replace mode, search for entry and replace values M, dM, R, D */ + if (REPLACE) { + int found; + found = FALSE; + m = catalog[0].average[n].offset; + for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) { + if (catalog[0].measure[m+k].source == thiscode[0].code) { + catalog[0].measure[m+k].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); + catalog[0].measure[m+k].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); + catalog[0].measure[m+k].M = 1000.0*stars[N].M; + catalog[0].measure[m+k].dM = 1000.0*MIN (255, stars[N].dM); /* error in input files stored in thousandths of mag */ + found = TRUE; + stars[N].found = m+k; /* save first match, in case coincidences are found */ + } + } + if (found) continue; + } + + /** insert star in measurement list */ + /* find last measurement of this star */ + m = catalog[0].average[n].offset; + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; + /* set up references */ + next[Nmeas] = next[m]; + next[m] = Nmeas; + /* last just was moved */ + if (next[Nmeas] == -1) last = Nmeas; + + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); + catalog[0].measure[Nmeas].M = 1000.0*stars[N].M; + catalog[0].measure[Nmeas].Mcal = 0; + catalog[0].measure[Nmeas].dM = 1000.0*MIN (255, stars[N].dM); /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = TIMEREF; /** careful : time_t vs e_time **/ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = thiscode[0].code; + catalog[0].measure[Nmeas].dophot = stars[N].dophot; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].flags = 0; + catalog[0].measure[Nmeas].dt = 0xffff; + + catalog[0].measure[Nmeas].Mgal = NO_MAG; + catalog[0].measure[Nmeas].airmass = 0; + catalog[0].measure[Nmeas].FWx = NO_MAG; + catalog[0].measure[Nmeas].fwy = 0xff; + catalog[0].measure[Nmeas].theta = 0xff; + + /* check for entries in the secfilt lists */ + Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; + if (*Mval == NO_MAG) *Mval = catalog[0].measure[Nmeas].M; /*** no airmass correction needed ***/ + /*** the problem here is that the ref photcode should never equate to a pri/sec photcode, right? ***/ + /* just skip assigning an average magnitude for this object? */ + + /*** handle multiple stars */ + /* this image star matches more than one catalog star */ + if (stars[N].found > -1) { + catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE; + catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; + } + if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ + catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; + } + if (stars[N].found == -1) { /* this image star matches only this catalog star */ + stars[N].found = Nmeas; /* save first match, in case coincidences are found */ + } + /* this catalog star matches more than one image star */ + if (catalog[0].found[n] > -1) { + catalog[0].measure[catalog[0].found[n]].flags |= BLEND_CATALOG; + catalog[0].measure[Nmeas].flags |= BLEND_CATALOG; + } else { + catalog[0].found[n] = Nmeas; + } + + catalog[0].average[n].Nm ++; + Nmeas ++; + if (Nmeas == NMEAS) { + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + + /* update values of Ro, Do */ + if (catalog[0].average[n].Nm > 2) { + double R, D, r, d, r2, d2, Npt; + Npt = r = d = r2 = d2 = 0; + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm; k++) { + if (catalog[0].measure[m].t == 0) { + m = next[m]; + continue; + } + R = catalog[0].measure[m].dR; + D = catalog[0].measure[m].dD; + r += R; + d += D; + r2 += R*R; + d2 += D*D; + m = next[m]; + Npt += 1.0; + } + if (Npt > 2) { + r = r / Npt; /* these are corrections in 1/100 arcsec to RA and DEC */ + d = d / Npt; + R = r2 / Npt - r*r; + D = d2 / Npt - d*d; + /* Xp is scatter in position in hundredths of arcsec */ + catalog[0].average[n].Xp = sqrt (D + R / SQ(cos(catalog[0].average[n].D*RAD_DEG))); + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm; k++) { + catalog[0].measure[m].dR = catalog[0].measure[m].dR - r; + catalog[0].measure[m].dD = catalog[0].measure[m].dD - d; + m = next[m]; + } + catalog[0].average[n].R = catalog[0].average[n].R - r / 360000.0; + catalog[0].average[n].D = catalog[0].average[n].D - d / 360000.0; + } + } + } + } + j = first_j; + i++; + } + + /* incorporate unmatched image stars, if this star is in field of this catalog */ + /* skip if we want to force matches - combined with -replace, this lets us keep + the reference up-to-date with known stars only */ + for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { + N = N1[i]; + if (stars[N].found < 0) { + catalog[0].average[Nave].R = stars[N].R; + catalog[0].average[Nave].D = stars[N].D; + catalog[0].average[Nave].M = NO_MAG; + for (j = 0; j < Nsecfilt; j++) { + catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; + catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; + } + + if (thiscode[0].type == PHOT_PRI) { + catalog[0].average[Nave].M = 1000.0*stars[N].M; + } + if (thiscode[0].type == PHOT_SEC) { + catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = 1000.0*stars[N].M; + } + catalog[0].average[Nave].Nm = 1; + catalog[0].average[Nave].Nn = 0; + catalog[0].average[Nave].Xp = NO_MAG; + catalog[0].average[Nave].Xm = NO_MAG; + catalog[0].average[Nave].dM = NO_MAG; + catalog[0].average[Nave].offset = Nmeas; + catalog[0].average[Nave].missing = -1; + catalog[0].average[Nave].code = 0; + + catalog[0].measure[Nmeas].dR = 0.0; + catalog[0].measure[Nmeas].dD = 0.0; + catalog[0].measure[Nmeas].M = 1000.0*stars[N].M; + catalog[0].measure[Nmeas].Mcal = 0; + catalog[0].measure[Nmeas].dM = MIN (255, stars[N].dM); + catalog[0].measure[Nmeas].t = 0; + catalog[0].measure[Nmeas].averef = Nave; + catalog[0].measure[Nmeas].source = thiscode[0].code; + catalog[0].measure[Nmeas].dophot = stars[N].dophot; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].flags = 0; + catalog[0].measure[Nmeas].dt = 0xffff; + + catalog[0].measure[Nmeas].Mgal = NO_MAG; + catalog[0].measure[Nmeas].airmass = 0; + catalog[0].measure[Nmeas].FWx = NO_MAG; + catalog[0].measure[Nmeas].fwy = 0xff; + catalog[0].measure[Nmeas].theta = 0xff; + + stars[N].found = Nmeas; + next[last] = Nmeas; + next[Nmeas] = -1; + last = Nmeas; + Nmeas ++; + if (Nmeas == NMEAS) { + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + Nave ++; + if (Nave == NAVE) { + NAVE = Nave + 1000; + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); + } + } + } + + free (catalog[0].found); + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + tmpmeasure[N].averef = i; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + /* fix order of Missing (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmissing, Missing, Nmiss); + for (i = 0; i < Nave; i++) { + if (catalog[0].average[i].Nn > 0) { + n = catalog[0].average[i].missing; + catalog[0].average[i].missing = N; + for (k = 0; k < catalog[0].average[i].Nn; k++, N++) { + tmpmissing[N] = catalog[0].missing[n]; + n = next_miss[n]; + } + } + } + free (catalog[0].missing); + catalog[0].missing = tmpmissing; + + /* note stars which have been found in this catalog */ + for (i = 0; i < Nstars; i++) { + if (stars[i].found > -1) { + stars[i].found = -2; + } else { + stars[i].found = -3; + } + } + + catalog[0].Naverage = Nave; + catalog[0].Nmeasure = Nmeas; + catalog[0].Nmissing = Nmiss; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch); +} + + +# if (0) + /* this is a bit of an issue. i don't want to include missing obs if the + ref star has not been previously detected (why?), but if i do this, it won't + get all missed entries if ref star is detected in the future... */ + /** now add references from all previous non-detection observations of this spot on the sky */ + for (j = 0; j < Nimage; j++) { + if (in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j])) { + if (catalog[0].average[Nave].Nn < 1) { + catalog[0].average[Nave].missing = Nmiss; + } + next_miss[last_miss] = Nmiss; + next_miss[Nmiss] = -1; + last_miss = Nmiss; + /* this can now be done exactly */ + RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j].coords); + catalog[0].missing[Nmiss].t = pimage[j].tzero + 1e-4*Y*pimage[j].trate; /* rough guess at time */ + catalog[0].average[Nave].Nn ++; + Nmiss ++; + if (Nmiss == NMISS) { + NMISS = Nmiss + 1000; + REALLOCATE (next_miss, int, NMISS); + REALLOCATE (catalog[0].missing, Missing, NMISS); + } + + } + } +# endif Index: /trunk/Ohana/src/addstar/src/find_proper.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_proper.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/find_proper.c (revision 3339) @@ -0,0 +1,158 @@ +# include "addusno.h" + +find_proper (catstats, catalog, usnostats, usno, Nusno) +CatStats catstats[]; +Catalog catalog[]; +USNOstats usnostats[]; +USNOdata usno[]; +int Nusno; +{ + + int i, j, k, n, m, N, first_j; + double RADIUS2, PROPER2; + float *X1, *Y1, *X2, *Y2, *usnodist; + float dX, dY, dR, dR2; + int *N1, *N2, *next, last; + int Nave, Nmeas, NMEAS, Nmatch; + unsigned int flags; + Measure *tmpmeasure; + Coords *tcoords; + int already_matched, far_enough; + + X1 = catstats[0].X; + Y1 = catstats[0].Y; + N1 = catstats[0].N; + Nave = catalog[0].Naverage; + + /* no need to do this twice!! */ + ALLOCATE (usnodist, float, Nave); + bzero (usnodist, Nave*sizeof(float)); + + X2 = usnostats[0].X; + Y2 = usnostats[0].Y; + N2 = usnostats[0].N; + + /* set up link listed pointers for new measurements */ + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + /* choose a radius for matches */ + PROPER2 = PROPER*PROPER; + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nave) && (j < Nusno); ) { + if (catalog[0].average[N1[i]].code & ID_MOVING) { + /* this is not a star, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].Nm < 3) { + /* may just be a noise spike, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].code & ID_USNO) { + /* already matched with USNO, skip this one */ + i++; + continue; + } + dX = X1[i] - X2[j]; + if (dX <= -2*PROPER) { + i++; + continue; + } + if (dX >= 2*PROPER) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*PROPER) && (j < Nusno); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < PROPER2) { + n = N1[i]; /* N1 refers to the average[] list */ + N = N2[j]; /* N2 refers to the usno[] list */ + if (usnostats[0].match[N] > -1) + continue; + if ((catalog[0].average[n].code & ID_USNO) && (usnodist[i] < dR)) { + /* existing USNO match is closer than this new one, skip this one */ + continue; + } + usnodist[i] = dR; + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) + m = next[m]; + next[Nmeas+1] = next[m]; /* insert 2 measurements in linked list */ + next[Nmeas] = Nmeas + 1; + next[m] = Nmeas; + if (next[Nmeas+1] == -1) { /* last just was moved */ + last = Nmeas+1; + } + Nmatch ++; + + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - usno[N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - usno[N].D); + catalog[0].measure[Nmeas].M = 1000.0*fabs(usno[N].r); + catalog[0].measure[Nmeas].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = USNO_RED; + catalog[0].measure[Nmeas+1].dR = catalog[0].measure[Nmeas].dR; + catalog[0].measure[Nmeas+1].dD = catalog[0].measure[Nmeas].dD; + catalog[0].measure[Nmeas+1].M = 1000.0*fabs(usno[N].b); + catalog[0].measure[Nmeas+1].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas+1].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas+1].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas+1].averef = n; + catalog[0].measure[Nmeas+1].source = USNO_BLUE; + /* add flag in average to mark as matched with the USNO catalog */ + catalog[0].average[n].code |= (ID_PROPER | ID_USNO); + + /* we add two measurement for each star: 1 in r, 1 in b */ + catalog[0].average[n].Nm +=2; + Nmeas +=2; + if (Nmeas == NMEAS - 2) { /* just to be safe... */ + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + } + } + j = first_j; + i++; + } + + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + catalog[0].Nmeasure = Nmeas; + if (VERBOSE) fprintf (stderr, "Nusno, Nave, Nmeas: %d %d %d, (%d matches)\n", Nusno, Nave, Nmeas, Nmatch); +} + Index: /trunk/Ohana/src/addstar/src/find_subset.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_subset.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/find_subset.c (revision 3339) @@ -0,0 +1,35 @@ +# include "addstar.h" + +Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { + + int i, N, NSUBSET; + Stars **subset; + double RA0, RA1, DEC0, DEC1; + + NSUBSET = 1000; + ALLOCATE (subset, Stars *, NSUBSET); + + RA0 = region[0].RA[0]; + RA1 = region[0].RA[1]; + DEC0 = region[0].DEC[0]; + DEC1 = region[0].DEC[1]; + + if (VERBOSE) fprintf (stderr, "subset for %f - %f, %f - %f\n", RA0, RA1, DEC0, DEC1); + + /* find stars within ra,dec region */ + for (i = N = 0; i < Nstars; i++) { + if (stars[i].R < RA0) continue; + if (stars[i].R >= RA1) continue; + if (stars[i].D < DEC0) continue; + if (stars[i].D >= DEC1) continue; + + subset[N] = &stars[i]; + N++; + if (N == NSUBSET - 1) { + NSUBSET += 1000; + REALLOCATE (subset, Stars, NSUBSET); + } + } + *Nsubset = N; + return (subset); +} Index: /trunk/Ohana/src/addstar/src/get2mass.c =================================================================== --- /trunk/Ohana/src/addstar/src/get2mass.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/get2mass.c (revision 3339) @@ -0,0 +1,35 @@ +# include "mosastro.h" + +StarData *get2mass (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Ngsc, Nregions, Nstars; + StarData *gsc; + StarData *stars; + CatStats *regions; + + Nstars = 0; + ALLOCATE (stars, StarData, 1); + + set_catalog (TWO_MASS_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + if (gsc == NULL) continue; + + Ns = Nstars; + Nstars += Ngsc; + + REALLOCATE (stars, StarData, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + stars[k].R = gsc[j].R; + stars[k].D = gsc[j].D; + stars[k].M = gsc[j].M; + } + free (gsc); + } + + if (VERBOSE) fprintf (stderr, "%d stars from 2MASS\n", Nstars); + *NSTARS = Nstars; + return (stars); +} Index: /trunk/Ohana/src/addstar/src/getgsc.c =================================================================== --- /trunk/Ohana/src/addstar/src/getgsc.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/getgsc.c (revision 3339) @@ -0,0 +1,36 @@ +# include "mosastro.h" + +StarData *getgsc (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Ngsc, Nregions, Nstars; + StarData *gsc; + StarData *stars; + CatStats *regions; + + Nstars = 0; + ALLOCATE (stars, StarData, 1); + + set_catalog (GSC_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + if (gsc == NULL) continue; + + Ns = Nstars; + Nstars += Ngsc; + + REALLOCATE (stars, StarData, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + stars[k].R = gsc[j].R; + stars[k].D = gsc[j].D; + stars[k].Mag = gsc[j].Mag; + } + free (gsc); + } + + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Nstars); + *NSTARS = Nstars; + return (stars); +} + Index: /trunk/Ohana/src/addstar/src/getusno.c =================================================================== --- /trunk/Ohana/src/addstar/src/getusno.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/getusno.c (revision 3339) @@ -0,0 +1,149 @@ +# include "mosastro.h" +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +StarData *getusno (CatStats *catstats, int *Nstars) { + + long int offset; + int i, bin, first, last, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128], c; + FILE *f; + int iRA0, iRA1, iDEC0, iDEC1; + double RA0, RA1, DEC0, DEC1, dec; + int spd, spd_start, spd_end, disk; + int NUSNO, Nusno; + StarData *stars; + + /** USNO-specific things **/ + USNO_RED = GetPhotcodebyName ("USNO_RED"); + if (!USNO_RED) { + fprintf (stderr, "ERROR: USNO_RED photcode not found in photcode table\n"); + exit (0); + } + USNO_BLUE = GetPhotcodebyName ("USNO_BLUE"); + if (!USNO_BLUE) { + fprintf (stderr, "ERROR: USNO_BLUE photcode not found in photcode table\n"); + exit (0); + } + + RA0 = catstats[0].RA[0]; + RA1 = catstats[0].RA[1]; + DEC0 = catstats[0].DEC[0]; + DEC1 = catstats[0].DEC[1]; + + /* identify ra & dec range of interest */ + iRA0 = RA0 * 360000.0; + iRA1 = RA1 * 360000.0; + iDEC0 = (DEC0 + 90.0) * 360000.0; + iDEC1 = (DEC1 + 90.0) * 360000.0; + + /* data is organized in south-pole distance zones */ + spd_start = (int)((DEC0 + 90) / 7.5) * 75.0; + dec = (DEC1 + 90) / 7.5; + if (dec > (int)(dec)) { + spd_end = (int)(1 + (DEC1 + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (DEC1 + 90) / 7.5) * 75.0; + } + + Nusno = 0; + NUSNO = 5000; + ALLOCATE (stars, StarData, NUSNO); + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (1); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = RA0 / 3.75; + if ((RA1 / 3.75) == (int) (RA1 / 3.75)) + last = RA1 / 3.75; + else + last = 1 + RA1 / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s\n", filename); + exit (1); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failure reading data from file %s\n", filename); + exit (1); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if ((buf[0] > iRA0) && (buf[0] < iRA1) && + (buf[1] > iDEC0) && (buf[1] < iDEC1)) { + bzero (&stars[Nusno], sizeof(StarData)); + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + stars[Nusno].Mag = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); + /* b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); */ + Nusno ++; + if (Nusno == NUSNO) { + NUSNO += 5000; + REALLOCATE (stars, StarData, NUSNO); + } + } + } + free (buffer); + } + fclose (f); + } + + *Nstars = Nusno; + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno); + return (stars); +} + + Index: /trunk/Ohana/src/addstar/src/gimages.c =================================================================== --- /trunk/Ohana/src/addstar/src/gimages.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/gimages.c (revision 3339) @@ -1,6 +1,6 @@ # include "addstar.h" -double opening_angle (); -Image *gimages (FILE *f, Image *image, int *Npimage) { +/* given image, find catalog images which overlap it */ +Image *gimages (Image *image, int *Npimage) { int i, j, k, Nimage, addtolist, size; @@ -15,17 +15,18 @@ struct stat filestatus; + f = GetDB (); fseek (f, 0, SEEK_SET); /* read header */ if (!fits_load_header (f, &header)) { - if (VERBOSE) fprintf (stderr, "ERROR: can't read image catalog %s\n", ImageCat); - exit (1); + if (VERBOSE) fprintf (stderr, "error: can't read image catalog %s\n", ImageCat); + Shutdown (); } fits_scan (&header, "ZERO_PT", "%lf", 1, &zeropt); if (fabs (ZeroPt - zeropt) > 1e-4) { - fprintf (stderr, "ERROR: zero point in image (%f:%s) inconsistent with zero point in catalog (%f)\n", + fprintf (stderr, "error: zero point in image (%f:%s) inconsistent with zero point in catalog (%f)\n", ZeroPt, image[0].name, zeropt); - exit (1); + Shutdown (); } @@ -56,15 +57,18 @@ /* check that file size makes sense */ if (stat (ImageCat, &filestatus) == -1) { - fprintf (stderr, "ERROR: failed to get status of image catalog\n"); - exit (1); + fprintf (stderr, "error: failed to get status of image catalog\n"); + Shutdown (); } size = Nimage*sizeof(Image) + header.size; if (size != filestatus.st_size) { int Ndata; + Ndata = (filestatus.st_size - header.size) / sizeof (Image); - fprintf (stderr, "ERROR: image catalog has inconsistent size\n"); - fprintf (stderr, "header: %d, data: %d\n", Nimage, Ndata); + fprintf (stderr, "error: image catalog has inconsistent size\n"); + fprintf (stderr, "header: %d images, data: %d images\n", Nimage, Ndata); if (!FORCE_READ) Shutdown (); - Nimage = Ndata; + + Nimage = MIN (Nimage, Ndata); + fprintf (stderr, "using %d images\n", Nimage); } @@ -76,4 +80,5 @@ ALLOCATE (pimage, Image, NPIMAGE); + /* run through image table, reading images in blocks */ for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) { ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image"); @@ -92,4 +97,5 @@ } /* check if any edges cross */ + /* not every robust - images must be close in size */ for (j = 0; (j < 4) && !addtolist; j++) { for (k = 0; (k < 4) && !addtolist; k++) { @@ -108,6 +114,6 @@ if (ntimage == 0) { - fprintf (stderr, "ERROR: expected %d images, only found %d\n", Nimage, Ntimage); - exit (1); + fprintf (stderr, "error: expected %d images, only found %d\n", Nimage, Ntimage); + Shutdown (); } } Index: /trunk/Ohana/src/addstar/src/greference.c =================================================================== --- /trunk/Ohana/src/addstar/src/greference.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/greference.c (revision 3339) @@ -0,0 +1,35 @@ +# include "addstar.h" + +Stars *grefcat (CatStats *catstats, int *Nrefcat) { + + int Nstars; + Stars *stars; + + if (VERBOSE) fprintf (stderr, "loading astrometric reference data from %s\n", REFCAT); + + Nstars = 0; + if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]); + + /* get stars from USNO for the given region */ + if (!strcmp (REFCAT, "USNO")) { + stars = getusno (catstats, &Nstars); + } + + /* get stars from HST GSC for the given region */ + if (!strcmp (REFCAT, "GSC")) { + stars = getgsc (catstats, &Nstars); + } + + /* get stars from 2MASS for the given region */ + if (!strcmp (REFCAT, "2MASS")) { + stars = get2mass (catstats, &Nstars); + } + + if (Nstars == 0) { + fprintf (stderr, "no ref objs: %s\n", REFCAT); + exit (1); + } + *Nrefcat = Nstars; + return (stars); + +} Index: /trunk/Ohana/src/addstar/src/grefstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/grefstars.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/grefstars.c (revision 3339) @@ -0,0 +1,44 @@ +# include "addstar.h" +# include +# include + +/* read ASCII file with ref star data */ +Stars *grefstars (char *file, int *Nstars) { + + FILE *f; + int N, NSTARS; + Stars *stars; + char line[256]; + + /* open file */ + f = fopen (file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't read data from %s\n", file); + Shutdown (); + } + + NSTARS = 100; + ALLOCATE (stars, Stars, NSTARS); + + /* read in stars line-by-line */ + for (N = 0; scan_line (f, line) != EOF; N++) { + stripwhite (line); + if (line[0] == 0) continue; + if (line[0] == '#') continue; + dparse (&stars[N].R, 1, line); + dparse (&stars[N].D, 2, line); + dparse (&stars[N].M, 3, line); + dparse (&stars[N].dM, 4, line); + while (stars[N].R < 0.0) stars[N].R += 360.0; + while (stars[N].R >= 360.0) stars[N].R -= 360.0; + stars[N].dophot = 0; + stars[N].found = -1; + if (N == NSTARS - 1) { + NSTARS += 100; + REALLOCATE (stars, Stars, NSTARS); + } + } + *Nstars = N; + return (stars); +} + Index: /trunk/Ohana/src/addstar/src/gregion_image.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_image.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/gregion_image.c (revision 3339) @@ -0,0 +1,55 @@ +# include "addstar.h" + +/* given image with coords, find regions which overlap image */ +GSCRegion *gregion_image (Image *image, int *Nregions) { + + GSCRegion *region; + FILE *f; + double x, y; + double dec, ra; + int i, j, done, nregion, NREGION; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC region file %s\n", GSCFILE); + exit (1); + } + + /* find regions at image corners */ + NREGION = 10; + ALLOCATE (region, GSCRegion, NREGION); + nregion = 0; + + /* this is not efficient, nor guaranteed - fix me */ + /* look for new regions on grid across image */ + for (x = 0.0; x <= 1.0; x+=0.1) { + for (y = 0.0; y <= 1.0; y+=0.1) { + XY_to_RD (&ra, &dec, image[0].NX*(1.1*x - 0.05), image[0].NY*(1.1*y - 0.05), &image[0].coords); + aregion (®ion[nregion], f, ra, dec); + done = FALSE; + for (j = 0; (j < nregion) && !done; j++) { + if (!strcmp (region[nregion].filename, region[j].filename)) { + nregion --; + done = TRUE; + } + } + nregion ++; + if (nregion == NREGION) { + NREGION += 10; + REALLOCATE (region, GSCRegion, NREGION); + } + } + } + + if (VERBOSE) { + fprintf (stderr, "found %d region files:\n", nregion); + for (i = 0; i < nregion; i++) { + fprintf (stderr, " %d %s\n", i, region[i].filename); + } + } + *Nregions = nregion; + + fclose (f); + return (region); + +} Index: /trunk/Ohana/src/addstar/src/gregion_match.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_match.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/gregion_match.c (revision 3339) @@ -0,0 +1,51 @@ +# include "addstar.h" +# include + +GSCRegion *gregion_match (int *nregions) { + + int i; + int Ngsc, Nregions; + glob_t pglob; + + GSCRegion *gsc, *regions; + + gsc = LoadRegions (&Ngsc); + + sprintf (globline, "%s/[nsNS]????/*.cpt", CATDIR); + pglob.gl_offs = 0; + glob (globline, 0, NULL, &pglob); + + Nregions = 0; + ALLOCATE (regions, GSCRegions, MAX (1, pglob.gl_pathc)); + + /* match the basename against the GSCRegion file names */ + for (i = 0; i < pglob.gl_pathc; i++) { + found = FALSE; + for (j = 0; !found && (j < Ngsc); j++) { + if (strcasecmp (pglob.gl_pathv, gsc[j].filename)) continue; + regions[Nregions] = gsc[j]; + found = TRUE; + } + } + *nregions = Nregions; + return (regions); + + /* should free things up here */ +} + +/* given a ra,dec bounded region, find all region files which overlap it + + - split region into two if overlapping 0,360 boundary + +*/ + +/* the goal here is to return a list of the existing region files + this means either: + + grabbing all files that are on disk and converting + them to GSCRegion files + + trying all GSCRegion files to see if they exist + (and are populated?) + +*/ Index: /trunk/Ohana/src/addstar/src/gregion_patch.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_patch.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/gregion_patch.c (revision 3339) @@ -0,0 +1,36 @@ +# include "addstar.h" + +GSCRegion *gregion_patch (GSCRegion *patch, int *nregions) { + + int i; + int Ngsc, Nregions; + GSCRegion *gsc, *regions; + + gsc = LoadRegions (&Ngsc); + + /* find lower and upper boundaries */ + Ns = FindRegionDecBandStart (gsc, Ngsc, patch[0].DEC[0]); + Ne = FindRegionDecBandStop (gsc, Ngsc, patch[0].DEC[1]); + + Nregions = 0; + ALLOCATE (regions, GSCRegions, MAX (1, (Ne - Ns))); + + /* find all regions in range within ra range */ + for (i = Ns; i < Ne; i++) { + if (gsc[i].RA[0] >= patch[0].RA[1]) continue; + if (gsc[i].RA[1] < patch[0].RA[0]) continue; + if (gsc[i].DEC[0] >= patch[0].DEC[1]) continue; + if (gsc[i].DEC[1] < patch[0].DEC[0]) continue; + /* add to the list */ + regions[Nregions] = gsc[i]; + Nregions ++; + } + + *nregions = Nregions; + return (regions); +} + +/* given a ra,dec bounded region, find all region files which overlap it + - split region into two if overlapping 0,360 boundary +*/ + Index: /trunk/Ohana/src/addstar/src/gregion_star.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_star.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/gregion_star.c (revision 3339) @@ -0,0 +1,43 @@ +# include "addstar.h" + +GSCRegion *gregion_stars (Stars *stars, int Nstars, int *Nregion) { + + int i; + int Ngsc; + GSCRegion *gsc; + + gsc = LoadRegions (&Ngsc); + + sort_stars (stars, Nstars); + for (i = 0; i < Nstars; i++) { + if (stars[0].found) continue; + region = FindRegionByPoint (star[0].R, star[0].D, gsc, Ngsc); + stars[0].found = TRUE; + /* search forward for all contained stars */ + for (j = i; j < Nstars; j++) { + if (stars[j].R >= region.RA[1]) break; + if (stars[j].D < region.DEC[0]) continue; + if (stars[j].D >= region.DEC[1]) continue; + stars[j].found = TRUE; + } + } + +} + + +/* given a list of stars, find all region files which contain them + + - sort by ra? + - loop over stars + - find region file which contains star + - go backwards until ra < RA[0] + - mark all stars in this region file + - go forwards until ra > RA[1] + - mark all stars in this region file + + - use DEC band information? + - would like to minimize the number of disk reads + +*/ + + Index: /trunk/Ohana/src/addstar/src/gstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/gstars.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/gstars.c (revision 3339) @@ -16,4 +16,7 @@ double tmp; int done, doneread, Nskip, Nextra, itmp; + int hour, min; + double sec, dMs; + char line[80]; if (!fits_read_header (file, &header)) { @@ -52,4 +55,11 @@ } + /* get ST (used for airmass calculation) */ + fits_scan (&header, STKeyword, "%s", 1, line); + /* remove ':' characters */ + for (c = strchr (line, ':'); c != (char *) NULL; c = strchr (line, ':')) { *c = ' '; } + sscanf (line, "%d %d %lf", &hour, &min, &sec); + SiderealTime = hour + min/60.0 + sec/3600.0; + /* CERROR in data file is in arcsec */ if (!fits_scan (&header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; @@ -120,6 +130,6 @@ /* secz is in units milli-airmass */ image[0].Mcal = 0.0; + image[0].Xm = NO_MAG; image[0].code = 0; - image[0].Xm = NO_MAG; bzero (image[0].dummy, sizeof(image[0].dummy)); @@ -131,4 +141,19 @@ } ALLOCATE (stars, Stars, image[0].nstar); + + /* temporary: get sky background from Flips data */ + if (SKYPROBE) { + char *p; + int sky; + + sky = 0; + fits_scan (&header, "HISTORY", "%S", 10, line); + p = strstr (line, "|I="); + if (p != (char *) NULL) { + p += 3; + sky = atoi (p); + } + image[0].Myyyy = MIN (sky - 0x8000, 0xffff); + } /* load in stars by blocks of 1000 */ @@ -200,4 +225,13 @@ while (stars[N].R >= 360.0) stars[N].R -= 360.0; stars[N].found = -1; + + if (SKYPROBE) { + dMs = get_subpix (stars[N].X, stars[N].Y); + stars[N].M -= dMs; + stars[N].Mgal -= dMs; + stars[N].Map -= dMs; + dMs = 1000.0*scat_subpix (stars[N].X, stars[N].Y); + stars[N].dM = hypot (stars[N].dM, dMs); + } } } @@ -210,7 +244,4 @@ return (stars); } - - - # if (0) Index: /trunk/Ohana/src/addstar/src/image-db.c =================================================================== --- /trunk/Ohana/src/addstar/src/image-db.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/image-db.c (revision 3339) @@ -0,0 +1,41 @@ +# include "addstar.h" + +/* these variables are needed by Shutdown */ +static FILE *f = (FILE *) NULL; + +/* clean up open / locked ImageCat before shutting down */ +int Shutdown () { + int dbstate; + + SetProtect (TRUE); + fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); + fprintf (stderr, "ERROR: addstar halted\n"); + exit (1); +} + +FILE *GetDB () { + return (f); +} + +void lock_image_db () { + /* lock the image catalog */ + check_permissions (ImageCat); + f = fsetlockfile (ImageCat, 3600.0, LCK_HARD, &dbstate); + if ((f == (FILE *) NULL) && (dbstate != LCK_EMPTY)) { + fprintf (stderr, "ERROR: can't lock image catalog\n"); + exit (1); + } + fseek (f, 0, SEEK_SET); +} + +void unlock_image_db () { + /* protect wimage from interrupt signals */ + SetProtect (TRUE); + wimage (f, dbstate, &image, Nstars); + SetProtect (FALSE); + + fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); + + mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + chmod (ImageCat, mode); +} Index: /trunk/Ohana/src/addstar/src/load_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/load_catalog.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/load_catalog.c (revision 3339) @@ -0,0 +1,30 @@ +# include "addstar.h" + +void load_catalog (Catalog *catalog, GSCRegion *region) { + + catalog[0].filename = region[0].filename; /* don't free region before catalog! */ + fprintf (stderr, "adding to %s\n", region[0].filename); + + switch (lock_catalog (catalog, LCK_XCLD)) { + case 0: + /* this is an unexpected error because ImageCat was NOT locked */ + fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); + exit (1); + case 1: + gcatalog (catalog); /* load from disk */ + break; + case 2: + mkcatalog (region, catalog); /* fills in new header info */ + break; + } +} + +void save_catalog (Catalog *catalog) { + + /* protect wcatalog from interrupt signals */ + Protect = TRUE; + if (!DUMP_MATCHES && !ONLY_IMAGES) wcatalog (catalog); + if (Trapped) Shutdown (); + Protect = FALSE; + unlock_catalog (catalog); +} Index: /trunk/Ohana/src/addstar/src/load_subpix.c =================================================================== --- /trunk/Ohana/src/addstar/src/load_subpix.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/load_subpix.c (revision 3339) @@ -0,0 +1,53 @@ +# include "addstar.h" + +typedef struct { + double Amp; + double Phase; + double dM; +} SubPixFix; + +static int Nsubpix; +static SubPixFix *Subpix; + +void load_subpix () { + + int i; + FILE *f; + + Nsubpix = 40; + ALLOCATE (Subpix, SubPixFix, Nsubpix); + + f = fopen (SubpixDatafile, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't load subpix datafile %s\n", SubpixDatafile); + Shutdown (); + } + for (i = 0; i < Nsubpix; i++) { + fscanf (f, "%*s %*s %lf %lf %lf %*s\n", + &Subpix[i].Amp, &Subpix[i].Phase, &Subpix[i].dM); + } + fclose (f); + +} + +double get_subpix (double x, double y) { + + int bin; + double dy, dM; + + dy = y - (int)(y); + bin = 5 * (int)(x/100) + (int)(y/100); + dM = Subpix[bin].Amp*sin(2*3.14159*dy + Subpix[bin].Phase); + return (dM); +} + +double scat_subpix (double x, double y) { + + int bin; + double dy, dM; + + dy = y - (int)(y); + bin = 5 * (int)(x/100) + (int)(y/100); + dM = Subpix[bin].dM; + return (dM); +} Index: /trunk/Ohana/src/addstar/src/match_refstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/match_refstars.c (revision 3339) +++ /trunk/Ohana/src/addstar/src/match_refstars.c (revision 3339) @@ -0,0 +1,68 @@ +# include "addstar.h" + +int match_refstars (Stars *stars, int Nstars) { + + int i, j, Nm, Nf, Nmissed, Nsubset; + int *subref; + Stars *subset; + GSCRegion region; + Catalog catalog; + + /* put stars in the appropriate db file */ + for (i = 0; i < Nstars; i++) { + if (stars[i].found != -1) continue; + gregion_star (&stars[i], ®ion); + check_permissions (region.filename); + make_backup (region.filename); + + catalog.filename = region.filename; /* don't free region before catalog! */ + fprintf (stderr, "adding to %s [star %d]\n", region.filename, i); + + switch (lock_catalog (&catalog, LCK_XCLD)) { + case 0: + fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); + exit (1); + case 1: + gcatalog (&catalog); /* load from disk */ + break; + case 2: + mkcatalog (®ion, &catalog); /* fills in new header info */ + break; + } + subset = find_subset (®ion, stars, Nstars, &subref, &Nsubset); + if (VERBOSE) fprintf (stderr, "%d stars in subset for %s\n", Nsubset, region.filename); + find_matches_refstars (®ion, subset, Nsubset, &catalog); + + /* protect wcatalog from interrupt signals */ + SetProtect (TRUE); + wcatalog (&catalog); + SetProtect (FALSE); + unlock_catalog (&catalog); + + Nm = 0; Nf = 0; + for (j = 0; j < Nsubset; j++) { + stars[subref[j]].found = subset[j].found; + if (subset[j].found == -2) Nm ++; + if (subset[j].found >= 0) Nf ++; + } + /* found: + -1 - not yet tried + -2 - found on previous catalog + -3 - not found on previous catalog + +N - found elsewhere + */ + if (VERBOSE) fprintf (stderr, "found %d stars, skipped %d stars\n", Nf, Nm); + free (subset); + free (subref); + } + + for (Nmissed = i = 0; i < Nstars; i++) { + if (stars[i].found == -1) { + fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM); + Nmissed ++; + } + } + if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed); + + return (TRUE); +} Index: /trunk/Ohana/src/addstar/src/parse_time.c =================================================================== --- /trunk/Ohana/src/addstar/src/parse_time.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/parse_time.c (revision 3339) @@ -75,4 +75,5 @@ for (c = strchr (line, 0x2d); c != (char *) NULL; c = strchr (line, 0x2d)) { *c = ' '; } + Nf = 0; switch (mode) { case 1: @@ -106,4 +107,6 @@ } + /* this should probably use localtime */ + /* convert yy.mm.dd hh.mm.ss to Nsec since 1970 (jd = 2440587.5) */ /* note that in this section, tm_mon has range 1-12, unlike for gmtime () */ Index: /trunk/Ohana/src/addstar/src/sort_lists.c =================================================================== --- /trunk/Ohana/src/addstar/src/sort_lists.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/sort_lists.c (revision 3339) @@ -1,2 +1,36 @@ + +void sort_GSCregions (GSCRegion *region, int Nregion) { + + int l,j,ir,i; + GSCRegion tmp; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tmp = region[l]; + } + else { + tmp = region[ir]; + region[ir] = region[0]; + if (--ir == 0) { + region[0] = tmp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && region[j].DEC < region[j+1].DEC) j++; + if (tmp.DEC < region[j].DEC) { + region[i] = region[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + region[i] = tmp; + } +} void sort_lists (float *X, float *Y, int *S, int N) { @@ -79,2 +113,42 @@ } } + +void fsort2 (float *X, float *Y, int N) { + + int l,j,ir,i; + float tX, tY; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + } + else { + tX = X[ir]; + tY = Y[ir]; + X[ir] = X[0]; + Y[ir] = Y[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + } +} Index: /trunk/Ohana/src/addstar/src/wimage.c =================================================================== --- /trunk/Ohana/src/addstar/src/wimage.c (revision 3338) +++ /trunk/Ohana/src/addstar/src/wimage.c (revision 3339) @@ -3,5 +3,5 @@ void wimage (FILE *f, int dbstate, Image *image, int Nstars) { - int Nimages, status; + int Nimages, status, offset; Header header; @@ -37,6 +37,8 @@ } - /* position to end of file for new image data */ - fseek (f, 0, SEEK_END); + /* position to end of data array */ + offset = (Nimages - 1)*sizeof(Image) + header.size; + fseek (f, offset, SEEK_SET); + status = Fwrite (image, sizeof(Image), 1, f, "image"); if (status != 1) {