Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/Makefile =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/Makefile (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/Makefile (revision 21919) @@ -0,0 +1,63 @@ +default: relastro +help: + @echo "make options: relastro (default)" + +include ../../Makefile.System +HOME = $(ROOT)/src/relastro +BIN = $(HOME)/bin +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +include ../../Makefile.Common + +# programs may add their own internal requirements here +FULL_CFLAGS = $(BASE_CFLAGS) +FULL_CPPFLAGS = $(BASE_CPPFLAGS) +FULL_LDFLAGS = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS) + +relastro: $(BIN)/relastro.$(ARCH) +install: $(DESTBIN)/relastro + +RELASTRO = \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/FitChip.$(ARCH).o \ +$(SRC)/FitMosaic.$(ARCH).o \ +$(SRC)/FitPM.$(ARCH).o \ +$(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitSimple.$(ARCH).o \ +$(SRC)/ImageOps.$(ARCH).o \ +$(SRC)/MosaicOps.$(ARCH).o \ +$(SRC)/ParFactor.$(ARCH).o \ +$(SRC)/SetSignals.$(ARCH).o \ +$(SRC)/Shutdown.$(ARCH).o \ +$(SRC)/UpdateChips.$(ARCH).o \ +$(SRC)/UpdateMosaic.$(ARCH).o \ +$(SRC)/UpdateObjects.$(ARCH).o \ +$(SRC)/UpdateSimple.$(ARCH).o \ +$(SRC)/args.$(ARCH).o \ +$(SRC)/bcatalog.$(ARCH).o \ +$(SRC)/dvo_astrom_ops.$(ARCH).o \ +$(SRC)/fitpoly.$(ARCH).o \ +$(SRC)/initialize.$(ARCH).o \ +$(SRC)/liststats.$(ARCH).o \ +$(SRC)/load_catalogs.$(ARCH).o \ +$(SRC)/load_images.$(ARCH).o \ +$(SRC)/mkpolyterm.$(ARCH).o \ +$(SRC)/plot_scatter.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/select_images.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o \ +$(SRC)/relastro.$(ARCH).o \ +$(SRC)/reload_catalogs.$(ARCH).o \ +$(SRC)/save_catalogs.$(ARCH).o \ +$(SRC)/setExclusions.$(ARCH).o \ +$(SRC)/write_coords.$(ARCH).o + +DROP = \ +$(SRC)/setMrelFinal.$(ARCH).o \ +$(SRC)/StarOps.$(ARCH).o + + +$(RELASTRO): $(INC)/relastro.h +$(BIN)/relastro.$(ARCH): $(RELASTRO) Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/bin/.cvsignore =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/bin/.cvsignore (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/bin/.cvsignore (revision 21919) @@ -0,0 +1,2 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix +*.linrh Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/notes.txt =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/notes.txt (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/notes.txt (revision 21919) @@ -0,0 +1,64 @@ + +2007.02.11 + + relastro major modes: + + - update the astrometry of objects in the images + - load objects within a region + - foreach object + - calculate average R,D + - update db tables (dR, dD) + - include external trends + - parallax factor + - atm trends + + - update the astrometric parameters of images + - use the average R,D for objects, update image terms + * include external refs with adjustable weighting + + - image parameter smoothing + + - simultaneous fit of objects and images + - identify parameters to constrain + - limited number of internal objects? + + +2006.04.08 : relastro contemplations + +relastro will perform the astrometry equivalent to relphot. For +region, it loads in the average objects, detections, and image +parameters. The goal is to determine improved astrometric positions +of objects in the images and to determine improved astrometric +parameters for the images. Just like relphot, this is an iterative +process in which the object postions are improved in one step, the +image parameters are improved in the next, then back to the images. + +Relastro needs to deal with outlier objects in a robust way. There +are a few reasons why objects should have poor astrometric solutions, +and these can be flagged up front: saturation and too many bad pixels. +For detections which come from psphot, we can use the quality +parameter to exclude objects which have excessively high bad pixel +counts. For detections from sextractor or other analysis tools, +this information is missing. In this case, we can use the distance +from the image boundaries as a trigger. + +An open question to be addressed is the issue of external calibration +sources. Unlike photometry, external astrometric measurements do not +need to have an unknown filter transformation applied. They can thus +be included as part of the solution automatically. It is not clear, +however, how to include *proper-motion* data in the analysis. If I +introduce, for example, USNO objects with proper motions, the +information provided consists in practice of a velocity vector and +errors. How does this get folded into the fit? One answer is to +construct additional artificial datapoints for each entry with the +proper errors and epochs. Seems like something of a hack. + +We may want to allow for different errors / weightings for the +different external sources. + +After the initial iterations, we will can solve for proper-motion and +parallax for (all?) objects + +- calculate the parallax normal vector for each image +- calculate the atmospheric / orbital motion effects for each image + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/parallax.pro =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/parallax.pro (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/doc/parallax.pro (revision 21919) @@ -0,0 +1,348 @@ + +macro sample + PMsim 15 30 20 10 $ANGLE + set dX = zero(X) + 10 + set dY = zero(Y) + 10 + + style -c black + plot -x 2 -pt 2 -sz 2.0 X Y -dx dX -dy dY + + # set x = (1000.0/$dS)*cos(t*$TP) + $vX*t + # set y = (1000.0/$dS)*sin(t*$TP)*dsin($ANGLE) + $vY*t + # plot -c grey60 -x 2 -sz 2.0 -pt 7 x y + + create t 0 1.5 0.01 + set x = (1000.0/$dS)*cos(t*$TP) + $vX*t + set y = (1000.0/$dS)*sin(t*$TP)*dsin($ANGLE) + $vY*t + plot -c grey60 -x 0 x y + + # label -x "offset east (milliarcsec)" + # label -y "offset north (milliarcsec)" +end + +macro PMsim + if ($0 != 6) + echo "USAGE: PMsim Do vX vY dP dT" + echo "Do : distance in parsec" + echo "vX, vY : velocity (milliarcsec/year)" + echo "dP : astrometric scatter (arcsec)" + echo "dT : ecliptic lattitude (degrees)" + break + end + + $TP = 2*3.1416 + + $dS = $1 + $vX = $2 + $vY = $3 + $dP = $4 + $dT = $5 + + mkbase + + set x = (1000.0/$dS)*cos(t*$TP) + $vX*t + set y = (1000.0/$dS)*sin(t*$TP)*dsin($dT) + $vY*t + + gaussdev dx x[] 0.0 $dP + gaussdev dy y[] 0.0 $dP + + set X = x + dx + set Y = y + dy + + if ($PLOT) + lim -200 200 -200 200; clear; box; + style -c black + plot -x 2 -pt 2 -sz 2.0 X Y + plot -x 0 X Y + end + + PMfit X Y t $dT + + if ($PLOT) + label -x "offset east (milliarcsec)" + label -y "offset north (milliarcsec)" + section a 0 0 1 1 + lim 0 1 0 1 + sprintf line "input: D: %5.1f (pc), V: %5.1f, %5.1f (mas/yr), dS: %4.1f" $1 $2 $3 $4 + textline -fn courier 14 0.15 0.95 "$line" + sprintf line "fit: D: %5.1f (pc), V: %5.1f, %5.1f (mas/yr)" $dS $vX $vY + textline -fn courier 14 0.15 0.90 "$line" + section default + end +end + +macro PMfit + if ($0 != 5) + echo "USAGE: PMfit X Y t (theta)" + echo "X,Y are offsets in milliarcsec" + echo "t is time in years" + break + end + + $TP = 2*3.1416 + set x = $1 + set y = $2 + set t = $3 + set cs = cos(t*$TP) + set snx = sin(t*$TP)*dsin($4) + + mcreate A 3 3 + create B 0 3 + + # define A values: + set tmp = cs^2 + snx^2 + vstat -q tmp + zap A 0 0 1 1 -v $TOTAL + + set tmp = t*cs + vstat -q tmp + zap A 1 0 1 1 -v $TOTAL + zap A 0 1 1 1 -v $TOTAL + + set tmp = t*snx + vstat -q tmp + zap A 2 0 1 1 -v $TOTAL + zap A 0 2 1 1 -v $TOTAL + + set tmp = t^2 + vstat -q tmp + zap A 1 1 1 1 -v $TOTAL + zap A 2 2 1 1 -v $TOTAL + + zap A 2 1 1 1 -v 0 + zap A 1 2 1 1 -v 0 + + # define B values: + set tmp = x*cs + y*snx + vstat -q tmp + B[0] = $TOTAL + + set tmp = x*t + vstat -q tmp + B[1] = $TOTAL + + set tmp = y*t + vstat -q tmp + B[2] = $TOTAL + + gaussj A B + + $dS = 1000/B[0] + $vX = B[1] + $vY = B[2] + + if ($PLOT) + echo "Do: $dS" + echo "vX: $vX" + echo "vY: $vY" + + set x = (1000.0/$dS)*cos(t*$TP) + $vX*t + set y = (1000.0/$dS)*sin(t*$TP)*dsin($4) + $vY*t + plot -c red -x 2 -sz 2.0 -pt 7 x y + end +end + +macro PMstats + + if ($0 != 7) + echo "USAGE: PMsim Do vX vY dP dT (Niter)" + echo "Do : distance in parsec" + echo "vX, vY : velocity (milliarcsec/year)" + echo "dP : astrometric scatter (arcsec)" + echo "dT : ecliptic lattitude (degrees)" + break + end + + delete -q vx vy ds + $PLOT = 0 + + for i 0 $6 + PMsim $1 $2 $3 $4 $5 + concat $dS ds + concat $vX vx + concat $vY vy + end + + vstat ds + vstat vx + vstat vy + + set gm = 1000/ds + vstat gm + + clear -s -n 1 + label -fn courier 14 + # section a 0.0 0.0 0.5 1.0 + lim {$2-30} {$2+30} {$3-30} {$3+30}; box + plot -x 2 -pt 2 -sz 0.5 vx vy + label -x "p.m. (east, mas/yr)" + label -y "p.m. (east, mas/yr)" + + clear -s -n 2 + # section b 0.5 0.0 0.5 1.0 + histogram gm Ng {1000/$1 - 30.0} {1000/$1 + 30.0} 1.0 + create dg {1000/$1 - 30.0} {1000/$1 + 30.0} 1.0 + lim dg Ng; box + plot -x 1 dg Ng + + label -x "parallax (mas)" + label -y "# of tests" + +end + +macro PMtrend + + PMstats 10 5.0 5.0 10 90 100 + echo 10 + PMstats 10 5.0 5.0 {750/2.38/25} 90 100 + echo {750/2.38/25} + PMstats 10 5.0 5.0 {750/2.38/10} 90 100 + echo {750/2.38/10} + PMstats 10 5.0 5.0 {750/2.38/5} 90 100 + echo {750/2.38/5} +end + +list dist + 5 + 10 + 15 + 20 + 30 + 40 + 50 + 60 + 70 + 80 +end + +macro ParDetect + + local i + delete -q xp yp dyp dym + + $ANGLE = 45 + + for i 0 $dist:n + PMstats $dist:$i 5.0 5.0 10 $ANGLE 100 + concat $dist:$i xp + concat {1000/$MEAN} yp + concat {1000/($MEAN - $SIGMA) - 1000/$MEAN} dyp + concat {1000/$MEAN - 1000/($MEAN + $SIGMA)} dym + end + + section a 0.0 0.0 0.5 1.0 + + lim xp -1.0 101; clear; box; + plot -c black -x 2 xp yp -dy dym +dy dyp + + delete -q xp yp + concat 0 xp + concat 100 xp + set yp = xp + plot -x 0 -c blue xp yp + + set yp = zero(xp) + 62.5 + plot -x 0 -c red -lt 1 xp yp + style -lt 0 + label -x "input dist (pc)" -y "output dist" -fn courier 14 + + $PLOT = 0 + section b 0.5 0.0 0.5 1.0 + lim -120 120 -120 120; box -labels 1001; + + for i 0 5 + sample + end + + label -x "offset east (milliarcsec)" + label +y "offset north" + ps -name ParallaxDist.ps +end + +macro fields + + $Dx = 0.25*4800*8/3600 + $Fo = 1.09 + + $dD = 2.65/1.09 + + $Nf = 0 + for i {-30+$dD/2} 90 $dD + $Np = int(dcos($i)*1.09*360/2.65) + 1 + $Nf = $Nf + $Np + fprintf "%5.1f %5.1f %4d %4d" $i {dcos($i)*1.09*360/2.65} $Np $Nf + end +end + +macro SNsim + $w = 0.76 + $m1 = 24.8 + $mu = 20.4 + + lim 13.5 26.5 0.0 4.0; clear; box + + $t = 1 + SN + + $t = 5 + SN + + $t = 30 + SN + + $t = 100 + SN + + set sn = 5 + zero (m) + set lsn = log(sn) + plot -lt 2 m lsn + + set sn = 25 + zero (m) + set lsn = log(sn) + plot -lt 2 m lsn + + set sn = 100 + zero (m) + set lsn = log(sn) + plot -lt 2 m lsn + style -lt 0 +end + +macro SN + create m 14 26 0.1 + set sn = ten (-0.2*(2*m - $mu - $m1)) * sqrt($t/$w^2/3.1416) + set lsn = log(sn) + plot m lsn +end + +# r: 30s @ 5sig = 22.65 (0.6") +# i: 30s @ 5sig = 22.42 (0.6") +# i: 30s @ 5sig = 22.27 (for 0.3 mag lower zp) + +macro mkbase + delete -q t + concat 0.0 t + concat {10.0/365.0} t + concat {20.0/365.0} t + concat { 6.0/52.0} t + + concat { 3.0/12.0} t + concat { 6.0/12.0} t + concat {12.0/12.0} t + concat {18.0/12.0} t +end + + +macro parfactor + + $RA = $1 + $Dec = $2 + $e = 23 + 27/60 + + create s 0 360 + set pR = dcos($e)*d + + set pR = +(dcos($e)*dsin(s)*dcos($RA) - dcos(s)*dsin($RA)) + set pD = -(dcos($e)*dsin(s)*dsin($RA) + dcos(s)*dcos($RA))*dsin($Dec) + dsin($e)*dsin(s)*dcos($Dec) + + lim -1.1 1.1 -1.1 1.1; clear; box; plot pR pD +end Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/include/relastro.h =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/include/relastro.h (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/include/relastro.h (revision 21919) @@ -0,0 +1,287 @@ +# include +# include +# 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; + +// structure to hold coordinate fitting terms +typedef struct { + int Npts; + int Nterms; + int Norder; + int Nsums; + int Nelems; + double **sum; + double **xsum; + double **ysum; + double **xfit; + double **yfit; +} CoordFit; + +typedef struct { + double Ro, dRo; + double Do, dDo; + + double uR, duR; + double uD, duD; + + double p, dp; + + double chisq; + int Nfit; +} PMFit; + +typedef struct { + unsigned int start; + unsigned int stop; + float Mcal; + float dMcal; + short Xm; + float secz; + char code; + Coords coords; +} Mosaic; + +typedef struct { + double median; + double mean; + double sigma; + double error; + double chisq; + double min; + double max; + double total; + int Nmeas; +} StatType; + +/* global variables set in parameter file */ +char ImageCat[256]; +char ImageTemplate[256]; +char CatTemplate[256]; +char GSCFILE[256]; +char CATDIR[256]; +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char CameraConfig[256]; +char SKY_TABLE[256]; +int SKY_DEPTH; /** XXX EAM : depth of catalog tables, fix usage */ + +double MAG_LIM; +double SIGMA_LIM; +double IMAGE_SCATTER; +double IMAGE_OFFSET; +double STAR_SCATTER; +double STAR_CHISQ; +double MIN_ERROR; + +int VERBOSE; + +int NLOOP; +int RESET; +int UPDATE; +int PLOTSTUFF; +int SAVEPLOT; +int SHOW_PARAMS; +char MOSAICNAME[256]; +char STATMODE[32]; +int STAR_BAD; +int MEAS_BAD; +int POS_TOOFEW; +int PM_TOOFEW; +double PM_DT_MIN; +int IMAGE_TOOFEW; +double IMAGE_GOOD_FRACTION; +int IMAGE_BAD; +int FREEZE_IMAGES; +int USE_GRID; +int PLOTDELAY; + +int RELPHOT_GRID_X; +int RELPHOT_GRID_Y; +int RELPHOT_GRID_BINNING; + +char *PHOTCODE_LIST; +int Nphotcodes; +PhotCode **photcodes; +int PhotPrimary; +int PhotNsec; + +int AreaSelect; +double AreaXmin, AreaXmax, AreaYmin, AreaYmax; + +int ImagSelect, ImagMin, ImagMax; + +double PlotMmin, PlotMmax, PlotdMmin, PlotdMmax; +enum {black, white, red, orange, yellow, green, blue, indigo, violet}; + +int DophotSelect, DophotValue; + +int TimeSelect; +time_t TSTART, TSTOP; + +int FIT_MODE; +enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR}; + +int FIT_TARGET; +enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS}; + +SkyRegion UserPatch; + +int DoUpdateObjects; +int DoUpdateSimple; +int DoUpdateChips; +int DoUpdateMosaics; + +/*** relphot prototypes ***/ +void ConfigInit PROTO((int *argc, char **argv)); +void GetConfig PROTO((char *config, char *field, char *format, int N, void *ptr)); +char *GetPhotnamebyCode PROTO((PhotCodeData *photcodes, int code)); +void InterpolateGrid PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords)); +int *SelectRefMosaic PROTO((Mosaic **refmosaic, int *Nimage)); +int args PROTO((int argc, char **argv)); +int bcatalog PROTO((Catalog *subcatalog, Catalog *catalog)); +void clean_images PROTO(()); +void clean_measures PROTO((Catalog *catalog, int Ncatalog, int final)); +void clean_mosaics PROTO(()); +void clean_stars PROTO((Catalog *catalog, int Ncatalog)); +int corner_check PROTO((double *x1, double *y1, double *x2, double *y2)); +void dumpGrid PROTO(()); +void dump_grid PROTO(()); +int edge_check PROTO((double *x1, double *y1, double *x2, double *y2)); +void findImages PROTO((Catalog *catalog, int Ncatalog)); +int findMosaics PROTO((Catalog *catalog, int Ncatalog)); +Image *find_images PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum)); +void set_db (FITS_DB *in); +int Shutdown (char *format, ...); +void TrapSignal (int sig); +void SetProtect (int mode); +int SetSignals (); + +GSCRegion *find_regions PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion)); +void freeGridBins PROTO((int Ncatalog)); +void freeImageBins PROTO((int Ncatalog)); +void freeMosaicBins PROTO((int Ncatalog)); +void free_catalogs PROTO((Catalog *catalog, int Ncatalog)); +int gcatalog PROTO((Catalog *catalog, int FINAL)); +Coords *getCoords PROTO((int meas, int cat)); +float getMcal PROTO((int meas, int cat)); +float getMgrid PROTO((int meas, int cat)); +float getMmos PROTO((int meas, int cat)); +float getMrel PROTO((Catalog *catalog, int meas, int cat)); +GSCRegion *get_regions PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions)); +void getfullregion PROTO((Image *image, int Nimage, GSCRegion *fullregion)); +Image *getimage PROTO((int N)); +Image *getimages PROTO((int *N)); +void global_stats PROTO((Catalog *catalog, int Ncatalog)); +void initGrid PROTO((int dX, int dY)); +void initGridBins PROTO((Catalog *catalog, int Ncatalog)); +void initImageBins PROTO((Catalog *catalog, int Ncatalog)); +void initImages PROTO((Image *input, int N)); +void initMosaicBins PROTO((Catalog *catalog, int Ncatalog)); +void initMosaicGrid PROTO((Image *image, int Nimage)); +void initMosaics PROTO((Image *image, int Nimage)); +void initMrel PROTO((Catalog *catalog, int Ncatalog)); +void initialize PROTO((int argc, char **argv)); +void initstats PROTO((char *mode)); +int liststats PROTO((double *value, double *dvalue, int N, StatType *stats)); +Catalog *load_catalogs PROTO((SkyList *skylist, int *Ncatalog, int subselect)); +SkyList *load_images PROTO((FITS_DB *db, SkyRegion *region)); +Image *select_images PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage)); + +void check_permissions (char *basefile); +void lock_image_db (FITS_DB *db, char *filename); +void unlock_image_db (FITS_DB *db); +void create_image_db (FITS_DB *db); + +int main PROTO((int argc, char **argv)); +void mark_images PROTO((Image *image, int Nimage, Image *timage, int Ntimage)); +void matchImage PROTO((Catalog *catalog, int meas, int cat)); +void matchMosaics PROTO((Catalog *catalog, int meas, int cat)); +GSCRegion *name_region PROTO((char *name, int *Nregions)); +double opening_angle PROTO((double x1, double y1, double x2, double y2, double x3, double y3)); +void plot_chisq PROTO((Catalog *catalog, int Ncatalog)); +void plot_defaults PROTO((Graphdata *graphdata)); +void plot_grid PROTO((Catalog *catalog)); +void plot_images PROTO(()); +void plot_list PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file)); +void plot_mosaic_fields PROTO((Catalog *catalog)); +void plot_mosaics PROTO(()); +void plot_scatter PROTO((Catalog *catalog, int Ncatalog)); +void plot_star_coords PROTO((Catalog *catalog, int Ncatalog)); +void plot_stars PROTO((Catalog *catalog, int Ncatalog)); +void reload_catalogs PROTO((SkyList *skylist)); +int setExclusions PROTO((Catalog *catalog, int Ncatalog)); +void setMcal PROTO((Catalog *catalog, int Poor)); +void setMcalFinal PROTO(()); +int setMcalOutput PROTO((Catalog *catalog, int Ncatalog)); +void setMgrid PROTO((Catalog *catalog)); +int setMmos PROTO((Catalog *catalog, int Poor)); +int setMrel PROTO((Catalog *catalog, int Ncatalog)); +void setMrelFinal PROTO((Catalog *catalog)); +int setMrelOutput PROTO((Catalog *catalog, int Ncatalog, int mark)); +void set_ZP PROTO((double ZERO)); +int setrefcode PROTO((Image *image, int Nimage)); +void skip_measurements PROTO((Catalog *catalog, int pass)); +void sortA PROTO((double *X, int N)); +void sortB PROTO((double *X, double *Y, int N)); +void sortC PROTO((double *X, double *Y, double *F1, double *F2, int N)); +void sortD PROTO((double *X, double *Y, double *Z, int N)); +StatType statsImageM PROTO((Catalog *catalog)); +StatType statsImageN PROTO((Catalog *catalog)); +StatType statsImageX PROTO((Catalog *catalog)); +StatType statsImagedM PROTO((Catalog *catalog)); +StatType statsMosaicM PROTO((Catalog *catalog)); +StatType statsMosaicN PROTO((Catalog *catalog)); +StatType statsMosaicX PROTO((Catalog *catalog)); +StatType statsMosaicdM PROTO((Catalog *catalog)); +StatType statsStarN PROTO((Catalog *catalog, int Ncatalog)); +StatType statsStarS PROTO((Catalog *catalog, int Ncatalog)); +StatType statsStarX PROTO((Catalog *catalog, int Ncatalog)); +void wcatalog PROTO((Catalog *catalog)); +void wimages PROTO(()); +void write_coords PROTO((Header *header, Coords *coords)); + + +double **array_init (int Nx, int Ny); +void array_free (double **array, int Nx); +CoordFit *fit_init (int order); +void fit_free (CoordFit *fit); +void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt); +void fit_eval (CoordFit *fit); +double **poly2d_dx (double **poly, int Nx, int Ny); +double **poly2d_dy (double **poly, int Nx, int Ny); +double **poly2d_copy (double **poly, int Nx, int Ny); +double poly2d_eval (double **poly, int Nx, int Ny, double x, double y); +void fit_apply_coords (CoordFit *fit, Coords *coords); +int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo); +CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo); +void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords); +void initObjectData (Catalog *catalog, int Ncatalog); +int UpdateObjects (Catalog *catalog, int Ncatalog); +int UpdateSimple (Catalog *catalog, int Ncatalog); +int UpdateChips (Catalog *catalog, int Ncatalog); +int UpdateMosaic (Catalog *catalog, int Ncatalog); +int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon); +int ParFactor (double *pR, double *pD, double R, double D, time_t T); +int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts); +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts); + +Image *getMosaicForImage (int N); + +StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic); +StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic); + +double getMeanR (Measure *measure, Average *average, SecFilt *secfilt); +double getMeanD (Measure *measure, Average *average, SecFilt *secfilt); +int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt); +int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt); Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ConfigInit.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ConfigInit.c (revision 21919) @@ -0,0 +1,80 @@ +# include "relastro.h" + +void ConfigInit (int *argc, char **argv) { + + double ZERO_POINT; + char *config, *file; + char CatdirPhotcodeFile[256]; + char MasterPhotcodeFile[256]; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + GetConfig (config, "MAG_LIM", "%lf", 0, &MAG_LIM); + GetConfig (config, "SIGMA_LIM", "%lf", 0, &SIGMA_LIM); + GetConfig (config, "STAR_SCATTER", "%lf", 0, &STAR_SCATTER); + GetConfig (config, "IMAGE_SCATTER", "%lf", 0, &IMAGE_SCATTER); + GetConfig (config, "IMAGE_OFFSET", "%lf", 0, &IMAGE_OFFSET); + + GetConfig (config, "STAR_CHISQ", "%lf", 0, &STAR_CHISQ); + GetConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN); + GetConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW); + GetConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW); + GetConfig (config, "IMAGE_TOOFEW", "%d", 0, &IMAGE_TOOFEW); + GetConfig (config, "IMAGE_GOOD_FRACTION", "%lf", 0, &IMAGE_GOOD_FRACTION); + + GetConfig (config, "GSCFILE", "%s", 0, GSCFILE); + GetConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig(config, "CATMODE", "%s", 0, CATMODE); + ScanConfig(config, "CATFORMAT", "%s", 0, CATFORMAT); + ScanConfig(config, "PHOTCODE_FILE", "%s", 0, MasterPhotcodeFile); + + sprintf (ImageCat, "%s/Images.dat", CATDIR); + + if (!ScanConfig (config, "SKY_DEPTH", "%d", 0, &SKY_DEPTH)) { + SKY_DEPTH = 2; + } + if (!ScanConfig (config, "SKY_TABLE", "%s", 0, SKY_TABLE)) { + SKY_TABLE[0] = 0; + } + + // GetConfig (config, "ZERO_PT", "%lf", 0, &ZERO_POINT); + // GetConfig (config, "RELPHOT_GRID_X", "%d", 0, &RELPHOT_GRID_X); + // GetConfig (config, "RELPHOT_GRID_Y", "%d", 0, &RELPHOT_GRID_Y); + // GetConfig (config, "RELPHOT_GRID_BINNING", "%d", 0, &RELPHOT_GRID_BINNING); + // GetConfig (config, "CAMERA_CONFIG", "%s", 0, CameraConfig); + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + /* XXX this does not yet write out the master photcode table */ + sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR); + if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) { + fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile); + exit (1); + } + SetZeroPoint (ZERO_POINT); + + free (config); + free (file); + +} + +void GetConfig (char *config, char *field, char *format, int N, void *ptr) { + + char *status; + + status = ScanConfig (config, field, format, N, ptr); + if (status == NULL) { + fprintf (stderr, "error in config, cannot find %s\n", field); + exit (1); + } + return; +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitChip.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitChip.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitChip.c (revision 21919) @@ -0,0 +1,36 @@ +# include "relastro.h" + +void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { + + int i; + CoordFit *fit; + + fit = fit_init (coords[0].Npolyterms); + for (i = 0; i < Nmatch; i++) { + if (raw[i].mask) continue; + fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0); + } + fit_eval (fit); + fit_apply_coords (fit, coords); + fit_free (fit); + /* FitChip and FitSimple update the coords in different ways? maybe not... */ +} + +/* XXX the supporting fitting functions are not quite right + I probably need to apply the projection and then do the fit + to L,M. I also need to incorporate the transformations + determined for psastro */ + + +/* in the simple case, we only have three coord systems of interest: + R,D : the sky + P,Q : the tangent plane + X,Y : the chip + + R,D -> P,Q (projection) + P,Q -> X,Y (polynomial transformation) +*/ + +/* XXX I'm not using the errors at all : this could at least be done with the dMag values */ + + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitMosaic.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitMosaic.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitMosaic.c (revision 21919) @@ -0,0 +1,36 @@ +# include "relastro.h" + +void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { + + int i; + CoordFit *fit; + + fit = fit_init (coords[0].Npolyterms); + for (i = 0; i < Nmatch; i++) { + if (raw[i].mask) continue; + fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0); + } + fit_eval (fit); + fit_apply_coords (fit, coords); + fit_free (fit); + /* FitChip and FitSimple update the coords in different ways? maybe not... */ +} + +/* XXX the supporting fitting functions are not quite right + I probably need to apply the projection and then do the fit + to L,M. I also need to incorporate the transformations + determined for psastro */ + + +/* in the simple case, we only have three coord systems of interest: + R,D : the sky + P,Q : the tangent plane + X,Y : the chip + + R,D -> P,Q (projection) + P,Q -> X,Y (polynomial transformation) +*/ + +/* XXX I'm not using the errors at all : this could at least be done with the dMag values */ + + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPM.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPM.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPM.c (revision 21919) @@ -0,0 +1,85 @@ +# include "relastro.h" + +/* do we want an init function which does the alloc and a clear function to free? */ +int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) { + + int i; + + double **A, **B; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + double chisq, Xf, Yf; + + /* do I need to do this as 2 2x2 matrix equations? */ + A = array_init (4, 4); + B = array_init (4, 1); + + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + for (i = 0; i < Npts; i++) { + /* handle case where dX or dY = 0.0 */ + wx = 1.0 / SQ(dX[i]); + wy = 1.0 / SQ(dY[i]); + + Wx += wx; + Wy += wy; + + Tx += T[i]*wx; + Ty += T[i]*wy; + + Tx2 += SQ(T[i])*wx; + Ty2 += SQ(T[i])*wy; + + Xs += X[i]*wx; + Ys += Y[i]*wy; + + XT += X[i]*T[i]*wx; + YT += Y[i]*T[i]*wy; + } + + A[0][0] = Wx; + A[0][1] = Tx; + + A[1][0] = Tx; + A[1][1] = Tx2; + + A[2][2] = Wy; + A[2][3] = Ty; + + A[3][2] = Ty; + A[3][3] = Ty2; + + B[0][0] = Xs; + B[1][0] = XT; + B[2][0] = Ys; + B[3][0] = YT; + + dgaussj (A, 4, B, 1); + + fit[0].Ro = B[0][0]; + fit[0].uR = B[1][0]; + fit[0].Do = B[2][0]; + fit[0].uD = B[3][0]; + fit[0].p = 0.0; + + fit[0].dRo = sqrt(A[0][0]); + fit[0].duR = sqrt(A[1][1]); + fit[0].dDo = sqrt(A[2][2]); + fit[0].duD = sqrt(A[3][3]); + fit[0].dp = 0.0; + + array_free (A, 4); + array_free (B, 4); + + // add up the chi square for the fit + chisq = 0.0; + for (i = 0; i < Npts; i++) { + Xf = fit[0].Ro + fit[0].uR*T[i]; + Yf = fit[0].Do + fit[0].uD*T[i]; + chisq += SQ(X[i] - Xf) / SQ(dX[i]); + chisq += SQ(Y[i] - Yf) / SQ(dY[i]); + } + fit[0].Nfit = Npts; + + // the reduced chisq is divided by (Ndof = 2*Npts - 4) + fit[0].chisq = chisq / (2.0*Npts - 4.0); + return (TRUE); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPMandPar.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitPMandPar.c (revision 21919) @@ -0,0 +1,98 @@ +# include "relastro.h" + +/* do we want an init function which does the alloc and a clear function to free? */ +int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) { + + int i; + + double **A, **B; + double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; + double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2; + + A = array_init (5, 5); + B = array_init (5, 1); + + PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0; + Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; + for (i = 0; i < Npts; i++) { + /* handle case where dX or dY = 0.0 */ + wx = 1.0 / SQ(dX[i]); + wy = 1.0 / SQ(dY[i]); + + Wx += wx; + Wy += wy; + + Tx += T[i]*wx; + Ty += T[i]*wy; + + Tx2 += SQ(T[i])*wx; + Ty2 += SQ(T[i])*wy; + + PR += pR[i]*wx; + PD += pD[i]*wy; + + PRT += pR[i]*T[i]*wx; + PDT += pD[i]*T[i]*wy; + + PRX += pR[i]*X[i]*wx; + PDY += pD[i]*Y[i]*wy; + + PR2 += SQ(pR[i])*wx; + PD2 += SQ(pD[i])*wy; + + Xs += X[i]*wx; + Ys += Y[i]*wy; + + XT += X[i]*T[i]*wx; + YT += Y[i]*T[i]*wy; + } + + A[0][0] = Wx; + A[0][1] = Tx; + A[0][4] = PR; + + A[1][0] = Tx; + A[1][1] = Tx2; + A[1][4] = PRT; + + A[2][2] = Wy; + A[2][3] = Ty; + A[2][4] = PD; + + A[3][2] = Ty; + A[3][3] = Ty2; + A[3][4] = PDT; + + A[4][0] = PR; + A[4][1] = PRT; + A[4][2] = PD; + A[4][3] = PDT; + A[4][4] = PR2 + PD2; + + B[0][0] = Xs; + B[1][0] = XT; + B[2][0] = Ys; + B[3][0] = YT; + B[4][0] = PRX + PDY; + + dgaussj (A, 5, B, 1); + + fit[0].Ro = B[0][0]; + fit[0].uR = B[1][0]; + fit[0].Do = B[2][0]; + fit[0].uD = B[3][0]; + fit[0].p = B[4][0]; + + fit[0].dRo = sqrt(A[0][0]); + fit[0].duR = sqrt(A[1][1]); + fit[0].dDo = sqrt(A[2][2]); + fit[0].duD = sqrt(A[3][3]); + fit[0].dp = sqrt(A[4][4]); + + array_free (A, 5); + array_free (B, 5); + + /* get the chisq from the matrix values */ + + return (TRUE); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitSimple.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitSimple.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/FitSimple.c (revision 21919) @@ -0,0 +1,35 @@ +# include "relastro.h" + +void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords) { + + int i; + CoordFit *fit; + + fit = fit_init (coords[0].Npolyterms); + for (i = 0; i < Nmatch; i++) { + if (raw[i].mask) continue; + fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0); + } + fit_eval (fit); + fit_apply_coords (fit, coords); + fit_free (fit); +} + +/* XXX the supporting fitting functions are not quite right + I probably need to apply the projection and then do the fit + to L,M. I also need to incorporate the transformations + determined for psastro */ + + +/* in the simple case, we only have three coord systems of interest: + R,D : the sky + P,Q : the tangent plane + X,Y : the chip + + R,D -> P,Q (projection) + P,Q -> X,Y (polynomial transformation) +*/ + +/* XXX I'm not using the errors at all : this could at least be done with the dMag values */ + + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ImageOps.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ImageOps.c (revision 21919) @@ -0,0 +1,261 @@ +# include "relastro.h" + +static unsigned int *start; +static unsigned int *stop; +static int **bin; + +static int **clist; +static int **mlist; +static int *Nlist; +static int *NLIST; + +static Image *image; +static int Nimage; + +Image *getimages (int *N) { + *N = Nimage; + return (image); +} + +Image *getimage (int N) { + return (&image[N]); +} + +void initImages (Image *input, int N) { + + int i; + + image = input; + Nimage = N; + + ALLOCATE (start, unsigned, Nimage); + ALLOCATE (stop, unsigned, Nimage); + + for (i = 0; i < Nimage; i++) { + start[i] = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1); + stop[i] = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1); + } +} + +void initImageBins (Catalog *catalog, int Ncatalog) { + + int i, j; + + ALLOCATE (bin, int *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; + } + + ALLOCATE (Nlist, int, Nimage); + ALLOCATE (NLIST, int, Nimage); + ALLOCATE (clist, int *, Nimage); + ALLOCATE (mlist, int *, Nimage); + + for (i = 0; i < Nimage; i++) { + Nlist[i] = 0; + NLIST[i] = 100; + ALLOCATE (clist[i], int, NLIST[i]); + ALLOCATE (mlist[i], int, NLIST[i]); + } +} + +void freeImageBins (int Ncatalog) { + + int i; + + for (i = 0; i < Ncatalog; i++) { + free (bin[i]); + } + free (bin); + for (i = 0; i < Nimage; i++) { + free (clist[i]); + free (mlist[i]); + } + free (clist); + free (mlist); +} + +/* match measurements to images */ +void findImages (Catalog *catalog, int Ncatalog) { + + int i, j; + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Nmeasure; j++) { + matchImage (catalog, j, i); + } + } +} + +/* modify this function to use the measure->imageID field */ +void matchImage (Catalog *catalog, int meas, int cat) { + + int i; + Measure *measure; + + measure = &catalog[cat].measure[meas]; + for (i = 0; i < Nimage; i++) { + if (image[0].photcode == -1) continue; + if (measure[0].photcode != image[i].photcode) continue; + if (measure[0].t < start[i]) continue; + if (measure[0].t > stop[i]) continue; + + bin[cat][meas] = i; + + clist[i][Nlist[i]] = cat; + mlist[i][Nlist[i]] = meas; + Nlist[i] ++; + + if (Nlist[i] == NLIST[i]) { + NLIST[i] += 100; + REALLOCATE (clist[i], int, NLIST[i]); + REALLOCATE (mlist[i], int, NLIST[i]); + } + return; + } +} + +Coords *getCoords (int meas, int cat) { + + int i; + + i = bin[cat][meas]; + if (i == -1) return (NULL); + return (&image[i].coords); +} + +void plot_images () { + + int i, bin; + double *xlist, *Mlist, *dlist; + Graphdata graphdata; + + if (FREEZE_IMAGES) return; + + ALLOCATE (xlist, double, Nimage); + ALLOCATE (dlist, double, Nimage); + ALLOCATE (Mlist, double, Nimage); + + /**** dMcal vs airmass ****/ + for (i = 0; i < Nimage; i++) { + Mlist[i] = image[i].Mcal; + dlist[i] = image[i].dMcal; + xlist[i] = image[i].secz; + } + + plot_defaults (&graphdata); + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png"); + plot_defaults (&graphdata); + plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL); + +# define NBIN 200 + REALLOCATE (xlist, double, NBIN); + REALLOCATE (Mlist, double, NBIN); + + /**** dMcal histgram ****/ + for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Nimage; i++) { + bin = image[i].dMcal / 0.00025; + bin = MAX (0, MIN (NBIN - 1, bin)); + Mlist[bin] += 1.0; + } + + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png"); + + free (dlist); + free (xlist); + free (Mlist); +} + +/* XXX re-write to use image->nstars if that is possible */ +StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) { + + int i, m, c; + + Image *mosaic; + Coords *moscoords; + StarData *raw; + + ALLOCATE (raw, StarData, Nlist[im]); + + mosaic = NULL; + moscoords = NULL; + if (isMosaic) { + mosaic = getMosaicForImage (im); + moscoords = &mosaic[0].coords; + } + + for (i = 0; i < Nlist[im]; i++) { + m = mlist[im][i]; + c = clist[im][i]; + + /* apply the current image transformation or use the current value of R+dR, D+dD? */ + raw[i].X = catalog[c].measure[m].Xccd; + raw[i].Y = catalog[c].measure[m].Yccd; + + raw[i].Mag = catalog[c].measure[m].M; + raw[i].dMag = catalog[c].measure[m].dM; + + /* note that for a Simple image, L,M = P,Q */ + XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords); + if (isMosaic) { + XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords); + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords); + } else { + raw[i].P = raw[i].L; + raw[i].Q = raw[i].M; + LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords); + } + } + + *Nstars = Nlist[im]; + return (raw); +} + +StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) { + + int i, m, c, n; + + Image *mosaic; + Coords *moscoords; + StarData *ref; + + ALLOCATE (ref, StarData, Nlist[im]); + + mosaic = NULL; + moscoords = NULL; + if (isMosaic) { + mosaic = getMosaicForImage (im); + moscoords = &mosaic[0].coords; + } + + for (i = 0; i < Nlist[im]; i++) { + m = mlist[im][i]; + c = clist[im][i]; + n = catalog[c].measure[m].averef; + + /* apply the current image transformation or use the current value of R+dR, D+dD? */ + ref[i].R = catalog[c].average[n].R; + ref[i].D = catalog[c].average[n].D; + + /* note that for a Simple image, L,M = P,Q */ + RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords); + if (isMosaic) { + LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords); + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, moscoords); + } else { + ref[i].L = ref[i].P; + ref[i].M = ref[i].Q; + LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords); + } + } + + *Nstars = Nlist[im]; + return (ref); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/MosaicOps.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/MosaicOps.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/MosaicOps.c (revision 21919) @@ -0,0 +1,302 @@ +# include "relastro.h" + +Image *getimages (int *N); + +static int Nmosaic; +static Mosaic *mosaic; + +static int *Nimlist; +static int **imlist; /* mosaic -> image[] */ +static int **bin; /* catalog, measure -> mosaic */ + +static int **clist; /* mosaic -> catalog[] */ +static int **mlist; /* mosiac -> measure[] */ +static int *Nlist; +static int *NLIST; + +/* find mosaic frames (unique time periods & photcode name matches mosaic) */ +/* XXX what is a mosaic? do we need to use the 'DIS/WRP' info to track? */ +void initMosaics (Image *image, int Nimage) { + + int i, j, status, found, NMOSAIC, *NIMLIST; + unsigned int start, stop; + char *pname; + + if (!MOSAICNAME[0]) return; + + Nmosaic = 0; + NMOSAIC = 10; + ALLOCATE (mosaic, Mosaic, NMOSAIC); + + ALLOCATE (imlist, int *, NMOSAIC); + ALLOCATE (Nimlist, int, NMOSAIC); + ALLOCATE (NIMLIST, int, NMOSAIC); + + /* generate list of unique mosaics */ + for (i = 0; i < Nimage; i++) { + + /* select valid mosaic images by photcode */ + pname = GetPhotcodeNamebyCode (image[i].photcode); + status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME)); + if (status) continue; + + /* set image time range */ + start = image[i].tzero - MAX(0.05*image[i].trate*image[i].NY, 1); + stop = image[i].tzero + MAX(1.05*image[i].trate*image[i].NY, 1); + + /* find existing mosaic with this time range */ + found = FALSE; + for (j = 0; !found && (j < Nmosaic); j++) { + if (stop < mosaic[j].start) continue; + if (start > mosaic[j].stop) continue; + found = TRUE; + + /* add image to mosaic image list */ + imlist[j][Nimlist[j]] = i; + Nimlist[j] ++; + if (Nimlist[j] == NIMLIST[j]) { + NIMLIST[j] += 10; + REALLOCATE (imlist[j], int, NIMLIST[j]); + } + + } + if (found) continue; + + /* a new mosaic, define ranges */ + mosaic[Nmosaic].start = start; + mosaic[Nmosaic].stop = stop; + mosaic[Nmosaic].Mcal = 0.0; + mosaic[Nmosaic].dMcal = 0.0; + mosaic[Nmosaic].Xm = 0.0; + mosaic[Nmosaic].code = image[i].code; + mosaic[Nmosaic].secz = image[i].secz; + + /* add image to mosaic image list */ + NIMLIST[Nmosaic] = 10; + Nimlist[Nmosaic] = 1; + ALLOCATE (imlist[Nmosaic], int, NIMLIST[Nmosaic]); + imlist[Nmosaic][0] = i; + + Nmosaic ++; + if (Nmosaic == NMOSAIC) { + NMOSAIC += 10; + REALLOCATE (mosaic, Mosaic, NMOSAIC); + REALLOCATE (imlist, int *, NMOSAIC); + REALLOCATE (Nimlist, int, NMOSAIC); + REALLOCATE (NIMLIST, int, NMOSAIC); + } + } + return; +} + +void initMosaicBins (Catalog *catalog, int Ncatalog) { + + int i, j; + + /* measure -> mosaic */ + if (!MOSAICNAME[0]) return; + ALLOCATE (bin, int *, Ncatalog); + for (i = 0; i < Ncatalog; i++) { + ALLOCATE (bin[i], int, MAX (catalog[i].Nmeasure, 1)); + for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1; + } + + /* mosaic -> measure */ + ALLOCATE (Nlist, int, Nmosaic); + ALLOCATE (NLIST, int, Nmosaic); + ALLOCATE (clist, int *, Nmosaic); + ALLOCATE (mlist, int *, Nmosaic); + + for (i = 0; i < Nmosaic; i++) { + Nlist[i] = 0; + NLIST[i] = 100; + ALLOCATE (clist[i], int, NLIST[i]); + ALLOCATE (mlist[i], int, NLIST[i]); + } +} + +void freeMosaicBins (int Ncatalog) { + + int i; + + /* measure -> mosaic */ + if (!MOSAICNAME[0]) return; + for (i = 0; i < Ncatalog; i++) { + free (bin[i]); + } + free (bin); + + /* mosaic -> measure */ + for (i = 0; i < Nmosaic; i++) { + free (clist[i]); + free (mlist[i]); + } + free (Nlist); + free (NLIST); + free (clist); + free (mlist); +} + +int findMosaics (Catalog *catalog, int Ncatalog) { + + int i, j; + + if (!MOSAICNAME[0]) return (FALSE); + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Nmeasure; j++) { + if (TimeSelect) { + if (catalog[i].measure[j].t < TSTART) continue; + if (catalog[i].measure[j].t > TSTOP) continue; + } + matchMosaics (catalog, j, i); + } + } + return (TRUE); +} + +/* modify this function to use the measure->imageID field ? */ +void matchMosaics (Catalog *catalog, int meas, int cat) { + + int i; + + for (i = 0; i < Nmosaic; i++) { + if (catalog[cat].measure[meas].t < mosaic[i].start) continue; + if (catalog[cat].measure[meas].t > mosaic[i].stop) continue; + + bin[cat][meas] = i; + + clist[i][Nlist[i]] = cat; + mlist[i][Nlist[i]] = meas; + Nlist[i] ++; + + if (Nlist[i] == NLIST[i]) { + NLIST[i] += 100; + REALLOCATE (clist[i], int, NLIST[i]); + REALLOCATE (mlist[i], int, NLIST[i]); + } + return; + } + return; +} + +void plot_mosaic_fields (Catalog *catalog) { + + int i, j, m, c, N, ave, Nimage; + double *xlist, *ylist; + double Xmin, Xmax, Ymin, Ymax; + char string[64]; + Image *image; + Graphdata graphdata; + + if (!MOSAICNAME[0]) return; + + image = getimages (&Nimage); + + N = 0; + for (i = 0; i < Nmosaic; i++) + N = MAX (N, Nlist[i]); + + ALLOCATE (xlist, double, N); + ALLOCATE (ylist, double, N); + + for (i = 0; i < Nmosaic; i++) { + N = 0; + Xmin = Ymin = +360.0; + Xmax = Ymax = -360.0; + for (j = 0; j < Nlist[i]; j++) { + + m = mlist[i][j]; + c = clist[i][j]; + + if (catalog[c].measure[m].flags & (ID_MEAS_AREA | ID_MEAS_NOCAL)) continue; + + ave = catalog[c].measure[m].averef; + xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0; + ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.0; + N++; + } + + sprintf (string, "Mosaic %d", i); + plot_defaults (&graphdata); + plot_list (&graphdata, xlist, ylist, N, string, NULL); + } + + free (ylist); + free (xlist); +} + +void plot_mosaics () { + + int i, bin; + double *xlist, *Mlist, *dlist; + Graphdata graphdata; + + if (!MOSAICNAME[0]) return; + + ALLOCATE (xlist, double, Nmosaic); + ALLOCATE (dlist, double, Nmosaic); + ALLOCATE (Mlist, double, Nmosaic); + + for (i = 0; i < Nmosaic; i++) { + Mlist[i] = mosaic[i].Mcal; + dlist[i] = mosaic[i].dMcal; + xlist[i] = mosaic[i].secz; + } + + plot_defaults (&graphdata); + graphdata.xmin = 0.95; + graphdata.xmax = 2.50; + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, Mlist, Nmosaic, "airmass vs Mcal", "airmass.png"); + plot_defaults (&graphdata); + graphdata.size = 1.5; + graphdata.ptype = 7; + plot_list (&graphdata, Mlist, dlist, Nmosaic, "Mcal vs dMcal", "MdM.png"); + +# define NBIN 200 + REALLOCATE (xlist, double, NBIN); + REALLOCATE (Mlist, double, NBIN); + + /**** dMcal histgram ****/ + for (i = 0; i < NBIN; i++) xlist[i] = 0.00005*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Nmosaic; i++) { + bin = mosaic[i].dMcal / 0.00005; + bin = MAX (0, MIN (NBIN - 1, bin)); + Mlist[bin] += 1.0; + } + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png"); + + free (dlist); + free (xlist); + free (Mlist); +} + +Image *getMosaicForImage (int Nim) { + + Image *mosaic; + + mosaic = NULL; + return mosaic; +} + +int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage) { + + int i, Imax, Nmax; + + Imax = 0; + Nmax = Nimlist[0]; + for (i = 0; i < Nmosaic; i++) { + if (Nimlist[i] > Nmax) { + Imax = i; + Nmax = Nimlist[i]; + } + } + + *refmosaic = &mosaic[Imax]; + *Nimage = Nmax; + return (imlist[Imax]); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ParFactor.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ParFactor.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/ParFactor.c (revision 21919) @@ -0,0 +1,62 @@ +# include "relastro.h" + +# if (0) +/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ +/* ra and dec are returned as decimal hours and decimal degrees. */ +void lpsun (double jd, double *ra, double *dec) { + + double n, L, g, lambda,epsilon,alpha,delta,x,y,z; + + n = jd - J2000; + L = 280.460 + 0.9856474 * n; + g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN; + lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN; + epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN; + + // this is the conversion from ecliptic to celestial coords + x = cos(lambda); + y = cos(epsilon)*sin(lambda); + z = sin(epsilon)*sin(lambda); + + *ra = (atan_circ(x,y))*HRS_IN_RADIAN; + *dec = (asin(z))*DEG_IN_RADIAN; +} +# endif + +/* code borrowed from Skycalc : fix this stuff XXX */ +/* Low precision formulae for the sun, from Almanac p. C24 (1990) */ +int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) { + + double n, L, g; + +# define J2000 2451545. /* Julian date at standard epoch */ + + n = jd - J2000; + L = 280.460 + 0.9856474 * n; + g = (357.528 + 0.9856003 * n)*RAD_DEG; + *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees + *beta = 0.0; // approx latitude + *epsilon = (23.439 - 0.0000004 * n); // obliquity of ecliptic in degrees + return TRUE; +} + +/* given RA, DEC, Time, calculate the parallax factor */ +int ParFactor (double *pR, double *pD, double R, double D, time_t T) { + + double jd; + double L, B, E, e, s, r, d; + + /* given a time T in UNIX seconds, determine the solar longitude S */ + + jd = sec_to_jd (T); + sun_ecliptic (jd, &L, &B, &E); + + e = E * RAD_DEG; + s = L * RAD_DEG; + r = R * RAD_DEG; + d = D * RAD_DEG; + + *pR = +(cos(e)*sin(s)*cos(r) - cos(s)*sin(r)); + *pD = -(cos(e)*sin(s)*sin(r) + cos(s)*cos(r))*sin(d) + sin(e)*sin(s)*cos(d); + return TRUE; +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/SetSignals.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/SetSignals.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/SetSignals.c (revision 21919) @@ -0,0 +1,118 @@ +# include "relastro.h" + +static int Protect = FALSE; +static int Trapped = FALSE; + +void TrapSignal (int sig) { + fprintf (stderr, "trapped signal %d\n", sig); + if (sig == 11) { + fprintf (stderr, "seg fault\n"); + exit (1); + } + if (Protect) { + Trapped = TRUE; + fprintf (stderr, "blocking until protected sections are clear\n"); + return; + } + Shutdown ("halted by signal (trapped)"); +} + +void SetProtect (int mode) { + Protect = mode; + if (Trapped && !Protect) Shutdown ("halted by signal (protect)"); +} + +int SetSignals () { + + int i; + + /* disable almost all signal interrupts */ + for (i = 0; i < 36; i++) { + switch (i) { + /* can't redirect this signals */ + case SIGKILL: /* kill -9: cannot be caught or ignored */ + case SIGSTOP: /* SIGSTOP: cannot be caught or ignored */ + /* ignore these signals */ + case SIGCHLD: /* child halted: ignore */ + case SIGPWR: /* power failure - why ignore this? */ + case SIGWINCH: /* window resized */ + case SIGCONT: /* continue - maintain this action */ + case SIGTSTP: /* stop signal sent from tty - why ignore? */ + case SIGURG: /* socket signal, ignore this */ + break; + + default: + signal (i, TrapSignal); + } + } + return (TRUE); +} +/* + + Signal Value Action Comment + ------------------------------------------------------------------------- + SIGHUP 1 A Hangup detected on controlling terminal + or death of controlling process + SIGINT 2 A Interrupt from keyboard + SIGQUIT 3 A Quit from keyboard + SIGILL 4 A Illegal Instruction + SIGABRT 6 C Abort signal from abort(3) + SIGFPE 8 C Floating point exception + SIGKILL 9 AEF Kill signal + SIGSEGV 11 C Invalid memory reference + SIGPIPE 13 A Broken pipe: write to pipe with no readers + SIGALRM 14 A Timer signal from alarm(2) + SIGTERM 15 A Termination signal + SIGUSR1 30,10,16 A User-defined signal 1 + SIGUSR2 31,12,17 A User-defined signal 2 + SIGCHLD 20,17,18 B Child stopped or terminated + SIGCONT 19,18,25 Continue if stopped + SIGSTOP 17,19,23 DEF Stop process + SIGTSTP 18,20,24 D Stop typed at tty + SIGTTIN 21,21,26 D tty input for background process + SIGTTOU 22,22,27 D tty output for background process + + Next various other signals. + + Signal Value Action Comment + --------------------------------------------------------------------- + SIGTRAP 5 CG Trace/breakpoint trap + SIGIOT 6 CG IOT trap. A synonym for SIGABRT + SIGEMT 7,-,7 G + SIGBUS 10,7,10 AG Bus error + SIGSYS 12,-,12 G Bad argument to routine (SVID) + SIGSTKFLT -,16,- AG Stack fault on coprocessor + SIGURG 16,23,21 BG Urgent condition on socket (4.2 BSD) + SIGIO 23,29,22 AG I/O now possible (4.2 BSD) + SIGPOLL AG A synonym for SIGIO (System V) + SIGCLD -,-,18 G A synonym for SIGCHLD + SIGXCPU 24,24,30 AG CPU time limit exceeded (4.2 BSD) + SIGXFSZ 25,25,31 AG File size limit exceeded (4.2 BSD) + SIGVTALRM 26,26,28 AG Virtual alarm clock (4.2 BSD) + SIGPROF 27,27,29 AG Profile alarm clock + SIGPWR 29,30,19 AG Power failure (System V) + SIGINFO 29,-,- G A synonym for SIGPWR + SIGLOST -,-,- AG File lock lost + SIGWINCH 28,28,20 BG Window resize signal (4.3 BSD, Sun) + SIGUNUSED -,31,- AG Unused signal + (Here - denotes that a signal is absent; there where three values are given, the first one is usually valid for alpha and + sparc, the middle one for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a + sparc.) + + The letters in the "Action" column have the following meanings: + + A Default action is to terminate the process. + + B Default action is to ignore the signal. + + C Default action is to dump core. + + D Default action is to stop the process. + + E Signal cannot be caught. + + F Signal cannot be ignored. + + G Not a POSIX.1 conformant signal. + +*/ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/Shutdown.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/Shutdown.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/Shutdown.c (revision 21919) @@ -0,0 +1,30 @@ +# include "relastro.h" + +static FITS_DB *db; + +void set_db (FITS_DB *in) { + db = in; +} + +/* clean up open / locked ImageCat before shutting down */ +int Shutdown (char *format, ...) { + va_list argp; + char *formatplus; + + ALLOCATE (formatplus, char, strlen(format)); + strcpy (formatplus, format); + strcat (formatplus, "\n"); + + va_start (argp, format); + vfprintf (stderr, formatplus, argp); + free (formatplus); + va_end (argp); + + SetProtect (TRUE); + gfits_db_close (db); + fprintf (stderr, "ERROR: addstar halted\n"); + exit (1); +} + + +/* XXX this is probably not needed anymore : just protect the write statements */ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/StarOps.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/StarOps.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/StarOps.c (revision 21919) @@ -0,0 +1,518 @@ +# include "relastro.h" + +static int Nmax; +static double *list; +static double *dlist; + +void initMrel (Catalog *catalog, int Ncatalog) { + + int i, j; + + Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nm); + } + } + + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); +} + +float getMrel (Catalog *catalog, int meas, int cat) { + + int ave; + float value; + + ave = catalog[cat].measure[meas].averef; + if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG); + + value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M; + return (value); +} + +int setMrel (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, N; + float Msys, Mcal, Mmos, Mgrid; + StatType stats; + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + m = catalog[i].average[j].offset; + + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + if ((Mcal = getMcal (m, i)) == NO_MAG) continue; + if ((Mmos = getMmos (m, i)) == NO_MAG) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + N++; + } + if (N < STAR_TOOFEW) { /* too few measurements */ + catalog[i].average[j].code |= ID_STAR_FEW; + } else { + catalog[i].average[j].code &= ~ID_STAR_FEW; + } + + liststats (list, dlist, N, &stats); + + catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma; + catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; + } + } + return (TRUE); +} + +int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) { + + int i, j, k, m, N, Nmax; + float Msys, Mcal, Mmos, Mgrid; + double *list, *dlist; + StatType stats; + + /* Nmeasure is now different, need to reallocate */ + Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nm); + } + } + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* skip stars already calibrated */ + if (catalog[i].found[j]) continue; + + N = 0; + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + if ((Mcal = getMcal (m, i)) == NO_MAG) continue; + if ((Mmos = getMmos (m, i)) == NO_MAG) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + N++; + } + if (N < 1) continue; + + liststats (list, dlist, N, &stats); + if (mark) catalog[i].found[j] = TRUE; + + /* use sigma or error in dM for output? */ + catalog[i].secfilt[PhotNsec*j+PhotSec].M = stats.mean; + catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma); + catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG; + } + } + + free (list); + free (dlist); + return (TRUE); +} + +/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_NOCAL */ +int setMcalOutput (Catalog *catalog, int Ncatalog) { + + int i, j, k, m; + float Mcal, Mmos, Mgrid; + + MEAS_BAD = ID_MEAS_NOCAL; + IMAGE_BAD = ID_IMAGE_NOCAL; + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + if ((Mcal = getMcal (m, i)) == NO_MAG) continue; + if ((Mmos = getMmos (m, i)) == NO_MAG) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid; + } + } + } + return (TRUE); +} + +void clean_stars (Catalog *catalog, int Ncatalog) { + + int i, j, Ndel, Nave, Ntot, mark; + float dM, Xm; + double Chisq, MaxScatter, MaxChisq; + double *xlist, *slist, *dlist; + StatType stats; + + if (VERBOSE) fprintf (stderr, "marking poor stars\n"); + + /* find Xm median -> ChiSq lim must be > median */ + for (i = Ntot = 0; i < Ncatalog; i++) { + Ntot += catalog[i].Naverage; + } + ALLOCATE (xlist, double, Ntot); + ALLOCATE (slist, double, Ntot); + ALLOCATE (dlist, double, Ntot); + for (i = Ntot = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].average[j].code & STAR_BAD) continue; + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + if (Xm == -1) continue; + Chisq = pow (10.0, 0.01*Xm); + xlist[Ntot] = Chisq; + slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + dlist[Ntot] = 1; + Ntot ++; + } + } + + initstats ("MEAN"); + liststats (xlist, dlist, Ntot, &stats); + MaxChisq = MAX (STAR_CHISQ, 2*stats.median); + liststats (slist, dlist, Ntot, &stats); + MaxScatter = MAX (STAR_SCATTER, 2*stats.median); + fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq); + + Ndel = Nave = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + Chisq = pow (10.0, 0.01*Xm); + mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq); + if (mark) { + catalog[i].average[j].code |= ID_STAR_POOR; + Ndel ++; + } else { + catalog[i].average[j].code &= ~ID_STAR_POOR; + } + Nave ++; + } + } + fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave); + initstats (STATMODE); + free (xlist); + free (slist); + free (dlist); +} + +void clean_measures (Catalog *catalog, int Ncatalog, int final) { + + int i, j, k, m, N, Ndel, Nave, Nmax, image_bad, TOOFEW; + int *ilist; + double *tlist, *list, *dlist, Ns; + float Msys, Mcal, Mmos, Mgrid; + StatType stats; + + if (VERBOSE) fprintf (stderr, "marking poor measures\n"); + /* Nmeasure is now different, need to reallocate */ + Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nm); + } + } + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); + ALLOCATE (ilist, int, Nmax); + ALLOCATE (tlist, double, Nmax); + + /* it makes no sense to mark 3-sigma outliers with <5 measurements */ + TOOFEW = MAX (5, STAR_TOOFEW); + + Ns = 3; + Ndel = Nave = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* skip bad stars to prevent them from becoming good (on inner sample) */ + if (catalog[i].average[j].code & STAR_BAD) continue; + + /* on final processing, skip stars already measured */ + if (final && catalog[i].found[j]) continue; + + /* accumulate list of valid measurements */ + m = catalog[i].average[j].offset; + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + /* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */ + if ((Mcal = getMcal (m, i)) == NO_MAG) continue; + if ((Mmos = getMmos (m, i)) == NO_MAG) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + N++; + } + if (N < TOOFEW) continue; + + /* 3-sigma clip based on stats of inner 50% */ + initstats ("INNER_MEAN"); + liststats (list, dlist, N, &stats); + stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */ + for (k = m = 0; k < N; k++) { + if (fabs (list[k] - stats.median) < Ns*stats.sigma) { + list[m] = list[k]; + m++; + } + } + initstats ("MEAN"); + liststats (list, dlist, m, &stats); + stats.sigma = MAX (MIN_ERROR, stats.sigma); + + /* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */ + image_bad = IMAGE_BAD; + IMAGE_BAD = ID_IMAGE_NOCAL; + m = catalog[i].average[j].offset; + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + /* if (catalog[i].measure[m].flags & MEAS_BAD) continue; */ + if ((Mcal = getMcal (m, i)) == NO_MAG) continue; + if ((Mmos = getMmos (m, i)) == NO_MAG) continue; + if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue; + + Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + list[N] = Msys - Mcal - Mmos - Mgrid; + dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + ilist[N] = m; + N++; + Nave ++; + } + if (N < TOOFEW) continue; + + /* mark bad measures */ + for (k = 0; k < N; k++) { + if (fabs (list[k] - stats.median) > Ns*stats.sigma) { + catalog[i].measure[ilist[k]].flags |= ID_MEAS_POOR; + Ndel ++; + } + } + IMAGE_BAD = image_bad; + } + } + initstats (STATMODE); + fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); + free (ilist); + free (tlist); +} + +StatType statsStarN (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, n, N, Ntot; + double *list, *dlist; + StatType stats; + + Ntot = 0; + for (i = 0; i < Ncatalog; i++) { + Ntot += catalog[i].Naverage; + } + + ALLOCATE (list, double, Ntot); + ALLOCATE (dlist, double, Ntot); + + n = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + m = catalog[i].average[j].offset; + + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (getMcal (m, i) == NO_MAG) continue; + if (getMmos (m, i) == NO_MAG) continue; + if (getMgrid (m, i) == NO_MAG) continue; + N++; + } + + list[n] = N; + dlist[n] = 1; + n++; + } + } + + liststats (list, dlist, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsStarX (Catalog *catalog, int Ncatalog) { + + int i, j, n, Ntot, Xm; + double *list, *dlist; + StatType stats; + + Ntot = 0; + for (i = 0; i < Ncatalog; i++) { + Ntot += catalog[i].Naverage; + } + + ALLOCATE (list, double, Ntot); + ALLOCATE (dlist, double, Ntot); + + n = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + + Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + if (Xm == NO_MAG) continue; + list[n] = pow (10.0, 0.01*Xm); + dlist[n] = 1; + n++; + } + } + + liststats (list, dlist, n, &stats); + free (list); + free (dlist); + return (stats); +} + +StatType statsStarS (Catalog *catalog, int Ncatalog) { + + int i, j, n, Ntot; + double *list, *dlist; + float dM; + StatType stats; + + Ntot = 0; + for (i = 0; i < Ncatalog; i++) { + Ntot += catalog[i].Naverage; + } + + ALLOCATE (list, double, Ntot); + ALLOCATE (dlist, double, Ntot); + + n = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + + dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + list[n] = dM; + dlist[n] = 1; + n++; + } + } + + liststats (list, dlist, n, &stats); + free (list); + free (dlist); + return (stats); +} + +void plot_stars (Catalog *catalog, int Ncatalog) { + + int i, j, bin; + float dMrel; + double *xlist, *Mlist; + Graphdata graphdata; + +# define NBIN 200 + ALLOCATE (xlist, double, NBIN); + ALLOCATE (Mlist, double, NBIN); + + for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i; + bzero (Mlist, NBIN*sizeof(double)); + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].average[j].code & STAR_BAD) continue; + dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM; + bin = dMrel / 0.00025; + bin = MAX (0, MIN (NBIN-1, bin)); + Mlist[bin] += 1.0; + } + } + + plot_defaults (&graphdata); + graphdata.style = 1; + plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "dMhist.png"); + + free (xlist); + free (Mlist); +} + +void plot_chisq (Catalog *catalog, int Ncatalog) { + + int i, j, N, Ntotal, value; + double *xlist, *ylist; + Graphdata graphdata; + + Ntotal = 0; + for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage; + + ALLOCATE (xlist, double, Ntotal); + ALLOCATE (ylist, double, Ntotal); + + N = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + if (catalog[i].average[j].code & STAR_BAD) continue; + xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M; + value = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm; + if (value == NO_MAG) continue; + ylist[N] = 0.01*value; + N++; + } + } + + plot_defaults (&graphdata); + graphdata.ymin = -3.0; + plot_list (&graphdata, xlist, ylist, N, "chisq", NULL); + free (xlist); + free (ylist); +} + +void plot_star_coords (Catalog *catalog, int Ncatalog) { + + int i, j, N; + double *xlist, *ylist; + double Xmin, Ymin, Xmax, Ymax; + Graphdata graphdata; + + N = 0; + for (i = 0; i < Ncatalog; i++) { + N += catalog[i].Naverage; + } + ALLOCATE (xlist, double, N); + ALLOCATE (ylist, double, N); + + N = 0; + Xmin = Ymin = +360.0; + Xmax = Ymax = -360.0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + xlist[N] = catalog[i].average[j].R; + ylist[N] = catalog[i].average[j].D; + N++; + } + } + plot_defaults (&graphdata); + plot_list (&graphdata, xlist, ylist, N, "coords", NULL); + + free (xlist); + free (ylist); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateChips.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateChips.c (revision 21919) @@ -0,0 +1,36 @@ +# include "relastro.h" + +int UpdateChips (Catalog *catalog, int Ncatalog) { + + /* we can measure new image parameters for each non-mosaic chip independently */ + int i, Nimage, Nstars; + Image *image; + StarData *raw, *ref; + + image = getimages (&Nimage); + + for (i = 0; i < Nimage; i++) { + + /* skip all except WRP images */ + if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + + /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */ + + /* convert measure coordinates to raw entries */ + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE); + + /* convert average coordinates to ref entries */ + ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE); + + /* XXX grab this code from mosastro */ + /* XXX check that it matches psastro */ + /* XXX this may need a 'mosaic' and 'non-mosaic' version */ + FitChip (raw, ref, Nstars, &image[i].coords); + + free (raw); + free (ref); + } + + return (TRUE); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateImages.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateImages.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateImages.c (revision 21919) @@ -0,0 +1,72 @@ +# include "relastro.h" + +static int Nmax; +static double *listR, *listD; +static double *dlistR, *dlistD; + +void initObjectData (Catalog *catalog, int Ncatalog) { + + int i, j; + + Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nm); + } + } + + ALLOCATE (list, double, MAX (1, Nmax)); + ALLOCATE (dlist, double, MAX (1, Nmax)); +} + +int UpdateObjects (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, N; + float chisq; + StatType statsR, statsD; + + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value of R,D for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + m = catalog[i].average[j].offset; + + N = 0; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + + listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + + /* the astrometric errors are not being carried yet (but should be!) */ + /* we use the photometric mag error as a weighting term */ + dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + N++; + } + if (N < STAR_TOOFEW) { /* too few measurements */ + catalog[i].average[j].code |= ID_STAR_FEW; + } else { + catalog[i].average[j].code &= ~ID_STAR_FEW; + } + + /* in here, we should be fitting the parallax and proper-motion components */ + + liststats (listR, dlistR, N, &statsR); + liststats (listR, dlistR, N, &statsD); + + catalog[i].average[j].R = statsR.mean; + catalog[i].average[j].dR = statsR.sigma; + + catalog[i].average[j].D = statsD.mean; + catalog[i].average[j].dD = statsD.sigma; + + chisq = 0.5*(statsR.chisq + statsD.chisq); + catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG; + + } + } + return (TRUE); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateMosaic.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateMosaic.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateMosaic.c (revision 21919) @@ -0,0 +1,36 @@ +# include "relastro.h" + +int UpdateMosaic (Catalog *catalog, int Ncatalog) { + + /* we can measure new image parameters for each mosaic independently */ + int i, Nimage, Nstars; + Image *image; + StarData *raw, *ref; + + image = getimages (&Nimage); + + for (i = 0; i < Nimage; i++) { + + /* skip all except DIS images */ + if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + + /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */ + + /* convert measure coordinates to raw entries */ + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE); + + /* convert average coordinates to ref entries */ + ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE); + + /* XXX grab this code from mosastro */ + /* XXX check that it matches psastro */ + /* XXX this may need a 'mosaic' and 'non-mosaic' version */ + FitMosaic (raw, ref, Nstars, &image[i].coords); + + free (raw); + free (ref); + } + + return (TRUE); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateObjects.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateObjects.c (revision 21919) @@ -0,0 +1,231 @@ +# include "relastro.h" + +static int Nmax; +static double *X, *dX; +static double *Y, *dY; +static double *R, *dR; +static double *D, *dD; +static double *pX; +static double *pY; +static double *T; +static double *dT; + +void initObjectData (Catalog *catalog, int Ncatalog) { + + int i, j; + + Nmax = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Nmax = MAX (Nmax, catalog[i].average[j].Nm); + } + } + + ALLOCATE (R, double, MAX (1, Nmax)); + ALLOCATE (D, double, MAX (1, Nmax)); + ALLOCATE (T, double, MAX (1, Nmax)); + ALLOCATE (X, double, MAX (1, Nmax)); + ALLOCATE (Y, double, MAX (1, Nmax)); + + ALLOCATE (dR, double, MAX (1, Nmax)); + ALLOCATE (dD, double, MAX (1, Nmax)); + ALLOCATE (dT, double, MAX (1, Nmax)); + ALLOCATE (dX, double, MAX (1, Nmax)); + ALLOCATE (dY, double, MAX (1, Nmax)); + + ALLOCATE (pX, double, MAX (1, Nmax)); + ALLOCATE (pY, double, MAX (1, Nmax)); +} + +int UpdateObjects (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, N; + StatType statsR, statsD; + Coords coords; + PMFit fit; + time_t To; + int mode, Nave, Npm, Npar, Nskip; + double Tmin, Tmax; + + initObjectData (catalog, Ncatalog); + + /* project coordinates to a plane centered on the object with units of arcsec */ + coords.crval1 = 0; + coords.crval2 = 0; + coords.crpix1 = 0; + coords.crpix2 = 0; + coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0; + coords.pc1_1 = coords.pc2_2 = 1.0; + coords.pc1_2 = coords.pc2_1 = 0.0; + coords.Npolyterms = 1; + strcpy (coords.ctype, "RA---SIN"); + + // use J2000 as a reference time + To = date_to_sec ("2000/01/01"); + Nave = Npar = Npm = 0; + + for (i = 0; i < Ncatalog; i++) { + + if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure); + + Nskip = 0; + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value of R,D for a single star */ + if (catalog[i].average[j].code & STAR_BAD) { + Nskip ++; + continue; + } + + N = 0; + m = catalog[i].average[j].offset; + + Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25); + mode = FIT_MODE; + + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + + R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years + + Tmin = MIN(Tmin, T[N]); + Tmax = MAX(Tmax, T[N]); + + /* the astrometric errors are not being carried yet (but should be!) */ + /* we use the photometric mag error as a weighting term */ + dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR); + dT[N] = catalog[i].measure[m].dt; + + N++; + } + + if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE; + if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE; + + // XXX This criterion needs to be better considered: adjust to match Ndof + if (N < POS_TOOFEW) { /* too few measurements */ + catalog[i].average[j].code |= ID_STAR_FEW; + continue; + } else { + catalog[i].average[j].code &= ~ID_STAR_FEW; + } + + /* we need to do the fit in a locally linear space; choose a ref coordinate */ + coords.crval1 = R[0]; + coords.crval2 = D[0]; + + /* project all of the R,D coordinates to a plane centered on this coordinate */ + for (k = 0; k < N; k++) { + RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords); + dX[k] = dR[k]; + dY[k] = dD[k]; + // fprintf (stderr, "%d %f %f %f %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]); + } + + /* fit the model components as needed */ + switch (mode) { + case FIT_AVERAGE: + liststats (R, dR, N, &statsR); + liststats (D, dD, N, &statsD); + + fit.Ro = statsR.mean; + fit.dRo = 3600.0*statsR.sigma; + + fit.Do = statsD.mean; + fit.dDo = 3600.0*statsD.sigma; + + fit.chisq = 0.5*(statsR.chisq + statsD.chisq); + fit.Nfit = N; + Nave ++; + break; + + case FIT_PM_ONLY: + FitPM (&fit, X, dX, Y, dY, T, N); + // fprintf (stderr, "fitted: %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); + // project Ro, Do back to RA,DEC + XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords); + // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p); + // continue; + Npm ++; + break; + + case FIT_PM_AND_PAR: + for (k = 0; k < N; k++) { + ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]); + } + FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N); + Npar ++; + break; + default: + fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__); + exit (2); + } + + if (0 && (j < 100)) { + fprintf (stderr, "%f %f -> %f %f (%f,%f)\n", + catalog[i].average[j].R, + catalog[i].average[j].D, + fit.Ro, fit.Do, + 3600*(catalog[i].average[j].R - fit.Ro), + 3600*(catalog[i].average[j].D - fit.Do)); + } + + // the measure fields must be updated before the average fields + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]); + } + + catalog[i].average[j].R = fit.Ro; // RA in degrees + catalog[i].average[j].D = fit.Do; // DEC in degrees + catalog[i].average[j].dR = fit.dRo; // RA scatter in arcsec + catalog[i].average[j].dD = fit.dDo; // DEC scatter in arcsec + + catalog[i].average[j].uR = fit.uR; // RA proper motion in arcsec/year + catalog[i].average[j].uD = fit.uD; // DEC proper motion in arcsec/year + catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year + catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year + + catalog[i].average[j].P = fit.p; // parallax in arcsec + catalog[i].average[j].dP = fit.dp; // parallax error in arcsec + + catalog[i].average[j].Xp = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG; + } + + if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip); + } + + if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar); + return (TRUE); +} + +/* fitting proper-motion and parallax: + +given a source at position r,d, at a time t, we need to calculate a vector (pr,pd) + +let x,y be the coordinate in the linearized frame with y parallel to DEC lines + +L,B are the ecliptic longitude and latitude of the object, +dL and dB are the offsets in the L and B directions + +dL = sin(t - topp) +dB = cos(t - topp)*sin(B) + +these need to be rotated to the R,D frame to yield pR,pD. Then, the equation of motion +for the source in the x,y frame is: + +x = Ro + uR * (t - to) + p * pR +y = Do + uD * (t - to) + p * pD + +the unknowns in these equations are Ro, uR, Do, uD, and p + +XXX think through the concepts for the pole a bit better. all objects near the pole +move the same way with the same phase. choose a projection center and define dL,dB relative +to that center point coordinate system? + +*/ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateSimple.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateSimple.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/UpdateSimple.c (revision 21919) @@ -0,0 +1,37 @@ +# include "relastro.h" + +int UpdateSimple (Catalog *catalog, int Ncatalog) { + + /* we can measure new image parameters for each non-mosaic chip independently */ + int i, Nimage, Nstars; + Image *image; + StarData *raw, *ref; + + image = getimages (&Nimage); + + for (i = 0; i < Nimage; i++) { + + /* skip WRP and DIS images */ + if (!strcmp(&image[i].coords.ctype[4], "-WRP")) continue; + if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue; + + /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */ + + /* convert measure coordinates to raw entries */ + raw = getImageRaw (catalog, Ncatalog, i, &Nstars, FALSE); + + /* convert average coordinates to ref entries */ + ref = getImageRef (catalog, Ncatalog, i, &Nstars, FALSE); + + /* XXX grab this code from mosastro */ + /* XXX check that it matches psastro */ + /* XXX this may need a 'mosaic' and 'non-mosaic' version */ + FitSimple (raw, ref, Nstars, &image[i].coords); + + free (raw); + free (ref); + } + + return (TRUE); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/args.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/args.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/args.c (revision 21919) @@ -0,0 +1,213 @@ +# include "relastro.h" +void usage (void); + +int args (int argc, char **argv) { + + int N; + double trange; + + /* possible operations */ + FIT_TARGET = TARGET_NONE; + FIT_MODE = FIT_AVERAGE; + if ((N = get_argument (argc, argv, "-update-objects"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_OBJECTS; + + // check for object fitting modes (not valid for images) + if ((N = get_argument (argc, argv, "-pm"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_ONLY; + } + if ((N = get_argument (argc, argv, "-pmpar"))) { + remove_argument (N, &argc, argv); + FIT_MODE = FIT_PM_AND_PAR; + } + } + if ((N = get_argument (argc, argv, "-update-simple"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_SIMPLE; + } + if ((N = get_argument (argc, argv, "-update-chips"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_CHIPS; + } + if ((N = get_argument (argc, argv, "-update-mosaics"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = TARGET_MOSAICS; + } + if (FIT_TARGET == TARGET_NONE) usage(); + + /* specify portion of the sky : allow default of all sky? */ + UserPatch.Rmin = 0; + UserPatch.Rmax = 360; + UserPatch.Dmin = -90; + UserPatch.Dmax = +90; + if ((N = get_argument (argc, argv, "-region"))) { + remove_argument (N, &argc, argv); + UserPatch.Rmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Rmax = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmin = atof (argv[N]); + remove_argument (N, &argc, argv); + UserPatch.Dmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } else { + usage (); + } + + /* define time */ + TimeSelect = FALSE; + if ((N = get_argument (argc, argv, "-time"))) { + remove_argument (N, &argc, argv); + if (!str_to_time (argv[N], &TSTART)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + remove_argument (N, &argc, argv); + if (!str_to_dtime (argv[N], &trange)) { + if (!str_to_time (argv[N], &TSTOP)) { + fprintf (stderr, "ERROR: syntax error\n"); + return (FALSE); + } + } else { + if (trange < 0) { + trange = fabs (trange); + TSTOP = TSTART; + TSTART -= trange; + } else { + TSTOP = TSTART + trange; + } + } + remove_argument (N, &argc, argv); + TimeSelect = TRUE; + } + + PHOTCODE_LIST = NULL; + if ((N = get_argument (argc, argv, "-photcode"))) { + remove_argument (N, &argc, argv); + PHOTCODE_LIST = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + + VERBOSE = FALSE; + if ((N = get_argument (argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTSTUFF = FALSE; + if ((N = get_argument (argc, argv, "-plot"))) { + PLOTSTUFF = TRUE; + remove_argument (N, &argc, argv); + } + + PLOTDELAY = 500000; + if ((N = get_argument (argc, argv, "-plotdelay"))) { + remove_argument (N, &argc, argv); + PLOTDELAY = 1e6*atof(argv[N]); + remove_argument (N, &argc, argv); + } + + strcpy (STATMODE, "CHI_INNER_WTMEAN"); + if ((N = get_argument (argc, argv, "-statmode"))) { + remove_argument (N, &argc, argv); + strcpy (STATMODE, argv[N]); + remove_argument (N, &argc, argv); + } + + RESET = FALSE; + if ((N = get_argument (argc, argv, "-reset"))) { + remove_argument (N, &argc, argv); + RESET = TRUE; + } + + UPDATE = FALSE; + if ((N = get_argument (argc, argv, "-update"))) { + remove_argument (N, &argc, argv); + UPDATE = TRUE; + } + + SHOW_PARAMS = FALSE; + if ((N = get_argument (argc, argv, "-params"))) { + remove_argument (N, &argc, argv); + SHOW_PARAMS = TRUE; + } + + /* XXX update these for relevant plots */ + PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0; + if ((N = get_argument (argc, argv, "-plrange"))) { + remove_argument (N, &argc, argv); + PlotMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmin = atof (argv[N]); + remove_argument (N, &argc, argv); + PlotdMmax = atof (argv[N]); + remove_argument (N, &argc, argv); + } + + /* XXX update this */ + MIN_ERROR = 0.001; + if ((N = get_argument (argc, argv, "-minerror"))) { + remove_argument (N, &argc, argv); + MIN_ERROR = atof (argv[N]); + remove_argument (N, &argc, argv); + /* require MIN_ERROR > 0 */ + } + + AreaSelect = FALSE; + if ((N = get_argument (argc, argv, "-area"))) { + remove_argument (N, &argc, argv); + AreaXmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaXmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmin = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaYmax = atof (argv[N]); + remove_argument (N, &argc, argv); + AreaSelect = TRUE; + } + + ImagSelect = FALSE; + if ((N = get_argument (argc, argv, "-instmag"))) { + remove_argument (N, &argc, argv); + ImagMin = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagMax = atof (argv[N]); + remove_argument (N, &argc, argv); + ImagSelect = TRUE; + } + + /* XXX drop this? */ + DophotSelect = FALSE; + if ((N = get_argument (argc, argv, "-dophot"))) { + remove_argument (N, &argc, argv); + DophotValue = atof (argv[N]); + remove_argument (N, &argc, argv); + DophotSelect = TRUE; + } + + if (argc != 1) usage (); + return TRUE; +} + +void usage () { + fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n"); + fprintf (stderr, " working options: \n"); + fprintf (stderr, " -time (start) (stop)\n"); + fprintf (stderr, " -v\n"); + fprintf (stderr, " -plot\n"); + fprintf (stderr, " -plotdelay (seconds)\n"); + fprintf (stderr, " -update\n"); + fprintf (stderr, " -params\n"); + fprintf (stderr, " -reset\n"); + fprintf (stderr, " -area Xmin Xmax Ymin Ymax\n"); + fprintf (stderr, " -instmag min max\n\n"); + fprintf (stderr, " planned options: \n"); + fprintf (stderr, " -photcode code[,code,code,..]\n"); + fprintf (stderr, " \n"); + exit (2); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/bcatalog.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/bcatalog.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/bcatalog.c (revision 21919) @@ -0,0 +1,112 @@ +# include "relastro.h" + +int bcatalog (Catalog *subcatalog, Catalog *catalog) { + + int i, j, k, offset, found; + int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm; + float mag; + + /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */ + NAVERAGE = 50; + NMEASURE = 1000; + ALLOCATE (subcatalog[0].average, Average, NAVERAGE); + ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec); + ALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + Nmeasure = Naverage = 0; + + /* exclude stars not in range or with too few measurements */ + for (i = 0; i < catalog[0].Naverage; i++) { + if (catalog[0].average[i].Nm < 2) continue; + + /* start with all stars good */ + subcatalog[0].average[Naverage] = catalog[0].average[i]; + subcatalog[0].average[Naverage].offset = Nmeasure; + for (j = 0; j < PhotNsec; j++) { + subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j]; + } + + if (RESET) { + // XXX reset the ra,dec coords? + subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW; + subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR; + } + + Nm = 0; + for (j = 0; j < catalog[0].average[i].Nm; j++) { + + offset = catalog[0].average[i].offset + j; + + /* select measurements by photcode, if specified */ + if (Nphotcodes > 0) { + found = FALSE; + for (k = 0; (k < Nphotcodes) && !found; k++) { + if (photcodes[k][0].code == catalog[0].measure[offset].photcode) found = TRUE; + } + if (!found) continue; + } + + /* select measurements by time */ + if (TimeSelect) { + if (catalog[0].measure[offset].t < TSTART) continue; + if (catalog[0].measure[offset].t > TSTOP) continue; + } + + /* select measurements by quality */ + if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue; + + /* select measurements by mag limit */ + mag = PhotCat (&catalog[0].measure[offset]); + if (mag > MAG_LIM) continue; + + /* select measurements by measurement error */ + if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue; + + /* select measurements by mag limit */ + if (ImagSelect) { + mag = PhotInst (&catalog[0].measure[offset]); + if (mag < ImagMin) continue; + if (mag > ImagMax) continue; + } + + subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_SKIP; + subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; + subcatalog[0].measure[Nmeasure].averef = Naverage; + if (RESET) { + subcatalog[0].measure[Nmeasure].Mcal = 0; + subcatalog[0].measure[Nmeasure].flags &= 0xff00; + subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_POOR; + subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_AREA; + subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_NOCAL; + } + Nmeasure ++; + Nm ++; + if (Nmeasure == NMEASURE) { + NMEASURE += 1000; + REALLOCATE (subcatalog[0].measure, Measure, NMEASURE); + } + } + if (Nm < 2) { /* enough measurements in band? */ + Nmeasure -= Nm; + continue; + } + subcatalog[0].average[Naverage].Nm = Nm; + Naverage ++; + if (Naverage == NAVERAGE) { + NAVERAGE += 50; + REALLOCATE (subcatalog[0].average, Average, NAVERAGE); + REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec); + } + } + REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); + REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1)); + REALLOCATE (subcatalog[0].secfilt, SecFilt, PhotNsec*MAX (Naverage, 1)); + subcatalog[0].Naverage = Naverage; + subcatalog[0].Nmeasure = Nmeasure; + subcatalog[0].Nsecfilt = catalog[0].Nsecfilt; + + if (VERBOSE) { + fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, + subcatalog[0].Naverage, subcatalog[0].Nmeasure); + } + return (TRUE); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/dvo_astrom_ops.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/dvo_astrom_ops.c (revision 21919) @@ -0,0 +1,91 @@ +# include "relastro.h" + +double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) { + + double ra; + + /* the measure carries the instantaneous mean position at the epoch t */ + ra = average[0].R - measure[0].dR / 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (ra); +} + +double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) { + + double dec; + + /* the measure carries the instantaneous mean position at the epoch t */ + dec = average[0].D - measure[0].dD / 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (dec); +} + +int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) { + + /* math to get from new fitted position to new measure offset + ra_obs = average[0].R - measure[0].dR / 3600.0; + measure[0].dR = (ra_fit - ra_obs) * 3600.0; + measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0 + measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR; + */ + + /* the measure carries the instantaneous mean position at the epoch t */ + measure[0].dR += (ra_fit - average[0].R) * 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (TRUE); +} + +int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) { + + /* math to get from new fitted position to new measure offset + dec_obs = average[0].D - measure[0].dD / 3600.0; + measure[0].dD = (dec_fit - dec_obs) * 3600.0; + measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0 + measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD; + */ + + /* the measure carries the instantaneous mean position at the epoch t */ + measure[0].dD += (dec_fit - average[0].D) * 3600.0; + + /* possible corrections to mean ra: + + - proper-motion and parallax + - abberation + - precession and nutation, etc + - refraction + - DCR + + */ + + return (TRUE); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/fitpoly.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/fitpoly.c (revision 21919) @@ -0,0 +1,295 @@ +# include "relastro.h" + +/* these functions support simultaneous 2D fits to + x2 = \sum a_i,j x1^i y1^j + y2 = \sum b_i,j x1^i y1^j + + the order of the fit (largest coefficient) is fixed to a single + value for both x1,y1 terms and for x2,y2 fits + + the code is currently confusing because we limit to i+j <= order. + this could be cleaner if we used masks and allowed i <= order, j <= order +*/ + +double **array_init (int Nx, int Ny) { + + int i; + double **array; + + ALLOCATE (array, double *, Nx); + for (i = 0; i < Nx; i++) { + ALLOCATE (array[i], double, Ny); + memset (array[i], 0, Ny*sizeof(double)); + } + return (array); +} + +void array_free (double **array, int Nx) { + + int i; + + for (i = 0; i < Nx; i++) { + free (array[i]); + } + free (array); +} + +// XXX define a fit structure and drop the file static variables? +CoordFit *fit_init (int order) { + + CoordFit *fit; + + ALLOCATE (fit, CoordFit, 1); + + fit[0].Npts = 0; + fit[0].Norder = order; + fit[0].Nterms = order + 1; + fit[0].Nsums = 2*order + 1; + fit[0].Nelems = SQ(order + 1); + + /* summing arrays for fit solution */ + + // xsum[i][j] holds \sum (x2 wt x1^i y1^j) + // ysum[i][j] holds \sum (y2 wt x1^i y1^j) + fit[0].xsum = array_init (fit[0].Nterms, fit[0].Nterms); + fit[0].ysum = array_init (fit[0].Nterms, fit[0].Nterms); + + // xfit[i][j] holds x2 coeff for x1^i y1^j + // yfit[i][j] holds y2 coeff for x1^i y1^j + fit[0].xfit = array_init (fit[0].Nterms, fit[0].Nterms); + fit[0].yfit = array_init (fit[0].Nterms, fit[0].Nterms); + + // sum[i][j] holds \sum (wt x1^i y1^j) + fit[0].sum = array_init (fit[0].Nsums, fit[0].Nsums); + + return (fit); +} + +void fit_free (CoordFit *fit) { + + array_free (fit[0].xfit, fit[0].Nterms); + array_free (fit[0].yfit, fit[0].Nterms); + + array_free (fit[0].sum, fit[0].Nsums); + array_free (fit[0].xsum, fit[0].Nterms); + array_free (fit[0].ysum, fit[0].Nterms); + + free (fit); +} + +// XXX use implicit masks as below or explicit masks (with function to set?) +// XXX eg, add a global mask to this file and +void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) { + + int ix, iy; + double xterm, yterm, term; + + xterm = 1; + for (ix = 0; ix < fit[0].Nsums; ix++) { + yterm = 1; + for (iy = 0; iy < fit[0].Nsums; iy++) { + term = xterm*yterm*wt; + fit[0].sum[ix][iy] += term; + if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) { + fit[0].xsum[ix][iy] += x2*term; + fit[0].ysum[ix][iy] += y2*term; + } + yterm *= y1; + } + xterm *= x1; + } + fit[0].Npts ++; +} + +/* convert the xsum,ysum,sum terms into vector,matrix and solve */ +void fit_eval (CoordFit *fit) { + + int i, j, ix, iy, jx, jy; + double **matrix, **vector; + + if (fit[0].Npts == 0) { + fprintf (stderr, "warning: no valid pts\n"); + } + + // matrix, vector hold the final linear system + matrix = array_init (fit[0].Nelems, fit[0].Nelems); + vector = array_init (fit[0].Nelems, 2); + + /* remap the xsum,ysum terms into the vector */ + for (i = 0; i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + vector[i][0] = fit[0].xsum[ix][iy]; + vector[i][1] = fit[0].ysum[ix][iy]; + + for (j = 0; j < fit[0].Nelems; j++) { + jx = j % fit[0].Nterms; + jy = j / fit[0].Nterms; + matrix[i][j] = fit[0].sum[ix+jx][iy+jy]; + } + } + + dgaussj (matrix, fit[0].Nelems, vector, 2); + + for (i = 0; i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + fprintf (stderr, "x2 : x^%dy^%d: %10.4g y2 : x^%dy^%d: %10.4g \n", + ix, iy, vector[i][0], ix, iy, vector[i][1]); + } + + /* remap the vector terms into xfit,yfit */ + for (i = 0; i < fit[0].Nelems; i++) { + ix = i % fit[0].Nterms; + iy = i / fit[0].Nterms; + fit[0].xfit[ix][iy] = vector[i][0]; + fit[0].yfit[ix][iy] = vector[i][1]; + } + + array_free (matrix, fit[0].Nelems); + array_free (vector, fit[0].Nelems); +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_dx (double **poly, int Nx, int Ny) { + + int i, j, Nxout, Nyout; + double **dpoly; + + Nxout = Nx - 1; + Nyout = Ny; + + // poly[i][j] holds coeff for x1^i y1^j + dpoly = array_init (Nxout, Nyout); + + for (i = 0; i < Nxout; i++) { + for (j = 0; j < Nyout; j++) { + dpoly[i][j] = poly[i+1][j] * (i+1); + } + } + return dpoly; +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_dy (double **poly, int Nx, int Ny) { + + int i, j, Nxout, Nyout; + double **dpoly; + + Nxout = Nx; + Nyout = Ny - 1; + + // poly[i][j] holds coeff for x1^i y1^j + dpoly = array_init (Nxout, Nyout); + + for (i = 0; i < Nxout; i++) { + for (j = 0; j < Nyout; j++) { + dpoly[i][j] = poly[i][j+1] * (j+1); + } + } + return dpoly; +} + +// Nx, Ny is the number of terms in x and in y +double **poly2d_copy (double **poly, int Nx, int Ny) { + + int i, j; + double **out; + + // poly[i][j] holds coeff for x1^i y1^j + out = array_init (Nx, Ny); + + for (i = 0; i < Nx; i++) { + for (j = 0; j < Ny; j++) { + out[i][j] = poly[i][j+1]; + } + } + return out; +} + +double poly2d_eval (double **poly, int Nx, int Ny, double x, double y) { + + int i, j; + double xsum, ysum, sum; + + sum = 0; + xsum = ysum = 1.0; + + for (i = 0; i < Nx; i++) { + ysum = xsum; + for (j = 0; j < Ny; j++) { + sum += poly[i][j] * ysum; + ysum *= y; + } + xsum *= x; + } + return (sum); +} + +/* linear portion of fit : NORDER is 1 */ +/* this should only apply to the polynomial, not the projection terms */ +/* compare with psastro supporting code */ +void fit_apply_coords (CoordFit *fit, Coords *coords) { + + double c11, c12; + double c21, c22; + double Xo, Yo, R; + CoordFit *modfit; + + /* I have L,M = fit(X,Y). set corresponding terms for coords */ + + // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ... + // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ... + + CoordsGetCenter (fit, 0.001, &Xo, &Yo); + coords[0].crpix1 = Xo; + coords[0].crpix2 = Yo; + + // resulting fit should have zero constant terms + modfit = CoordsSetCenter (fit, Xo, Yo); + + // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0) + coords[0].pc1_1 = modfit[0].xfit[1][0]; + coords[0].pc1_2 = modfit[0].xfit[0][1]; + coords[0].pc2_1 = modfit[0].yfit[1][0]; + coords[0].pc2_2 = modfit[0].yfit[0][1]; + + // set the polyterm elements + coords[0].polyterms[0][0] = modfit[0].xfit[2][0]; + coords[0].polyterms[1][0] = modfit[0].xfit[1][1]; + coords[0].polyterms[2][0] = modfit[0].xfit[0][2]; + + coords[0].polyterms[0][1] = modfit[0].yfit[2][0]; + coords[0].polyterms[1][1] = modfit[0].yfit[1][1]; + coords[0].polyterms[2][1] = modfit[0].yfit[0][2]; + + // I need to validate Norder + coords[0].polyterms[3][0] = modfit[0].xfit[3][0]; + coords[0].polyterms[4][0] = modfit[0].xfit[2][1]; + coords[0].polyterms[5][0] = modfit[0].xfit[1][2]; + coords[0].polyterms[6][0] = modfit[0].xfit[0][3]; + + coords[0].polyterms[3][1] = modfit[0].yfit[3][0]; + coords[0].polyterms[4][1] = modfit[0].yfit[2][1]; + coords[0].polyterms[5][1] = modfit[0].yfit[1][2]; + coords[0].polyterms[6][1] = modfit[0].yfit[0][3]; + + /* we do not modify crval1,2: these are kept at the default values */ + + /* normalize pc11,etc */ + c11 = coords[0].pc1_1; + c21 = coords[0].pc1_1; + c12 = coords[0].pc1_1; + c22 = coords[0].pc1_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; + + fit_free (modfit); + /* keep the order and type from initial values */ +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/global_stats.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/global_stats.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/global_stats.c (revision 21919) @@ -0,0 +1,44 @@ +# include "relastro.h" + +void global_stats (Catalog *catalog, int Ncatalog) { + + StatType stN, stX, stS, imN, imX, imM, imD, msM, msX, msN, msD; + + initstats ("MEAN"); + + stN = statsStarN (catalog, Ncatalog); + stX = statsStarX (catalog, Ncatalog); + stS = statsStarS (catalog, Ncatalog); + + imN = statsImageN (catalog); + imX = statsImageX (catalog); + imM = statsImageM (catalog); + imD = statsImagedM (catalog); + + msN = statsMosaicN (catalog); + msM = statsMosaicM (catalog); + msD = statsMosaicdM (catalog); + msX = statsMosaicX (catalog); + + fprintf (stderr, "STATS median mean sigma min max Nmeas\n"); + fprintf (stderr, "meas / image: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", imN.median, imN.mean, imN.sigma, imN.min, imN.max, imN.Nmeas); + fprintf (stderr, "meas / mosaic: %7.0f %7.0f %7.0f %7.0f %7.0f %6d\n", msN.median, msN.mean, msN.sigma, msN.min, msN.max, msN.Nmeas); + fprintf (stderr, "meas / star: %7.0f %7.1f %7.1f %7.0f %7.0f %6d\n", stN.median, stN.mean, stN.sigma, stN.min, stN.max, stN.Nmeas); + + fprintf (stderr, "chisq image: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", imX.median, imX.mean, imX.sigma, imX.min, imX.max, imX.Nmeas); + fprintf (stderr, "chisq mosaic: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas); + fprintf (stderr, "chisq star: %7.1f %7.1f %7.1f %7.1f %7.1f %6d\n", stX.median, stX.mean, stX.sigma, stX.min, stX.max, stX.Nmeas); + + fprintf (stderr, "Mcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imM.median, imM.mean, imM.sigma, imM.min, imM.max, imM.Nmeas); + fprintf (stderr, "dMcal image: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", imD.median, imD.mean, imD.sigma, imD.min, imD.max, imD.Nmeas); + + fprintf (stderr, "Mcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msM.median, msM.mean, msM.sigma, msM.min, msM.max, msM.Nmeas); + fprintf (stderr, "dMcal mosaic: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n", msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas); + + fprintf (stderr, "dMrel star: %7.4f %7.4f %7.4f %7.4f %7.4f %6d\n\n", stS.median, stS.mean, stS.sigma, stS.min, stS.max, stS.Nmeas); + + initstats (STATMODE); + +} + + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/initialize.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/initialize.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/initialize.c (revision 21919) @@ -0,0 +1,61 @@ +# include "relastro.h" + +void initialize (int argc, char **argv) { + + int NPHOTCODES; + char *codename, *ptr; + + ConfigInit (&argc, argv); + args (argc, argv); + + /* XXX we need to build a list of accepted photcodes. these will be used by bcatalog to accept or + reject loaded data */ + + Nphotcodes = 0; + photcodes = NULL; + if (PHOTCODE_LIST != NULL) { + NPHOTCODES = 10; + ALLOCATE (photcodes, PhotCode *, NPHOTCODES); + + /* parse the comma-separated list of photcodes */ + strtok_r (PHOTCODE_LIST, ",", &ptr); + while ((codename = strtok_r (NULL, ",", &ptr)) != NULL) { + fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_LIST); + fprintf (stderr, "codename: %s\n", codename); + if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) { + fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename); + exit (1); + } + Nphotcodes ++; + CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10); + } + } + PhotNsec = GetPhotcodeNsecfilt (); + + initstats (STATMODE); + + IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP; + STAR_BAD = ID_STAR_POOR | ID_STAR_FEW; + MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA; + + /* XXX drop irrelevant entries */ + if (SHOW_PARAMS) { + fprintf (stderr, "current parameter settings:\n"); + if (TimeSelect) { + fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", sec_to_date (TSTART), sec_to_date (TSTOP)); + } else { + fprintf (stderr, "TimeSelect: FALSE\n"); + } + fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF); + + fprintf (stderr, "MAG_LIM %lf\n", MAG_LIM); + fprintf (stderr, "STAR_SCATTER %lf\n", STAR_SCATTER); + fprintf (stderr, "IMAGE_SCATTER %lf\n", IMAGE_SCATTER); + fprintf (stderr, "IMAGE_OFFSET %lf\n", IMAGE_OFFSET); + fprintf (stderr, "IMAGE_CATALOG %s\n", ImageCat); + fprintf (stderr, "GSCFILE %s\n", GSCFILE); + fprintf (stderr, "CATDIR %s\n", CATDIR); + exit (0); + } +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/liststats.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/liststats.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/liststats.c (revision 21919) @@ -0,0 +1,119 @@ +# include "relastro.h" + +enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, + M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN}; + +static int statmode; + +void initstats (char *mode) { + + statmode = -1; + if (!strcmp (mode, "MEAN")) statmode = M_MEAN; + if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN; + if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN; + if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN; + if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN; + if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN; + if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN; + + if (statmode == -1) { + fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode); + exit (1); + } +} + +int liststats (double *value, double *dvalue, int N, StatType *stats) { + + int i, ks, ke, Nm; + double Mo, dMo, M, dM, X2, dS, *chi; + + ke = ks = dMo = 0; + + stats[0].Nmeas = N; + stats[0].mean = 0; + stats[0].sigma = 0; + stats[0].error = 0; + stats[0].chisq = 0; + if (N < 1) return (FALSE); + + sortB (value, dvalue, N); + stats[0].median = value[(int)(0.5*N)]; + stats[0].min = value[0]; + stats[0].max = value[N-1]; + + switch (statmode) { + case M_MEDIAN: + ks = 0; + ke = N; + Mo = stats[0].median; + Nm = N; + goto chisq; + break; + case M_MEAN: + case M_WT_MEAN: + ks = 0; + ke = N; + break; + case M_INNER_MEAN: + case M_INNER_WTMEAN: + case M_CHI_INNER_MEAN: + case M_CHI_INNER_WTMEAN: + ks = 0.25*N + 0.50; + ke = 0.75*N + 0.25; + if (N <= 3) { + ks = 0; + ke = N; + } + break; + } + + if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) { + ALLOCATE (chi, double, N); + for (i = 0; i < N; i++) { + chi[i] = (value[i] - stats[0].median) / dvalue[i]; + } + sortD (chi, value, dvalue, N); + free (chi); + } + + /* calculating the per-star offset based on the weighted average */ + M = dM = Nm = 0; + if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) { + for (i = ks; i < ke; i++) { + M += value[i] / SQ (dvalue[i]); + dM += 1.0 / SQ (dvalue[i]); + Nm ++; + } + Mo = M / dM; + dMo = sqrt (1.0 / dM); + } else { + for (i = ks; i < ke; i++) { + M += value[i]; + dM += SQ (dvalue[i]); + Nm ++; + } + Mo = M / (double) Nm; + dMo = sqrt (dM / (double) Nm); + } + + chisq: + /* find sigma and chisq */ + X2 = dS = 0; + for (i = ks; i < ke; i++) { + M = SQ (value[i] - Mo); + dM = SQ (dvalue[i]); + X2 += M / dM; + dS += M; + } + X2 = X2 / Nm; + dS = sqrt (dS / Nm); + + stats[0].mean = Mo; + stats[0].Nmeas = Nm; + stats[0].chisq = X2; + stats[0].sigma = dS; + stats[0].error = dMo; + + return (TRUE); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_catalogs.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_catalogs.c (revision 21919) @@ -0,0 +1,66 @@ +# include "relastro.h" + +Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) { + + int i, j, k, m, Nstar; + Catalog *catalog, *pcatalog, tcatalog; + + if (VERBOSE) fprintf (stderr, "loading catalog data\n"); + + ALLOCATE (catalog, Catalog, skylist[0].Nregions); + + // load data from each region file, only use bright stars + for (i = 0; i < skylist[0].Nregions; i++) { + + pcatalog = subselect ? &tcatalog : &catalog[i]; + + // set up the basic catalog info + pcatalog[0].filename = skylist[0].filename[i]; + pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + pcatalog[0].catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + pcatalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point + pcatalog[0].Nsecfilt = GetPhotcodeNsecfilt (); + + if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename); + exit (1); + } + if (VERBOSE && !pcatalog[0].Nave_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename); + + // select only the brighter stars + if (subselect) { + bcatalog (&catalog[i], &tcatalog); + dvo_catalog_unlock (&tcatalog); + dvo_catalog_free (&tcatalog); + } else { + if (RESET) { + for (j = 0; j < catalog[i].Naverage; j++) { + catalog[i].average[j].code = 0; + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++) { + catalog[i].measure[m+k].flags = 0; + } + } + } + } + } + + // XXX keep this test? + Nstar = 0; + for (i = 0; i < skylist[0].Nregions; i++) { + Nstar += catalog[i].Naverage; + } + if (Nstar < 2) { + fprintf (stderr, "warning: insufficient stars %d\n", Nstar); + } + + // XXX consider only returning the populated catalogs + *Ncatalog = skylist[0].Nregions; + return (catalog); +} + +/* this function loads all relevant catalog files for the first pass. it currently loads the data + read only (SOFT lock) since it assumes the image table has been locked. if we go to the new + addstar locking paradigm, in which the images and catalogs are updated independently, then we may + need to use an XCLD lock here. +*/ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_images.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_images.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/load_images.c (revision 21919) @@ -0,0 +1,34 @@ +# include "relastro.h" + +SkyList *load_images (FITS_DB *db, SkyRegion *region) { + + Image *image, *subset; + int Nimage, Nsubset; + int *LineNumber; + + SkyTable *sky = NULL; + SkyList *skylist = NULL; + + // load the current sky table (layout of all SkyRegions) + sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, SKY_DEPTH, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + + // determine the populated SkyRegions overlapping the requested area + skylist = SkyListByPatch (sky, -1, region); + + // convert database table to internal structure + image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped); + + // select the images which overlap the selected sky regions + subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset); + + gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset); + + initImages (subset, Nsubset); + initMosaics (subset, Nsubset); + + /* unlock, if we can (else, unlocked below) */ + if (!UPDATE) dvo_image_unlock (db); + + return (skylist); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/mkpolyterm.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/mkpolyterm.c (revision 21919) @@ -0,0 +1,159 @@ +# include "relastro.h" + +double factorial (int N) { + + int i; + double F; + + F = 1; + for (i = N; i > 1; i--) { + F *= i; + } + return F; +} + +// XXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ? + +// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y +// for which L,M = 0,0. tol is the allowed error on x,y. +int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo) { + + int i, Nx, Ny; + double R, Xo, Yo, dPos; + double **XdX, **XdY, **YdX, **YdY, **alpha, **beta; + double **xfit, **yfit; + + xfit = fit[0].xfit; + yfit = fit[0].yfit; + + // solve for Xo,Yo where L,M(Xo,Yo) = 0,0 + // start with linear solution for Xo,Yo: + R = (xfit[1][0]*yfit[0][1] - xfit[0][1]*yfit[1][0]); + Xo = (yfit[0][0]*xfit[0][1] - xfit[0][0]*yfit[0][1])/R; + Yo = (xfit[0][0]*yfit[1][0] - yfit[0][0]*xfit[1][0])/R; + + Nx = fit[0].Nterms; + Ny = fit[0].Nterms; + + // iterate to actual solution: requires small non-linear terms + if (fit[0].Norder > 1) { + XdX = poly2d_dx (xfit, Nx, Ny); + XdY = poly2d_dy (xfit, Nx, Ny); + YdX = poly2d_dx (yfit, Nx, Ny); + YdY = poly2d_dy (yfit, Nx, Ny); + + alpha = array_init (2, 2); + beta = array_init (2, 1); + + /* this loop uses the Newton-Raphson method to solve for Xo,Yo + * it needs the high order terms to be small + * Xo,Yo are in pixels; + */ + dPos = tol + 1; + for (i = 0; (dPos > tol) && (i < 20); i++) { + // NOTE: order for alpha is: [y][x] + alpha[0][0] = poly2d_eval (XdX, Nx-1, Ny, Xo, Yo); + alpha[1][0] = poly2d_eval (XdY, Nx, Ny-1, Xo, Yo); + alpha[0][1] = poly2d_eval (YdX, Nx-1, Ny, Xo, Yo); + alpha[1][1] = poly2d_eval (YdY, Nx, Ny-1, Xo, Yo); + + beta[0][0] = poly2d_eval (xfit, Nx, Ny, Xo, Yo); + beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo); + + dgaussj (alpha, 2, beta, 1); + + Xo -= beta[0][0]; + Yo -= beta[1][0]; + dPos = hypot(beta[0][0], beta[1][0]); + } + array_free (alpha, 2); + array_free (beta, 2); + array_free (XdX, Nx - 1); + array_free (XdY, Nx); + array_free (YdX, Nx - 1); + array_free (YdY, Nx); + } + *xo = Xo; + *yo = Yo; + return TRUE; +} + +// convert a transformation L(x,y) to L'(x-xo,y-yo) +CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo) { + + int i, j, Nx, Ny; + double **xPx, **yPx, **xPy, **yPy, **tmp; + + CoordFit *output; + output = fit_init (input->Norder); + + /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j + * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and + * noting that the constant term in each is the coefficient in the case of L(x,y) and is the + * value of L'(-xo,-yo) in the second case. + */ + + Nx = input->Nterms; + Ny = input->Nterms; + + xPx = poly2d_copy (input->xfit, Nx, Ny); + yPx = poly2d_copy (input->yfit, Nx, Ny); + + for (i = 0; i <= input->Nterms; i++) { + xPy = poly2d_copy (xPx, Nx, Ny); + yPy = poly2d_copy (yPx, Nx, Ny); + for (j = 0; j <= input->Nterms; j++) { + output->xfit[i][j] = poly2d_eval (xPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j); + output->yfit[i][j] = poly2d_eval (yPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j); + + // take the next derivative wrt y, catch output (is NULL on last pass) + if (Ny > 0) { + tmp = poly2d_dy(xPy, Nx, Ny); + array_free (xPy, Nx); + xPy = tmp; + + tmp = poly2d_dy(yPy, Nx, Ny); + array_free (yPy, Nx); + yPy = tmp; + } else { + array_free (xPy, Nx); + array_free (yPy, Nx); + } + + Ny --; + } + Ny = input->Nterms; + // take the next derivative wrt x, catch output (is NULL on last pass) + if (Nx > 0) { + tmp = poly2d_dx(xPx, Nx, Ny); + array_free (xPx, Nx); + xPx = tmp; + tmp = poly2d_dx(yPx, Nx, Ny); + array_free (yPx, Nx); + yPx = tmp; + Nx --; + } else { + array_free (xPx, Nx); + array_free (yPx, Nx); + } + } + return output; +} + +/* + Coords uses a rigid sequence for the coefficients: + + x^0 y^0 : crpix1,2 (but note this is applied before higher order terms) + + x^1 y^0 : pc1_1, pc2_1 + x^0 y^1 : pc1_2, pc2_2 + + x^2 y^0 : polyterm[0][0,1] + x^1 y^1 : polyterm[1][0,1] + x^0 y^2 : polyterm[2][0,1] + + x^3 y^0 : polyterm[3][0,1] + x^2 y^1 : polyterm[4][0,1] + x^1 y^2 : polyterm[5][0,1] + x^0 y^3 : polyterm[6][0,1] +*/ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plot_scatter.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plot_scatter.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plot_scatter.c (revision 21919) @@ -0,0 +1,59 @@ +# include "relastro.h" + +void plot_scatter (Catalog *catalog, int Ncatalog) { + +# if (0) + int i, j, k, m, N, Ntot; + float Mrel, Mcal, Mmos, Mgrid; + double *xlist, *ylist, *ilist; + Graphdata graphdata; + + Ntot = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + Ntot += catalog[i].average[j].Nm; + } + } + ALLOCATE (xlist, double, Ntot); + ALLOCATE (ylist, double, Ntot); + ALLOCATE (ilist, double, Ntot); + + N = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + + /* calculate the average value for a single star */ + if (catalog[i].average[j].code & STAR_BAD) continue; + m = catalog[i].average[j].offset; + + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + if (catalog[i].measure[m].flags & MEAS_BAD) continue; + + Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M; + xlist[N] = Mrel; + ylist[N] = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel; + ilist[N] = PhotInst (&catalog[i].measure[m]); + N++; + } + } + } + + plot_defaults (&graphdata); + graphdata.xmin = PlotMmin; + graphdata.xmax = PlotMmax; + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "Mag.png"); + + plot_defaults (&graphdata); + graphdata.ymin = PlotdMmin; + graphdata.ymax = PlotdMmax; + plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "iMag.png"); + free (xlist); + free (ylist); + free (ilist); + +# endif +} + +/* XXX this should become astrometrically relevant */ Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plotstuff.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plotstuff.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/plotstuff.c (revision 21919) @@ -0,0 +1,154 @@ +# include "relastro.h" +# include + +static int Xgraph[5] = {0,0,0,0,0}; +static int active; + +/* +static union { unsigned char c[4]; float f; } f_undef = { {0xff, 0xff, 0xff, 0xfe} }; +# define fUNDEF (f_undef.f) +*/ + +static union { unsigned char c[8]; float d; } d_undef = { {0xff, 0xff, 0xff, 0xfe, 0x00, 0x00, 0x00, 0x00} }; +# define dUNDEF (d_undef.d) + +void XDead () { + signal (SIGPIPE, XDead); + fprintf (stderr, "kapa is dead, must restart\n"); + Xgraph[active] = -1; +} + +int open_graph (int N) { + + Xgraph[N] = KiiOpen ("kapa", NULL); + if (Xgraph[N] < 0) { + fprintf (stderr, "error starting kapa\n"); + return (FALSE); + } + return (TRUE); +} + +void DonePlotting (Graphdata *graphmode, int N) { + + if (Xgraph[N] == 0) return; + + KapaBox (Xgraph[N], graphmode); + return; +} + +void JpegPlot (Graphdata *graphmode, int N, char *filename) { + + if (Xgraph[N] == 0) return; + + KapaPNG (Xgraph[N], filename); + return; +} + +void PSPlot (Graphdata *graphmode, int N, char *filename) { + + if (Xgraph[N] == 0) return; + + KiiPS (Xgraph[N], filename, TRUE, KAPA_PS_NEWPLOT, "default"); + return; +} + +void PrepPlotting (int Npts, Graphdata *graphmode, int N) { + + if (Npts < 1) return; + + if (Xgraph[N] < 1) if (!open_graph(N)) return; + + KapaClear (Xgraph[N], TRUE); + + KiiSendCommand (Xgraph[N], 4, "PLOT"); + + /* send Xgraph the plot details */ + KiiSendMessage (Xgraph[N], "%8d %8d %d %d %d %d %d %f %f", + Npts, graphmode[0].style, + graphmode[0].ptype, graphmode[0].ltype, + graphmode[0].etype, graphmode[0].ebar, graphmode[0].color, + graphmode[0].lweight, graphmode[0].size); + KiiSendMessage (Xgraph[N], "%g %g %g %g", + graphmode[0].xmin, graphmode[0].xmax, + graphmode[0].ymin, graphmode[0].ymax); + + return; +} + +void PlotLabel (char *string, int N) { + + if (Xgraph[N] == 0) return; + + KapaSendLabel (Xgraph[N], string, 2); +} + +void PlotVector (int Npts, double *vect, int mode, int N) { + + float *values; + int i, Nbytes; + + if (Npts < 1) return; + + ALLOCATE (values, float, Npts); + for (i = 0; i < Npts; i++) { + values[i] = vect[i]; + } + + Nbytes = Npts * sizeof (float); + write (Xgraph[N], values, Nbytes); + free (values); + return; +} + +void plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file) { + + int i; + StatType stats; + + stats.min = stats.max = xlist[0]; + for (i = 0; i < N; i++) { + stats.min = MIN (stats.min, xlist[i]); + stats.max = MAX (stats.max, xlist[i]); + } + if (graphdata[0].xmin == dUNDEF) graphdata[0].xmin = 1.05*stats.min - 0.05*stats.max; + if (graphdata[0].xmax == dUNDEF) graphdata[0].xmax = 1.05*stats.max - 0.05*stats.min; + + stats.min = stats.max = ylist[0]; + for (i = 0; i < N; i++) { + stats.min = MIN (stats.min, ylist[i]); + stats.max = MAX (stats.max, ylist[i]); + } + if (graphdata[0].ymin == dUNDEF) graphdata[0].ymin = 1.05*stats.min - 0.05*stats.max; + if (graphdata[0].ymax == dUNDEF) graphdata[0].ymax = 1.05*stats.max - 0.05*stats.min; + + PrepPlotting (N, graphdata, 0); + PlotVector (N, xlist, 0, 0); + PlotVector (N, ylist, 1, 0); + if (label != NULL) PlotLabel (label, 0); + DonePlotting (graphdata, 0); + + if ((file != NULL) && SAVEPLOT) JpegPlot (graphdata, 0, file); + if (PLOTDELAY > 500000) { + fprintf (stdout, "press return\n"); + fscanf (stdin, "%*c"); + } else { + usleep (PLOTDELAY); + } +} + +void plot_defaults (Graphdata *graphdata) { + + graphdata[0].style = 2; + graphdata[0].ptype = 2; + graphdata[0].ltype = 0; + graphdata[0].etype = 0; + graphdata[0].color = black; + graphdata[0].lweight = 0; + graphdata[0].size = 0.5; + + graphdata[0].xmin = dUNDEF; + graphdata[0].xmax = dUNDEF; + graphdata[0].ymin = dUNDEF; + graphdata[0].ymax = dUNDEF; + +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/relastro.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/relastro.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/relastro.c (revision 21919) @@ -0,0 +1,82 @@ +# include "relastro.h" + +int main (int argc, char **argv) { + + int status, Ncatalog; + Catalog *catalog; + FITS_DB db; + + SkyList *skylist = NULL; + + /* get configuration info, args */ + initialize (argc, argv); + + /* register database handle with shutdown procedure */ + set_db (&db); + + /* lock and load the image db table */ + status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT)); + if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename); + if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename); + if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename); + + /* load regions and images based on specified sky patch */ + skylist = load_images (&db, &UserPatch); + + /* load catalog data from region files : subselect only if we are not doing the objects */ + catalog = load_catalogs (skylist, &Ncatalog, (FIT_TARGET != TARGET_OBJECTS)); + + /* match measurements with images, mosaics */ + initImageBins (catalog, Ncatalog); + // initMosaicBins (catalog, Ncatalog); + + findImages (catalog, Ncatalog); + // findMosaics (catalog, Ncatalog); /* also sets Grid values */ + + if (PLOTSTUFF) { + // plot_star_coords (catalog, Ncatalog); + plot_mosaic_fields (catalog); + } + + /* major modes */ + switch (FIT_TARGET) { + case TARGET_OBJECTS: + UpdateObjects (catalog, Ncatalog); + break; + case TARGET_SIMPLE: + UpdateSimple (catalog, Ncatalog); + break; + + case TARGET_CHIPS: + UpdateChips (catalog, Ncatalog); + break; + + case TARGET_MOSAICS: + UpdateMosaic (catalog, Ncatalog); + break; + default: + fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__); + exit (2); + } + + /* put these types of functions inside the update functions + plot_scatter (catalog, Ncatalog); + plot_grid (catalog); + plot_mosaics (); + plot_images (); + plot_stars (catalog, Ncatalog); + plot_chisq (catalog, Ncatalog); + */ + + if (!UPDATE) exit (0); + + if (FIT_TARGET == TARGET_OBJECTS) { + save_catalogs (catalog, Ncatalog); + } else { + dvo_image_update (&db, VERBOSE); + dvo_image_unlock (&db); + } + + exit (0); +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/reload_catalogs.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/reload_catalogs.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/reload_catalogs.c (revision 21919) @@ -0,0 +1,51 @@ +# include "relastro.h" + +void reload_catalogs (SkyList *skylist) { + + int i; + int status; + struct stat filestat; + Catalog catalog; + + if (VERBOSE) fprintf (stderr, "re-loading catalog data\n"); + + /* load data from each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + catalog.filename = skylist[0].filename[i]; + + // only update existing db tables + status = stat (catalog.filename, &filestat); + if ((status == -1) && (errno == ENOENT)) { + if (VERBOSE) fprintf (stderr, "no file %s, skipping\n", catalog.filename); + continue; + } + + catalog.catformat = dvo_catalog_catformat (CATFORMAT); // set the default catformat from config data + catalog.catmode = dvo_catalog_catmode (CATMODE); // set the default catmode from config data + catalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF; + + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename); + exit (1); + } + if (VERBOSE && (catalog.Nave_disk == 0)) { + fprintf (stderr, "no data in %s, skipping\n", catalog.filename); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + continue; + } + + initImageBins (&catalog, 1); + initMosaicBins (&catalog, 1); + + findImages (&catalog, 1); + findMosaics (&catalog, 1); + + dvo_catalog_save (&catalog, VERBOSE); + dvo_catalog_unlock (&catalog); + dvo_catalog_free (&catalog); + + freeImageBins (1); + freeMosaicBins (1); + } +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/save_catalogs.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/save_catalogs.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/save_catalogs.c (revision 21919) @@ -0,0 +1,15 @@ +# include "relastro.h" + +void save_catalogs (Catalog *catalog, int Ncatalog) { + + int i; + + /* load data from each region file */ + for (i = 0; i < Ncatalog; i++) { + + if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename); + dvo_catalog_save (&catalog[i], VERBOSE); + dvo_catalog_unlock (&catalog[i]); + dvo_catalog_free (&catalog[i]); + } +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/select_images.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/select_images.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/select_images.c (revision 21919) @@ -0,0 +1,244 @@ +# include "relastro.h" + +/* this function returns a list of all images which overlap the given SkyList (set of + SkyRegions). All images in the image catalog are tested once, so there is no check that an + image already has been included. LineNum stores the locations in the Image database of the + list of images */ + +typedef struct { + double Xc[5]; + double Yc[5]; + double Rc; + double Dc; +} SkyRegionCoords; + +Image *select_images (SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage) { + + Image *image; + int i, j, k, m, found, nimage, NIMAGE; + int InRange; + double Ri[5], Di[5], Xi[5], Yi[5], dx, dy; + int *line_number; + Coords tcoords; + SkyRegionCoords *skycoords; + + if (skylist[0].Nregions < 1) { + *Nimage = 0; + *LineNumber = NULL; + if (VERBOSE) fprintf (stderr, "no matching sky regions\n"); + return NULL; + } + + // the comparison is made in the catalog local projection. below we set crval1,2 + tcoords.crpix1 = tcoords.crpix2 = 0.0; + tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0; + tcoords.pc1_1 = tcoords.pc2_2 = 1.0; + tcoords.pc1_2 = tcoords.pc2_1 = 0.0; + strcpy (tcoords.ctype, "RA---TAN"); + + ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions); + + /* compare with each region file */ + for (i = 0; i < skylist[0].Nregions; i++) { + + /* we make positional comparisons in the projection of catalog */ + skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin); + skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin); + tcoords.crval1 = skycoords[i].Rc; + tcoords.crval2 = skycoords[i].Dc; + + /* define catalog corners */ + RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords); + RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords); + RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords); + skycoords[i].Xc[4] = skycoords[i].Xc[0]; + skycoords[i].Yc[4] = skycoords[i].Yc[0]; + + dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]); + dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]); + skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy; + skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy; + skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy; + skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy; + skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy; + } + + if (VERBOSE) fprintf (stderr, "finding images\n"); + BuildChipMatch (timage, Ntimage); + + nimage = 0; + NIMAGE = 100; + ALLOCATE (image, Image, NIMAGE); + ALLOCATE (line_number, int, NIMAGE); + + // go through the complete list of images, selecting ones which overlap any region + for (i = 0; i < Ntimage; i++) { + + /* exclude images by photcode, if selected */ + if (Nphotcodes > 0) { + found = FALSE; + for (j = 0; (j < Nphotcodes) && !found; j++) { + if (photcodes[j][0].code == timage[i].photcode) found = TRUE; + } + if (!found) continue; + } + + /* exclude images by time */ + if (TimeSelect) { + if (timage[i].tzero < TSTART) continue; + if (timage[i].tzero > TSTOP) continue; + } + + if (!FindMosaicForImage (timage, Ntimage, i)) continue; + + /* define image corners */ + Xi[0] = 0; Yi[0] = 0; + Xi[1] = timage[i].NX; Yi[1] = 0; + Xi[2] = timage[i].NX; Yi[2] = timage[i].NY; + Xi[3] = 0; Yi[3] = timage[i].NY; + Xi[4] = 0; Yi[4] = 0; + found = FALSE; + + /* transform to ra,dec */ + for (j = 0; j < 5; j++) { + XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords); + } + + /* compare with each region file */ + for (m = 0; (m < skylist[0].Nregions) && !found; m++) { + + /* we make positional comparisons in the projection of catalog */ + tcoords.crval1 = skycoords[m].Rc; + tcoords.crval2 = skycoords[m].Dc; + + /* transform to ra,dec */ + InRange = TRUE; + for (j = 0; (j < 5) && InRange; j++) { + InRange = RD_to_XY (&Xi[j], &Yi[j], Ri[j], Di[j], &tcoords); + } + if (!InRange) continue; + + /* check if image corner inside catalog */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]); + } + /* check if catalog corner inside image */ + for (j = 0; (j < 4) && !found; j++) { + found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]); + } + /* check if edges cross */ + for (j = 0; (j < 4) && !found; j++) { + for (k = 0; (k < 4) && !found; k++) { + found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]); + } + } + if (!found) continue; + + image[nimage] = timage[i]; + /* always allow 'few' images to succeed, if possible */ + if (image[nimage].code & ID_IMAGE_FEW) { + image[nimage].code &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR); + } + if (RESET) { + assignMcal (&image[nimage], (double *) NULL, -1); + image[nimage].dMcal = NO_MAG; + image[nimage].code &= ~ID_IMAGE_POOR; + } + line_number[nimage] = i; + nimage ++; + if (nimage == NIMAGE) { + NIMAGE += 100; + REALLOCATE (image, Image, NIMAGE); + REALLOCATE (line_number, int, NIMAGE); + } + } + } + + if (VERBOSE) fprintf (stderr, "found %d images\n", nimage); + + REALLOCATE (image, Image, MAX (nimage, 1)); + REALLOCATE (line_number, int, MAX (nimage, 1)); + free (skycoords); + + *Nimage = nimage; + *LineNumber = line_number; + return (image); +} + +/* check if line between points 0 and 1 of x1 + crosses line between points 0 and 1 of x2 */ +int edge_check (double *x1, double *y1, double *x2, double *y2) { + + double theta1, theta2; + double Theta1, Theta2; + + theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); + theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + Theta1 = theta1; + Theta2 = theta2; + theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); + theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); + + + if (theta1*theta2 < 0.0) { + return (FALSE); + } + + if (fabs(theta1) < fabs(theta2)) { + return (FALSE); + } + + return (TRUE); + +} + +/* check if point x1,y1 is in box formed by x2[0-4] */ +int corner_check (double *x1, double *y1, double *x2, double *y2) { + + int i; + double theta; + + theta = 0; + + for (i = 0; i < 4; i++) { + theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); + } + if (fabs(theta) > 6) { + return (TRUE); + } else { + return (FALSE); + } +} + +/* returns the opening angle between the three points (2 is in middle) + in range -pi to pi */ + +double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) { + + double dx1, dy1, dx2, dy2, ct, st, theta; + + dx1 = x1 - x2; + dy1 = y1 - y2; + + dx2 = x3 - x2; + dy2 = y3 - y2; + + ct = (dx1*dx2 + dy1*dy2); + st = (dx1*dy2 - dx2*dy1); + + theta = atan2 (st, ct); + + return (theta); + +} + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setExclusions.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setExclusions.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setExclusions.c (revision 21919) @@ -0,0 +1,55 @@ +# include "relastro.h" + +/* XXX I think all of these, except for X,Y selection, are done / can be done in bcatalog */ +int setExclusions (Catalog *catalog, int Ncatalog) { + + int i, j, k, m, Narea, Nnocal; + Coords *coords; + double r, d, x, y; + + Nnocal = Narea = 0; + for (i = 0; i < Ncatalog; i++) { + for (j = 0; j < catalog[i].Naverage; j++) { + m = catalog[i].average[j].offset; + for (k = 0; k < catalog[i].average[j].Nm; k++, m++) { + + /* select measurements by photcode */ + // XXXX is this done when loaded, or is it needed for reload? + // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode); + // if (ecode != photcode[0].code) goto mark_nocal; + + /* select measurements by time */ + if (TimeSelect) { + if (catalog[i].measure[m].t < TSTART) goto mark_nocal; + if (catalog[i].measure[m].t > TSTOP) goto mark_nocal; + } + + /* select measurements by mag limit */ + if (AreaSelect) { + r = catalog[i].average[j].R + catalog[i].measure[m].dR / 3600.0; + d = catalog[i].average[j].D + catalog[i].measure[m].dD / 3600.0; + if ((coords = getCoords (m, i)) == NULL) goto markbad; + RD_to_XY (&x, &y, r, d, coords); + if (x < AreaXmin) goto markbad; + if (x > AreaXmax) goto markbad; + if (y < AreaYmin) goto markbad; + if (y > AreaYmax) goto markbad; + } + continue; + + markbad: + catalog[i].measure[m].flags |= ID_MEAS_AREA; + Narea ++; + continue; + + mark_nocal: + catalog[i].measure[m].flags |= ID_MEAS_NOCAL; + Nnocal ++; + continue; + } + } + } + if (VERBOSE) fprintf (stderr, "%d measurements marked by area\n", Narea); + if (VERBOSE) fprintf (stderr, "%d measurements marked nocal\n", Nnocal); + return (TRUE); +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setMrelFinal.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setMrelFinal.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/setMrelFinal.c (revision 21919) @@ -0,0 +1,154 @@ +# include "relastro.h" + +void setMrelFinal (Catalog *catalog) { + + int i, j, m, ecode; + float *p; + short *q; + + /* if we reset the catalog, reset all the current measurements */ + if (RESET) { + + for (i = 0; i < catalog[0].Naverage; i++) { + p = (PhotPrimary) ? &catalog[0].average[i].M : &catalog[0].secfilt[PhotNsec*i+PhotSec].M; + *p = NO_MAG; + p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM; + *p = NO_MAG; + q = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm; + *q = NO_MAG; + + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + + /* select measurements by photcode */ + ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); + if (ecode != photcode[0].code) continue; + + /* select measurements by time */ + if (TimeSelect) { + if (catalog[0].measure[m].t < TSTART) continue; + if (catalog[0].measure[m].t > TSTOP) continue; + } + + catalog[0].measure[m].Mcal = 0; + catalog[0].measure[m].flags &= 0xff00; + catalog[0].measure[m].flags &= ~ID_MEAS_POOR; + catalog[0].measure[m].flags &= ~ID_MEAS_SKIP; + catalog[0].measure[m].flags &= ~ID_MEAS_AREA; + catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL; + } + } + } + + setExclusions (catalog, 1); /* mark by area */ + + /* set catalog[0].found[i] = FALSE */ + ALLOCATE (catalog[0].found, int, MAX (1, catalog[0].Naverage)); + for (i = 0; i < catalog[0].Naverage; i++) { + catalog[0].found[i] = FALSE; + } + + for (i = 0; i < 5; i++) { + skip_measurements (catalog, i); /* set ID_MEAS_SKIP for measures to be skipped */ + setMrelOutput (catalog, 1, FALSE); /* set Mrel using all measures not skipped */ + clean_measures (catalog, 1, TRUE); /* mark outliers ID_MEAS_POOR */ + setMrelOutput (catalog, 1, TRUE); /* set Mrel using remaining measures */ + } + setMcalOutput (catalog, 1); + + /* clear ID_STAR_POOR, ID_STAR_FEW, ID_MEAS_NOCAL values before writing ??? */ + for (i = 0; i < catalog[0].Naverage; i++) { + catalog[0].average[i].code &= ~ID_STAR_FEW; + catalog[0].average[i].code &= ~ID_STAR_POOR; + m = catalog[0].average[i].offset; + for (j = 0; j < catalog[0].average[i].Nm; j++, m++) { + catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL; + } + } +} + +/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */ +void skip_measurements (Catalog *catalog, int pass) { + + int i, k, m, ecode, d1, d2; + int Ntot, Ntry, Nkeep, Nskip; + float mag; + + Ntot = Ntry = Nskip = Nkeep = 0; + + /* allow measures from images marked POOR and FEW */ + if (pass >= 3) IMAGE_BAD = ID_IMAGE_NOCAL; + + /* allow measures marked as outliers (POOR) and off image region (AREA) */ + if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP; + + /* mark measures which should be ignored on second pass */ + for (i = 0; i < catalog[0].Naverage; i++) { + Ntot += catalog[0].average[i].Nm; + if (catalog[0].found[i]) continue; + + m = catalog[0].average[i].offset; + for (k = 0; k < catalog[0].average[i].Nm; k++, m++) { + Ntry++; + + /* clear SKIP for all measures at first */ + catalog[0].measure[m].flags &= ~ID_MEAS_SKIP; + + /** never use these measurements (wrong photcode, bad time range) */ + /* skipped via NOCAL, don't mark as skipped */ + ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode); + if (ecode != photcode[0].code) continue; + + /* skip measurements by time range */ + if (TimeSelect) { + if (catalog[0].measure[m].t < TSTART) goto skip; + if (catalog[0].measure[m].t > TSTOP) goto skip; + } + + /* skip measurements with sat. dophot values */ + if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip; + + /* skip measurements from BAD images and mosaics */ + if ((getMcal (m, 0)) == NO_MAG) goto skip; + if ((getMmos (m, 0)) == NO_MAG) goto skip; + + /* skip measurements by inst mag limit */ + if ((pass < 4) && ImagSelect) { + mag = PhotInst (&catalog[0].measure[m]); + if (mag < ImagMin) goto skip; + if (mag > ImagMax) goto skip; + } + + d1 = (catalog[0].measure[m].dophot == 1); + d2 = (catalog[0].measure[m].dophot == 2); + + /* skip measurements by measurement error */ + if ((pass < 2) && !(d1 || d2)) goto skip; + + /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */ + if ((pass < 1) && !d1) goto skip; + Nkeep ++; + + continue; + + skip: + catalog[0].measure[m].flags |= ID_MEAS_SKIP; + Nskip ++; + } + } + fprintf (stderr, "pass %d, Ntot: %d, Ntry: %d, Nskip: %d, Nkeep: %d\n", + pass, Ntot, Ntry, Nskip, Nkeep); +} + +/* dophot types: + +1 - star +2 - galaxy +3 - star +4 - weak +5 - edge +7 - weak +9 - weak +10 - sat +*/ + Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/sort.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/sort.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/sort.c (revision 21919) @@ -0,0 +1,182 @@ +# include "relastro.h" + +void sortA (double *X, int N) { + + int l,j,ir,i; + double tX; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + if (--ir == 0) { + X[0] = tX; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + } +} + +void sortB (double *X, double *Y, int N) { + + int l,j,ir,i; + double tX, tY; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + } +} + +void sortC (double *X, double *Y, double *F1, double *F2, int N) { + + int l,j,ir,i; + double tX, tY; + double t1, t2; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + t1 = F1[l]; + t2 = F2[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + t1 = F1[ir]; + F1[ir] = F1[0]; + t2 = F2[ir]; + F2[ir] = F2[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + F1[0] = t1; + F2[0] = t2; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + F1[i] = F1[j]; + F2[i] = F2[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + F1[i] = t1; + F2[i] = t2; + } +} + +void sortD (double *X, double *Y, double *Z, int N) { + + int l,j,ir,i; + double tX, tY, tZ; + + if (N < 2) return; + + 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; + } + } + 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; + } +} Index: /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/write_coords.c =================================================================== --- /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/write_coords.c (revision 21919) +++ /tags/ohana/relastro-1-6-0/Ohana/src/relastro/src/write_coords.c (revision 21919) @@ -0,0 +1,29 @@ +# include "relastro.h" + +# define CD_COORDS 1 + +void write_coords (Header *header, Coords *coords) { + + gfits_modify (header, "CTYPE1", "%s", 1, "RA---TAN"); + gfits_modify (header, "CTYPE2", "%s", 1, "DEC--TAN"); + + gfits_modify (header, "CRVAL1", "%lf", 1, coords[0].crval1); + gfits_modify (header, "CRVAL2", "%lf", 1, coords[0].crval2); + + gfits_modify (header, "CRPIX1", "%lf", 1, coords[0].crpix1); + gfits_modify (header, "CRPIX2", "%lf", 1, coords[0].crpix2); + +# if (CD_COORDS) + gfits_modify (header, "CD1_1", "%le", 1, coords[0].pc1_1 * coords[0].cdelt1); + gfits_modify (header, "CD2_1", "%le", 1, coords[0].pc2_1 * coords[0].cdelt1); + gfits_modify (header, "CD1_2", "%le", 1, coords[0].pc1_2 * coords[0].cdelt2); + gfits_modify (header, "CD2_2", "%le", 1, coords[0].pc2_2 * coords[0].cdelt2); +# else + gfits_modify (header, "CDELT1", "%le", 1, coords[0].cdelt1); + gfits_modify (header, "CDELT2", "%le", 1, coords[0].cdelt2); + gfits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1); + gfits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2); + gfits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1); + gfits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2); +# endif +}