Index: /tags/addstar-1-2/Ohana/src/addstar/Makefile =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/Makefile (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/Makefile (revision 4773) @@ -0,0 +1,133 @@ +default: addstar +help: + @echo "make options: addstar (default)" + +include ../../Configure + +HOME = $(ROOT)/src/addstar +BIN = $(HOME)/bin +INC = $(HOME)/include +SRC = $(HOME)/src +MAN = $(HOME)/doc +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +# +INCS = -I$(INC) -I$(LINC) -I$(XINC) +LIBS = -L$(LLIB) -lohana -lFITS -lz -lm +CFLAGS = $(INCS) -g +LFLAGS = $(LIBS) + +ADDSTAR = \ +$(SRC)/addstar.$(ARCH).o \ +$(SRC)/airmass.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/calibrate.$(ARCH).o \ +$(SRC)/check_permissions.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/conversions.$(ARCH).o \ +$(SRC)/dump.$(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)/get2mass_as.$(ARCH).o \ +$(SRC)/get2mass_dr2.$(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)/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 \ +$(SRC)/rfits.$(ARCH).o \ +$(SRC)/rtext.$(ARCH).o + +MKACC-2MASS = \ +$(SRC)/mkacc-2mass.$(ARCH).o \ +$(SRC)/image-db.$(ARCH).o \ +$(SRC)/check_permissions.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/conversions.$(ARCH).o \ +$(SRC)/wimage.$(ARCH).o \ +$(SRC)/get2mass_as.$(ARCH).o + +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 : $(BIN)/addstar.$(ARCH) +$(BIN)/addstar.$(ARCH) : $(ADDSTAR) +# $(ADDSTAR): $(INC)/addstar.h + +gztest : $(BIN)/gztest.$(ARCH) +$(BIN)/gztest.$(ARCH) : $(SRC)/gztest.$(ARCH).o + +mkacc-2mass : $(BIN)/mkacc-2mass.$(ARCH) +$(BIN)/mkacc-2mass.$(ARCH) : $(MKACC-2MASS) +mkacc-2mass.install : $(DESTBIN)/mkacc-2mass + +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) $^ -o $@ $(LFLAGS) + +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* + +$(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) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` Index: /tags/addstar-1-2/Ohana/src/addstar/bin/.cvsignore =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/bin/.cvsignore (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/bin/.cvsignore (revision 4773) @@ -0,0 +1,1 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix Index: /tags/addstar-1-2/Ohana/src/addstar/doc/Changes.txt =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/doc/Changes.txt (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/doc/Changes.txt (revision 4773) @@ -0,0 +1,25 @@ + +2005.08.15: + cleanup of the minor Wall,Werror messages + +2005.02.24 : I am adding the mosaic two-level astrometry information to the Image table. + +2005.07.06 : current release is addstar-1.1 + I have made a variety of fairly substantial changes since + v1.0. I have added autocode for the basic DVO data types. I + have merged the functionality of addspphot, addrefs, and + addusno into the single addstar function. There are now + options to add an ascii text table or stars from known + external reference catalogs (HST GSC, USNO, 2MASS). Also, the + functions added for skyprobe (in-flight Mcal calculation, + color terms) are also incorporated. I also substantially + re-organized the basic code to make the different modes more + sensible. There are still a few changes needed for addstar + v2.0 (which will be able to handle the pan-starrs bandwidth). + specifically, addstar needs to handle the catalog files using + the TBD sky.db interfaces, it needs to interact with a DVO + server in a client/server way (to divide up the sky), it needs + to handle multiple versions of the table formats, and it might + need to know about the mysql backend interface. + + Index: /tags/addstar-1-2/Ohana/src/addstar/doc/database.txt =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/doc/database.txt (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/doc/database.txt (revision 4773) @@ -0,0 +1,91 @@ + +\begin{verbatim} + +typedef struct { + float R; /* RA in decimal degrees */ + float D; /* DEC in decimal degrees */ + short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */ + unsigned short int Nm; /* number of measurements */ + unsigned short int Xp, Xm; /* chisq values in tenths */ + unsigned int offset; /* offset to first Measure-ment */ +} Average; /* = 20 bytes / average */ + +typedef struct { + char dR, dD; /* tenths of arcsec (-12.7 to +12.7 valid range) */ + short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */ + unsigned char dM; /* thousandths of mag (0.000 -- 0.255 valid range) */ + float t; /* time in seconds (what is reference?) */ + unsigned int average; /* reference to corresponding Average entry, + upper byte of value contains flags. + limit of 16,777,215 stars (Naverage) + in a file (=0xFFFFFF). + flags = average & 0x1000000 */ +} Measure; /* = 13 bytes / measure */ + +# define BLEND_IMAGE 0X01000000 +# define BLEND_CATALOG 0X02000000 +# define UPPER_LIMIT 0X04000000 +# define CALIBRATED 0X08000000 + +\end{verbatim} + +The above two structures define the entries in the photometry +database. The database consists of a large number of files +representing a small patch on the sky (roughly 1.5 degree$^2$ in most +places). These files are organized into directories representing +bands of Declination. A reference file determines the coordinate +boundaries for each of the files so that a given point on the sky can +unambiguously be associated with a specific file in a specific +directory. The sky coordinates for each file is the same as those +used by the HST Guide Star catalog, except for the region around the +North celestial pole, for which all stars are included in a single +file. + +Within a given file, the data are stored in a binary format, with an +ASCII FITS-like header. The header is examply in the format of a +normal FITS header, but with the exception that all files have a fixed +number of header blocks (for now 3 blocks = 8640 bytes). This is done +to speed loading the header and finding the beginning of the binary +data. The number of 3 blocks seems quite generous, as currently only +a few FITS keywords have been defined for each file, basically +keywords to define the number of stars and the total number of +measurements stored in the file, as well as values to define the RA +and DEC range of the file. + +The first section of data following the header blocks consists of +average measurements for each uniquely observed star. Each star +occupies 20 bytes, the size of the Average structure defined above. +The Average structure contains the average Ra, Dec, and Magnitude for +the star, as well as the number of measurements, and \chisq\ values +for the magnitude and position. Finally, there is a 32 bit integer +which defines the offset to the first measurement for this star. This +offset is defined as the number of Measure records from the start of +the Measure structure. + +The second section of data, following the Average data contains all +measurements for each star listed in Average. Each measurement +occupies 13 bytes, the size of the Measure structure. This structure +contains the difference of this position from the average RA and DEC, +and the instrumental magnitude of this measurement (in the units +defined by the fstat program, which give m = -2.5*log(cts) + Mo, where +Mo is currently 24.5 [10/15]). There is also the magnitude error for +this measurement, the time of the measurement (in seconds relative to +a to-be-determined zero point), and a reference to the entry in the +Average structure so we can relate a given measurement with a given star. +This last entry also includes a byte of flags, of which only 4 have +currently been defined. This means the Average offset can only be as +large as 16,777,215 (0xffffff), limiting the possible number of stars +allowed in a given file. This does not seem like a long term problem, +though: aside from the fact that this number is very large and we +only expect in the vicinity of 20,000 stars per file, the file can +easily be divided into pieces at a later date if needed. this last +step is trivial, consisting of splitting the data up into smaller +RA,DEC regions, and updating the reference catalog. With the above +definitions for the Average and Measure structures, we will typically +expect 20,000 * 15 measurements per year and 20,000 average entries in +a given file. This implies a file size of about 4.3 MB at the end of +a year. It is also possible that we will choose to split the files up +in the future if the number of measurements makes their size +unweildy. 4.3 MB is sufficiently small that this is not a problem, +but after only 5 years, the files will be over 21 MB each, getting to +be fairly significant. Index: /tags/addstar-1-2/Ohana/src/addstar/doc/description.txt =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/doc/description.txt (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/doc/description.txt (revision 4773) @@ -0,0 +1,94 @@ + +ADDSTAR: + +This program takes a photometry file, with astrometry, provided by +earlier routines and incorporates the stellar measurements into the +photometry database. This routine does NOT try to calibrate or check +the photometry, nor does it try to improve the coordinate +determinations or calculate chisquares for position or photometry - +those tasks are performed by other routines. + +Addstar starts by loading in the photometry file and converting the +pixel coordinates to sky coordinates using the astrometry information +in the header of the file. It also gets the detection limits and +stores the image astrometry information. + +Next, addstar determines which of all the previous images overlaps +with this image and also which of the catalog regions overlaps with +this image. + +Addstar considers each catalog region in turn, and compares the +positions of stars in the catalog with stars in the image. If a +catalog star is detected, a new measurement is added to the catalog. +If a catalog star is not detected, but is in the field of view of this +image, the detection limit for this image is added to the list of +measurements as an upper limit. If the image contains a star which is +not already included in the catalog, the star is added to the catalog, +and all previous images are checked to see if they included this +location. If so, upper limits are added to the list of measurements. + +To deal with blended images, we have taken the philosophy that all +measurements compatible with the coordinates of a given star should be +included in the list of measurements. This assumes that programs or people +downstream will figure out which is the "correct" measurement. To +this end, if a catalog star is consistent with multiple measurements, +they are all added to the catalog measurement and the measurements are +given a flag to say that there was blending in the catalog (ie, it is +probably the catalog star which is a blend). If multiple catalog +stars are consistent with the same image star, that measurement is +added to each catalog star, and given a flag to say that there was +blending on the image. + +After all catalogs have been checked, updated, and written to disk, +the image is added to the database of images. + +The format for a catalog region file is as follows: + +-- +Header +-- +Averages +-- +Measurements +-- + +The Header follows the format of a FITS header but always stored with +3 blocks to speed up read time: 2880 x 3 bytes. Two important +keywords in the header: NSTARS & NMEAS + +All average values are stored consecutively in a single block. They +are written as an array with NSTARS elements of the Average structure +defined below. The individual measurements follow the averages as a +block. To find a measurement for a given star, you need to get the +value of the offset for the star. this is the element number for the +first measurement of this star, and are followed by next measurements +for this star, until the value Nm is reached. + +Here are the structures being used for the average and measurement +values: + +typedef struct { + float R; /* RA in decimal degrees */ + float D; /* DEC in decimal degrees */ + short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */ + unsigned short int Nm; /* number of measurements */ + unsigned short int Xp, Xm; /* chisq values in tenths */ + unsigned int offset; /* offset to first measurement */ +} Average; /* = 20 bytes / average */ + +typedef struct { + char dR, dD; /* tenths of arcsec (-12.7 to +12.7 valid range) */ + short int M; /* thousandths of mag (-32.000 to 32.000 valid range) */ + unsigned char dM; /* thousandths of mag (0.000 -- 0.255 valid range) */ + float t; /* time in seconds (what is reference?) */ + unsigned int average; + /* reference to corresponding average entry, upper byte stores flags: + limit of 16,777,215 stars (Naverage) in a file (=0xFFFFFF) + flags: average & 0x1000000 */ +} Measure; /* = 13 bytes / measure */ + +/* flags for measurements: */ +# define BLEND_IMAGE 0X01000000 +# define BLEND_CATALOG 0X02000000 +# define UPPER_LIMIT 0X04000000 +# define CALIBRATED 0X08000000 Index: /tags/addstar-1-2/Ohana/src/addstar/doc/merge.txt =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/doc/merge.txt (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/doc/merge.txt (revision 4773) @@ -0,0 +1,30 @@ + +the goal of this project is to merge the functionality of the four +addstar-like programs in DVO: addstar, addrefs, addspphot, addusno. + +files which are identical or could be in all four: + +addstar/src/aregion.c +addstar/src/addstar.c +addstar/src/args.c +addstar/src/calibrate.c +addstar/src/check_permissions.c +addstar/src/ConfigInit.c +addstar/src/find_matches.c +addstar/src/gcatalog.c +addstar/src/gimages.c +addstar/src/gregions.c +addstar/src/gstars.c +addstar/src/mkcatalog.c +addstar/src/parse_time.c +addstar/src/sort_lists.c +addstar/src/wcatalog.c +addstar/src/wimage.c + +files to watch: +[eugene@kawelu src]$ ll addstar/src/ +-rw-rw-r-- 1 eugene eugene 3658 Nov 24 14:13 calibrate.c + +[eugene@kawelu src]$ ll addspphot/src/ +total 116 +-rw-rw-r-- 1 eugene eugene 4480 Nov 24 15:18 calibrate.c Index: /tags/addstar-1-2/Ohana/src/addstar/doc/notes.txt =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/doc/notes.txt (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/doc/notes.txt (revision 4773) @@ -0,0 +1,41 @@ + +2005.03.07 : notes related to new version of addstar + +- USNO seems to work fine now: + + addstar -region 9.0 12.0 9.0 12.0 -cat usno -p USNO_RED + addstar -region 9.0 12.0 9.0 12.0 -cat usno -p USNO_BLUE + +- 2MASS seems to work as well: + + addstar -region 3.0 12.0 12.0 14.0 -cat 2mass -p 2MASS_J + addstar -region 3.0 12.0 12.0 14.0 -cat 2mass -p 2MASS_H + addstar -region 3.0 12.0 12.0 14.0 -cat 2mass -p 2MASS_K + +- both of these are somewhat inefficient since they need to search + through the (larger) USNO/2MASS data files for each output ptolemy + region file. + +- GSC seems to work fine now: + + addstar -region 3.0 12.0 12.0 14.0 -cat gsc + (photcode is not needed here since there is only one) + +- images seems to work fine now + +- general concerns: + + - RA 0,360 boundary is likely to have problems, especially with + region selection + - mosaic astrometry concepts still need to be added. + + - average magnitudes are being set incorrectly in addstar (not + taking account of the zero-points). I don't really want to apply + the colors here, but I need to do something, perhaps? + + one option is to only apply the correction for the first + measurement of an object, unlike the way astrometry works + currently. both of these should be updated based on + relphot/uniphot and eqiuvalent relastro/uniastro independently + from the addstar step. + Index: /tags/addstar-1-2/Ohana/src/addstar/etc/phottemp.cat =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/etc/phottemp.cat (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/etc/phottemp.cat (revision 4773) @@ -0,0 +1,296 @@ +SIMPLE = F / ASTROCAM PHOTOMETRY FILE \ +BITPIX = 16 / \ +NAXIS = 2 / NUMBER OF AXES \ +NAXIS1 = 1106 / NUMBER OF COLUMNS \ +NAXIS2 = 1024 / NUMBER OF ROWS \ +BSCALE = 1.000000 / \ +BZERO = 0.000000 / \ +DATE = '19/12/94' / UT Date of file creation (DD/MM/YY) \ +ORIGIN = 'MDM Observatory' / Michigan-Dartmouth-MIT \ +LATITUDE= 31.9500 / Latitude (degrees N) \ +LONGITUD= -111.6150 / Longitude (degrees E) \ +OBSERVER= 'Metzger ' / Name of observer \ +TELESCOP= '1.3m McGraw-Hill' / Telescope used for observation \ +INSTRUME= 'Charlotte Direct' / Instrument used for observation \ +DETECTOR= 'Charlotte/Tek 1024^2 CCD' / Detector used for observation \ +FRAME = 32 / Frame number of observation \ +CCDPICNO= 32 / Frame number of observation \ +OBJECT = 'ocl0327 ' / Name of object \ +IMAGETYP= 'OBJECT ' / Type of observation \ +EXPTIME = 250.000 / Integration time (seconds) \ +DARKTIME= 250.067 / Dark current time (seconds) \ +DATE-OBS= '19/12/94' / UT Date of observation (DD/MM/YY) \ +UT = ' 05:25:44.00' / Universal time (UTC) at exposure start \ +JD = 2449705.726204 \ +RA = 50.700 / Right Ascension \ +DEC = 89.000 / Declination \ +DIRECTN = -1 / Moving South \ +EQUINOX = 1950.000 / Equinox of RA and DEC \ +HA = ' 02:14:06.72' / Hour angle at start \ +ST = ' 03:49:36.85' / Sidereal time at start \ +ZD = ' 36:26:07.17' / Zenith distance (degrees) \ +AIRMASS = 1.243 / Airmass at start \ +FILTER = 'I KP ' / Filter description \ +GAIN = 3.350 / Nominal gain (e-/ADU) \ +SECPIX1 = 0.508 / Arcseconds per pixel in fast dir \ +SECPIX2 = 0.508 / Arcseconds per pixel in slow dir \ +CCDBIN1 = 1 / On-chip column binning (fast dir) \ +CCDBIN2 = 1 / On-chip row binning (slow dir) \ +GPROBE = ' 5000.00 7000.00' / Guide probe X Y \ +DATASEC = '[51:1074,1:1024]' / Image area of frame \ +CCDSEC = '[1:1074,1:1024]' / Image area relative to full chip \ +BIASSEC = '[1080:1106,1:1023]' / Overscan area of frame \ +UNSIGN = T / \ +NSTARS = 0 / NUMBER OF stars \ +CTYPE1 = 'RA---SIN ' / +CTYPE2 = 'DEC--SIN ' / +CDELT1 = 0.000733 / +CDELT2 = 0.000733 / +CRVAL1 = 50.191418 / +CRVAL2 = 88.998663 / +CRPIX1 = 1028.644173 / +CRPIX2 = 503.129830 / +PC001001= 0.999906 / +PC001002= 0.012081 / +PC002001= -0.011936 / +PC002002= 0.999950 / +END \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + 246.5 352.1 13.000 020 1 3.2 + 246.5 352.1 14.100 020 1 3.2 + 279.4 328.9 13.800 020 1 3.2 + 187.3 238.8 9.800 020 1 3.2 + 187.3 238.8 9.800 020 1 3.2 + 436.6 413.5 14.300 020 1 3.2 + 9.3 856.5 13.000 020 1 3.2 + 9.3 856.2 13.600 020 1 3.2 + 319.7 697.9 12.000 020 1 3.2 + 319.7 697.7 12.100 020 1 3.2 + 244.3 9.5 14.100 020 1 3.2 + 284.7 334.4 14.000 020 1 3.2 + 127.1 122.8 14.200 020 1 3.2 + 140.2 791.1 13.700 020 1 3.2 + 140.2 790.8 14.400 020 1 3.2 + 131.7 636.3 13.800 020 1 3.2 + 423.8 700.6 13.800 020 1 3.2 + 423.7 700.7 12.500 020 1 3.2 + 266.2 290.7 13.800 020 1 3.2 + 266.0 290.7 13.100 020 1 3.2 + 570.7 721.6 13.600 020 1 3.2 + 570.7 721.6 13.100 020 1 3.2 + 548.7 673.0 13.700 020 1 3.2 + 264.2 219.7 12.800 020 1 3.2 + 264.0 219.7 12.100 020 1 3.2 + 161.9 797.4 11.900 020 1 3.2 + 162.0 797.2 12.200 020 1 3.2 + 505.8 851.8 13.700 020 1 3.2 + 88.2 419.0 11.400 020 1 3.2 + 88.2 419.0 11.800 020 1 3.2 + 426.9 613.2 9.400 020 1 3.2 + 427.4 613.5 7.600 020 1 3.2 + 326.1 595.2 12.100 020 1 3.2 + 326.2 595.1 12.600 020 1 3.2 + 444.4 845.2 14.200 020 1 3.2 + 254.1 916.8 14.500 020 1 3.2 + 502.6 882.4 14.200 020 1 3.2 + 128.3 275.7 12.300 020 1 3.2 + 128.3 275.7 12.900 020 1 3.2 + 230.9 554.8 13.300 020 1 3.2 + 230.9 554.7 14.200 020 1 3.2 + 47.5 685.2 13.500 020 1 3.2 + 47.5 684.9 14.200 020 1 3.2 + 382.6 658.7 14.000 020 1 3.2 + 182.2 330.3 14.200 020 1 3.2 + 245.1 548.3 12.600 020 1 3.2 + 245.3 548.0 13.200 020 1 3.2 + 50.3 849.0 13.700 020 1 3.2 + 50.5 848.8 14.300 020 1 3.2 + 17.9 772.7 14.200 020 1 3.2 + 499.8 712.8 13.500 020 1 3.2 + 499.8 712.6 14.400 020 1 3.2 + 294.5 349.4 14.000 020 1 3.2 + 497.7 747.8 11.000 020 1 3.2 + 497.5 747.8 10.900 020 1 3.2 + 109.7 323.5 14.000 020 1 3.2 + 59.5 917.6 14.300 020 1 3.2 + 59.4 917.8 13.700 020 1 3.2 + 82.2 880.5 13.800 020 1 3.2 + 145.1 341.2 14.500 020 1 3.2 + 452.7 592.4 14.400 020 1 3.2 + 623.4 820.0 13.700 020 1 3.2 + 609.7 829.3 13.500 020 1 3.2 + 481.8 938.3 14.400 020 1 3.2 + 617.0 792.7 13.500 020 1 3.2 + 183.2 539.7 14.400 020 1 3.2 + 133.9 679.6 14.400 020 1 3.2 + 423.7 667.8 6.500 020 1 3.2 + 434.7 111.9 13.200 020 1 3.2 + 434.5 111.9 13.200 020 1 3.2 + 434.6 111.9 14.200 020 1 3.2 +1222.3 684.2 13.900 020 1 3.2 +1077.6 282.5 14.000 020 1 3.2 +1077.8 282.5 13.400 020 1 3.2 + 720.8 562.1 14.100 020 1 3.2 + 721.1 561.9 14.400 020 1 3.2 +1182.1 944.4 12.100 020 1 3.2 +1181.9 944.4 11.800 020 1 3.2 +1181.9 944.4 12.300 020 1 3.2 + 761.2 423.1 14.300 020 1 3.2 + 761.0 423.1 14.100 020 1 3.2 + 926.7 436.8 13.400 020 1 3.2 + 926.6 436.8 14.200 020 1 3.2 +1209.1 45.0 13.200 020 1 3.2 +1209.2 43.6 12.800 020 1 3.2 +1135.2 653.1 8.100 020 1 3.2 +1135.6 653.4 8.200 020 1 3.2 +1123.5 735.8 12.000 020 1 3.2 +1123.6 735.8 11.800 020 1 3.2 +1123.7 735.7 11.600 020 1 3.2 + 730.7 28.6 12.100 020 1 3.2 + 730.7 28.6 12.200 020 1 3.2 + 783.1 267.5 14.200 020 1 3.2 + 783.2 267.5 13.400 020 1 3.2 + 876.2 462.7 13.800 020 1 3.2 + 876.4 462.7 14.200 020 1 3.2 + 917.2 399.9 13.700 020 1 3.2 + 917.3 399.9 13.200 020 1 3.2 +1227.5 218.4 10.400 020 1 3.2 +1227.7 218.4 10.500 020 1 3.2 +1022.7 146.0 13.300 020 1 3.2 +1022.5 146.0 14.100 020 1 3.2 +1324.7 72.3 12.900 020 1 3.2 +1324.6 72.3 13.200 020 1 3.2 + 743.2 379.4 14.300 020 1 3.2 + 742.9 379.4 13.900 020 1 3.2 +1255.9 544.0 12.600 020 1 3.2 +1256.0 543.9 12.300 020 1 3.2 + 848.9 251.1 13.200 020 1 3.2 + 849.0 251.1 12.700 020 1 3.2 + 495.6 147.4 11.400 020 1 3.2 + 495.6 147.4 11.300 020 1 3.2 + 495.5 147.4 11.000 020 1 3.2 +1016.4 735.5 14.200 020 1 3.2 +1016.3 735.7 13.900 020 1 3.2 + 795.2 675.9 14.400 020 1 3.2 + 887.2 968.8 10.800 020 1 3.2 + 626.8 726.9 14.400 020 1 3.2 + 718.6 800.1 9.700 020 1 3.2 + 675.4 301.6 14.400 020 1 3.2 + 801.6 436.8 14.200 020 1 3.2 +1066.4 870.0 13.700 020 1 3.2 +1066.5 869.9 13.800 020 1 3.2 + 644.5 662.9 9.900 020 1 3.2 + 639.1 760.4 13.000 020 1 3.2 + 822.7 663.1 12.800 020 1 3.2 + 822.7 663.1 13.100 020 1 3.2 + 677.2 716.2 12.900 020 1 3.2 + 765.6 641.5 13.200 020 1 3.2 + 765.7 641.4 13.500 020 1 3.2 + 998.3 425.8 14.400 020 1 3.2 + 713.7 716.2 13.800 020 1 3.2 + 814.7 887.6 2.000 020 1 3.2 + 808.7 884.1 8.200 020 1 3.2 +1812.3 259.3 14.500 020 1 3.2 +1655.9 655.4 12.200 020 1 3.2 +1655.9 655.4 12.100 020 1 3.2 +1577.2 839.6 11.500 020 1 3.2 +1577.2 839.6 11.400 020 1 3.2 +1643.7 597.0 14.100 020 1 3.2 +1552.8 251.1 10.100 020 1 3.2 +1552.8 251.1 10.300 020 1 3.2 +1410.0 999.9 13.900 020 1 3.2 +1335.6 687.0 13.000 020 1 3.2 +1335.8 686.8 12.700 020 1 3.2 +1471.4 262.0 14.200 020 1 3.2 +1479.6 608.3 14.200 020 1 3.2 +1478.1 211.5 12.000 020 1 3.2 +1478.3 211.5 12.000 020 1 3.2 +1552.5 46.4 14.500 020 1 3.2 +2042.4 548.3 10.800 020 1 3.2 +2042.6 548.3 11.000 020 1 3.2 +1448.0 915.7 14.100 020 1 3.2 +1462.3 989.2 13.000 020 1 3.2 +1462.3 989.0 12.700 020 1 3.2 +1316.9 787.3 14.100 020 1 3.2 +1787.6 883.5 12.200 020 1 3.2 +1787.6 883.5 12.200 020 1 3.2 +1979.1 300.3 12.900 020 1 3.2 +1979.1 300.3 12.900 020 1 3.2 +1677.2 106.4 13.100 020 1 3.2 +1677.0 106.4 13.300 020 1 3.2 +1324.2 429.9 9.200 020 1 3.2 +1324.2 429.9 9.000 020 1 3.2 +1462.8 214.3 13.700 020 1 3.2 +1462.8 214.3 14.200 020 1 3.2 +1405.4 95.5 11.700 020 1 3.2 +1405.2 95.5 11.600 020 1 3.2 +1604.8 90.0 11.200 020 1 3.2 +1604.8 90.0 11.300 020 1 3.2 +1921.5 738.1 13.500 020 1 3.2 +1921.5 738.1 13.500 020 1 3.2 +1459.8 473.6 13.400 020 1 3.2 +1459.7 473.6 13.800 020 1 3.2 +1438.7 77.8 12.200 020 1 3.2 +1438.6 77.8 12.200 020 1 3.2 +1355.0 114.6 13.000 020 1 3.2 +1354.9 114.6 13.200 020 1 3.2 +1941.4 940.7 13.000 020 1 3.2 +1941.6 940.7 12.900 020 1 3.2 +1994.0 800.8 13.900 020 1 3.2 +1993.8 800.8 14.000 020 1 3.2 +2024.4 907.0 13.700 020 1 3.2 +2024.4 907.0 13.800 020 1 3.2 +2035.7 881.6 13.900 020 1 3.2 +2035.9 881.6 13.800 020 1 3.2 +1994.0 924.2 12.100 020 1 3.2 +1994.0 924.2 12.000 020 1 3.2 Index: /tags/addstar-1-2/Ohana/src/addstar/etc/template.cat =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/etc/template.cat (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/etc/template.cat (revision 4773) @@ -0,0 +1,296 @@ +SIMPLE = F / ASTROCAM PHOTOMETRY FILE \ +BITPIX = 16 / \ +NAXIS = 2 / NUMBER OF AXES \ +NAXIS1 = 1106 / NUMBER OF COLUMNS \ +NAXIS2 = 1024 / NUMBER OF ROWS \ +BSCALE = 1.000000 / \ +BZERO = 0.000000 / \ +DATE = '19/12/94' / UT Date of file creation (DD/MM/YY) \ +ORIGIN = 'MDM Observatory' / Michigan-Dartmouth-MIT \ +LATITUDE= 31.9500 / Latitude (degrees N) \ +LONGITUD= -111.6150 / Longitude (degrees E) \ +OBSERVER= 'Metzger ' / Name of observer \ +TELESCOP= '1.3m McGraw-Hill' / Telescope used for observation \ +INSTRUME= 'Charlotte Direct' / Instrument used for observation \ +DETECTOR= 'Charlotte/Tek 1024^2 CCD' / Detector used for observation \ +FRAME = 32 / Frame number of observation \ +CCDPICNO= 32 / Frame number of observation \ +OBJECT = 'ocl0327 ' / Name of object \ +IMAGETYP= 'OBJECT ' / Type of observation \ +EXPTIME = 250.000 / Integration time (seconds) \ +DARKTIME= 250.067 / Dark current time (seconds) \ +DATE-OBS= '19/12/94' / UT Date of observation (DD/MM/YY) \ +UT = ' 05:25:44.00' / Universal time (UTC) at exposure start \ +JD = 2449705.726204 \ +RA = 50.700 / Right Ascension \ +DEC = 89.000 / Declination \ +DIRECTN = -1 / Moving South \ +EQUINOX = 1950.000 / Equinox of RA and DEC \ +HA = ' 02:14:06.72' / Hour angle at start \ +ST = ' 03:49:36.85' / Sidereal time at start \ +ZD = ' 36:26:07.17' / Zenith distance (degrees) \ +AIRMASS = 1.243 / Airmass at start \ +FILTER = 'I KP ' / Filter description \ +GAIN = 3.350 / Nominal gain (e-/ADU) \ +SECPIX1 = 0.508 / Arcseconds per pixel in fast dir \ +SECPIX2 = 0.508 / Arcseconds per pixel in slow dir \ +CCDBIN1 = 1 / On-chip column binning (fast dir) \ +CCDBIN2 = 1 / On-chip row binning (slow dir) \ +GPROBE = ' 5000.00 7000.00' / Guide probe X Y \ +DATASEC = '[51:1074,1:1024]' / Image area of frame \ +CCDSEC = '[1:1074,1:1024]' / Image area relative to full chip \ +BIASSEC = '[1080:1106,1:1023]' / Overscan area of frame \ +UNSIGN = T / \ +NSTARS = 0 / NUMBER OF stars \ +CTYPE1 = 'RA---SIN ' / +CTYPE2 = 'DEC--SIN ' / +CDELT1 = 0.000733 / +CDELT2 = 0.000733 / +CRVAL1 = 50.191418 / +CRVAL2 = 88.998663 / +CRPIX1 = 1028.644173 / +CRPIX2 = 503.129830 / +PC001001= 0.999906 / +PC001002= 0.012081 / +PC002001= -0.011936 / +PC002002= 0.999950 / +END \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + \ + 246.5 352.1 13.000 020 1 3.2 + 246.5 352.1 14.100 020 1 3.2 + 279.4 328.9 13.800 020 1 3.2 + 187.3 238.8 9.800 020 1 3.2 + 187.3 238.8 9.800 020 1 3.2 + 436.6 413.5 14.300 020 1 3.2 + 9.3 856.5 13.000 020 1 3.2 + 9.3 856.2 13.600 020 1 3.2 + 319.7 697.9 12.000 020 1 3.2 + 319.7 697.7 12.100 020 1 3.2 + 244.3 9.5 14.100 020 1 3.2 + 284.7 334.4 14.000 020 1 3.2 + 127.1 122.8 14.200 020 1 3.2 + 140.2 791.1 13.700 020 1 3.2 + 140.2 790.8 14.400 020 1 3.2 + 131.7 636.3 13.800 020 1 3.2 + 423.8 700.6 13.800 020 1 3.2 + 423.7 700.7 12.500 020 1 3.2 + 266.2 290.7 13.800 020 1 3.2 + 266.0 290.7 13.100 020 1 3.2 + 570.7 721.6 13.600 020 1 3.2 + 570.7 721.6 13.100 020 1 3.2 + 548.7 673.0 13.700 020 1 3.2 + 264.2 219.7 12.800 020 1 3.2 + 264.0 219.7 12.100 020 1 3.2 + 161.9 797.4 11.900 020 1 3.2 + 162.0 797.2 12.200 020 1 3.2 + 505.8 851.8 13.700 020 1 3.2 + 88.2 419.0 11.400 020 1 3.2 + 88.2 419.0 11.800 020 1 3.2 + 426.9 613.2 9.400 020 1 3.2 + 427.4 613.5 7.600 020 1 3.2 + 326.1 595.2 12.100 020 1 3.2 + 326.2 595.1 12.600 020 1 3.2 + 444.4 845.2 14.200 020 1 3.2 + 254.1 916.8 14.500 020 1 3.2 + 502.6 882.4 14.200 020 1 3.2 + 128.3 275.7 12.300 020 1 3.2 + 128.3 275.7 12.900 020 1 3.2 + 230.9 554.8 13.300 020 1 3.2 + 230.9 554.7 14.200 020 1 3.2 + 47.5 685.2 13.500 020 1 3.2 + 47.5 684.9 14.200 020 1 3.2 + 382.6 658.7 14.000 020 1 3.2 + 182.2 330.3 14.200 020 1 3.2 + 245.1 548.3 12.600 020 1 3.2 + 245.3 548.0 13.200 020 1 3.2 + 50.3 849.0 13.700 020 1 3.2 + 50.5 848.8 14.300 020 1 3.2 + 17.9 772.7 14.200 020 1 3.2 + 499.8 712.8 13.500 020 1 3.2 + 499.8 712.6 14.400 020 1 3.2 + 294.5 349.4 14.000 020 1 3.2 + 497.7 747.8 11.000 020 1 3.2 + 497.5 747.8 10.900 020 1 3.2 + 109.7 323.5 14.000 020 1 3.2 + 59.5 917.6 14.300 020 1 3.2 + 59.4 917.8 13.700 020 1 3.2 + 82.2 880.5 13.800 020 1 3.2 + 145.1 341.2 14.500 020 1 3.2 + 452.7 592.4 14.400 020 1 3.2 + 623.4 820.0 13.700 020 1 3.2 + 609.7 829.3 13.500 020 1 3.2 + 481.8 938.3 14.400 020 1 3.2 + 617.0 792.7 13.500 020 1 3.2 + 183.2 539.7 14.400 020 1 3.2 + 133.9 679.6 14.400 020 1 3.2 + 423.7 667.8 6.500 020 1 3.2 + 434.7 111.9 13.200 020 1 3.2 + 434.5 111.9 13.200 020 1 3.2 + 434.6 111.9 14.200 020 1 3.2 +1222.3 684.2 13.900 020 1 3.2 +1077.6 282.5 14.000 020 1 3.2 +1077.8 282.5 13.400 020 1 3.2 + 720.8 562.1 14.100 020 1 3.2 + 721.1 561.9 14.400 020 1 3.2 +1182.1 944.4 12.100 020 1 3.2 +1181.9 944.4 11.800 020 1 3.2 +1181.9 944.4 12.300 020 1 3.2 + 761.2 423.1 14.300 020 1 3.2 + 761.0 423.1 14.100 020 1 3.2 + 926.7 436.8 13.400 020 1 3.2 + 926.6 436.8 14.200 020 1 3.2 +1209.1 45.0 13.200 020 1 3.2 +1209.2 43.6 12.800 020 1 3.2 +1135.2 653.1 8.100 020 1 3.2 +1135.6 653.4 8.200 020 1 3.2 +1123.5 735.8 12.000 020 1 3.2 +1123.6 735.8 11.800 020 1 3.2 +1123.7 735.7 11.600 020 1 3.2 + 730.7 28.6 12.100 020 1 3.2 + 730.7 28.6 12.200 020 1 3.2 + 783.1 267.5 14.200 020 1 3.2 + 783.2 267.5 13.400 020 1 3.2 + 876.2 462.7 13.800 020 1 3.2 + 876.4 462.7 14.200 020 1 3.2 + 917.2 399.9 13.700 020 1 3.2 + 917.3 399.9 13.200 020 1 3.2 +1227.5 218.4 10.400 020 1 3.2 +1227.7 218.4 10.500 020 1 3.2 +1022.7 146.0 13.300 020 1 3.2 +1022.5 146.0 14.100 020 1 3.2 +1324.7 72.3 12.900 020 1 3.2 +1324.6 72.3 13.200 020 1 3.2 + 743.2 379.4 14.300 020 1 3.2 + 742.9 379.4 13.900 020 1 3.2 +1255.9 544.0 12.600 020 1 3.2 +1256.0 543.9 12.300 020 1 3.2 + 848.9 251.1 13.200 020 1 3.2 + 849.0 251.1 12.700 020 1 3.2 + 495.6 147.4 11.400 020 1 3.2 + 495.6 147.4 11.300 020 1 3.2 + 495.5 147.4 11.000 020 1 3.2 +1016.4 735.5 14.200 020 1 3.2 +1016.3 735.7 13.900 020 1 3.2 + 795.2 675.9 14.400 020 1 3.2 + 887.2 968.8 10.800 020 1 3.2 + 626.8 726.9 14.400 020 1 3.2 + 718.6 800.1 9.700 020 1 3.2 + 675.4 301.6 14.400 020 1 3.2 + 801.6 436.8 14.200 020 1 3.2 +1066.4 870.0 13.700 020 1 3.2 +1066.5 869.9 13.800 020 1 3.2 + 644.5 662.9 9.900 020 1 3.2 + 639.1 760.4 13.000 020 1 3.2 + 822.7 663.1 12.800 020 1 3.2 + 822.7 663.1 13.100 020 1 3.2 + 677.2 716.2 12.900 020 1 3.2 + 765.6 641.5 13.200 020 1 3.2 + 765.7 641.4 13.500 020 1 3.2 + 998.3 425.8 14.400 020 1 3.2 + 713.7 716.2 13.800 020 1 3.2 + 814.7 887.6 2.000 020 1 3.2 + 808.7 884.1 8.200 020 1 3.2 +1812.3 259.3 14.500 020 1 3.2 +1655.9 655.4 12.200 020 1 3.2 +1655.9 655.4 12.100 020 1 3.2 +1577.2 839.6 11.500 020 1 3.2 +1577.2 839.6 11.400 020 1 3.2 +1643.7 597.0 14.100 020 1 3.2 +1552.8 251.1 10.100 020 1 3.2 +1552.8 251.1 10.300 020 1 3.2 +1410.0 999.9 13.900 020 1 3.2 +1335.6 687.0 13.000 020 1 3.2 +1335.8 686.8 12.700 020 1 3.2 +1471.4 262.0 14.200 020 1 3.2 +1479.6 608.3 14.200 020 1 3.2 +1478.1 211.5 12.000 020 1 3.2 +1478.3 211.5 12.000 020 1 3.2 +1552.5 46.4 14.500 020 1 3.2 +2042.4 548.3 10.800 020 1 3.2 +2042.6 548.3 11.000 020 1 3.2 +1448.0 915.7 14.100 020 1 3.2 +1462.3 989.2 13.000 020 1 3.2 +1462.3 989.0 12.700 020 1 3.2 +1316.9 787.3 14.100 020 1 3.2 +1787.6 883.5 12.200 020 1 3.2 +1787.6 883.5 12.200 020 1 3.2 +1979.1 300.3 12.900 020 1 3.2 +1979.1 300.3 12.900 020 1 3.2 +1677.2 106.4 13.100 020 1 3.2 +1677.0 106.4 13.300 020 1 3.2 +1324.2 429.9 9.200 020 1 3.2 +1324.2 429.9 9.000 020 1 3.2 +1462.8 214.3 13.700 020 1 3.2 +1462.8 214.3 14.200 020 1 3.2 +1405.4 95.5 11.700 020 1 3.2 +1405.2 95.5 11.600 020 1 3.2 +1604.8 90.0 11.200 020 1 3.2 +1604.8 90.0 11.300 020 1 3.2 +1921.5 738.1 13.500 020 1 3.2 +1921.5 738.1 13.500 020 1 3.2 +1459.8 473.6 13.400 020 1 3.2 +1459.7 473.6 13.800 020 1 3.2 +1438.7 77.8 12.200 020 1 3.2 +1438.6 77.8 12.200 020 1 3.2 +1355.0 114.6 13.000 020 1 3.2 +1354.9 114.6 13.200 020 1 3.2 +1941.4 940.7 13.000 020 1 3.2 +1941.6 940.7 12.900 020 1 3.2 +1994.0 800.8 13.900 020 1 3.2 +1993.8 800.8 14.000 020 1 3.2 +2024.4 907.0 13.700 020 1 3.2 +2024.4 907.0 13.800 020 1 3.2 +2035.7 881.6 13.900 020 1 3.2 +2035.9 881.6 13.800 020 1 3.2 +1994.0 924.2 12.100 020 1 3.2 +1994.0 924.2 12.000 020 1 3.2 Index: /tags/addstar-1-2/Ohana/src/addstar/include/2mass.h =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/include/2mass.h (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/include/2mass.h (revision 4773) @@ -0,0 +1,13 @@ + +/* structure for data on a catalog region */ +typedef struct { + char filename[256]; + double RA[2]; + int Nrec; +} TM_Region; + +short TM_J, TM_H, TM_K; + +GSCRegion *get2mass_acc (GSCRegion *patch, int *nfiles, char *path, char *accel); +Stars *get2mass_2DR_data (GSCRegion *region, GSCRegion *patch, int *nstars); +Stars *get2mass_AS_data (GSCRegion *region, GSCRegion *patch, int *nstars); Index: /tags/addstar-1-2/Ohana/src/addstar/include/addstar.h =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/include/addstar.h (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/include/addstar.h (revision 4773) @@ -0,0 +1,216 @@ +# include +# include +# 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 = GetPhotcodeCodebyName (NAME); \ + if (!CODE) { \ + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", NAME); \ + exit (0); } + +# if (0) +/* this structure is used for loading the data from the FITS smp file */ +typedef struct { + float X; + float Y; + float M; + float dM; + float Mgal; + float Map; + float sky; + float fx; + float fy; + float df; + char dophot; + char dummy[3]; +} StarData; +# endif + +typedef struct { + double X; + double Y; + double R; + double D; + double M, dM; + char dophot; + double sky; + double fx, fy, df; + double Mgal, Map; + int found; + short int code; + e_time t; +} Stars; + +enum {M_IMAGE, M_REFLIST, M_REFCAT}; + +/* global variables set in parameter file */ +char ImageCat[256]; +char ImageTemplate[256]; +char CatTemplate[256]; +char GSCFILE[256]; +char CATDIR[256]; +char TWO_MASS_DIR_AS[256]; +char TWO_MASS_DIR_DR2[256]; +char GSCDIR[256]; +char CDROM[256]; +double DEFAULT_RADIUS, NSIGMA, SNLIMIT; +double ZeroPt; +double Latitude, SiderealTime; +int VERBOSE; +int FITS_INPUT; +int SKIP_MISSED; +int XOVERSCAN, YOVERSCAN; +int XMIN, XMAX, YMIN, YMAX; +int ACCEPT_ASTROM; + +PhotCode *thiscode; +GSCRegion UserPatch; +Coords *MOSAIC; + +char DateKeyword[64]; +char DateMode[64]; +char UTKeyword[64]; +char MJDKeyword[64]; +char JDKeyword[64]; +char ExptimeKeyword[64]; +char AirmassKeyword[64]; +char CCDNumKeyword[64]; +char STKeyword[64]; +char SubpixDatafile[256]; + +int CalReference; +int CalColor; +int CALIBRATE; +int DUMP_MATCHES; +int ONLY_IMAGES; +int EXISTING_REGIONS; +int FORCE_READ; +int SKYPROBE; +int ADDREFS; +int REPLACE; +int ONLY_MATCH; + +int MODE; +char *DUMP; +char *SELECT_2MASS_QUALITY; + +time_t TIMEREF; + +int main (int argc, char **argv); +int Shutdown (char *message, ...); +void TrapSignal (int sig); +int SetSignals (); +void help (); +int args (int argc, char **argv); +void SetProtect (int mode); + +int gcatalog (Catalog *catalog); +void gregion_star (Stars *star, GSCRegion *region); +GSCRegion *gregion_image (Image *image, int *Nregions); +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); +Stars *gstars (char *file, int *NSTARS, Image *image); +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); +short airmass (short secz_image, double ra, double dec, double st, double latitude); +void InitCalibration (); +void FindCalibration (Image *image); +double get_subpix (double x, double y); +void load_subpix (); +double scat_subpix (double x, double y); +void mkcatalog (GSCRegion *region, Catalog *catalog); +Stars *grefstars (char *file, int *Nstars); +Stars **find_subset (GSCRegion *region, Stars *stars, int Nstars, int *NSTARS); +int match_refstars (Stars *stars, int Nstars); +int parse_time (Header *header); + +void check_permissions (char *basefile); +void make_backup (char *filename); +void uppercase (char *string); + +void sort_lists (float *X, float *Y, int *S, int N); +void fsort (float *X, int N); +void fsort2 (float *X, float *Y, int N); +void sort_stars (Stars *stars, int N); +void sort_regions (GSCRegion *region, int N); + +char *sec_to_date (time_t second); +double sec_to_jd (time_t second); +int chk_time (char *line); +int dms_to_ddd (double *Value, char *string); +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string); +int str_to_dtime (char *line, double *second); +int str_to_radec (double *ra, double *dec, char *str1, char *str2); +int str_to_time (char *line, time_t *second); +time_t jd_to_sec (double jd); +time_t date_to_sec (char *date); +time_t short_date_to_sec (char *date); + + +/** + there is an inconsistency to be resolved: fixed structures (like Image) + need a fixed bit-length time (e_time), but these functions all use the + UNIX time_t types, which may be 32 or 64 bits, depending on the machine. + This can be resolved by using time_t with these functions, but casting + between e_time and time_t when necessary (ie, cannot return data to an + e_time pointer from one of these functions) +**/ + +Stars *get2mass (GSCRegion *patch, int *NSTARS, int mode); +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 FindDecBand (double dec, double *DEC0, double *DEC1); + +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); +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); +int dump_rawstars (Stars *stars, int Nstars); +Stars *rfits (FILE *f, int *nstars); +Stars *rtext (FILE *f, int *nstars); + Index: /tags/addstar-1-2/Ohana/src/addstar/src/ConfigInit.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/ConfigInit.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/ConfigInit.c (revision 4773) @@ -0,0 +1,83 @@ +# include "addstar.h" + +void ConfigInit (int *argc, char **argv) { + + char *config, *file; + char RadiusWord[80]; + char PhotCodeFile[256]; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (1); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + ScanConfig (config, "RADIUS", "%s", 0, RadiusWord); + if (!strcasecmp (RadiusWord, "header")) { + DEFAULT_RADIUS = 0; + } else { + DEFAULT_RADIUS = atof (RadiusWord); + } + + XMIN = XMAX = YMIN = YMAX = 0; + ScanConfig (config, "NSIGMA", "%lf", 0, &NSIGMA); + ScanConfig (config, "XOVERSCAN", "%d", 0, &XOVERSCAN); + ScanConfig (config, "YOVERSCAN", "%d", 0, &YOVERSCAN); + ScanConfig (config, "ADDSTAR_XMIN", "%d", 0, &XMIN); + ScanConfig (config, "ADDSTAR_XMAX", "%d", 0, &XMAX); + ScanConfig (config, "ADDSTAR_YMIN", "%d", 0, &YMIN); + ScanConfig (config, "ADDSTAR_YMAX", "%d", 0, &YMAX); + + ScanConfig (config, "2MASS_DIR_AS", "%s", 0, TWO_MASS_DIR_AS); + ScanConfig (config, "2MASS_DIR_DR2", "%s", 0, TWO_MASS_DIR_DR2); + ScanConfig (config, "GSCDIR", "%s", 0, GSCDIR); + ScanConfig (config, "USNO_CDROM", "%s", 0, CDROM); + + ScanConfig (config, "IMAGE_CATALOG", "%s", 0, ImageCat); + ScanConfig (config, "IMAGE_CATALOG_TEMPLATE", "%s", 0, ImageTemplate); + ScanConfig (config, "CATALOG_TEMPLATE", "%s", 0, CatTemplate); + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "MIN_SN_FSTAT", "%lf", 0, &SNLIMIT); + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + + /* used by parse_time to find time-related keywords */ + ScanConfig (config, "DATE-KEYWORD", "%s", 0, DateKeyword); + ScanConfig (config, "DATE-MODE", "%s", 0, DateMode); + ScanConfig (config, "UT-KEYWORD", "%s", 0, UTKeyword); + ScanConfig (config, "MJD-KEYWORD", "%s", 0, MJDKeyword); + ScanConfig (config, "JD-KEYWORD", "%s", 0, JDKeyword); + + ScanConfig (config, "EXPTIME-KEYWORD", "%s", 0, ExptimeKeyword); + ScanConfig (config, "AIRMASS-KEYWORD", "%s", 0, AirmassKeyword); + ScanConfig (config, "CCDNUM-KEYWORD", "%s", 0, CCDNumKeyword); + + ScanConfig (config, "ST-KEYWORD", "%s", 0, STKeyword); + ScanConfig (config, "OBSERVATORY-LATITUDE", "%lf", 0, &Latitude); + ScanConfig (config, "SUBPIX_DATAFILE", "%s", 0, SubpixDatafile); + + if (!LoadPhotcodes (PhotCodeFile)) { + fprintf (stderr, "error loading photcode file %s\n", PhotCodeFile); + exit (1); + } + + free (config); + free (file); + +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/ConvertStars.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/ConvertStars.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/ConvertStars.c (revision 4773) @@ -0,0 +1,58 @@ +# include "addstar.h" + +# define SWAP_BYTE(X) \ + tmp = byte[X+0]; byte[X+0] = byte[X+1]; byte[X+1] = tmp; +# define SWAP_WORD(X) \ + tmp = byte[X+0]; byte[X+0] = byte[X+3]; byte[X+3] = tmp; \ + tmp = byte[X+1]; byte[X+1] = byte[X+2]; byte[X+2] = tmp; +# define SWAP_DBLE(X) \ + tmp = byte[X+0]; byte[X+0] = byte[X+7]; byte[X+7] = tmp; \ + tmp = byte[X+1]; byte[X+1] = byte[X+6]; byte[X+6] = tmp; \ + tmp = byte[X+2]; byte[X+2] = byte[X+5]; byte[X+5] = tmp; \ + tmp = byte[X+3]; byte[X+3] = byte[X+4]; byte[X+4] = tmp; + +# ifdef linux +# define BYTE_SWAP +# endif + +# ifdef sid +# define BYTE_SWAP +# endif + +# ifdef dec +# define BYTE_SWAP +# endif + +# define STARS_SIZE 44 + +int ConvertStars (StarData *data, int size, int nitems) { + + int i; + unsigned char *byte, tmp; + +# ifdef BYTE_SWAP + + if (size != STARS_SIZE) { + fprintf (stderr, "mismatch in type sizes (StarData) %d vs %d\n", size, STARS_SIZE); + return (FALSE); + } + + byte = (char *) data; + for (i = 0; i < nitems; i++, byte += size) { + SWAP_WORD (0); /* R */ + SWAP_WORD (4); /* Y */ + SWAP_WORD (8); /* M */ + SWAP_WORD (12); /* dM */ + SWAP_WORD (16); /* Mgal */ + SWAP_WORD (20); /* Map */ + SWAP_WORD (24); /* sky */ + SWAP_WORD (28); /* fx */ + SWAP_WORD (32); /* fy */ + SWAP_WORD (36); /* df */ + } + return (TRUE); + +# else + return (TRUE); +# endif +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/RegionOps.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/RegionOps.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/RegionOps.c (revision 4773) @@ -0,0 +1,233 @@ +# include "addstar.h" + +# define NBANDS 24 + +double DecBands[] = { + -90.0, -82.5, -75.0, -67.5, -60.0, + -52.5, -45.0, -37.5, -30.0, -22.5, + -15.0, -07.5, +00.0, +07.5, +15.0, + +22.5, +30.0, +37.5, +45.0, +52.5, + +60.0, +67.5, +75.0, +82.5, +90.0 +}; + +char *DecNames[] = { + "s8230", "s7500", "s6730", "s6000", "s5230", + "s4500", "s3730", "s3000", "s2230", "s1500", + "s0730", "s0000", "n0000", "n0730", "n1500", + "n2230", "n3000", "n3730", "n4500", "n5230", + "n6000", "n6730", "n7500", "n8230", +}; + +GSCRegion *LoadRegions (int *nregions) { + + int Nregions, i, NDec, Nx, Ny; + char temp[50], file[256]; + double RA0, RA1, DEC0, DEC1, D0, D1; + FILE *f; + GSCRegion *regions; + Header header; + Header theader; + Matrix matrix; + FTable table; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC region file %s\n", GSCFILE); + exit (1); + } + + /* init & load in table data */ + table.header = &theader; + if (!fits_fread_header (f, &header)) Shutdown ("can't read primary header for GSC Region table"); + if (!fits_fread_matrix (f, &matrix, &header)) Shutdown ("can't read primary matrix for GSC Region table"); + if (!fits_fread_ftable (f, &table, "REGIONS")) Shutdown ("can't read GSC Region table"); + + fits_scan (table.header, "NAXIS1", "%d", 1, &Nx); + fits_scan (table.header, "NAXIS2", "%d", 1, &Ny); + + Nregions = 0; + ALLOCATE (regions, GSCRegion, Ny + 2); /* allow room for N & S poles */ + + /* pole region is artificial - not in table */ + sprintf (file, "%s/n8230/pole.cpt", CATDIR); + regions[Nregions].DEC[0] = 86.25; + regions[Nregions].DEC[1] = 93.75; + regions[Nregions].RA[0] = 0.0; + regions[Nregions].RA[1] = 360.0; + strcpy (regions[Nregions].filename, file); + Nregions ++; + + /* fill out GSCRegion structure */ + for (i = 0; i < Ny; i++) { + strncpy (temp, &table.buffer[i*48], 48); + temp[49] = 0; + hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); + if (RA1 < RA0) RA1 += 360.0; + + if (DEC0 < DEC1) { + regions[Nregions].DEC[0] = DEC0; + regions[Nregions].DEC[1] = DEC1; + } else { + regions[Nregions].DEC[0] = DEC1; + regions[Nregions].DEC[1] = DEC0; + } + /* skip the pole entries */ + if (regions[Nregions].DEC[1] > 86.25) continue; + + regions[Nregions].RA[0] = RA0; + regions[Nregions].RA[1] = RA1; + + /** convert DEC0, DEC1 to NBigDec **/ + NDec = FindDecBand (0.5*(DEC0 + DEC1), &D0, &D1); + if (NDec == -1) Shutdown ("programming / table error"); + + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt", CATDIR, DecNames[NDec], &temp[1]); + + strcpy (regions[Nregions].filename, file); + Nregions ++; + } + if (Nregions > Ny + 2) { + fprintf (stderr, "ERROR: too many regions\n"); + exit (2); + } + + sort_regions (regions, Nregions); + *nregions = Nregions; + return (regions); +} + +/* given dec value, find upper and lower bounds of enclosing dec band */ +int FindDecBand (double dec, double *DEC0, double *DEC1) { + + int i; + + for (i = 0; i < NBANDS; i++) { + if (DecBands[i + 0] > dec) continue; + if (DecBands[i + 1] <= dec) continue; + *DEC0 = DecBands[i+0]; + *DEC1 = DecBands[i+1]; + return (i); + } + return (-1); +} + +/* note that region.DEC[0] is sorted, but not necessarily DEC[1] */ +int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi; + double D0, D1; + + if (dec <= 90) return (0); + + /* D0, D1 are upper lower band boundaries */ + /* all regions in this band have DEC[0] >= D0 */ + N = FindDecBand (dec, &D0, &D1); + if (N == -1) Shutdown ("programming / table error"); + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[0] < D0) { + Nlo = N; + } else { + Nhi = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] < v <= v[Nhi] */ + + /* find the last entry < dec */ + for (N = Nhi; N >= Nlo; N--) { + if (region[N].DEC[0] < D0) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +/* note that region.DEC[0] is sorted, but not necessarily DEC[1] */ +int FindRegionDecBandStop (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi; + double D0, D1; + + if (dec >= 90) return (Nregion - 1); + + /* D0, D1 are upper lower band boundaries */ + /* all regions in this band have DEC[0] < D1 */ + N = FindDecBand (dec, &D0, &D1); + if (N == -1) Shutdown ("programming / table error"); + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[0] > D1) { + Nhi = N; + } else { + Nlo = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] <= v < v[Nhi] */ + + /* find the last entry < dec */ + for (N = Nlo; N < Nhi; N++) { + if (region[N].DEC[0] > D1) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +int FindRegionByPoint (GSCRegion *region, double ra, double dec, GSCRegion *gsc, int Ngsc) { + + int i, Ns; + + Ns = FindRegionDecBandStart (gsc, Ngsc, dec); + if (Ns < 0) return (FALSE); + + for (i = Ns; i < Ngsc; i++) { + if (gsc[i].RA[0] >= ra) continue; + if (gsc[i].RA[1] < ra) continue; + if (gsc[i].DEC[0] >= dec) continue; + if (gsc[i].DEC[1] < dec) continue; + region[0] = gsc[i]; + return (TRUE); + } + return (FALSE); +} + +/* bracket search: + + Nlo T N Nhi + X X X X X X X V X X X X X X X X X X X X X X + Nhi + + Nlo N T Nhi + X X X X X X X X X X X X X X X V X X X X X X + Nlo + +*/ + +/* + if (buffer[i*48 + 19] == ' ') continue; + strncpy (temp, &buffer[i*48 + 19], 20); temp[20] = 0; + hms_to_deg (&RA, &DEC, temp, ' ', 3); + + + dBigRA = 360.0 / (int)(0.5 + 48*cos(3.1415927*0.5*(BigDecBounds[NBigDec] + BigDecBounds[NBigDec + 1])/180.0)); + + NBigRA = ra / dBigRA; + + NBig = NBigRA; + for (i = 0; (i < 12) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + for (i = 13; (i < 24) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + + + fprintf (stderr, "%d %d %d %d %f %f -> %f %f %f\n", NBigDec, NBigRA, NBigRASections[NBigDec], NBig, ra, dec, BigDecBounds[NBigDec], BigDecBounds[NBigDec + 1], dBigRA); +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/SetSignals.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/SetSignals.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/SetSignals.c (revision 4773) @@ -0,0 +1,118 @@ +# include "addstar.h" + +static int Protect = FALSE; +static int Trapped = FALSE; + +void TrapSignal (int sig) { + fprintf (stderr, "trapped signal %d\n", sig); + if (sig == 11) { + fprintf (stderr, "seg fault\n"); + exit (1); + } + if (Protect) { + Trapped = TRUE; + fprintf (stderr, "blocking until protected sections are clear\n"); + return; + } + Shutdown ("halted by signal (trapped)"); +} + +void SetProtect (int mode) { + Protect = mode; + if (Trapped && !Protect) Shutdown ("halted by signal (protect)"); +} + +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: /tags/addstar-1-2/Ohana/src/addstar/src/addstar.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/addstar.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/addstar.c (revision 4773) @@ -0,0 +1,90 @@ +# include "addstar.h" + +int main (int argc, char **argv) { + + int i, Nstars, Noverlap, Nregions, Nsubset, Naverage, Nmeasure; + Stars *stars, **subset; + Image image, *overlap; + GSCRegion *regions; + Catalog catalog; + + double dtime; + struct timeval start, stop; + + gettimeofday (&start, NULL); + + SetSignals (); + ConfigInit (&argc, argv); + args (argc, argv); + + stars = NULL; + overlap = NULL; + regions = NULL; + + if (SKYPROBE) load_subpix (); + + lock_image_db (); + + switch (MODE) { + case M_IMAGE: + stars = gstars (argv[1], &Nstars, &image); + if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars (stars, Nstars); + regions = gregion_image (&image, &Nregions); + overlap = gimages (&image, &Noverlap); + break; + case M_REFLIST: + stars = grefstars (argv[1], &Nstars); + regions = gregion_stars (stars, Nstars, &Nregions); + break; + case M_REFCAT: + regions = gregion_patch (&UserPatch, &Nregions); + break; + } + if (ONLY_MATCH || EXISTING_REGIONS) { + regions = gregion_match (regions, &Nregions); + } + if (VERBOSE) fprintf (stderr, "writing to %d regions\n", Nregions); + + Naverage = Nmeasure = 0; + for (i = 0; i < Nregions; i++) { + if (!load_pt_catalog (&catalog, ®ions[i])) continue; + + Naverage += catalog.Naverage; + Nmeasure += catalog.Nmeasure; + + switch (MODE) { + case M_IMAGE: + Nsubset = Nstars; + find_matches (®ions[i], stars, Nstars, &catalog, &image, overlap, Noverlap); + break; + case M_REFCAT: + stars = grefcat (argv[1], ®ions[i], &Nstars); + case M_REFLIST: + subset = find_subset (®ions[i], stars, Nstars, &Nsubset); + find_matches_refstars (®ions[i], subset, Nsubset, &catalog); + if (Nsubset) free (subset); + break; + } + if (Nsubset == 0) { + unlock_catalog (&catalog); + } else { + save_pt_catalog (&catalog); + } + } + + if (CALIBRATE) { FindCalibration (&image); } + unlock_image_db (&image); + + gettimeofday (&stop, NULL); + dtime = DTIME (stop, start); + fprintf (stderr, "SUCCESS: elapsed time %8.3f sec for %5d stars, %6d average, %7d measure\n", dtime, Nstars, Naverage, Nmeasure); + exit (0); +} + +/* names: + + catalog - existing object db table + regions - sky area which may or may contain data + patch - RA,DEC bounded portion of sky + +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/airmass.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/airmass.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/airmass.c (revision 4773) @@ -0,0 +1,19 @@ +# include "addstar.h" + +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... ***/ + /* ra, dec, latitude in dec deg; st in dec hours */ + + /* hour : hour angle in degrees */ + hour = 15.0*st - ra; + cosz = sin (RAD_DEG*dec) * sin (RAD_DEG*latitude) + cos (RAD_DEG*dec) * cos (RAD_DEG*hour) * cos (RAD_DEG*latitude); + + /* modify to match storage units (millimags) */ + secz = 1000.0 / cosz; + return (secz); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/aregion.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/aregion.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/aregion.c (revision 4773) @@ -0,0 +1,171 @@ +# include "addstar.h" + +double BigDecBounds[] = {0.0, 7.5, 15.0, 22.5, 30.0, 37.5, 45.0, + 52.5, 60.0, 67.5, 75.0, 82.5, 90.0, + 0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, + -52.5, -60.0, -67.5, -75.0, -82.5, -90.0}; +char *DecSections[] = {"N0000", "N0730", "N1500", "N2230", "N3000", "N3730", "N4500", + "N5230", "N6000", "N6730", "N7500", "N8230", "weirdness", + "S0000", "S0730", "S1500", "S2230", "S3000", "S3730", "S4500", + "S5230", "S6000", "S6730", "S7500", "S8230", "weirdness"}; + +char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", + "n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", + "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", + "s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"}; + +char *disk[] = {"disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", + "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "weirdness", + "disk 1", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", + "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "weirdness"}; + +int NBigRASections [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3, 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3}; + +int NDecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, + 0, 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48}; + +# define DEBUG 1 + +/* find region file which contains ra, dec */ +void aregion (GSCRegion *region, FILE *f, double ra, double dec) { + + char buffer[28800], temp[50], file[256]; + double RA0, RA1, DEC0, DEC1; + int i, NBigDec, NLINES, done; + + while (ra < 0) { ra += 360.0; } + while (ra >= 360.0) { ra -= 360.0; } + + if (dec >= 86.25) { + sprintf (file, "%s/n8230/pole.cpt", CATDIR); + region[0].DEC[0] = 86.25; + region[0].DEC[1] = 93.75; + region[0].RA[0] = 0.0; + region[0].RA[1] = 360.0; + strcpy (region[0].filename, file); + return; + } + + NBigDec = -1; + for (i = 0; i < 12; i++) { + if ((dec >= BigDecBounds[i]) && (dec < BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + if (NBigDec < 0) { + for (i = 13; i < 24; i++) { + if ((dec < BigDecBounds[i]) && (dec >= BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + } + if (NBigDec < 0) { + fprintf (stderr, "dec out of range: %f\n", dec); + } + + NLINES = 0; + for (i = 0; i < NBigDec; i++) { + NLINES += NDecLines[i]; + } + fseek (f, 5*2880 + 48*NLINES, SEEK_SET); + + done = FALSE; + Fread (buffer, 1, 48*NDecLines[NBigDec], f, "char"); + for (i = 0; !done && (i < NDecLines[NBigDec]); i++) { + strncpy (temp, &buffer[i*48], 48); + temp[49] = 0; + hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); + if (RA1 < RA0) RA1 += 360.0; + if ((dec >= 0) && (dec >= DEC0) && (dec < DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + if ((dec < 0) && (dec < DEC0) && (dec >= DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + } + + if (!done) { + fprintf (stderr, "ERROR in search: %f %f\n", ra, dec); + exit (1); + } + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt", CATDIR, Dec2Sections[NBigDec], &temp[1]); + if (DEC0 < DEC1) { + region[0].DEC[0] = DEC0; + region[0].DEC[1] = DEC1; + } else { + region[0].DEC[0] = DEC1; + region[0].DEC[1] = DEC0; + } + region[0].RA[0] = RA0; + region[0].RA[1] = RA1; + strcpy (region[0].filename, file); + return; +} + + +/**********/ +int hms_to_deg (double *h0, double *h1, double *d0, double *d1, char *string) { + + int flag_d0, flag_d1, flag_h0, flag_h1; + double tmp; + + *d0 = *h0 = *d1 = *h1 = 0; + + flag_h0 = dparse (h0, 1, string); + flag_h1 = dparse (h1, 4, string); + flag_d0 = dparse (d0, 7, string); + flag_d1 = dparse (d1, 9, string); + *h0 *= flag_h0; + *h1 *= flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + dparse (&tmp, 2, string); + *h0 += tmp/60.0; + dparse (&tmp, 3, string); + *h0 += tmp/3600.0; + + dparse (&tmp, 5, string); + *h1 += tmp/60.0; + dparse (&tmp, 6, string); + *h1 += tmp/3600.0; + + dparse (&tmp, 8, string); + *d0 += tmp/60.0; + + dparse (&tmp, 10, string); + *d1 += tmp/60.0; + + *h0 *= 15*flag_h0; + *h1 *= 15*flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + return (TRUE); +} + + +/* + if (buffer[i*48 + 19] == ' ') continue; + strncpy (temp, &buffer[i*48 + 19], 20); temp[20] = 0; + hms_to_deg (&RA, &DEC, temp, ' ', 3); + + + dBigRA = 360.0 / (int)(0.5 + 48*cos(3.1415927*0.5*(BigDecBounds[NBigDec] + BigDecBounds[NBigDec + 1])/180.0)); + + NBigRA = ra / dBigRA; + + NBig = NBigRA; + for (i = 0; (i < 12) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + for (i = 13; (i < 24) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + + + fprintf (stderr, "%d %d %d %d %f %f -> %f %f %f\n", NBigDec, NBigRA, NBigRASections[NBigDec], NBig, ra, dec, BigDecBounds[NBigDec], BigDecBounds[NBigDec + 1], dBigRA); +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/args.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/args.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/args.c (revision 4773) @@ -0,0 +1,223 @@ +# include "addstar.h" +# define NARGS 2 /* minimum is: addstar (filename) */ + +void help () { + + fprintf (stderr, "USAGE\n"); + fprintf (stderr, " addstar (filename)\n"); + fprintf (stderr, " add specified image (cmp format) to database\n\n"); + fprintf (stderr, " addstar -ref (filename)"); + fprintf (stderr, " add ASCII data (ra dec mag dmag) to database\n\n"); + fprintf (stderr, " addstar -cat (catalog)"); + fprintf (stderr, " add data from catalog (USNO/2MASS/GSC) to database\n\n"); + + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -region ra ra dec dec : only add data in specified region (-ref mode only)\n"); + fprintf (stderr, " -p (photcode) : specify photcode (override header)\n"); + fprintf (stderr, " -time (YYYY/MM/DD,HH:MM:SS) : specify date/time (override header)\n"); + fprintf (stderr, " -mosaic (filename) : identify associated mosaic frame for chip image\n"); + fprintf (stderr, " -fits : input file is FITS table, not TEXT table\n"); + fprintf (stderr, " -existing-regions : only add measurements to existing catalog files\n"); + fprintf (stderr, " -only-match : only add measurements to existing objects\n"); + fprintf (stderr, " -missed : skipped 'missed' entries\n"); + fprintf (stderr, " -replace : replace time/photcode measurements (no duplication)\n"); + fprintf (stderr, " -image : only insert image data\n"); + fprintf (stderr, " -cal : perform zero-point calibration\n"); + fprintf (stderr, " -skyprobe : specify skyprobe mode\n"); + fprintf (stderr, " -accept : accept bad astrometry from header\n"); + fprintf (stderr, " -force : force read of database with inconsistent info\n"); + fprintf (stderr, " -v : verbose mode\n"); + fprintf (stderr, " -dump (mode) : output test data\n"); + fprintf (stderr, " -help : this list\n"); + fprintf (stderr, " -h : this list\n\n"); + exit (2); +} + +int args (int argc, char **argv) { + + int N; + + /* check for help request */ + if (get_argument (argc, argv, "-help") || + get_argument (argc, argv, "-h")) { + help (); + } + + /*** check for command line options ***/ + + /* basic mode: image, list, refcat */ + MODE = M_IMAGE; + if ((N = get_argument (argc, argv, "-ref"))) { + MODE = M_REFLIST; + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-cat"))) { + MODE = M_REFCAT; + remove_argument (N, &argc, argv); + } + + /*** provide additional data ***/ + /* restrict to a portion of the sky? (REFCAT only) */ + UserPatch.RA[0] = 0; + UserPatch.RA[1] = 360; + UserPatch.DEC[0] = -90; + UserPatch.DEC[1] = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.RA[0] = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.RA[1] = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.DEC[0] = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.DEC[1] = atof (argv[N]); + remove_argument (N, &argc, argv); + } + /* 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 ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!str_to_time (argv[N], &TIMEREF)) { + fprintf (stderr, "syntax error in time\n"); + exit (1); + } + remove_argument (N, &argc, argv); + } + /* provide a mosaic for distortion */ + MOSAIC = NULL; + if ((N = get_argument (argc, argv, "-mosaic"))) { + Header header; + ALLOCATE (MOSAIC, Coords, 1); + + remove_argument (N, &argc, argv); + if (!fits_read_header (argv[N], &header)) { + fprintf (stderr, "ERROR: can't read header for mosaic %s\n", argv[N]); + exit (1); + } + if (!GetCoords (MOSAIC, &header)) { + fprintf (stderr, "ERROR: no astrometric solution in header\n"); + exit (1); + } + if (strcmp(&MOSAIC[0].ctype[4], "-DIS")) { + fprintf (stderr, "ERROR: not a mosaic distortion header\n"); + exit (1); + } + remove_argument (N, &argc, argv); + fits_free_header (&header); + } + FITS_INPUT = FALSE; + if ((N = get_argument (argc, argv, "-fits"))) { + FITS_INPUT = TRUE; + remove_argument (N, &argc, argv); + } + + /*** modify behavior ***/ + /* only add to existing objects */ + EXISTING_REGIONS = FALSE; + if ((N = get_argument (argc, argv, "-existing-regions"))) { + EXISTING_REGIONS = TRUE; + remove_argument (N, &argc, argv); + } + /* only add to existing objects */ + ONLY_MATCH = FALSE; + if ((N = get_argument (argc, argv, "-only-match"))) { + ONLY_MATCH = TRUE; + remove_argument (N, &argc, argv); + } + /* don't add missed pts to Missed table (image only) */ + SKIP_MISSED = FALSE; + if ((N = get_argument (argc, argv, "-missed"))) { + SKIP_MISSED = TRUE; + remove_argument (N, &argc, argv); + } + /* replace measurement, don't duplicate (ref/cat only) */ + REPLACE = FALSE; + if ((N = get_argument (argc, argv, "-replace"))) { + REPLACE = TRUE; + remove_argument (N, &argc, argv); + } + /* only add image potion to image table */ + ONLY_IMAGES = FALSE; + if ((N = get_argument (argc, argv, "-image"))) { + ONLY_IMAGES = TRUE; + remove_argument (N, &argc, argv); + } + /* apply average zpt offset calibration (image only) */ + CALIBRATE = FALSE; + if ((N = get_argument (argc, argv, "-cal"))) { + CALIBRATE = TRUE; + remove_argument (N, &argc, argv); + } + /* define 2MASS quality flags to keep */ + SELECT_2MASS_QUALITY = NULL; + if ((N = get_argument (argc, argv, "-2massquality"))) { + remove_argument (N, &argc, argv); + SELECT_2MASS_QUALITY = strcreate (argv[N]); + remove_argument (N, &argc, argv); + } + + /*** optional situations ***/ + /* treat data specially for skyprobe (calibration, subpix) */ + SKYPROBE = FALSE; + if ((N = get_argument (argc, argv, "-skyprobe"))) { + SKYPROBE = TRUE; + remove_argument (N, &argc, argv); + } + /* accept bad header astrometry */ + ACCEPT_ASTROM = FALSE; + if ((N = get_argument (argc, argv, "-accept"))) { + ACCEPT_ASTROM = 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"))) { + FORCE_READ = TRUE; + remove_argument (N, &argc, argv); + } + /* extra error messages */ + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + DUMP = NULL; + if ((N = get_argument (argc, argv, "-dump"))) { + remove_argument (N, &argc, argv); + DUMP = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + + if (argc != 2) { + fprintf (stderr, "USAGE: addstar (filename)\n"); + fprintf (stderr, "USAGE: addstar -cat (catalog)\n"); + fprintf (stderr, "USAGE: addstar -ref (filename)\n"); + exit (2); + } + return (TRUE); +} + +/** addstar modes: + + addstar (image.smp) - add cmp/smp image data to db + addstar -ref (file.dat) (photcode) + addstar -cat (USNO/2MASS/GSC) -region (ra dec - ra dec) + + -replace : ref/cat - replace existing match (photcode/time) + -match : ref/cat - only add measures to existing averages + + ref types: + ASCII - RA,DEC,M,dM in a table + + addstar + +**/ + Index: /tags/addstar-1-2/Ohana/src/addstar/src/bracket.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/bracket.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/bracket.c (revision 4773) @@ -0,0 +1,37 @@ +int FindRegionDecBandStart (GSCRegion *region, int Nregion, double dec) { + + int N, Nlo, Nhi, N1, N2; + + /* bracket value of interest */ + Nlo = 0; Nhi = Nregion; + while (Nhi - Nlo > 10) { + N = 0.5*(Nlo + Nhi); + if (region[N].DEC[0] < dec) { + Nlo = N; + } else { + Nhi = N; + } + } + /* Nlo - Nhi bracket the transition across dec */ + /* v[Nlo] < v <= v[Nhi] */ + + /* find the first entry >= dec */ + for (N = N1; N < N2; N++) { + if (region[N].DEC[0] >= T) return (N); + } + return (-1); + /* this will fail if transition value is not in the data range */ +} + +/* bracket search: + + Nlo T N Nhi + X X X X X X X V X X X X X X X X X X X X X X + Nhi + + Nlo N T Nhi + X X X X X X X X X X X X X X X V X X X X X X + Nlo + +*/ + Index: /tags/addstar-1-2/Ohana/src/addstar/src/calibrate.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/calibrate.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/calibrate.c (revision 4773) @@ -0,0 +1,179 @@ +# include "addstar.h" + +static int Ncal, NCAL, *Nstar; +static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm; + +void InitCalibration () { + + fprintf (stderr, "calibrating the image... %d\n", thiscode[0].equiv); + Ncal = 0; + NCAL = 1000; + ALLOCATE (Mobs, float, NCAL); + ALLOCATE (dMobs, float, NCAL); + ALLOCATE (Mref, float, NCAL); + ALLOCATE (dMref, float, NCAL); + ALLOCATE (Cref, float, NCAL); + ALLOCATE (Airm, float, NCAL); + ALLOCATE (Nstar, int, NCAL); +} + +void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N) { + + Mobs[Ncal] = M; + dMobs[Ncal] = dM; + Mref[Ncal] = Mr; + dMref[Ncal] = dMr; + Cref[Ncal] = Mc; + Airm[Ncal] = A; + Nstar[Ncal] = N; + Ncal ++; + + if (Ncal == NCAL) { + NCAL += 1000; + REALLOCATE (dMobs, float, NCAL); + REALLOCATE (Mobs, float, NCAL); + REALLOCATE (Mref, float, NCAL); + REALLOCATE (dMref, float, NCAL); + REALLOCATE (Cref, float, NCAL); + REALLOCATE (Airm, float, NCAL); + REALLOCATE (Nstar, int, 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; +} + +void FindCalibration (Image *image) { + + 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; + + /* reject multiple matched-stars */ + /* find maximum value of Nstar[] */ + MaxN = -1; + for (i = 0; i < Ncal; i++) { + MaxN = MAX (Nstar[i], MaxN); + } + if (MaxN == -1) { + fprintf (stderr, "no clean stars\n"); + image[0].Mcal = 10000; + image[0].dMcal = 10000; + return; + } + /* create a hash array from Nstar[] entries */ + ALLOCATE (Nlist, int, MaxN + 1); + bzero (Nlist, MAX (0, MaxN*sizeof(int))); + for (i = 0; i < Ncal; i++) { + Nlist[Nstar[i]] ++; + } + + /* accumulate delta mags */ + ALLOCATE (Dmag, float, Ncal); + ALLOCATE (dDmag, float, Ncal); + Nkeep = 0; + Klam = thiscode[0].K; + Clam = thiscode[0].C; + Xlam = thiscode[0].X[0]; + for (i = 0; i < Ncal; i++) { + /* if this entry has too many (or two few?) matches, skip it */ + if (Nlist[Nstar[i]] != 1) continue; + Mint = Clam + Mobs[i] + Xlam*Cref[i] + Klam*(Airm[i] - 1000); + Mabs = 0.001*Mint - ZeroPt; + /* note: subpix correction is applied in gstars */ + + /* skip stars brighter than 8.0 */ + if (SKYPROBE) { + if (Mabs > 9.0) continue; + if (Mabs < 5.0) continue; + } + + dMr = MAX (0.005, 0.001*dMref[i]); + dMo = MAX (0.005, 0.001*dMobs[i]); + + Dmag[Nkeep] = (Mabs - 0.001*Mref[i]); + dDmag[Nkeep] = (dMr*dMr + dMo*dMo); + Nkeep ++; + } + + if (Nkeep < 5) { + fprintf (stderr, "too few stars\n"); + image[0].Mcal = 10000; + image[0].dMcal = 10000; + return; + } + fsort2 (Dmag, dDmag, Nkeep); + + /* take sort list of Dmag, find median */ + Dmed = Dmag[(int)(0.5*Nkeep)]; + + /* exclude points with abs(Dmag - Dmed) / dDmag > 2.5 */ + + /* accumulate delta mags (25% - 75% of clipped range) */ + W1 = 0.0; + W2 = 0.0; + M1 = 0.0; + M2 = 0.0; + N = 0.0; + for (i = 0; i < Nkeep; i++) { + NSigma = fabs (Dmag[i] - Dmed) / sqrt (dDmag[i]); + if (NSigma > 2.5) continue; + W1 += Dmag[i] / dDmag[i]; + W2 += 1 / dDmag[i]; + M1 += Dmag[i]; + M2 += SQ (Dmag[i]); + N += 1.0; + } + + if (N > 1) { + M1 = M1 / N; + M2 = sqrt (fabs(M2/N - M1*M1)); + Mw = W1 / W2; + fprintf (stderr, "N: %.0f, mean: %f, wt mean: %f, stdev: %f, precision: %f\n", N, M1, Mw, M2, M2 / sqrt (N)); + image[0].Mcal = 1000 * M1; + image[0].dMcal = 1000 * M2 / sqrt (N); + image[0].Mxxxx = N; + } else { + fprintf (stderr, "too few stars\n"); + image[0].Mcal = 10000; + image[0].dMcal = 10000; + image[0].Mxxxx = 0; + } +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/check_permissions.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/check_permissions.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/check_permissions.c (revision 4773) @@ -0,0 +1,67 @@ +# include "addstar.h" + +void check_permissions (char *basefile) { + + char *c, dir[256], filename[256]; + struct stat filestat; + uid_t uid; + gid_t gid; + int status, cmode; + + uid = getuid(); + gid = getgid(); + + /* check permission to write to directory */ + sprintf (filename, "%s", basefile); + c = strrchr (filename, '/'); + if (c == (char *) NULL) { + strcpy (dir, "."); + } else { + *c = 0; + strcpy (dir, filename); + } + status = stat (dir, &filestat); + if (status == -1) { + fprintf (stderr, "directory %s does not exist, creating...\n", dir); + cmode = S_IRWXU | S_IRWXG | S_IRWXO; + status = mkdir (dir, cmode); + if (status == -1) { + fprintf (stderr, "ERROR: can't create %s\n", dir); + exit (1); + } + } + status = stat (dir, &filestat); + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRWXU)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRWXG)) || + (filestat.st_mode & S_IRWXO)) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", dir); + exit (1); + } + + /* check permission to write to file */ + sprintf (filename, "%s", basefile); + status = stat (filename, &filestat); + if (status == 0) { /* file exists, are permissions OK? */ + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || + ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", filename); + exit (1); + } + } + + /* check permission to write to backup file */ + sprintf (filename, "%s~", basefile); + status = stat (filename, &filestat); + if (status == 0) { /* file exists, are permissions OK? */ + if (((uid == filestat.st_uid) && (filestat.st_mode & S_IRUSR) && (filestat.st_mode & S_IWUSR)) || + ((gid == filestat.st_gid) && (filestat.st_mode & S_IRGRP) && (filestat.st_mode & S_IWGRP)) || + ((filestat.st_mode & S_IROTH) && (filestat.st_mode & S_IWOTH))) { + } else { + fprintf (stderr, "ERROR: can't write to %s\n", filename); + exit (1); + } + } +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/conversions.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/conversions.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/conversions.c (revision 4773) @@ -0,0 +1,364 @@ +# include "addstar.h" +enum {TIME_NONE, TIME_DATE, TIME_HOURS, TIME_DAYS, TIME_JD, TIME_MJD, TIME_SEC}; + +/***** convert [-]00:00:00 to 0.0000 ****/ +int dms_to_ddd (double *Value, char *string) { + + int valid, neg, status; + double tmp, value; + char *p1, *p2, *px; + + valid = FALSE; + neg = FALSE; + stripwhite (string); + p1 = string; + px = string + strlen(string); + + if (string[0] == '-') { + valid = TRUE; + neg = TRUE; + p1 = &string[1]; + } + if (string[0] == '+') { + valid = TRUE; + neg = FALSE; + p1 = &string[1]; + } + if (isdigit(string[0])) { + valid = TRUE; + p1 = &string[0]; + } + if (!valid) { return (FALSE); } + + status = 1; + tmp = strtod (p1, &p2); + if (p2 == p1) return (FALSE); /* entry not a number: +fred */ + value = tmp; + if (p2 == px) goto escape; /* entry only number: +1.0 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:fred */ + status = 2; + value += tmp / 60.0; + if (p2 == px) goto escape; /* entry only number: +1:1 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */ + value += tmp / 3600.0; + + escape: + if (neg) { + value *= -1; + } + *Value = value; + + return (status); +} + +/**********/ +int str_to_radec (double *ra, double *dec, char *str1, char *str2) { + + double Ra, Dec; + + *ra = *dec = 0; + switch (dms_to_ddd (&Ra, str1)) { + case 0: + fprintf (stderr, "syntax error in RA\n"); + return (FALSE); + case 1: + break; + case 2: + Ra = Ra * 15; + break; + } + switch (dms_to_ddd (&Dec, str2)) { + case 0: + fprintf (stderr, "syntax error in DEC\n"); + return (FALSE); + case 1: + case 2: + break; + } + *ra = Ra; + *dec = Dec; + return (TRUE); +} + +# define FORMAT_DAYS 1 +# define FORMAT_HOURS 2 +# define FORMAT_MINUTES 3 +# define FORMAT_SECONDS 4 +# define FORMAT_JD 5 +# define FORMAT_DATE 6 + +/**********/ +int chk_time (char *line) { + + char *p1, *p2; + double tmp; + int mode; + + p1 = line; + tmp = strtod (p1, &p2); + mode = FORMAT_DATE; + if (p2 == p1 + strlen (p1) - 1) { + if (*p2 == 'd') { + mode = FORMAT_DAYS; + } + if (*p2 == 'h') { + mode = FORMAT_HOURS; + } + if (*p2 == 'm') { + mode = FORMAT_MINUTES; + } + if (*p2 == 's') { + mode = FORMAT_SECONDS; + } + if (*p2 == 'j') { + mode = FORMAT_JD; + } + } + return (mode); +} + +/**********/ +int str_to_time (char *line, time_t *second) { + + char *tmpline; + struct tm *gmt; + struct timeval now; + double jd; + time_t tsec; + + if (!strcasecmp (line, "NOW")) { + gettimeofday (&now, (struct timezone *) NULL); + *second = now.tv_sec; + return (TRUE); + } + + if (!strncasecmp (line, "TODAY", 5)) { + gettimeofday (&now, (struct timezone *) NULL); + if (line[5]) { /* line has extra data (ie, hh:mm:ss) */ + tsec = now.tv_sec; + ALLOCATE (tmpline, char, 64); + gmt = gmtime (&tsec); + sprintf (tmpline, "%04d/%02d/%02d,%s", 1900 + gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, &line[6]); + *second = date_to_sec (tmpline); + free (tmpline); + return (TRUE); + } else { + *second = 86400 * ((int)(now.tv_sec / 86400)); + return (TRUE); + } + } + + switch (chk_time (line)) { + case 0: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + case FORMAT_JD: + jd = strtod (line, 0); + *second = jd_to_sec (jd); + return (TRUE); + case FORMAT_DATE: + *second = date_to_sec (line); + return (TRUE); + } + return (FALSE); +} + + +/**********/ +int str_to_dtime (char *line, double *second) { + + switch (chk_time (line)) { + case 0: + case FORMAT_JD: + case FORMAT_DATE: + return (FALSE); + case FORMAT_DAYS: + *second = strtod (line, 0) * 86400.0; + return (TRUE); + case FORMAT_HOURS: + *second = strtod (line, 0) * 3600.0; + return (TRUE); + case FORMAT_MINUTES: + *second = strtod (line, 0) * 60.0; + return (TRUE); + case FORMAT_SECONDS: + *second = strtod (line, 0); + return (TRUE); + } + return (FALSE); +} + +/**********/ +double sec_to_jd (time_t second) { + + double jd; + + jd = second/86400.0 + 2440587.5; + return (jd); +} + +/**********/ +time_t jd_to_sec (double jd) { + + time_t second; + + second = (jd - 2440587.5)*86400; + return (second); +} + +/**********/ +char *sec_to_date (time_t second) { + + struct tm *gmt; + char *line; + + ALLOCATE (line, char, 64); + gmt = gmtime (&second); + sprintf (line, "%04d/%02d/%02d,%02d:%02d:%02d", + 1900 + gmt[0].tm_year, gmt[0].tm_mon+1, gmt[0].tm_mday, + gmt[0].tm_hour, gmt[0].tm_min, gmt[0].tm_sec); + return (line); + +} + +/***** date in format yyyy/mm/dd,hh:mm:ss *****/ +time_t date_to_sec (char *date) { + + time_t second; + double tmp, jd; + struct tm now; + char *p1, *p2, *px; + + p1 = date; + px = date + strlen(date); + bzero (&now, sizeof(now)); + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_year = tmp; + if (now.tm_year > 1000) now.tm_year -= 1900; + if (now.tm_year < 50) now.tm_year += 100; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mon = tmp - 1; /* mon runs from 0 - 11 */ + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_mday = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + p1 = p2 + 1; + now.tm_hour = tmp; + if (p2 == px) goto escape; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_min = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; + now.tm_sec = tmp; + if (p2 == px) goto escape; + p1 = p2 + 1; + + escape: + jd = now.tm_mday - 32075 + (int)(1461*(1900 + now.tm_year + 4800 + (int)(((now.tm_mon+1)-14)/12))/4) + + (int)(367*((now.tm_mon+1) - 2 - (int)(((now.tm_mon+1) - 14)/12)*12)/12) + - (int)(3*(int)((1900 + now.tm_year + 4900 + (int)(((now.tm_mon+1) - 14)/12))/100)/4) - 0.5; + + second = (jd - 2440587.5)*86400 + 3600.0*now.tm_hour + now.tm_min*60.0 + now.tm_sec; + + return (second); +} + +/***** 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)); + + sscanf (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: /tags/addstar-1-2/Ohana/src/addstar/src/dump.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/dump.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/dump.c (revision 4773) @@ -0,0 +1,21 @@ +# include "addstar.h" + +int dump_rawstars (Stars *stars, int Nstars) { + + int i; + FILE *f; + + f = fopen ("stars.dat", "w"); + + for (i = 0; i < Nstars; i++) { + fprintf (f, "%4d %10.6f %10.6f %8.2f %8.2f %7.2f %7.2f\n", + i, + stars[i].R, stars[i].D, + stars[i].X, stars[i].Y, + stars[i].M, stars[i].dM); + } + + fclose (f); + exit (1); +} + Index: /tags/addstar-1-2/Ohana/src/addstar/src/edge_check.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/edge_check.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/edge_check.c (revision 4773) @@ -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: /tags/addstar-1-2/Ohana/src/addstar/src/find_matches.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/find_matches.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/find_matches.c (revision 4773) @@ -0,0 +1,398 @@ +# include "addstar.h" + +void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *overlap, int Noverlap) { + + 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; + int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch; + Measure *tmpmeasure; + Missing *tmpmissing; + Coords tcoords; + int Nsecfilt, Nsec, MTIME; + short int Mrel, *Mval; + + /* photcode data - must by of type DEP, (PRI, SEC) - probably should restrict to DEP */ + Nsecfilt = GetPhotcodeNsecfilt (); + Nsec = (thiscode[0].type == PHOT_DEP) ? + GetPhotcodeNsec (thiscode[0].equiv) : + GetPhotcodeNsec (thiscode[0].code); + /* this function requires incoming stars to have the same photcode */ + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, float, Nstars); + ALLOCATE (Y1, float, Nstars); + ALLOCATE (N1, int, Nstars); + + /** allocate local arrays (catalog) **/ + Nave = catalog[0].Naverage; + NAVE = Nave + 1000; + ALLOCATE (X2, float, NAVE); + ALLOCATE (Y2, float, NAVE); + ALLOCATE (N2, int, NAVE); + ALLOCATE (catalog[0].found, int, NAVE); + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); + /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */ + + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + + Nmiss = catalog[0].Nmissing; + NMISS = Nmiss + 1000; + ALLOCATE (next_miss, int, NMISS); + REALLOCATE (catalog[0].missing, Missing, NMISS); + + /* project onto rectilinear grid with 1 arcsec pixels */ + /* we keep the original crpix1,2 and crref1,2 */ + /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + tcoords = MOSAIC[0]; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 1; + strcpy (tcoords.ctype, "RA---TAN"); + } else { + tcoords = image[0].coords; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + } + + for (i = 0; i < Nstars; i++) { + fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].D, &tcoords); + N1[i] = i; + } + if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars); + + for (i = 0; i < Nave; i++) { + fRD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + catalog[0].found[N2[i]] = -1; + } + if (Nave > 1) sort_lists (X2, Y2, N2, Nave); + + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + for (i = 0; i < Nmiss - 1; i++) { + next_miss[i] = i+1; + } + next_miss[i] = -1; + last_miss = i; + + /* choose a radius for matches */ + if (DEFAULT_RADIUS == 0) { + RADIUS = NSIGMA * 0.02 * image[0].cerror; /* 0.02 corrects cerror to arcsec from storage units */ + } else { + RADIUS = DEFAULT_RADIUS; /* provided by config */ + } + RADIUS2 = RADIUS*RADIUS; + + /* correct instrumental mags for exposure time */ + if (image[0].exptime > 0) { + MTIME = 2500*log10(image[0].exptime); + } else { + MTIME = 0; + } + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + + dX = X1[i] - X2[j]; + if (dX <= -2*RADIUS) { + i++; + continue; + } + if (dX >= 2*RADIUS) { + j++; + continue; + } + + /* negative dX: j is too large, positive dX, i is too large */ + 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) 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 = 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); + + /*** this should probably use the mag lookup functions - zeropoints are incorrectly applied ***/ + /* 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; + + /** 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. **/ + if (CALIBRATE) { + 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++; + } + + /* add reference for undetected catalog stars */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) RegisterMosaic (MOSAIC); + for (j = 0; (j < Nave) && !SKIP_MISSED; j++) { + n = N2[j]; + 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)) 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; + } + } + /* 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); + } + } + } + + /* incorporate unmatched image stars, if this star is in field of this catalog */ + for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) { + N = N1[i]; + 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 (!FindMosaicForImage (overlap, Noverlap, j)) continue; + 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); + } + } + + free (catalog[0].found); + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + tmpmeasure[N].averef = i; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + /* fix order of Missing (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmissing, Missing, Nmiss); + for (i = 0; i < Nave; i++) { + if (catalog[0].average[i].Nn > 0) { + n = catalog[0].average[i].missing; + catalog[0].average[i].missing = N; + for (k = 0; k < catalog[0].average[i].Nn; k++, N++) { + tmpmissing[N] = catalog[0].missing[n]; + n = next_miss[n]; + } + } + } + free (catalog[0].missing); + catalog[0].missing = tmpmissing; + + /* note stars which have been found in this catalog */ + for (i = 0; i < Nstars; i++) { + if (stars[i].found > -1) { + stars[i].found = -2; + } else { + stars[i].found = -3; + } + } + + /* check if the catalog has changed? if no change, no need to write */ + catalog[0].Naverage = Nave; + catalog[0].Nmeasure = Nmeas; + catalog[0].Nmissing = Nmiss; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch); +} + +/* + 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 + +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/find_matches_refstars.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/find_matches_refstars.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/find_matches_refstars.c (revision 4773) @@ -0,0 +1,297 @@ +# include "addstar.h" + +void find_matches_refstars (GSCRegion *region, Stars **stars, int Nstars, Catalog *catalog) { + + int i, j, k, n, m, N, J; + double RADIUS, RADIUS2; + float *X1, *Y1, *X2, *Y2; + float dX, dY, dR; + int *N1, *N2, *next, *next_miss, last, last_miss; + int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch; + Measure *tmpmeasure; + Missing *tmpmissing; + Coords tcoords; + int Nsecfilt; + + /* photcode data -- should not have to modify secfilt / average */ + Nsecfilt = GetPhotcodeNsecfilt (); + + /** allocate local arrays (stars) **/ + ALLOCATE (X1, float, Nstars); + ALLOCATE (Y1, float, Nstars); + ALLOCATE (N1, int, Nstars); + + /** allocate local arrays (catalog) **/ + Nave = catalog[0].Naverage; + NAVE = Nave + 1000; + ALLOCATE (X2, float, NAVE); + ALLOCATE (Y2, float, NAVE); + ALLOCATE (N2, int, NAVE); + ALLOCATE (catalog[0].found, int, NAVE); + REALLOCATE (catalog[0].average, Average, NAVE); + REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt); + /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */ + + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + + Nmiss = catalog[0].Nmissing; + NMISS = Nmiss + 1000; + ALLOCATE (next_miss, int, NMISS); + REALLOCATE (catalog[0].missing, Missing, NMISS); + + /* project onto rectilinear grid with 1 arcsec pixels, sort by X */ + /* reference for coords is catalog center */ + tcoords.crval1 = 0.5*(region[0].RA[0] + region[0].RA[1]); + tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]); + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + tcoords.Npolyterms = 0; + + for (i = 0; i < Nstars; i++) { + fRD_to_XY (&X1[i], &Y1[i], stars[i][0].R, stars[i][0].D, &tcoords); + N1[i] = i; + } + if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars); + + for (i = 0; i < Nave; i++) { + fRD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords); + N2[i] = i; + catalog[0].found[N2[i]] = -1; + } + if (Nave > 1) sort_lists (X2, Y2, N2, Nave); + + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + for (i = 0; i < Nmiss - 1; i++) { + next_miss[i] = i+1; + } + next_miss[i] = -1; + last_miss = i; + + /* choose a radius for matches */ + if (DEFAULT_RADIUS == 0) { + RADIUS = 2.0; /* hardwired default for refstars */ + } else { + RADIUS = DEFAULT_RADIUS; /* provided by config */ + } + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nstars) && (j < Nave); ) { + + dX = X1[i] - X2[j]; + if (dX <= -2*RADIUS) { + i++; + continue; + } + if (dX >= 2*RADIUS) { + j++; + continue; + } + + /* negative dX: j is too large, positive dX, i is too large */ + 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) 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][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].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; + + 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; + + /** 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) { + 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++; + } + + /* 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); + } + } + + free (catalog[0].found); + REALLOCATE (catalog[0].average, Average, Nave); + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + tmpmeasure[N].averef = i; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + /* fix order of Missing (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmissing, Missing, Nmiss); + for (i = 0; i < Nave; i++) { + if (catalog[0].average[i].Nn > 0) { + n = catalog[0].average[i].missing; + catalog[0].average[i].missing = N; + for (k = 0; k < catalog[0].average[i].Nn; k++, N++) { + tmpmissing[N] = catalog[0].missing[n]; + n = next_miss[n]; + } + } + } + free (catalog[0].missing); + catalog[0].missing = tmpmissing; + + /* note stars which have been found in this catalog */ + for (i = 0; i < Nstars; i++) { + if (stars[i][0].found > -1) { + stars[i][0].found = -2; + } else { + stars[i][0].found = -3; + } + } + + catalog[0].Naverage = Nave; + catalog[0].Nmeasure = Nmeas; + catalog[0].Nmissing = Nmiss; + if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/find_proper.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/find_proper.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/find_proper.c (revision 4773) @@ -0,0 +1,158 @@ +# include "addusno.h" + +find_proper (catstats, catalog, usnostats, usno, Nusno) +CatStats catstats[]; +Catalog catalog[]; +USNOstats usnostats[]; +USNOdata usno[]; +int Nusno; +{ + + int i, j, k, n, m, N, first_j; + double RADIUS2, PROPER2; + float *X1, *Y1, *X2, *Y2, *usnodist; + float dX, dY, dR, dR2; + int *N1, *N2, *next, last; + int Nave, Nmeas, NMEAS, Nmatch; + unsigned int flags; + Measure *tmpmeasure; + Coords *tcoords; + int already_matched, far_enough; + + X1 = catstats[0].X; + Y1 = catstats[0].Y; + N1 = catstats[0].N; + Nave = catalog[0].Naverage; + + /* no need to do this twice!! */ + ALLOCATE (usnodist, float, Nave); + bzero (usnodist, Nave*sizeof(float)); + + X2 = usnostats[0].X; + Y2 = usnostats[0].Y; + N2 = usnostats[0].N; + + /* set up link listed pointers for new measurements */ + Nmatch = 0; + Nmeas = catalog[0].Nmeasure; + NMEAS = Nmeas + 1000; + ALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + /* set up pointers for linked list of measurements */ + for (i = 0; i < Nmeas - 1; i++) { + next[i] = i+1; + } + next[i] = -1; + last = i; + + /* choose a radius for matches */ + PROPER2 = PROPER*PROPER; + RADIUS2 = RADIUS*RADIUS; + + /** find matched stars **/ + for (i = j = 0; (i < Nave) && (j < Nusno); ) { + if (catalog[0].average[N1[i]].code & ID_MOVING) { + /* this is not a star, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].Nm < 3) { + /* may just be a noise spike, skip */ + i++; + continue; + } + if (catalog[0].average[N1[i]].code & ID_USNO) { + /* already matched with USNO, skip this one */ + i++; + continue; + } + dX = X1[i] - X2[j]; + if (dX <= -2*PROPER) { + i++; + continue; + } + if (dX >= 2*PROPER) { + j++; + continue; + } + /* negative dX: j is too large, positive dX, i is too large */ + first_j = j; + for (; (dX > -2*PROPER) && (j < Nusno); j++) { + dX = X1[i] - X2[j]; + dY = Y1[i] - Y2[j]; + dR = dX*dX + dY*dY; + if (dR < PROPER2) { + n = N1[i]; /* N1 refers to the average[] list */ + N = N2[j]; /* N2 refers to the usno[] list */ + if (usnostats[0].match[N] > -1) + continue; + if ((catalog[0].average[n].code & ID_USNO) && (usnodist[i] < dR)) { + /* existing USNO match is closer than this new one, skip this one */ + continue; + } + usnodist[i] = dR; + m = catalog[0].average[n].offset; /* first measurement of this star */ + for (k = 0; k < catalog[0].average[n].Nm - 1; k++) + m = next[m]; + next[Nmeas+1] = next[m]; /* insert 2 measurements in linked list */ + next[Nmeas] = Nmeas + 1; + next[m] = Nmeas; + if (next[Nmeas+1] == -1) { /* last just was moved */ + last = Nmeas+1; + } + Nmatch ++; + + /** add measurements for this star **/ + catalog[0].measure[Nmeas].dR = 360000.0*(catalog[0].average[n].R - usno[N].R); + catalog[0].measure[Nmeas].dD = 360000.0*(catalog[0].average[n].D - usno[N].D); + catalog[0].measure[Nmeas].M = 1000.0*fabs(usno[N].r); + catalog[0].measure[Nmeas].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas].averef = n; + catalog[0].measure[Nmeas].source = USNO_RED; + catalog[0].measure[Nmeas+1].dR = catalog[0].measure[Nmeas].dR; + catalog[0].measure[Nmeas+1].dD = catalog[0].measure[Nmeas].dD; + catalog[0].measure[Nmeas+1].M = 1000.0*fabs(usno[N].b); + catalog[0].measure[Nmeas+1].Mcal= 0; /* above measurement is exact */ + catalog[0].measure[Nmeas+1].dM = 100; /* error in input files stored in thousandths of mag */ + catalog[0].measure[Nmeas+1].t = 0; /* a flag: if 0, image is not in database */ + catalog[0].measure[Nmeas+1].averef = n; + catalog[0].measure[Nmeas+1].source = USNO_BLUE; + /* add flag in average to mark as matched with the USNO catalog */ + catalog[0].average[n].code |= (ID_PROPER | ID_USNO); + + /* we add two measurement for each star: 1 in r, 1 in b */ + catalog[0].average[n].Nm +=2; + Nmeas +=2; + if (Nmeas == NMEAS - 2) { /* just to be safe... */ + NMEAS = Nmeas + 1000; + REALLOCATE (next, int, NMEAS); + REALLOCATE (catalog[0].measure, Measure, NMEAS); + } + } + } + j = first_j; + i++; + } + + REALLOCATE (catalog[0].measure, Measure, Nmeas); + + /* fix order of Measure (memory intensive, but fast) */ + N = 0; + ALLOCATE (tmpmeasure, Measure, Nmeas); + for (i = 0; i < Nave; i++) { + n = catalog[0].average[i].offset; + catalog[0].average[i].offset = N; + for (k = 0; k < catalog[0].average[i].Nm; k++, N++) { + tmpmeasure[N] = catalog[0].measure[n]; + n = next[n]; + } + } + free (catalog[0].measure); + catalog[0].measure = tmpmeasure; + + catalog[0].Nmeasure = Nmeas; + if (VERBOSE) fprintf (stderr, "Nusno, Nave, Nmeas: %d %d %d, (%d matches)\n", Nusno, Nave, Nmeas, Nmatch); +} + Index: /tags/addstar-1-2/Ohana/src/addstar/src/find_subset.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/find_subset.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/find_subset.c (revision 4773) @@ -0,0 +1,36 @@ +# include "addstar.h" + +/* find stars within this region */ +Stars **find_subset (GSCRegion *region, Stars *stars, int Nstars, int *Nsubset) { + + int i, N, NSUBSET; + Stars **subset; + double RA0, RA1, DEC0, DEC1; + + NSUBSET = 1000; + ALLOCATE (subset, Stars *, NSUBSET); + + RA0 = region[0].RA[0]; + RA1 = region[0].RA[1]; + DEC0 = region[0].DEC[0]; + DEC1 = region[0].DEC[1]; + + if (VERBOSE) fprintf (stderr, "subset for %f - %f, %f - %f\n", RA0, RA1, DEC0, DEC1); + + /* find stars within ra,dec region */ + for (i = N = 0; i < Nstars; i++) { + if (stars[i].R < RA0) continue; + if (stars[i].R >= RA1) continue; + if (stars[i].D < DEC0) continue; + if (stars[i].D >= DEC1) continue; + + subset[N] = &stars[i]; + N++; + if (N == NSUBSET - 1) { + NSUBSET += 1000; + REALLOCATE (subset, Stars *, NSUBSET); + } + } + *Nsubset = N; + return (subset); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/gcatalog.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gcatalog.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gcatalog.c (revision 4773) @@ -0,0 +1,47 @@ +# include "addstar.h" + +int gcatalog (Catalog *catalog) { + + int Nsecfilt; + char mode; + + /* read catalog header */ + mode = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + if (!load_catalog (catalog, mode, VERBOSE)) { + fprintf (stderr, "ERROR: failure loading catalog\n"); + exit (1); + } + + /* check Nsecfile value, update if needed */ + Nsecfilt = GetPhotcodeNsecfilt (); + if (catalog[0].Nsecfilt < Nsecfilt) { + + int i, j, Nextra, in, out; + SecFilt *insec, *outsec; + + Nextra = Nsecfilt - catalog[0].Nsecfilt; + insec = catalog[0].secfilt; + ALLOCATE (outsec, SecFilt, catalog[0].Naverage * Nsecfilt); + 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].Xm = insec[in].Xm; + outsec[out].dM = insec[in].dM; + } + for (j = 0; j < Nextra; j++, out++) { + outsec[out].M = NO_MAG; + outsec[out].Xm = NO_MAG; + outsec[out].dM = NO_MAG; + } + } + free (catalog[0].secfilt); + catalog[0].secfilt = outsec; + catalog[0].Nsecfilt = Nsecfilt; + } + + if (catalog[0].Nsecfilt > Nsecfilt) { + fprintf (stderr, "ERROR: can't reduce number of secondary filters\n"); + exit (1); + } + return (TRUE); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/get2mass.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/get2mass.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/get2mass.c (revision 4773) @@ -0,0 +1,95 @@ +# include "addstar.h" +# include "2mass.h" + +Stars *get2mass (GSCRegion *patch, int *NSTARS, int mode) { + + char *path; + int i, Nstars, Nrefcat, Nregions; + Stars *stars; + Stars *refcat; + GSCRegion *regions; + + NAMED_PHOTCODE (TM_J, "2MASS_J"); + NAMED_PHOTCODE (TM_H, "2MASS_H"); + NAMED_PHOTCODE (TM_K, "2MASS_K"); + if (thiscode == NULL) Shutdown ("photcode not specified"); + if (thiscode[0].code == TM_J) goto good_code; + if (thiscode[0].code == TM_H) goto good_code; + if (thiscode[0].code == TM_K) goto good_code; + Shutdown ("2MASS photcode not specified"); + +good_code: + path = TWO_MASS_DIR_AS; + if (mode == 1) path = TWO_MASS_DIR_DR2; + + regions = get2mass_acc (patch, &Nregions, path, "accel.dat"); + + Nstars = 0; + ALLOCATE (stars, Stars, 1); + + for (i = 0; i < Nregions; i++) { + switch (mode) { + case 0: + refcat = get2mass_AS_data (®ions[i], patch, &Nrefcat); + if (VERBOSE) fprintf (stderr, "loaded %d stars from 2MASS (allsky)\n", Nrefcat); + break; + case 1: + refcat = get2mass_2DR_data (®ions[i], patch, &Nrefcat); + if (VERBOSE) fprintf (stderr, "loaded %d stars from 2MASS (dr2)\n", Nrefcat); + break; + } + + REALLOCATE (stars, Stars, MAX (1, Nstars + Nrefcat)); + memcpy (&stars[Nstars], refcat, Nrefcat*sizeof(Stars)); + Nstars += Nrefcat; + + free (refcat); + } + + if (VERBOSE) fprintf (stderr, "loaded total %d stars from 2MASS\n", Nstars); + *NSTARS = Nstars; + return (stars); +} + +/* watch for patch which cross 0,360 boundary */ +GSCRegion *get2mass_acc (GSCRegion *patch, int *nregions, char *path, char *accel) { + + int Nregions, NREGIONS; + GSCRegion *regions; + char accelfile[1024], line[256], filename[128]; + FILE *f; + double Rs, Re, Ds, De; + + sprintf (accelfile, "%s/%s", path, accel); + f = fopen (accelfile, "r"); + if (f == NULL) Shutdown ("can't read data from accelerator %s", accelfile); + + Nregions = 0; + NREGIONS = 200; + ALLOCATE (regions, GSCRegion, NREGIONS); + + /* read in stars line-by-line */ + while (scan_line (f, line) != EOF) { + stripwhite (line); + if (line[0] == 0) continue; + if (line[0] == '#') continue; + sscanf (line, "%s %lf %lf %lf %lf", filename, &Rs, &Re, &Ds, &De); + Rs *= 15.0; + Re *= 15.0; + if (Rs > patch[0].RA[1]) continue; + if (Re < patch[0].RA[0]) continue; + if (Ds > patch[0].DEC[1]) continue; + if (De < patch[0].DEC[0]) continue; + sprintf (regions[Nregions].filename, "%s/%s", path, filename); + regions[Nregions].RA[0] = Rs; + regions[Nregions].RA[1] = Re; + regions[Nregions].DEC[0] = Ds; + regions[Nregions].DEC[1] = De; + Nregions ++; + CHECK_REALLOCATE (regions, GSCRegion, NREGIONS, Nregions, 20); + } + fclose (f); + + *nregions = Nregions; + return (regions); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_as.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_as.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_as.c (revision 4773) @@ -0,0 +1,205 @@ +# include "addstar.h" +# include "2mass.h" + +/* unlike the DR2 data, the AS data is NOT fixed bytes/row + * we need to handle fractional lines at the end of each read block + */ + +/* read in chunks of ~16MB */ +# define NBYTE 0x1000000 + +char *skipNbounds (char *line, char bound, int Nbound, int Nbyte); +void *memrchr(const void *s, int c, size_t n); +e_time get2mass_time (char *ptr, char *buffer, int Nbound, int Nbyte); + +Stars *get2mass_AS_data (GSCRegion *region, GSCRegion *patch, int *nstars) { + + int FilterSkip, TimeSkip; + int Nstars, NSTARS, Nbyte, Nextra, code, Qentry; + Stars *stars; + gzFile gf; + char qc; + char *buffer, *filename; + char *p, *q, *ptr, *tmp; + double RA, DEC, M, dM; + double RA0, RA1, DEC0, DEC1; + e_time time; + + ALLOCATE (buffer, char, NBYTE); + + RA0 = MAX (patch[0].RA[0], UserPatch.RA[0]); + RA1 = MIN (patch[0].RA[1], UserPatch.RA[1]); + DEC0 = MAX (patch[0].DEC[0], UserPatch.DEC[0]); + DEC1 = MIN (patch[0].DEC[1], UserPatch.DEC[1]); + + code = thiscode[0].code; + FilterSkip = TimeSkip = 0; + if (thiscode[0].code == TM_J) { + FilterSkip = 6; + TimeSkip = 23; + Qentry = 0; + } + if (thiscode[0].code == TM_H) { + FilterSkip = 10; + TimeSkip = 19; + Qentry = 1; + } + if (thiscode[0].code == TM_K) { + FilterSkip = 14; + TimeSkip = 15; + Qentry = 2; + } + if (!FilterSkip) Shutdown ("invalid photcode %d", thiscode[0].code); + + filename = region[0].filename; + gf = gzopen (filename, "rb"); + if (gf == NULL) Shutdown ("can't read 2mass data file: %s", filename); + + Nstars = 0; + NSTARS = 10000; + ALLOCATE (stars, Stars, NSTARS); + + Nextra = 0; + while ((Nbyte = gzread (gf, &buffer[Nextra], NBYTE-Nextra)) != 0) { + if (Nbyte == -1) Shutdown ("error reading from gzipped file %s", filename); + Nbyte += Nextra; + + if (VERBOSE) fprintf (stderr, "."); + + /* find bounds on first complete line */ + p = buffer; + q = memchr (p, '\n', Nbyte); + if (q == NULL) Shutdown ("incomplete line at end of file\n"); + + while (1) { + + RA = strtod (p, NULL); + ptr = skipNbounds (p, '|', 1, Nbyte - (p - buffer)); + DEC = strtod (ptr, NULL); + if (DEC > 90) Shutdown ("weird DEC value: something is wrong"); + + /* skip stars which are outside desired region */ + if (DEC > DEC1) goto skip_star; + if (DEC < DEC0) goto skip_star; + if (RA < RA0) goto skip_star; + if (RA > RA1) goto skip_star; + + ptr = skipNbounds (p, '|', FilterSkip, Nbyte - (p - buffer)); + if (ptr == NULL) Shutdown ("format error in 2mass"); + M = strtod (ptr, NULL); + ptr = skipNbounds (ptr, '|', 1, Nbyte - (ptr - buffer)); + dM = strtod (ptr, NULL); + time = get2mass_time (ptr, buffer, TimeSkip, Nbyte); + + /* filter on the ph_qual flag for this filter (field 19) */ + if (SELECT_2MASS_QUALITY != NULL) { + ptr = skipNbounds (p, '|', 18, Nbyte - (p - buffer)); + qc = ptr[Qentry]; + if (strchr (SELECT_2MASS_QUALITY, qc) == NULL) goto skip_star; + } + + stars[Nstars].R = RA; + stars[Nstars].D = DEC; + stars[Nstars].M = M; + stars[Nstars].dM = dM; + stars[Nstars].code = code; + stars[Nstars].t = time; + stars[Nstars].found = -1; + + Nstars ++; + CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 5000); + + skip_star: + /* start of the next line */ + tmp = p; + p = q + 1; + if (p - buffer == Nbyte) { + Nextra = 0; + break; + } + /* end of the next line */ + q = memchr (p, '\n', Nbyte - (p - buffer)); + if (q == NULL) { + Nextra = Nbyte - (p - buffer); + memmove (buffer, p, Nextra); + break; + } + } + } + if (VERBOSE) fprintf (stderr, "\n"); + + gzclose (gf); + free (buffer); + *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. +*/ +/* don't bother to seek ahead : position is not sufficiently predictable + and gzseek is as expensive as gzread */ +/* + Noffset = region[0].Nrec * (patch[0].DEC[0] + 90) / 180.0; + gzseek (gf, Noffset * NBYTE, SEEK_SET); + Nbyte = gzread (gf, buffer, NLINE*NBYTE); +*/ + +/* return a pointer to the first char after Nbound of value bound */ +char *skipNbounds (char *line, char bound, int Nbound, int Nbyte) { + + int i; + char *p, *q; + + p = line; + for (i = 0; i < Nbound; i++) { + q = memchr (p, bound, Nbyte - (p - line)); + if (q == NULL) return (NULL); + p = q + 1; + if (p - line == Nbyte) return (NULL); + } + return (p); +} + +e_time get2mass_time (char *ptr, char *buffer, int Nbound, int Nbyte) { + + e_time time; + char *p, *end; + + p = skipNbounds (ptr, '|', Nbound, Nbyte - (ptr - buffer)); + if (p == NULL) Shutdown ("format error in 2mass"); + end = memchr (p, '|', Nbyte - (p - buffer)); + if (end == NULL) Shutdown ("format error in 2mass"); + *end = 0; + time = date_to_sec (ptr); + *end = '|'; + + return (time); +} + +# if (0) + /** need to re-think this test **/ + if (0) { + /* search for end of last complete line */ + p = memrchr (buffer, '\n', Nbyte); + if (p == NULL) Shutdown ("incomplete line in at end of file\n"); + + /* search for start of last complete line */ + /* last block may be only one line */ + q = memrchr (buffer, '\n', (p - buffer - 1)); + if (q == NULL) { + q = buffer; + } else { + q ++; + } + + /* skip past block not yet in range */ + RA = strtod (q, NULL); + ptr = skipNbounds (q, '|', 1, Nbyte - (q - buffer)); + DEC = strtod (ptr, NULL); + if (DEC < DEC0) continue; + } + +# endif Index: /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_dr2.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_dr2.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/get2mass_dr2.c (revision 4773) @@ -0,0 +1,111 @@ +# include "addstar.h" +# include "2mass.h" + +# define NBYTE 302 +# define NLINE 30000 + +Stars *get2mass_2DR_data (GSCRegion *region, GSCRegion *patch, int *nstars) { + + int i, Nstars, NSTARS, Nbyte, Nline; + Stars *stars; + gzFile gf; + char *buffer, *filename; + char line[303]; + double RA, DEC, J, H, K, dJ, dH, dK; + double RA0, RA1, DEC0, DEC1; + + ALLOCATE (buffer, char, (NBYTE*NLINE)); + + if (region == NULL) exit (2); + if (patch == NULL) exit (3); + + RA0 = MAX (patch[0].RA[0], UserPatch.RA[0]); + RA1 = MIN (patch[0].RA[1], UserPatch.RA[1]); + DEC0 = MAX (patch[0].DEC[0], UserPatch.DEC[0]); + DEC1 = MIN (patch[0].DEC[1], UserPatch.DEC[1]); + + fprintf (stderr, "overlap: %f - %f, %f - %f\n", RA0, RA1, DEC0, DEC1); + + filename = region[0].filename; + gf = gzopen (filename, "rb"); + if (gf == NULL) Shutdown ("can't read 2mass data file: %s", filename); + + Nstars = 0; + NSTARS = 10000; + ALLOCATE (stars, Stars, NSTARS); + + while ((Nbyte = gzread (gf, buffer, NLINE*NBYTE)) != 0) { + if (Nbyte == 0) Shutdown ("error reading from gzipped file %s", filename); + if (Nbyte == -1) Shutdown ("error reading from gzipped file %s", filename); + if (Nbyte % NBYTE) Shutdown ("error reading complete line from gzipped file %s", filename); + Nline = Nbyte / NBYTE; + + /* skip past block not yet in range */ + sscanf (&buffer[NBYTE*(Nline - 1)], "%lf %lf", &RA, &DEC); + if (DEC < DEC0) continue; + + memcpy (line, &buffer[NBYTE*(Nline-1)], NBYTE); + line[302] = 0; + + for (i = 0; i < Nline; i++) { + + dparse (&RA, 1, &buffer[NBYTE*i+ 0]); + dparse (&DEC, 2, &buffer[NBYTE*i+ 0]); + + /* dr2 is nicely sorted in dec order */ + if (DEC > DEC1) goto finished; + if (DEC < DEC0) continue; + if (RA < RA0) continue; + if (RA > RA1) continue; + + stars[Nstars].R = RA; + stars[Nstars].D = DEC; + stars[Nstars].t = short_date_to_sec (&buffer[NBYTE*i + 164]); + stars[Nstars].found = -1; + + if (thiscode[0].code == TM_J) { + dparse (&J, 1, &buffer[NBYTE*i + 53]); + dparse (&dJ, 2, &buffer[NBYTE*i + 53]); + stars[Nstars].M = J; + stars[Nstars].dM = dJ; + stars[Nstars].code = TM_J; + } + if (thiscode[0].code == TM_H) { + dparse (&H, 1, &buffer[NBYTE*i + 72]); + dparse (&dH, 2, &buffer[NBYTE*i + 72]); + stars[Nstars].M = H; + stars[Nstars].dM = dH; + stars[Nstars].code = TM_H; + } + if (thiscode[0].code == TM_K) { + dparse (&K, 1, &buffer[NBYTE*i + 91]); + dparse (&dK, 2, &buffer[NBYTE*i + 91]); + stars[Nstars].M = K; + stars[Nstars].dM = dK; + stars[Nstars].code = TM_K; + } + Nstars ++; + CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 5000); + } + } +finished: + gzclose (gf); + free (buffer); + + *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. +*/ +/* don't bother to seek ahead : position is not sufficiently predictable + and gzseek is as expensive as gzread */ +/* + Noffset = region[0].Nrec * (patch[0].DEC[0] + 90) / 180.0; + gzseek (gf, Noffset * NBYTE, SEEK_SET); + Nbyte = gzread (gf, buffer, NLINE*NBYTE); +*/ + Index: /tags/addstar-1-2/Ohana/src/addstar/src/getgsc.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/getgsc.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/getgsc.c (revision 4773) @@ -0,0 +1,87 @@ +# include "addstar.h" +# define BYTES_STAR 23 +# define BLOCK 1000 +# define DNSTARS 1000 + +static short GSC_M; + +Stars *getgsc (GSCRegion *patch, int *NSTARS) { + + int i, Ngsc, Nregions, Nstars; + Stars *gsc; + Stars *stars; + GSCRegion *regions; + char *tmp; + + NAMED_PHOTCODE (GSC_M, "GSC_M"); + + /* sort of sleaze way to get the right GSC files */ + tmp = strcreate (CATDIR); + strcpy (CATDIR, GSCDIR); + regions = gregion_patch (patch, &Nregions); + strcpy (CATDIR, tmp); + free (tmp); + + Nstars = 0; + ALLOCATE (stars, Stars, 1); + + for (i = 0; i < Nregions; i++) { + gsc = rd_gsc (regions[i].filename, &Ngsc); + + REALLOCATE (stars, Stars, MAX (1, Nstars + Ngsc)); + memcpy (&stars[Nstars], gsc, Ngsc*sizeof(Stars)); + Nstars += Ngsc; + + free (gsc); + } + + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Nstars); + *NSTARS = Nstars; + return (stars); +} + +Stars *rd_gsc (char *filename, int *nstars) { + + int i, NSTARS, Nstars, Nbytes, nbytes, Nline, Nbyte; + 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]); + if (stars[Nstars].R < UserPatch.RA[0]) continue; + if (stars[Nstars].R > UserPatch.RA[1]) continue; + if (stars[Nstars].D < UserPatch.DEC[0]) continue; + if (stars[Nstars].D > UserPatch.DEC[1]) continue; + + dparse (&stars[Nstars].M, 3, &buffer[Nbyte]); + stars[Nstars].dM = NO_ERR; + stars[Nstars].t = 0; + stars[Nstars].code = GSC_M; + stars[Nstars].found = -1; + Nstars ++; + CHECK_REALLOCATE (stars, Stars, NSTARS, Nstars, 1000); + } + } + free (buffer); + REALLOCATE (stars, Stars, Nstars); + *nstars = Nstars; + return (stars); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/getusno.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/getusno.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/getusno.c (revision 4773) @@ -0,0 +1,153 @@ +# include "addstar.h" +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +Stars *getusno (GSCRegion *catstats, int *Nstars) { + + long int offset; + int i, bin, first, last, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128], c; + FILE *f; + int iRA0, iRA1, iDEC0, iDEC1; + double dec; + int spd, spd_start, spd_end, disk; + int NUSNO, Nusno; + Stars *stars; + short int USNO_RED, USNO_BLUE; + + /* require photcode */ + NAMED_PHOTCODE (USNO_RED, "USNO_RED"); + NAMED_PHOTCODE (USNO_BLUE, "USNO_BLUE"); + if (thiscode == NULL) Shutdown ("photcode not specified"); + if (thiscode[0].code == USNO_RED) goto good_code; + if (thiscode[0].code == USNO_BLUE) goto good_code; + Shutdown ("USNO photcode not specified"); +good_code: + + /* identify ra & dec range of interest */ + /* note: the use of UserPatch to restrict here limits general utility of function */ + iRA0 = MAX (catstats[0].RA[0], UserPatch.RA[0]) * 360000.0; + iRA1 = MIN (catstats[0].RA[1], UserPatch.RA[1]) * 360000.0; + iDEC0 = (MAX (catstats[0].DEC[0], UserPatch.DEC[0]) + 90.0) * 360000.0; + iDEC1 = (MIN (catstats[0].DEC[1], UserPatch.DEC[1]) + 90.0) * 360000.0; + + /* data is organized in south-pole distance zones */ + 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 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } + + Nusno = 0; + NUSNO = 5000; + ALLOCATE (stars, Stars, NUSNO); + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (1); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = 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 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s\n", filename); + exit (1); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failure reading data from file %s\n", filename); + exit (1); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if (buf[0] < iRA0) 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].dM = NO_ERR; + stars[Nusno].t = 0; + stars[Nusno].found = -1; + + /* one pass of addstar does either r or b */ + if (thiscode[0].code == USNO_RED) { + stars[Nusno].code = USNO_RED; + stars[Nusno].M = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); + } + if (thiscode[0].code == USNO_BLUE) { + stars[Nusno].code = USNO_BLUE; + stars[Nusno].M = fabs (0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000)); + } + Nusno ++; + CHECK_REALLOCATE (stars, Stars, NUSNO, Nusno, 5000); + } + free (buffer); + } + fclose (f); + } + + *Nstars = Nusno; + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno); + return (stars); +} + + Index: /tags/addstar-1-2/Ohana/src/addstar/src/gimages.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gimages.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gimages.c (revision 4773) @@ -0,0 +1,160 @@ +# include "addstar.h" + +/* given image, find catalog images which overlap it */ +Image *gimages (Image *image, int *Npimage) { + + int i, j, k, Nimage, addtolist, size, dbstate; + int NTIMAGE, Ntimage, ntimage; + int NPIMAGE, npimage; + Image *timage, *pimage; + Coords tcoords; + Header header; + double r, d; + double Xi[5], Yi[5], Xo[5], Yo[5]; /* image and original corners */ + double zeropt; + struct stat filestatus; + FILE *f; + + /* if image.db is opened successfully, f will not be NULL. + if image.db is empty, dbstate will be LCK_EMPTY + if image.db is not locked & opened successfully, we will not reach here! + */ + f = GetDB (&dbstate); + if (dbstate == LCK_EMPTY) { + *Npimage = 0; + return (NULL); + } + fseek (f, 0, SEEK_SET); + + /* read header */ + if (!fits_load_header (f, &header)) { + Shutdown ("can't read image catalog %s", ImageCat); + } + + fits_scan (&header, "ZERO_PT", "%lf", 1, &zeropt); + if (fabs (ZeroPt - zeropt) > 1e-4) { + Shutdown ("zero point in image (%f:%s) inconsistent with zero point in catalog (%f)", + ZeroPt, image[0].name, zeropt); + } + + /* project onto rectilinear grid with 1 arcsec pixels */ + /* we keep the original crpix1,2 and crref1,2 */ + /* for mosaic astrometry, the grid should be w.r.t. the tangent-plane, not chip coords */ + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + tcoords = MOSAIC[0]; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + tcoords.Npolyterms = 1; + strcpy (tcoords.ctype, "RA---TAN"); + RegisterMosaic (MOSAIC); + /* register so image->sky conversions below have correct mosaic */ + } else { + tcoords = image[0].coords; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + } + + /* define original corners */ + Xo[0] = 0; Yo[0] = 0; + Xo[1] = image[0].NX; Yo[1] = 0; + Xo[2] = image[0].NX; Yo[2] = image[0].NY; + Xo[3] = 0; Yo[3] = image[0].NY; + Xo[4] = 0; Yo[4] = 0; /* so we can make a loop easily */ + for (j = 0; j < 5; j++) { + /* XY-to_RD is two-level if ctype == WRP */ + XY_to_RD (&r, &d, Xo[j], Yo[j], &image[0].coords); + RD_to_XY (&Xo[j], &Yo[j], r, d, &tcoords); + } + + /* advance to data segment */ + fseek (f, header.size, SEEK_SET); + + Nimage = 0; + fits_scan (&header, "NIMAGES", "%d", 1, &Nimage); + /* check that file size makes sense */ + if (stat (ImageCat, &filestatus) == -1) { + Shutdown ("failed to get status of image catalog"); + } + size = Nimage*sizeof(Image) + header.size; + if (size != filestatus.st_size) { + int Ndata; + + Ndata = (filestatus.st_size - header.size) / sizeof (Image); + fprintf (stderr, "error: image catalog has inconsistent size\n"); + fprintf (stderr, "header: %d images, data: %d images\n", Nimage, Ndata); + if (!FORCE_READ) Shutdown ("refuse to read image table"); + + Nimage = MIN (Nimage, Ndata); + fprintf (stderr, "using %d images\n", Nimage); + } + + /* kind of dangerous, see note below */ + NTIMAGE = MAX (MIN (Nimage, 1000), 1); + ALLOCATE (timage, Image, NTIMAGE); + + npimage = 0; + NPIMAGE = 20; + ALLOCATE (pimage, Image, NPIMAGE); + + /* run through image table, reading images in blocks */ + for (Ntimage = 0; Ntimage < Nimage; Ntimage += ntimage) { + ntimage = Fread (timage, sizeof(Image), NTIMAGE, f, "image"); + if (ntimage == 0) Shutdown ("expected %d images, only found %d", Nimage, Ntimage); + BuildChipMatch (timage, ntimage); + + for (i = 0; i < ntimage; i++) { + + /* if any of these images are WRP images, need to find matching DIS */ + if (!FindMosaicForImage (timage, ntimage, i)) continue; + + addtolist = FALSE; + /* define image corners */ + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; /* so we can make a loop easily */ + /* transform to tcoords */ + for (j = 0; j < 5; j++) { + XY_to_RD (&r, &d, Xi[j], Yi[j], &timage[i].coords); + RD_to_XY (&Xi[j], &Yi[j], r, d, &tcoords); + } + /* check if any edges cross */ + /* not every robust - images must be close in size */ + for (j = 0; (j < 4) && !addtolist; j++) { + for (k = 0; (k < 4) && !addtolist; k++) { + addtolist |= edge_check (&Xi[j], &Yi[j], &Xo[k], &Yo[k]); + } + } + if (addtolist) { + pimage[npimage] = timage[i]; + npimage ++; + if (npimage == NPIMAGE) { + NPIMAGE += 20; + REALLOCATE (pimage, Image, NPIMAGE); + } + } + } + } + + BuildChipMatch (pimage, npimage); + + if (VERBOSE) fprintf (stderr, "found %d overlapping images\n", npimage); + + *Npimage = npimage; + return (pimage); + +} + +/* this function tries to find overlapping images with the current image. + it loads the image table in blocks (up to 1000 at a time) and looks + for overlaps within the set. the search is not very robust. in + addition, it needs to look for mosaic DIS images matching any WRP images. + if the WRP image is loaded near the beginning of a block, then the corresponding DIS + may be in the previous block, and thus missed + + we are not adding the matching DIS image to the table. +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/greference.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/greference.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/greference.c (revision 4773) @@ -0,0 +1,44 @@ +# include "addstar.h" + +Stars *grefcat (char *Refcat, GSCRegion *catstats, int *nstars) { + + int Nstars; + Stars *stars; + + 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; + stars = NULL; + + /* get stars from USNO for the given region */ + if (!strcasecmp (Refcat, "USNO")) { + stars = getusno (catstats, &Nstars); + } + + /* get stars from HST GSC for the given region */ + if (!strcasecmp (Refcat, "GSC")) { + stars = getgsc (catstats, &Nstars); + } + + /* get stars from 2MASS for the given region */ + if (!strcasecmp (Refcat, "2MASS")) { + stars = get2mass (catstats, &Nstars, 0); + } + + /* get stars from 2MASS for the given region */ + if (!strcasecmp (Refcat, "2MASS-ALLSKY")) { + stars = get2mass (catstats, &Nstars, 0); + } + + /* get stars from 2MASS for the given region */ + if (!strcasecmp (Refcat, "2MASS-DR2")) { + stars = get2mass (catstats, &Nstars, 1); + } + + if (VERBOSE && (Nstars == 0)) fprintf (stderr, "no ref objects in region %s\n", catstats[0].filename); + + *nstars = Nstars; + return (stars); + +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/grefstars.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/grefstars.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/grefstars.c (revision 4773) @@ -0,0 +1,41 @@ +# include "addstar.h" + +/* read ASCII file with ref star data */ +Stars *grefstars (char *file, int *Nstars) { + + FILE *f; + int N, NSTARS; + Stars *stars; + char line[256]; + + /* require photcode */ + if (thiscode == NULL) Shutdown ("photcode not specified"); + + /* open file */ + f = fopen (file, "r"); + if (f == NULL) Shutdown ("can't read data from %s", file); + + NSTARS = 100; + ALLOCATE (stars, Stars, NSTARS); + + /* read in stars line-by-line */ + for (N = 0; scan_line (f, line) != EOF; N++) { + stripwhite (line); + if (line[0] == 0) continue; + if (line[0] == '#') continue; + dparse (&stars[N].R, 1, line); + dparse (&stars[N].D, 2, line); + dparse (&stars[N].M, 3, line); + dparse (&stars[N].dM, 4, line); + while (stars[N].R < 0.0) stars[N].R += 360.0; + while (stars[N].R >= 360.0) stars[N].R -= 360.0; + stars[N].t = 0; + stars[N].code = thiscode[0].code; + stars[N].found = FALSE; + CHECK_REALLOCATE (stars, Stars, NSTARS, N+1, 100); + } + *Nstars = N; + return (stars); +} + +/* add the photcode here */ Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregion_image.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregion_image.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregion_image.c (revision 4773) @@ -0,0 +1,33 @@ +# include "addstar.h" + +/* given image with coords, find regions which overlap image */ +GSCRegion *gregion_image (Image *image, int *Nregions) { + + GSCRegion *regions, impatch; + double x, y, X, Y, ra, dec; + + /* this has trouble for images which overlap 0,360 */ + impatch.RA[0] = 360; + impatch.RA[1] = 0; + impatch.DEC[0] = +90; + impatch.DEC[1] = -90; + + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + if (MOSAIC == NULL) Shutdown ("no mosaic for WRP image (use -mosaic)"); + RegisterMosaic (MOSAIC); + } + + for (x = 0; x <= 1; x++) { + for (y = 0; y <= 1; y++) { + X = image[0].NX*(1.1*x - 0.05); + Y = image[0].NY*(1.1*y - 0.05); + XY_to_RD (&ra, &dec, X, Y, &image[0].coords); + impatch.RA[0] = MIN (ra, impatch.RA[0]); + impatch.RA[1] = MAX (ra, impatch.RA[1]); + impatch.DEC[0] = MIN (dec, impatch.DEC[0]); + impatch.DEC[1] = MAX (dec, impatch.DEC[1]); + } + } + regions = gregion_patch (&impatch, Nregions); + return (regions); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match.c (revision 4773) @@ -0,0 +1,29 @@ +# include "addstar.h" + +GSCRegion *gregion_match (GSCRegion *regions, int *nregions) { + + int i, status, Nregions, Nsubset, NSUBSET; + GSCRegion *subset; + struct stat filestats; + + Nregions = *nregions; + + Nsubset = 0; + NSUBSET = 100; + ALLOCATE (subset, GSCRegion, NSUBSET); + + /* match the basename against the GSCRegion file names */ + for (i = 0; i < Nregions; i++) { + status = stat (regions[i].filename, &filestats); + if ((status == -1) && (errno == ENOENT)) continue; + /* give an error for other conditions? */ + + subset[Nsubset] = regions[i]; + Nsubset ++; + CHECK_REALLOCATE (subset, GSCRegion, NSUBSET, Nsubset, 100); + } + free (regions); + + *nregions = Nsubset; + return (subset); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match_glob.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match_glob.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregion_match_glob.c (revision 4773) @@ -0,0 +1,51 @@ +# include "addstar.h" +# include + +GSCRegion *gregion_match (int *nregions) { + + int i; + int Ngsc, Nregions; + glob_t pglob; + + GSCRegion *gsc, *regions; + + gsc = LoadRegions (&Ngsc); + + sprintf (globline, "%s/[nsNS]????/*.cpt", CATDIR); + pglob.gl_offs = 0; + glob (globline, 0, NULL, &pglob); + + Nregions = 0; + ALLOCATE (regions, GSCRegion, MAX (1, pglob.gl_pathc)); + + /* match the basename against the GSCRegion file names */ + for (i = 0; i < pglob.gl_pathc; i++) { + found = FALSE; + for (j = 0; !found && (j < Ngsc); j++) { + if (strcasecmp (pglob.gl_pathv, gsc[j].filename)) continue; + regions[Nregions] = gsc[j]; + found = TRUE; + } + } + *nregions = Nregions; + return (regions); + + /* should free things up here */ +} + +/* given a ra,dec bounded region, find all region files which overlap it + + - split region into two if overlapping 0,360 boundary + +*/ + +/* the goal here is to return a list of the existing region files + this means either: + + grabbing all files that are on disk and converting + them to GSCRegion files + + trying all GSCRegion files to see if they exist + (and are populated?) + +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregion_patch.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregion_patch.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregion_patch.c (revision 4773) @@ -0,0 +1,43 @@ +# include "addstar.h" + +GSCRegion *gregion_patch (GSCRegion *patch, int *nregions) { + + int i, Ns, Ne, Ngsc, Nregions; + GSCRegion *gsc, *regions; + + gsc = LoadRegions (&Ngsc); + + /* find lower and upper boundaries */ + Ns = FindRegionDecBandStart (gsc, Ngsc, patch[0].DEC[0]); + Ne = FindRegionDecBandStop (gsc, Ngsc, patch[0].DEC[1]); + if (Ns < 0) Shutdown ("invalid patch lower bound %f", patch[0].DEC[0]); + if (Ne < 0) Shutdown ("invalid patch upper bound %f", patch[0].DEC[1]); + + Nregions = 0; + ALLOCATE (regions, GSCRegion, Ne - Ns); + + /* find all regions in range within ra range */ + for (i = Ns; i < Ne; i++) { + if (gsc[i].RA[0] >= patch[0].RA[1]) continue; + if (gsc[i].RA[1] < patch[0].RA[0]) continue; + if (gsc[i].DEC[0] >= patch[0].DEC[1]) continue; + if (gsc[i].DEC[1] < patch[0].DEC[0]) continue; + /* add to the list */ + regions[Nregions] = gsc[i]; + Nregions ++; + } + + if (VERBOSE) fprintf (stderr, "Nregions in patch: %d, Ne: %d, Ns: %d\n", Nregions, Ne, Ns); + if (Nregions > (Ne - Ns)) { + fprintf (stderr, "ERROR: too many regions found\n"); + exit (1); + } + + *nregions = Nregions; + return (regions); +} + +/* given a ra,dec bounded region, find all region files which overlap it + - split region into two if overlapping 0,360 boundary +*/ + Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregion_star.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregion_star.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregion_star.c (revision 4773) @@ -0,0 +1,47 @@ +# include "addstar.h" + +GSCRegion *gregion_stars (Stars *stars, int Nstars, int *nregions) { + + int i, j, Ngsc, Nregions, NREGIONS; + GSCRegion *gsc, *regions; + + gsc = LoadRegions (&Ngsc); + + Nregions = 0; + NREGIONS = 100; + ALLOCATE (regions, GSCRegion, NREGIONS); + + sort_stars (stars, Nstars); /* sort by RA */ + for (i = 0; i < Nstars; i++) { + if (stars[i].found) continue; + FindRegionByPoint (®ions[Nregions], stars[i].R, stars[i].D, gsc, Ngsc); + stars[i].found = TRUE; + /* search forward for all contained stars */ + for (j = i; j < Nstars; j++) { + if (stars[j].R >= regions[Nregions].RA[1]) break; + if (stars[j].D < regions[Nregions].DEC[0]) continue; + if (stars[j].D >= regions[Nregions].DEC[1]) continue; + stars[j].found = TRUE; + } + Nregions ++; + CHECK_REALLOCATE (regions, GSCRegion, Nregions, NREGIONS, 100); + + } + *nregions = Nregions; + return (regions); +} + +/* given a list of stars, find all region files which contain them + + - sort by ra? + - loop over stars + - find region file which contains star + - go backwards until ra < RA[0] + - mark all stars in this region file + - go forwards until ra > RA[1] + - mark all stars in this region file + + - use DEC band information? + - would like to minimize the number of disk reads + +*/ Index: /tags/addstar-1-2/Ohana/src/addstar/src/gregions.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gregions.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gregions.c (revision 4773) @@ -0,0 +1,53 @@ +# include "addstar.h" + +GSCRegion *gregions (Image *image, int *Nregions) { + + GSCRegion *region; + FILE *f; + double x, y; + double dec, ra; + int i, j, done, nregion, NREGION; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC region file %s\n", GSCFILE); + exit (1); + } + + /* find regions at image corners */ + NREGION = 10; + ALLOCATE (region, GSCRegion, NREGION); + nregion = 0; + + /* look for new regions on grid across image */ + for (x = 0.0; x <= 1.0; x+=0.25) { + for (y = 0.0; y <= 1.0; y+=0.25) { + XY_to_RD (&ra, &dec, image[0].NX*(1.1*x - 0.05), image[0].NY*(1.1*y - 0.05), &image[0].coords); + aregion (®ion[nregion], f, ra, dec); + done = FALSE; + for (j = 0; (j < nregion) && !done; j++) { + if (!strcmp (region[nregion].filename, region[j].filename)) { + nregion --; + done = TRUE; + } + } + nregion ++; + if (nregion == NREGION) { + NREGION += 10; + REALLOCATE (region, GSCRegion, NREGION); + } + } + } + + if (VERBOSE) { + fprintf (stderr, "found %d region files:\n", nregion); + for (i = 0; i < nregion; i++) { + fprintf (stderr, " %d %s\n", i, region[i].filename); + } + } + *Nregions = nregion; + + fclose (f); + return (region); + +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/gstars.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gstars.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gstars.c (revision 4773) @@ -0,0 +1,207 @@ +# include "addstar.h" + +Stars *gstars (char *file, int *NSTARS, Image *image) { + + FILE *f; + int j, N, Nbytes; + int itmp, hour, min; + char *c, photcode[64], line[80]; + double tmp, sec, dMs; + Stars *stars, *rdstars; + Header header; + + if (!fits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't read header for %s\n", file); + exit (1); + } + + /* find image rootname */ + c = strrchr (file, 0x2f); + if (c == (char *) NULL) { + strcpy (image[0].name, file); + } else { + strcpy (image[0].name, (c+1)); + } + + /* get astrometry information */ + if (!GetCoords (&image[0].coords, &header)) { + fprintf (stderr, "ERROR: no astrometric solution in header\n"); + exit (1); + } + itmp = 0; + fits_scan (&header, "NASTRO", "%d", 1, &itmp); + if ((itmp == 0) && !ACCEPT_ASTROM) { + fprintf (stderr, "ERROR: bad astrometric solution in header\n"); + exit (1); + } + if (!strcmp (&image[0].coords.ctype[4], "-WRP")) { + if (MOSAIC == NULL) { + fprintf (stderr, "ERROR: no mosaic for WRP image (use -mosaic)\n"); + exit (1); + } + RegisterMosaic (MOSAIC); + } else { + /* force image to lie in 0-360 range */ + while (image[0].coords.crval1 < 0) image[0].coords.crval1 += 360.0; + while (image[0].coords.crval1 > 360.0) image[0].coords.crval1 -= 360.0; + } + + /* get ST (used for airmass calculation) */ + fits_scan (&header, STKeyword, "%s", 1, line); + /* remove ':' characters */ + for (c = strchr (line, ':'); c != (char *) NULL; c = strchr (line, ':')) { *c = ' '; } + sscanf (line, "%d %d %lf", &hour, &min, &sec); + SiderealTime = hour + min/60.0 + sec/3600.0; + + /* CERROR in data file is in arcsec */ + if (!fits_scan (&header, "CERROR", "%lf", 1, &tmp)) tmp = 1.0; + image[0].cerror = tmp * 50.0; + fits_scan (&header, "NAXIS1", "%hd", 1, &image[0].NX); + fits_scan (&header, "NAXIS2", "%hd", 1, &image[0].NY); + + /* get photcode from header */ + if (thiscode == NULL) { + fits_scan (&header, "PHOTCODE", "%s", 1, photcode); + thiscode = GetPhotcodebyName (photcode); + } + if (thiscode == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", photcode); + exit (1); + } + image[0].source = thiscode[0].code; + + image[0].NX -= XOVERSCAN; + image[0].NY -= YOVERSCAN; + fits_scan (&header, ExptimeKeyword, "%lf", 1, &tmp); + image[0].exptime = tmp; + + tmp = 0; + /* fits_scan (&header, "APMIFIT", "%lf", 1, &tmp); */ + image[0].apmifit = tmp * 1000.0; + + tmp = 0; + /* fits_scan (&header, "dAPMIFIT", "%lf", 1, &tmp); */ + image[0].dapmifit = tmp * 1000.0; + + tmp = 0; + fits_scan (&header, "FLIMIT", "%lf", 1, &tmp); + image[0].detection_limit = tmp * 10.0; + + tmp = 0; + fits_scan (&header, "FSATUR", "%lf", 1, &tmp); + image[0].saturation_limit = tmp * 10.0; + + tmp = 0; + fits_scan (&header, "FWHM_X", "%lf", 1, &tmp); + image[0].fwhm_x = tmp * 25.0 * image[0].coords.cdelt1 * 3600.0; + + tmp = 0; + fits_scan (&header, "FWHM_Y", "%lf", 1, &tmp); + image[0].fwhm_y = tmp * 25.0 * image[0].coords.cdelt1 * 3600.0; + + if (!fits_scan (&header, "TZERO", "%d", 1, &image[0].tzero)) { + image[0].tzero = parse_time (&header); + } + + tmp = 0; + fits_scan (&header, "TRATE", "%lf", 1, &tmp); + image[0].trate = 10000 * tmp; + + tmp = 0; + fits_scan (&header, AirmassKeyword, "%lf", 1, &tmp); + image[0].secz = MIN (NO_MAG, 1000*tmp); + + if (!fits_scan (&header, CCDNumKeyword, "%d", 1, &itmp)) { + image[0].ccdnum = 0xff; + } else { + image[0].ccdnum = itmp; + } + + fits_scan (&header, "ZERO_PT", "%lf", 1, &ZeroPt); + + /* secz is in units milli-airmass */ + image[0].Mcal = 0.0; + image[0].Xm = NO_MAG; + image[0].code = 0; + bzero (image[0].dummy, sizeof(image[0].dummy)); + + /* find expected number of stars */ + if (!fits_scan (&header, "NSTARS", "%d", 1, &image[0].nstar)) { + fprintf (stderr, "ERROR: can't get NSTARS from header\n"); + exit (1); + } + + /* 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; + int sky; + + sky = 0; + fits_scan (&header, "HISTORY", "%S", 10, line); + p = strstr (line, "|I="); + if (p != (char *) NULL) { + p += 3; + sky = atoi (p); + } + image[0].Myyyy = MIN (sky - 0x8000, 0xffff); + } + + /* read in data section */ + f = fopen (file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't read data from %s\n", file); + exit (1); + } + fseek (f, header.size, SEEK_SET); + + if (FITS_INPUT) { + Nbytes = fits_matrix_size (&header); + fseek (f, Nbytes, SEEK_CUR); + rdstars = rfits (f, &image[0].nstar); + } else { + rdstars = rtext (f, &image[0].nstar); + } + fclose (f); + + /* modify resulting star list */ + ALLOCATE (stars, Stars, image[0].nstar); + for (N = j = 0; j < image[0].nstar; j++) { + /* allow for some dynamic filtering of star list */ + if (rdstars[j].dM > 1000.0 / SNLIMIT) continue; + if (XMAX && (rdstars[j].X > XMAX)) continue; + if (XMIN && (rdstars[j].X < XMIN)) continue; + if (YMAX && (rdstars[j].Y > YMAX)) continue; + if (YMIN && (rdstars[j].Y < YMIN)) continue; + stars[N] = rdstars[j]; + + XY_to_RD (&stars[N].R, &stars[N].D, stars[N].X, stars[N].Y, &image[0].coords); + while (stars[N].R < 0.0) stars[N].R += 360.0; + while (stars[N].R >= 360.0) stars[N].R -= 360.0; + stars[N].found = -1; + stars[N].code = thiscode[0].code; + + if (SKYPROBE) { + dMs = get_subpix (stars[N].X, stars[N].Y); + stars[N].M -= dMs; + stars[N].Mgal -= dMs; + stars[N].Map -= dMs; + dMs = 1000.0*scat_subpix (stars[N].X, stars[N].Y); + stars[N].dM = hypot (stars[N].dM, dMs); + } + N ++; + } + image[0].nstar = N; + REALLOCATE (stars, Stars, image[0].nstar); + free (rdstars); + + if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar); + *NSTARS = image[0].nstar; + + return (stars); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/gztest.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/gztest.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/gztest.c (revision 4773) @@ -0,0 +1,99 @@ +# include "ohana.h" +# include "zlib.h" +# define DTIME(A,B) ((A.tv_sec - B.tv_sec) + 1e-6*(A.tv_usec - B.tv_usec)) + +main (int argc, char **argv) { + + int N, Nread, Nseek, Nback; + char *filename, *buffer; + gzFile gf; + FILE *f; + struct timeval start, stop; + + if (argc != 6) { + fprintf (stderr, "USAGE: gztest (mode) (file) (Nseek) (Nread) (Nback)\n"); + exit (2); + } + + filename = argv[2]; + Nseek = atoi (argv[3]); + Nread = atoi (argv[4]); + Nback = -1 * atoi (argv[5]); + ALLOCATE (buffer, char, Nread); + + if (!strcmp (argv[1], "gz")) { + gf = gzopen (filename, "rb"); + if (gf == NULL) { + fprintf (stderr, "can't read data file: %s", filename); + exit (1); + } + + gettimeofday (&start, NULL); + N = gzseek (gf, Nseek, SEEK_SET); + gettimeofday (&stop, NULL); + if (N != Nseek) { + fprintf (stderr, "error seeking\n"); + exit (1); + } + fprintf (stdout, "seek: %f\n", DTIME (stop, start)); + + gettimeofday (&start, NULL); + N = gzread (gf, buffer, Nread); + if (N != Nread) { + fprintf (stderr, "error reading\n"); + exit (1); + } + gettimeofday (&stop, NULL); + fprintf (stdout, "read: %f\n", DTIME (stop, start)); + + gettimeofday (&start, NULL); + N = gzseek (gf, Nback, SEEK_CUR); + gettimeofday (&stop, NULL); + if (N == -1) { + fprintf (stderr, "error seeking\n"); + exit (1); + } + fprintf (stdout, "back: %f\n", DTIME (stop, start)); + exit (0); + } + + if (!strcmp (argv[1], "raw")) { + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "can't read data file: %s", filename); + exit (1); + } + + gettimeofday (&start, NULL); + N = fseek (f, Nseek, SEEK_SET); + gettimeofday (&stop, NULL); + if (N) { + fprintf (stderr, "error seeking\n"); + exit (1); + } + fprintf (stdout, "seek: %f\n", DTIME (stop, start)); + + gettimeofday (&start, NULL); + N = fread (buffer, 1, Nread, f); + gettimeofday (&stop, NULL); + if (N != Nread) { + fprintf (stderr, "error reading\n"); + exit (1); + } + fprintf (stdout, "read: %f\n", DTIME (stop, start)); + + gettimeofday (&start, NULL); + N = fseek (f, Nback, SEEK_CUR); + gettimeofday (&stop, NULL); + if (N) { + fprintf (stderr, "error seeking\n"); + exit (1); + } + fprintf (stdout, "back: %f\n", DTIME (stop, start)); + exit (0); + } + + + fprintf (stderr, "unknown mode\n"); + exit (1); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/image-db.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/image-db.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/image-db.c (revision 4773) @@ -0,0 +1,59 @@ +# include "addstar.h" + +# define LOCK LCK_XCLD + +/* 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 (char *format, ...) { + va_list argp; + char *formatplus; + + ALLOCATE (formatplus, char, strlen(format)); + strcpy (formatplus, format); + strcat (formatplus, "\n"); + + va_start (argp, format); + vfprintf (stderr, formatplus, argp); + free (formatplus); + va_end (argp); + + SetProtect (TRUE); + fclearlockfile (ImageCat, f, LOCK, &dbstate); + fprintf (stderr, "ERROR: addstar halted\n"); + exit (1); +} + +void lock_image_db () { + /* lock the image catalog */ + check_permissions (ImageCat); + f = fsetlockfile (ImageCat, 3600.0, LOCK, &dbstate); + if (f == NULL) { + fprintf (stderr, "ERROR: can't lock image catalog\n"); + exit (1); + } + fseek (f, 0, SEEK_SET); +} + +void unlock_image_db (Image *image) { + + mode_t mode; + + /* protect wimage from interrupt signals */ + if (MODE == M_IMAGE) { + SetProtect (TRUE); + wimage (image); + SetProtect (FALSE); + } + + fclearlockfile (ImageCat, f, LOCK, &dbstate); + mode = S_IRUSR | S_IWUSR | S_IRGRP | S_IWGRP | S_IROTH | S_IWOTH; + chmod (ImageCat, mode); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/in_image.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/in_image.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/in_image.c (revision 4773) @@ -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: /tags/addstar-1-2/Ohana/src/addstar/src/load_pt_catalog.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/load_pt_catalog.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/load_pt_catalog.c (revision 4773) @@ -0,0 +1,26 @@ +# include "addstar.h" + +int load_pt_catalog (Catalog *catalog, GSCRegion *region) { + + check_permissions (region[0].filename); + + /* 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[0].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: /tags/addstar-1-2/Ohana/src/addstar/src/load_subpix.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/load_subpix.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/load_subpix.c (revision 4773) @@ -0,0 +1,51 @@ +# include "addstar.h" + +typedef struct { + double Amp; + double Phase; + double dM; +} SubPixFix; + +static int Nsubpix; +static SubPixFix *Subpix; + +void load_subpix () { + + int i; + FILE *f; + + Nsubpix = 40; + ALLOCATE (Subpix, SubPixFix, Nsubpix); + + f = fopen (SubpixDatafile, "r"); + if (f == NULL) Shutdown ("can't load subpix datafile %s", SubpixDatafile); + + for (i = 0; i < Nsubpix; i++) { + fscanf (f, "%*s %*s %lf %lf %lf %*s\n", + &Subpix[i].Amp, &Subpix[i].Phase, &Subpix[i].dM); + } + fclose (f); + +} + +double get_subpix (double x, double y) { + + int bin; + double dy, dM; + + dy = y - (int)(y); + bin = 5 * (int)(x/100) + (int)(y/100); + dM = Subpix[bin].Amp*sin(2*3.14159*dy + Subpix[bin].Phase); + return (dM); +} + +double scat_subpix (double x, double y) { + + int bin; + double dy, dM; + + dy = y - (int)(y); + bin = 5 * (int)(x/100) + (int)(y/100); + dM = Subpix[bin].dM; + return (dM); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/make_backup.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/make_backup.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/make_backup.c (revision 4773) @@ -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: /tags/addstar-1-2/Ohana/src/addstar/src/match_refstars.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/match_refstars.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/match_refstars.c (revision 4773) @@ -0,0 +1,68 @@ +# include "addstar.h" + +int match_refstars (Stars *stars, int Nstars) { + + int i, j, Nm, Nf, Nmissed, Nsubset; + int *subref; + Stars *subset; + GSCRegion region; + Catalog catalog; + + /* put stars in the appropriate db file */ + for (i = 0; i < Nstars; i++) { + if (stars[i].found != -1) continue; + gregion_star (&stars[i], ®ion); + check_permissions (region.filename); + make_backup (region.filename); + + catalog.filename = region.filename; /* don't free region before catalog! */ + fprintf (stderr, "adding to %s [star %d]\n", region.filename, i); + + switch (lock_catalog (&catalog, LCK_XCLD)) { + case 0: + fprintf (stderr, "ERROR: can't lock file %s\n", catalog.filename); + exit (1); + case 1: + gcatalog (&catalog); /* load from disk */ + break; + case 2: + mkcatalog (®ion, &catalog); /* fills in new header info */ + break; + } + subset = find_subset (®ion, stars, Nstars, &subref, &Nsubset); + if (VERBOSE) fprintf (stderr, "%d stars in subset for %s\n", Nsubset, region.filename); + find_matches_refstars (®ion, subset, Nsubset, &catalog); + + /* protect wcatalog from interrupt signals */ + SetProtect (TRUE); + wcatalog (&catalog); + SetProtect (FALSE); + unlock_catalog (&catalog); + + Nm = 0; Nf = 0; + for (j = 0; j < Nsubset; j++) { + stars[subref[j]].found = subset[j].found; + if (subset[j].found == -2) Nm ++; + if (subset[j].found >= 0) Nf ++; + } + /* found: + -1 - not yet tried + -2 - found on previous catalog + -3 - not found on previous catalog + +N - found elsewhere + */ + if (VERBOSE) fprintf (stderr, "found %d stars, skipped %d stars\n", Nf, Nm); + free (subset); + free (subref); + } + + for (Nmissed = i = 0; i < Nstars; i++) { + if (stars[i].found == -1) { + fprintf (stderr, "%d %f %f %f %f\n", i, stars[i].R, stars[i].D, stars[i].M, stars[i].dM); + Nmissed ++; + } + } + if (Nmissed) fprintf (stderr, "WARNING: %d stars in image were missed!\n", Nmissed); + + return (TRUE); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/mkacc-2mass.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/mkacc-2mass.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/mkacc-2mass.c (revision 4773) @@ -0,0 +1,68 @@ +# include "addstar.h" +# include "2mass.h" + +int main (int argc, char **argv) { + + int i, N, Nrefcat; + Stars *refcat; + GSCRegion *regions, accregion; + + ConfigInit (&argc, argv); + + /* 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); + } + + if (argc != 2) { + fprintf (stderr, "USAGE: mkacc.2mass (file)\n"); + exit (2); + } + + NAMED_PHOTCODE (TM_J, "2MASS_J"); + NAMED_PHOTCODE (TM_H, "2MASS_H"); + NAMED_PHOTCODE (TM_K, "2MASS_K"); + if (thiscode == NULL) Shutdown ("photcode not specified"); + if (thiscode[0].code == TM_J) goto good_code; + if (thiscode[0].code == TM_H) goto good_code; + if (thiscode[0].code == TM_K) goto good_code; + Shutdown ("2MASS photcode not specified"); + +good_code: + UserPatch.RA[0] = 0; + UserPatch.RA[1] = 360; + UserPatch.DEC[0] = -90; + UserPatch.DEC[1] = +90; + + ALLOCATE (regions, GSCRegion, 1); + strcpy (regions[0].filename, argv[1]); + regions[0].RA[0] = 0; + regions[0].RA[1] = 360; + regions[0].DEC[0] = -90; + regions[0].DEC[1] = +90; + + refcat = get2mass_AS_data (regions, &UserPatch, &Nrefcat); + + /* find upper and lower file limits in RA and DEC */ + + strcpy (accregion.filename, regions[0].filename); + accregion.RA[0] = 360; + accregion.RA[1] = 0; + accregion.DEC[0] = +90; + accregion.DEC[1] = -90; + + for (i = 0; i < Nrefcat; i++) { + accregion.RA[0] = MIN (refcat[i].R, accregion.RA[0]); + accregion.RA[1] = MAX (refcat[i].R, accregion.RA[1]); + accregion.DEC[0] = MIN (refcat[i].D, accregion.DEC[0]); + accregion.DEC[1] = MAX (refcat[i].D, accregion.DEC[1]); + } + + fprintf (stderr, "%s %10.6f %10.6f %10.6f %10.6f %d\n", + accregion.filename, accregion.RA[0]/15.0, accregion.RA[1]/15.0, accregion.DEC[0], accregion.DEC[1], Nrefcat); + + exit (0); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/mkcatalog.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/mkcatalog.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/mkcatalog.c (revision 4773) @@ -0,0 +1,37 @@ +# include "addstar.h" + +void mkcatalog (GSCRegion *region, Catalog *catalog) { + + char filename[64], line[64]; + struct tm *local; + struct timeval now; + + if (VERBOSE) fprintf (stderr, "new catalog file: %s\n", region[0].filename); + if (!fits_read_header (CatTemplate, &catalog[0].header)) { + fprintf (stderr, "ERROR: can't find template header %s\n", CatTemplate); + exit (1); + } + + /* assign header values for RA, DEC, DATE, etc */ + fits_modify (&catalog[0].header, "RA0", "%lf", 1, region[0].RA[0]); + fits_modify (&catalog[0].header, "DEC0", "%lf", 1, region[0].DEC[0]); + fits_modify (&catalog[0].header, "RA1", "%lf", 1, region[0].RA[1]); + fits_modify (&catalog[0].header, "DEC1", "%lf", 1, region[0].DEC[1]); + gettimeofday (&now, NULL); + local = localtime (&now.tv_sec); + sprintf (line, "%02d/%02d/%02d", local[0].tm_year, local[0].tm_mon + 1, local[0].tm_mday); + fits_modify (&catalog[0].header, "DATE", "%s", 1, line); + strcpy (filename, CatTemplate); + ALLOCATE (catalog[0].average, Average, 1); + ALLOCATE (catalog[0].measure, Measure, 1); + ALLOCATE (catalog[0].missing, Missing, 1); + catalog[0].Naverage = catalog[0].Nmeasure = catalog[0].Nmissing = 0; + + /* setup secondary filters to match photcodes */ + ALLOCATE (catalog[0].secfilt, SecFilt, 1); + catalog[0].Nsecfilt = GetPhotcodeNsecfilt (); + + /* Nsecfilt is number of filters. Number of entries in array is + Nsecfilt * Naverage. At this point, N entries == 0 */ + +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/opening_angle.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/opening_angle.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/opening_angle.c (revision 4773) @@ -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: /tags/addstar-1-2/Ohana/src/addstar/src/parse_time.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/parse_time.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/parse_time.c (revision 4773) @@ -0,0 +1,121 @@ +# include "addstar.h" + +int parse_time (Header *header) { + + double jd; + int Ny, Nf, mode; + int Nsec, hour, min, sec, year, month, day; + char *py, *pm, *pd, *c; + char line[256]; + + /* we want to find JD or MJD to get Nsec (seconds since 01/01/1970) */ + + /* try JD first */ + if (strcasecmp (JDKeyword, "NONE")) { + uppercase (JDKeyword); + fits_scan (header, JDKeyword, "%lf", 1, &jd); + Nsec = (jd - 2440587.5)*86400; + return (Nsec); + } + + /* try MJD next */ + if (strcasecmp (MJDKeyword, "NONE")) { + uppercase (MJDKeyword); + fits_scan (header, MJDKeyword, "%lf", 1, &jd); + Nsec = (jd - 40587.0)*86400; + return (Nsec); + } + + /* get UT and DATE */ + uppercase (UTKeyword); + fits_scan (header, UTKeyword, "%s", 1, line); + /* remove ':' characters */ + for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; } + sscanf (line, "%d %d %d", &hour, &min, &sec); + + /* parse mode line */ + uppercase (DateMode); + for (Ny = 0, c = strchr (DateMode, 'Y'); c != (char ) NULL; c = strchr (c + 1, 'Y'), Ny++); + if ((Ny != 2) && (Ny != 4)) { + fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode); + exit (1); + } + py = strchr (DateMode, 'Y'); + pm = strchr (DateMode, 'M'); + pd = strchr (DateMode, 'D'); + if ((py == (char *) NULL) || (pm == (char *) NULL) || (pd == (char *) NULL)) { + fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode); + exit (1); + } + if ((py > pm) && (py < pd)) { + fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode); + exit (1); + } + if ((py > pd) && (py < pm)) { + fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode); + exit (1); + } + mode = 0; + if ((py < pm) && (pm < pd)) { mode = 1; } /* yyyy-mm-dd */ + if ((py < pm) && (pm > pd)) { mode = 2; } /* yyyy-dd-mm */ + if ((py > pm) && (pm < pd)) { mode = 3; } /* mm-dd-yyyy */ + if ((py > pm) && (pm > pd)) { mode = 4; } /* dd-mm-yyyy */ + if (!mode) { + fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode); + exit (1); + } + + /* parse date entry */ + uppercase (DateKeyword); + fits_scan (header, DateKeyword, "%s", 1, line); + /* remove possible separators: ':', '/' '.', '-' */ + for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; } + for (c = strchr (line, 0x2f); c != (char *) NULL; c = strchr (line, 0x2f)) { *c = ' '; } + for (c = strchr (line, 0x2e); c != (char *) NULL; c = strchr (line, 0x2e)) { *c = ' '; } + for (c = strchr (line, 0x2d); c != (char *) NULL; c = strchr (line, 0x2d)) { *c = ' '; } + + Nf = 0; + switch (mode) { + case 1: + Nf = sscanf (line, "%d %d %d", &year, &month, &day); + break; + case 2: + Nf = sscanf (line, "%d %d %d", &year, &day, &month); + break; + case 3: + Nf = sscanf (line, "%d %d %d", &month, &day, &year); + break; + case 4: + Nf = sscanf (line, "%d %d %d", &day, &month, &year); + break; + } + if (Nf != 3) { + fprintf (stderr, "error in date entry (%s) or DATE-MODE format (%s)\n", line, DateMode); + exit (1); + } + + if (year > 1000) { + if (Ny == 2) { + fprintf (stderr, "warning: mode line claims 2 digit year, but 4 digit year found\n"); + } + } else { + if (Ny == 4) { + fprintf (stderr, "warning: mode line claims 4 digit year, but 2 digit year found\n"); + } + if (year < 50) year += 100; + year += 1900; + } + + /* this should probably use localtime */ + + /* convert yy.mm.dd hh.mm.ss to Nsec since 1970 (jd = 2440587.5) */ + /* note that in this section, tm_mon has range 1-12, unlike for gmtime () */ + jd = day - 32075 + (int)(1461*(year + 4800 + (int)(((month)-14)/12))/4) + + (int)(367*((month) - 2 - (int)(((month) - 14)/12)*12)/12) + - (int)(3*(int)((year + 4900 + (int)(((month) - 14)/12))/100)/4) - 0.5; + /* jd is the julian day of the whole day only not the time */ + Nsec = (jd - 2440587.5)*86400 + 3600.0*hour + min*60.0 + sec; + + return (Nsec); + +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/replace_match.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/replace_match.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/replace_match.c (revision 4773) @@ -0,0 +1,20 @@ +# include "addstar.h" + +int replace_match (Average *average, Measure *measure, Stars *star) { + + int i; + + if (!REPLACE) return (FALSE); + + /** in replace mode, search for entry and replace values M, dM, R, D */ + for (i = 0; i < average[0].Nm; i++) { + if (measure[i].source != thiscode[0].code) continue; + measure[i].dR = 360000.0*(average[0].R - star[0].R); + measure[i].dD = 360000.0*(average[0].D - star[0].D); + measure[i].M = MIN (1000*star[0].M, NO_MAG); + measure[i].dM = MIN (1000*star[0].dM, NO_ERR); + star[0].found = average[0].offset + i; + return (TRUE); + } + return (FALSE); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/rfits.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/rfits.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/rfits.c (revision 4773) @@ -0,0 +1,41 @@ +# include "addstar.h" + +Stars *rfits (FILE *f, int *nstars) { + + int i, Nstars; + Header theader; + FTable table; + Stars *stars; + SMPData *smpdata; + + /* if no stars, no table */ + if (*nstars == 0) return (NULL); + + /* init & load in table data */ + table.header = &theader; + if (!fits_fread_ftable (f, &table, "SMPFILE")) goto escape; + + smpdata = fits_table_get_SMPData (&table, &Nstars); + + ALLOCATE (stars, Stars, Nstars); + for (i = 0; i < Nstars; i++) { + stars[i].X = smpdata[i].X; + stars[i].Y = smpdata[i].Y; + stars[i].M = smpdata[i].M; + stars[i].dM = smpdata[i].dM; + stars[i].dophot = smpdata[i].dophot; + + stars[i].Mgal = smpdata[i].M; + stars[i].Map = smpdata[i].dM; + stars[i].fx = smpdata[i].fx; + stars[i].fy = smpdata[i].fy; + stars[i].df = smpdata[i].df; + } + *nstars = Nstars; + return (stars); + +escape: + fprintf (stderr, "error reading file\n"); + *nstars = 0; + return (NULL); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/rtext.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/rtext.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/rtext.c (revision 4773) @@ -0,0 +1,85 @@ +# include "addstar.h" +# define D_NSTARS 1000 +# define BYTES_STAR 66 +# define BLOCK 1000 + +Stars *rtext (FILE *f, int *nstars) { + + int j, N, Nextra, Ninstar, Nskip, Nbytes, nbytes; + int done; + char *buffer, *c, *c2; + double tmp; + Stars *stars; + + /* load in stars by blocks of 1000 */ + N = 0; + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR) + 1); + buffer[BLOCK*BYTES_STAR] = 0; + Nextra = 0; + + ALLOCATE (stars, Stars, *nstars); + + while (N < *nstars) { + /* load next data block */ + Nbytes = BYTES_STAR * BLOCK - Nextra; + nbytes = Fread (&buffer[Nextra], 1, Nbytes, f, "char"); + if (nbytes == 0) { + *nstars = N; + return (stars); + } + nbytes += Nextra; + + /* check line-by-line integrity */ + c = buffer; + done = FALSE; + while ((c < buffer + nbytes) && (!done)) { + for (c2 = c; *c2 == '\n'; c2++); + if (c2 > c) { /* extra return chars */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + if (VERBOSE) fprintf (stderr, "deleted %d extra return chars\n", Nskip); + } + c2 = strchr (c, '\n'); + if (c2 == (char *) NULL) { + done = TRUE; + continue; + } + c2++; + if ((c2 - c) != BYTES_STAR) { /* bad line, delete it */ + memmove (c, c2, (int)(buffer + nbytes - c2)); + Nskip = c2 - c; + nbytes -= Nskip; + bzero (buffer + nbytes, Nskip); + if (VERBOSE) fprintf (stderr, "deleted line, %d extra chars\n", Nskip); + } else { + c = c2; + } + } + + /* extract data for stars */ + Ninstar = nbytes / BYTES_STAR; + Nextra = nbytes % BYTES_STAR; + for (j = 0; (j < Ninstar) && (N < *nstars); j++, N++) { + dparse (&stars[N].X, 1, &buffer[j*BYTES_STAR]); + dparse (&stars[N].Y, 2, &buffer[j*BYTES_STAR]); + dparse (&stars[N].M, 3, &buffer[j*BYTES_STAR]); + + /* cmp files carry dM in millimags */ + dparse (&tmp, 4, &buffer[j*BYTES_STAR]); + stars[N].dM = 0.001*tmp; + + 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]); + dparse (&stars[N].fx, 9, &buffer[j*BYTES_STAR]); + dparse (&stars[N].fy, 10, &buffer[j*BYTES_STAR]); + dparse (&stars[N].df, 11, &buffer[j*BYTES_STAR]); + } + } + *nstars = N; + return (stars); +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/save_pt_catalog.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/save_pt_catalog.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/save_pt_catalog.c (revision 4773) @@ -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: /tags/addstar-1-2/Ohana/src/addstar/src/sort_lists.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/sort_lists.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/sort_lists.c (revision 4773) @@ -0,0 +1,186 @@ +# include "addstar.h" + +void sort_stars (Stars *stars, int N) { + + int l,j,ir,i; + Stars tmp; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tmp = stars[l]; + } + else { + tmp = stars[ir]; + stars[ir] = stars[0]; + if (--ir == 0) { + stars[0] = tmp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && stars[j].R < stars[j+1].R) j++; + if (tmp.R < stars[j].R) { + stars[i] = stars[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + stars[i] = tmp; + } +} + +void sort_regions (GSCRegion *region, int N) { + + int l,j,ir,i; + GSCRegion tmp; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tmp = region[l]; + } else { + tmp = region[ir]; + region[ir] = region[0]; + if (--ir == 0) { + region[0] = tmp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && region[j].DEC[0] < region[j+1].DEC[0]) j++; + if (tmp.DEC[0] < region[j].DEC[0]) { + region[i] = region[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + region[i] = tmp; + } +} + +void sort_lists (float *X, float *Y, int *S, int N) { + + int l,j,ir,i; + double tX, tY, tS; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tS = S[l]; + } else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tS = S[ir]; + S[ir] = S[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + S[0] = tS; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + S[i] = S[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + S[i] = tS; + } +} + +void fsort (float *X, int N) { + + int l,j,ir,i; + float tX; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + } else { + tX = X[ir]; + X[ir] = X[0]; + if (--ir == 0) { + X[0] = tX; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + } +} + +void fsort2 (float *X, float *Y, int N) { + + int l,j,ir,i; + float tX, tY; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + } + else { + tX = X[ir]; + tY = Y[ir]; + X[ir] = X[0]; + Y[ir] = Y[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + } +} Index: /tags/addstar-1-2/Ohana/src/addstar/src/update_coords.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/update_coords.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/update_coords.c (revision 4773) @@ -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) { + + int i, m, Npt; + double R, D, r, d; + double r2, d2, dR2, dD2; + + 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 (i = 0; i < average[0].Nm; i++) { + 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 (i = 0; i < average[0].Nm; i++) { + 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: /tags/addstar-1-2/Ohana/src/addstar/src/wcatalog.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/wcatalog.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/wcatalog.c (revision 4773) @@ -0,0 +1,35 @@ +# include "addstar.h" + +void wcatalog (Catalog *catalog) { + + /* we'll leave these for now, but they are LONEOS specific */ + fits_modify (&catalog[0].header, "MARKSTAR", "%t", 1, FALSE); + fits_modify (&catalog[0].header, "ADDUSNO", "%t", 1, FALSE); + fits_modify (&catalog[0].header, "MARKROCK", "%t", 1, FALSE); + + if (!save_catalog (catalog, VERBOSE)) { + fprintf (stderr, "ERROR: failure saving catalog\n"); + exit (1); + } + + /* free, initialize data structures */ + + if (catalog[0].Naverage) { + free (catalog[0].average); + catalog[0].Naverage = 0; + } + if (catalog[0].Nmeasure) { + free (catalog[0].measure); + catalog[0].Nmeasure = 0; + } + if (catalog[0].Nmissing) { + free (catalog[0].missing); + catalog[0].Nmissing = 0; + } + if (catalog[0].Nsecfilt) { + free (catalog[0].secfilt); + catalog[0].Nsecfilt = 0; + } + fits_free_header (&catalog[0].header); +} + Index: /tags/addstar-1-2/Ohana/src/addstar/src/wimage.c =================================================================== --- /tags/addstar-1-2/Ohana/src/addstar/src/wimage.c (revision 4773) +++ /tags/addstar-1-2/Ohana/src/addstar/src/wimage.c (revision 4773) @@ -0,0 +1,57 @@ +# include "addstar.h" + +/* add image to image db */ +void wimage (Image *image) { + + FILE *f; + int Nimages, status, offset, dbstate; + Header header; + + f = GetDB (&dbstate); + fseek (f, 0, SEEK_SET); + + if (dbstate == LCK_EMPTY) { + if (VERBOSE) fprintf (stderr, "can't find %s, creating a new one\n", ImageCat); + if (!fits_read_header (ImageTemplate, &header)) { + fprintf (stderr, "ERROR: can't find template header %s\n", ImageTemplate); + exit (1); + } + /* 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)) { + fprintf (stderr, "ERROR: can't read image header %s\n", ImageCat); + exit (1); + } + } + + Nimages = 0; + fits_scan (&header, "NIMAGES", "%d", 1, &Nimages); + Nimages ++; + fits_modify (&header, "NIMAGES", "%d", 1, Nimages); + + /* position to begining of file to write header */ + fseek (f, 0, SEEK_SET); + status = Fwrite (header.buffer, 1, header.size, f, "char"); + if (status != header.size) { + fprintf (stderr, "ERROR: failed writing data to image header\n"); + exit (1); + } + + /* position to end of data array */ + offset = (Nimages - 1)*sizeof(Image) + header.size; + fseek (f, offset, SEEK_SET); + + status = Fwrite (image, sizeof(Image), 1, f, "image"); + if (status != 1) { + fprintf (stderr, "ERROR: failed writing data to image catalog\n"); + exit (1); + } +} + +/* the image we add in this routine is always a new image. + We only need to do two things: + 1) append the data for the image to the end of the file + 2) update the NIMAGES field in the header (which means, read in header, re-write) + */