Index: /trunk/Ohana/src/addstar/Makefile =================================================================== --- /trunk/Ohana/src/addstar/Makefile (revision 3360) +++ /trunk/Ohana/src/addstar/Makefile (revision 3361) @@ -1,6 +1,9 @@ +default: addstar +help: + @echo "make options: addstar (default)" + include ../../Configure + HOME = $(ROOT)/src/addstar -PROGRAM = addstar - BIN = $(HOME)/bin INC = $(HOME)/include @@ -15,43 +18,95 @@ INCS = -I$(INC) -I$(LINC) -I$(XINC) LIBS = -L$(LLIB) -lFITS -lohana -lm -CFLAGS = -o $*.$(ARCH).o $(INCS) -CCFLAGS = $(INCS) $(LIBS) +CFLAGS = $(INCS) +LFLAGS = $(LIBS) ADDSTAR = \ -$(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 +$(SRC)/addstar.$(ARCH).o \ +$(SRC)/airmass.$(ARCH).o \ +$(SRC)/aregion.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/calibrate.$(ARCH).o \ +$(SRC)/check_permissions.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/conversions.$(ARCH).o \ +$(SRC)/edge_check.$(ARCH).o \ +$(SRC)/find_matches.$(ARCH).o \ +$(SRC)/find_matches_refstars.$(ARCH).o \ +$(SRC)/find_subset.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/get2mass.$(ARCH).o \ +$(SRC)/getgsc.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/gimages.$(ARCH).o \ +$(SRC)/greference.$(ARCH).o \ +$(SRC)/grefstars.$(ARCH).o \ +$(SRC)/gregion_image.$(ARCH).o \ +$(SRC)/gregion_match.$(ARCH).o \ +$(SRC)/gregion_patch.$(ARCH).o \ +$(SRC)/gregion_star.$(ARCH).o \ +$(SRC)/gstars.$(ARCH).o \ +$(SRC)/image-db.$(ARCH).o \ +$(SRC)/in_image.$(ARCH).o \ +$(SRC)/load_pt_catalog.$(ARCH).o \ +$(SRC)/load_subpix.$(ARCH).o \ +$(SRC)/make_backup.$(ARCH).o \ +$(SRC)/mkcatalog.$(ARCH).o \ +$(SRC)/opening_angle.$(ARCH).o \ +$(SRC)/parse_time.$(ARCH).o \ +$(SRC)/RegionOps.$(ARCH).o \ +$(SRC)/replace_match.$(ARCH).o \ +$(SRC)/save_pt_catalog.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/sort_lists.$(ARCH).o \ +$(SRC)/update_coords.$(ARCH).o \ +$(SRC)/wcatalog.$(ARCH).o \ +$(SRC)/wimage.$(ARCH).o -OBJ = $(ADDSTAR) +EXTRA = \ +$(SRC)/bracket.$(ARCH).o \ +$(SRC)/find_proper.$(ARCH).o \ +$(SRC)/gregion_match_glob.$(ARCH).o \ +$(SRC)/gregions.$(ARCH).o \ +$(SRC)/match_refstars.$(ARCH).o \ default: $(PROGRAM) -$(ADDSTAR): $(INC)/addstar.h +addstar : $(BIN)/addstar.$(ARCH) -# dependancy rules for binary code ########################## -$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) +$(BIN)/addstar.$(ARCH) : $(ADDSTAR) -$(BIN)/$(PROGRAM).$(ARCH): $(OBJ) +# $(ADDSTAR): $(INC)/addstar.h + +INSTALL = addstar + +# dependancy rules for binary code ######################### +.PRECIOUS: %.$(ARCH).o +.PRECIOUS: $(BIN)/%.$(ARCH) + +%.$(ARCH).o : %.c + $(CC) $(CFLAGS) -c $< -o $@ + +$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(OBJ) -o $(BIN)/$(PROGRAM).$(ARCH) $(CCFLAGS) + $(CC) $^ -o $@ $(LFLAGS) -install: $(DESTBIN)/$(PROGRAM) +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* -$(DESTBIN)/$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) - @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi - rm -f $(DESTBIN)/$(PROGRAM) - cp $(BIN)/$(PROGRAM).$(ARCH) $(DESTBIN)/$(PROGRAM) +$(INSTALL) $(DEVEL): % : $(BIN)/%.$(ARCH) + +%.clean : + rm -f $(BIN)/$*.$(ARCH) + +%.install: + make $(DESTBIN)/$* # utilities ################################################# + +install: + for i in $(INSTALL); do make $$i.install; done + clean: rm -f $(BIN)/*.$(ARCH) @@ -59,9 +114,2 @@ rm -f `find . -name "*~"` rm -f `find . -name "#*"` - -.SUFFIXES: .$(ARCH).o - -.c.$(ARCH).o: - $(CC) $(CFLAGS) -c $< - - Index: /trunk/Ohana/src/addstar/include/addstar.h =================================================================== --- /trunk/Ohana/src/addstar/include/addstar.h (revision 3360) +++ /trunk/Ohana/src/addstar/include/addstar.h (revision 3361) @@ -2,4 +2,18 @@ # include # include +# include +# include + +/* used in find_matches, find_matches_refstars */ +# 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])) + +/* grab named photcode */ +# define NAMED_PHOTCODE(CODE,NAME) \ + CODE = GetPhotcodebyName (NAME); \ + if (!CODE) { \ + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", NAME); \ + exit (0); } typedef struct { @@ -14,5 +28,14 @@ double Mgal, Map; int found; + short int code; + e_time t; } Stars; + +/* structure for data on a catalog region */ +typedef struct { + char filename[256]; + double RA[2]; + int Nrec; +} TM_Region; enum {M_IMAGE, M_REFLIST, M_REFCAT}; @@ -34,4 +57,5 @@ PhotCode *thiscode; +GSCRegion patch; char DateKeyword[64]; @@ -57,4 +81,6 @@ int ONLY_MATCH; +int MODE; + time_t TIMEREF; @@ -70,19 +96,13 @@ 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); double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); -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); +short airmass (short secz_image, 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); @@ -91,5 +111,5 @@ void mkcatalog (GSCRegion *region, Catalog *catalog); Stars *grefstars (char *file, int *Nstars); -Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, 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); @@ -113,4 +133,6 @@ time_t jd_to_sec (double jd); time_t date_to_sec (char *date); +time_t short_date_to_sec (char *date); + /** @@ -122,2 +144,47 @@ e_time pointer from one of these functions) **/ + +Stars *get2mass (GSCRegion *patch, int *NSTARS); +TM_Region *get2mass_acc (GSCRegion *patch, int *nfiles); +Stars *get2mass_data (char *filename, GSCRegion *patch, int *nstars); +GSCRegion *gregion_stars (Stars *stars, int Nstars, int *Nregion); +Stars *grefstars (char *file, int *Nstars); +GSCRegion *LoadRegions (int *nregions); +int FindRegionByPoint (GSCRegion *regions, double ra, double dec, GSCRegion *gsc, int Ngsc); +int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec); +int FindRegionDecBandStop (GSCRegion *region, int Nregion, double dec); +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string); +void find_matches_refstars (GSCRegion *region, Stars **stars, int Nstars, Catalog *catalog); +Stars *gstars (char *file, int *NSTARS, Image *image); +void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap); +Stars *getusno (GSCRegion *catstats, int *Nstars); +Stars *rd_gsc (char *filename, int *nstars); +Stars *getgsc (GSCRegion *patch, int *NSTARS); +Stars *grefcat (char *Refcat, GSCRegion *catstats, int *nstars); +void load_subpix (); +double get_subpix (double x, double y); +double scat_subpix (double x, double y); +void help (); +int args (int argc, char **argv); +void InitCalibration (); +void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N); +void AddToCalibration (Average *average, Measure *measure, Measure *new, int Nstar); +void FindCalibration (Image *image); +int load_pt_catalog (Catalog *catalog, GSCRegion *region); /*** choose new name ***/ +void save_pt_catalog (Catalog *catalog); /*** choose new name ***/ +void wcatalog (Catalog *catalog); +GSCRegion *gregions (Image *image, int *Nregions); +int gcatalog (Catalog *catalog); +GSCRegion *gregion_image (Image *image, int *Nregions); +void update_coords (Average *average, Measure *measure, int *next); +int replace_match (Average *average, Measure *measure, Stars *star); +FILE *GetDB (int *state); +int Shutdown (); +void lock_image_db (); +void unlock_image_db (Image *image); +void wimage (Image *image); +Image *gimages (Image *image, int *Npimage); +double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); +int edge_check (double *x1, double *y1, double *x2, double *y2); +GSCRegion *gregion_match (GSCRegion *regions, int *nregions); +GSCRegion *gregion_patch (GSCRegion *patch, int *nregions); Index: /trunk/Ohana/src/addstar/src/RegionOps.c =================================================================== --- /trunk/Ohana/src/addstar/src/RegionOps.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/RegionOps.c (revision 3361) @@ -177,46 +177,4 @@ */ -/**********/ -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; Index: unk/Ohana/src/addstar/src/addstar-old.c =================================================================== --- /trunk/Ohana/src/addstar/src/addstar-old.c (revision 3360) +++ (revision ) @@ -1,186 +1,0 @@ -# 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 3360) +++ /trunk/Ohana/src/addstar/src/addstar.c (revision 3361) @@ -3,9 +3,8 @@ int main (int argc, char **argv) { - int i, mode, dbstate; - int Nstars, Nimage, Nregions, Nmissed; - Stars *stars; - Image image, *image_db; - GSCRegion *region; + int i, Nstars, Noverlap, Nregions, Nsubset; + Stars *stars, **subset; + Image image, *overlap; + GSCRegion *regions; Catalog catalog; @@ -14,48 +13,50 @@ args (argc, argv); + stars = NULL; + overlap = NULL; + regions = NULL; + if (SKYPROBE) load_subpix (); lock_image_db (); - switch (mode) { + switch (MODE) { case M_IMAGE: stars = gstars (argv[1], &Nstars, &image); - region = gregion_image (&image, &Nregions); + regions = gregion_image (&image, &Nregions); overlap = gimages (&image, &Noverlap); break; case M_REFLIST: stars = grefstars (argv[1], &Nstars); - region = gregion_stars (stars, &Nregion); + regions = gregion_stars (stars, Nstars, &Nregions); break; case M_REFCAT: - regions = gregion_patch (patch, &Nregions); - if (ONLY_MATCH) { - regions = gregion_match (regions, &Nregions); - } + regions = gregion_patch (&patch, &Nregions); break; + } + if (ONLY_MATCH) { + regions = gregion_match (regions, &Nregions); } for (i = 0; i < Nregions; i++) { - load_catalog (&catalog, ®ion[i]); + if (!load_pt_catalog (&catalog, ®ions[i])) continue; - switch (mode) { + switch (MODE) { case M_IMAGE: - find_matches (®ion[i], stars, Nstars, &catalog, &image, overlap, Noverlap); + find_matches (®ions[i], stars, Nstars, &catalog, &image, overlap, Noverlap); break; case M_REFCAT: - stars = grefcat (argv[1], region, &Nstars); + stars = grefcat (argv[1], regions, &Nstars); case M_REFLIST: - subset = find_subset (®ion, stars, Nstars, &subref, &Nsubset); - find_matches_refstars (®ion, subset, Nsubset, &catalog); + subset = find_subset (®ions[i], stars, Nstars, &Nsubset); + find_matches_refstars (®ions[i], subset, Nsubset, &catalog); + if (Nsubset) free (subset); break; } - wcatalog (&catalog); + save_pt_catalog (&catalog); } if (CALIBRATE) { FindCalibration (&image); } - if (DUMP_MATCHES) Shutdown (); - - if (mode == IMAGE) wimage (f, dbstate, &image, Nstars); - unlock_image_db (); + unlock_image_db (&image); exit (0); } @@ -64,5 +65,5 @@ catalog - existing object db table - region - sky area which may or may contain data + regions - 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 3360) +++ /trunk/Ohana/src/addstar/src/airmass.c (revision 3361) @@ -1,7 +1,9 @@ # include "addstar.h" -double airmass (double ra, double dec, double st, double latitude) { +short airmass (short secz_image, double ra, double dec, double st, double latitude) { double hour, cosz, secz; + + if (!SKYPROBE) return (secz_image); /*** make this optional? we may not have ST... ***/ @@ -12,7 +14,6 @@ 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; - + /* modify to match storage units (millimags) */ + secz = 1000.0 / cosz; return (secz); - } Index: /trunk/Ohana/src/addstar/src/args.c =================================================================== --- /trunk/Ohana/src/addstar/src/args.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/args.c (revision 3361) @@ -10,5 +10,4 @@ fprintf (stderr, "\n"); exit (2); - } @@ -26,14 +25,15 @@ /* basic mode: image, list, refcat */ - mode = M_IMAGE; + MODE = M_IMAGE; if ((N = get_argument (argc, argv, "-ref"))) { - mode = M_REFLIST; + MODE = M_REFLIST; remove_argument (N, &argc, argv); } if ((N = get_argument (argc, argv, "-cat"))) { - mode = M_REFCAT; + MODE = M_REFCAT; remove_argument (N, &argc, argv); } + /*** provide additional data ***/ /* restrict to a portion of the sky? (REFCAT only) */ patch.RA[0] = 0; @@ -48,6 +48,24 @@ patch.DEC[1] = atof (argv[N]); } + /* override any header PHOTCODE values */ + thiscode = NULL; + if ((N = get_argument (argc, argv, "-p"))) { + remove_argument (N, &argc, argv); + thiscode = GetPhotcodebyName (argv[N]); + remove_argument (N, &argc, argv); + } + /* provide a time for dataset */ + 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 add to existing objects (REFCAT only) */ + /*** modify behavior ***/ + /* only add to existing objects */ ONLY_MATCH = FALSE; if ((i = get_argument (argc, argv, "-only-match"))) { @@ -55,5 +73,5 @@ remove_argument (i, &argc, argv); } - /* don't add missed pts to Missed table */ + /* don't add missed pts to Missed table (image only) */ SKIP_MISSED = FALSE; if ((N = get_argument (argc, argv, "-missed"))) { @@ -61,5 +79,5 @@ remove_argument (N, &argc, argv); } - /* replace measurement, don't duplicate */ + /* replace measurement, don't duplicate (ref/cat only) */ REPLACE = FALSE; if ((i = get_argument (argc, argv, "-replace"))) { @@ -73,23 +91,12 @@ remove_argument (N, &argc, argv); } - /* apply average zpt offset calibration */ + /* apply average zpt offset calibration (image only) */ CALIBRATE = FALSE; if ((N = get_argument (argc, argv, "-cal"))) { CALIBRATE = TRUE; remove_argument (N, &argc, argv); - CalReference = GetPhotcodeCodebyName (argv[N]); - if (!CalReference) { - fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]); - exit (1); - } - remove_argument (N, &argc, argv); - CalColor = GetPhotcodeCodebyName (argv[N]); - if (!CalColor) { - fprintf (stderr, "ERROR: photcode color %s not found in photcode table\n", argv[N]); - exit (1); - } - remove_argument (N, &argc, argv); } + /*** optional situations ***/ /* treat data specially for skyprobe (calibration, subpix) */ SKYPROBE = FALSE; @@ -98,4 +105,5 @@ remove_argument (N, &argc, argv); } + /* accept bad header astrometry */ ACCEPT_ASTROM = FALSE; if ((N = get_argument (argc, argv, "-accept"))) { @@ -103,9 +111,5 @@ remove_argument (N, &argc, argv); } - DUMP_MATCHES = FALSE; - if ((N = get_argument (argc, argv, "-dump"))) { - DUMP_MATCHES = TRUE; - remove_argument (N, &argc, argv); - } + /* force read of image database with mismatched NSTARS & size */ FORCE_READ = FALSE; if ((N = get_argument (argc, argv, "-force"))) { @@ -113,28 +117,9 @@ remove_argument (N, &argc, argv); } - + /* extra error messages */ VERBOSE = FALSE; if ((N = get_argument (argc, argv, "-v"))) { VERBOSE = TRUE; remove_argument (N, &argc, argv); - } - thiscode = NULL; - if ((N = get_argument (argc, argv, "-p"))) { - remove_argument (N, &argc, argv); - thiscode = GetPhotcodebyName (argv[N]); - remove_argument (N, &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); - } - if (DUMP_MATCHES && !CALIBRATE) { - fprintf (stderr, "-dump only valid -with -cal\n"); - exit (1); } Index: /trunk/Ohana/src/addstar/src/calibrate.c =================================================================== --- /trunk/Ohana/src/addstar/src/calibrate.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/calibrate.c (revision 3361) @@ -1,12 +1,10 @@ # include "addstar.h" -static int Ncal, NCAL; +static int Ncal, NCAL, *Nstar; static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm; -static float *RA, *DEC, *X, *Y; -int *Nstar; void InitCalibration () { - fprintf (stderr, "calibrating the image... %d\n", CalReference); + fprintf (stderr, "calibrating the image... %d\n", thiscode[0].equiv); Ncal = 0; NCAL = 1000; @@ -18,11 +16,7 @@ ALLOCATE (Airm, float, NCAL); ALLOCATE (Nstar, int, NCAL); - ALLOCATE (RA, float, NCAL); - ALLOCATE (DEC, float, NCAL); - ALLOCATE (X, float, NCAL); - ALLOCATE (Y, float, NCAL); } -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 SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N) { Mobs[Ncal] = M; @@ -33,9 +27,4 @@ Airm[Ncal] = A; Nstar[Ncal] = N; - - RA[Ncal] = ra; - DEC[Ncal] = dec; - X[Ncal] = x; - Y[Ncal] = y; Ncal ++; @@ -49,9 +38,44 @@ REALLOCATE (Airm, float, NCAL); REALLOCATE (Nstar, int, NCAL); - REALLOCATE (RA, float, NCAL); - REALLOCATE (DEC, float, NCAL); - REALLOCATE (X, float, NCAL); - REALLOCATE (Y, float, NCAL); } +} + +/* this does not use the linked list to navigate the measure structure + this is OK for old measures, but will break for entries which have + already added a new measure to the list. */ +void AddToCalibration (Average *average, Measure *measure, Measure *new, int Nstar) { + + int i, found0, found1, found2; + float CalM0, CalM1, CalM2, dCalM; + short CalC0, CalC1, CalC2; + + found0 = found1 = found2 = FALSE; + CalM0 = CalM1 = CalM2 = dCalM = NO_MAG; + + CalC0 = thiscode[0].equiv; + CalC1 = thiscode[0].c1; + CalC2 = thiscode[0].c2; + + for (i = 0; i < average[0].Nm; i++) { + if (measure[0].source == CalC0) { + found0 = TRUE; + CalM0 = measure[i].M; + dCalM = measure[i].dM; + } + if (measure[i].source == CalC1) { + found1 = TRUE; + CalM1 = measure[i].M; + } + if (measure[i].source == CalC2) { + found2 = TRUE; + CalM2 = measure[i].M; + } + if (found0 && found1 && found2) { + SaveCalibration (new[0].M, new[0].dM, CalM0, dCalM, CalM1-CalM2, new[0].airmass, Nstar); + return; + } + /* n = next[n]; - if we use this, we need to pass in measure[0] above */ + } + return; } @@ -59,4 +83,5 @@ int i, MaxN, *Nlist, Nkeep; + short int Mint; float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag, *dDmag; float dMo, dMr, Mw, Dmed, W1, W2, NSigma; @@ -91,11 +116,7 @@ /* if this entry has too many (or two few?) matches, skip it */ if (Nlist[Nstar[i]] != 1) continue; - Mabs = 0.001*(Mobs[i] + Klam*(Airm[i] - 1000) + Clam + Xlam*(Mref[i] - Cref[i])) - ZeroPt; - + Mint = Clam + Mobs[i] + Xlam*Cref[i] + Klam*(Airm[i] - 1000); + Mabs = 0.001*Mint - ZeroPt; /* 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 %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 */ @@ -157,57 +178,2 @@ } } - -# 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"); - image[0].Mcal = 10000; - image[0].dMcal = 10000; - return; - } - fsort2 (Dmag, dDmag, Nkeep); - - /* accumulate delta mags (25% - 75% of list) */ - M1 = 0.0; - M2 = 0.0; - N = 0.0; - for (i = 0.25*Nkeep; i <= 0.75*Nkeep; i++) { - - /* a straight mean works better, given the bad photometry of skyprobe */ - M1 += Dmag[i]; - M2 += SQ (Dmag[i]); - N += 1.0; - - } - - if (N > 1) { - M1 = M1 / N; - M2 = sqrt (fabs(M2/N - M1*M1)); - fprintf (stderr, "N: %.0f, mean: %f, stdev: %f, precision: %f\n", N, M1, M2, M2 / sqrt (N)); - image[0].Mcal = 1000 * M1; - image[0].dMcal = 1000 * M2 / sqrt (N); - } else { - fprintf (stderr, "too few stars\n"); - image[0].Mcal = 10000; - 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/check_permissions.c =================================================================== --- /trunk/Ohana/src/addstar/src/check_permissions.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/check_permissions.c (revision 3361) @@ -65,26 +65,3 @@ } } - } - -/* uses cp only */ -void make_backup (char *filename) { - - int status, cmode; - struct stat filestat; - char line [256]; - - status = stat (filename, &filestat); - if (status == 0) { /* file exists, make backup copy */ - sprintf (line, "cp %s %s~", filename, filename); - status = system (line); - if (status) { - fprintf (stderr, "ERROR: unable to create %s~, exiting\n", filename); - exit (1); - } - cmode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; - sprintf (line, "%s~", filename); - chmod (line, cmode); - } - -} Index: /trunk/Ohana/src/addstar/src/conversions.c =================================================================== --- /trunk/Ohana/src/addstar/src/conversions.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/conversions.c (revision 3361) @@ -298,2 +298,67 @@ return (second); } + +/***** short date in format yymmdd *****/ +time_t short_date_to_sec (char *date) { + + time_t second; + double jd; + struct tm now; + + bzero (&now, sizeof(now)); + + scanf (date, "%2d%2d%2d", &now.tm_year, &now.tm_mon, &now.tm_mday); + + if (now.tm_year > 51) now.tm_year += 0; + if (now.tm_year < 50) now.tm_year += 100; + now.tm_mon --; /* tm_mon runs from 0 - 11 */ + + 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); +} + +/**********/ +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); +} Index: /trunk/Ohana/src/addstar/src/edge_check.c =================================================================== --- /trunk/Ohana/src/addstar/src/edge_check.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/edge_check.c (revision 3361) @@ -0,0 +1,36 @@ +# include "addstar.h" + +int edge_check (double *x1, double *y1, double *x2, double *y2) { + + double theta1, theta2; + double Theta1, Theta2; + + theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); + theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + Theta1 = theta1; + Theta2 = theta2; + theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); + theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); + + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + return (TRUE); + +} + Index: /trunk/Ohana/src/addstar/src/find_matches.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/find_matches.c (revision 3361) @@ -1,14 +1,9 @@ # 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 (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap) { - 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, dCalM, *Xs, *Ys; + int i, j, k, n, m, N, J; + double X, Y, RADIUS, RADIUS2, secz; + float *X1, *Y1, *X2, *Y2; float dX, dY, dR; int *N1, *N2, *next, *next_miss, last, last_miss; @@ -21,7 +16,7 @@ /* photcode data */ - CalM = CalC = dCalM = 0; Nsecfilt = GetPhotcodeNsecfilt (); Nsec = GetPhotcodeNsec (thiscode[0].code); + /* this function requires incoming stars to have the same photcode */ /** allocate local arrays (stars) **/ @@ -36,6 +31,4 @@ ALLOCATE (Y2, float, NAVE); ALLOCATE (N2, int, NAVE); - ALLOCATE (Xs, float, NAVE); - ALLOCATE (Ys, float, NAVE); ALLOCATE (catalog[0].found, int, NAVE); REALLOCATE (catalog[0].average, Average, NAVE); @@ -63,5 +56,5 @@ for (i = 0; i < Nstars; i++) { - fRD_to_XY (&X1[i], &Y1[i], stars[0][i].R, stars[0][i].D, &tcoords); + fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); N1[i] = i; } @@ -69,5 +62,4 @@ for (i = 0; i < Nave; i++) { - fRD_to_XY (&Xs[i], &Ys[i], catalog[0].average[i].R, catalog[0].average[i].D, &image[0].coords); fRD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); N2[i] = i; @@ -93,5 +85,5 @@ RADIUS = NSIGMA * 0.02 * image[0].cerror; /* 0.02 corrects cerror to arcsec from storage units */ } else { - RADIUS = DEFAULT_RADIUS; + RADIUS = DEFAULT_RADIUS; /* provided by config */ } RADIUS2 = RADIUS*RADIUS; @@ -108,5 +100,4 @@ dX = X1[i] - X2[j]; - if (dX <= -2*RADIUS) { i++; @@ -119,147 +110,89 @@ /* 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]; + for (J = j; (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]; - - /** 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; + if (dR < RADIUS2) continue; + + Nmatch ++; + n = N2[J]; + N = N1[i]; + m = catalog[0].average[n].offset; + + /** insert star in measurement list */ + /* find last measurement of this star */ + 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; - /* 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; - found0 = FALSE; - found1 = FALSE; - 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; - 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[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]); - } - } - } - - /** add measurements for this star **/ - 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[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[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[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 */ + /* calculate accurate per-star airmass */ + secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude); + + + /** 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 (1000*stars[N].M + MTIME, NO_MAG); + catalog[0].measure[Nmeas].dM = MIN (1000*stars[N].dM, NO_ERR); /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].Mcal = image[0].Mcal; + 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].averef = n; + catalog[0].measure[Nmeas].source = stars[N].code; /* photcode */ + catalog[0].measure[Nmeas].dophot = stars[N].dophot; + catalog[0].measure[Nmeas].flags = 0; + catalog[0].measure[Nmeas].dt = MTIME; + catalog[0].measure[Nmeas].airmass = secz; + + catalog[0].measure[Nmeas].Mgal = MIN (1000*stars[N].Mgal + MTIME, NO_MAG); + catalog[0].measure[Nmeas].FWx = MIN (100*stars[N].fx, NO_MAG); + catalog[0].measure[Nmeas].fwy = MIN (100*(stars[N].fy / stars[N].fx), NO_ERR); + catalog[0].measure[Nmeas].theta = MIN ((255/360)*stars[N].df, NO_ERR); - /* check for entries in the secfilt lists */ - Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K; - Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; - if (*Mval == NO_MAG) *Mval = Mrel; - - /*** handle multiple stars */ - /* this image star matches more than one catalog star */ - 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[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[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 */ - 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; + /* it is not valid to pass PRI/SEC/REF photcodes to this routine */ + /* check for entries in the secfilt lists */ + Mrel = catalog[0].measure[Nmeas].M + secz * thiscode[0].K; + Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; + if (*Mval == NO_MAG) *Mval = Mrel; + + if (CALIBRATE) { + /** this is a little tricky: the new measures are not in the measure sequence + this call can run into the next star if we have had two measures found at this location. **/ + AddToCalibration (&catalog[0].average[n], &catalog[0].measure[m], &catalog[0].measure[Nmeas], N); + } + + /*** flag 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_coords (&catalog[0].average[n], &catalog[0].measure[0], next); + } i++; } @@ -270,30 +203,30 @@ if (catalog[0].found[n] < 0) { /* should the catalog star be on this image? project into image coords */ - if (in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) { - /* find last missing meas */ - if (catalog[0].average[n].Nn < 1) { /* no previous missing obs */ - catalog[0].average[n].missing = Nmiss; - next_miss[last_miss] = Nmiss; - next_miss[Nmiss] = -1; + if (!in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) continue; + /* find last missing meas */ + if (catalog[0].average[n].Nn < 1) { /* no previous missing obs */ + catalog[0].average[n].missing = Nmiss; + next_miss[last_miss] = Nmiss; + next_miss[Nmiss] = -1; + last_miss = Nmiss; + } else { + m = catalog[0].average[n].missing; + for (k = 0; k < catalog[0].average[n].Nn - 1; k++) + m = next_miss[m]; + next_miss[Nmiss] = next_miss[m]; + next_miss[m] = Nmiss; + if (next_miss[Nmiss] == -1) { /* last just was moved */ last_miss = Nmiss; - } else { - m = catalog[0].average[n].missing; - for (k = 0; k < catalog[0].average[n].Nn - 1; k++) - m = next_miss[m]; - next_miss[Nmiss] = next_miss[m]; - next_miss[m] = Nmiss; - if (next_miss[Nmiss] == -1) { /* last just was moved */ - last_miss = Nmiss; - } } - RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords); - catalog[0].missing[Nmiss].t = image[0].tzero + 1e-4*Y*image[0].trate; /* trate is in 0.1 msec / row */ - catalog[0].average[n].Nn ++; - Nmiss ++; - if (Nmiss == NMISS) { - NMISS = Nmiss + 1000; - REALLOCATE (next_miss, int, NMISS); - REALLOCATE (catalog[0].missing, Missing, NMISS); - } + } + /* calculate time of exposure for this coordinate in the image */ + RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords); + catalog[0].missing[Nmiss].t = image[0].tzero + 1e-4*Y*image[0].trate; /* trate is in 0.1 msec / row */ + catalog[0].average[n].Nn ++; + Nmiss ++; + if (Nmiss == NMISS) { + NMISS = Nmiss + 1000; + REALLOCATE (next_miss, int, NMISS); + REALLOCATE (catalog[0].missing, Missing, NMISS); } } @@ -301,90 +234,88 @@ /* incorporate unmatched image stars, if this star is in field of this catalog */ - for (i = 0; i < Nstars; i++) { + for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { N = N1[i]; - 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++) { - catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; - catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; - } - - 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; - 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 = 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[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[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[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; - last = Nmeas; - Nmeas ++; - if (Nmeas == NMEAS) { - NMEAS = Nmeas + 1000; - REALLOCATE (next, int, NMEAS); - 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 < Noverlap) && !SKIP_MISSED; j++) { - if (in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[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, &overlap[j].coords); - catalog[0].missing[Nmiss].t = overlap[j].tzero + 1e-4*Y*overlap[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); - } - - } - } - Nave ++; - if (Nave == NAVE) { - NAVE = Nave + 1000; - REALLOCATE (catalog[0].average, Average, NAVE); - REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); - } + if (stars[N].found >= 0) continue; + if (!IN_CATALOG (stars[N].R, stars[N].D)) continue; + + secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude); + + catalog[0].average[Nave].R = stars[N].R; + catalog[0].average[Nave].D = stars[N].D; + catalog[0].average[Nave].M = NO_MAG; + catalog[0].average[Nave].dM = NO_MAG; + 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].Xg = NO_MAG; + catalog[0].average[Nave].offset = Nmeas; + catalog[0].average[Nave].missing = -1; + catalog[0].average[Nave].code = 0; + + for (j = 0; j < Nsecfilt; j++) { + catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; + catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; + catalog[0].secfilt[Nave*Nsecfilt+j].dM = NO_MAG; + } + + catalog[0].measure[Nmeas].dR = 0.0; + catalog[0].measure[Nmeas].dD = 0.0; + catalog[0].measure[Nmeas].M = MIN (1000*stars[N].M + MTIME, NO_MAG); + catalog[0].measure[Nmeas].dM = MIN (1000*stars[N].dM, NO_ERR); + catalog[0].measure[Nmeas].Mcal = image[0].Mcal; + 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].averef = Nave; + catalog[0].measure[Nmeas].source = stars[N].code; /* photcode */ + catalog[0].measure[Nmeas].dophot = stars[N].dophot; + catalog[0].measure[Nmeas].flags = 0; + catalog[0].measure[Nmeas].dt = MTIME; + catalog[0].measure[Nmeas].airmass = secz; + + catalog[0].measure[Nmeas].Mgal = MIN (1000*stars[N].Mgal + MTIME, NO_MAG); + catalog[0].measure[Nmeas].FWx = MIN (100*stars[N].fx, NO_MAG); + catalog[0].measure[Nmeas].fwy = MIN (100*(stars[N].fy / stars[N].fx), NO_ERR); + catalog[0].measure[Nmeas].theta = MIN ((255/360)*stars[N].df, NO_ERR); + + Mrel = catalog[0].measure[Nmeas].M + secz * thiscode[0].K; + Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M; + if (*Mval == NO_MAG) *Mval = Mrel; + + /** now add references from all previous non-detection observations of this spot on the sky */ + for (j = 0; (j < Noverlap) && !SKIP_MISSED; j++) { + if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue; + 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; + /* get time of exposure of this portion of the image */ + RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j].coords); + catalog[0].missing[Nmiss].t = overlap[j].tzero + 1e-4*Y*overlap[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); + } + } + + 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); } } @@ -393,5 +324,5 @@ REALLOCATE (catalog[0].average, Average, Nave); REALLOCATE (catalog[0].measure, Measure, Nmeas); - + /* fix order of Measure (memory intensive, but fast) */ N = 0; @@ -427,6 +358,8 @@ /* note stars which have been found in this catalog */ for (i = 0; i < Nstars; i++) { - if (stars[0][i].found > -1) { - stars[0][i].found = -2; + if (stars[i].found > -1) { + stars[i].found = -2; + } else { + stars[i].found = -3; } } @@ -438,95 +371,13 @@ } -int in_image (double r, double d, Image *image) { - - double X, Y; - - RD_to_XY (&X, &Y, r, d, &image[0].coords); - if (((X) >= 0) && ((X) < image[0].NX) && ((Y) >= 0) && ((Y) < image[0].NY)) - return (TRUE); - else - return (FALSE); - -} - /* notes: - - for finding if a catalog star is in an image or an image star is in the catalog: - - catalogs have boundaries defined by RA and DEC, but they may curve in projection - images have boundaries which are lines in pixels coords, but curve in RA and DEC - - catalog[0].found[Ncat] but stars[Nstar].found - - */ - - -# ifdef OLDJUNKHERE - - - /* fix measurement list order */ - for (i = 0; i < Nmeas - 1; i++) { - if (next[i] != i + 1) { - n0 = next[i]; - n1 = next[i+1]; - n2 = next[n0]; - t = i; - while ((t != -1) && (next[t] != i+1)) - t = next[t]; - if ((t == -1) || (n2 == i+1)) { - t = n2; - } - next[i] = i+1; - next[i+1] = n2; - next[n0] = n1; - next[t] = n0; - tmpmeasure = catalog[0].measure[n0]; - catalog[0].measure[n0] = catalog[0].measure[i+1]; - catalog[0].measure[i+1] = tmpmeasure; - } - } - - /* fix offset references */ - n = (0x00ffffff & catalog[0].measure[0].average); - catalog[0].average[n].offset = 0; - for (i = 1; i < Nmeas; i++) { - if (n != (0x00ffffff & catalog[0].measure[i].average)) { - n = (0x00ffffff & catalog[0].measure[i].average); - catalog[0].average[n].offset = i; - } - } - - /* fix missing list order */ - for (i = 0; i < Nmiss - 1; i++) { - if (next_miss[i] != i + 1) { - n0 = next_miss[i]; - n1 = next_miss[i+1]; - n2 = next_miss[n0]; - t = i; - while ((t != -1) && (next_miss[t] != i+1)) - t = next_miss[t]; - if ((t == -1) || (n2 == i+1)) { - t = n2; - } - next_miss[i] = i+1; - next_miss[i+1] = n2; - next_miss[n0] = n1; - next_miss[t] = n0; - tmpmissing = catalog[0].missing[n0]; - catalog[0].missing[n0] = catalog[0].missing[i+1]; - catalog[0].missing[i+1] = tmpmissing; - } - } - - /* fix offset references */ - n = (0x00ffffff & catalog[0].measure[0].average); - catalog[0].average[n].offset = 0; - for (i = 1; i < Nmeas; i++) { - if (n != (0x00ffffff & catalog[0].measure[i].average)) { - n = (0x00ffffff & catalog[0].measure[i].average); - catalog[0].average[n].offset = i; - } - } - -# endif + + for finding if a catalog star is in an image or an image star is in the catalog: + + catalogs have boundaries defined by RA and DEC, but they may curve in projection + images have boundaries which are lines in pixels coords, but curve in RA and DEC + + catalog[0].found[Ncat] but stars[Nstar].found + +*/ Index: /trunk/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/find_matches_refstars.c (revision 3361) @@ -1,11 +1,7 @@ # 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; + int i, j, k, n, m, N, J; double RADIUS, RADIUS2; float *X1, *Y1, *X2, *Y2; @@ -16,9 +12,8 @@ Missing *tmpmissing; Coords tcoords; - int Nsecfilt, Nsec; + int Nsecfilt; /* photcode data -- should not have to modify secfilt / average */ Nsecfilt = GetPhotcodeNsecfilt (); - Nsec = GetPhotcodeNsec (thiscode[0].code); /** allocate local arrays (stars) **/ @@ -88,7 +83,7 @@ /* choose a radius for matches */ if (DEFAULT_RADIUS == 0) { - RADIUS = 2.0; + RADIUS = 2.0; /* hardwired default for refstars */ } else { - RADIUS = DEFAULT_RADIUS; + RADIUS = DEFAULT_RADIUS; /* provided by config */ } RADIUS2 = RADIUS*RADIUS; @@ -98,5 +93,4 @@ dX = X1[i] - X2[j]; - if (dX <= -2*RADIUS) { i++; @@ -109,175 +103,37 @@ /* 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]; + for (J = j; (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; + if (dR > RADIUS2) continue; + + Nmatch ++; + n = N2[J]; + N = N1[i]; + m = catalog[0].average[n].offset; + + /** in replace mode, search for entry and replace values M, dM, R, D */ + if (REPLACE && replace_match (&catalog[0].average[n], &catalog[0].measure[m], stars[N])) continue; + + /** insert star in measurement list */ + /* find last measurement of this star */ + 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; + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N][0].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N][0].D); + catalog[0].measure[Nmeas].M = MIN (1000*stars[N][0].M, NO_MAG); + catalog[0].measure[Nmeas].dM = MIN (1000*stars[N][0].dM, NO_ERR); 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].t = (TIMEREF == 0) ? stars[N][0].t : TIMEREF; /** careful : time_t vs e_time **/ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = stars[N][0].code; + catalog[0].measure[Nmeas].dophot = 0; catalog[0].measure[Nmeas].flags = 0; catalog[0].measure[Nmeas].dt = 0xffff; @@ -288,9 +144,28 @@ catalog[0].measure[Nmeas].fwy = 0xff; catalog[0].measure[Nmeas].theta = 0xff; - - stars[N].found = Nmeas; - next[last] = Nmeas; - next[Nmeas] = -1; - last = Nmeas; + + /** don't update average / secfilt values for REF photcodes **/ + + /*** handle multiple stars */ + /* this image star matches more than one catalog star */ + if (stars[N][0].found > -1) { + catalog[0].measure[stars[N][0].found].flags |= BLEND_IMAGE; + catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; + } + if (stars[N][0].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ + catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; + } + if (stars[N][0].found == -1) { /* this image star matches only this catalog star */ + stars[N][0].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) { @@ -299,10 +174,74 @@ 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); - } + + update_coords (&catalog[0].average[n], &catalog[0].measure[0], next); + } + i++; + } + + /* we don't add missed entries for refcat + (already in database, not refcat) */ + + /* incorporate unmatched refcat stars */ + /* skip if we want to require 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][0].found >= 0) continue; + + catalog[0].average[Nave].R = stars[N][0].R; + catalog[0].average[Nave].D = stars[N][0].D; + catalog[0].average[Nave].M = NO_MAG; + catalog[0].average[Nave].dM = NO_MAG; + 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].Xg = NO_MAG; + catalog[0].average[Nave].offset = Nmeas; + catalog[0].average[Nave].missing = -1; + catalog[0].average[Nave].code = 0; + + for (j = 0; j < Nsecfilt; j++) { + catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; + catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; + catalog[0].secfilt[Nave*Nsecfilt+j].dM = NO_MAG; + } + + catalog[0].measure[Nmeas].dR = 0.0; + catalog[0].measure[Nmeas].dD = 0.0; + catalog[0].measure[Nmeas].M = MIN (1000*stars[N][0].M, NO_MAG); + catalog[0].measure[Nmeas].dM = MIN (1000*stars[N][0].dM, NO_ERR); + catalog[0].measure[Nmeas].Mcal = 0; + catalog[0].measure[Nmeas].t = (stars[N][0].t == 0) ? TIMEREF : stars[N][0].t; /** careful : time_t vs e_time **/ + catalog[0].measure[Nmeas].averef = Nave; + catalog[0].measure[Nmeas].source = stars[N][0].code; + catalog[0].measure[Nmeas].dophot = 0; + 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][0].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); } } @@ -345,8 +284,8 @@ /* 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[i][0].found > -1) { + stars[i][0].found = -2; } else { - stars[i].found = -3; + stars[i][0].found = -3; } } @@ -357,31 +296,2 @@ 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_subset.c =================================================================== --- /trunk/Ohana/src/addstar/src/find_subset.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/find_subset.c (revision 3361) @@ -1,5 +1,6 @@ # include "addstar.h" -Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { +/* find stars within this region */ +Stars **find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { int i, N, NSUBSET; @@ -28,5 +29,5 @@ if (N == NSUBSET - 1) { NSUBSET += 1000; - REALLOCATE (subset, Stars, NSUBSET); + REALLOCATE (subset, Stars *, NSUBSET); } } Index: /trunk/Ohana/src/addstar/src/gcatalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/gcatalog.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/gcatalog.c (revision 3361) @@ -25,10 +25,12 @@ for (in = out = i = 0; i < catalog[0].Naverage; i++) { for (j = 0; j < catalog[0].Nsecfilt; j++, in++, out++) { - outsec[out].M = insec[in].M; + outsec[out].M = insec[in].M; outsec[out].Xm = insec[in].Xm; + outsec[out].dM = insec[in].dM; } for (j = 0; j < Nextra; j++, out++) { - outsec[out].M = NO_MAG; + outsec[out].M = NO_MAG; outsec[out].Xm = NO_MAG; + outsec[out].dM = NO_MAG; } } Index: /trunk/Ohana/src/addstar/src/get2mass.c =================================================================== --- /trunk/Ohana/src/addstar/src/get2mass.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/get2mass.c (revision 3361) @@ -1,31 +1,29 @@ -# include "mosastro.h" +# include "addstar.h" -StarData *get2mass (CatStats *catstats, int *NSTARS) { +Stars *get2mass (GSCRegion *patch, int *NSTARS) { - int i, j, k, Ns, Ngsc, Nregions, Nstars; - StarData *gsc; - StarData *stars; - CatStats *regions; + int i, j, k, Ns, Nstars, Nrefcat, Nregions; + Stars *stars; + Stars *refcat; + TM_Region *regions; + short int TM_J, TM_H, TM_K; + NAMED_PHOTCODE (TM_J, "2MASS_J"); + NAMED_PHOTCODE (TM_H, "2MASS_H"); + NAMED_PHOTCODE (TM_K, "2MASS_K"); + + regions = get2mass_acc (patch, &Nregions); + Nstars = 0; - ALLOCATE (stars, StarData, 1); - - set_catalog (TWO_MASS_DIR); - regions = gregions (catstats, &Nregions); + ALLOCATE (stars, Stars, 1); for (i = 0; i < Nregions; i++) { - gsc = gcatalog (regions[i].name, &Ngsc); - if (gsc == NULL) continue; + refcat = get2mass_data (regions[i].name, patch, &Nrefcat); - Ns = Nstars; - Nstars += Ngsc; + REALLOCATE (stars, Stars, MAX (1, Nstars + Nrefcat)); + memcpy (&stars[Nstars], refcat, Nrefcat*sizeof(Stars)); + Nstars += Nrefcat; - REALLOCATE (stars, StarData, 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); + free (refcat); } @@ -34,2 +32,132 @@ return (stars); } + +/* watch for patch which cross 0,360 boundary */ +TM_Region *get2mass_acc (GSCRegion *patch, int *nfiles) { + + int Nregions + TM_Region *regions; + + sprintf (accel, "%s/idr2.acc", TWO_MASS_DIR); + f = fopen (accel, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't read data from %s\n", file); + Shutdown (); + } + + Nregions = 0; + NREGIONS = 20; + ALLOCATE (regions, TM_Region, NREGIONS); + + /* 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; + scanf (line, "%s %lf %lf %d", filename, &Rs, &Re, &Nrec); + if (Rs > patch[0].RA[1]) continue; + if (Re < patch[0].RA[0]) continue; + strcpy (regions[Nregions].filename, "%s/%s", TWO_MASS_DIR, filename); + regions[Nregions].RA[0] = Rs; + regions[Nregions].RA[1] = Re; + regions[Nregions].Nrec = Nrec; + Nregions ++; + CHECK_REALLOCATE (regions, TM_Region, NREGIONS, Nregions, 20); + } + + *nregions = Nregions; + return (regions); +} + +# define NBYTE 302 +# define NLINE 100 + +Stars *get2mass_data (char *filename, GSCRegion *patch, int *nstars) { + + int Nstars; + Stars *stars; + gzFile gf; + char buffer[NLINE*NBYTE]; + + gf = gzopen (filename, "rb"); + if (gf == NULL) { + fprintf (stderr, "ERROR: can't read data from %s\n", file); + Shutdown (); + } + + Nstars = 0; + NSTARS = 10000; + ALLOCATE (stars, Stars, NSTARS); + + /* z_off_t gzseek (gzFile file, z_off_t offset, int whence); */ + + while ((Nbyte = gzread (gf, buffer, NLINE*NBYTE)) != 0) { + + if (Nbyte == -1) { + fprintf (stderr, "error reading from gzipped file %s\n", filename); + Shutdown (); + } + + Nfrac = Nbyte % NBYTE; + if (Nfrac) { + fprintf (stderr, "error reading complete line from gzipped file %s\n", filename); + Shutdown (); + } + + Nline = Nbyte / NBYTE; + + for (i = 0; i < Nline; i++) { + Nval = scanf (&buffer[NBYTE*i + 0], "%lf %lf", &RA, &DEC); + if (RA < patch[0].RA[0]) continue; + if (RA > patch[0].RA[1]) continue; + if (DEC < patch[0].DEC[0]) continue; + if (DEC < patch[0].DEC[1]) continue; + + Nval += scanf (&buffer[NBYTE*i + 53], "%lf %lf", &J, &dJ); + Nval += scanf (&buffer[NBYTE*i + 72], "%lf %lf", &H, &dH); + Nval += scanf (&buffer[NBYTE*i + 91], "%lf %lf", &K, &dK); + Nval += scanf (&buffer[NBYTE*i + 164], "%s", datestr); + if (Nval != 9) { + fprintf (stderr, "error in 2mass data line\n"); + exit (1); + } + + time = short_date_to_set (datestr); + + stars[Nstars].R = RA; + stars[Nstars].D = DEC; + stars[Nstars].M = J; + stars[Nstars].dM = dJ; + stars[Nstars].t = time; + stars[Nstars].code = TM_J; + Nstars ++; + + stars[Nstars].R = RA; + stars[Nstars].D = DEC; + stars[Nstars].M = H; + stars[Nstars].dM = dH; + stars[Nstars].t = time; + stars[Nstars].code = TM_H; + Nstars ++; + + stars[Nstars].R = RA; + stars[Nstars].D = DEC; + stars[Nstars].M = K; + stars[Nstars].dM = dK; + stars[Nstars].t = time; + stars[Nstars].code = TM_K; + Nstars ++; + + CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars - 2, 300); + } + } + + *nstars = Nstars; + return (stars); +} + + /* this just scans along in the file. file is sorted by dec, so we + should be skipping large chunks - but we would need to have + the size from the accel file (won't fit in GSCRegion) and need + to use gzseek, if it exists. + */ Index: /trunk/Ohana/src/addstar/src/getgsc.c =================================================================== --- /trunk/Ohana/src/addstar/src/getgsc.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/getgsc.c (revision 3361) @@ -1,30 +1,29 @@ -# include "mosastro.h" +# include "addstar.h" +# define BYTES_STAR 23 +# define BLOCK 1000 +# define DNSTARS 1000 -StarData *getgsc (CatStats *catstats, int *NSTARS) { +Stars *getgsc (GSCRegion *patch, int *NSTARS) { int i, j, k, Ns, Ngsc, Nregions, Nstars; StarData *gsc; StarData *stars; - CatStats *regions; + GSCRegions *regions; + short int GSC_M; + NAMED_PHOTCODE (GSC_M, "GSC_M"); + + regions = gregion_patch (patch, &Nregions); + Nstars = 0; - ALLOCATE (stars, StarData, 1); - - set_catalog (GSC_DIR); - regions = gregions (catstats, &Nregions); + ALLOCATE (stars, Stars, 1); for (i = 0; i < Nregions; i++) { - gsc = gcatalog (regions[i].name, &Ngsc); - if (gsc == NULL) continue; + gsc = rd_gsc (regions[i].name, &Ngsc); - Ns = Nstars; + REALLOCATE (stars, Stars, MAX (1, Nstars + Ngsc)); + memcpy (&stars[Nstars], gsc, Ngsc*sizeof(Stars)); Nstars += Ngsc; - REALLOCATE (stars, StarData, 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); } @@ -35,2 +34,41 @@ } +Stars *rd_gsc (char *filename, int *nstars) { + + int i, NSTARS, Nstars, Nbytes, nbytes; + char *buffer; + FILE *f; + Stars *stars; + + Nstars = 0; + NSTARS = 1000; + ALLOCATE (stars, Stars, NSTARS); + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find catalog file %s\n", filename); + *nstars = Nstars; + return (stars); + } + + Nbytes = BLOCK*BYTES_STAR; + ALLOCATE (buffer, char, Nbytes); + while ((nbytes = fread (buffer, 1, Nbytes, f)) > 0) { + Nline = nbytes / BYTES_STAR; + for (i = 0; i < Nline; i++) { + Nbyte = i*BYTES_STAR; + dparse (&stars[Nstars].R, 1, &buffer[Nbyte]); + dparse (&stars[Nstars].D, 2, &buffer[Nbyte]); + dparse (&stars[Nstars].M, 3, &buffer[Nbyte]); + stars[Nstars].dM = NO_ERR; + stars[Nstars].t = 0; + stars[Nstars].code = GSC_M; + Nstars ++; + CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 1000); + } + } + free (buffer); + REALLOCATE (stars, Stars, Nstars); + *nstars = Nstars; + return (stars); +} Index: /trunk/Ohana/src/addstar/src/getusno.c =================================================================== --- /trunk/Ohana/src/addstar/src/getusno.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/getusno.c (revision 3361) @@ -1,3 +1,3 @@ -# include "mosastro.h" +# include "addstar.h" # define NZONE 24 @@ -9,5 +9,5 @@ 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; -Stars *getusno (CatStats *catstats, int *Nstars) { +Stars *getusno (GSCRegion *catstats, int *Nstars) { long int offset; @@ -18,39 +18,26 @@ FILE *f; int iRA0, iRA1, iDEC0, iDEC1; - double RA0, RA1, DEC0, DEC1, dec; + double dec; int spd, spd_start, spd_end, disk; int NUSNO, Nusno; Stars *stars; + short int USNO_RED, USNO_BLUE; - /** 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]; + NAMED_PHOTCODE (USNO_RED, "USNO_RED"); + NAMED_PHOTCODE (USNO_BLUE, "USNO_BLUE"); /* 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; + iRA0 = catstats[0].RA[0] * 360000.0; + iRA1 = catstats[0].RA[1] * 360000.0; + iDEC0 = (catstats[0].DEC[0] + 90.0) * 360000.0; + iDEC1 = (catstats[0].DEC[1] + 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; + spd_start = (int)((catstats[0].DEC[0] + 90) / 7.5) * 75.0; + dec = (catstats[0].DEC[1] + 90) / 7.5; if (dec > (int)(dec)) { - spd_end = (int)(1 + (DEC1 + 90) / 7.5) * 75.0; + spd_end = (int)(1 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; } else { - spd_end = (int)(0 + (DEC1 + 90) / 7.5) * 75.0; + spd_end = (int)(0 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; } @@ -89,9 +76,9 @@ fclose (f); - first = RA0 / 3.75; - if ((RA1 / 3.75) == (int) (RA1 / 3.75)) - last = RA1 / 3.75; + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; else - last = 1 + RA1 / 3.75; + last = 1 + catstats[0].RA[1] / 3.75; if ((first > Nbins) || (last > Nbins)) { @@ -123,21 +110,28 @@ /* 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(Stars)); - 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)))); - stars[Nusno].code = USNO_RED; - /* r data */ - Nusno ++; - stars[Nusno].R = buf[0]/360000.0; - stars[Nusno].D = buf[1]/360000.0 - 90.0; - stars[Nusno].Mag = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000)); - stars[Nusno].code = USNO_BLUE; - /* b data */ - Nusno ++; - CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno-1, 5000); - } + if (buf[0] < iRA0) continue; + if (buf[0] > iRA1) continue; + if (buf[1] < iDEC0) continue; + if (buf[1] > iDEC1) continue; + + bzero (&stars[Nusno], sizeof(Stars)); + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + stars[Nusno].M = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); + stars[Nusno].dM = NO_ERR; + stars[Nusno].t = 0; + stars[Nusno].code = USNO_RED; + /* r data */ + Nusno ++; + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + stars[Nusno].M = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000)); + stars[Nusno].dM = NO_ERR; + stars[Nusno].t = 0; + stars[Nusno].code = USNO_BLUE; + /* b data */ + Nusno ++; + + CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno-1, 5000); } free (buffer); Index: /trunk/Ohana/src/addstar/src/gimages.c =================================================================== --- /trunk/Ohana/src/addstar/src/gimages.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/gimages.c (revision 3361) @@ -126,55 +126,2 @@ } -int edge_check (double *x1, double *y1, double *x2, double *y2) { - - double theta1, theta2; - double Theta1, Theta2; - - theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); - theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); - - if (theta1*theta2 < 0.0) { - return (FALSE); - } - - if (fabs(theta1) < fabs(theta2)) { - return (FALSE); - } - - Theta1 = theta1; - Theta2 = theta2; - theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); - theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); - - - if (theta1*theta2 < 0.0) { - return (FALSE); - } - - if (fabs(theta1) < fabs(theta2)) { - return (FALSE); - } - - return (TRUE); - -} - -double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) { - - double dx1, dy1, dx2, dy2, ct, st, theta; - - dx1 = x1 - x2; - dy1 = y1 - y2; - - dx2 = x3 - x2; - dy2 = y3 - y2; - - ct = (dx1*dx2 + dy1*dy2); - st = (dx1*dy2 - dx2*dy1); - - theta = atan2 (st, ct); - - return (theta); - -} - Index: /trunk/Ohana/src/addstar/src/greference.c =================================================================== --- /trunk/Ohana/src/addstar/src/greference.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/greference.c (revision 3361) @@ -1,5 +1,5 @@ # include "addstar.h" -Stars *grefcat (char *Refcat, CatStats *catstats, int *Nrefcat) { +Stars *grefcat (char *Refcat, GSCRegion *catstats, int *nstars) { int Nstars; @@ -7,7 +7,7 @@ if (VERBOSE) fprintf (stderr, "loading reference catalog data from %s\n", Refcat); + 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]); 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 */ @@ -30,5 +30,5 @@ exit (1); } - *Nrefcat = Nstars; + *nstars = Nstars; return (stars); Index: /trunk/Ohana/src/addstar/src/grefstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/grefstars.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/grefstars.c (revision 3361) @@ -1,5 +1,3 @@ # include "addstar.h" -# include -# include /* read ASCII file with ref star data */ @@ -38,7 +36,7 @@ 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; - stars[N].code = thiscode[0].code;; + stars[N].t = 0; + stars[N].code = thiscode[0].code; + stars[N].found = FALSE; CHECK_REALLOCATE (stars, Stars, NSTARS, N+1, 100); } Index: /trunk/Ohana/src/addstar/src/gregion_star.c =================================================================== --- /trunk/Ohana/src/addstar/src/gregion_star.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/gregion_star.c (revision 3361) @@ -1,5 +1,5 @@ # include "addstar.h" -GSCRegion *gregion_stars (Stars *stars, int Nstars, int *Nregion) { +GSCRegion *gregion_stars (Stars *stars, int Nstars, int *nregions) { int i; @@ -11,5 +11,5 @@ Nregions = 0; NREGIONS = 100; - ALLOCATE (region, GSCRegion, NREGIONS); + ALLOCATE (regions, GSCRegion, NREGIONS); sort_stars (stars, Nstars); @@ -27,6 +27,8 @@ Nregions ++; CHECK_REALLOCATE (regions, GSCRegion, Nregions, NREGIONS, 100); + } - + *nregions = Nregions; + return (regions); } Index: /trunk/Ohana/src/addstar/src/gstars.c =================================================================== --- /trunk/Ohana/src/addstar/src/gstars.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/gstars.c (revision 3361) @@ -3,6 +3,4 @@ # define BYTES_STAR 66 # define BLOCK 1000 -#include -#include Stars *gstars (char *file, int *NSTARS, Image *image) { @@ -142,5 +140,11 @@ ALLOCATE (stars, Stars, image[0].nstar); - /* temporary: get sky background from Flips data */ + /* this is really lame. get sky background from Flips data */ + /*** the problems are: + 1) need to define sky entry in image table + 2) need to have sky measurement (in reg.db) + (could be derived from star data?) + ***/ + if (SKYPROBE) { char *p; @@ -205,5 +209,9 @@ dparse (&stars[N].Y, 2, &buffer[j*BYTES_STAR]); dparse (&stars[N].M, 3, &buffer[j*BYTES_STAR]); - dparse (&stars[N].dM, 4, &buffer[j*BYTES_STAR]); + + /* cmp files carry dM in millimags */ + dparse (&tmp, 4, &buffer[j*BYTES_STAR]); + stars[N].dM = 0.001*tmp; + /* allow for some dynamic filtering of star list */ if (stars[N].dM > 1000.0 / SNLIMIT) { N--; continue; } @@ -215,4 +223,5 @@ dparse (&tmp, 5, &buffer[j*BYTES_STAR]); stars[N].dophot = tmp; + dparse (&stars[N].Mgal, 7, &buffer[j*BYTES_STAR]); dparse (&stars[N].Map, 8, &buffer[j*BYTES_STAR]); @@ -225,4 +234,5 @@ while (stars[N].R >= 360.0) stars[N].R -= 360.0; stars[N].found = -1; + stars[N].code = thiscode[0].code; if (SKYPROBE) { Index: /trunk/Ohana/src/addstar/src/image-db.c =================================================================== --- /trunk/Ohana/src/addstar/src/image-db.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/image-db.c (revision 3361) @@ -1,19 +1,21 @@ # include "addstar.h" -/* these variables are needed by Shutdown */ +/* file-pointer to Image.db table. if image db does not yet exist, + +*/ static FILE *f = (FILE *) NULL; +static int dbstate = LCK_UNLOCK; + +FILE *GetDB (int *state) { + *state = dbstate; + return (f); +} /* 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); } @@ -29,12 +31,14 @@ } -void unlock_image_db () { +void unlock_image_db (Image *image) { + /* protect wimage from interrupt signals */ - SetProtect (TRUE); - wimage (f, dbstate, &image, Nstars); - SetProtect (FALSE); + if (mode == M_IMAGE) { + SetProtect (TRUE); + wimage (image); + 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/in_image.c =================================================================== --- /trunk/Ohana/src/addstar/src/in_image.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/in_image.c (revision 3361) @@ -0,0 +1,14 @@ +# include "addstar.h" + +int in_image (double r, double d, Image *image) { + + double X, Y; + + RD_to_XY (&X, &Y, r, d, &image[0].coords); + if (X < 0) return (FALSE); + if (Y < 0) return (FALSE); + if (X >= image[0].NX) return (FALSE); + if (Y >= image[0].NY) return (FALSE); + return (TRUE); +} + Index: unk/Ohana/src/addstar/src/load_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/load_catalog.c (revision 3360) +++ (revision ) @@ -1,30 +1,0 @@ -# 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_pt_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/load_pt_catalog.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/load_pt_catalog.c (revision 3361) @@ -0,0 +1,24 @@ +# include "addstar.h" + +int load_pt_catalog (Catalog *catalog, GSCRegion *region) { + + /* don't free region before catalog! */ + catalog[0].filename = region[0].filename; + if (VERBOSE) fprintf (stderr, "adding to %s\n", region[0].filename); + + switch (lock_catalog (catalog, LCK_XCLD)) { + case 0: + /* catalog file is locked. this is an unexpected error because + ImageCat was NOT locked so we should be the only locking entity */ + fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); + exit (1); + case 1: + gcatalog (catalog); /* load from disk */ + break; + case 2: + if (ONLY_MATCH) return (FALSE); + mkcatalog (region, catalog); /* fills in new header info */ + break; + } + return (TRUE); +} Index: /trunk/Ohana/src/addstar/src/make_backup.c =================================================================== --- /trunk/Ohana/src/addstar/src/make_backup.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/make_backup.c (revision 3361) @@ -0,0 +1,23 @@ +# include "addstar.h" + +/* uses cp only */ +void make_backup (char *filename) { + + int status, cmode; + struct stat filestat; + char line [256]; + + status = stat (filename, &filestat); + if (status == 0) { /* file exists, make backup copy */ + sprintf (line, "cp %s %s~", filename, filename); + status = system (line); + if (status) { + fprintf (stderr, "ERROR: unable to create %s~, exiting\n", filename); + exit (1); + } + cmode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + sprintf (line, "%s~", filename); + chmod (line, cmode); + } + +} Index: /trunk/Ohana/src/addstar/src/mkcatalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/mkcatalog.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/mkcatalog.c (revision 3361) @@ -1,4 +1,3 @@ # include "addstar.h" -# include void mkcatalog (GSCRegion *region, Catalog *catalog) { Index: /trunk/Ohana/src/addstar/src/opening_angle.c =================================================================== --- /trunk/Ohana/src/addstar/src/opening_angle.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/opening_angle.c (revision 3361) @@ -0,0 +1,20 @@ +# include "addstar.h" + +double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) { + + double dx1, dy1, dx2, dy2, ct, st, theta; + + dx1 = x1 - x2; + dy1 = y1 - y2; + + dx2 = x3 - x2; + dy2 = y3 - y2; + + ct = (dx1*dx2 + dy1*dy2); + st = (dx1*dy2 - dx2*dy1); + + theta = atan2 (st, ct); + + return (theta); + +} Index: /trunk/Ohana/src/addstar/src/replace_match.c =================================================================== --- /trunk/Ohana/src/addstar/src/replace_match.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/replace_match.c (revision 3361) @@ -0,0 +1,18 @@ +# include "addstar.h" + +int replace_match (Average *average, Measure *measure, Stars *star) { + + if (!REPLACE) return (FALSE); + + /** in replace mode, search for entry and replace values M, dM, R, D */ + for (k = 0; k < average[0].Nm; k++) { + if (measure[k].source != thiscode[0].code) continue; + measure[k].dR = 360000.0*(average[0].R - stars[0].R); + measure[k].dD = 360000.0*(average[0].D - stars[0].D); + measure[k].M = MIN (1000*stars[0].M, NO_MAG); + measure[k].dM = MIN (1000*stars[0].dM, NO_ERR); + stars[N].found = average[0].offset + k; + return (TRUE); + } + return (FALSE); +} Index: /trunk/Ohana/src/addstar/src/save_pt_catalog.c =================================================================== --- /trunk/Ohana/src/addstar/src/save_pt_catalog.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/save_pt_catalog.c (revision 3361) @@ -0,0 +1,10 @@ +# include "addstar.h" + +void save_pt_catalog (Catalog *catalog) { + + /* protect wcatalog from interrupt signals */ + SetProtect (TRUE); + if (!DUMP_MATCHES && !ONLY_IMAGES) wcatalog (catalog); + SetProtect (FALSE); + unlock_catalog (catalog); +} Index: /trunk/Ohana/src/addstar/src/update_coords.c =================================================================== --- /trunk/Ohana/src/addstar/src/update_coords.c (revision 3361) +++ /trunk/Ohana/src/addstar/src/update_coords.c (revision 3361) @@ -0,0 +1,53 @@ +# include "addstar.h" + +/* find the average and scatter for R and D - no outlier rejection */ +void update_coords (Average *average, Measure *measure, int *next) { + + double R, D, r, d; + double r2, d2, dR2, dD2; + int Npt; + + Npt = r = d = r2 = d2 = 0; + + if (average[0].Nm < 2) return; + + /* find the average & sum-square */ + m = average[0].offset; /* first measurement of this star */ + for (k = 0; k < average[0].Nm; k++) { + if (measure[m].t == 0) { + m = next[m]; + continue; + } + R = measure[m].dR; + D = measure[m].dD; + r += R; + d += D; + r2 += R*R; + d2 += D*D; + m = next[m]; + Npt ++; + } + if (Npt < 1) return; + + /* apply average offset */ + r = r / Npt; /* these are corrections in 1/100 arcsec to RA and DEC */ + d = d / Npt; + average[0].R -= r / 360000.0; + average[0].D -= d / 360000.0; + m = average[0].offset; /* first measurement of this star */ + for (k = 0; k < average[0].Nm; k++) { + measure[m].dR -= r; + measure[m].dD -= d; + m = next[m]; + } + + /* measure scatter, if possible */ + if (Npt < 2) return; + + dR2 = r2 / Npt - r*r; + dD2 = d2 / Npt - d*d; + average[0].Xp = sqrt (dD2 + dR2 / SQ(cos(d*RAD_DEG))); + /* Xp is scatter in position in hundredths of arcsec */ + + return; +} Index: /trunk/Ohana/src/addstar/src/wimage.c =================================================================== --- /trunk/Ohana/src/addstar/src/wimage.c (revision 3360) +++ /trunk/Ohana/src/addstar/src/wimage.c (revision 3361) @@ -1,10 +1,12 @@ # include "addstar.h" -void wimage (FILE *f, int dbstate, Image *image, int Nstars) { +/* add image to image db */ +void wimage (Image *image) { - int Nimages, status, offset; + FILE *f; + int Nimages, status, offset, dbstate; Header header; - /* place file pointer at beginning of file */ + f = GetDB (&dbstate); fseek (f, 0, SEEK_SET); @@ -17,4 +19,5 @@ /* assign Zero Point provided by first image in dB */ fits_modify (&header, "ZERO_PT", "%lf", 1, ZeroPt); + fits_modify (&header, "NIMAGES", "%d", 1, 0); } else { if (!fits_load_header (f, &header)) {