Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/Makefile =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/Makefile (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/Makefile (revision 21807) @@ -0,0 +1,147 @@ +default: mosastro +help: + @echo "make options: mosastro (default) mkstandards mkobs warptest" + +include ../../Configure + +HOME = $(ROOT)/src/mosastro +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) -ldvo -lFITS -lohana -lm +CFLAGS = $(INCS) +LFLAGS = $(LIBS) + +MOS = \ +$(SRC)/mosastro.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/LoadStars.$(ARCH).o \ +$(SRC)/match.$(ARCH).o \ +$(SRC)/mkheader.$(ARCH).o \ +$(SRC)/fakefield.$(ARCH).o \ +$(SRC)/greference.$(ARCH).o \ +$(SRC)/project.$(ARCH).o \ +$(SRC)/transforms.$(ARCH).o \ +$(SRC)/dump.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/getstone.$(ARCH).o \ +$(SRC)/get2mass.$(ARCH).o \ +$(SRC)/getgsc.$(ARCH).o \ +$(SRC)/getptolemy.$(ARCH).o \ +$(SRC)/gptolemy.$(ARCH).o \ +$(SRC)/gregions.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/wstars.$(ARCH).o \ +$(SRC)/field.$(ARCH).o \ +$(SRC)/chips.$(ARCH).o \ +$(SRC)/GetGradients.$(ARCH).o \ +$(SRC)/FitGradients.$(ARCH).o \ +$(SRC)/FitChips.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/output.$(ARCH).o \ +$(SRC)/wstars.$(ARCH).o \ +$(SRC)/gaussj.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/mkmosaic.$(ARCH).o \ +$(SRC)/GetScatter.$(ARCH).o \ +$(SRC)/clip.$(ARCH).o \ +$(SRC)/parse_time.$(ARCH).o \ +$(SRC)/testcoords.$(ARCH).o \ +$(SRC)/rfits.$(ARCH).o \ +$(SRC)/wfits.$(ARCH).o \ +$(SRC)/rtext.$(ARCH).o \ +$(SRC)/ConvertMatch.$(ARCH).o \ +$(SRC)/SaveResiduals.$(ARCH).o + +mosastro: $(BIN)/mosastro.$(ARCH) +$(BIN)/mosastro.$(ARCH) : $(MOS) + +STD = \ +$(SRC)/mkstandards.$(ARCH).o \ +$(SRC)/mkheader.$(ARCH).o \ +$(SRC)/random.$(ARCH).o \ +$(SRC)/wstars.$(ARCH).o \ +$(SRC)/wfits.$(ARCH).o + +mkstandards: $(BIN)/mkstandards.$(ARCH) +$(BIN)/mkstandards.$(ARCH) : $(STD) + +OBS = \ +$(SRC)/args_obs.$(ARCH).o \ +$(SRC)/mkobs.$(ARCH).o \ +$(SRC)/mkheader.$(ARCH).o \ +$(SRC)/fakefield.$(ARCH).o \ +$(SRC)/greference.$(ARCH).o \ +$(SRC)/project.$(ARCH).o \ +$(SRC)/transforms.$(ARCH).o \ +$(SRC)/dump.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/get2mass.$(ARCH).o \ +$(SRC)/getgsc.$(ARCH).o \ +$(SRC)/getptolemy.$(ARCH).o \ +$(SRC)/getstone.$(ARCH).o \ +$(SRC)/gptolemy.$(ARCH).o \ +$(SRC)/gregions.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/random.$(ARCH).o \ +$(SRC)/wstars.$(ARCH).o + +mkobs: $(BIN)/mkobs.$(ARCH) +$(BIN)/mkobs.$(ARCH) : $(OBS) + +warptest: $(BIN)/warptest.(ARCH) +$(BIN)/warptest.(ARCH) : $(SRC)/warptest.$(ARCH).o + +$(MOS) : $(INC)/mosastro.h +$(OBS) : $(INC)/mosastro.h +$(STD) : $(INC)/mosastro.h + +# INSTALL = mosastro mkstandards mkobs +INSTALL = mosastro + +# 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/ohana/mosastro-1-0/Ohana/src/mosastro/bin/.cvsignore =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/bin/.cvsignore (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/bin/.cvsignore (revision 21807) @@ -0,0 +1,1 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/ChangeLog.txt =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/ChangeLog.txt (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/ChangeLog.txt (revision 21807) @@ -0,0 +1,5 @@ + +2005.10.20 : mosastro-1-0 + + added support for new dvo load methods + minor changes to use new libohana (v1.5) / libfits (v1.4) Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/astrom.txt =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/astrom.txt (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/astrom.txt (revision 21807) @@ -0,0 +1,99 @@ +mosaic astrometry keywords + +terms to model: + +boresite + projection +Ro, Do, Theta, PS +CTYPE, CRVAL, CRPIX, CD1_1 + +distortion +A(i,j), D(i,j) +WLi_j +WMi_j + +chips +WXi_j +WYi_j + +-- these are not general since they don't allow for more than 2 dimensions +-- these have the advantage of specifying which terms are available (WX0_2 is unambiguously X term for L^0M^2) + +- load the astrometry from a table (not very efficient unless all chip + terms are in the same table, which is an option) + +X = \Sum (WX_i_j) L^i M^j +Y = \Sum (WY_i_j) L^i M^j + +L = \Sum (WL_i_j) P^i Q^j +M = \Sum (WM_i_j) P^i Q^j + +P,Q = proj (R,D)[CTYPE, CRVAL, CRPIX, CDi_j] + + +Implementing the mosaic astrometry information in Elixir systems vs IPP. + +Elixir uses libohana/coordops to handle all coordinate +transformations. The coordinate structure (Coords) is part of the +DVO/Images.h structure. I cannot add an additional set of entries to +that structure without breaking old DVO or other elixir code. It is +also not strictly relevant since the Coords structure is used to +define the transformation of a single chip (and is sufficient to carry +those terms up to 3rd order). + +For elixir, we will need a two stage process, with one stage +transforming the celestial coords to the focal plane, with PLY terms +and a second to represent the chip to FP transform, with a varient on +PLY terms. (we can't use PLY since those use a projection, while we +want a transform that only corrsponds to an Nth order polynomial onto +a plane. We should call this transform WRP, with ctype PLN-WRP and +terms: + +crval1,2 = L,M coordinates for reference pixel +crpix1,2 = X,Y coordinates for reference pixel +cdelt1,2 = plate scale (should be 1,1 by default) +etc. + +A third order warp couple to a third order distortion leads to 6th +order terms in general. + +-- + +My plan is now to allow the dvo image table to represent chips in a mosaic +with either the current (per-chip) entry or an entry per chip coupled +to an entry per mosaic. The chip image entries would have astrometric +terms of up-to third order which transform from chip to focal plane, +coupled in turn to a mosaic entry in the table giving astrometry from +focal plane to the sky. The chip-to-fp entries would be identified +with astrometry CTYPE value of "-WRP" and the fp-to-sky would have +entries with CTYPE value of "-MOS". + +- does this information go in the header of each *.smp file? +- does this require us to package all chips in a single MEF smp +- package? +- mosastro would create the astrometry terms (per-chip and +- per-mosaic), where does it store them? +- minimal impact on addstar? + +- what are the implications for dvo? + - image I/O functions need to know about the mosaic + +- how do we match a chip to a mosaic? + chip is unique based on time and photcode. + - we could require the mosaic photcode entries to have values of + PRI/SEC + +: mosastro should construct a PHU element 654321oXX.smp which contains + the mosaic astrometry terms + +: alternatively: each .smp file could have a collection of header + keywords which specify the mosaic astrometry terms. + +: alternatively again, the + +: photcode for the mosaic entry is CAMERA.FILTER.XX +: addstar would have to have access to this file for each chip +: addstar should not create a new mosaic entry + +: how do we handle GetCoords (which needs to construct TWO coords + structures for output in the case of a mosaic?) + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/notes.txt =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/notes.txt (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/notes.txt (revision 21807) @@ -0,0 +1,59 @@ + +- need to update header astrometry + - determine higher-order polynomial terms from fit? + - insert all distortion & warping terms into header? + +- better input parameters: + - use camera to define chips we expect to find + - camera -> default warp + - add dispersion stats to header in output (CERROR, CPRECISE) + +- allow chips to have full linear terms (after warp) + +- load 2mass, USNO SA data + +----------------------------------------- +- load *.smp / *.cmp : + chip[i].pt[j].x + chip[i].pt[j].y + chip[i].coords + (global Ro, Do = RA_DEG, DEC_DEG) + +- load usno in region + stars[i].r, stars[i].d + +- match usno:smp: + chip[i].pt[j].r, chip[i].pt[j].d + +- project ra,dec to p,q (arcsec on sky relative to Ro, Do) + chip[i].pt[j].p, chip[i].pt[j].q + +- init chip paramters: + chip[i].Px, chip[i].Py - parities + chip[i].Po, chip[i].Qo, chip[i].To + + * p' = Px*(p - Po), q' = Py*(q - Qo) + * x = cos(To) p' + sin(To) q' + * y = cos(To) p' + sin(To) q' + +- other parameter: + * field.Ro, field.Do - not fitted + field.A[0] - field.A[5] + field.D[0] - field.D[5] + + + (r,d) -> field.coords -> (p,q) [arcsec from field center] + (p,q) -> field.A,D -> (L,M) [pixels from field center] + (L,M) -> Px,Py,Po,Qo -> (P,Q) [pixels from chip corner] + (P,Q) -> To -> (x,y) [pixels on chip] + + L = A[0] + A[1]*p + A[2]*q + A[3]*p^2 + A[4]*q^2 + A[5]*p*q + M = D[0] + D[1]*p + D[2]*q + D[3]*p^2 + D[4]*q^2 + D[5]*p*q + + p' = Px*(L - Po) + q' = Py*(M - Qo) + + x = cos(To) p' + sin(To) q' + y = cos(To) p' + sin(To) q' + + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/update.txt =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/update.txt (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/doc/update.txt (revision 21807) @@ -0,0 +1,8 @@ + +mosastro USAGE: + + - mosastro (inglob) (outext) [-dist (distort)] [-chips (chipset)] + + * given a collection of chip data, measure the complete model + parameters relative to reference. output is individual chip files + with a new extension (and updated header contents) Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/include/mosastro.h =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/include/mosastro.h (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/include/mosastro.h (revision 21807) @@ -0,0 +1,209 @@ +# include +# include + +typedef struct { + double R, D; /* Sky Coords - degrees */ + double P, Q; /* Tangent Plane - pixels */ + double L, M; /* Focal Plane - pixels */ + double X, Y; /* Chip Coords - pixels */ + double Mag, dMag; + int mask; +} StarData; + +typedef struct { + double RA[2], DEC[2]; + double Area; + char *name; +} CatStats; + +typedef struct { + char *file; + int Nstars; + int NX, NY; + StarData *stars; + Header header; + Header theader; + Matrix matrix; + + char *buffer; + int Nbuffer; + int FITS; + + int Nmatch; + StarData *raw, *ref; + + Coords coords; + Coords map; +} Chip; + +typedef struct { + double Ro, Do, PSx, PSy, To; + double Rmin, Rmax, Dmin, Dmax; + int Norder; + double **A; + double **D; + Coords project; + Coords distort; + int fit; +} Field; + +typedef struct { + double *dPdL; + double *dPdM; + double *dQdL; + double *dQdM; + double *Lo; + double *Mo; + int Npts; +} Gradients; + +# if (0) +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]; +} Stars; +# endif + +typedef struct { + double Rraw, Draw; /* Sky Coords - degrees */ + float Praw, Qraw; /* Tangent Plane - pixels */ + float Lraw, Mraw; /* Focal Plane - pixels */ + float Xraw, Yraw; /* Chip Coords - pixels */ + + double Rref, Dref; /* Sky Coords - degrees */ + float Pref, Qref; /* Tangent Plane - pixels */ + float Lref, Mref; /* Focal Plane - pixels */ + float Xref, Yref; /* Chip Coords - pixels */ + + float Mcat, dMcat; + float Minst, dMinst; + char mask; +} MatchData; + +int ChipOrder; +int Nchip; +Chip *chip; +Field field; +double Year; /** carried for precession - probably put this in chip data **/ +double RADIUS; /** raw / ref matching radius (pixels on Focal Plane) **/ +double SIGMA; +double SIGMA_LIM; +double IMAG_MIN; +double IMAG_MAX; +double INST_BRIGHT; +double ZERO_POINT; + +char StoneRegions[256]; +char CDROM[256]; +char REFCAT[256]; +char CATDIR[256]; +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char ExptimeKeyword[256]; +char DateKeyword[256]; +char DateMode[256]; +char UTKeyword[256]; +char MJDKeyword[256]; +char JDKeyword[256]; +char TWO_MASS_DIR[256]; +int VERBOSE; +int NO_CHIPS; +int SAVE_RESID; + +char GSCFILE[256]; +char GSC_DIR[256]; + +char *FIELD; +char *CHIPS; +char *OUTPUT; +char *DUMP; +char *FOCAL_PLANE; + +/*** mosastro prototypes ***/ +void ChipToFP PROTO((StarData *stars, int Nstars, Coords *coords)); +void ChipToSky PROTO((StarData *stars, int Nstars, Coords *coords)); +int ClipOnFP PROTO((double Nsigma)); +void ConfigInit PROTO((int *argc, char **argv)); +int ConvertMatch PROTO((MatchData *data, int size, int nitems)); +void FPtoChip PROTO((StarData *stars, int Nstars, Coords *coords)); +void FPtoTP PROTO((StarData *stars, int Nstars, Coords *coords)); +void FitChip PROTO((StarData *raw, StarData *ref, int Nmatch, Coords *coords)); +void FitChipLinear PROTO((StarData *raw, StarData *ref, int Nmatch, Coords *coords)); +void FitChipResid PROTO((StarData *raw, StarData *ref, int Nmatch, Coords *coords)); +void FitChips PROTO((int Norder)); +void FitGradients PROTO((Gradients *grad)); +void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); +Gradients *GetGradients PROTO(()); +double GetScatter PROTO((int *Nscatter, double *DL, double *DM, int bright)); +int LoadStars PROTO((int Nfile, char **file)); +void SaveResiduals PROTO((FILE *f, Header *header)); +void SkyToTP PROTO((StarData *stars, int Nstars, Coords *coords)); +void TPtoFP PROTO((StarData *stars, int Nstars, Coords *coords)); +void TPtoSky PROTO((StarData *stars, int Nstars, Coords *coords)); +void add_to_regions PROTO((CatStats *area)); +void area_of_region PROTO((CatStats *region)); +void args PROTO((int *argc, char **argv)); +int deproject_raw PROTO(()); +int deproject_stars PROTO(()); +int dump_grads PROTO((Gradients *grad, char *filename)); +int dump_match PROTO(()); +int dump_rawstars PROTO(()); +int dump_refcat PROTO((StarData *refcat, int Nrefcat)); +int dump_stars PROTO((FILE *f, StarData *stars, int Nstars)); +int fake_field PROTO((double RA, double DEC)); +void field_combine PROTO(()); +void field_stats PROTO(()); +int find_dec_bands PROTO((CatStats *area)); +void fit_add PROTO((double x1, double y1, double x2, double y2)); +void fit_apply_coords PROTO((Coords *coords)); +void fit_apply_grads PROTO((Coords *distort, Coords *project, int term)); +void fit_correct_grads PROTO((Gradients *in, Gradients *out, int term)); +void fit_eval PROTO(()); +void fit_free PROTO(()); +void fit_init PROTO((int order)); +int gaussj PROTO((double **a, int n, double **b, int m)); +StarData *gcatalog PROTO((char *filename, int *Nstars)); +StarData *get2mass PROTO((CatStats *catstats, int *NSTARS)); +StarData *getgsc PROTO((CatStats *catstats, int *NSTARS)); +StarData *getptolemy PROTO((CatStats *catstats, int *NSTARS)); +StarData *getstone PROTO((CatStats *input, int *nstars)); +StarData *getusno PROTO((CatStats *catstats, int *Nstars)); +StarData *gptolemy PROTO((char *filename, int *NSTARS)); +StarData *greference PROTO((int *Nrefcat)); +CatStats *gregions PROTO((CatStats *patch, int *nregion)); +int init_chips PROTO(()); +int init_field PROTO(()); +void init_regions PROTO(()); +int load_chips PROTO((char *filename)); +int load_field PROTO((char *filename)); +int load_ra_blocks PROTO((int Ndec, CatStats *area)); +int match PROTO((StarData *refcat, int Nrefcat)); +Header *mkheader PROTO((int Nx, int Ny, int Nstars, Coords *coords)); +Header *mkmosaic PROTO((int Nx, int Ny, int Nstars, Coords *coords)); +int mkpolyterm PROTO((int n, int m)); +int mkvector PROTO((int n, int m, int norder)); +void output PROTO((char *ext, char *phu)); +int parse_GSC_line PROTO((CatStats *tregion, char *line)); +e_time parse_time PROTO((Header *header)); +void print_help PROTO(()); +int project_ref PROTO(()); +int project_refcat PROTO((StarData *refcat, int Nrefcat)); +int project_stars PROTO(()); +int rfits PROTO((Chip *mychip)); +int rtext PROTO((Chip *mychip)); +void set_catalog PROTO((char *catdir)); +int sortthree PROTO((double *X, double *Y, int *Z, int N)); +void uppercase PROTO((char *string)); +void wchip PROTO((char *filename, Chip *data)); +void wfits PROTO((char *filename, SMPData *stars, int Nstars, Header *header)); +void wstars PROTO((char *filename, SMPData *stars, int Nstars, Header *header)); Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConfigInit.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConfigInit.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConfigInit.c (revision 21807) @@ -0,0 +1,64 @@ +# include "mosastro.h" + +void ConfigInit (int *argc, char **argv) { + + char *config, *file; + + VERBOSE = TRUE; + + /*** 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 (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + GetConfig (config, "USNO_CDROM", "%s", 0, CDROM); + ScanConfig (config, "2MASS_DIR", "%s", 0, TWO_MASS_DIR); + + GetConfig (config, "ASTRO_REFCAT", "%s", 0, REFCAT); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig (config, "GSCDIR", "%s", 0, GSC_DIR); + ScanConfig (config, "STONE_DIR", "%s", 0, StoneRegions); + + /* abstracted header 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); + + GetConfig (config, "RADIUS", "%lf", 0, &RADIUS); + GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + ScanConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + + ScanConfig (config, "INST_MAG_MIN", "%lf", 0, &IMAG_MIN); + ScanConfig (config, "INST_MAG_MAX", "%lf", 0, &IMAG_MAX); + ScanConfig (config, "INST_BRIGHT", "%lf", 0, &INST_BRIGHT); + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + free (config); + free (file); + + return; +} + +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/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertMatch.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertMatch.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertMatch.c (revision 21807) @@ -0,0 +1,70 @@ +# include "mosastro.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 MATCH_SIZE 104 + +int ConvertMatch (MatchData *data, int size, int nitems) { + + int i; + unsigned char *byte, tmp; + +# ifdef BYTE_SWAP + + if (size != MATCH_SIZE) { + fprintf (stderr, "mismatch in type sizes (MatchData) %d vs %d\n", size, MATCH_SIZE); + return (FALSE); + } + + byte = (char *) data; + for (i = 0; i < nitems; i++, byte += size) { + SWAP_DBLE (0); /* R */ + SWAP_DBLE (8); /* D */ + SWAP_WORD (16); /* P */ + SWAP_WORD (20); /* Q */ + SWAP_WORD (24); /* L */ + SWAP_WORD (28); /* M */ + SWAP_WORD (32); /* X */ + SWAP_WORD (36); /* Y */ + + SWAP_DBLE (40); /* R */ + SWAP_DBLE (40); /* D */ + SWAP_WORD (56); /* P */ + SWAP_WORD (60); /* Q */ + SWAP_WORD (64); /* L */ + SWAP_WORD (68); /* M */ + SWAP_WORD (72); /* X */ + SWAP_WORD (76); /* Y */ + + SWAP_WORD (80); /* Mcat */ + SWAP_WORD (84); /* dMcat */ + SWAP_WORD (88); /* Minst */ + SWAP_WORD (92); /* dMinst */ + } + return (TRUE); + +# else + return (TRUE); +# endif +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertStars.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertStars.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/ConvertStars.c (revision 21807) @@ -0,0 +1,103 @@ +# include "mosastro.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 +# define MATCH_SIZE 104 + +int ConvertStars (Stars *data, int size, int nitems) { + + int i; + unsigned char *byte, tmp; + +# ifdef BYTE_SWAP + + if (size != STARS_SIZE) { + fprintf (stderr, "mismatch in type sizes (Stars) %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 +} + +int ConvertMatch (MatchData *data, int size, int nitems) { + + int i; + unsigned char *byte, tmp; + +# ifdef BYTE_SWAP + + if (size != MATCH_SIZE) { + fprintf (stderr, "mismatch in type sizes (MatchData) %d vs %d\n", size, MATCH_SIZE); + return (FALSE); + } + + byte = (char *) data; + for (i = 0; i < nitems; i++, byte += size) { + SWAP_DBLE (0); /* R */ + SWAP_DBLE (8); /* D */ + SWAP_WORD (16); /* P */ + SWAP_WORD (20); /* Q */ + SWAP_WORD (24); /* L */ + SWAP_WORD (28); /* M */ + SWAP_WORD (32); /* X */ + SWAP_WORD (36); /* Y */ + + SWAP_DBLE (40); /* R */ + SWAP_DBLE (40); /* D */ + SWAP_WORD (56); /* P */ + SWAP_WORD (60); /* Q */ + SWAP_WORD (64); /* L */ + SWAP_WORD (68); /* M */ + SWAP_WORD (72); /* X */ + SWAP_WORD (76); /* Y */ + + SWAP_WORD (80); /* Mcat */ + SWAP_WORD (84); /* dMcat */ + SWAP_WORD (88); /* Minst */ + SWAP_WORD (92); /* dMinst */ + } + return (TRUE); + +# else + return (TRUE); +# endif +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitChips.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitChips.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitChips.c (revision 21807) @@ -0,0 +1,26 @@ +# include "mosastro.h" + +void FitChips (int Norder) { + + int i; + + for (i = 0; i < Nchip; i++) { + chip[i].map.Npolyterms = Norder; + FitChip (chip[i].raw, chip[i].ref, chip[i].Nmatch, &chip[i].map); + } + deproject_raw (); + project_ref (); +} + +void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { + + int i; + + fit_init (coords[0].Npolyterms); + for (i = 0; i < Nmatch; i++) { + if (raw[i].mask) continue; + fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M); + } + fit_eval (); + fit_apply_coords (coords); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitGradients.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitGradients.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/FitGradients.c (revision 21807) @@ -0,0 +1,47 @@ +# include "mosastro.h" + +/* fit dPdL, etc vs L,M */ +void FitGradients (Gradients *grad) { + + int i, Norder; + Gradients grfix; + + ALLOCATE (grfix.dPdL, double, grad[0].Npts); + ALLOCATE (grfix.dPdM, double, grad[0].Npts); + ALLOCATE (grfix.dQdL, double, grad[0].Npts); + ALLOCATE (grfix.dQdM, double, grad[0].Npts); + + ALLOCATE (grfix.Lo, double, grad[0].Npts); + ALLOCATE (grfix.Mo, double, grad[0].Npts); + grfix.Npts = grad[0].Npts; + + /* where do we set field.distort.Npolyterms? */ + Norder = 3; + + fit_init (Norder - 1); + for (i = 0; i < grad[0].Npts; i++) { + fit_add (grad[0].Lo[i], grad[0].Mo[i], grad[0].dPdL[i], grad[0].dPdM[i]); + } + fit_eval (); + fit_apply_grads (&field.distort, &field.project, 0); + fit_correct_grads (grad, &grfix, 0); + fit_free (); + + fit_init (Norder - 1); + for (i = 0; i < grad[0].Npts; i++) { + fit_add (grad[0].Lo[i], grad[0].Mo[i], grad[0].dQdL[i], grad[0].dQdM[i]); + } + fit_eval (); + fit_apply_grads (&field.distort, &field.project, 1); + fit_correct_grads (grad, &grfix, 1); + fit_free (); + + /* use new model of field to get TP & FP coords */ + deproject_raw (); + project_ref (); + + if ((DUMP != NULL) && !strcmp (DUMP, "grads")) { + dump_grads (&grfix, "gradfix.dat"); + exit (0); + } +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetGradients.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetGradients.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetGradients.c (revision 21807) @@ -0,0 +1,108 @@ +# include "mosastro.h" + +/* measure local gradient vs field position : dP, dQ = f(L,M) + gradient is measured in small boxes across the field */ +Gradients *GetGradients () { + + int i, j, nx, ny, Nval, Nx, Ny, Npts, DX, DY; + double **a, **b; + double L, M, dP, dQ; + double Xmin, Xmax, Ymin, Ymax; + + Gradients *grad; + + ALLOCATE (grad, Gradients, 1); + + /** this should not be hard-wired **/ + Nx = 2; + Ny = 2; + + ALLOCATE (grad[0].dPdL, double, Nchip*Nx*Ny); + ALLOCATE (grad[0].dPdM, double, Nchip*Nx*Ny); + ALLOCATE (grad[0].dQdL, double, Nchip*Nx*Ny); + ALLOCATE (grad[0].dQdM, double, Nchip*Nx*Ny); + + ALLOCATE (grad[0].Lo, double, Nchip*Nx*Ny); + ALLOCATE (grad[0].Mo, double, Nchip*Nx*Ny); + + ALLOCATE (a, double *, 3); + ALLOCATE (b, double *, 3); + for (i = 0; i < 3; i++) { + ALLOCATE (a[i], double, 3); + ALLOCATE (b[i], double, 2); + } + + Nval = 0; + for (i = 0; i < Nchip; i++) { + DX = chip[i].NX / Nx; + DY = chip[i].NY / Ny; + for (nx = 0; nx < Nx; nx++) { + for (ny = 0; ny < Ny; ny++) { + Xmin = nx*DX; + Xmax = Xmin + DX; + Ymin = ny*DY; + Ymax = Ymin + DY; + + for (j = 0; j < 3; j++) { + bzero (a[j], 3*sizeof(double)); + bzero (b[j], 2*sizeof(double)); + } + + /* find local gradient in limited box */ + for (j = 0; j < chip[i].Nmatch; j++) { + if (chip[i].raw[j].X < Xmin) continue; + if (chip[i].raw[j].X > Xmax) continue; + if (chip[i].raw[j].Y < Ymin) continue; + if (chip[i].raw[j].Y > Ymax) continue; + + L = chip[i].ref[j].L; + M = chip[i].ref[j].M; + dP = chip[i].ref[j].P - chip[i].raw[j].P; + dQ = chip[i].ref[j].Q - chip[i].raw[j].Q; + + a[0][0] += 1; + a[1][0] += L; + a[2][0] += M; + a[1][1] += L*L; + a[2][1] += L*M; + a[2][2] += M*M; + + b[0][0] += dP; + b[1][0] += dP*L; + b[2][0] += dP*M; + + b[0][1] += dQ; + b[1][1] += dQ*L; + b[2][1] += dQ*M; + + } + /* fitting dP = a[0][0] + a[1][0]*L + a[2][0]*M + dQ = a[0][1] + a[1][1]*L + a[2][1]*M + */ + a[0][1] = a[1][0]; + a[0][2] = a[2][0]; + a[1][2] = a[2][1]; + + Npts = a[0][0]; + grad[0].Lo[Nval] = a[1][0] / a[0][0]; + grad[0].Mo[Nval] = a[2][0] / a[0][0]; + /* point-weighted average coordinate */ + + if (Npts < 5) continue; + if (!gaussj (a, 3, b, 2)) continue; + + /* we only care about the slopes, not the offsets */ + grad[0].dPdL[Nval] = b[1][0]; + grad[0].dPdM[Nval] = b[2][0]; + grad[0].dQdL[Nval] = b[1][1]; + grad[0].dQdM[Nval] = b[2][1]; + + Nval ++; + + } + } + } + grad[0].Npts = Nval; + if ((DUMP != NULL) && !strcmp (DUMP, "grads")) dump_grads (grad, "grads.dat"); + return (grad); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetScatter.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetScatter.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/GetScatter.c (revision 21807) @@ -0,0 +1,37 @@ +# include "mosastro.h" + +/* measure scatter on the focal plane */ +double GetScatter (int *Nscatter, double *DL, double *DM, int bright) { + + int i, j, Ntotal; + double dL, dM, dL2, dM2, dl, dm, dR; + StarData *raw, *ref; + + Ntotal = 0.0; + dL = dL2 = dM = dM2 = 0; + for (i = 0; i < Nchip; i++) { + raw = chip[i].raw; + ref = chip[i].ref; + for (j = 0; j < chip[i].Nmatch; j++) { + if (raw[j].mask) continue; + if (bright && (raw[j].Mag > INST_BRIGHT)) continue; + dl = raw[j].L - ref[j].L; + dm = raw[j].M - ref[j].M; + dL += dl; + dL2 += SQ(dl); + dM += dm; + dM2 += SQ(dm); + Ntotal ++; + } + } + dL = sqrt (fabs(dL2 / Ntotal - (dL*dL) / (Ntotal*Ntotal))); + dM = sqrt (fabs(dM2 / Ntotal - (dM*dM) / (Ntotal*Ntotal))); + dR = hypot(dL, dM) * 3600.0 * field.project.cdelt1; + *Nscatter = Ntotal; + *DL = dL; + *DM = dM; + return (dR); +} + +/* sigma is returned in arcsec, dL, dM are in pixels */ +/* return dL and dM independently as well */ Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/LoadStars.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/LoadStars.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/LoadStars.c (revision 21807) @@ -0,0 +1,77 @@ +# include "mosastro.h" + +int LoadStars (int Nfile, char **file) { + + int i, itmp, status, extend; + time_t tsval; + struct tm *tmval; + + ALLOCATE (chip, Chip, Nfile); + + /* Nchip is number of valid chips */ + + Nchip = 0; + for (i = 0; i < Nfile; i++) { + chip[Nchip].file = strcreate (file[i]); + + /* load header */ + if (!fits_read_header (chip[Nchip].file, &chip[Nchip].header)) { + fprintf (stderr, "ERROR: can't read header for %s\n", file[i]); + free (chip[Nchip].file); + continue; + } + + /* get image dimensions */ + fits_scan (&chip[Nchip].header, "NAXIS1", "%d", 1, &chip[Nchip].NX); + fits_scan (&chip[Nchip].header, "NAXIS2", "%d", 1, &chip[Nchip].NY); + + /* get astrometry information */ + if (!GetCoords (&chip[Nchip].coords, &chip[Nchip].header)) { + fprintf (stderr, "ERROR: no astrometric solution in header\n"); + free (chip[Nchip].file); + fits_free_header (&chip[Nchip].header); + continue; + } + while (chip[Nchip].coords.crval1 < 0) chip[Nchip].coords.crval1 += 360.0; + while (chip[Nchip].coords.crval1 > 360.0) chip[Nchip].coords.crval1 -= 360.0; + + itmp = 0; + fits_scan (&chip[Nchip].header, "NASTRO", "%d", 1, &itmp); + if (itmp == 0) { + fprintf (stderr, "ERROR: bad astrometric solution in header %s\n", file[i]); + free (chip[Nchip].file); + fits_free_header (&chip[Nchip].header); + continue; + } + + /* get time info */ + tsval = parse_time (&chip[Nchip].header); + tmval = gmtime (&tsval); + Year = tmval[0].tm_year; + + extend = FALSE; + fits_scan (&chip[Nchip].header, "EXTEND", "%t", 1, &extend); + if (extend) { + status = rfits (&chip[Nchip]); + } else { + status = rtext (&chip[Nchip]); + } + if (!status) { + /* skip on failure */ + free (chip[Nchip].file); + fits_free_header (&chip[Nchip].header); + continue; + } + + if (VERBOSE) fprintf (stderr, "loaded %d stars from %s\n", chip[Nchip].Nstars, chip[Nchip].file); + + Nchip ++; + } + + if (Nchip == 0) { + fprintf (stderr, "no valid data in chip files, exiting\n"); + exit (1); + } + + return (TRUE); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/SaveResiduals.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/SaveResiduals.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/SaveResiduals.c (revision 21807) @@ -0,0 +1,97 @@ +# include "mosastro.h" + +void SaveResiduals (FILE *f, Header *header) { + + int i, j, N, Nmatch; + MatchData *match; + Matrix matrix; + Header theader; + FTable table; + + header[0].extend = TRUE; + header[0].Naxes = 0; + fits_modify (header, "NAXIS", "%d", 1, 0); + fits_modify (header, "EXTEND", "%t", 1, TRUE); + fits_modify (header, "NEXTEND", "%d", 1, 1); + + /* add in some keywords to specify the datatype & software version? */ + + /* create (empty) data matrix */ + fits_create_matrix (header, &matrix); + + /* create bintable header */ + fits_create_table_header (&theader, "BINTABLE", "MOSASTRO_RESIDUALS"); + + /* define bintable layout */ + fits_define_bintable_column (&theader, "D", "R_RAW", "ra (raw)", "degrees", 1.0, 0.0); + fits_define_bintable_column (&theader, "D", "D_RAW", "dec (raw)", "degrees", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "P_RAW", "P coord (raw)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "Q_RAW", "Q coord (raw)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "L_RAW", "L coord (raw)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "M_RAW", "M coord (raw)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "X_RAW", "X coord (raw)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "Y_RAW", "Y coord (raw)", "pixels", 1.0, 0.0); + + fits_define_bintable_column (&theader, "D", "R_REF", "ra (ref)", "degrees", 1.0, 0.0); + fits_define_bintable_column (&theader, "D", "D_REF", "dec (ref)", "degrees", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "P_REF", "P coord (ref)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "Q_REF", "Q coord (ref)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "L_REF", "L coord (ref)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "M_REF", "M coord (ref)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "X_REF", "X coord (ref)", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "Y_REF", "Y coord (ref)", "pixels", 1.0, 0.0); + + fits_define_bintable_column (&theader, "E", "MAG_REF", "catalog mag", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "DMAG_REF", "catalog mag err", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "MAG_RAW", "instrum mag", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "E", "DMAG_RAW", "instrum mag err", "pixels", 1.0, 0.0); + fits_define_bintable_column (&theader, "B", "MASK", "excluded from fit?", "", 1.0, 0.0); + fits_define_bintable_column (&theader, "7A", "DUMMY", "padding", "", 1.0, 0.0); + + /* create table, add data values */ + fits_create_table (&theader, &table); + + /* create output data block and assign */ + Nmatch = 0; + for (i = 0; i < Nchip; i++) Nmatch += chip[i].Nmatch; + ALLOCATE (match, MatchData, Nmatch); + + for (i = N = 0; i < Nchip; i++) { + for (j = 0; j < chip[i].Nmatch; j++, N++) { + match[N].Rraw = chip[i].raw[j].R; + match[N].Draw = chip[i].raw[j].D; + match[N].Praw = chip[i].raw[j].P; + match[N].Qraw = chip[i].raw[j].Q; + match[N].Lraw = chip[i].raw[j].L; + match[N].Mraw = chip[i].raw[j].M; + match[N].Xraw = chip[i].raw[j].X; + match[N].Yraw = chip[i].raw[j].Y; + + match[N].Rref = chip[i].ref[j].R; + match[N].Dref = chip[i].ref[j].D; + match[N].Pref = chip[i].ref[j].P; + match[N].Qref = chip[i].ref[j].Q; + match[N].Lref = chip[i].ref[j].L; + match[N].Mref = chip[i].ref[j].M; + match[N].Xref = chip[i].ref[j].X; + match[N].Yref = chip[i].ref[j].Y; + + match[N].Mcat = chip[i].ref[j].Mag; + match[N].dMcat = chip[i].ref[j].dMag; + match[N].Minst = chip[i].raw[j].Mag; + match[N].dMinst = chip[i].raw[j].dMag; + match[N].mask = chip[i].raw[j].mask; + } + } + + /* fix byte order issues */ + ConvertMatch (match, sizeof (MatchData), Nmatch); + fits_add_rows (&table, (char *) match, Nmatch, sizeof (MatchData)); + + fits_fwrite_header (f, header); + fits_fwrite_matrix (f, &matrix); + fits_fwrite_Theader (f, &theader); + fits_fwrite_table (f, &table); + fclose (f); + return; +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args.c (revision 21807) @@ -0,0 +1,68 @@ +# include "mosastro.h" + +void print_help () { + + fprintf (stderr, "mosastro -- mosaic astrometry\n"); + fprintf (stderr, "\n"); + exit (0); + +} + +void args (int *argc, char **argv) { + + int N; + + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + print_help (); + } + + 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); + } + + SAVE_RESID = FALSE; + if ((N = get_argument (*argc, argv, "-save-residuals"))) { + remove_argument (N, argc, argv); + SAVE_RESID = TRUE; + } + + CHIPS = (char *) NULL; + if ((N = get_argument (*argc, argv, "-chips"))) { + remove_argument (N, argc, argv); + CHIPS = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + FIELD = (char *) NULL; + if ((N = get_argument (*argc, argv, "-field"))) { + remove_argument (N, argc, argv); + FIELD = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + /** currently unused **/ + field.Norder = 0; + if ((N = get_argument (*argc, argv, "-order"))) { + remove_argument (N, argc, argv); + field.Norder = atoi (argv[N]); + remove_argument (N, argc, argv); + } + + /** currently unused **/ + ChipOrder = 1; + if ((N = get_argument (*argc, argv, "-chiporder"))) { + remove_argument (N, argc, argv); + ChipOrder = atoi (argv[N]); + remove_argument (N, argc, argv); + } +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args_obs.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args_obs.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/args_obs.c (revision 21807) @@ -0,0 +1,91 @@ +# include "mosastro.h" + +void print_help () { + fprintf (stderr, "USAGE: mkobs (RA) (DEC) (output) [-p param value]\n"); + exit (1); +} + +void args (int *argc, char **argv) { + + int N, No, Np, Nx, Ny; + double theta; + char line[500]; + + if (get_argument (*argc, argv, "--help")) print_help (); + if (get_argument (*argc, argv, "-h")) print_help (); + + FOCAL_PLANE = NULL; + if ((N = get_argument (*argc, argv, "-fp"))) { + remove_argument (N, argc, argv); + FOCAL_PLANE = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + + NO_CHIPS = FALSE; + if ((N = get_argument (*argc, argv, "-nochips"))) { + remove_argument (N, argc, argv); + NO_CHIPS = TRUE; + } + + SIGMA = 0; + if ((N = get_argument (*argc, argv, "-sigma"))) { + remove_argument (N, argc, argv); + SIGMA = atof (argv[N]); + remove_argument (N, argc, argv); + } + + fake_field_defaults (); + + while ((N = get_argument (*argc, argv, "-p"))) { + remove_argument (N, argc, argv); + + if (!strcmp (argv[N], "help")) { + fprintf (stderr, "valid parameters:\n"); + fprintf (stderr, "cdelt (arcsec/pix)\n"); + fprintf (stderr, "crpix (Xo) (Yo)\n"); + fprintf (stderr, "theta (angle)\n"); + exit (2); + } + + if (!strcmp (argv[N], "cdelt")) { + remove_argument (N, argc, argv); + field.project.cdelt2 = field.project.cdelt1 = atof(argv[N])/3600.0; + remove_argument (N, argc, argv); + continue; + } + if (!strcmp (argv[N], "crpix")) { + remove_argument (N, argc, argv); + field.project.crpix1 = atof(argv[N]); + remove_argument (N, argc, argv); + field.project.crpix2 = atof(argv[N]); + remove_argument (N, argc, argv); + continue; + } + if (!strcmp (argv[N], "theta")) { + remove_argument (N, argc, argv); + theta = atof(argv[N]); + remove_argument (N, argc, argv); + field.project.pc1_1 = +cos (RAD_DEG*theta); + field.project.pc1_2 = -sin (RAD_DEG*theta); + field.project.pc2_1 = +sin (RAD_DEG*theta); + field.project.pc2_2 = +cos (RAD_DEG*theta); + continue; + } + if (!strncmp (argv[N], "pca", 3)) { + No = argv[N][3] - '0'; + Nx = argv[N][5] - '0'; + Ny = argv[N][7] - '0'; + if ((Nx + Ny > 3) || (No > 1)) { + fprintf (stderr, "PCA out of range\n"); + exit (1); + } + Np = mkpolyterm (Nx, Ny); + remove_argument (N, argc, argv); + field.distort.polyterms[Np][No] = atof(argv[N]); + remove_argument (N, argc, argv); + continue; + } + } + + if (*argc != 4) print_help (); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/chips.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/chips.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/chips.c (revision 21807) @@ -0,0 +1,61 @@ +# include "mosastro.h" + +/* set chip model based on initial header info */ + +/** when the chips and field terms are treated independently, they must be of + ctype PLY. When we reach a solution, and write out a coupled set of chip + and field terms, then we need to transition to WRP (chip) & DIS (field +**/ + +int init_chips () { + + int i, j; + double R, D, P, Q, L, M, Scale, Det; + + if (CHIPS != (char *) NULL) { + load_chips (CHIPS); + return (1); + } + + for (i = 0; i < Nchip; i++) { + + /* bore site center guess */ + strcpy (chip[i].map.ctype, "DEC--PLY"); + chip[i].map.crpix1 = 0.0; + chip[i].map.crpix2 = 0.0; + chip[i].map.cdelt1 = 1.0; + chip[i].map.cdelt2 = 1.0; + + /* find (L,M) coords of reference pixel (0,0) */ + XY_to_RD (&R, &D, 0.0, 0.0, &chip[i].coords); + RD_to_XY (&P, &Q, R, D, &field.project); + RD_to_XY (&L, &M, P, Q, &field.distort); + chip[i].map.crval1 = L; + chip[i].map.crval2 = M; + + /** we preserve the rotation and parity of coords.pc_ij, but renormalize to unity scale **/ + Det = chip[i].coords.pc1_1*chip[i].coords.pc2_2 - chip[i].coords.pc1_2*chip[i].coords.pc2_1; + Scale = 1.0 / sqrt(fabs(chip[i].coords.cdelt1*chip[i].coords.cdelt2*Det)); + + /** test for NaN Scale **/ + + chip[i].map.pc1_1 = Scale * chip[i].coords.pc1_1 * chip[i].coords.cdelt1; + chip[i].map.pc2_2 = Scale * chip[i].coords.pc2_2 * chip[i].coords.cdelt2; + chip[i].map.pc1_2 = Scale * chip[i].coords.pc1_2 * chip[i].coords.cdelt2; + chip[i].map.pc2_1 = Scale * chip[i].coords.pc2_1 * chip[i].coords.cdelt1; + + chip[i].map.Npolyterms = 1; + for (j = 0; j < 7; j++) { + chip[i].map.polyterms[j][0] = 0; + chip[i].map.polyterms[j][1] = 0; + } + } + return (1); +} + +int load_chips (char *filename) { + + fprintf (stderr, "not ready yet\n"); + exit (1); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/clip.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/clip.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/clip.c (revision 21807) @@ -0,0 +1,35 @@ +# include "mosastro.h" + +int ClipOnFP (double Nsigma) { + + int i, j, Nscatter, Nmask, Nkeep; + double DL, DM, dL, dM; + double sigma; + StarData *raw, *ref; + + Nmask = Nkeep = 0; + sigma = GetScatter (&Nscatter, &DL, &DM, FALSE); + + for (i = 0; i < Nchip; i++) { + raw = chip[i].raw; + ref = chip[i].ref; + for (j = 0; j < chip[i].Nmatch; j++) { + dL = raw[j].L - ref[j].L; + dM = raw[j].M - ref[j].M; + if ((fabs(dL) > Nsigma*DL) || (fabs(dM) > Nsigma*DM)) { + raw[j].mask = TRUE; + Nmask ++; + } else { + raw[j].mask = FALSE; + Nkeep ++; + } + } + } + + fprintf (stderr, "Nmask: %d, Nkeep: %d\n", Nmask, Nkeep); + return (TRUE); +} + +/* + sigma = sigma / (3600.0 * field.project.cdelt1); +*/ Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/dump.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/dump.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/dump.c (revision 21807) @@ -0,0 +1,73 @@ +# include "mosastro.h" + +int dump_stars (FILE *f, StarData *stars, int Nstars) { + + int i; + + for (i = 0; i < Nstars; i++) { + fprintf (f, "%4d %10.6f %10.6f %8.2f %8.2f %8.2f %8.2f %7.2f %7.2f %7.2f %7.2f %d\n", + i, + stars[i].R, stars[i].D, + stars[i].P, stars[i].Q, + stars[i].L, stars[i].M, + stars[i].X, stars[i].Y, + stars[i].Mag, stars[i].dMag, stars[i].mask); + } + return (1); +} + +int dump_grads (Gradients *grad, char *filename) { + + int i; + FILE *f; + + fprintf (stderr, "printing to file %s\n", filename); + + f = fopen (filename, "w"); + + for (i = 0; i < grad[0].Npts; i++) { + fprintf (f, "%4d %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n", + i, + grad[0].dPdL[i], grad[0].dPdM[i], + grad[0].dQdL[i], grad[0].dQdM[i], + grad[0].Lo[i], grad[0].Mo[i]); + } + fclose (f); + return (1); +} + +int dump_match () { + + int i; + FILE *f, *g; + f = fopen ("raw.dat", "w"); + g = fopen ("ref.dat", "w"); + for (i = 0; i < Nchip; i++) { + dump_stars (f, chip[i].raw, chip[i].Nmatch); + dump_stars (g, chip[i].ref, chip[i].Nmatch); + } + fclose (f); + fclose (g); + exit (1); +} + +int dump_rawstars () { + + int i; + FILE *f; + f = fopen ("stars.dat", "w"); + for (i = 0; i < Nchip; i++) { + dump_stars (f, chip[i].stars, chip[i].Nstars); + } + fclose (f); + exit (1); +} + +int dump_refcat (StarData *refcat, int Nrefcat) { + + FILE *f; + f = fopen ("refcat.dat", "w"); + dump_stars (f, refcat, Nrefcat); + fclose (f); + exit (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fakefield.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fakefield.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fakefield.c (revision 21807) @@ -0,0 +1,62 @@ +# include "mosastro.h" + +/* determine an initial guess to field parameters from data */ +int fake_field_center (double RA, double DEC) { + + field.Rmin = RA - 0.5; + field.Rmax = RA + 0.5; + field.Dmin = DEC - 0.5; + field.Dmax = DEC + 0.5; + + /* bore site center guess */ + field.project.crval1 = 0.5*(field.Rmin + field.Rmax); + field.project.crval2 = 0.5*(field.Dmin + field.Dmax); + + return (1); +} + +/* set default field parameters (overridden in args_obs) */ +int fake_field_defaults () { + + int i; + + /* bore site center guess */ + strcpy (field.project.ctype, "DEC--TAN"); + field.project.crpix1 = 0; + field.project.crpix2 = 0; + + /* set TP plate scale */ + field.project.cdelt1 = 1.0/3600.0; + field.project.cdelt2 = 1.0/3600.0; + + /** allow guess at field rotation?? **/ + field.project.pc1_1 = 1; + field.project.pc2_2 = 1; + field.project.pc1_2 = 0; + field.project.pc2_1 = 0; + field.project.Npolyterms = 1; + + /** distort only has power in polyterms **/ + strcpy (field.distort.ctype, "DEC--PLY"); + field.distort.crval1 = 0.0; + field.distort.crval2 = 0.0; + field.distort.crpix1 = 0.0; + field.distort.crpix2 = 0.0; + field.distort.cdelt1 = 1.0; + field.distort.cdelt2 = 1.0; + + field.distort.pc1_1 = 1; + field.distort.pc2_2 = 1; + field.distort.pc1_2 = 0; + field.distort.pc2_1 = 0; + + /* how do we handle renormalization? (fixed at 1000 pixels??) */ + field.distort.Npolyterms = 3; + for (i = 0; i < 7; i++) { + field.distort.polyterms[i][0] = 0; + field.distort.polyterms[i][1] = 0; + } + + return (1); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/field.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/field.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/field.c (revision 21807) @@ -0,0 +1,181 @@ +# include "mosastro.h" + +/* determine an initial guess to field parameters from data */ +int init_field () { + + int i; + + if (FIELD != (char *) NULL) { + load_field (FIELD); + return (TRUE); + } + + /* bore site center guess */ + strcpy (field.project.ctype, "DEC--TAN"); + field.project.crval1 = 0.5*(field.Rmin + field.Rmax); + field.project.crval2 = 0.5*(field.Dmin + field.Dmax); + field.project.crpix1 = 0; + field.project.crpix2 = 0; + + /* measure average plate scale - would be better using parabolic min... */ + field.project.cdelt1 = 0; + field.project.cdelt2 = 0; + for (i = 0; i < Nchip; i++) { + field.project.cdelt1 += chip[i].coords.cdelt1; + field.project.cdelt2 += chip[i].coords.cdelt2; + } + field.project.cdelt1 /= Nchip; + field.project.cdelt2 /= Nchip; + /* force starting guess to have equal x & y plate scales? */ + + /** allow guess at field rotation?? **/ + field.project.pc1_1 = 1; + field.project.pc2_2 = 1; + field.project.pc1_2 = 0; + field.project.pc2_1 = 0; + field.project.Npolyterms = 1; + + /* bore site center guess */ + strcpy (field.distort.ctype, "DEC--PLY"); + field.distort.crval1 = 0.0; + field.distort.crval2 = 0.0; + field.distort.crpix1 = 0.0; + field.distort.crpix2 = 0.0; + field.distort.cdelt1 = 1.0; + field.distort.cdelt2 = 1.0; + + /** allow guess at field rotation?? **/ + field.distort.pc1_1 = 1; + field.distort.pc2_2 = 1; + field.distort.pc1_2 = 0; + field.distort.pc2_1 = 0; + + /* allow 2nd and 3rd order? */ + /* how do we handle renormalization? (fixed at 1000 pixels??) */ + field.distort.Npolyterms = 1; + for (i = 0; i < 7; i++) { + field.distort.polyterms[i][0] = 0; + field.distort.polyterms[i][1] = 0; + } + return (TRUE); +} + +void field_stats () { + + int i, j; + double Rmin, Rmax, Dmin, Dmax; + + /* find range for single image */ + Rmin = Dmin = 360.0; + Rmax = Dmax = -90.0; + + for (i = 0; i < Nchip; i++) { + for (j = 0; j < chip[i].Nstars; j++) { + Rmin = MIN (Rmin, chip[i].stars[j].R); + Dmin = MIN (Dmin, chip[i].stars[j].D); + Rmax = MAX (Rmax, chip[i].stars[j].R); + Dmax = MAX (Dmax, chip[i].stars[j].D); + } + } + field.Rmin = Rmin; + field.Rmax = Rmax; + field.Dmin = Dmin; + field.Dmax = Dmax; + +} + +void field_combine () { + + int i; + double cd1, cd2, pc11, pc12, pc21, pc22; + + /* combine boresite & distortion parameters: ctype DIS */ + strcpy (field.project.ctype, "DEC--DIS"); + + cd1 = field.project.cdelt1; + cd2 = field.project.cdelt2; + pc11 = field.project.pc1_1; + pc12 = field.project.pc1_2; + pc21 = field.project.pc2_1; + pc22 = field.project.pc2_2; + + field.project.Npolyterms = field.distort.Npolyterms; + + for (i = 0; i < 7; i++) { + field.project.polyterms[i][0] = (pc11*cd1*field.distort.polyterms[i][0] + pc12*cd2*field.distort.polyterms[i][1]); + field.project.polyterms[i][1] = (pc21*cd1*field.distort.polyterms[i][0] + pc22*cd2*field.distort.polyterms[i][1]); + } + for (i = 0; i < 2; i++) { + field.project.polyterms[0][i] = field.project.polyterms[0][i] / (cd1*cd1); + field.project.polyterms[1][i] = field.project.polyterms[1][i] / (cd1*cd2); + field.project.polyterms[2][i] = field.project.polyterms[2][i] / (cd2*cd2); + + field.project.polyterms[3][i] = field.project.polyterms[3][i] / (cd1*cd1*cd1); + field.project.polyterms[4][i] = field.project.polyterms[4][i] / (cd1*cd1*cd2); + field.project.polyterms[5][i] = field.project.polyterms[5][i] / (cd1*cd2*cd2); + field.project.polyterms[6][i] = field.project.polyterms[6][i] / (cd2*cd2*cd2); + } +} + +int load_field (char *filename) { + + int i; + Coords coords; + Header header; + + /* load header */ + if (!fits_read_header (filename, &header)) { + fprintf (stderr, "ERROR: can't read header for %s\n", filename); + exit (1); + } + /* get astrometry information */ + if (!GetCoords (&coords, &header)) { + fprintf (stderr, "ERROR: no astrometric solution in field %s\n", filename); + exit (1); + } + + /* separate field into boresite + distortion terms */ + + /* bore site center guess */ + strcpy (field.project.ctype, "DEC--TAN"); + field.project.crval1 = 0.5*(field.Rmin + field.Rmax); + field.project.crval2 = 0.5*(field.Dmin + field.Dmax); + field.project.crpix1 = 0; + field.project.crpix2 = 0; + + /* measure average plate scale - would be better using parabolic min... */ + field.project.cdelt1 = coords.cdelt1; + field.project.cdelt2 = coords.cdelt2; + + /** allow guess at field rotation?? **/ + field.project.pc1_1 = coords.pc1_1; + field.project.pc2_2 = coords.pc2_2; + field.project.pc1_2 = coords.pc1_2; + field.project.pc2_1 = coords.pc2_1; + field.project.Npolyterms = 1; + + /* bore site center guess */ + strcpy (field.distort.ctype, "DEC--PLY"); + field.distort.crval1 = 0.0; + field.distort.crval2 = 0.0; + field.distort.crpix1 = 0.0; + field.distort.crpix2 = 0.0; + field.distort.cdelt1 = 1.0; + field.distort.cdelt2 = 1.0; + + /** allow guess at field rotation?? **/ + field.distort.pc1_1 = 1; + field.distort.pc2_2 = 1; + field.distort.pc1_2 = 0; + field.distort.pc2_1 = 0; + + /* allow 2nd and 3rd order? */ + /* how do we handle renormalization? (fixed at 1000 pixels??) */ + field.distort.Npolyterms = coords.Npolyterms; + for (i = 0; i < 7; i++) { + field.distort.polyterms[i][0] = coords.polyterms[i][0]; + field.distort.polyterms[i][1] = coords.polyterms[i][1]; + } + return (TRUE); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fitpoly.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fitpoly.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/fitpoly.c (revision 21807) @@ -0,0 +1,284 @@ +# include "mosastro.h" + +static int NTERM, NPOWER, NPARS, NORDER, Npts; +static double **sum, **xsum, **ysum; +static double **matrix, **vector; + +void fit_init (int order) { + + int i; + + Npts = 0; + NORDER = order; + NPOWER = NORDER + 1; + NTERM = 2*NORDER + 1; + NPARS = (NORDER + 1)*(NORDER + 2) / 2; + + /* allocate arrays for fit solution */ + ALLOCATE (sum, double *, NTERM); + ALLOCATE (xsum, double *, NTERM); + ALLOCATE (ysum, double *, NTERM); + for (i = 0; i < NTERM; i++) { + ALLOCATE (sum[i], double, NTERM); + bzero (sum[i], NTERM*sizeof(double)); + ALLOCATE (xsum[i], double, NTERM); + bzero (xsum[i], NTERM*sizeof(double)); + ALLOCATE (ysum[i], double, NTERM); + bzero (ysum[i], NTERM*sizeof(double)); + } + ALLOCATE (matrix, double *, NPARS); + ALLOCATE (vector, double *, NPARS); + for (i = 0; i < NPARS; i++) { + ALLOCATE (matrix[i], double, NPARS); + ALLOCATE (vector[i], double, 2); + bzero (vector[i], 2*sizeof(double)); + bzero (matrix[i], NPARS*sizeof(double)); + } + +} + +# define SCALE 1.0 +void fit_add (double x1, double y1, double x2, double y2) { + + int n, m; + double xterm, yterm, term; + + xterm = 1; + for (n = 0; n < NTERM; n++) { + yterm = 1; + for (m = 0; m < NTERM; m++) { + term = xterm*yterm; + if (n+m < NTERM) { + sum[n][m] += term; + } + if (n+m < NPOWER) { + xsum[n][m] += x2*term; + ysum[n][m] += y2*term; + } + yterm *= y1/SCALE; + } + xterm *= x1/SCALE; + } + Npts ++; +} + +/** I am renormalizing here by the max pivots to keep gaussj sane ** + ** would not be needed if the fit used scaled ind. variables **/ +void fit_eval () { + + int i, j, n, m, M, N; + + if (Npts == 0) { + fprintf (stderr, "warning: no valid pts\n"); + } + + i = 0; + for (m = 0; m < NPOWER; m++) { + for (n = 0; n < NPOWER - m; n++, i++) { + vector[i][0] = xsum[n][m]; + vector[i][1] = ysum[n][m]; + } + } + j = 0; + for (M = 0; M < NPOWER; M++) { + for (N = 0; N < NPOWER - M; N++, j++) { + i = 0; + for (m = 0; m < NPOWER; m++) { + for (n = 0; n < NPOWER - m; n++, i++) { + matrix[i][j] = sum[n+N][m+M]; + } + } + } + } +# if (0) + max = 0.0; + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + max = MAX (max, fabs(matrix[i][j])); + } + max = MAX (max, fabs(vector[i][0])); + max = MAX (max, fabs(vector[i][1])); + } + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + matrix[i][j] /= max; + } + vector[i][0] /= max; + vector[i][1] /= max; + } +# endif + + gaussj (matrix, NPARS, vector, 2); + +# if (0) + i = 0; + for (m = 0; m < NPOWER; m++) { + for (n = 0; n < NPOWER - m; n++, i++) { + fprintf (stderr, "RA x^%dy^%d: %10.4g DEC x^%dy^%d: %10.4g \n", + n, m, vector[i][0], n, m, vector[i][1]); + } + } +# endif +} + +/* linear portion of fit : NORDER is 1 */ +void fit_apply_coords (Coords *coords) { + + int i, j, Np, Nv, N; + double c11, c12; + double c21, c22; + double R; + + /* update the higher order terms */ + if (NORDER > 1) { + for (i = 0; i < NPOWER; i++) { + for (j = 0; j < (NPOWER - i); j++) { + if (i + j < 2) continue; + Np = mkpolyterm (i, j); + Nv = mkvector (i, j, NORDER); + coords[0].polyterms[Np][0] = vector[Nv][0]; + coords[0].polyterms[Np][1] = vector[Nv][1]; + } + } + } + + /* get the correct vector entries for the linear terms */ + N = mkvector (0, 0, NORDER); + coords[0].crval1 = vector[N][0]; + coords[0].crval2 = vector[N][1]; + + N = mkvector (1, 0, NORDER); + c11 = vector[N][0]; + c21 = vector[N][1]; + N = mkvector (0, 1, NORDER); + c12 = vector[N][0]; + c22 = vector[N][1]; + coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21)); + R = 1 / coords[0].cdelt1; + + coords[0].pc1_1 = c11*R; + coords[0].pc2_1 = c21*R; + coords[0].pc1_2 = c12*R; + coords[0].pc2_2 = c22*R; + + coords[0].crpix1 = 0; + coords[0].crpix2 = 0; + + coords[0].Npolyterms = NORDER; + strcpy (coords[0].ctype, "DEC--PLY"); +} + +/* + if we have just linear terms, the following holds for crpix1,2: + D = R / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + coords[0].crpix1 = D * (coords[0].pc1_2*c20 - coords[0].pc2_2*c10); + coords[0].crpix2 = D * (coords[0].pc2_1*c10 - coords[0].pc1_1*c20); +*/ + +/* NORDER is order of gradient fit : Npolyterms is Norder + 1 */ +void fit_apply_grads (Coords *distort, Coords *project, int term) { + + int i, j, Np, Nv1, Nv2; + + distort[0].Npolyterms = NORDER + 1; + if (distort[0].Npolyterms > 1) strcpy (distort[0].ctype, "DEC--PLY"); + + for (i = 0; i < NPOWER + 1; i++) { + for (j = 0; j < (NPOWER + 1 - i); j++) { + + if (i + j < 2) continue; + Np = mkpolyterm (i, j); + Nv1 = mkvector (i-1, j, NORDER); + Nv2 = mkvector (i, j-1, NORDER); + + /** why do we have the negative sign? **/ + if (j == 0) { + distort[0].polyterms[Np][term] = vector[Nv1][0] / i; + } + if (i == 0) { + distort[0].polyterms[Np][term] = vector[Nv2][1] / j; + } + if ((i > 0) && (j > 0)) { + distort[0].polyterms[Np][term] = 0.5*(vector[Nv1][0] / i + vector[Nv2][1] / j); + } + } + } + + Nv1 = mkvector (0, 0, NORDER); + if (term == 0) { + project[0].pc1_1 = project[0].pc1_1 * (1 + vector[Nv1][0]); + project[0].pc1_2 = project[0].pc1_2 * (1 + vector[Nv1][1]); + } else { + project[0].pc2_1 = project[0].pc2_1 * (1 + vector[Nv1][0]); + project[0].pc2_2 = project[0].pc2_2 * (1 + vector[Nv1][1]); + } +} + +void fit_correct_grads (Gradients *in, Gradients *out, int term) { + + int i, k, m, n; + double x, y, dx, dy, dz1, dz2; + + for (i = 0; i < in[0].Npts; i++) { + + dx = in[0].Lo[i]; + dy = in[0].Mo[i]; + dz1 = dz2 = 0.0; + + k = 0; + x = y = 1; + for (m = 0; m < NPOWER; m++) { + x = y; + for (n = 0; n < NPOWER - m; n++, k++) { + dz1 += vector[k][0]*x; + dz2 += vector[k][1]*x; + x = x * dx / SCALE; + } + y = y * dy / SCALE; + } + + out[0].Lo[i] = dx; + out[0].Mo[i] = dy; + if (term == 0) { + out[0].dPdL[i] = in[0].dPdL[i] - dz1; + out[0].dPdM[i] = in[0].dPdM[i] - dz2; + } else { + out[0].dQdL[i] = in[0].dQdL[i] - dz1; + out[0].dQdM[i] = in[0].dQdM[i] - dz2; + } + } +} + +int mkvector (int n, int m, int norder) { + + int i, N; + + N = 0; + for (i = 0; i < m; i++) { + N += (norder - i + 1); + } + N += n; + return (N); +} + +void fit_free () { + + int i; + + for (i = 0; i < NTERM; i++) { + free (sum[i]); + free (xsum[i]); + free (ysum[i]); + } + free (sum); + free (xsum); + free (ysum); + + for (i = 0; i < NPARS; i++) { + free (matrix[i]); + free (vector[i]); + } + free (matrix); + free (vector); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gaussj.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gaussj.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gaussj.c (revision 21807) @@ -0,0 +1,76 @@ +# include "mosastro.h" + +int gaussj (double **a, int n, double **b, int m) { + + int *indxc,*indxr,*ipiv; + int i, icol, irow, j, k, l, ll; + double big,dum,pivinv; + + ALLOCATE (indxc, int, n); + ALLOCATE (indxr, int, n); + ALLOCATE (ipiv, int, n); + for (j = 0; j < n; j++) + ipiv[j] = 0; + + icol = irow = 0; + for (i = 0; i < n; i++) { + big = 0.0; + for (j = 0; j < n; j++) { + if (ipiv[j] != 1) { + for (k = 0; k < n; k++) { + if (ipiv[k] == 0) { + if (fabs (a[j][k]) >= big) { + big = fabs (a[j][k]); + irow = j; + icol = k; + } + } + else + if (ipiv[k] > 1) { + /* fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); */ + return (0); + } + } + } + } + ipiv[icol]++; + if (irow != icol) { + for (l = 0; l < n; l++) + SWAP (a[irow][l], a[icol][l]); + for (l = 0; l < m; l++) + SWAP (b[irow][l], b[icol][l]); + } + indxr[i] = irow; + indxc[i] = icol; + if (a[icol][icol] == 0.0) { + /* fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); */ + return (0); + } + pivinv = 1.0 / a[icol][icol]; + a[icol][icol] = 1.0; + for (l = 0; l < n; l++) + a[icol][l] *= pivinv; + for (l = 0; l < m; l++) + b[icol][l] *= pivinv; + for (ll = 0; ll < n; ll++) { + if (ll != icol) { + dum = a[ll][icol]; + a[ll][icol] = 0.0; + for (l = 0; l < n; l++) + a[ll][l] -= a[icol][l]*dum; + for (l = 0; l < m; l++) + b[ll][l] -= b[icol][l]*dum; + } + } + } + + for (l = n - 1; l >= 0; l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < n; k++) + SWAP (a[k][indxr[l]], a[k][indxc[l]]); + } + free (ipiv); + free (indxr); + free (indxc); + return (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gcatalog.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gcatalog.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gcatalog.c (revision 21807) @@ -0,0 +1,48 @@ +# include "mosastro.h" + +# define BYTES_STAR 23 +# define BLOCK 1000 + +StarData *gcatalog (char *filename, int *Nstars) { + + StarData *stars; + int i, NSTAR, nstar, Nbytes, nbytes; + char *buffer; + FILE *f; + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find catalog file %s\n", filename); + return (NULL); + } + + nstar = 0; + NSTAR = 1000; + ALLOCATE (stars, StarData, NSTAR); + + Nbytes = BLOCK*BYTES_STAR; + ALLOCATE (buffer, char, Nbytes); + while ((nbytes = fread (buffer, 1, Nbytes, f)) > 0) { + for (i = 0; i < nbytes / BYTES_STAR; i++) { + bzero (&stars[nstar], sizeof(StarData)); + dparse (&stars[nstar].R, 1, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].D, 2, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].Mag, 3, &buffer[i*BYTES_STAR]); + nstar++; + if (nstar == NSTAR) { + NSTAR += 1000; + REALLOCATE (stars, StarData, NSTAR); + } + } + } + free (buffer); + + *Nstars = nstar; + return (stars); +} + +/** this function reads a file in the format of my extracted GSC files + this is a fairly weak way of storing photometry / astrometry data + eventually, this should be moved to ptolemy / DVO or something else + or else, the existing 2mass and HST-GSC files should be parsed +**/ Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/get2mass.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/get2mass.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/get2mass.c (revision 21807) @@ -0,0 +1,35 @@ +# include "mosastro.h" + +StarData *get2mass (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Ngsc, Nregions, Nstars; + StarData *gsc; + StarData *stars; + CatStats *regions; + + Nstars = 0; + ALLOCATE (stars, StarData, 1); + + set_catalog (TWO_MASS_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + if (gsc == NULL) continue; + + Ns = Nstars; + Nstars += Ngsc; + + REALLOCATE (stars, StarData, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + stars[k].R = gsc[j].R; + stars[k].D = gsc[j].D; + stars[k].M = gsc[j].M; + } + free (gsc); + } + + if (VERBOSE) fprintf (stderr, "%d stars from 2MASS\n", Nstars); + *NSTARS = Nstars; + return (stars); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getgsc.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getgsc.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getgsc.c (revision 21807) @@ -0,0 +1,36 @@ +# include "mosastro.h" + +StarData *getgsc (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Ngsc, Nregions, Nstars; + StarData *gsc; + StarData *stars; + CatStats *regions; + + Nstars = 0; + ALLOCATE (stars, StarData, 1); + + set_catalog (GSC_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + if (gsc == NULL) continue; + + Ns = Nstars; + Nstars += Ngsc; + + REALLOCATE (stars, StarData, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + stars[k].R = gsc[j].R; + stars[k].D = gsc[j].D; + stars[k].Mag = gsc[j].Mag; + } + free (gsc); + } + + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Nstars); + *NSTARS = Nstars; + return (stars); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getptolemy.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getptolemy.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getptolemy.c (revision 21807) @@ -0,0 +1,32 @@ +# include "mosastro.h" + +StarData *getptolemy (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Nptolemy, Nregions, Nstars; + CatStats *regions; + StarData *stars, *ptolemy; + + Nstars = 0; + ALLOCATE (stars, StarData, 1); + + set_catalog (CATDIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + ptolemy = gptolemy (regions[i].name, &Nptolemy); + if (ptolemy == NULL) continue; + + Ns = Nstars; + Nstars += Nptolemy; + + REALLOCATE (stars, StarData, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Nptolemy; k++, j++) { + stars[k] = ptolemy[j]; + } + free (ptolemy); + } + + if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", Nstars); + *NSTARS = Nstars; + return (stars); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getstone.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getstone.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getstone.c (revision 21807) @@ -0,0 +1,89 @@ +# include "mosastro.h" + +StarData *getstone (CatStats *input, int *nstars) { + + FILE *f; + int i, Nregion, NREGION, Nstars, NSTARS; + double r, d, R, D, M, Tr, Td, dRdT, dDdT; + char **regname, filename[1024], line[1024]; + StarData *stars; + CatStats *regions; + + sprintf (filename, "%s/Regions.dat", StoneRegions); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open stone %s\n", filename); + exit (1); + } + + Nregion = 0; + NREGION = 20; + ALLOCATE (regions, CatStats, NREGION); + ALLOCATE (regname, char *, NREGION); + + /* strip off first commented line */ + scan_line (f, line); + while (scan_line (f, line) != EOF) { + sscanf (line, "%s %lf %lf %lf %lf", filename, + ®ions[Nregion].RA[0], ®ions[Nregion].RA[1], + ®ions[Nregion].DEC[0], ®ions[Nregion].DEC[1]); + sprintf (line, "%s/%s", StoneRegions, filename); + regname[Nregion] = strcreate (line); + Nregion ++; + if (Nregion == NREGION) { + NREGION += 50; + REALLOCATE (regions, CatStats, NREGION); + REALLOCATE (regname, char *, NREGION); + } + } + fclose (f); + + Nstars = 0; + NSTARS = 1000; + ALLOCATE (stars, StarData, NSTARS); + + /* find the region(s) that overlap the input region */ + for (i = 0; i < Nregion; i++) { + if (input[0].RA[0] > regions[i].RA[1]) continue; + if (input[0].RA[1] < regions[i].RA[0]) continue; + if (input[0].DEC[0] > regions[i].DEC[1]) continue; + if (input[0].DEC[1] < regions[i].DEC[0]) continue; + /* this will fail on 0,360 boundary */ + fprintf (stderr, "loading data from %s\n", regname[i]); + + f = fopen (regname[i], "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open stone %s\n", regname[i]); + exit (1); + } + + while (scan_line (f, line) != EOF) { + + dparse (&R, 3, line); + dparse (&D, 4, line); + dparse (&M, 5, line); + + dparse (&Tr, 12, line); + dparse (&Td, 13, line); + dparse (&dRdT, 14, line); + dparse (&dDdT, 15, line); + + r = R + dRdT*(Year - Tr)/(100.0*3600.0); + d = D + dDdT*(Year - Td)/(100.0*3600.0); + + stars[Nstars].M = M; + stars[Nstars].R = r*15.0; + stars[Nstars].D = d; + Nstars ++; + if (Nstars == NSTARS) { + NSTARS += 5000; + REALLOCATE (stars, StarData, NSTARS); + } + } + fclose (f); + } + + *nstars = Nstars; + return (stars); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getusno.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getusno.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/getusno.c (revision 21807) @@ -0,0 +1,137 @@ +# include "mosastro.h" +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +StarData *getusno (CatStats *catstats, int *Nstars) { + + long int offset; + int i, bin, first, last, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128], c; + FILE *f; + int iRA0, iRA1, iDEC0, iDEC1; + double RA0, RA1, DEC0, DEC1, dec; + int spd, spd_start, spd_end, disk; + int NUSNO, Nusno; + StarData *stars; + + RA0 = catstats[0].RA[0]; + RA1 = catstats[0].RA[1]; + DEC0 = catstats[0].DEC[0]; + DEC1 = catstats[0].DEC[1]; + + /* identify ra & dec range of interest */ + iRA0 = RA0 * 360000.0; + iRA1 = RA1 * 360000.0; + iDEC0 = (DEC0 + 90.0) * 360000.0; + iDEC1 = (DEC1 + 90.0) * 360000.0; + + /* data is organized in south-pole distance zones */ + spd_start = (int)((DEC0 + 90) / 7.5) * 75.0; + dec = (DEC1 + 90) / 7.5; + if (dec > (int)(dec)) { + spd_end = (int)(1 + (DEC1 + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (DEC1 + 90) / 7.5) * 75.0; + } + + Nusno = 0; + NUSNO = 5000; + ALLOCATE (stars, StarData, NUSNO); + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (1); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = RA0 / 3.75; + if ((RA1 / 3.75) == (int) (RA1 / 3.75)) + last = RA1 / 3.75; + else + last = 1 + RA1 / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", CDROM, spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s\n", filename); + exit (1); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failure reading data from file %s\n", filename); + exit (1); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if ((buf[0] > iRA0) && (buf[0] < iRA1) && + (buf[1] > iDEC0) && (buf[1] < iDEC1)) { + bzero (&stars[Nusno], sizeof(StarData)); + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + stars[Nusno].Mag = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); + /* b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); */ + Nusno ++; + if (Nusno == NUSNO) { + NUSNO += 5000; + REALLOCATE (stars, StarData, NUSNO); + } + } + } + free (buffer); + } + fclose (f); + } + + *Nstars = Nusno; + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno); + return (stars); +} + + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gptolemy.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gptolemy.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gptolemy.c (revision 21807) @@ -0,0 +1,43 @@ +# include "mosastro.h" + +StarData *gptolemy (char *filename, int *NSTARS) { + + int i, Nstars; + Catalog catalog; + StarData *stars; + + Nstars = 0; + catalog.filename = filename; + switch (lock_catalog (&catalog, LCK_SOFT)) { + case 0: + case 2: + fprintf (stderr, "can't lock catalog data %s\n", filename); + *NSTARS = Nstars; + return (NULL); + case 1: + break; + } + + /* CATMODE and CATFORMAT determined from catalog */ + catalog.catflags = LOAD_AVES | LOAD_MEAS; + if (!load_catalog (&catalog, FALSE)) { + fprintf (stderr, "can't load catalog data %s\n", filename); + *NSTARS = Nstars; + return (NULL); + } + unlock_catalog (&catalog); + + Nstars = catalog.Naverage; + ALLOCATE (stars, StarData, MAX (1, Nstars)); + for (i = 0; i < catalog.Naverage; i++) { + bzero (&stars[i], sizeof(StarData)); + stars[i].R = catalog.average[i].R_PS; + stars[i].D = catalog.average[i].D_PS; + stars[i].Mag = catalog.measure[catalog.average[i].offset].M_PS; + } + free (catalog.average); + free (catalog.measure); + + *NSTARS = Nstars; + return (stars); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/greference.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/greference.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/greference.c (revision 21807) @@ -0,0 +1,54 @@ +# include "mosastro.h" + +StarData *greference (int *Nrefcat) { + + int Nstars; + StarData *stars; + CatStats catstats; + + if (VERBOSE) fprintf (stderr, "loading astrometric reference data from %s\n", REFCAT); + + stars = NULL; + Nstars = 0; + catstats.RA[0] = field.Rmin; + catstats.RA[1] = field.Rmax; + catstats.DEC[0] = field.Dmin; + catstats.DEC[1] = field.Dmax; + + if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", catstats.RA[0], catstats.RA[1], catstats.DEC[0], catstats.DEC[1]); + + /* get stars from the Stone et al catalog for the given region */ + if (!strcmp (REFCAT, "STONE")) { + stars = getstone (&catstats, &Nstars); + } + + /* get stars from the USNO A catalog for the given region */ + if (!strcmp (REFCAT, "USNO")) { + stars = getusno (&catstats, &Nstars); + } + + /* get stars from the HST GSC catalog for the given region */ + if (!strcmp (REFCAT, "GSC")) { + stars = getgsc (&catstats, &Nstars); + } + + /* get stars from 2MASS for the given region -- add PHOTCODE check? */ + if (!strcmp (REFCAT, "2MASS")) { + strcpy (CATDIR, TWO_MASS_DIR); + stars = getptolemy (&catstats, &Nstars); + /* stars = get2mass (&catstats, &Nstars); */ + } + + /* get stars from the DVO CATDIR for the given region */ + if (!strcmp (REFCAT, "PTOLEMY")) { + stars = getptolemy (&catstats, &Nstars); + } + + if (Nstars == 0) { + fprintf (stderr, "no ref objs: %s\n", REFCAT); + exit (1); + } + *Nrefcat = Nstars; + return (stars); + +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gregions.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gregions.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/gregions.c (revision 21807) @@ -0,0 +1,255 @@ +# include "mosastro.h" +# define MAX_NAME 256 + +static int Nregion, NREGION; +static CatStats *region; +static char *CatDir; + +/* data which defines the layout of the regions */ +# define NDECBANDS 12 +double DecBand[] = {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}; + +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, 0}; + +int NRaBlocks [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0, + 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0}; + +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"}; + +/* given patch on sky with RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */ +CatStats *gregions (CatStats *patch, int *nregion) { + + CatStats area; + double dtmp; + char file[MAX_NAME]; + + init_regions (); + + /* force RA[i] to be in range 0 - 360 , DEC[0] < DEC[1] */ + while (patch[0].RA[0] > 360.0) { patch[0].RA[0] -= 360.0; } + while (patch[0].RA[0] < 0.0) { patch[0].RA[0] += 360.0; } + while (patch[0].RA[1] > 360.0) { patch[0].RA[1] -= 360.0; } + while (patch[0].RA[1] < 0.0) { patch[0].RA[1] += 360.0; } + if (patch[0].DEC[0] > patch[0].DEC[1]) { + dtmp = patch[0].DEC[0]; + patch[0].DEC[0] = patch[0].DEC[1]; + patch[0].DEC[1] = dtmp; + } + + /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!) */ + if (patch[0].DEC[1] > 86.25) { + area.RA[0] = 0.00; + area.RA[1] = 360.00; + area.DEC[0] = 86.25; + area.DEC[1] = 90.00; + sprintf (file, "n8230/pole"); + area.name = strcreate (file); + add_to_regions (&area); + free (area.name); + } + + /* if RA[0] > RA[1], split in two passes */ + if (patch[0].RA[0] > patch[0].RA[1]) { + area = patch[0]; + area.RA[0] = 0.0; + find_dec_bands (&area); + + area = patch[0]; + area.RA[1] = 360.0; + find_dec_bands (&area); + } else { + area = patch[0]; + find_dec_bands (&area); + } + + *nregion = Nregion; + return (region); +} + +/* given an area [RA[0] - RA[1], DEC[0] - DEC[1], + where 0.0 < RA[0] < RA[1] < 360.0 + and -90.0 < (DEC[0], DEC[1]) < 90.0 + and DEC[0] < DEC[1] + search through the dec bands to identify ones which overlap the area */ + +int find_dec_bands (CatStats *area) { + + int i; + + /* for NORTH region (DecBand[i] < DecBand[i+1]) */ + for (i = 0; i < NDECBANDS; i++) { + if ((area[0].DEC[1] > DecBand[i]) || (area[0].DEC[0] < DecBand[i+1])) { + load_ra_blocks (i, area); + } + } + + /* for SOUTH region (DecBand[i] > DecBand[i+1]) */ + for (i = NDECBANDS+1; i < 2*NDECBANDS; i++) { + if ((area[0].DEC[0] < DecBand[i]) || (area[0].DEC[1] > DecBand[i+1])) { + load_ra_blocks (i, area); + } + } + return (TRUE); +} + +/* given the number of a DecBlock and a defined area, return all area names in that block */ +/* area is defined so that RA[0] < RA[1], DEC[0] < DEC[1] */ +int load_ra_blocks (int Ndec, CatStats *area) { + + CatStats tregion; + int i, Nlines, Nskip; + char *buffer, file[MAX_NAME]; + FILE *f; + + /* count lines before Dec Band, and in Dec Band */ + Nskip = 0; + for (i = 0; i < Ndec; i++) Nskip += NDecLines[i]; + Nlines = NDecLines[Ndec]; + + ALLOCATE (buffer, char, 48*Nlines); + ALLOCATE (tregion.name, char, MAX_NAME); + + /* get data from GSC file */ + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE); + exit (1); + } + fseek (f, 5*2880 + 48*Nskip, SEEK_SET); + fread (buffer, 48*Nlines, 1, f); + fclose (f); + + /* find overlapping regions */ + for (i = 0; i < Nlines; i++) { + parse_GSC_line (&tregion, &buffer[i*48]); + if (area[0].RA[0] > tregion.RA[1]) continue; + if (area[0].RA[1] < tregion.RA[0]) continue; + if (area[0].DEC[0] > tregion.DEC[1]) continue; + if (area[0].DEC[1] < tregion.DEC[0]) continue; + sprintf (file, "%s/%s", Dec2Sections[Ndec], tregion.name); + strcpy (tregion.name, file); + add_to_regions (&tregion); + } + + free (buffer); + free (tregion.name); + return (TRUE); +} + +/* parse the GSC line, placing values in 'region'. force RA[0] < RA[1], DEC[0] < DEC[1] */ +int parse_GSC_line (CatStats *tregion, char *line) { + + double v0, s, dtmp; + char tmp[10]; + + /* copy name (4 chars) */ + bzero (tregion[0].name, 10); + strncpy (tregion[0].name, &line[1], 4); + + /* get RA[0] */ + bzero (tmp, 10); + strncpy (tmp, &line[7], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[10], 2); + v0 += atof (tmp) / 60.0; + + bzero (tmp, 10); + strncpy (tmp, &line[13], 5); + v0 += atof (tmp) / 3600.0; + tregion[0].RA[0] = v0 * 15.0; + + /* get RA[1] */ + bzero (tmp, 10); + strncpy (tmp, &line[19], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[22], 2); + v0 += atof (tmp) / 60.0; + + bzero (tmp, 10); + strncpy (tmp, &line[25], 5); + v0 += atof (tmp) / 3600.0; + tregion[0].RA[1] = v0 * 15.0; + + /* this only happens at the 0,360 boundary */ + if (tregion[0].RA[1] < tregion[0].RA[0]) tregion[0].RA[1] += 360.0; + + /* get DEC[0] */ + if (line[31] == '-') s = -1; + else s = +1; + + bzero (tmp, 10); + strncpy (tmp, &line[32], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[35], 2); + v0 += atof (tmp) / 60.0; + tregion[0].DEC[0] = s * v0; + + /* get DEC[1] */ + if (line[40] == '-') s = -1; + else s = +1; + + bzero (tmp, 10); + strncpy (tmp, &line[41], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[44], 2); + v0 += atof (tmp) / 60.0; + tregion[0].DEC[1] = s * v0; + + if (tregion[0].DEC[0] > tregion[0].DEC[1]) { + dtmp = tregion[0].DEC[1]; + tregion[0].DEC[1] = tregion[0].DEC[0]; + tregion[0].DEC[0] = dtmp; + } + return (TRUE); +} + +void init_regions () { + + Nregion = 0; + NREGION = 100; + + ALLOCATE (region, CatStats, NREGION); + +} + +void add_to_regions (CatStats *area) { + + region[Nregion] = area[0]; + ALLOCATE (region[Nregion].name, char, MAX_NAME); + sprintf (region[Nregion].name, "%s/%s.%s", CatDir, area[0].name, "cpt"); + + area_of_region (®ion[Nregion]); + + Nregion ++; + + if (Nregion == NREGION) { + NREGION += 100; + REALLOCATE (region, CatStats, NREGION); + } + +} + +void area_of_region (CatStats *region) { + + double area; + + area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG)); + region[0].Area = area; +} + +void set_catalog (char *catdir) { + CatDir = catdir; +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/match.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/match.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/match.c (revision 21807) @@ -0,0 +1,117 @@ +# include "mosastro.h" +/* match is done on tangent plane : should be focal plane in some cases? */ + +int match (StarData *refcat, int Nrefcat) { + + int i, j, k, K, Ntotal; + int Nmatch, NMATCH; + double dp, dq, radius, Radius; + double *p, *q, *P, *Q; + int *U, *u; + + /* requested radius is in arcsec ; internally radius is in pixels */ + Radius = RADIUS / (3600.0 * field.project.cdelt2); + + /* sort the REFCAT data by P in tangent plane */ + ALLOCATE (P, double, Nrefcat); + ALLOCATE (Q, double, Nrefcat); + ALLOCATE (U, int, Nrefcat); + for (i = 0; i < Nrefcat; i++) { + U[i] = i; + P[i] = refcat[i].P; + Q[i] = refcat[i].Q; + } + sortthree (P, Q, U, Nrefcat); + + + for (i = 0; i < Nchip; i++) { + + Nmatch = 0; + NMATCH = 1000; + ALLOCATE (chip[i].raw, StarData, NMATCH); + ALLOCATE (chip[i].ref, StarData, NMATCH); + + /* sort the star data by P in tangent plane */ + ALLOCATE (p, double, chip[i].Nstars); + ALLOCATE (q, double, chip[i].Nstars); + ALLOCATE (u, int, chip[i].Nstars); + for (j = 0; j < chip[i].Nstars; j++) { + u[j] = j; + p[j] = chip[i].stars[j].P; + q[j] = chip[i].stars[j].Q; + } + sortthree (p, q, u, chip[i].Nstars); + + /* find star matches in the tangent plane coord system */ + for (j = k = 0; (j < chip[i].Nstars) && (k < Nrefcat); ) { + + /** instrumental magnitude limits for raw data **/ + if (IMAG_MIN && (chip[i].stars[u[j]].Mag - ZERO_POINT < IMAG_MIN)) { + j++; + continue; + } + if (IMAG_MAX && (chip[i].stars[u[j]].Mag - ZERO_POINT > IMAG_MAX)) { + j++; + continue; + } + /* skip anything with dMag too large */ + if (SIGMA_LIM > 0.0) { + if (chip[i].stars[u[j]].dMag > SIGMA_LIM) { + j++; + continue; + } + } + + dp = p[j] - P[k]; + + if (dp <= -2*Radius) { + j++; + continue; + } + if (dp >= 2*Radius) { + k++; + continue; + } + + K = k; + for (; (dp > -2*Radius) && (k < Nrefcat); k++) { + dp = p[j] - P[k]; + dq = q[j] - Q[k]; + radius = hypot (dp, dq); + if (radius < Radius) { + chip[i].ref[Nmatch] = refcat[U[k]]; + chip[i].raw[Nmatch] = chip[i].stars[u[j]]; + chip[i].raw[Nmatch].Mag -= ZERO_POINT; /* raw in instrumental mags */ + Nmatch ++; + if (Nmatch == NMATCH) { + NMATCH += 1000; + REALLOCATE (chip[i].raw, StarData, NMATCH); + REALLOCATE (chip[i].ref, StarData, NMATCH); + } + goto done; + } + } + done: + k = K; + j++; + } + REALLOCATE (chip[i].raw, StarData, MAX (1, Nmatch)); + REALLOCATE (chip[i].ref, StarData, MAX (1, Nmatch)); + chip[i].Nmatch = Nmatch; + free (p); + free (q); + free (u); + } + free (P); + free (Q); + free (U); + + Ntotal = 0; + for (i = 0; i < Nchip; i++) { + Ntotal += chip[i].Nmatch; + if (VERBOSE) fprintf (stderr, "chip %d: %d of %d stars\n", + i, chip[i].Nmatch, chip[i].Nstars); + } + fprintf (stderr, "Nchips: %d Nmatch: %d\n", Nchip, Ntotal); + return (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkheader.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkheader.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkheader.c (revision 21807) @@ -0,0 +1,32 @@ +# include "mosastro.h" + +Header *mkheader (int Nx, int Ny, int Nstars, Coords *coords) { + + Header *header; + + ALLOCATE (header, Header, 1); + + header[0].simple = TRUE; + header[0].bitpix = -32; + header[0].unsign = FALSE; + header[0].bscale = 1.0; + header[0].bzero = 0.0; + header[0].Naxes = 2; + header[0].Naxis[0] = Nx; + header[0].Naxis[1] = Ny; + + fits_create_header (header); + + fits_modify (header, "NSTARS", "%d", 1, Nstars); + fits_modify (header, "PHOTCODE", "%s", 1, "STD.R"); + fits_modify (header, "DATE-OBS", "%s", 1, "2004-04-22"); + fits_modify (header, "UTC-OBS", "%s", 1, "14:27:45.30"); + fits_modify (header, "ZERO_PT", "%lf", 1, 25.0); + fits_modify (header, "EXPTIME", "%lf", 1, 2.0); + + PutCoords (coords, header); + + fits_modify (header, "NASTRO", "%d", 1, 1); + + return (header); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkmosaic.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkmosaic.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkmosaic.c (revision 21807) @@ -0,0 +1,90 @@ +# include "mosastro.h" + +Header *mkmosaic (int Nx, int Ny, int Nstars, Coords *coords) { + + int Nastro; + double Xmin, Xmax, Ymin, Ymax, tmp; + double DL, DM, Cerror, Cprecise; + Header *header; + char line[80]; + + /* mosaic header has info to define the outline of the mosaic */ + RD_to_XY (&Xmin, &Ymin, field.Rmin, field.Dmin, coords); + RD_to_XY (&Xmax, &Ymax, field.Rmax, field.Dmax, coords); + + ALLOCATE (header, Header, 1); + + header[0].bitpix = -32; + header[0].unsign = FALSE; + header[0].extend = FALSE; + header[0].bscale = 1.0; + header[0].bzero = 0.0; + header[0].Naxes = 2; + header[0].Naxis[0] = Xmax - Xmin; + header[0].Naxis[1] = Ymax - Ymin; + + fits_create_header (header); + + /* calculate image-wide astrometry scatter */ + Cerror = GetScatter (&Nastro, &DL, &DM, FALSE); + Cprecise = Cerror / sqrt (Nastro); + fits_modify (header, "CERROR", "%lf", 1, Cerror); + fits_modify (header, "CERR_L", "%lf", 1, DL*3600.0 * field.project.cdelt1); + fits_modify (header, "CERR_M", "%lf", 1, DM*3600.0 * field.project.cdelt1); + + fits_modify (header, "CPRECISE", "%lf", 1, Cprecise); + fits_modify (header, "NSTARS", "%d", 1, 0); /*** modify addstar to allow NSTARS = 0 ***/ + fits_modify (header, "NASTRO", "%d", 1, Nastro); + + GetScatter (&Nastro, &DL, &DM, TRUE); + fits_modify (header, "CERR_LB", "%lf", 1, DL*3600.0 * field.project.cdelt1); + fits_modify (header, "CERR_MB", "%lf", 1, DM*3600.0 * field.project.cdelt1); + fits_modify (header, "NBRIGHT", "%d", 1, Nastro); + + fprintf (stderr, "bright: %f %f %d\n", DL*3600.0 * field.project.cdelt1, DM*3600.0 * field.project.cdelt1, Nastro); + + /* make a better selection for the mosaic photcode? */ + fits_scan (&chip[0].header, "PHOTCODE", "%s", 1, line); + fits_modify (header, "PHOTCODE", "%s", 1, line); + + fits_scan (&chip[0].header, "ZERO_PT", "%lf", 1, &tmp); + fits_modify (header, "ZERO_PT", "%lf", 1, tmp); + + fits_scan (&chip[0].header, ExptimeKeyword, "%lf", 1, &tmp); + fits_modify (header, ExptimeKeyword, "%lf", 1, tmp); + + /* insert time data from chip[0] */ + /* try JD first */ + if (strcasecmp (JDKeyword, "NONE")) { + uppercase (JDKeyword); + fits_scan (&chip[0].header, JDKeyword, "%lf", 1, &tmp); + fits_modify (header, JDKeyword, "%lf", 1, tmp); + goto got_time; + } + + /* try MJD next */ + if (strcasecmp (MJDKeyword, "NONE")) { + uppercase (MJDKeyword); + fits_scan (&chip[0].header, MJDKeyword, "%lf", 1, &tmp); + fits_modify (header, MJDKeyword, "%lf", 1, tmp); + goto got_time; + } + + /* get UT and DATE */ + if (strcasecmp (UTKeyword, "NONE") && strcasecmp (DateKeyword, "NONE")) { + uppercase (UTKeyword); + fits_scan (&chip[0].header, UTKeyword, "%s", 1, line); + fits_modify (header, UTKeyword, "%s", 1, line); + uppercase (DateKeyword); + fits_scan (&chip[0].header, DateKeyword, "%s", 1, line); + fits_modify (header, DateKeyword, "%s", 1, line); + goto got_time; + } + fprintf (stderr, "ERROR: missing time abstraction in config\n"); + exit (1); + +got_time: + + PutCoords (coords, header); + return (header); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkobs.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkobs.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkobs.c (revision 21807) @@ -0,0 +1,142 @@ +# include "mosastro.h" + +int main (int argc, char **argv) { + + int i, j, N, Nrefcat, Nchip, Nstars, NSTARS; + double Ro, Do, Po, Qo, Lo, Mo; + double dX, dY, Dist; + double RA, DEC; + char filename[64]; + Coords map, coords; + StarData *refcat; + SMPData *stars; + Header *header; + + ConfigInit (&argc, argv); + args (&argc, argv); + init_random (); + + RA = atof(argv[1]); + DEC = atof(argv[2]); + OUTPUT = argv[3]; + fake_field_center (RA, DEC); + + /* chip size in pixels (output is 5x5 grid of chips */ + dX = 0.1 / field.project.cdelt1; + dY = 0.1 / field.project.cdelt2; + + refcat = greference (&Nrefcat); + project_refcat (refcat, Nrefcat); + if (FOCAL_PLANE != NULL) { + FILE *f; + f = fopen (FOCAL_PLANE, "w"); + dump_stars (f, refcat, Nrefcat); + fclose (f); + } + if (NO_CHIPS) exit (0); + + Nchip = 0; + for (i = -2; i < 3; i++) { + for (j = -2; j < 3; j++) { + + /* find chip center in TP, FP and Sky */ + Po = 2*dX*i; + Qo = 2*dY*j; + RD_to_XY (&Lo, &Mo, Po, Qo, &field.distort); + XY_to_RD (&Ro, &Do, Po, Qo, &field.project); + + /* FP-Chip terms */ + strcpy (map.ctype, "DEC--WRP"); + map.crval1 = Lo; + map.crval2 = Mo; + map.crpix1 = dX; + map.crpix2 = dX; + + map.cdelt1 = 1.0; + map.cdelt2 = 1.0; + + map.pc1_1 = 1; + map.pc2_2 = 1; + map.pc1_2 = 0; + map.pc2_1 = 0; + + map.Npolyterms = 2; + for (N = 0; N < 7; N++) { + map.polyterms[N][0] = 0; + map.polyterms[N][1] = 0; + } + map.polyterms[2][1] = 1e-5; + + /* project catalog stars to chip */ + FPtoChip (refcat, Nrefcat, &map); + + Nstars = 0; + NSTARS = 1000; + ALLOCATE (stars, SMPData, NSTARS); + + /* find catalog stars on chip */ + for (N = 0; N < Nrefcat; N++) { + if (refcat[N].X < 0) continue; + if (refcat[N].Y < 0) continue; + if (refcat[N].X >= 2*dX) continue; + if (refcat[N].Y >= 2*dY) continue; + + /* add random noise - not gaussian noise */ + stars[Nstars].X = refcat[N].X + SIGMA*(drand48() - 0.5); + stars[Nstars].Y = refcat[N].Y + SIGMA*(drand48() - 0.5); + stars[Nstars].M = 16.0; + stars[Nstars].dM = 0.02; + stars[Nstars].dophot = 1; + stars[Nstars].sky = 1.0; + stars[Nstars].Mgal = 16.0; + stars[Nstars].Map = 16.0; + stars[Nstars].fx = 1.0; + stars[Nstars].fy = 1.0; + stars[Nstars].df = 0.0; + + Nstars ++; + if (Nstars >= NSTARS) { + NSTARS += 1000; + REALLOCATE (stars, SMPData, NSTARS); + } + } + + + /* Chip-Sky terms */ + strcpy (coords.ctype, "DEC--WRP"); + coords.crval1 = Ro; + coords.crval2 = Do; + coords.crpix1 = dX; + coords.crpix2 = dY; + + { + double dP, dQ, dL, dM; + double Mx, Lx, scale; + + dP = 10; + dQ = 10; + RD_to_XY (&Lx, &Mx, Po+dP, Qo+dQ, &field.distort); + dL = Lx - Lo; + dM = Mx - Mo; + scale = hypot (10.0, 10.0) / hypot (dL, dM); + coords.cdelt1 = scale/3600.0; + coords.cdelt2 = scale/3600.0; + } + + coords.pc1_1 = 1; + coords.pc2_2 = 1; + coords.pc1_2 = 0; + coords.pc2_1 = 0; + coords.Npolyterms = 1; + + header = mkheader (2*dX, 2*dY, Nstars, &coords); + + sprintf (filename, "%s.%02d.fits", OUTPUT, Nchip); + wstars (filename, stars, Nstars, header); + fits_free_header (header); + free (stars); + Nchip ++; + } + } + exit (0); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkpolyterm.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkpolyterm.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkpolyterm.c (revision 21807) @@ -0,0 +1,14 @@ +# include "mosastro.h" + +int mkpolyterm (int n, int m) { + + int i, nt, N; + + N = 0; + nt = n + m; + for (i = 2; i < nt; i++) { + N += i + 1; + } + N += m; + return (N); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkstandards.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkstandards.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mkstandards.c (revision 21807) @@ -0,0 +1,99 @@ +# include "mosastro.h" +# define DX 30.0 +# define DY 30.0 + +extern double drand48(); +void init_random (); + +/* build a grid of reference stars */ +int main (int argc, char **argv) { + + int i, N, Nstars, NSTARS, Random; + double x, y, m, dX, dY; + Header *header; + Coords coords; + SMPData *stars; + + Random = FALSE; + if ((N = get_argument (argc, argv, "-random"))) { + Random = TRUE; + remove_argument (N, &argc, argv); + Nstars = atoi (argv[N]); + remove_argument (N, &argc, argv); + } + + if (argc != 5) { + fprintf (stderr, "USAGE: mkstandards (RA) (DEC) (half-width) (output) [-random Nstars]\n"); + exit (2); + } + + /* generate grid of stars */ + dX = dY = 3600*atof (argv[3]); + + /* bore site center guess */ + strcpy (coords.ctype, "DEC--TAN"); + coords.crval1 = atof(argv[1]); + coords.crval2 = atof(argv[2]); + coords.crpix1 = dX; + coords.crpix2 = dX; + + coords.cdelt1 = 1/3600.0; + coords.cdelt2 = 1/3600.0; + + coords.pc1_1 = 1; + coords.pc2_2 = 1; + coords.pc1_2 = 0; + coords.pc2_1 = 0; + coords.Npolyterms = 1; + + if (Random) { + + init_random (); + + ALLOCATE (stars, SMPData, Nstars); + for (i = 0; i < Nstars; i++) { + x = 2*dX*drand48(); + y = 2*dY*drand48(); + m = 14.0 + 4.0*drand48(); + stars[i].X = x; + stars[i].Y = y; + stars[i].M = m; + stars[i].dM = 0.02; + } + } else { + + Nstars = 0; + NSTARS = 1000; + ALLOCATE (stars, SMPData, NSTARS); + + for (x = 0; x < 2*dX; x += DX) { + for (y = 0; y < 2*dY; y += DY) { + stars[Nstars].X = x; + stars[Nstars].Y = y; + stars[Nstars].M = 16.0; + + Nstars ++; + if (Nstars >= NSTARS) { + NSTARS += 1000; + REALLOCATE (stars, SMPData, NSTARS); + } + } + } + } + + for (i = 0; i < Nstars; i++) { + stars[i].dophot = 1; + stars[i].sky = 1.0; + stars[i].Mgal = 16.0; + stars[i].Map = 16.0; + stars[i].fx = 1.0; + stars[i].fy = 1.0; + stars[i].df = 0.0; + } + + header = mkheader (2*dX, 2*dY, Nstars, &coords); + wstars (argv[4], stars, Nstars, header); + fits_free_header (header); + + exit (0); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mosastro.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mosastro.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/mosastro.c (revision 21807) @@ -0,0 +1,97 @@ +# include "mosastro.h" +# include + +int main (int argc, char **argv) { + + int Nrefcat, Nastro; + double Cerror, DL, DM; + double Scale; + glob_t pglob; + StarData *refcat; + Gradients *grad; + + + ConfigInit (&argc, argv); + args (&argc, argv); + + if (argc != 4) { + fprintf (stderr, "USAGE: mosastro (inglob) (ext) (mosaic.phu)\n"); + exit (2); + } + + pglob.gl_offs = 0; + glob (argv[1], 0, NULL, &pglob); + + LoadStars (pglob.gl_pathc, pglob.gl_pathv); + + /* use per-chip astrometry to find ra,dec range */ + deproject_stars (); + field_stats (); /** needs coords from deproject_stars **/ + init_field (); /** needs results from field stats **/ + init_chips (); /** needs results from init_field **/ + Scale = 3600.0 * field.project.cdelt1; + + /* use field model to get TP & FP coords */ + project_stars (); + if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars(); + + refcat = greference (&Nrefcat); + project_refcat (refcat, Nrefcat); + if ((DUMP != NULL) && !strcmp (DUMP, "refcat")) dump_refcat(refcat, Nrefcat); + + match (refcat, Nrefcat); + deproject_raw (); + project_ref (); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (raw) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + if ((DUMP != NULL) && !strcmp (DUMP, "rawmatch")) dump_match(); + + grad = GetGradients (); + FitGradients (grad); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (grad) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + if ((DUMP != NULL) && !strcmp (DUMP, "fitgrads")) dump_match(); + Scale = 3600.0 * field.project.cdelt1; + + FitChips (ChipOrder); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (chip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + if ((DUMP != NULL) && !strcmp (DUMP, "fitchips_unclip")) dump_match(); + + ClipOnFP (2.5); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (clip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + + FitChips (ChipOrder); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (chip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + + ClipOnFP (2.5); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (clip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + + FitChips (ChipOrder); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (chip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + + ClipOnFP (2.5); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (clip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + + FitChips (ChipOrder); + Cerror = GetScatter (&Nastro, &DL, &DM, FALSE); + fprintf (stderr, "scatter (brit) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + Cerror = GetScatter (&Nastro, &DL, &DM, TRUE); + fprintf (stderr, "scatter (chip) : %6.4f for %d stars (%6.4f, %6.4f)\n", Cerror, Nastro, DL*Scale, DM*Scale); + if ((DUMP != NULL) && !strcmp (DUMP, "fitchips")) dump_match(); + + /* testcoords (); */ + output (argv[2], argv[3]); + exit (0); + +/* + FitField (3); + FitChips (3); +*/ + +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/output.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/output.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/output.c (revision 21807) @@ -0,0 +1,45 @@ +# include "mosastro.h" + +void output (char *ext, char *phu) { + + int i, N; + char *p, outname[256]; + FILE *f; + Header *header; + + for (i = 0; i < Nchip; i++) { + + /* convert chip[i].file to output name */ + + p = strrchr (chip[i].file, '.'); + N = p - chip[i].file + 1; + if (p == NULL) N = strlen (chip[i].file); + + bzero (outname, 256); + strncpy (outname, chip[i].file, N); + + strcat (outname, ext); + + /* convert chip from PLY to WRP (linked to field model) */ + strcpy (chip[i].map.ctype, "DEC--WRP"); + PutCoords (&chip[i].map, &chip[i].header); + wchip (outname, &chip[i]); + } + + field_combine (); + + f = fopen (phu, "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file %s\n", phu); + exit (1); + } + + header = mkmosaic (1, 1, 0, &field.project); + if (SAVE_RESID) { + SaveResiduals (f, header); + return; + } + + fwrite (header[0].buffer, 1, header[0].size, f); + fclose (f); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/parse_time.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/parse_time.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/parse_time.c (revision 21807) @@ -0,0 +1,120 @@ +# include "mosastro.h" + +e_time parse_time (Header *header) { + + double jd; + int Ny, Nf, mode; + int hour, min, sec, year, month, day; + char *py, *pm, *pd, *c; + char line[256]; + e_time Nsec; + + /* 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; + } + + /* 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/ohana/mosastro-1-0/Ohana/src/mosastro/src/project.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/project.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/project.c (revision 21807) @@ -0,0 +1,54 @@ +# include "mosastro.h" + +int project_refcat (StarData *refcat, int Nrefcat) { + + SkyToTP (refcat, Nrefcat, &field.project); + TPtoFP (refcat, Nrefcat, &field.distort); + + return (1); +} + +int deproject_stars () { + + int i; + + for (i = 0; i < Nchip; i++) { + ChipToSky (chip[i].stars, chip[i].Nstars, &chip[i].coords); + } + return (1); +} + +int project_stars () { + + int i; + + for (i = 0; i < Nchip; i++) { + SkyToTP (chip[i].stars, chip[i].Nstars, &field.project); + TPtoFP (chip[i].stars, chip[i].Nstars, &field.distort); + } + return (1); +} + +int deproject_raw () { + + int i; + + for (i = 0; i < Nchip; i++) { + ChipToFP (chip[i].raw, chip[i].Nmatch, &chip[i].map); + FPtoTP (chip[i].raw, chip[i].Nmatch, &field.distort); + TPtoSky (chip[i].raw, chip[i].Nmatch, &field.project); + } + return (1); +} + +int project_ref () { + + int i; + + for (i = 0; i < Nchip; i++) { + SkyToTP (chip[i].ref, chip[i].Nmatch, &field.project); + TPtoFP (chip[i].ref, chip[i].Nmatch, &field.distort); + FPtoChip (chip[i].ref, chip[i].Nmatch, &chip[i].map); + } + return (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/random.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/random.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/random.c (revision 21807) @@ -0,0 +1,12 @@ +# include +# include +# include + +void init_random () { + + long A, B; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rfits.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rfits.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rfits.c (revision 21807) @@ -0,0 +1,58 @@ +# include "mosastro.h" + +int rfits (Chip *mychip) { + + int i, Nx; + FILE *f; + FTable table; + SMPData *stars; + + /* open file for stars */ + f = fopen (mychip[0].file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't open file to load stars\n"); + exit (1); + } + fseek (f, mychip[0].header.size, SEEK_SET); + + /* init & load in table data */ + table.header = &mychip[0].theader; + if (!fits_fread_matrix (f, &mychip[0].matrix, &mychip[0].header)) { + fprintf (stderr, "error reading file\n"); + fclose (f); + return (FALSE); + } + if (!fits_fread_ftable (f, &table, "SMPFILE")) { + fprintf (stderr, "error reading file\n"); + fits_free_matrix (&mychip[0].matrix); + fclose (f); + return (FALSE); + } + + stars = fits_table_get_SMPData (&table, &mychip[0].Nstars, NULL); + fits_scan (table.header, "NAXIS1", "%d", 1, &Nx); + + /* save raw data for output */ + mychip[0].FITS = TRUE; + mychip[0].buffer = (char *) stars; + mychip[0].Nbuffer = Nx*mychip[0].Nstars; + + if (mychip[0].Nstars < 5) { + fprintf (stderr, "Too few stars for reliable solution, only %d\n", mychip[0].Nstars); + fits_free_matrix (&mychip[0].matrix); + fits_free_table (&table); + fclose (f); + return (FALSE); + } + + /* note the different structures for mychip[0].stars and stars */ + ALLOCATE (mychip[0].stars, StarData, mychip[0].Nstars); + for (i = 0; i < mychip[0].Nstars; i++) { + mychip[0].stars[i].X = stars[i].X; + mychip[0].stars[i].Y = stars[i].Y; + mychip[0].stars[i].Mag = stars[i].M; + mychip[0].stars[i].dMag = stars[i].dM; + } + return (TRUE); +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rtext.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rtext.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/rtext.c (revision 21807) @@ -0,0 +1,40 @@ +# include "mosastro.h" +# define BYTES_STAR 66 + +int rtext (Chip *mychip) { + + int i, Nbytes, nbytes; + FILE *f; + + f = fopen (mychip[0].file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't read data from %s\n", mychip[0].file); + return (FALSE); + } + fseek (f, mychip[0].header.size, SEEK_SET); + + /* find expected number of stars */ + fits_scan (&mychip[0].header, "NSTARS", "%d", 1, &mychip[0].Nstars); + if (mychip[0].Nstars == 0) { + fprintf (stderr, "ERROR: can't get NSTARS from header for %s\n", mychip[0].file); + return (FALSE); + } + ALLOCATE (mychip[0].stars, StarData, mychip[0].Nstars); + + Nbytes = BYTES_STAR * mychip[0].Nstars; + ALLOCATE (mychip[0].buffer, char, Nbytes + 1); + nbytes = fread (mychip[0].buffer, 1, Nbytes, f); + if (nbytes != Nbytes) { exit (1); } + mychip[0].Nbuffer = Nbytes; + mychip[0].FITS = FALSE; + + for (i = 0; i < mychip[0].Nstars; i++) { + bzero (&mychip[0].stars[i], sizeof(StarData)); + dparse (&mychip[0].stars[i].X, 1, &mychip[0].buffer[i*BYTES_STAR]); + dparse (&mychip[0].stars[i].Y, 2, &mychip[0].buffer[i*BYTES_STAR]); + dparse (&mychip[0].stars[i].Mag, 3, &mychip[0].buffer[i*BYTES_STAR]); + dparse (&mychip[0].stars[i].dMag, 4, &mychip[0].buffer[i*BYTES_STAR]); + mychip[0].stars[i].dMag *= 0.001; /* millimag errors stored in file */ + } + return (TRUE); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/sort.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/sort.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/sort.c (revision 21807) @@ -0,0 +1,50 @@ +# include "mosastro.h" + +int sortthree (double *X, double *Y, int *Z, int N) { + + int l,j,ir,i; + double tX, tY; + int tZ; + + if (N < 2) return (1); + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tZ = Z[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tZ = Z[ir]; + Z[ir] = Z[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + Z[0] = tZ; + return (1); + } + } + 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]; + Z[i] = Z[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + Z[i] = tZ; + } + return (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/testcoords.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/testcoords.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/testcoords.c (revision 21807) @@ -0,0 +1,66 @@ +# include "mosastro.h" + +void testcoords () { + + int i; + FILE *f; + double R, D, P, Q, L, M, X, Y; + + fprintf (stderr, "starting test.dat\n"); + + /* version 1: transform with three separate stages */ + /* X,Y -> L,M -> P,Q -> R,D */ + f = fopen ("test.1.dat", "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file \n"); + exit (1); + } + for (i = 0; i < chip[0].Nmatch; i++) { + X = chip[0].raw[i].X; + Y = chip[0].raw[i].Y; + XY_to_RD (&L, &M, X, Y, &chip[0].map); + XY_to_RD (&P, &Q, L, M, &field.distort); + XY_to_RD (&R, &D, P, Q, &field.project); + fprintf (f, "%4d %10.6f %10.6f %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f\n", + i, R, D, P, Q, L, M, X, Y); + } + fclose (f); + + /* version 2: transform with two separate stages */ + /* X,Y -> L,M -> R,D */ + f = fopen ("test.2.dat", "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file \n"); + exit (1); + } + field_combine (); + for (i = 0; i < chip[0].Nmatch; i++) { + X = chip[0].raw[i].X; + Y = chip[0].raw[i].Y; + XY_to_RD (&L, &M, X, Y, &chip[0].map); + XY_to_RD (&R, &D, L, M, &field.project); + fprintf (f, "%4d %10.6f %10.6f %8.2f %8.2f %8.2f %8.2f\n", + i, R, D, L, M, X, Y); + } + fclose (f); + + + /* version 3: transform with automatic two stages in coordops */ + /* X,Y -> R,D */ + f = fopen ("test.3.dat", "w"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file \n"); + exit (1); + } + RegisterMosaic (&field.project); + strcpy (chip[0].map.ctype, "DEC--WRP"); + for (i = 0; i < chip[0].Nmatch; i++) { + X = chip[0].raw[i].X; + Y = chip[0].raw[i].Y; + XY_to_RD (&R, &D, X, Y, &chip[0].map); + fprintf (f, "%4d %10.6f %10.6f %8.2f %8.2f\n", + i, R, D, X, Y); + } + fclose (f); + exit (1); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/transforms.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/transforms.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/transforms.c (revision 21807) @@ -0,0 +1,54 @@ +# include "mosastro.h" + +void ChipToSky (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + XY_to_RD (&stars[i].R, &stars[i].D, stars[i].X, stars[i].Y, coords); + while (stars[i].R < 0.0) stars[i].R += 360.0; + while (stars[i].R >= 360.0) stars[i].R -= 360.0; + } +} + +void ChipToFP (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + XY_to_RD (&stars[i].L, &stars[i].M, stars[i].X, stars[i].Y, coords); + } +} + +void FPtoTP (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + XY_to_RD (&stars[i].P, &stars[i].Q, stars[i].L, stars[i].M, coords); + } +} + +void TPtoSky (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + XY_to_RD (&stars[i].R, &stars[i].D, stars[i].P, stars[i].Q, coords); + while (stars[i].R < 0.0) stars[i].R += 360.0; + while (stars[i].R >= 360.0) stars[i].R -= 360.0; + } +} + +void SkyToTP (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + RD_to_XY (&stars[i].P, &stars[i].Q, stars[i].R, stars[i].D, coords); + } +} + +void TPtoFP (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + RD_to_XY (&stars[i].L, &stars[i].M, stars[i].P, stars[i].Q, coords); + } +} + +void FPtoChip (StarData *stars, int Nstars, Coords *coords) { + int i; + for (i = 0; i < Nstars; i++) { + RD_to_XY (&stars[i].X, &stars[i].Y, stars[i].L, stars[i].M, coords); + } +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/warptest.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/warptest.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/warptest.c (revision 21807) @@ -0,0 +1,55 @@ +# include "mosastro.h" +int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords); +int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords); + +main (int argc, char **argv) { + + int i, status; + double X, Y, x, y, p, q; + Coords coords; + + if (argc != 2) { + fprintf (stderr, "USAGE: warptest (system)\n"); + exit (2); + } + + /* bore site center guess */ + sprintf (coords.ctype, "RA---%s", argv[1]); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + coords.crpix1 = 0.0; + coords.crpix2 = 0.0; + coords.cdelt1 = 1.0; + coords.cdelt2 = 1.0; + + /* + coords.cdelt1 = 1.0/3600.0; + coords.cdelt2 = 1.0/3600.0; + */ + + /** allow guess at field rotation?? **/ + coords.pc1_1 = 1; + coords.pc2_2 = 1; + coords.pc1_2 = 0; + coords.pc2_1 = 0; + + /* allow 2nd and 3rd order? */ + /* how do we handle renormalization? (fixed at 1000 pixels??) */ + coords.Npolyterms = 3; + for (i = 0; i < 7; i++) { + coords.polyterms[i][0] = 0; + coords.polyterms[i][1] = 0; + } + coords.polyterms[3][0] = 1e-10; + coords.polyterms[6][1] = 1e-10; + + for (x = -10000; x <= 10000; x+= 500) { + for (y = -10000; y <= 10000; y+= 500) { + status = XY_to_RD (&p, &q, x, y, &coords); + status = RD_to_XY (&X, &Y, p, q, &coords); + fprintf (stdout, "%f %f %f %f %f %f\n", x, y, p, q, X, Y); + if (!status) exit (1); + } + } + exit (0); +} Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wfits.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wfits.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wfits.c (revision 21807) @@ -0,0 +1,29 @@ +# include "mosastro.h" + +void wfits (char *filename, SMPData *stars, int Nstars, Header *header) { + + Matrix matrix; + Header theader; + FTable table; + + header[0].extend = TRUE; + header[0].Naxes = 0; + fits_modify (header, "NAXIS", "%d", 1, 0); + fits_modify (header, "EXTEND", "%t", 1, TRUE); + fits_modify (header, "NEXTEND", "%d", 1, 1); + + /* add in some keywords to specify the datatype & software version? */ + + /* create (empty) data matrix */ + fits_create_matrix (header, &matrix); + + table.header = &theader; + fits_table_set_SMPData (&table, stars, Nstars); + + fits_write_header (filename, header); + fits_write_matrix (filename, &matrix); + fits_write_Theader (filename, &theader); + fits_write_table (filename, &table); + return; +} + Index: /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wstars.c =================================================================== --- /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wstars.c (revision 21807) +++ /tags/ohana/mosastro-1-0/Ohana/src/mosastro/src/wstars.c (revision 21807) @@ -0,0 +1,63 @@ +# include "mosastro.h" +# define NCHAR 66 /* 65 char EXCLUDING return */ + +void wchip (char *filename, Chip *data) { + + FILE *g; + + if (data[0].FITS) { + wfits (filename, (SMPData *) data[0].buffer, data[0].Nstars, &data[0].header); + } else { + + g = fopen (filename, "w"); + if (g == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file %s\n", filename); + exit (1); + } + + fwrite (data[0].header.buffer, 1, data[0].header.size, g); + fwrite (data[0].buffer, 1, data[0].Nbuffer, g); + fclose (g); + } +} + +void wstars (char *filename, SMPData *stars, int Nstars, Header *header) { + + int i, Nchar; + FILE *g; + char line[NCHAR + 3]; + + g = fopen (filename, "w"); + if (g == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't create output file %s\n", filename); + exit (1); + } + + fwrite (header[0].buffer, 1, header[0].size, g); + + for (i = 0; i < Nstars; i++) { + if (stars[i].dophot == 0) continue; + Nchar = snprintf (line, NCHAR, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f", + stars[i].X, stars[i].Y, stars[i].M, + (int)(1000*stars[i].dM), stars[i].dophot, stars[i].sky, + stars[i].Mgal, stars[i].Map, stars[i].fx, stars[i].fy, stars[i].df); + + /* this is just a little funny. NCHAR (in) includes the trailing NULL, Nchar (out) excludes it */ + if (Nchar != NCHAR - 1) { + fprintf (stderr, "funny line %d (%d)\n%s\n", i, Nchar, line); + } else { + fprintf (g, "%s\n", line); + } + } + + fclose (g); + +} + +/* + + 63.6 2869.5 17.568 157 17.568 17.568 25.01 25.00 360.0 7 2.9 + + 63.6 2869.5 17.568 157 7 2.9 + +*/