Changeset 3361
- Timestamp:
- Mar 2, 2005, 7:54:17 AM (21 years ago)
- Location:
- trunk/Ohana/src/addstar
- Files:
-
- 8 added
- 2 deleted
- 24 edited
-
Makefile (modified) (3 diffs)
-
include/addstar.h (modified) (8 diffs)
-
src/RegionOps.c (modified) (1 diff)
-
src/addstar-old.c (deleted)
-
src/addstar.c (modified) (3 diffs)
-
src/airmass.c (modified) (2 diffs)
-
src/args.c (modified) (9 diffs)
-
src/calibrate.c (modified) (7 diffs)
-
src/check_permissions.c (modified) (1 diff)
-
src/conversions.c (modified) (1 diff)
-
src/edge_check.c (added)
-
src/find_matches.c (modified) (13 diffs)
-
src/find_matches_refstars.c (modified) (9 diffs)
-
src/find_subset.c (modified) (2 diffs)
-
src/gcatalog.c (modified) (1 diff)
-
src/get2mass.c (modified) (2 diffs)
-
src/getgsc.c (modified) (2 diffs)
-
src/getusno.c (modified) (5 diffs)
-
src/gimages.c (modified) (1 diff)
-
src/greference.c (modified) (3 diffs)
-
src/grefstars.c (modified) (2 diffs)
-
src/gregion_star.c (modified) (3 diffs)
-
src/gstars.c (modified) (5 diffs)
-
src/image-db.c (modified) (2 diffs)
-
src/in_image.c (added)
-
src/load_catalog.c (deleted)
-
src/load_pt_catalog.c (added)
-
src/make_backup.c (added)
-
src/mkcatalog.c (modified) (1 diff)
-
src/opening_angle.c (added)
-
src/replace_match.c (added)
-
src/save_pt_catalog.c (added)
-
src/update_coords.c (added)
-
src/wimage.c (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Ohana/src/addstar/Makefile
r3339 r3361 1 default: addstar 2 help: 3 @echo "make options: addstar (default)" 4 1 5 include ../../Configure 6 2 7 HOME = $(ROOT)/src/addstar 3 PROGRAM = addstar4 5 8 BIN = $(HOME)/bin 6 9 INC = $(HOME)/include … … 15 18 INCS = -I$(INC) -I$(LINC) -I$(XINC) 16 19 LIBS = -L$(LLIB) -lFITS -lohana -lm 17 CFLAGS = -o $*.$(ARCH).o$(INCS)18 CCFLAGS = $(INCS)$(LIBS)20 CFLAGS = $(INCS) 21 LFLAGS = $(LIBS) 19 22 20 23 ADDSTAR = \ 21 $(SRC)/addstar.$(ARCH).o $(SRC)/gcatalog.$(ARCH).o \ 22 $(SRC)/gregion_image.$(ARCH).o $(SRC)/gregion_star.$(ARCH).o \ 23 $(SRC)/gimages.$(ARCH).o $(SRC)/wimage.$(ARCH).o \ 24 $(SRC)/gstars.$(ARCH).o $(SRC)/grefstars.$(ARCH).o \ 25 $(SRC)/find_matches.$(ARCH).o $(SRC)/find_matches_refstars.$(ARCH).o \ 26 $(SRC)/wcatalog.$(ARCH).o $(SRC)/conversions.$(ARCH).o \ 27 $(SRC)/sort_lists.$(ARCH).o $(SRC)/ConfigInit.$(ARCH).o \ 28 $(SRC)/aregion.$(ARCH).o $(SRC)/find_subset.$(ARCH).o \ 29 $(SRC)/load_subpix.$(ARCH).o $(SRC)/airmass.$(ARCH).o \ 30 $(SRC)/mkcatalog.$(ARCH).o $(SRC)/calibrate.$(ARCH).o \ 31 $(SRC)/args.$(ARCH).o $(SRC)/parse_time.$(ARCH).o \ 32 $(SRC)/match_refstars.$(ARCH).o \ 33 $(SRC)/check_permissions.$(ARCH).o 24 $(SRC)/addstar.$(ARCH).o \ 25 $(SRC)/airmass.$(ARCH).o \ 26 $(SRC)/aregion.$(ARCH).o \ 27 $(SRC)/args.$(ARCH).o \ 28 $(SRC)/calibrate.$(ARCH).o \ 29 $(SRC)/check_permissions.$(ARCH).o \ 30 $(SRC)/ConfigInit.$(ARCH).o \ 31 $(SRC)/conversions.$(ARCH).o \ 32 $(SRC)/edge_check.$(ARCH).o \ 33 $(SRC)/find_matches.$(ARCH).o \ 34 $(SRC)/find_matches_refstars.$(ARCH).o \ 35 $(SRC)/find_subset.$(ARCH).o \ 36 $(SRC)/gcatalog.$(ARCH).o \ 37 $(SRC)/get2mass.$(ARCH).o \ 38 $(SRC)/getgsc.$(ARCH).o \ 39 $(SRC)/getusno.$(ARCH).o \ 40 $(SRC)/gimages.$(ARCH).o \ 41 $(SRC)/greference.$(ARCH).o \ 42 $(SRC)/grefstars.$(ARCH).o \ 43 $(SRC)/gregion_image.$(ARCH).o \ 44 $(SRC)/gregion_match.$(ARCH).o \ 45 $(SRC)/gregion_patch.$(ARCH).o \ 46 $(SRC)/gregion_star.$(ARCH).o \ 47 $(SRC)/gstars.$(ARCH).o \ 48 $(SRC)/image-db.$(ARCH).o \ 49 $(SRC)/in_image.$(ARCH).o \ 50 $(SRC)/load_pt_catalog.$(ARCH).o \ 51 $(SRC)/load_subpix.$(ARCH).o \ 52 $(SRC)/make_backup.$(ARCH).o \ 53 $(SRC)/mkcatalog.$(ARCH).o \ 54 $(SRC)/opening_angle.$(ARCH).o \ 55 $(SRC)/parse_time.$(ARCH).o \ 56 $(SRC)/RegionOps.$(ARCH).o \ 57 $(SRC)/replace_match.$(ARCH).o \ 58 $(SRC)/save_pt_catalog.$(ARCH).o \ 59 $(SRC)/SetSignals.$(ARCH).o \ 60 $(SRC)/sort_lists.$(ARCH).o \ 61 $(SRC)/update_coords.$(ARCH).o \ 62 $(SRC)/wcatalog.$(ARCH).o \ 63 $(SRC)/wimage.$(ARCH).o 34 64 35 OBJ = $(ADDSTAR) 65 EXTRA = \ 66 $(SRC)/bracket.$(ARCH).o \ 67 $(SRC)/find_proper.$(ARCH).o \ 68 $(SRC)/gregion_match_glob.$(ARCH).o \ 69 $(SRC)/gregions.$(ARCH).o \ 70 $(SRC)/match_refstars.$(ARCH).o \ 36 71 37 72 default: $(PROGRAM) 38 73 39 $(ADDSTAR): $(INC)/addstar.h 74 addstar : $(BIN)/addstar.$(ARCH) 40 75 41 # dependancy rules for binary code ########################## 42 $(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) 76 $(BIN)/addstar.$(ARCH) : $(ADDSTAR) 43 77 44 $(BIN)/$(PROGRAM).$(ARCH): $(OBJ) 78 # $(ADDSTAR): $(INC)/addstar.h 79 80 INSTALL = addstar 81 82 # dependancy rules for binary code ######################### 83 .PRECIOUS: %.$(ARCH).o 84 .PRECIOUS: $(BIN)/%.$(ARCH) 85 86 %.$(ARCH).o : %.c 87 $(CC) $(CFLAGS) -c $< -o $@ 88 89 $(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o 45 90 @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi 46 $(CC) $ (OBJ) -o $(BIN)/$(PROGRAM).$(ARCH) $(CCFLAGS)91 $(CC) $^ -o $@ $(LFLAGS) 47 92 48 install: $(DESTBIN)/$(PROGRAM) 93 $(DESTBIN)/%: $(BIN)/%.$(ARCH) 94 @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi 95 rm -f $(DESTBIN)/$* 96 cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* 49 97 50 $(DESTBIN)/$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) 51 @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi 52 rm -f $(DESTBIN)/$(PROGRAM) 53 cp $(BIN)/$(PROGRAM).$(ARCH) $(DESTBIN)/$(PROGRAM) 98 $(INSTALL) $(DEVEL): % : $(BIN)/%.$(ARCH) 99 100 %.clean : 101 rm -f $(BIN)/$*.$(ARCH) 102 103 %.install: 104 make $(DESTBIN)/$* 54 105 55 106 # utilities ################################################# 107 108 install: 109 for i in $(INSTALL); do make $$i.install; done 110 56 111 clean: 57 112 rm -f $(BIN)/*.$(ARCH) … … 59 114 rm -f `find . -name "*~"` 60 115 rm -f `find . -name "#*"` 61 62 .SUFFIXES: .$(ARCH).o63 64 .c.$(ARCH).o:65 $(CC) $(CFLAGS) -c $<66 67 -
trunk/Ohana/src/addstar/include/addstar.h
r3347 r3361 2 2 # include <loneos.h> 3 3 # include <signal.h> 4 # include <sys/time.h> 5 # include <time.h> 6 7 /* used in find_matches, find_matches_refstars */ 8 # define IN_CATALOG(R,D) ( \ 9 ((D) >= region[0].DEC[0]) && ((D) < region[0].DEC[1]) && \ 10 ((R) >= region[0].RA[0]) && ((R) < region[0].RA[1])) 11 12 /* grab named photcode */ 13 # define NAMED_PHOTCODE(CODE,NAME) \ 14 CODE = GetPhotcodebyName (NAME); \ 15 if (!CODE) { \ 16 fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", NAME); \ 17 exit (0); } 4 18 5 19 typedef struct { … … 14 28 double Mgal, Map; 15 29 int found; 30 short int code; 31 e_time t; 16 32 } Stars; 33 34 /* structure for data on a catalog region */ 35 typedef struct { 36 char filename[256]; 37 double RA[2]; 38 int Nrec; 39 } TM_Region; 17 40 18 41 enum {M_IMAGE, M_REFLIST, M_REFCAT}; … … 34 57 35 58 PhotCode *thiscode; 59 GSCRegion patch; 36 60 37 61 char DateKeyword[64]; … … 57 81 int ONLY_MATCH; 58 82 83 int MODE; 84 59 85 time_t TIMEREF; 60 86 … … 70 96 void gregion_star (Stars *star, GSCRegion *region); 71 97 GSCRegion *gregion_image (Image *image, int *Nregions); 72 Image *gimages (FILE *f, Image *image, int *Npimage);73 98 int edge_check (double *x1, double *y1, double *x2, double *y2); 74 99 double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); 75 void wimage (FILE *f, int dbstate, Image *image, int Nstars);76 100 Stars *gstars (char *file, int *NSTARS, Image *image); 77 void find_matches_refstars (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog);78 void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage);79 101 int in_image (double r, double d, Image *image); 80 102 void wcatalog (Catalog *catalog); 81 103 void ConfigInit (int *argc, char **argv); 82 104 void aregion (GSCRegion *region, FILE *f, double ra, double dec); 83 void wimage (FILE *f, int dbstate, Image *image, int Nstars); 84 double airmass (double ra, double dec, double st, double latitude); 105 short airmass (short secz_image, double ra, double dec, double st, double latitude); 85 106 void InitCalibration (); 86 void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y);87 107 void FindCalibration (Image *image); 88 108 double get_subpix (double x, double y); … … 91 111 void mkcatalog (GSCRegion *region, Catalog *catalog); 92 112 Stars *grefstars (char *file, int *Nstars); 93 Stars * find_subset (GSCRegion *region, Stars *stars, int Nstars, int *NSTARS);113 Stars **find_subset (GSCRegion *region, Stars *stars, int Nstars, int *NSTARS); 94 114 int match_refstars (Stars *stars, int Nstars); 95 115 int parse_time (Header *header); … … 113 133 time_t jd_to_sec (double jd); 114 134 time_t date_to_sec (char *date); 135 time_t short_date_to_sec (char *date); 136 115 137 116 138 /** … … 122 144 e_time pointer from one of these functions) 123 145 **/ 146 147 Stars *get2mass (GSCRegion *patch, int *NSTARS); 148 TM_Region *get2mass_acc (GSCRegion *patch, int *nfiles); 149 Stars *get2mass_data (char *filename, GSCRegion *patch, int *nstars); 150 GSCRegion *gregion_stars (Stars *stars, int Nstars, int *Nregion); 151 Stars *grefstars (char *file, int *Nstars); 152 GSCRegion *LoadRegions (int *nregions); 153 int FindRegionByPoint (GSCRegion *regions, double ra, double dec, GSCRegion *gsc, int Ngsc); 154 int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec); 155 int FindRegionDecBandStop (GSCRegion *region, int Nregion, double dec); 156 int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string); 157 void find_matches_refstars (GSCRegion *region, Stars **stars, int Nstars, Catalog *catalog); 158 Stars *gstars (char *file, int *NSTARS, Image *image); 159 void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap); 160 Stars *getusno (GSCRegion *catstats, int *Nstars); 161 Stars *rd_gsc (char *filename, int *nstars); 162 Stars *getgsc (GSCRegion *patch, int *NSTARS); 163 Stars *grefcat (char *Refcat, GSCRegion *catstats, int *nstars); 164 void load_subpix (); 165 double get_subpix (double x, double y); 166 double scat_subpix (double x, double y); 167 void help (); 168 int args (int argc, char **argv); 169 void InitCalibration (); 170 void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N); 171 void AddToCalibration (Average *average, Measure *measure, Measure *new, int Nstar); 172 void FindCalibration (Image *image); 173 int load_pt_catalog (Catalog *catalog, GSCRegion *region); /*** choose new name ***/ 174 void save_pt_catalog (Catalog *catalog); /*** choose new name ***/ 175 void wcatalog (Catalog *catalog); 176 GSCRegion *gregions (Image *image, int *Nregions); 177 int gcatalog (Catalog *catalog); 178 GSCRegion *gregion_image (Image *image, int *Nregions); 179 void update_coords (Average *average, Measure *measure, int *next); 180 int replace_match (Average *average, Measure *measure, Stars *star); 181 FILE *GetDB (int *state); 182 int Shutdown (); 183 void lock_image_db (); 184 void unlock_image_db (Image *image); 185 void wimage (Image *image); 186 Image *gimages (Image *image, int *Npimage); 187 double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3); 188 int edge_check (double *x1, double *y1, double *x2, double *y2); 189 GSCRegion *gregion_match (GSCRegion *regions, int *nregions); 190 GSCRegion *gregion_patch (GSCRegion *patch, int *nregions); -
trunk/Ohana/src/addstar/src/RegionOps.c
r3347 r3361 177 177 */ 178 178 179 /**********/180 int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) {181 182 int flag_d0, flag_d1, flag_h0, flag_h1;183 double tmp;184 185 *d0 = *h0 = *d1 = *h1 = 0;186 187 flag_h0 = dparse (h0, 1, string);188 flag_h1 = dparse (h1, 4, string);189 flag_d0 = dparse (d0, 7, string);190 flag_d1 = dparse (d1, 9, string);191 *h0 *= flag_h0;192 *h1 *= flag_h1;193 *d0 *= flag_d0;194 *d1 *= flag_d1;195 196 dparse (&tmp, 2, string);197 *h0 += tmp/60.0;198 dparse (&tmp, 3, string);199 *h0 += tmp/3600.0;200 201 dparse (&tmp, 5, string);202 *h1 += tmp/60.0;203 dparse (&tmp, 6, string);204 *h1 += tmp/3600.0;205 206 dparse (&tmp, 8, string);207 *d0 += tmp/60.0;208 209 dparse (&tmp, 10, string);210 *d1 += tmp/60.0;211 212 *h0 *= 15*flag_h0;213 *h1 *= 15*flag_h1;214 *d0 *= flag_d0;215 *d1 *= flag_d1;216 217 return (TRUE);218 }219 220 221 179 /* 222 180 if (buffer[i*48 + 19] == ' ') continue; -
trunk/Ohana/src/addstar/src/addstar.c
r3347 r3361 3 3 int main (int argc, char **argv) { 4 4 5 int i, mode, dbstate; 6 int Nstars, Nimage, Nregions, Nmissed; 7 Stars *stars; 8 Image image, *image_db; 9 GSCRegion *region; 5 int i, Nstars, Noverlap, Nregions, Nsubset; 6 Stars *stars, **subset; 7 Image image, *overlap; 8 GSCRegion *regions; 10 9 Catalog catalog; 11 10 … … 14 13 args (argc, argv); 15 14 15 stars = NULL; 16 overlap = NULL; 17 regions = NULL; 18 16 19 if (SKYPROBE) load_subpix (); 17 20 18 21 lock_image_db (); 19 22 20 switch ( mode) {23 switch (MODE) { 21 24 case M_IMAGE: 22 25 stars = gstars (argv[1], &Nstars, &image); 23 region = gregion_image (&image, &Nregions);26 regions = gregion_image (&image, &Nregions); 24 27 overlap = gimages (&image, &Noverlap); 25 28 break; 26 29 case M_REFLIST: 27 30 stars = grefstars (argv[1], &Nstars); 28 region = gregion_stars (stars, &Nregion);31 regions = gregion_stars (stars, Nstars, &Nregions); 29 32 break; 30 33 case M_REFCAT: 31 regions = gregion_patch (patch, &Nregions); 32 if (ONLY_MATCH) { 33 regions = gregion_match (regions, &Nregions); 34 } 34 regions = gregion_patch (&patch, &Nregions); 35 35 break; 36 } 37 if (ONLY_MATCH) { 38 regions = gregion_match (regions, &Nregions); 36 39 } 37 40 38 41 for (i = 0; i < Nregions; i++) { 39 load_catalog (&catalog, ®ion[i]);42 if (!load_pt_catalog (&catalog, ®ions[i])) continue; 40 43 41 switch ( mode) {44 switch (MODE) { 42 45 case M_IMAGE: 43 find_matches (®ion [i], stars, Nstars, &catalog, &image, overlap, Noverlap);46 find_matches (®ions[i], stars, Nstars, &catalog, &image, overlap, Noverlap); 44 47 break; 45 48 case M_REFCAT: 46 stars = grefcat (argv[1], region , &Nstars);49 stars = grefcat (argv[1], regions, &Nstars); 47 50 case M_REFLIST: 48 subset = find_subset (®ion, stars, Nstars, &subref, &Nsubset); 49 find_matches_refstars (®ion, subset, Nsubset, &catalog); 51 subset = find_subset (®ions[i], stars, Nstars, &Nsubset); 52 find_matches_refstars (®ions[i], subset, Nsubset, &catalog); 53 if (Nsubset) free (subset); 50 54 break; 51 55 } 52 wcatalog (&catalog);56 save_pt_catalog (&catalog); 53 57 } 54 58 55 59 if (CALIBRATE) { FindCalibration (&image); } 56 if (DUMP_MATCHES) Shutdown (); 57 58 if (mode == IMAGE) wimage (f, dbstate, &image, Nstars); 59 unlock_image_db (); 60 unlock_image_db (&image); 60 61 exit (0); 61 62 } … … 64 65 65 66 catalog - existing object db table 66 region - sky area which may or may contain data67 regions - sky area which may or may contain data 67 68 patch - RA,DEC bounded portion of sky 68 69 -
trunk/Ohana/src/addstar/src/airmass.c
r3339 r3361 1 1 # include "addstar.h" 2 2 3 double airmass (double ra, double dec, double st, double latitude) {3 short airmass (short secz_image, double ra, double dec, double st, double latitude) { 4 4 5 5 double hour, cosz, secz; 6 7 if (!SKYPROBE) return (secz_image); 6 8 7 9 /*** make this optional? we may not have ST... ***/ … … 12 14 cosz = sin (RAD_DEG*dec) * sin (RAD_DEG*latitude) + cos (RAD_DEG*dec) * cos (RAD_DEG*hour) * cos (RAD_DEG*latitude); 13 15 14 secz = 1.0 / cosz;15 16 /* modify to match storage units (millimags) */ 17 secz = 1000.0 / cosz; 16 18 return (secz); 17 18 19 } -
trunk/Ohana/src/addstar/src/args.c
r3347 r3361 10 10 fprintf (stderr, "\n"); 11 11 exit (2); 12 13 12 } 14 13 … … 26 25 27 26 /* basic mode: image, list, refcat */ 28 mode= M_IMAGE;27 MODE = M_IMAGE; 29 28 if ((N = get_argument (argc, argv, "-ref"))) { 30 mode= M_REFLIST;29 MODE = M_REFLIST; 31 30 remove_argument (N, &argc, argv); 32 31 } 33 32 if ((N = get_argument (argc, argv, "-cat"))) { 34 mode= M_REFCAT;33 MODE = M_REFCAT; 35 34 remove_argument (N, &argc, argv); 36 35 } 37 36 37 /*** provide additional data ***/ 38 38 /* restrict to a portion of the sky? (REFCAT only) */ 39 39 patch.RA[0] = 0; … … 48 48 patch.DEC[1] = atof (argv[N]); 49 49 } 50 /* override any header PHOTCODE values */ 51 thiscode = NULL; 52 if ((N = get_argument (argc, argv, "-p"))) { 53 remove_argument (N, &argc, argv); 54 thiscode = GetPhotcodebyName (argv[N]); 55 remove_argument (N, &argc, argv); 56 } 57 /* provide a time for dataset */ 58 TIMEREF = 0; 59 if ((i = get_argument (argc, argv, "-time"))) { 60 remove_argument (i, &argc, argv); 61 if (!str_to_time (argv[i], &TIMEREF)) { 62 fprintf (stderr, "syntax error in time\n"); 63 exit (1); 64 } 65 remove_argument (i, &argc, argv); 66 } 50 67 51 /* only add to existing objects (REFCAT only) */ 68 /*** modify behavior ***/ 69 /* only add to existing objects */ 52 70 ONLY_MATCH = FALSE; 53 71 if ((i = get_argument (argc, argv, "-only-match"))) { … … 55 73 remove_argument (i, &argc, argv); 56 74 } 57 /* don't add missed pts to Missed table */75 /* don't add missed pts to Missed table (image only) */ 58 76 SKIP_MISSED = FALSE; 59 77 if ((N = get_argument (argc, argv, "-missed"))) { … … 61 79 remove_argument (N, &argc, argv); 62 80 } 63 /* replace measurement, don't duplicate */81 /* replace measurement, don't duplicate (ref/cat only) */ 64 82 REPLACE = FALSE; 65 83 if ((i = get_argument (argc, argv, "-replace"))) { … … 73 91 remove_argument (N, &argc, argv); 74 92 } 75 /* apply average zpt offset calibration */93 /* apply average zpt offset calibration (image only) */ 76 94 CALIBRATE = FALSE; 77 95 if ((N = get_argument (argc, argv, "-cal"))) { 78 96 CALIBRATE = TRUE; 79 97 remove_argument (N, &argc, argv); 80 CalReference = GetPhotcodeCodebyName (argv[N]);81 if (!CalReference) {82 fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]);83 exit (1);84 }85 remove_argument (N, &argc, argv);86 CalColor = GetPhotcodeCodebyName (argv[N]);87 if (!CalColor) {88 fprintf (stderr, "ERROR: photcode color %s not found in photcode table\n", argv[N]);89 exit (1);90 }91 remove_argument (N, &argc, argv);92 98 } 93 99 100 /*** optional situations ***/ 94 101 /* treat data specially for skyprobe (calibration, subpix) */ 95 102 SKYPROBE = FALSE; … … 98 105 remove_argument (N, &argc, argv); 99 106 } 107 /* accept bad header astrometry */ 100 108 ACCEPT_ASTROM = FALSE; 101 109 if ((N = get_argument (argc, argv, "-accept"))) { … … 103 111 remove_argument (N, &argc, argv); 104 112 } 105 DUMP_MATCHES = FALSE; 106 if ((N = get_argument (argc, argv, "-dump"))) { 107 DUMP_MATCHES = TRUE; 108 remove_argument (N, &argc, argv); 109 } 113 /* force read of image database with mismatched NSTARS & size */ 110 114 FORCE_READ = FALSE; 111 115 if ((N = get_argument (argc, argv, "-force"))) { … … 113 117 remove_argument (N, &argc, argv); 114 118 } 115 119 /* extra error messages */ 116 120 VERBOSE = FALSE; 117 121 if ((N = get_argument (argc, argv, "-v"))) { 118 122 VERBOSE = TRUE; 119 123 remove_argument (N, &argc, argv); 120 }121 thiscode = NULL;122 if ((N = get_argument (argc, argv, "-p"))) {123 remove_argument (N, &argc, argv);124 thiscode = GetPhotcodebyName (argv[N]);125 remove_argument (N, &argc, argv);126 }127 TIMEREF = 0;128 if ((i = get_argument (argc, argv, "-time"))) {129 remove_argument (i, &argc, argv);130 if (!str_to_time (argv[i], &TIMEREF)) {131 fprintf (stderr, "syntax error in time\n");132 exit (1);133 }134 remove_argument (i, &argc, argv);135 }136 if (DUMP_MATCHES && !CALIBRATE) {137 fprintf (stderr, "-dump only valid -with -cal\n");138 exit (1);139 124 } 140 125 -
trunk/Ohana/src/addstar/src/calibrate.c
r3347 r3361 1 1 # include "addstar.h" 2 2 3 static int Ncal, NCAL ;3 static int Ncal, NCAL, *Nstar; 4 4 static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm; 5 static float *RA, *DEC, *X, *Y;6 int *Nstar;7 5 8 6 void InitCalibration () { 9 7 10 fprintf (stderr, "calibrating the image... %d\n", CalReference);8 fprintf (stderr, "calibrating the image... %d\n", thiscode[0].equiv); 11 9 Ncal = 0; 12 10 NCAL = 1000; … … 18 16 ALLOCATE (Airm, float, NCAL); 19 17 ALLOCATE (Nstar, int, NCAL); 20 ALLOCATE (RA, float, NCAL);21 ALLOCATE (DEC, float, NCAL);22 ALLOCATE (X, float, NCAL);23 ALLOCATE (Y, float, NCAL);24 18 } 25 19 26 void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N , float ra, float dec, float x, float y) {20 void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N) { 27 21 28 22 Mobs[Ncal] = M; … … 33 27 Airm[Ncal] = A; 34 28 Nstar[Ncal] = N; 35 36 RA[Ncal] = ra;37 DEC[Ncal] = dec;38 X[Ncal] = x;39 Y[Ncal] = y;40 29 Ncal ++; 41 30 … … 49 38 REALLOCATE (Airm, float, NCAL); 50 39 REALLOCATE (Nstar, int, NCAL); 51 REALLOCATE (RA, float, NCAL);52 REALLOCATE (DEC, float, NCAL);53 REALLOCATE (X, float, NCAL);54 REALLOCATE (Y, float, NCAL);55 40 } 41 } 42 43 /* this does not use the linked list to navigate the measure structure 44 this is OK for old measures, but will break for entries which have 45 already added a new measure to the list. */ 46 void AddToCalibration (Average *average, Measure *measure, Measure *new, int Nstar) { 47 48 int i, found0, found1, found2; 49 float CalM0, CalM1, CalM2, dCalM; 50 short CalC0, CalC1, CalC2; 51 52 found0 = found1 = found2 = FALSE; 53 CalM0 = CalM1 = CalM2 = dCalM = NO_MAG; 54 55 CalC0 = thiscode[0].equiv; 56 CalC1 = thiscode[0].c1; 57 CalC2 = thiscode[0].c2; 58 59 for (i = 0; i < average[0].Nm; i++) { 60 if (measure[0].source == CalC0) { 61 found0 = TRUE; 62 CalM0 = measure[i].M; 63 dCalM = measure[i].dM; 64 } 65 if (measure[i].source == CalC1) { 66 found1 = TRUE; 67 CalM1 = measure[i].M; 68 } 69 if (measure[i].source == CalC2) { 70 found2 = TRUE; 71 CalM2 = measure[i].M; 72 } 73 if (found0 && found1 && found2) { 74 SaveCalibration (new[0].M, new[0].dM, CalM0, dCalM, CalM1-CalM2, new[0].airmass, Nstar); 75 return; 76 } 77 /* n = next[n]; - if we use this, we need to pass in measure[0] above */ 78 } 79 return; 56 80 } 57 81 … … 59 83 60 84 int i, MaxN, *Nlist, Nkeep; 85 short int Mint; 61 86 float N, M1, M2, Klam, Clam, Xlam, Mabs, *Dmag, *dDmag; 62 87 float dMo, dMr, Mw, Dmed, W1, W2, NSigma; … … 91 116 /* if this entry has too many (or two few?) matches, skip it */ 92 117 if (Nlist[Nstar[i]] != 1) continue; 93 M abs = 0.001*(Mobs[i] + Klam*(Airm[i] - 1000) + Clam + Xlam*(Mref[i] - Cref[i])) - ZeroPt;94 118 Mint = Clam + Mobs[i] + Xlam*Cref[i] + Klam*(Airm[i] - 1000); 119 Mabs = 0.001*Mint - ZeroPt; 95 120 /* note: subpix correction is applied in gstars */ 96 97 if (DUMP_MATCHES)98 fprintf (stdout, "%d %6.3f %6.3f %6.3f %10.6f %10.6f %7.2f %7.2f %7.2f %7.2f\n",99 i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i], Airm[i], 0.001*dMobs[i]);100 121 101 122 /* skip stars brighter than 8.0 */ … … 157 178 } 158 179 } 159 160 # if (0)161 162 /* old correction - superceeded by subpix grid */163 164 /* find lo (Dmed - 0.1) + hi (Dmed + 0.1)165 lo = -1; hi = -1;166 for (i = 0; (lo == -1) && (i < Nkeep); i++) { if (Dmed - Dmag[i] < 0.1) lo = i; }167 for (i = Nkeep - 1; (hi == -1) && (i >= 0); i--) { if (Dmag[i] - Dmed < 0.1) hi = i; }168 if (lo == -1) lo = 0;169 if (hi == -1) hi = Nkeep - 1;170 */171 172 if (Nkeep < 3) {173 fprintf (stderr, "too few stars\n");174 image[0].Mcal = 10000;175 image[0].dMcal = 10000;176 return;177 }178 fsort2 (Dmag, dDmag, Nkeep);179 180 /* accumulate delta mags (25% - 75% of list) */181 M1 = 0.0;182 M2 = 0.0;183 N = 0.0;184 for (i = 0.25*Nkeep; i <= 0.75*Nkeep; i++) {185 186 /* a straight mean works better, given the bad photometry of skyprobe */187 M1 += Dmag[i];188 M2 += SQ (Dmag[i]);189 N += 1.0;190 191 }192 193 if (N > 1) {194 M1 = M1 / N;195 M2 = sqrt (fabs(M2/N - M1*M1));196 fprintf (stderr, "N: %.0f, mean: %f, stdev: %f, precision: %f\n", N, M1, M2, M2 / sqrt (N));197 image[0].Mcal = 1000 * M1;198 image[0].dMcal = 1000 * M2 / sqrt (N);199 } else {200 fprintf (stderr, "too few stars\n");201 image[0].Mcal = 10000;202 image[0].dMcal = 10000;203 }204 205 206 /* subpix correction : make parameters global! */207 dAs = 0;208 if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i];209 dYs = Y[i] - (int) (Y[i]) + 0.55;210 dMs = dAs * sin (6.28*dYs);211 /* Mabs -= dMs; */212 213 # endif -
trunk/Ohana/src/addstar/src/check_permissions.c
r2457 r3361 65 65 } 66 66 } 67 68 67 } 69 70 /* uses cp only */71 void make_backup (char *filename) {72 73 int status, cmode;74 struct stat filestat;75 char line [256];76 77 status = stat (filename, &filestat);78 if (status == 0) { /* file exists, make backup copy */79 sprintf (line, "cp %s %s~", filename, filename);80 status = system (line);81 if (status) {82 fprintf (stderr, "ERROR: unable to create %s~, exiting\n", filename);83 exit (1);84 }85 cmode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH;86 sprintf (line, "%s~", filename);87 chmod (line, cmode);88 }89 90 } -
trunk/Ohana/src/addstar/src/conversions.c
r3339 r3361 298 298 return (second); 299 299 } 300 301 /***** short date in format yymmdd *****/ 302 time_t short_date_to_sec (char *date) { 303 304 time_t second; 305 double jd; 306 struct tm now; 307 308 bzero (&now, sizeof(now)); 309 310 scanf (date, "%2d%2d%2d", &now.tm_year, &now.tm_mon, &now.tm_mday); 311 312 if (now.tm_year > 51) now.tm_year += 0; 313 if (now.tm_year < 50) now.tm_year += 100; 314 now.tm_mon --; /* tm_mon runs from 0 - 11 */ 315 316 jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4) 317 + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12) 318 - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5; 319 320 second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec; 321 322 return (second); 323 } 324 325 /**********/ 326 int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) { 327 328 int flag_d0, flag_d1, flag_h0, flag_h1; 329 double tmp; 330 331 *d0 = *h0 = *d1 = *h1 = 0; 332 333 flag_h0 = dparse (h0, 1, string); 334 flag_h1 = dparse (h1, 4, string); 335 flag_d0 = dparse (d0, 7, string); 336 flag_d1 = dparse (d1, 9, string); 337 *h0 *= flag_h0; 338 *h1 *= flag_h1; 339 *d0 *= flag_d0; 340 *d1 *= flag_d1; 341 342 dparse (&tmp, 2, string); 343 *h0 += tmp/60.0; 344 dparse (&tmp, 3, string); 345 *h0 += tmp/3600.0; 346 347 dparse (&tmp, 5, string); 348 *h1 += tmp/60.0; 349 dparse (&tmp, 6, string); 350 *h1 += tmp/3600.0; 351 352 dparse (&tmp, 8, string); 353 *d0 += tmp/60.0; 354 355 dparse (&tmp, 10, string); 356 *d1 += tmp/60.0; 357 358 *h0 *= 15*flag_h0; 359 *h1 *= 15*flag_h1; 360 *d0 *= flag_d0; 361 *d1 *= flag_d1; 362 363 return (TRUE); 364 } -
trunk/Ohana/src/addstar/src/find_matches.c
r3347 r3361 1 1 # include "addstar.h" 2 2 3 # define IN_CATALOG(R,D) ( \4 ((D) >= region[0].DEC[0]) && ((D) < region[0].DEC[1]) && \5 ((R) >= region[0].RA[0]) && ((R) < region[0].RA[1]))6 7 3 void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap) { 8 4 9 double secz; 10 int i, j, k, n, m, N, first_j, found, found0, found1; 11 double X, Y, RADIUS, RADIUS2; 12 float *X1, *Y1, *X2, *Y2, CalM, CalC, dCalM, *Xs, *Ys; 5 int i, j, k, n, m, N, J; 6 double X, Y, RADIUS, RADIUS2, secz; 7 float *X1, *Y1, *X2, *Y2; 13 8 float dX, dY, dR; 14 9 int *N1, *N2, *next, *next_miss, last, last_miss; … … 21 16 22 17 /* photcode data */ 23 CalM = CalC = dCalM = 0;24 18 Nsecfilt = GetPhotcodeNsecfilt (); 25 19 Nsec = GetPhotcodeNsec (thiscode[0].code); 20 /* this function requires incoming stars to have the same photcode */ 26 21 27 22 /** allocate local arrays (stars) **/ … … 36 31 ALLOCATE (Y2, float, NAVE); 37 32 ALLOCATE (N2, int, NAVE); 38 ALLOCATE (Xs, float, NAVE);39 ALLOCATE (Ys, float, NAVE);40 33 ALLOCATE (catalog[0].found, int, NAVE); 41 34 REALLOCATE (catalog[0].average, Average, NAVE); … … 63 56 64 57 for (i = 0; i < Nstars; i++) { 65 fRD_to_XY (&X1[i], &Y1[i], stars[ 0][i].R, stars[0][i].D, &tcoords);58 fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); 66 59 N1[i] = i; 67 60 } … … 69 62 70 63 for (i = 0; i < Nave; i++) { 71 fRD_to_XY (&Xs[i], &Ys[i], catalog[0].average[i].R, catalog[0].average[i].D, &image[0].coords);72 64 fRD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); 73 65 N2[i] = i; … … 93 85 RADIUS = NSIGMA * 0.02 * image[0].cerror; /* 0.02 corrects cerror to arcsec from storage units */ 94 86 } else { 95 RADIUS = DEFAULT_RADIUS; 87 RADIUS = DEFAULT_RADIUS; /* provided by config */ 96 88 } 97 89 RADIUS2 = RADIUS*RADIUS; … … 108 100 109 101 dX = X1[i] - X2[j]; 110 111 102 if (dX <= -2*RADIUS) { 112 103 i++; … … 119 110 120 111 /* negative dX: j is too large, positive dX, i is too large */ 121 first_j = j; 122 for (; (dX > -2*RADIUS) && (j < Nave); j++) { 123 dX = X1[i] - X2[j]; 124 dY = Y1[i] - Y2[j]; 112 for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) { 113 dX = X1[i] - X2[J]; 114 dY = Y1[i] - Y2[J]; 125 115 dR = dX*dX + dY*dY; 126 if (dR < RADIUS2) { /* new measurement of this star */ 127 Nmatch ++; 128 n = N2[j]; 129 N = N1[i]; 130 131 /** insert star in measurement list */ 132 /* find last measurement of this star */ 133 m = catalog[0].average[n].offset; 134 for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; 135 /* set up references */ 136 next[Nmeas] = next[m]; 137 next[m] = Nmeas; 138 /* last just was moved */ 139 if (next[Nmeas] == -1) last = Nmeas; 116 if (dR < RADIUS2) continue; 117 118 Nmatch ++; 119 n = N2[J]; 120 N = N1[i]; 121 m = catalog[0].average[n].offset; 122 123 /** insert star in measurement list */ 124 /* find last measurement of this star */ 125 for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; 126 /* set up references */ 127 next[Nmeas] = next[m]; 128 next[m] = Nmeas; 129 /* last just was moved */ 130 if (next[Nmeas] == -1) last = Nmeas; 140 131 141 /* calculate accurate per-star airmass */ 142 secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude); 143 144 /* find ref measurements & add to calibration table */ 145 /* is this search right? it does not use the next[] list */ 146 if (CALIBRATE) { 147 found = FALSE; 148 found0 = FALSE; 149 found1 = FALSE; 150 m = catalog[0].average[n].offset; /* first measurement of this star */ 151 for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) { 152 if (catalog[0].measure[m+k].source == CalReference) { 153 found0 = TRUE; 154 CalM = catalog[0].measure[m+k].M; 155 dCalM = catalog[0].measure[m+k].dM; 156 } 157 if (catalog[0].measure[m+k].source == CalColor) { 158 found1 = TRUE; 159 CalC = catalog[0].measure[m+k].M; 160 } 161 if (found0 && found1) { 162 found = TRUE; 163 SaveCalibration (MIN (NO_MAG, 1000*stars[0][N].M + MTIME), 164 MIN (255, stars[0][N].dM), 165 CalM, dCalM, CalC, 166 secz, N, 167 stars[0][N].R, stars[0][N].D, Xs[n], Ys[n]); 168 } 169 } 170 } 171 172 /** add measurements for this star **/ 173 catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[0][N].R); 174 catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[0][N].D); 175 catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000*stars[0][N].M + MTIME); 176 catalog[0].measure[Nmeas].Mcal = image[0].Mcal; 177 catalog[0].measure[Nmeas].dM = MIN (255, stars[0][N].dM); /* error in input files stored in thousandths of mag */ 178 catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate; /* trate is in 0.1 msec / row */ 179 catalog[0].measure[Nmeas].averef = n; 180 catalog[0].measure[Nmeas].source = image[0].source; /* photometry source */ 181 catalog[0].measure[Nmeas].dophot = stars[0][N].dophot; 182 catalog[0].measure[Nmeas].flags = 0; 183 catalog[0].measure[Nmeas].dt = MTIME; 184 185 catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000*stars[0][N].Mgal + MTIME); 186 catalog[0].measure[Nmeas].airmass = secz; 187 catalog[0].measure[Nmeas].FWx = MIN (NO_MAG, 100*stars[0][N].fx); 188 catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx)); 189 catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[0][N].df); 190 /* refers to same number as first measurement */ 132 /* calculate accurate per-star airmass */ 133 secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude); 134 135 136 /** add measurements for this star **/ 137 catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); 138 catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); 139 catalog[0].measure[Nmeas].M = MIN (1000*stars[N].M + MTIME, NO_MAG); 140 catalog[0].measure[Nmeas].dM = MIN (1000*stars[N].dM, NO_ERR); /* error in input files stored in thousandths of mag */ 141 catalog[0].measure[Nmeas].Mcal = image[0].Mcal; 142 catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */ 143 catalog[0].measure[Nmeas].averef = n; 144 catalog[0].measure[Nmeas].source = stars[N].code; /* photcode */ 145 catalog[0].measure[Nmeas].dophot = stars[N].dophot; 146 catalog[0].measure[Nmeas].flags = 0; 147 catalog[0].measure[Nmeas].dt = MTIME; 148 catalog[0].measure[Nmeas].airmass = secz; 149 150 catalog[0].measure[Nmeas].Mgal = MIN (1000*stars[N].Mgal + MTIME, NO_MAG); 151 catalog[0].measure[Nmeas].FWx = MIN (100*stars[N].fx, NO_MAG); 152 catalog[0].measure[Nmeas].fwy = MIN (100*(stars[N].fy / stars[N].fx), NO_ERR); 153 catalog[0].measure[Nmeas].theta = MIN ((255/360)*stars[N].df, NO_ERR); 191 154 192 /* check for entries in the secfilt lists */ 193 Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K; 194 Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; 195 if (*Mval == NO_MAG) *Mval = Mrel; 196 197 /*** handle multiple stars */ 198 /* this image star matches more than one catalog star */ 199 if (stars[0][N].found > -1) { 200 catalog[0].measure[stars[0][N].found].flags |= BLEND_IMAGE; 201 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; 202 } 203 if (stars[0][N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ 204 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; 205 } 206 if (stars[0][N].found == -1) { /* this image star matches only this catalog star */ 207 stars[0][N].found = Nmeas; /* save first match, in case coincidences are found */ 208 } 209 /* this catalog star matches more than one image star */ 210 if (catalog[0].found[n] > -1) { 211 catalog[0].measure[catalog[0].found[n]].flags |= BLEND_CATALOG; 212 catalog[0].measure[Nmeas].flags |= BLEND_CATALOG; 213 } else { 214 catalog[0].found[n] = Nmeas; 215 } 216 217 catalog[0].average[n].Nm ++; 218 Nmeas ++; 219 if (Nmeas == NMEAS) { 220 NMEAS = Nmeas + 1000; 221 REALLOCATE (next, int, NMEAS); 222 REALLOCATE (catalog[0].measure, Measure, NMEAS); 223 } 224 225 /* update values of Ro, Do */ 226 if (catalog[0].average[n].Nm > 2) { 227 double R, D, r, d, r2, d2, Npt; 228 Npt = r = d = r2 = d2 = 0; 229 m = catalog[0].average[n].offset; /* first measurement of this star */ 230 for (k = 0; k < catalog[0].average[n].Nm; k++) { 231 if (catalog[0].measure[m].t == 0) { 232 m = next[m]; 233 continue; 234 } 235 R = catalog[0].measure[m].dR; 236 D = catalog[0].measure[m].dD; 237 r += R; 238 d += D; 239 r2 += R*R; 240 d2 += D*D; 241 m = next[m]; 242 Npt += 1.0; 243 } 244 if (Npt > 2) { 245 r = r / Npt; /* these are corrections in 1/100 arcsec to RA and DEC */ 246 d = d / Npt; 247 R = r2 / Npt - r*r; 248 D = d2 / Npt - d*d; 249 /* Xp is scatter in position in hundredths of arcsec */ 250 catalog[0].average[n].Xp = sqrt (D + R / SQ(cos(catalog[0].average[n].D*RAD_DEG))); 251 m = catalog[0].average[n].offset; /* first measurement of this star */ 252 for (k = 0; k < catalog[0].average[n].Nm; k++) { 253 catalog[0].measure[m].dR = catalog[0].measure[m].dR - r; 254 catalog[0].measure[m].dD = catalog[0].measure[m].dD - d; 255 m = next[m]; 256 } 257 catalog[0].average[n].R = catalog[0].average[n].R - r / 360000.0; 258 catalog[0].average[n].D = catalog[0].average[n].D - d / 360000.0; 259 } 260 } 261 } 262 } 263 j = first_j; 155 /* it is not valid to pass PRI/SEC/REF photcodes to this routine */ 156 /* check for entries in the secfilt lists */ 157 Mrel = catalog[0].measure[Nmeas].M + secz * thiscode[0].K; 158 Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; 159 if (*Mval == NO_MAG) *Mval = Mrel; 160 161 if (CALIBRATE) { 162 /** this is a little tricky: the new measures are not in the measure sequence 163 this call can run into the next star if we have had two measures found at this location. **/ 164 AddToCalibration (&catalog[0].average[n], &catalog[0].measure[m], &catalog[0].measure[Nmeas], N); 165 } 166 167 /*** flag multiple stars */ 168 /* this image star matches more than one catalog star */ 169 if (stars[N].found > -1) { 170 catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE; 171 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; 172 } 173 if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ 174 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; 175 } 176 if (stars[N].found == -1) { /* this image star matches only this catalog star */ 177 stars[N].found = Nmeas; /* save first match, in case coincidences are found */ 178 } 179 /* this catalog star matches more than one image star */ 180 if (catalog[0].found[n] > -1) { 181 catalog[0].measure[catalog[0].found[n]].flags |= BLEND_CATALOG; 182 catalog[0].measure[Nmeas].flags |= BLEND_CATALOG; 183 } else { 184 catalog[0].found[n] = Nmeas; 185 } 186 187 catalog[0].average[n].Nm ++; 188 Nmeas ++; 189 if (Nmeas == NMEAS) { 190 NMEAS = Nmeas + 1000; 191 REALLOCATE (next, int, NMEAS); 192 REALLOCATE (catalog[0].measure, Measure, NMEAS); 193 } 194 195 update_coords (&catalog[0].average[n], &catalog[0].measure[0], next); 196 } 264 197 i++; 265 198 } … … 270 203 if (catalog[0].found[n] < 0) { 271 204 /* should the catalog star be on this image? project into image coords */ 272 if (in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) { 273 /* find last missing meas */ 274 if (catalog[0].average[n].Nn < 1) { /* no previous missing obs */ 275 catalog[0].average[n].missing = Nmiss; 276 next_miss[last_miss] = Nmiss; 277 next_miss[Nmiss] = -1; 205 if (!in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) continue; 206 /* find last missing meas */ 207 if (catalog[0].average[n].Nn < 1) { /* no previous missing obs */ 208 catalog[0].average[n].missing = Nmiss; 209 next_miss[last_miss] = Nmiss; 210 next_miss[Nmiss] = -1; 211 last_miss = Nmiss; 212 } else { 213 m = catalog[0].average[n].missing; 214 for (k = 0; k < catalog[0].average[n].Nn - 1; k++) 215 m = next_miss[m]; 216 next_miss[Nmiss] = next_miss[m]; 217 next_miss[m] = Nmiss; 218 if (next_miss[Nmiss] == -1) { /* last just was moved */ 278 219 last_miss = Nmiss; 279 } else {280 m = catalog[0].average[n].missing;281 for (k = 0; k < catalog[0].average[n].Nn - 1; k++)282 m = next_miss[m];283 next_miss[Nmiss] = next_miss[m];284 next_miss[m] = Nmiss;285 if (next_miss[Nmiss] == -1) { /* last just was moved */286 last_miss = Nmiss;287 }288 220 } 289 RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords); 290 catalog[0].missing[Nmiss].t = image[0].tzero + 1e-4*Y*image[0].trate; /* trate is in 0.1 msec / row */ 291 catalog[0].average[n].Nn ++; 292 Nmiss ++; 293 if (Nmiss == NMISS) { 294 NMISS = Nmiss + 1000; 295 REALLOCATE (next_miss, int, NMISS); 296 REALLOCATE (catalog[0].missing, Missing, NMISS); 297 } 221 } 222 /* calculate time of exposure for this coordinate in the image */ 223 RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords); 224 catalog[0].missing[Nmiss].t = image[0].tzero + 1e-4*Y*image[0].trate; /* trate is in 0.1 msec / row */ 225 catalog[0].average[n].Nn ++; 226 Nmiss ++; 227 if (Nmiss == NMISS) { 228 NMISS = Nmiss + 1000; 229 REALLOCATE (next_miss, int, NMISS); 230 REALLOCATE (catalog[0].missing, Missing, NMISS); 298 231 } 299 232 } … … 301 234 302 235 /* incorporate unmatched image stars, if this star is in field of this catalog */ 303 for (i = 0; i < Nstars; i++) {236 for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { 304 237 N = N1[i]; 305 if ((stars[0][N].found < 0) && IN_CATALOG (stars[0][N].R, stars[0][N].D)) { 306 catalog[0].average[Nave].R = stars[0][N].R; 307 catalog[0].average[Nave].D = stars[0][N].D; 308 catalog[0].average[Nave].M = NO_MAG; 309 for (j = 0; j < Nsecfilt; j++) { 310 catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; 311 catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; 312 } 313 314 if (SKYPROBE) 315 secz = 1000 * airmass (stars[0][N].R, stars[0][N].D, SiderealTime, Latitude); 316 else 317 secz = image[0].secz; 318 319 catalog[0].average[Nave].Nm = 1; 320 catalog[0].average[Nave].Nn = 0; 321 catalog[0].average[Nave].Xp = NO_MAG; 322 catalog[0].average[Nave].Xm = NO_MAG; 323 catalog[0].average[Nave].dM = NO_MAG; 324 catalog[0].average[Nave].offset = Nmeas; 325 catalog[0].average[Nave].missing = -1; 326 catalog[0].average[Nave].code = 0; 327 328 catalog[0].measure[Nmeas].dR = 0.0; 329 catalog[0].measure[Nmeas].dD = 0.0; 330 catalog[0].measure[Nmeas].M = MIN (NO_MAG, 1000.0*stars[0][N].M + MTIME); 331 catalog[0].measure[Nmeas].Mcal = image[0].Mcal; 332 catalog[0].measure[Nmeas].dM = MIN (255, stars[0][N].dM); 333 catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[0][N].Y*image[0].trate; /* trate is in 0.1 msec / row */ 334 catalog[0].measure[Nmeas].averef = Nave; 335 catalog[0].measure[Nmeas].source = image[0].source; /* photometry source */ 336 catalog[0].measure[Nmeas].dophot = stars[0][N].dophot; 337 catalog[0].measure[Nmeas].flags = 0; 338 catalog[0].measure[Nmeas].dt = MTIME; 339 340 catalog[0].measure[Nmeas].Mgal = MIN (NO_MAG, 1000.0*stars[0][N].Mgal + MTIME); 341 catalog[0].measure[Nmeas].airmass = secz; 342 catalog[0].measure[Nmeas].FWx = MIN (0x7fff, 100*stars[0][N].fx); 343 catalog[0].measure[Nmeas].fwy = MIN (255, 100*(stars[0][N].fy / stars[0][N].fx)); 344 catalog[0].measure[Nmeas].theta = MIN (255, (255/360)*stars[0][N].df); 345 346 Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K; 347 Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M; 348 if (*Mval == NO_MAG) *Mval = Mrel; 349 350 stars[0][N].found = Nmeas; 351 next[last] = Nmeas; 352 next[Nmeas] = -1; 353 last = Nmeas; 354 Nmeas ++; 355 if (Nmeas == NMEAS) { 356 NMEAS = Nmeas + 1000; 357 REALLOCATE (next, int, NMEAS); 358 REALLOCATE (catalog[0].measure, Measure, NMEAS); 359 } 360 361 /** now add references from all previous non-detection observations of this spot on the sky */ 362 for (j = 0; (j < Noverlap) && !SKIP_MISSED; j++) { 363 if (in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) { 364 if (catalog[0].average[Nave].Nn < 1) { 365 catalog[0].average[Nave].missing = Nmiss; 366 } 367 next_miss[last_miss] = Nmiss; 368 next_miss[Nmiss] = -1; 369 last_miss = Nmiss; 370 /* this can now be done exactly */ 371 RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j].coords); 372 catalog[0].missing[Nmiss].t = overlap[j].tzero + 1e-4*Y*overlap[j].trate; /* rough guess at time */ 373 catalog[0].average[Nave].Nn ++; 374 Nmiss ++; 375 if (Nmiss == NMISS) { 376 NMISS = Nmiss + 1000; 377 REALLOCATE (next_miss, int, NMISS); 378 REALLOCATE (catalog[0].missing, Missing, NMISS); 379 } 380 381 } 382 } 383 Nave ++; 384 if (Nave == NAVE) { 385 NAVE = Nave + 1000; 386 REALLOCATE (catalog[0].average, Average, NAVE); 387 REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); 388 } 238 if (stars[N].found >= 0) continue; 239 if (!IN_CATALOG (stars[N].R, stars[N].D)) continue; 240 241 secz = airmass (image[0].secz, stars[N].R, stars[N].D, SiderealTime, Latitude); 242 243 catalog[0].average[Nave].R = stars[N].R; 244 catalog[0].average[Nave].D = stars[N].D; 245 catalog[0].average[Nave].M = NO_MAG; 246 catalog[0].average[Nave].dM = NO_MAG; 247 catalog[0].average[Nave].Nm = 1; 248 catalog[0].average[Nave].Nn = 0; 249 catalog[0].average[Nave].Xp = NO_MAG; 250 catalog[0].average[Nave].Xm = NO_MAG; 251 catalog[0].average[Nave].Xg = NO_MAG; 252 catalog[0].average[Nave].offset = Nmeas; 253 catalog[0].average[Nave].missing = -1; 254 catalog[0].average[Nave].code = 0; 255 256 for (j = 0; j < Nsecfilt; j++) { 257 catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; 258 catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; 259 catalog[0].secfilt[Nave*Nsecfilt+j].dM = NO_MAG; 260 } 261 262 catalog[0].measure[Nmeas].dR = 0.0; 263 catalog[0].measure[Nmeas].dD = 0.0; 264 catalog[0].measure[Nmeas].M = MIN (1000*stars[N].M + MTIME, NO_MAG); 265 catalog[0].measure[Nmeas].dM = MIN (1000*stars[N].dM, NO_ERR); 266 catalog[0].measure[Nmeas].Mcal = image[0].Mcal; 267 catalog[0].measure[Nmeas].t = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */ 268 catalog[0].measure[Nmeas].averef = Nave; 269 catalog[0].measure[Nmeas].source = stars[N].code; /* photcode */ 270 catalog[0].measure[Nmeas].dophot = stars[N].dophot; 271 catalog[0].measure[Nmeas].flags = 0; 272 catalog[0].measure[Nmeas].dt = MTIME; 273 catalog[0].measure[Nmeas].airmass = secz; 274 275 catalog[0].measure[Nmeas].Mgal = MIN (1000*stars[N].Mgal + MTIME, NO_MAG); 276 catalog[0].measure[Nmeas].FWx = MIN (100*stars[N].fx, NO_MAG); 277 catalog[0].measure[Nmeas].fwy = MIN (100*(stars[N].fy / stars[N].fx), NO_ERR); 278 catalog[0].measure[Nmeas].theta = MIN ((255/360)*stars[N].df, NO_ERR); 279 280 Mrel = catalog[0].measure[Nmeas].M + secz * thiscode[0].K; 281 Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M; 282 if (*Mval == NO_MAG) *Mval = Mrel; 283 284 /** now add references from all previous non-detection observations of this spot on the sky */ 285 for (j = 0; (j < Noverlap) && !SKIP_MISSED; j++) { 286 if (!in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j])) continue; 287 if (catalog[0].average[Nave].Nn < 1) { 288 catalog[0].average[Nave].missing = Nmiss; 289 } 290 next_miss[last_miss] = Nmiss; 291 next_miss[Nmiss] = -1; 292 last_miss = Nmiss; 293 /* get time of exposure of this portion of the image */ 294 RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &overlap[j].coords); 295 catalog[0].missing[Nmiss].t = overlap[j].tzero + 1e-4*Y*overlap[j].trate; /* rough guess at time */ 296 catalog[0].average[Nave].Nn ++; 297 Nmiss ++; 298 if (Nmiss == NMISS) { 299 NMISS = Nmiss + 1000; 300 REALLOCATE (next_miss, int, NMISS); 301 REALLOCATE (catalog[0].missing, Missing, NMISS); 302 } 303 } 304 305 stars[N].found = Nmeas; 306 next[last] = Nmeas; 307 next[Nmeas] = -1; 308 last = Nmeas; 309 Nmeas ++; 310 if (Nmeas == NMEAS) { 311 NMEAS = Nmeas + 1000; 312 REALLOCATE (next, int, NMEAS); 313 REALLOCATE (catalog[0].measure, Measure, NMEAS); 314 } 315 Nave ++; 316 if (Nave == NAVE) { 317 NAVE = Nave + 1000; 318 REALLOCATE (catalog[0].average, Average, NAVE); 319 REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); 389 320 } 390 321 } … … 393 324 REALLOCATE (catalog[0].average, Average, Nave); 394 325 REALLOCATE (catalog[0].measure, Measure, Nmeas); 395 326 396 327 /* fix order of Measure (memory intensive, but fast) */ 397 328 N = 0; … … 427 358 /* note stars which have been found in this catalog */ 428 359 for (i = 0; i < Nstars; i++) { 429 if (stars[0][i].found > -1) { 430 stars[0][i].found = -2; 360 if (stars[i].found > -1) { 361 stars[i].found = -2; 362 } else { 363 stars[i].found = -3; 431 364 } 432 365 } … … 438 371 } 439 372 440 int in_image (double r, double d, Image *image) {441 442 double X, Y;443 444 RD_to_XY (&X, &Y, r, d, &image[0].coords);445 if (((X) >= 0) && ((X) < image[0].NX) && ((Y) >= 0) && ((Y) < image[0].NY))446 return (TRUE);447 else448 return (FALSE);449 450 }451 452 373 /* 453 374 notes: 454 455 for finding if a catalog star is in an image or an image star is in the catalog: 456 457 catalogs have boundaries defined by RA and DEC, but they may curve in projection 458 images have boundaries which are lines in pixels coords, but curve in RA and DEC 459 460 catalog[0].found[Ncat] but stars[Nstar].found 461 462 */ 463 464 465 # ifdef OLDJUNKHERE 466 467 468 /* fix measurement list order */ 469 for (i = 0; i < Nmeas - 1; i++) { 470 if (next[i] != i + 1) { 471 n0 = next[i]; 472 n1 = next[i+1]; 473 n2 = next[n0]; 474 t = i; 475 while ((t != -1) && (next[t] != i+1)) 476 t = next[t]; 477 if ((t == -1) || (n2 == i+1)) { 478 t = n2; 479 } 480 next[i] = i+1; 481 next[i+1] = n2; 482 next[n0] = n1; 483 next[t] = n0; 484 tmpmeasure = catalog[0].measure[n0]; 485 catalog[0].measure[n0] = catalog[0].measure[i+1]; 486 catalog[0].measure[i+1] = tmpmeasure; 487 } 488 } 489 490 /* fix offset references */ 491 n = (0x00ffffff & catalog[0].measure[0].average); 492 catalog[0].average[n].offset = 0; 493 for (i = 1; i < Nmeas; i++) { 494 if (n != (0x00ffffff & catalog[0].measure[i].average)) { 495 n = (0x00ffffff & catalog[0].measure[i].average); 496 catalog[0].average[n].offset = i; 497 } 498 } 499 500 /* fix missing list order */ 501 for (i = 0; i < Nmiss - 1; i++) { 502 if (next_miss[i] != i + 1) { 503 n0 = next_miss[i]; 504 n1 = next_miss[i+1]; 505 n2 = next_miss[n0]; 506 t = i; 507 while ((t != -1) && (next_miss[t] != i+1)) 508 t = next_miss[t]; 509 if ((t == -1) || (n2 == i+1)) { 510 t = n2; 511 } 512 next_miss[i] = i+1; 513 next_miss[i+1] = n2; 514 next_miss[n0] = n1; 515 next_miss[t] = n0; 516 tmpmissing = catalog[0].missing[n0]; 517 catalog[0].missing[n0] = catalog[0].missing[i+1]; 518 catalog[0].missing[i+1] = tmpmissing; 519 } 520 } 521 522 /* fix offset references */ 523 n = (0x00ffffff & catalog[0].measure[0].average); 524 catalog[0].average[n].offset = 0; 525 for (i = 1; i < Nmeas; i++) { 526 if (n != (0x00ffffff & catalog[0].measure[i].average)) { 527 n = (0x00ffffff & catalog[0].measure[i].average); 528 catalog[0].average[n].offset = i; 529 } 530 } 531 532 # endif 375 376 for finding if a catalog star is in an image or an image star is in the catalog: 377 378 catalogs have boundaries defined by RA and DEC, but they may curve in projection 379 images have boundaries which are lines in pixels coords, but curve in RA and DEC 380 381 catalog[0].found[Ncat] but stars[Nstar].found 382 383 */ -
trunk/Ohana/src/addstar/src/find_matches_refstars.c
r3347 r3361 1 1 # include "addstar.h" 2 2 3 # define IN_CATALOG(R,D) ( \4 ((D) >= region[0].DEC[0]) && ((D) < region[0].DEC[1]) && \5 ((R) >= region[0].RA[0]) && ((R) < region[0].RA[1]))6 7 3 void find_matches_refstars (GSCRegion *region, Stars **stars, int Nstars, Catalog *catalog) { 8 4 9 int i, j, k, n, m, N, first_j;5 int i, j, k, n, m, N, J; 10 6 double RADIUS, RADIUS2; 11 7 float *X1, *Y1, *X2, *Y2; … … 16 12 Missing *tmpmissing; 17 13 Coords tcoords; 18 int Nsecfilt , Nsec;14 int Nsecfilt; 19 15 20 16 /* photcode data -- should not have to modify secfilt / average */ 21 17 Nsecfilt = GetPhotcodeNsecfilt (); 22 Nsec = GetPhotcodeNsec (thiscode[0].code);23 18 24 19 /** allocate local arrays (stars) **/ … … 88 83 /* choose a radius for matches */ 89 84 if (DEFAULT_RADIUS == 0) { 90 RADIUS = 2.0; 85 RADIUS = 2.0; /* hardwired default for refstars */ 91 86 } else { 92 RADIUS = DEFAULT_RADIUS; 87 RADIUS = DEFAULT_RADIUS; /* provided by config */ 93 88 } 94 89 RADIUS2 = RADIUS*RADIUS; … … 98 93 99 94 dX = X1[i] - X2[j]; 100 101 95 if (dX <= -2*RADIUS) { 102 96 i++; … … 109 103 110 104 /* negative dX: j is too large, positive dX, i is too large */ 111 first_j = j; 112 for (; (dX > -2*RADIUS) && (j < Nave); j++) { 113 dX = X1[i] - X2[j]; 114 dY = Y1[i] - Y2[j]; 105 for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) { 106 dX = X1[i] - X2[J]; 107 dY = Y1[i] - Y2[J]; 115 108 dR = dX*dX + dY*dY; 116 if (dR < RADIUS2) { /* new measurement of this star */ 117 Nmatch ++; 118 n = N2[j]; 119 N = N1[i]; 120 121 /** in replace mode, search for entry and replace values M, dM, R, D */ 122 if (REPLACE) { 123 int found; 124 found = FALSE; 125 m = catalog[0].average[n].offset; 126 for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) { 127 if (catalog[0].measure[m+k].source == thiscode[0].code) { 128 catalog[0].measure[m+k].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); 129 catalog[0].measure[m+k].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); 130 catalog[0].measure[m+k].M = 1000.0*stars[N].M; 131 catalog[0].measure[m+k].dM = 1000.0*MIN (255, stars[N].dM); /* error in input files stored in thousandths of mag */ 132 found = TRUE; 133 stars[N].found = m+k; /* save first match, in case coincidences are found */ 134 } 135 } 136 if (found) continue; 137 } 138 139 /** insert star in measurement list */ 140 /* find last measurement of this star */ 141 m = catalog[0].average[n].offset; 142 for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; 143 /* set up references */ 144 next[Nmeas] = next[m]; 145 next[m] = Nmeas; 146 /* last just was moved */ 147 if (next[Nmeas] == -1) last = Nmeas; 109 if (dR > RADIUS2) continue; 110 111 Nmatch ++; 112 n = N2[J]; 113 N = N1[i]; 114 m = catalog[0].average[n].offset; 115 116 /** in replace mode, search for entry and replace values M, dM, R, D */ 117 if (REPLACE && replace_match (&catalog[0].average[n], &catalog[0].measure[m], stars[N])) continue; 118 119 /** insert star in measurement list */ 120 /* find last measurement of this star */ 121 for (k = 0; k < catalog[0].average[n].Nm - 1; k++) m = next[m]; 122 /* set up references */ 123 next[Nmeas] = next[m]; 124 next[m] = Nmeas; 125 /* last just was moved */ 126 if (next[Nmeas] == -1) last = Nmeas; 148 127 149 /** add measurements for this star **/ 150 catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N].R); 151 catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N].D); 152 catalog[0].measure[Nmeas].M = 1000.0*stars[N].M; 153 catalog[0].measure[Nmeas].Mcal = 0; 154 catalog[0].measure[Nmeas].dM = 1000.0*MIN (255, stars[N].dM); /* error in input files stored in thousandths of mag */ 155 catalog[0].measure[Nmeas].t = TIMEREF; /** careful : time_t vs e_time **/ 156 catalog[0].measure[Nmeas].averef = n; 157 catalog[0].measure[Nmeas].source = thiscode[0].code; 158 catalog[0].measure[Nmeas].dophot = stars[N].dophot; /* error in input files stored in thousandths of mag */ 159 catalog[0].measure[Nmeas].flags = 0; 160 catalog[0].measure[Nmeas].dt = 0xffff; 161 162 catalog[0].measure[Nmeas].Mgal = NO_MAG; 163 catalog[0].measure[Nmeas].airmass = 0; 164 catalog[0].measure[Nmeas].FWx = NO_MAG; 165 catalog[0].measure[Nmeas].fwy = 0xff; 166 catalog[0].measure[Nmeas].theta = 0xff; 167 168 /* check for entries in the secfilt lists */ 169 Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M; 170 if (*Mval == NO_MAG) *Mval = catalog[0].measure[Nmeas].M; /*** no airmass correction needed ***/ 171 /*** the problem here is that the ref photcode should never equate to a pri/sec photcode, right? ***/ 172 /* just skip assigning an average magnitude for this object? */ 173 174 /*** handle multiple stars */ 175 /* this image star matches more than one catalog star */ 176 if (stars[N].found > -1) { 177 catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE; 178 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; 179 } 180 if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ 181 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; 182 } 183 if (stars[N].found == -1) { /* this image star matches only this catalog star */ 184 stars[N].found = Nmeas; /* save first match, in case coincidences are found */ 185 } 186 /* this catalog star matches more than one image star */ 187 if (catalog[0].found[n] > -1) { 188 catalog[0].measure[catalog[0].found[n]].flags |= BLEND_CATALOG; 189 catalog[0].measure[Nmeas].flags |= BLEND_CATALOG; 190 } else { 191 catalog[0].found[n] = Nmeas; 192 } 193 194 catalog[0].average[n].Nm ++; 195 Nmeas ++; 196 if (Nmeas == NMEAS) { 197 NMEAS = Nmeas + 1000; 198 REALLOCATE (next, int, NMEAS); 199 REALLOCATE (catalog[0].measure, Measure, NMEAS); 200 } 201 202 /* update values of Ro, Do */ 203 if (catalog[0].average[n].Nm > 2) { 204 double R, D, r, d, r2, d2, Npt; 205 Npt = r = d = r2 = d2 = 0; 206 m = catalog[0].average[n].offset; /* first measurement of this star */ 207 for (k = 0; k < catalog[0].average[n].Nm; k++) { 208 if (catalog[0].measure[m].t == 0) { 209 m = next[m]; 210 continue; 211 } 212 R = catalog[0].measure[m].dR; 213 D = catalog[0].measure[m].dD; 214 r += R; 215 d += D; 216 r2 += R*R; 217 d2 += D*D; 218 m = next[m]; 219 Npt += 1.0; 220 } 221 if (Npt > 2) { 222 r = r / Npt; /* these are corrections in 1/100 arcsec to RA and DEC */ 223 d = d / Npt; 224 R = r2 / Npt - r*r; 225 D = d2 / Npt - d*d; 226 /* Xp is scatter in position in hundredths of arcsec */ 227 catalog[0].average[n].Xp = sqrt (D + R / SQ(cos(catalog[0].average[n].D*RAD_DEG))); 228 m = catalog[0].average[n].offset; /* first measurement of this star */ 229 for (k = 0; k < catalog[0].average[n].Nm; k++) { 230 catalog[0].measure[m].dR = catalog[0].measure[m].dR - r; 231 catalog[0].measure[m].dD = catalog[0].measure[m].dD - d; 232 m = next[m]; 233 } 234 catalog[0].average[n].R = catalog[0].average[n].R - r / 360000.0; 235 catalog[0].average[n].D = catalog[0].average[n].D - d / 360000.0; 236 } 237 } 238 } 239 } 240 j = first_j; 241 i++; 242 } 243 244 /* incorporate unmatched image stars, if this star is in field of this catalog */ 245 /* skip if we want to force matches - combined with -replace, this lets us keep 246 the reference up-to-date with known stars only */ 247 for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { 248 N = N1[i]; 249 if (stars[N].found < 0) { 250 catalog[0].average[Nave].R = stars[N].R; 251 catalog[0].average[Nave].D = stars[N].D; 252 catalog[0].average[Nave].M = NO_MAG; 253 for (j = 0; j < Nsecfilt; j++) { 254 catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; 255 catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; 256 } 257 258 if (thiscode[0].type == PHOT_PRI) { 259 catalog[0].average[Nave].M = 1000.0*stars[N].M; 260 } 261 if (thiscode[0].type == PHOT_SEC) { 262 catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = 1000.0*stars[N].M; 263 } 264 catalog[0].average[Nave].Nm = 1; 265 catalog[0].average[Nave].Nn = 0; 266 catalog[0].average[Nave].Xp = NO_MAG; 267 catalog[0].average[Nave].Xm = NO_MAG; 268 catalog[0].average[Nave].dM = NO_MAG; 269 catalog[0].average[Nave].offset = Nmeas; 270 catalog[0].average[Nave].missing = -1; 271 catalog[0].average[Nave].code = 0; 272 273 catalog[0].measure[Nmeas].dR = 0.0; 274 catalog[0].measure[Nmeas].dD = 0.0; 275 catalog[0].measure[Nmeas].M = 1000.0*stars[N].M; 128 /** add measurements for this star **/ 129 catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - stars[N][0].R); 130 catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - stars[N][0].D); 131 catalog[0].measure[Nmeas].M = MIN (1000*stars[N][0].M, NO_MAG); 132 catalog[0].measure[Nmeas].dM = MIN (1000*stars[N][0].dM, NO_ERR); 276 133 catalog[0].measure[Nmeas].Mcal = 0; 277 catalog[0].measure[Nmeas].dM = MIN (255, stars[N].dM); 278 catalog[0].measure[Nmeas].t = 0; 279 catalog[0].measure[Nmeas].averef = Nave; 280 catalog[0].measure[Nmeas].source = thiscode[0].code; 281 catalog[0].measure[Nmeas].dophot = stars[N].dophot; /* error in input files stored in thousandths of mag */ 134 catalog[0].measure[Nmeas].t = (TIMEREF == 0) ? stars[N][0].t : TIMEREF; /** careful : time_t vs e_time **/ 135 catalog[0].measure[Nmeas].averef = n; 136 catalog[0].measure[Nmeas].source = stars[N][0].code; 137 catalog[0].measure[Nmeas].dophot = 0; 282 138 catalog[0].measure[Nmeas].flags = 0; 283 139 catalog[0].measure[Nmeas].dt = 0xffff; … … 288 144 catalog[0].measure[Nmeas].fwy = 0xff; 289 145 catalog[0].measure[Nmeas].theta = 0xff; 290 291 stars[N].found = Nmeas; 292 next[last] = Nmeas; 293 next[Nmeas] = -1; 294 last = Nmeas; 146 147 /** don't update average / secfilt values for REF photcodes **/ 148 149 /*** handle multiple stars */ 150 /* this image star matches more than one catalog star */ 151 if (stars[N][0].found > -1) { 152 catalog[0].measure[stars[N][0].found].flags |= BLEND_IMAGE; 153 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE; 154 } 155 if (stars[N][0].found == -2) { /* this image star matches a catalog star on a neighboring catalog */ 156 catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR; 157 } 158 if (stars[N][0].found == -1) { /* this image star matches only this catalog star */ 159 stars[N][0].found = Nmeas; /* save first match, in case coincidences are found */ 160 } 161 /* this catalog star matches more than one image star */ 162 if (catalog[0].found[n] > -1) { 163 catalog[0].measure[catalog[0].found[n]].flags |= BLEND_CATALOG; 164 catalog[0].measure[Nmeas].flags |= BLEND_CATALOG; 165 } else { 166 catalog[0].found[n] = Nmeas; 167 } 168 169 catalog[0].average[n].Nm ++; 295 170 Nmeas ++; 296 171 if (Nmeas == NMEAS) { … … 299 174 REALLOCATE (catalog[0].measure, Measure, NMEAS); 300 175 } 301 Nave ++; 302 if (Nave == NAVE) { 303 NAVE = Nave + 1000; 304 REALLOCATE (catalog[0].average, Average, NAVE); 305 REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); 306 } 176 177 update_coords (&catalog[0].average[n], &catalog[0].measure[0], next); 178 } 179 i++; 180 } 181 182 /* we don't add missed entries for refcat 183 (already in database, not refcat) */ 184 185 /* incorporate unmatched refcat stars */ 186 /* skip if we want to require matches 187 combined with -replace, this lets us keep 188 the reference up-to-date with known stars only */ 189 190 for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { 191 N = N1[i]; 192 if (stars[N][0].found >= 0) continue; 193 194 catalog[0].average[Nave].R = stars[N][0].R; 195 catalog[0].average[Nave].D = stars[N][0].D; 196 catalog[0].average[Nave].M = NO_MAG; 197 catalog[0].average[Nave].dM = NO_MAG; 198 catalog[0].average[Nave].Nm = 1; 199 catalog[0].average[Nave].Nn = 0; 200 catalog[0].average[Nave].Xp = NO_MAG; 201 catalog[0].average[Nave].Xm = NO_MAG; 202 catalog[0].average[Nave].Xg = NO_MAG; 203 catalog[0].average[Nave].offset = Nmeas; 204 catalog[0].average[Nave].missing = -1; 205 catalog[0].average[Nave].code = 0; 206 207 for (j = 0; j < Nsecfilt; j++) { 208 catalog[0].secfilt[Nave*Nsecfilt+j].M = NO_MAG; 209 catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG; 210 catalog[0].secfilt[Nave*Nsecfilt+j].dM = NO_MAG; 211 } 212 213 catalog[0].measure[Nmeas].dR = 0.0; 214 catalog[0].measure[Nmeas].dD = 0.0; 215 catalog[0].measure[Nmeas].M = MIN (1000*stars[N][0].M, NO_MAG); 216 catalog[0].measure[Nmeas].dM = MIN (1000*stars[N][0].dM, NO_ERR); 217 catalog[0].measure[Nmeas].Mcal = 0; 218 catalog[0].measure[Nmeas].t = (stars[N][0].t == 0) ? TIMEREF : stars[N][0].t; /** careful : time_t vs e_time **/ 219 catalog[0].measure[Nmeas].averef = Nave; 220 catalog[0].measure[Nmeas].source = stars[N][0].code; 221 catalog[0].measure[Nmeas].dophot = 0; 222 catalog[0].measure[Nmeas].flags = 0; 223 catalog[0].measure[Nmeas].dt = 0xffff; 224 225 catalog[0].measure[Nmeas].Mgal = NO_MAG; 226 catalog[0].measure[Nmeas].airmass = 0; 227 catalog[0].measure[Nmeas].FWx = NO_MAG; 228 catalog[0].measure[Nmeas].fwy = 0xff; 229 catalog[0].measure[Nmeas].theta = 0xff; 230 231 stars[N][0].found = Nmeas; 232 next[last] = Nmeas; 233 next[Nmeas] = -1; 234 last = Nmeas; 235 Nmeas ++; 236 if (Nmeas == NMEAS) { 237 NMEAS = Nmeas + 1000; 238 REALLOCATE (next, int, NMEAS); 239 REALLOCATE (catalog[0].measure, Measure, NMEAS); 240 } 241 Nave ++; 242 if (Nave == NAVE) { 243 NAVE = Nave + 1000; 244 REALLOCATE (catalog[0].average, Average, NAVE); 245 REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); 307 246 } 308 247 } … … 345 284 /* note stars which have been found in this catalog */ 346 285 for (i = 0; i < Nstars; i++) { 347 if (stars[i] .found > -1) {348 stars[i] .found = -2;286 if (stars[i][0].found > -1) { 287 stars[i][0].found = -2; 349 288 } else { 350 stars[i] .found = -3;289 stars[i][0].found = -3; 351 290 } 352 291 } … … 357 296 if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch); 358 297 } 359 360 361 # if (0)362 /* this is a bit of an issue. i don't want to include missing obs if the363 ref star has not been previously detected (why?), but if i do this, it won't364 get all missed entries if ref star is detected in the future... */365 /** now add references from all previous non-detection observations of this spot on the sky */366 for (j = 0; j < Nimage; j++) {367 if (in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j])) {368 if (catalog[0].average[Nave].Nn < 1) {369 catalog[0].average[Nave].missing = Nmiss;370 }371 next_miss[last_miss] = Nmiss;372 next_miss[Nmiss] = -1;373 last_miss = Nmiss;374 /* this can now be done exactly */375 RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j].coords);376 catalog[0].missing[Nmiss].t = pimage[j].tzero + 1e-4*Y*pimage[j].trate; /* rough guess at time */377 catalog[0].average[Nave].Nn ++;378 Nmiss ++;379 if (Nmiss == NMISS) {380 NMISS = Nmiss + 1000;381 REALLOCATE (next_miss, int, NMISS);382 REALLOCATE (catalog[0].missing, Missing, NMISS);383 }384 385 }386 }387 # endif -
trunk/Ohana/src/addstar/src/find_subset.c
r3339 r3361 1 1 # include "addstar.h" 2 2 3 Stars *find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { 3 /* find stars within this region */ 4 Stars **find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { 4 5 5 6 int i, N, NSUBSET; … … 28 29 if (N == NSUBSET - 1) { 29 30 NSUBSET += 1000; 30 REALLOCATE (subset, Stars , NSUBSET);31 REALLOCATE (subset, Stars *, NSUBSET); 31 32 } 32 33 } -
trunk/Ohana/src/addstar/src/gcatalog.c
r2457 r3361 25 25 for (in = out = i = 0; i < catalog[0].Naverage; i++) { 26 26 for (j = 0; j < catalog[0].Nsecfilt; j++, in++, out++) { 27 outsec[out].M = insec[in].M;27 outsec[out].M = insec[in].M; 28 28 outsec[out].Xm = insec[in].Xm; 29 outsec[out].dM = insec[in].dM; 29 30 } 30 31 for (j = 0; j < Nextra; j++, out++) { 31 outsec[out].M = NO_MAG;32 outsec[out].M = NO_MAG; 32 33 outsec[out].Xm = NO_MAG; 34 outsec[out].dM = NO_MAG; 33 35 } 34 36 } -
trunk/Ohana/src/addstar/src/get2mass.c
r3347 r3361 1 # include " mosastro.h"1 # include "addstar.h" 2 2 3 Star Data *get2mass (CatStats *catstats, int *NSTARS) {3 Stars *get2mass (GSCRegion *patch, int *NSTARS) { 4 4 5 int i, j, k, Ns, Ngsc, Nregions, Nstars; 6 StarData *gsc; 7 StarData *stars; 8 CatStats *regions; 5 int i, j, k, Ns, Nstars, Nrefcat, Nregions; 6 Stars *stars; 7 Stars *refcat; 8 TM_Region *regions; 9 short int TM_J, TM_H, TM_K; 9 10 11 NAMED_PHOTCODE (TM_J, "2MASS_J"); 12 NAMED_PHOTCODE (TM_H, "2MASS_H"); 13 NAMED_PHOTCODE (TM_K, "2MASS_K"); 14 15 regions = get2mass_acc (patch, &Nregions); 16 10 17 Nstars = 0; 11 ALLOCATE (stars, StarData, 1); 12 13 set_catalog (TWO_MASS_DIR); 14 regions = gregions (catstats, &Nregions); 18 ALLOCATE (stars, Stars, 1); 15 19 16 20 for (i = 0; i < Nregions; i++) { 17 gsc = gcatalog (regions[i].name, &Ngsc); 18 if (gsc == NULL) continue; 21 refcat = get2mass_data (regions[i].name, patch, &Nrefcat); 19 22 20 Ns = Nstars; 21 Nstars += Ngsc; 23 REALLOCATE (stars, Stars, MAX (1, Nstars + Nrefcat)); 24 memcpy (&stars[Nstars], refcat, Nrefcat*sizeof(Stars)); 25 Nstars += Nrefcat; 22 26 23 REALLOCATE (stars, StarData, Nstars); 24 for (k = Ns, j = 0; j < Ngsc; k++, j++) { 25 stars[k].R = gsc[j].R; 26 stars[k].D = gsc[j].D; 27 stars[k].M = gsc[j].M; 28 } 29 free (gsc); 27 free (refcat); 30 28 } 31 29 … … 34 32 return (stars); 35 33 } 34 35 /* watch for patch which cross 0,360 boundary */ 36 TM_Region *get2mass_acc (GSCRegion *patch, int *nfiles) { 37 38 int Nregions 39 TM_Region *regions; 40 41 sprintf (accel, "%s/idr2.acc", TWO_MASS_DIR); 42 f = fopen (accel, "r"); 43 if (f == NULL) { 44 fprintf (stderr, "ERROR: can't read data from %s\n", file); 45 Shutdown (); 46 } 47 48 Nregions = 0; 49 NREGIONS = 20; 50 ALLOCATE (regions, TM_Region, NREGIONS); 51 52 /* read in stars line-by-line */ 53 for (N = 0; scan_line (f, line) != EOF; N++) { 54 stripwhite (line); 55 if (line[0] == 0) continue; 56 if (line[0] == '#') continue; 57 scanf (line, "%s %lf %lf %d", filename, &Rs, &Re, &Nrec); 58 if (Rs > patch[0].RA[1]) continue; 59 if (Re < patch[0].RA[0]) continue; 60 strcpy (regions[Nregions].filename, "%s/%s", TWO_MASS_DIR, filename); 61 regions[Nregions].RA[0] = Rs; 62 regions[Nregions].RA[1] = Re; 63 regions[Nregions].Nrec = Nrec; 64 Nregions ++; 65 CHECK_REALLOCATE (regions, TM_Region, NREGIONS, Nregions, 20); 66 } 67 68 *nregions = Nregions; 69 return (regions); 70 } 71 72 # define NBYTE 302 73 # define NLINE 100 74 75 Stars *get2mass_data (char *filename, GSCRegion *patch, int *nstars) { 76 77 int Nstars; 78 Stars *stars; 79 gzFile gf; 80 char buffer[NLINE*NBYTE]; 81 82 gf = gzopen (filename, "rb"); 83 if (gf == NULL) { 84 fprintf (stderr, "ERROR: can't read data from %s\n", file); 85 Shutdown (); 86 } 87 88 Nstars = 0; 89 NSTARS = 10000; 90 ALLOCATE (stars, Stars, NSTARS); 91 92 /* z_off_t gzseek (gzFile file, z_off_t offset, int whence); */ 93 94 while ((Nbyte = gzread (gf, buffer, NLINE*NBYTE)) != 0) { 95 96 if (Nbyte == -1) { 97 fprintf (stderr, "error reading from gzipped file %s\n", filename); 98 Shutdown (); 99 } 100 101 Nfrac = Nbyte % NBYTE; 102 if (Nfrac) { 103 fprintf (stderr, "error reading complete line from gzipped file %s\n", filename); 104 Shutdown (); 105 } 106 107 Nline = Nbyte / NBYTE; 108 109 for (i = 0; i < Nline; i++) { 110 Nval = scanf (&buffer[NBYTE*i + 0], "%lf %lf", &RA, &DEC); 111 if (RA < patch[0].RA[0]) continue; 112 if (RA > patch[0].RA[1]) continue; 113 if (DEC < patch[0].DEC[0]) continue; 114 if (DEC < patch[0].DEC[1]) continue; 115 116 Nval += scanf (&buffer[NBYTE*i + 53], "%lf %lf", &J, &dJ); 117 Nval += scanf (&buffer[NBYTE*i + 72], "%lf %lf", &H, &dH); 118 Nval += scanf (&buffer[NBYTE*i + 91], "%lf %lf", &K, &dK); 119 Nval += scanf (&buffer[NBYTE*i + 164], "%s", datestr); 120 if (Nval != 9) { 121 fprintf (stderr, "error in 2mass data line\n"); 122 exit (1); 123 } 124 125 time = short_date_to_set (datestr); 126 127 stars[Nstars].R = RA; 128 stars[Nstars].D = DEC; 129 stars[Nstars].M = J; 130 stars[Nstars].dM = dJ; 131 stars[Nstars].t = time; 132 stars[Nstars].code = TM_J; 133 Nstars ++; 134 135 stars[Nstars].R = RA; 136 stars[Nstars].D = DEC; 137 stars[Nstars].M = H; 138 stars[Nstars].dM = dH; 139 stars[Nstars].t = time; 140 stars[Nstars].code = TM_H; 141 Nstars ++; 142 143 stars[Nstars].R = RA; 144 stars[Nstars].D = DEC; 145 stars[Nstars].M = K; 146 stars[Nstars].dM = dK; 147 stars[Nstars].t = time; 148 stars[Nstars].code = TM_K; 149 Nstars ++; 150 151 CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars - 2, 300); 152 } 153 } 154 155 *nstars = Nstars; 156 return (stars); 157 } 158 159 /* this just scans along in the file. file is sorted by dec, so we 160 should be skipping large chunks - but we would need to have 161 the size from the accel file (won't fit in GSCRegion) and need 162 to use gzseek, if it exists. 163 */ -
trunk/Ohana/src/addstar/src/getgsc.c
r3347 r3361 1 # include "mosastro.h" 1 # include "addstar.h" 2 # define BYTES_STAR 23 3 # define BLOCK 1000 4 # define DNSTARS 1000 2 5 3 Star Data *getgsc (CatStats *catstats, int *NSTARS) {6 Stars *getgsc (GSCRegion *patch, int *NSTARS) { 4 7 5 8 int i, j, k, Ns, Ngsc, Nregions, Nstars; 6 9 StarData *gsc; 7 10 StarData *stars; 8 CatStats *regions; 11 GSCRegions *regions; 12 short int GSC_M; 9 13 14 NAMED_PHOTCODE (GSC_M, "GSC_M"); 15 16 regions = gregion_patch (patch, &Nregions); 17 10 18 Nstars = 0; 11 ALLOCATE (stars, StarData, 1); 12 13 set_catalog (GSC_DIR); 14 regions = gregions (catstats, &Nregions); 19 ALLOCATE (stars, Stars, 1); 15 20 16 21 for (i = 0; i < Nregions; i++) { 17 gsc = gcatalog (regions[i].name, &Ngsc); 18 if (gsc == NULL) continue; 22 gsc = rd_gsc (regions[i].name, &Ngsc); 19 23 20 Ns = Nstars; 24 REALLOCATE (stars, Stars, MAX (1, Nstars + Ngsc)); 25 memcpy (&stars[Nstars], gsc, Ngsc*sizeof(Stars)); 21 26 Nstars += Ngsc; 22 27 23 REALLOCATE (stars, StarData, Nstars);24 for (k = Ns, j = 0; j < Ngsc; k++, j++) {25 stars[k].R = gsc[j].R;26 stars[k].D = gsc[j].D;27 stars[k].Mag = gsc[j].Mag;28 }29 28 free (gsc); 30 29 } … … 35 34 } 36 35 36 Stars *rd_gsc (char *filename, int *nstars) { 37 38 int i, NSTARS, Nstars, Nbytes, nbytes; 39 char *buffer; 40 FILE *f; 41 Stars *stars; 42 43 Nstars = 0; 44 NSTARS = 1000; 45 ALLOCATE (stars, Stars, NSTARS); 46 47 f = fopen (filename, "r"); 48 if (f == NULL) { 49 fprintf (stderr, "ERROR: can't find catalog file %s\n", filename); 50 *nstars = Nstars; 51 return (stars); 52 } 53 54 Nbytes = BLOCK*BYTES_STAR; 55 ALLOCATE (buffer, char, Nbytes); 56 while ((nbytes = fread (buffer, 1, Nbytes, f)) > 0) { 57 Nline = nbytes / BYTES_STAR; 58 for (i = 0; i < Nline; i++) { 59 Nbyte = i*BYTES_STAR; 60 dparse (&stars[Nstars].R, 1, &buffer[Nbyte]); 61 dparse (&stars[Nstars].D, 2, &buffer[Nbyte]); 62 dparse (&stars[Nstars].M, 3, &buffer[Nbyte]); 63 stars[Nstars].dM = NO_ERR; 64 stars[Nstars].t = 0; 65 stars[Nstars].code = GSC_M; 66 Nstars ++; 67 CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 1000); 68 } 69 } 70 free (buffer); 71 REALLOCATE (stars, Stars, Nstars); 72 *nstars = Nstars; 73 return (stars); 74 } -
trunk/Ohana/src/addstar/src/getusno.c
r3347 r3361 1 # include " mosastro.h"1 # include "addstar.h" 2 2 # define NZONE 24 3 3 … … 9 9 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; 10 10 11 Stars *getusno ( CatStats*catstats, int *Nstars) {11 Stars *getusno (GSCRegion *catstats, int *Nstars) { 12 12 13 13 long int offset; … … 18 18 FILE *f; 19 19 int iRA0, iRA1, iDEC0, iDEC1; 20 double RA0, RA1, DEC0, DEC1,dec;20 double dec; 21 21 int spd, spd_start, spd_end, disk; 22 22 int NUSNO, Nusno; 23 23 Stars *stars; 24 short int USNO_RED, USNO_BLUE; 24 25 25 /** USNO-specific things **/ 26 USNO_RED = GetPhotcodebyName ("USNO_RED"); 27 if (!USNO_RED) { 28 fprintf (stderr, "ERROR: USNO_RED photcode not found in photcode table\n"); 29 exit (0); 30 } 31 USNO_BLUE = GetPhotcodebyName ("USNO_BLUE"); 32 if (!USNO_BLUE) { 33 fprintf (stderr, "ERROR: USNO_BLUE photcode not found in photcode table\n"); 34 exit (0); 35 } 36 37 RA0 = catstats[0].RA[0]; 38 RA1 = catstats[0].RA[1]; 39 DEC0 = catstats[0].DEC[0]; 40 DEC1 = catstats[0].DEC[1]; 26 NAMED_PHOTCODE (USNO_RED, "USNO_RED"); 27 NAMED_PHOTCODE (USNO_BLUE, "USNO_BLUE"); 41 28 42 29 /* identify ra & dec range of interest */ 43 iRA0 = RA0* 360000.0;44 iRA1 = RA1* 360000.0;45 iDEC0 = ( DEC0+ 90.0) * 360000.0;46 iDEC1 = ( DEC1+ 90.0) * 360000.0;30 iRA0 = catstats[0].RA[0] * 360000.0; 31 iRA1 = catstats[0].RA[1] * 360000.0; 32 iDEC0 = (catstats[0].DEC[0] + 90.0) * 360000.0; 33 iDEC1 = (catstats[0].DEC[1] + 90.0) * 360000.0; 47 34 48 35 /* data is organized in south-pole distance zones */ 49 spd_start = (int)(( DEC0+ 90) / 7.5) * 75.0;50 dec = ( DEC1+ 90) / 7.5;36 spd_start = (int)((catstats[0].DEC[0] + 90) / 7.5) * 75.0; 37 dec = (catstats[0].DEC[1] + 90) / 7.5; 51 38 if (dec > (int)(dec)) { 52 spd_end = (int)(1 + ( DEC1+ 90) / 7.5) * 75.0;39 spd_end = (int)(1 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; 53 40 } else { 54 spd_end = (int)(0 + ( DEC1+ 90) / 7.5) * 75.0;41 spd_end = (int)(0 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; 55 42 } 56 43 … … 89 76 fclose (f); 90 77 91 first = RA0/ 3.75;92 if (( RA1 / 3.75) == (int) (RA1/ 3.75))93 last = RA1/ 3.75;78 first = catstats[0].RA[0] / 3.75; 79 if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) 80 last = catstats[0].RA[1] / 3.75; 94 81 else 95 last = 1 + RA1/ 3.75;82 last = 1 + catstats[0].RA[1] / 3.75; 96 83 97 84 if ((first > Nbins) || (last > Nbins)) { … … 123 110 /* print out data from slice within RA and DEC range */ 124 111 for (i = 0; i < number[bin]; i++, buf+=3) { 125 if ((buf[0] > iRA0) && (buf[0] < iRA1) && 126 (buf[1] > iDEC0) && (buf[1] < iDEC1)) { 127 bzero (&stars[Nusno], sizeof(Stars)); 128 stars[Nusno].R = buf[0]/360000.0; 129 stars[Nusno].D = buf[1]/360000.0 - 90.0; 130 stars[Nusno].Mag = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); 131 stars[Nusno].code = USNO_RED; 132 /* r data */ 133 Nusno ++; 134 stars[Nusno].R = buf[0]/360000.0; 135 stars[Nusno].D = buf[1]/360000.0 - 90.0; 136 stars[Nusno].Mag = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000)); 137 stars[Nusno].code = USNO_BLUE; 138 /* b data */ 139 Nusno ++; 140 CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno-1, 5000); 141 } 112 if (buf[0] < iRA0) continue; 113 if (buf[0] > iRA1) continue; 114 if (buf[1] < iDEC0) continue; 115 if (buf[1] > iDEC1) continue; 116 117 bzero (&stars[Nusno], sizeof(Stars)); 118 stars[Nusno].R = buf[0]/360000.0; 119 stars[Nusno].D = buf[1]/360000.0 - 90.0; 120 stars[Nusno].M = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); 121 stars[Nusno].dM = NO_ERR; 122 stars[Nusno].t = 0; 123 stars[Nusno].code = USNO_RED; 124 /* r data */ 125 Nusno ++; 126 stars[Nusno].R = buf[0]/360000.0; 127 stars[Nusno].D = buf[1]/360000.0 - 90.0; 128 stars[Nusno].M = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000)); 129 stars[Nusno].dM = NO_ERR; 130 stars[Nusno].t = 0; 131 stars[Nusno].code = USNO_BLUE; 132 /* b data */ 133 Nusno ++; 134 135 CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno-1, 5000); 142 136 } 143 137 free (buffer); -
trunk/Ohana/src/addstar/src/gimages.c
r3339 r3361 126 126 } 127 127 128 int edge_check (double *x1, double *y1, double *x2, double *y2) {129 130 double theta1, theta2;131 double Theta1, Theta2;132 133 theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]);134 theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]);135 136 if (theta1*theta2 < 0.0) {137 return (FALSE);138 }139 140 if (fabs(theta1) < fabs(theta2)) {141 return (FALSE);142 }143 144 Theta1 = theta1;145 Theta2 = theta2;146 theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]);147 theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]);148 149 150 if (theta1*theta2 < 0.0) {151 return (FALSE);152 }153 154 if (fabs(theta1) < fabs(theta2)) {155 return (FALSE);156 }157 158 return (TRUE);159 160 }161 162 double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) {163 164 double dx1, dy1, dx2, dy2, ct, st, theta;165 166 dx1 = x1 - x2;167 dy1 = y1 - y2;168 169 dx2 = x3 - x2;170 dy2 = y3 - y2;171 172 ct = (dx1*dx2 + dy1*dy2);173 st = (dx1*dy2 - dx2*dy1);174 175 theta = atan2 (st, ct);176 177 return (theta);178 179 }180 -
trunk/Ohana/src/addstar/src/greference.c
r3347 r3361 1 1 # include "addstar.h" 2 2 3 Stars *grefcat (char *Refcat, CatStats *catstats, int *Nrefcat) {3 Stars *grefcat (char *Refcat, GSCRegion *catstats, int *nstars) { 4 4 5 5 int Nstars; … … 7 7 8 8 if (VERBOSE) fprintf (stderr, "loading reference catalog data from %s\n", Refcat); 9 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]); 9 10 10 11 Nstars = 0; 11 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]);12 12 13 13 /* get stars from USNO for the given region */ … … 30 30 exit (1); 31 31 } 32 * Nrefcat= Nstars;32 *nstars = Nstars; 33 33 return (stars); 34 34 -
trunk/Ohana/src/addstar/src/grefstars.c
r3347 r3361 1 1 # include "addstar.h" 2 # include <sys/time.h>3 # include <time.h>4 2 5 3 /* read ASCII file with ref star data */ … … 38 36 while (stars[N].R < 0.0) stars[N].R += 360.0; 39 37 while (stars[N].R >= 360.0) stars[N].R -= 360.0; 40 stars[N]. dophot = 0;41 stars[N]. found = -1;42 stars[N]. code = thiscode[0].code;;38 stars[N].t = 0; 39 stars[N].code = thiscode[0].code; 40 stars[N].found = FALSE; 43 41 CHECK_REALLOCATE (stars, Stars, NSTARS, N+1, 100); 44 42 } -
trunk/Ohana/src/addstar/src/gregion_star.c
r3347 r3361 1 1 # include "addstar.h" 2 2 3 GSCRegion *gregion_stars (Stars *stars, int Nstars, int * Nregion) {3 GSCRegion *gregion_stars (Stars *stars, int Nstars, int *nregions) { 4 4 5 5 int i; … … 11 11 Nregions = 0; 12 12 NREGIONS = 100; 13 ALLOCATE (region , GSCRegion, NREGIONS);13 ALLOCATE (regions, GSCRegion, NREGIONS); 14 14 15 15 sort_stars (stars, Nstars); … … 27 27 Nregions ++; 28 28 CHECK_REALLOCATE (regions, GSCRegion, Nregions, NREGIONS, 100); 29 29 30 } 30 31 *nregions = Nregions; 32 return (regions); 31 33 } 32 34 -
trunk/Ohana/src/addstar/src/gstars.c
r3347 r3361 3 3 # define BYTES_STAR 66 4 4 # define BLOCK 1000 5 #include <sys/time.h>6 #include <time.h>7 5 8 6 Stars *gstars (char *file, int *NSTARS, Image *image) { … … 142 140 ALLOCATE (stars, Stars, image[0].nstar); 143 141 144 /* temporary: get sky background from Flips data */ 142 /* this is really lame. get sky background from Flips data */ 143 /*** the problems are: 144 1) need to define sky entry in image table 145 2) need to have sky measurement (in reg.db) 146 (could be derived from star data?) 147 ***/ 148 145 149 if (SKYPROBE) { 146 150 char *p; … … 205 209 dparse (&stars[N].Y, 2, &buffer[j*BYTES_STAR]); 206 210 dparse (&stars[N].M, 3, &buffer[j*BYTES_STAR]); 207 dparse (&stars[N].dM, 4, &buffer[j*BYTES_STAR]); 211 212 /* cmp files carry dM in millimags */ 213 dparse (&tmp, 4, &buffer[j*BYTES_STAR]); 214 stars[N].dM = 0.001*tmp; 215 208 216 /* allow for some dynamic filtering of star list */ 209 217 if (stars[N].dM > 1000.0 / SNLIMIT) { N--; continue; } … … 215 223 dparse (&tmp, 5, &buffer[j*BYTES_STAR]); 216 224 stars[N].dophot = tmp; 225 217 226 dparse (&stars[N].Mgal, 7, &buffer[j*BYTES_STAR]); 218 227 dparse (&stars[N].Map, 8, &buffer[j*BYTES_STAR]); … … 225 234 while (stars[N].R >= 360.0) stars[N].R -= 360.0; 226 235 stars[N].found = -1; 236 stars[N].code = thiscode[0].code; 227 237 228 238 if (SKYPROBE) { -
trunk/Ohana/src/addstar/src/image-db.c
r3339 r3361 1 1 # include "addstar.h" 2 2 3 /* these variables are needed by Shutdown */ 3 /* file-pointer to Image.db table. if image db does not yet exist, 4 5 */ 4 6 static FILE *f = (FILE *) NULL; 7 static int dbstate = LCK_UNLOCK; 8 9 FILE *GetDB (int *state) { 10 *state = dbstate; 11 return (f); 12 } 5 13 6 14 /* clean up open / locked ImageCat before shutting down */ 7 15 int Shutdown () { 8 int dbstate;9 10 16 SetProtect (TRUE); 11 17 fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); 12 18 fprintf (stderr, "ERROR: addstar halted\n"); 13 19 exit (1); 14 }15 16 FILE *GetDB () {17 return (f);18 20 } 19 21 … … 29 31 } 30 32 31 void unlock_image_db () { 33 void unlock_image_db (Image *image) { 34 32 35 /* protect wimage from interrupt signals */ 33 SetProtect (TRUE); 34 wimage (f, dbstate, &image, Nstars); 35 SetProtect (FALSE); 36 if (mode == M_IMAGE) { 37 SetProtect (TRUE); 38 wimage (image); 39 SetProtect (FALSE); 40 } 36 41 37 42 fclearlockfile (ImageCat, f, LCK_HARD, &dbstate); 38 39 43 mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; 40 44 chmod (ImageCat, mode); -
trunk/Ohana/src/addstar/src/mkcatalog.c
r2457 r3361 1 1 # include "addstar.h" 2 # include <time.h>3 2 4 3 void mkcatalog (GSCRegion *region, Catalog *catalog) { -
trunk/Ohana/src/addstar/src/wimage.c
r3339 r3361 1 1 # include "addstar.h" 2 2 3 void wimage (FILE *f, int dbstate, Image *image, int Nstars) { 3 /* add image to image db */ 4 void wimage (Image *image) { 4 5 5 int Nimages, status, offset; 6 FILE *f; 7 int Nimages, status, offset, dbstate; 6 8 Header header; 7 9 8 /* place file pointer at beginning of file */10 f = GetDB (&dbstate); 9 11 fseek (f, 0, SEEK_SET); 10 12 … … 17 19 /* assign Zero Point provided by first image in dB */ 18 20 fits_modify (&header, "ZERO_PT", "%lf", 1, ZeroPt); 21 fits_modify (&header, "NIMAGES", "%d", 1, 0); 19 22 } else { 20 23 if (!fits_load_header (f, &header)) {
Note:
See TracChangeset
for help on using the changeset viewer.
