Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/Makefile =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/Makefile (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/Makefile (revision 21563) @@ -0,0 +1,64 @@ +default: gastro2 +help: + @echo "make options: gastro2 (default)" + +include ../../Configure +HOME = $(ROOT)/src/gastro2 +BIN = $(HOME)/bin +LIB = $(HOME)/lib +SRC = $(HOME)/src +MAN = $(HOME)/doc +INC = $(HOME)/include +include ../../Makefile.Common + +LDFLAGS := -ldvo -lFITS -lohana $(LDFLAGS) + +gastro2: $(BIN)/gastro2.$(ARCH) +coordtest: $(BIN)/coordtest.$(ARCH) +extr2mass: $(BIN)/extr2mass.$(ARCH) +install: $(DESTBIN)/gastro2 + +EXTR2MASS = \ +$(SRC)/extr2mass.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/gregions2.$(ARCH).o + +COORDTEST = \ +$(SRC)/coordtest.$(ARCH).o \ +$(SRC)/gpairs.$(ARCH).o \ +$(SRC)/polyfit.$(ARCH).o + +GASTRO2 = \ +$(SRC)/plots.$(ARCH).o \ +$(SRC)/gheader2.$(ARCH).o \ +$(SRC)/gfit2.$(ARCH).o \ +$(SRC)/gpairs.$(ARCH).o \ +$(SRC)/polyfit.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/rotate2.$(ARCH).o \ +$(SRC)/gcenter2.$(ARCH).o \ +$(SRC)/gproject2.$(ARCH).o \ +$(SRC)/grid.$(ARCH).o \ +$(SRC)/lumfunc.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o \ +$(SRC)/misc.$(ARCH).o \ +$(SRC)/gargs.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/gastro2.$(ARCH).o \ +$(SRC)/gstars2.$(ARCH).o \ +$(SRC)/greference2.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/getusnob.$(ARCH).o \ +$(SRC)/getptolemy.$(ARCH).o \ +$(SRC)/getgsc.$(ARCH).o \ +$(SRC)/remove_clumps.$(ARCH).o \ +$(SRC)/rfits.$(ARCH).o \ +$(SRC)/rtext.$(ARCH).o + +$(GASTRO2): $(INC)/gastro2.h +$(COORDTEST): $(INC)/gastro2.h +$(EXTR2MASS): $(INC)/gastro2.h + +$(BIN)/gastro2.$(ARCH): $(GASTRO2) +$(BIN)/coordtest.$(ARCH): $(COORDTEST) +$(BIN)/extr2mass.$(ARCH): $(EXTR2MASS) Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/bin/.cvsignore =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/bin/.cvsignore (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/bin/.cvsignore (revision 21563) @@ -0,0 +1,2 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix +*.linrh Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/doc/ChangeLog.txt =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/doc/ChangeLog.txt (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/doc/ChangeLog.txt (revision 21563) @@ -0,0 +1,20 @@ + + gastro-2-3 2006.10.04 + * converted to gfits APIs (forces libfits 1.6) + * converted to new DVO APIs (forces libdvo 1.3) + * removed unused code + * better error checks in ConfigInit + * fixed luminosity function matching + * added ASCA mode + * added mode to measure ptolemy fill-factor + * cleaned up polynomial fitting / header conversions + * cleaned up fitting code + * added error test for two few / nan fits + + gastro-2-2: + added USNO-B + + gastro-2-1: + added support for dvo table modes / formats (libdvo v1.0) + minor changes to sync with libfits modes (libfits v1.4) + added support for both fits and text input data files Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/include/gastro2.h =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/include/gastro2.h (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/include/gastro2.h (revision 21563) @@ -0,0 +1,207 @@ +# include +# include + +typedef struct { + double R, D; + double P, Q; + double X, Y; + double M, dM; + int type; +} StarData; + +typedef struct { + double dNdM; + double Mo; + double Mmin; + double Mmax; + double Mz; +} LumStats; + +typedef struct { + double angle; + double Xoff; + double Yoff; + double Chi; + double dR; + int N; +} Answer; + +typedef struct { + Header header; /* cmp file header */ + LumStats lum; + Coords coords; /* current best guess for astrometry */ + Answer answer; + + double Area; + StarData *stars; /* array with all star data */ + int N; /* number of stars */ +} CmpCatalog; + +typedef struct { + LumStats lum; + + double Area; + double Moff; + double R0, R1; + double D0, D1; + int N; /* number of stars */ + + StarData *stars; /* array with all star data */ +} RefCatalog; + +typedef struct { + double RA[2], DEC[2]; + double Area; + char *name; +} CatStats; + +typedef struct { + double R, D; + double r, b; +} USNOdata; + +/* global variables, from ConfigInit or args */ +double DEFAULT_RADIUS; +double MINIMUM_RADIUS; +double MAX_ERROR, MAX_NONLINEAR; +double MIN_PRECISE; +double CCD_PC1_1; +double CCD_PC2_2; +double CCD_PC1_2; +double CCD_PC2_1; +double NFIELD; +double SEARCH_RADIUS; +double ROT_ZERO; +double dROT; +double RA_OFFSET, DEC_OFFSET; +double POLE_RA, POLE_DEC; +int POLAR_ALIGNMENT; +int NROT; +int VERBOSE; +int LONEOS_COORDS; +int CATDUMP; +int MATCHDUMP; +int NOMATCHDUMP; +int NEWPHOTCODE; +int MIN_MATCHES; +char *PHOTCODE; +int FLIPX, FLIPY; +int NPOLYTERMS; +char CATDIR[256]; +char CATMODE[16]; /* raw, mef, split, mysql */ +char CATFORMAT[16]; /* internal, elixir, loneos, panstarrs */ +char REFCAT[256]; +char HEADER[256]; +int PLOTSTUFF; +int MAGLIMS; +int NMAX_STARS; +char PhotCodeFile[256]; +int GASTRO_MAX_NSTARS; +int TEXTMODE; +int PTOLEMY_FILL_FACTOR; +int MAGMANUAL; +double MAGLIM_MIN; +double MAGLIM_MAX; +int NGRID_PIX; + +int ASCA; +int FORCE; +double F_RA; +double F_DEC; + +double ASEC_PIX; +char ROUGH_ASTROMETRY[64]; + +/* locations for reference data */ +char GSCFILE[256]; +char GSCDIR[256]; +char CATDIR[256]; +char USNO_A_DIR[256]; +char USNO_B_DIR[256]; +char TWO_MASS_DIR[256]; +char ASTROM_CATDIR[256]; +char LONEOS_REGION_FILE[256]; + +typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, ebar, color; + double lweight, size; +} Graphdata; + +StarData *rtext (FILE *f, int *nstars); +StarData *rfits (FILE *f, int *nstars); + +StarData *remove_clumps (StarData *instars, int *nstars, int NX, int NY); +void ConfigInit (int *argc, char **argv); +void DonePlotting (Graphdata *graphmode, int N); +void PlotReset (int N); +void PlotVector (int Npts, float *vect, int mode, int N); +void PrepPlotting (int Npts, Graphdata *graphmode, int N); +void XDead (); +void add_to_regions (CatStats *area); +void ahelp (); +void area_of_region (CatStats *region); +double area_of_skyregion (SkyRegion *region); +void args (int *argc, char **argv, Coords *coords); +void define_region (CatStats *catstats, CmpCatalog *Target); +int dms_to_ddd (double *Value, char *string); +void dump_coords (CmpCatalog *Target); +void fill_lumfunc (StarData *stars, int N, float *lbin, float *bin, int *nb); +int find_dec_bands (CatStats *area); +void fit_add (double x1, double y1, double x2, double y2, double wt); +int fit_adjust (Coords *coords); +void fit_apply (double *x, double *y, double X, double Y); +void fit_eval (); +void fit_init (int order); +void fit_lum_bin (double *x, double *y, int N, double *C0, double *C1); +void fit_norm (); +double fit_scat (StarData *st, StarData *sr, Coords *coords); +int gaussj (double **a, int n, double **b, int m); +void gcenter (CmpCatalog *Target, RefCatalog *Ref); +int get_luminosity_func (StarData *stars, int N, LumStats *lum); +int getptolemy (CatStats *catstats, RefCatalog *Ref); +int getusno (CatStats *catstats, RefCatalog *Ref); +int getusnob (CatStats *catstats, RefCatalog *Ref, double epoch); +int getgsc (CatStats *catstats, RefCatalog *Ref); + +void gfit (CmpCatalog *Target, RefCatalog *Ref, int order); +void gheader (char *file, CmpCatalog *Target); +void gproject (CmpCatalog *Target, RefCatalog *Ref, RefCatalog *Subset); +void greference (CmpCatalog *Target, RefCatalog *Ref); +void grid (CmpCatalog *Target, RefCatalog *Subset, Answer *answer); +int gridbin (double dX, double dY); +void gridfree (); +void gridinit (double XMIN, double XMAX, double YMIN, double YMAX, int Nr, int Nt); +void gstars (char *filename, CmpCatalog *Target); +void hh_hms (double hh, int *hr, int *mn, double *sc); +void hms_format (char *line, double value); +void init_regions (); +int load_ra_blocks (int Ndec, CatStats *area); +int mk_polyterm (int n, int m, int norder); +int mk_vector (int n, int m, int norder); +int open_graph (int N); + +void pair_add (int i1, int i2); +void pair_init (); +int pair_lists (int **index1, int **index2); + +int parse_GSC_line (CatStats *tregion, char *line); +int plot_addpt_gridplot (double x, double y); +void plot_done_gridplot (); +void plot_fullfield (CmpCatalog *Target, RefCatalog *Ref); +void plot_fullfield_pairs (float *x, float *y, int n); +void plot_gridpts (double *pts, int Npts); +void plot_init_gridplot (); +void plot_lumfunc (CmpCatalog *Target, RefCatalog *Ref); +void plot_resid (StarData *st, StarData *sr, Coords *coords); +void plot_resid_init (int version, double xmax); +void plot_resid_plot (int version, float *xvect, float *yvect, int Nvect); + +void rotate (RefCatalog *Subset, RefCatalog *Ref, double angle); +void set_catalog (char *catdir); +void sort (double *X, int N); +void sort_lists (double *X, double *Y, int *S, int N); +void sort_lum (double *R, double *X, double *Y, int N); +void sort_stars_X (StarData *stars, int N); +void sort_stars_mag (StarData *stars, int N); +int str_to_radec (double *ra, double *dec, char *str1, char *str2); Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/ConfigInit.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/ConfigInit.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/ConfigInit.c (revision 21563) @@ -0,0 +1,71 @@ +# include "gastro2.h" + +void ConfigInit (int *argc, char **argv) { + + char *config, *file; + + /*** load configuration info ***/ + file = SelectConfigFile (argc, argv, "ptolemy"); + config = LoadConfigFile (file); + if (config == (char *) NULL) { + fprintf (stderr, "ERROR: can't find configuration file %s\n", file); + if (file != (char *) NULL) free (file); + exit (1); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + /* default values for config variables: used if key is missing from the config file */ + strcpy (ROUGH_ASTROMETRY, "header"); + + ScanConfig (config, "CCD_PC1_1", "%lf", 0, &CCD_PC1_1); // guess if WCS is missing + ScanConfig (config, "CCD_PC2_2", "%lf", 0, &CCD_PC2_2); // guess if WCS is missing + ScanConfig (config, "CCD_PC1_2", "%lf", 0, &CCD_PC1_2); // guess if WCS is missing + ScanConfig (config, "CCD_PC2_1", "%lf", 0, &CCD_PC2_1); // guess if WCS is missing + ScanConfig (config, "ASEC_PIX", "%lf", 0, &ASEC_PIX); // guess if WCS is missing + ScanConfig (config, "NFIELD", "%lf", 0, &NFIELD); // search region *padding* in field units + ScanConfig (config, "NGRID_PIX", "%d", 0, &NGRID_PIX); // resolution of grid search (pixels) + ScanConfig (config, "NPOLYTERMS", "%d", 0, &NPOLYTERMS); // high-order fit terms (2 or 3) + ScanConfig (config, "ROT_ZERO", "%lf", 0, &ROT_ZERO); // rotation search region + ScanConfig (config, "dROT", "%lf", 0, &dROT); // rotation search region + ScanConfig (config, "NROT", "%d", 0, &NROT); // rotation search region + ScanConfig (config, "POLAR_ALIGNMENT", "%d", 0, &POLAR_ALIGNMENT); // apply polar alignment correction + ScanConfig (config, "POLAR_AXIS_RA", "%lf", 0, &POLE_RA); // true coords of pole (should be HA, not RA) + ScanConfig (config, "POLAR_AXIS_DEC", "%lf", 0, &POLE_DEC); // true coords of pole + ScanConfig (config, "RA_OFFSET", "%lf", 0, &RA_OFFSET); // ?? not well defined (should be euler angle) + ScanConfig (config, "DEC_OFFSET", "%lf", 0, &DEC_OFFSET); // ?? not well defined (should be euler angle) + + /* possible sources of astrometric reference data */ + if (!ScanConfig (config, "USNO_A_DIR", "%s", 0, USNO_A_DIR)) { // location of USNO A data (USNO_CDROM in gastro) + ScanConfig (config, "USNO_CDROM", "%s", 0, USNO_A_DIR); // alternate location of USNO A data + } + ScanConfig (config, "USNO_B_DIR", "%s", 0, USNO_B_DIR); // location of USNO B ref data + ScanConfig (config, "GSCDIR", "%s", 0, GSCDIR); // location of HST GSC ref data + ScanConfig (config, "2MASS_DIR", "%s", 0, TWO_MASS_DIR); // location of 2MASS ref data + ScanConfig (config, "ASTROM_CATDIR", "%s", 0, ASTROM_CATDIR); // location of ptolemy-format ref data + + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); // location of sky table + ScanConfig (config, "ASTRO_REFCAT", "%s", 0, REFCAT); // which astrometry catalog to use + ScanConfig (config, "ROUGH_ASTROMETRY", "%s", 0, ROUGH_ASTROMETRY); // where to get initial guess (header, config) + ScanConfig (config, "GASTRO_MAX_NSTARS", "%d", 0, &GASTRO_MAX_NSTARS); // max number of stars from image to fit + ScanConfig (config, "GASTRO_MAX_MAG_ERROR", "%lf", 0, &MAX_ERROR); // S/N limit on image stars used in fit + + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); // not used + + if (NFIELD <= 0.0) { + fprintf (stderr, "NFIELD is not sensible: choose a non-zero number\n"); + exit (1); + } + + if (!GASTRO_MAX_NSTARS) GASTRO_MAX_NSTARS = 300; + if (!MAX_ERROR) MAX_ERROR = 0.2; + if (!NGRID_PIX) NGRID_PIX = 50.0; + if (!NFIELD) NFIELD = 0.1; + + if (strcasecmp (ROUGH_ASTROMETRY, "header") && + strcasecmp (ROUGH_ASTROMETRY, "config")) { + fprintf (stderr, "ROUGH_ASTROMETRY must be one of: header, config\n"); + exit (1); + } + free (config); + free (file); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/coordtest.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/coordtest.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/coordtest.c (revision 21563) @@ -0,0 +1,142 @@ +# include "gastro2.h" + +# define A00 +500.00 +# define A10 +4.00 +# define A01 +1.00 +# define A20 -0.001 +# define A11 +0.000 +# define A02 -0.001 +# define A30 +0.00002 +# define A21 +0.00001 +# define A12 -0.00001 +# define A03 -0.00002 + +# define B00 400.00 +# define B10 -1.00 +# define B01 4.00 +# define B20 -0.001 +# define B11 0.000 +# define B02 -0.001 +# define B30 -0.00002 +# define B21 -0.00001 +# define B12 +0.00001 +# define B03 +0.00002 + +int main (int argc, char **argv) { + + /* generate a set of fake data and send to fitter */ + int i, N, order; + double x, y, z; + double L[50000], M[50000], X[50000], Y[50000]; + Coords coords; + FILE *f; + + if (argc != 2) { + fprintf (stderr, "USAGE: coordtest (order)\n"); + exit (2); + } + + order = atoi (argv[1]); + fit_init (order); + N = 0; + + switch (order) { + case 1: + for (x = -10.0; x < 10.1; x += 0.1) { + for (y = -10.0; y < 10.1; y += 0.1) { + + L[N] = A00 + A10*x + A01*y; + M[N] = B00 + B10*x + B01*y; + X[N] = x; + Y[N] = y; + fit_add (x, y, L[N], M[N], 1.0); + N++; + } + } + break; + case 2: + for (x = -10.0; x < 10.1; x += 0.1) { + for (y = -10.0; y < 10.1; y += 0.1) { + + L[N] = A00 + A10*x + A01*y + A20*x*x + A11*x*y + A02*y*y; + M[N] = B00 + B10*x + B01*y + B20*x*x + B11*x*y + B02*y*y; + X[N] = x; + Y[N] = y; + fit_add (x, y, L[N], M[N], 1.0); + N++; + } + } + break; + case 3: + for (x = -10.0; x < 10.1; x += 0.1) { + for (y = -10.0; y < 10.1; y += 0.1) { + + L[N] = A00 + A10*x + A01*y + A20*x*x + A11*x*y + A02*y*y + A30*x*x*x + A21*x*x*y + A12*x*y*y + A03*y*y*y; + M[N] = B00 + B10*x + B01*y + B20*x*x + B11*x*y + B02*y*y + B30*x*x*x + B21*x*x*y + B12*x*y*y + B03*y*y*y; + X[N] = x; + Y[N] = y; + fit_add (x, y, L[N], M[N], 1.0); + N++; + } + } + break; + } + + fit_eval (); + + strcpy (coords.ctype, "RA---PLY"); + coords.crval1 = 0.0; + coords.crval2 = 0.0; + + fit_adjust (&coords); + + fprintf (stderr, "CTYPE: %s\n", coords.ctype); + fprintf (stderr, "CRVAL: %f %f\n", coords.crval1, coords.crval2); + fprintf (stderr, "CRPIX: %f %f\n", coords.crpix1, coords.crpix2); + fprintf (stderr, "CDELT: %f %f\n", coords.cdelt1, coords.cdelt2); + fprintf (stderr, "PC1_j: %f %f\n", coords.pc1_1, coords.pc1_2); + fprintf (stderr, "PC2_j: %f %f\n", coords.pc2_1, coords.pc2_2); + + fprintf (stderr, "Npolyterms: %d\n", coords.Npolyterms); + for (i = 0; i < 7; i++) { + fprintf (stderr, "%f %f\n", coords.polyterms[i][0], coords.polyterms[i][1]); + } + + f = fopen ("test.dat", "w"); + { + double Lo, Mo, dL, dL2, dM, dM2; + double Xo, Yo, dX, dX2, dY, dY2; + double Lx, Mx; + + dL = dL2 = dM = dM2 = 0.0; + dX = dX2 = dY = dY2 = 0.0; + for (i = 0; i < N; i++) { + XY_to_RD (&Lo, &Mo, X[i], Y[i], &coords); + dL += (L[i] - Lo); + dM += (M[i] - Mo); + dL2 += SQ(L[i] - Lo); + dM2 += SQ(M[i] - Mo); + + RD_to_XY (&Xo, &Yo, L[i], M[i], &coords); + dX += (X[i] - Xo); + dY += (Y[i] - Yo); + dX2 += SQ(X[i] - Xo); + dY2 += SQ(Y[i] - Yo); + + // fit_apply (&Lx, &Mx, X[i], Y[i]); + // fprintf (stderr, "%f,%f -> %f,%f | %f,%f : %f, %f\n", X[i], Y[i], Lx, Mx, L[i], M[i], Lo, Mo); + + fprintf (f, "%f %f : %f %f :: %f %f : %f %f\n", + X[i], Y[i], L[i], M[i], Lo, Mo, Xo, Yo); + } + fclose (f); + + fprintf (stderr, "dL: %f\n", sqrt(fabs(dL2/N - SQ(dL/N)))); + fprintf (stderr, "dM: %f\n", sqrt(fabs(dM2/N - SQ(dM/N)))); + fprintf (stderr, "dX: %f\n", sqrt(fabs(dX2/N - SQ(dX/N)))); + fprintf (stderr, "dY: %f\n", sqrt(fabs(dY2/N - SQ(dY/N)))); + } + + exit (0); +} + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gargs.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gargs.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gargs.c (revision 21563) @@ -0,0 +1,156 @@ +# include "gastro2.h" +# define NARGS 2 /* minimum is: gastro catalog */ + +void ahelp () { + + fprintf (stderr, "gastro -- astrometry for LONEOS\n"); + + fprintf (stderr, " USAGE: gastro pixscale filename"); + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -v (verbose mode)\n"); + fprintf (stderr, " -dump (dump catalog stars, don't complete astrometry)\n"); + fprintf (stderr, " -mdmp (dump matched catalog stars)\n"); + fprintf (stderr, "\n"); + exit (0); + +} + +void args (int *argc, char **argv, Coords *coords) { + + int N; + + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + ahelp (); + } + + VERBOSE = FALSE; + if ((N = get_argument (*argc, argv, "-v"))) { + VERBOSE = TRUE; + remove_argument (N, argc, argv); + } + + /* force read of image database with mismatched NSTARS & size */ + TEXTMODE = FALSE; + if ((N = get_argument (*argc, argv, "-textmode"))) { + TEXTMODE = TRUE; + remove_argument (N, argc, argv); + } + + /* force read of image database with mismatched NSTARS & size */ + PTOLEMY_FILL_FACTOR = FALSE; + if ((N = get_argument (*argc, argv, "-ptolemy-fill-factor"))) { + PTOLEMY_FILL_FACTOR = TRUE; + remove_argument (N, argc, argv); + } + + PLOTSTUFF = FALSE; + if ((N = get_argument (*argc, argv, "-plot"))) { + PLOTSTUFF = TRUE; + remove_argument (N, argc, argv); + } + + MAGLIMS = TRUE; + if ((N = get_argument (*argc, argv, "-maglims"))) { + MAGLIMS = FALSE; + remove_argument (N, argc, argv); + } + if ((N = get_argument (*argc, argv, "-magrange"))) { + MAGLIMS = TRUE; + MAGMANUAL = TRUE; + remove_argument (N, argc, argv); + MAGLIM_MIN = atof (argv[N]); + remove_argument (N, argc, argv); + MAGLIM_MAX = atof (argv[N]); + remove_argument (N, argc, argv); + } + + NMAX_STARS = 300; + if ((N = get_argument (*argc, argv, "-nstars"))) { + remove_argument (N, argc, argv); + NMAX_STARS = atof (argv[N]); + remove_argument (N, argc, argv); + } + + HEADER[0] = 0; + if ((N = get_argument (*argc, argv, "-header"))) { + remove_argument (N, argc, argv); + strcpy (HEADER, argv[N]); + remove_argument (N, argc, argv); + } + HEADER[0] = 0; + if ((N = get_argument (*argc, argv, "-head"))) { + remove_argument (N, argc, argv); + strcpy (HEADER, argv[N]); + remove_argument (N, argc, argv); + } + + FLIPX = FALSE; + if ((N = get_argument (*argc, argv, "-fx"))) { + FLIPX = TRUE; + remove_argument (N, argc, argv); + } + + FLIPY = FALSE; + if ((N = get_argument (*argc, argv, "-fy"))) { + FLIPY = TRUE; + remove_argument (N, argc, argv); + } + + FORCE = FALSE; + if ((N = get_argument (*argc, argv, "-coords"))) { + FORCE = TRUE; + remove_argument (N, argc, argv); + F_RA = atof (argv[N]); + remove_argument (N, argc, argv); + F_DEC = atof (argv[N]); + remove_argument (N, argc, argv); + } + + /** XXX temporary trick to deal with the very wide-field ASCA images + this alters the definition of the reference field boundaries */ + ASCA = FALSE; + if ((N = get_argument (*argc, argv, "-asca"))) { + ASCA = TRUE; + remove_argument (N, argc, argv); + } + + CATDUMP = FALSE; + if ((N = get_argument (*argc, argv, "-dump"))) { + CATDUMP = TRUE; + remove_argument (N, argc, argv); + } + + MATCHDUMP = FALSE; + if ((N = get_argument (*argc, argv, "-mdmp"))) { + MATCHDUMP = TRUE; + remove_argument (N, argc, argv); + } + + NOMATCHDUMP = FALSE; + if ((N = get_argument (*argc, argv, "-cdmp"))) { + NOMATCHDUMP = TRUE; + remove_argument (N, argc, argv); + } + + // XXX these options are being ignored + NEWPHOTCODE = FALSE; + if ((N = get_argument (*argc, argv, "-p"))) { + NEWPHOTCODE = TRUE; + remove_argument (N, argc, argv); + PHOTCODE = strcreate (argv[N]); + remove_argument (N, argc, argv); + } + LONEOS_COORDS = FALSE; + if ((N = get_argument (*argc, argv, "-loneos"))) { + LONEOS_COORDS = TRUE; + remove_argument (N, argc, argv); + } + + if (*argc != NARGS) { + fprintf (stderr, "USAGE: gastro filename\n"); + exit (0); + } + +} + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gastro2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gastro2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gastro2.c (revision 21563) @@ -0,0 +1,99 @@ +# include "gastro2.h" + +int main (int argc, char **argv) { + + int i; + RefCatalog Ref; + CmpCatalog Target; + + /* start_timer (); */ + + ConfigInit (&argc, argv); + args (&argc, argv, &Target.coords); + + gstars (argv[1], &Target); + greference (&Target, &Ref); + gcenter (&Target, &Ref); + + for (i = 0; i < 3; i++) { + gfit (&Target, &Ref, 1); + fprintf (stderr, "precision: %f\n", Target.answer.dR / sqrt(Target.answer.N)); + if ((Target.answer.N < 2) || isnan(Target.answer.dR)) { + fprintf (stderr, "ERROR: bad fit\n"); + exit (1); + } + if (VERBOSE) { + fprintf (stderr, "%s\n", Target.coords.ctype); + fprintf (stderr, "%f %f\n", Target.coords.crval1, Target.coords.crval2); + fprintf (stderr, "%f %f\n", Target.coords.crpix1, Target.coords.crpix2); + fprintf (stderr, "%f %f\n", Target.coords.pc1_1, Target.coords.pc1_2); + fprintf (stderr, "%f %f\n", Target.coords.pc2_1, Target.coords.pc2_2); + fprintf (stderr, "%f %f\n", Target.coords.cdelt1, Target.coords.cdelt2); + } + + } + + for (i = 0; i < 5; i++) { + gfit (&Target, &Ref, MIN (MAX (1, NPOLYTERMS), 3)); + fprintf (stderr, "precision: %f\n", Target.answer.dR / sqrt(Target.answer.N)); + if ((Target.answer.N < 2) || isnan(Target.answer.dR)) { + fprintf (stderr, "ERROR: bad fit\n"); + exit (1); + } + } + + if (VERBOSE) { + fprintf (stderr, "%s\n", Target.coords.ctype); + fprintf (stderr, "%f %f\n", Target.coords.crval1, Target.coords.crval2); + fprintf (stderr, "%f %f\n", Target.coords.crpix1, Target.coords.crpix2); + fprintf (stderr, "%f %f\n", Target.coords.pc1_1, Target.coords.pc1_2); + fprintf (stderr, "%f %f\n", Target.coords.pc2_1, Target.coords.pc2_2); + fprintf (stderr, "%f %f\n", Target.coords.cdelt1, Target.coords.cdelt2); + } + + gheader (argv[1], &Target); + + /* print_timer (); */ + + fprintf (stderr, "SUCCESS\n"); + exit (0); +} + + /* + load config & args + + load stars & header from cmp file + (this should include filtering based on CONFIG data: + limit number of stars, limit types, etc) + + identify and load reference catalog(s) + - load all catalogs available? + - keep ra, dec, mag + + project catalog to initial guess + + find simple x, y offset in limited number of rotations + + reload reference catalog if dx, dy large + + project to new guess + + fit on star-by-star basis + - include weighting by dmag + - downweight by Nmatch to each star + - iterate a few times? + + project to new guess (why is this not part of the routines?) + + try a higher order fit + + adjust header, re-write file + + TODO: + + 1) downweight points by Nmatch + 2) better criterion for success + 3) better criterion for failure + 4) polyterms are broken + + */ Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gcenter2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gcenter2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gcenter2.c (revision 21563) @@ -0,0 +1,67 @@ +# include "gastro2.h" + +void gcenter (CmpCatalog *Target, RefCatalog *Ref) { + + int i, N, Imin; + double angle, ChiMin; + Answer *answer; + RefCatalog Subset; + + gproject (Target, Ref, &Subset); + if (PLOTSTUFF) plot_fullfield (Target, &Subset); + + fprintf (stderr, "Target: %d, Ref: %d\n", Target[0].N, Subset.N); + dump_coords (Target); + + ALLOCATE (answer, Answer, 2*NROT + 1); + + N = 0; + for (angle = ROT_ZERO-dROT*NROT; angle <= ROT_ZERO+dROT*NROT; angle += dROT, N++) { + + if (N == 2*NROT + 1) { + fprintf (stderr, "ERROR in logic: Nanswer > 2*NROT+1 (%d, %d)\n", N, 2*NROT+1); + exit (1); + } + answer[N].angle = angle; + grid (Target, &Subset, &answer[N]); + } + + Imin = 0; + ChiMin = answer[0].Chi; + for (i = 0; i < N; i++) { + if (answer[i].Chi < ChiMin) { + Imin = i; + ChiMin = answer[i].Chi; + } + } + + fprintf (stderr, "best solution: angle: %6.1f, (%6.1f,%6.1f) - %10.8f for %d pairs\n", + answer[Imin].angle, answer[Imin].Xoff, answer[Imin].Yoff, answer[Imin].Chi, answer[Imin].N); + + Target[0].answer = answer[Imin]; + + /* adjust original coordinates for new center */ + { + + double cs, sn; + double pc11, pc12, pc21, pc22; + double Xo, Yo; + + cs = cos(RAD_DEG*answer[Imin].angle); sn = sin(RAD_DEG*answer[Imin].angle); + + pc11 = Target[0].coords.pc1_1; + pc12 = Target[0].coords.pc1_2; + pc21 = Target[0].coords.pc2_1; + pc22 = Target[0].coords.pc2_2; + + Target[0].coords.pc1_1 = pc11*cs - pc21*sn; + Target[0].coords.pc1_2 = pc11*sn + pc12*cs; + Target[0].coords.pc2_1 = pc21*cs - pc22*sn; + Target[0].coords.pc2_2 = pc21*sn + pc22*cs; + + Xo = Target[0].coords.crpix1; + Yo = Target[0].coords.crpix2; + Target[0].coords.crpix1 = answer[Imin].Xoff + cs*Xo - sn*Yo; + Target[0].coords.crpix2 = answer[Imin].Yoff + sn*Xo + cs*Yo; + } +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getgsc.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getgsc.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getgsc.c (revision 21563) @@ -0,0 +1,95 @@ +# include "gastro2.h" + +StarData *rd_gsc (char *filename, int *Nstars); + +int getgsc (CatStats *catstats, RefCatalog *Ref) { + + int i, j, k, Ns, Ngsc; + StarData *gsc; + SkyList *skylist; + SkyTable *sky; + SkyRegion patch; + + Ref[0].N = 0; + Ref[0].Area = 0; + ALLOCATE (Ref[0].stars, StarData, 1); + + patch.Rmin = catstats[0].RA[0]; + patch.Rmax = catstats[0].RA[1]; + patch.Dmin = catstats[0].DEC[0]; + patch.Dmax = catstats[0].DEC[1]; + + /* load regions from GSC table, restrict to patch */ + sky = SkyTableFromGSC (GSCFILE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (sky, GSCDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &patch); + + for (i = 0; i < skylist[0].Nregions; i++) { + gsc = rd_gsc (skylist[0].filename[i], &Ngsc); + + Ns = Ref[0].N; + Ref[0].N += Ngsc; + Ref[0].Area += area_of_skyregion (skylist[0].regions[i]); + + REALLOCATE (Ref[0].stars, StarData, MAX (1, Ref[0].N)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + Ref[0].stars[k].R = gsc[j].R; + Ref[0].stars[k].D = gsc[j].D; + Ref[0].stars[k].M = gsc[j].M; + } + free (gsc); + } + SkyTableFree (sky); + + Ref[0].R0 = catstats[0].RA[0]; + Ref[0].R1 = catstats[0].RA[1]; + Ref[0].D0 = catstats[0].DEC[0]; + Ref[0].D1 = catstats[0].DEC[1]; + + /* calculate luminosity function of stars */ + get_luminosity_func (Ref[0].stars, Ref[0].N, &Ref[0].lum); + + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Ref[0].N); + return (TRUE); +} + +# define BYTES_STAR 23 +# define BLOCK 1000 +StarData *rd_gsc (char *filename, int *Nstars) { + + StarData *stars; + int i, NSTAR, nstar, Nbytes, nbytes; + char *buffer; + FILE *f; + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find catalog file %s\n", filename); + exit (1); + } + + nstar = 0; + NSTAR = 1000; + ALLOCATE (stars, StarData, NSTAR); + + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR)); + Nbytes = BLOCK*BYTES_STAR; + + while ((nbytes = fread (buffer, 1, Nbytes, f)) > 0) { + for (i = 0; i < nbytes / BYTES_STAR; i++) { + dparse (&stars[nstar].R, 1, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].D, 2, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].M, 3, &buffer[i*BYTES_STAR]); + nstar++; + if (nstar == NSTAR) { + NSTAR += 1000; + REALLOCATE (stars, StarData, NSTAR); + } + } + } + + free (buffer); + + *Nstars = nstar; + return (stars); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getptolemy.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getptolemy.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getptolemy.c (revision 21563) @@ -0,0 +1,115 @@ +# include "gastro2.h" +double catalog_area (Average *average, int Naverage, SkyRegion *region); + +int getptolemy (CatStats *catstats, RefCatalog *Ref) { + + int i, j, k, Ns; + double FracArea; + Catalog catalog; + SkyList *skylist; + SkyTable *sky; + SkyRegion patch; + + Ref[0].N = 0; + Ref[0].Area = 0; + ALLOCATE (Ref[0].stars, StarData, 1); + + patch.Rmin = catstats[0].RA[0]; + patch.Rmax = catstats[0].RA[1]; + patch.Dmin = catstats[0].DEC[0]; + patch.Dmax = catstats[0].DEC[1]; + + /* load regions from GSC table, restrict to patch */ + sky = SkyTableLoadOptimal (CATDIR, NULL, GSCFILE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &patch); + + for (i = 0; i < skylist[0].Nregions; i++) { + // set the parameters which guide catalog open/load/create + catalog.filename = skylist[0].filename[i]; + catalog.catflags = LOAD_AVES | LOAD_MEAS; + catalog.Nsecfilt = 0; + + // an error exit status here is a significant error + if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) { + fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename); + exit (2); + } + dvo_catalog_unlock (&catalog); + + // Nave_disk == 0 implies an empty catalog file + // for only_match, skip empty catalogs + if (!catalog.Nave_disk) { + dvo_catalog_free (&catalog); + continue; + } + + // this measurement adjusts a DVO database for partial coverage + // XXX the correction is bogus if the sky density of the catalog + // is too low (<< 100 stars per field) + FracArea = 1.0; + if (PTOLEMY_FILL_FACTOR) { + FracArea = catalog_area (catalog.average, catalog.Naverage, skylist[0].regions[i]); + } + + Ns = Ref[0].N; + Ref[0].N += catalog.Naverage; + Ref[0].Area += FracArea * area_of_skyregion (skylist[0].regions[i]); + + REALLOCATE (Ref[0].stars, StarData, MAX (1, Ref[0].N)); + for (k = Ns, j = 0; j < catalog.Naverage; k++, j++) { + Ref[0].stars[k].R = catalog.average[j].R; + Ref[0].stars[k].D = catalog.average[j].D; + Ref[0].stars[k].M = catalog.measure[catalog.average[j].offset].M_PS; + } + dvo_catalog_free (&catalog); + } + + Ref[0].R0 = catstats[0].RA[0]; + Ref[0].R1 = catstats[0].RA[1]; + Ref[0].D0 = catstats[0].DEC[0]; + Ref[0].D1 = catstats[0].DEC[1]; + + /* calculate luminosity function of stars */ + get_luminosity_func (Ref[0].stars, Ref[0].N, &Ref[0].lum); + + if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", Ref[0].N); + return (TRUE); +} + +double catalog_area (Average *average, int Naverage, SkyRegion *region) { + + int i, xb, yb, Nb; + int bin[10][10]; + double frac, Rmin, Rmax, Dmin, Dmax, dR, dD; + + Rmin = region[0].Rmin; + Rmax = region[0].Rmax; + Dmin = region[0].Dmin; + Dmax = region[0].Dmax; + dR = Rmax - Rmin; + dD = Dmax - Dmin; + + for (xb = 0; xb < 10; xb++) { + for (yb = 0; yb < 10; yb++) { + bin[xb][yb] = 0; + } + } + + for (i = 0; i < Naverage; i++) { + xb = MAX (MIN (0, 10 * (average[i].R - Rmin) / dR), 9); + yb = MAX (MIN (0, 10 * (average[i].D - Dmin) / dD), 9); + bin[xb][yb] ++; + } + + Nb = 0; + for (xb = 0; xb < 10; xb++) { + for (yb = 0; yb < 10; yb++) { + if (bin[xb][yb]) Nb ++; + } + } + + frac = Nb / 100.0; + + return (frac); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusno.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusno.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusno.c (revision 21563) @@ -0,0 +1,138 @@ +# include "gastro2.h" +# define NZONE 24 + +int SPDzone[] = { + 0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, + 675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350}; + +int USNOdisk[] = { + 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10}; + +int getusno (CatStats *catstats, RefCatalog *Ref) { + + long int offset; + int i, bin, first, last, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128]; + FILE *f; + double DEC1; + int iDEC0, iDEC1, iRA0, iRA1; + int spd, spd_start, spd_end, disk; + int NUSNO, Nusno; + StarData *stars; + + /* identify ra & dec range of interest */ + iRA0 = catstats[0].RA[0] * 360000.0; + iRA1 = catstats[0].RA[1] * 360000.0; + iDEC0 = (catstats[0].DEC[0] + 90.0) * 360000.0; + iDEC1 = (catstats[0].DEC[1] + 90.0) * 360000.0; + + /* data is organized in south-pole distance zones */ + spd_start = (int)((catstats[0].DEC[0] + 90) / 7.5) * 75.0; + DEC1 = (catstats[0].DEC[1] + 90) / 7.5; + if (DEC1 > (int)(DEC1)) { + spd_end = (int)(1 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } else { + spd_end = (int)(0 + (catstats[0].DEC[1] + 90) / 7.5) * 75.0; + } + + Nusno = 0; + NUSNO = 5000; + ALLOCATE (stars, StarData, NUSNO); + + for (spd = spd_start; spd < spd_end; spd += 75) { + disk = -1; + for (i = 0; i < NZONE; i++) { + if (spd == SPDzone[i]) + disk = USNOdisk[i]; + } + if (disk < 0) { + fprintf (stderr, "ERROR: can't find USNO zone for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc", USNO_A_DIR, spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open accelerator file %s\n", filename); + exit (1); + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; + else + last = 1 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", USNO_A_DIR, spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s\n", filename); + exit (1); + } + /* advance file pointer to first slice */ + offset = 3*sizeof(int)*(start[first] - 1); + fseek (f, offset, SEEK_SET); + /* on each loop, load data from an RA slice of the catalog */ + for (bin = first; bin < last; bin++) { + Nitems = 3*number[bin]; + ALLOCATE (buffer, int, Nitems); + nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failure reading data from file %s\n", filename); + exit (1); + } + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < number[bin]; i++, buf+=3) { + if ((buf[0] > iRA0) && (buf[0] < iRA1) && + (buf[1] > iDEC0) && (buf[1] < iDEC1)) { + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + /* note that this is the RED mag */ + stars[Nusno].M = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000)))); + /* b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); */ + Nusno ++; + if (Nusno == NUSNO) { + NUSNO += 5000; + REALLOCATE (stars, StarData, NUSNO); + } + } + } + free (buffer); + } + fclose (f); + } + + area_of_region (catstats); + + REALLOCATE (stars, StarData, MAX (1, Nusno)); + + Ref[0].stars = stars; + Ref[0].N = Nusno; + Ref[0].R0 = catstats[0].RA[0]; + Ref[0].R1 = catstats[0].RA[1]; + Ref[0].D0 = catstats[0].DEC[0]; + Ref[0].D1 = catstats[0].DEC[1]; + Ref[0].Area = catstats[0].Area; + + get_luminosity_func (Ref[0].stars, Ref[0].N, &Ref[0].lum); + + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno); + return (TRUE); +} + + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusnob.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusnob.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/getusnob.c (revision 21563) @@ -0,0 +1,163 @@ +# include "gastro2.h" +# define NBYTE 4 +# define NELEM 20 + +int getusnob (CatStats *catstats, RefCatalog *Ref, double epoch) { + + long int offset; + int i, bin, first, last, nitems, Nitems, Nbins; + float hours[100]; + int start[100], number[100], *buffer, *buf; + char filename[128]; + FILE *f; + double DEC1; + double uR, uD; + float mB1, mB2, mR1, mR2, mB, mR; + int iDEC0, iDEC1, iRA0, iRA1; + int spd, spd_start, spd_end; + int NUSNO, Nusno, Nstars; + StarData *stars; + + /* identify ra & dec range of interest */ + iRA0 = catstats[0].RA[0] * 360000.0; + iRA1 = catstats[0].RA[1] * 360000.0; + iDEC0 = (catstats[0].DEC[0] + 90.0) * 360000.0; + iDEC1 = (catstats[0].DEC[1] + 90.0) * 360000.0; + /* note that DEC is in SPD, while both have units to 0.01 degrees */ + + /* data is organized in south-pole distance zones, 1 deg per direction, 0.1 deg per file */ + spd_start = (int)(10*(catstats[0].DEC[0] + 90)); + DEC1 = 10*(catstats[0].DEC[1] + 90); + if (DEC1 > (int)(DEC1)) { + spd_end = (int)(1 + 10*(catstats[0].DEC[1] + 90)); + } else { + spd_end = (int)(0 + 10*(catstats[0].DEC[1] + 90)); + } + + Nusno = 0; + NUSNO = 5000; + ALLOCATE (stars, StarData, NUSNO); + + for (spd = spd_start; spd < spd_end; spd ++) { + + /* load accelerator file */ + sprintf (filename, "%s/%03d/b%04d.acc", USNO_B_DIR, (int)(spd/10), spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open accelerator file %s\n", filename); + exit (1); + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; + else + last = 1 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/%03d/b%04d.cat", USNO_B_DIR, (int)(spd/10), spd); + if (VERBOSE) fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s\n", filename); + exit (1); + } + + /** USNO-B consists of 20 x 4byte (int) records **/ + /* advance file pointer to first slice */ + offset = NELEM*NBYTE*(start[first] - 1); + fseek (f, offset, SEEK_SET); + + /* sum the number of stars in data segment of interest */ + Nstars = 0; + for (bin = first; bin < last; bin++) { + Nstars += number[bin]; + } + Nitems = NELEM*Nstars; /* number of integer blocks; need to use Fread for byte-swapping read */ + + /* allocate space for stars in segment */ + ALLOCATE (buffer, int, Nitems); + // data has the WRONG byte order? + // nitems = Fread (buffer, sizeof(int), Nitems, f, "int"); + nitems = fread (buffer, sizeof(int), Nitems, f); + if (nitems != Nitems) { + fprintf (stderr, "ERROR: failure reading data from file %s\n", filename); + exit (1); + } + + buf = buffer; + /* print out data from slice within RA and DEC range */ + for (i = 0; i < Nstars; i++, buf += NELEM) { + if (buf[0] < iRA0) continue; + if (buf[0] > iRA1) continue; + if (buf[1] < iDEC0) continue; + if (buf[1] > iDEC1) continue; + + bzero (&stars[Nusno], sizeof(StarData)); + stars[Nusno].R = buf[0]/360000.0; + stars[Nusno].D = buf[1]/360000.0 - 90.0; + + uR = (buf[2] % 10000); + uR = (uR - 5000.0) * 0.002 / 3600.0; + uD = ((buf[2] / 10000) % 10000); + uD = (uD - 5000.0) * 0.002 / 3600.0; + + /* 1st blue mag */ + mB1 = 0.01 * (buf[5] % 10000); + /* 1st blue mag */ + mB2 = 0.01 * (buf[6] % 10000); + /* 1st blue mag */ + mR1 = 0.01 * (buf[7] % 10000); + /* 1st blue mag */ + mR2 = 0.01 * (buf[8] % 10000); + + if (mB1 && mB2) { + mB = 0.5*(mB1 + mB2); + } else { + mB = (mB1) ? mB1 : mB2; + } + + if (mR1 && mR2) { + mR = 0.5*(mR1 + mR2); + } else { + mR = (mR1) ? mR1 : mR2; + } + + stars[Nusno].M = mB; + stars[Nusno].R += uR*(epoch - 2000.0); + stars[Nusno].D += uD*(epoch - 2000.0); + Nusno ++; + CHECK_REALLOCATE (stars, StarData, NUSNO, Nusno, 5000); + } + free (buffer); + fclose (f); + } + + area_of_region (catstats); + + REALLOCATE (stars, StarData, MAX (1, Nusno)); + + Ref[0].stars = stars; + Ref[0].N = Nusno; + Ref[0].R0 = catstats[0].RA[0]; + Ref[0].R1 = catstats[0].RA[1]; + Ref[0].D0 = catstats[0].DEC[0]; + Ref[0].D1 = catstats[0].DEC[1]; + Ref[0].Area = catstats[0].Area; + + get_luminosity_func (Ref[0].stars, Ref[0].N, &Ref[0].lum); + + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno); + return (TRUE); +} + + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gfit2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gfit2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gfit2.c (revision 21563) @@ -0,0 +1,83 @@ +# include "gastro2.h" + +void gfit (CmpCatalog *Target, RefCatalog *Ref, int order) { + + int i, j, j0; + int Npair, *idx1, *idx2; + double Radius, Radius2; + double dX, dY, dR; + RefCatalog Subset; + StarData *st, *sr; + + /* XXX why is this hardwired here? */ + NFIELD = 0.1; + gproject (Target, Ref, &Subset); + if (PLOTSTUFF) plot_fullfield (Target, &Subset); + + if (Subset.N < 3) { + fprintf (stderr, "ERROR: solution off target\n"); + exit (1); + } + + /* need the stars sorted in X */ + sort_stars_X (Target[0].stars, Target[0].N); + sort_stars_X (Subset.stars, Subset.N); + + Radius = MAX (2.5 * Target[0].answer.dR, 0.5); + Radius2 = Radius*Radius; + + st = Target[0].stars; + sr = Subset.stars; + + /* find the matched pairs of stars within the radius */ + pair_init (); + for (i = j = 0; (i < Target[0].N) && (j < Subset.N);) { + /* get in right X range */ + dX = st[i].X - sr[j].X; + if (dX < -Radius) { + i++; + continue; + } + if (dX > Radius) { + j++; + continue; + } + + /* check for pairs in this X range */ + j0 = j; + for (; (dX > -Radius) && (j < Subset.N); j++) { + + dX = st[i].X - sr[j].X; + dY = st[i].Y - sr[j].Y; + + dR = dX*dX + dY*dY; + if (dR > Radius2) { + j++; + continue; + } + pair_add (i, j); + } + j = j0; + i ++; + } + + Npair = pair_lists (&idx1, &idx2); + /* find fit for matched pairs */ + fit_init (order); + for (i = 0; i < Npair; i++) { + fit_add (st[idx1[i]].X, st[idx1[i]].Y, sr[idx2[i]].P, sr[idx2[i]].Q, 1.0); + } + fit_eval (); + + /* XXX this is weak: the fit_scat call requires the coords from the fit_adjust call */ + Target[0].answer.N = fit_adjust (&Target[0].coords); + Target[0].answer.dR = fit_scat (st, sr, &Target[0].coords); + + if (PLOTSTUFF) fprintf (stderr, "ploting resid (2)\n"); + plot_resid_init (0, (double) Target[0].header.Naxis[0]); + // XXX test: plotting fit_apply and RD_to_XY results plot_resid_init (1, (double) Target[0].header.Naxis[1]); + plot_resid_init (1, (double) Target[0].header.Naxis[0]); + if (PLOTSTUFF) plot_resid (st, sr, &Target[0].coords); + free (idx1); + free (idx2); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gheader2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gheader2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gheader2.c (revision 21563) @@ -0,0 +1,106 @@ +# include "gastro2.h" + +void gheader (char *file, CmpCatalog *Target) { + + double dR; + Header header; + FILE *f, *g; + int i, oldsize, nbytes, status; + char line[1024]; + + if (!gfits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (3)\n", file); + exit(0); + } + oldsize = header.size; + + /* check for insufficient number of stars */ + switch (Target[0].coords.Npolyterms) { + case 0: + case 1: + if (Target[0].answer.N < 6) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + case 2: + if (Target[0].answer.N < 12) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + case 3: + if (Target[0].answer.N < 20) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + break; + default: + fprintf (stderr, "invalid order\n"); + exit (2); + } + + gfits_modify (&header, "NASTRO", "%d", 1, Target[0].answer.N); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + + /*** use PutCoords to update header ***/ + PutCoords (&Target[0].coords, &header); + + dR = fabs (Target[0].answer.dR*Target[0].coords.cdelt1*3600.0); + gfits_modify (&header, "CERROR", "%lf", 1, dR); + gfits_modify (&header, "CERROR", "%C", 1, "scatter in astrometry soln (arcsec)"); + gfits_modify (&header, "CPRECISE", "%lf", 1, dR / sqrt(1.0*Target[0].answer.N)); + gfits_modify (&header, "CPRECISE", "%C", 1, "precision of astrometry soln (arcsec)"); + gfits_modify (&header, "EQUINOX", "%lf", 1, 2000.0); + /* we force equinox to be 2000.0 for all images */ + +skipstuff: + if (header.size > oldsize) { + if (VERBOSE) fprintf (stderr, "header expanded, creating new copy\n"); + sprintf (line, "mv %s %s~", file, file); + status = system (line); + if (status) { + fprintf (stderr, "ERROR: unable to create %s~, exiting\n", file); + exit (0); + } + sprintf (line, "%s~", file); + f = fopen (line, "r"); + g = fopen (file, "w"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find image file %s (4)\n", line); + exit(0); + } + if (g == NULL) { + fprintf (stderr, "ERROR: can't open output image file %s (4)\n", file); + exit(0); + } + nbytes = fwrite (header.buffer, 1, header.size, g); + fseek (f, oldsize, SEEK_SET); + for (i = 0; (nbytes = fread (header.buffer, 1, header.size, f)) > 0; i++) { + if (nbytes != fwrite (header.buffer, 1, nbytes, g)) { + fprintf (stderr, "ERROR: failure writing output data file\n"); + exit (0); + } + } + fclose (f); + fclose (g); + } else { + f = fopen (file, "r+"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find image file %s (4)\n", file); + exit(0); + } + + fseek (f, 0, SEEK_SET); + nbytes = fwrite (header.buffer, 1, header.size, f); + + fclose (f); + } + free (header.buffer); + +} + + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gpairs.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gpairs.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gpairs.c (revision 21563) @@ -0,0 +1,35 @@ +# include "gastro2.h" + +static int NPAIR, Npair; +static int *idx1, *idx2; + +void pair_init () { + + Npair = 0; + NPAIR = 100; + ALLOCATE (idx1, int, NPAIR); + ALLOCATE (idx2, int, NPAIR); + +} + +void pair_add (int i1, int i2) { + + idx1[Npair] = i1; + idx2[Npair] = i2; + + Npair ++; + + if (Npair == NPAIR) { + NPAIR += 100; + REALLOCATE (idx1, int, NPAIR); + REALLOCATE (idx2, int, NPAIR); + } + +} + +int pair_lists (int **index1, int **index2) { + + *index1 = idx1; + *index2 = idx2; + return (Npair); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gproject2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gproject2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gproject2.c (revision 21563) @@ -0,0 +1,134 @@ +# include "gastro2.h" + +void gproject (CmpCatalog *Target, RefCatalog *Ref, RefCatalog *Subset) { + + int i, N; + double X, Y, P, Q, M, Moff; + double XMIN, XMAX, YMIN, YMAX, MMIN, MMAX; + Coords TPtoSky, FPtoTP, *coords; + StarData *in, *out; + + Subset[0] = Ref[0]; + + ALLOCATE (Subset[0].stars, StarData, MAX (1, Ref[0].N)); + + in = Ref[0].stars; + out = Subset[0].stars; + + coords = &Target[0].coords; + + /* create Tangent Plane to Sky transformation from input coords */ + strcpy (TPtoSky.ctype, coords[0].ctype); + TPtoSky.crval1 = coords[0].crval1; + TPtoSky.crval2 = coords[0].crval2; + + TPtoSky.cdelt1 = TPtoSky.cdelt2 = 1; + TPtoSky.crpix1 = TPtoSky.crpix2 = 0; + TPtoSky.pc1_1 = TPtoSky.pc2_2 = 1; + TPtoSky.pc1_2 = TPtoSky.pc2_1 = 0; + TPtoSky.Npolyterms = 0; + for (i = 0; i < 7; i++) { + TPtoSky.polyterms[i][0] = 0; + TPtoSky.polyterms[i][1] = 0; + } + + /* create Focal Plane to Tangent Plane transformation from input coords */ + strcpy (FPtoTP.ctype, "FP---PLY"); + FPtoTP.crval1 = FPtoTP.crval2 = 0; + + FPtoTP.cdelt1 = coords[0].cdelt1; + FPtoTP.cdelt2 = coords[0].cdelt2; + FPtoTP.crpix1 = coords[0].crpix1; + FPtoTP.crpix2 = coords[0].crpix2; + FPtoTP.pc1_1 = coords[0].pc1_1; + FPtoTP.pc1_2 = coords[0].pc1_2; + FPtoTP.pc2_1 = coords[0].pc2_1; + FPtoTP.pc2_2 = coords[0].pc2_2; + + FPtoTP.Npolyterms = coords[0].Npolyterms; + for (i = 0; i < 7; i++) { + FPtoTP.polyterms[i][0] = coords[0].polyterms[i][0]; + FPtoTP.polyterms[i][1] = coords[0].polyterms[i][1]; + } + + Moff = Ref[0].Moff; + + XMIN = -0.5*NFIELD*Target[0].header.Naxis[0]; + XMAX = 0.5*NFIELD*Target[0].header.Naxis[0] + Target[0].header.Naxis[0]; + YMIN = -0.5*NFIELD*Target[0].header.Naxis[1]; + YMAX = 0.5*NFIELD*Target[0].header.Naxis[1] + Target[0].header.Naxis[1]; + + /* need to allow some leeway? use a fixed +/- 0.5 mag for now */ + if (MAGLIMS && !MAGMANUAL) { + MMAX = Target[0].lum.Mmax + 1.0; + MMIN = Target[0].lum.Mmin - 1.0; + + if (MMAX < Ref[0].lum.Mmin + Moff) + fprintf (stderr, "warning: reference catalog probably too faint: %5.3f < %5.3f\n", Target[0].lum.Mmax, Ref[0].lum.Mmin + Moff); + + if (MMIN > Ref[0].lum.Mmax + Moff) + fprintf (stderr, "warning: reference catalog probably too bright: %5.3f > %5.3f\n", Target[0].lum.Mmin, Ref[0].lum.Mmax + Moff); + } + if (MAGMANUAL) { + MMIN = MAGLIM_MIN; + MMAX = MAGLIM_MAX; + Moff = 0; + } + + if (VERBOSE) fprintf (stderr, "limited reference stars to mag range %f - %f\n", MMIN - Moff, MMAX - Moff); + + for (N = i = 0; i < Ref[0].N; i++) { + RD_to_XY (&X, &Y, in[i].R, in[i].D, coords); + M = in[i].M + Moff; + + if (X < XMIN) continue; + if (X > XMAX) continue; + if (Y < YMIN) continue; + if (Y > YMAX) continue; + + if (MAGLIMS) { + if (M < MMIN) continue; + if (M > MMAX) continue; + } + + /* get tangent-plane coordinates as well */ + RD_to_XY (&P, &Q, in[i].R, in[i].D, &TPtoSky); + + out[N] = in[i]; + out[N].X = X; + out[N].Y = Y; + out[N].P = P; + out[N].Q = Q; + out[N].M = M; + N++; + } + + if (N < 3) { + fprintf (stderr, "ERROR: too few reference stars accepted\n"); + exit (1); + } + + Subset[0].N = N; + sort_stars_mag (Subset[0].stars, N); + if (GASTRO_MAX_NSTARS && (GASTRO_MAX_NSTARS < Subset[0].N)) { + Subset[0].N = GASTRO_MAX_NSTARS; + REALLOCATE (Subset[0].stars, StarData, Subset[0].N); + } + if (VERBOSE) fprintf (stderr, "using %d stars from ref catalog\n", Subset[0].N); + + REALLOCATE (Subset[0].stars, StarData, MAX (1, Subset[0].N)); + + /* get tangent-plane coords for target as well */ + for (i = 0; i < Target[0].N; i++) { + XY_to_RD (&P, &Q, Target[0].stars[i].X, Target[0].stars[i].Y, &FPtoTP); + Target[0].stars[i].P = P; + Target[0].stars[i].Q = Q; + } + +} + +/* in this function, we convert the Ra & Dec coords to the rough X, Y coords + we also convert the magnitudes to the approximate system with Ref[0].Moff + also, reduce the domain to those within X, Y, M limits +*/ + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/greference2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/greference2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/greference2.c (revision 21563) @@ -0,0 +1,131 @@ +# include "gastro2.h" + +void greference (CmpCatalog *Target, RefCatalog *Ref) { + + CatStats catstats; + + if (VERBOSE) fprintf (stderr, "loading astrometric reference data from %s\n", REFCAT); + + define_region (&catstats, Target); + + Ref[0].N = 0; + /* get stars from the USNO A catalog for the given region */ + if (!strcasecmp (REFCAT, "USNO")) { + getusno (&catstats, Ref); + /* calculate Ref[0].Moff from Target & Ref dMdN, Mo */ + } + + /* get stars from the USNO B catalog for the given region */ + if (!strcasecmp (REFCAT, "USNOB")) { + getusnob (&catstats, Ref, 2000.0); + /* calculate Ref[0].Moff from Target & Ref dMdN, Mo */ + } + + /* get stars from the HST GSC catalog for the given region */ + if (!strcasecmp (REFCAT, "GSC")) { + getgsc (&catstats, Ref); + } + + /* get stars from 2MASS for the given region -- add PHOTCODE check? */ + if (!strcasecmp (REFCAT, "2MASS")) { + strcpy (CATDIR, TWO_MASS_DIR); + getptolemy (&catstats, Ref); + } + + /* get stars from the DVO CATDIR for the given region */ + if (!strcasecmp (REFCAT, "PTOLEMY")) { + strcpy (CATDIR, ASTROM_CATDIR); + getptolemy (&catstats, Ref); + } + + if (Ref[0].N == 0) { + fprintf (stderr, "no ref objs: %s\n", REFCAT); + exit (1); + } + + { + double Mref, Mtar, logRho; + + /* what is the offset between the two lines at the average magnitude? */ + Mref = 0.5*(Ref[0].lum.Mmin + Ref[0].lum.Mmax); + logRho = Mref * Ref[0].lum.dNdM + Ref[0].lum.Mo - log10(Ref[0].Area); + Mtar = (logRho + log10(Target[0].Area) - Target[0].lum.Mo) / Target[0].lum.dNdM; + + Ref[0].lum.Mz = Mref; + Target[0].lum.Mz = Mtar; + Ref[0].Moff = Target[0].lum.Mz - Ref[0].lum.Mz; + fprintf (stderr, "mag offset: %f (Areas: %f vs %f; log(Rho): %f @ %f mags)\n", Ref[0].Moff, Target[0].Area, Ref[0].Area, logRho, Mref); + } + + if (PLOTSTUFF) plot_lumfunc (Target, Ref); + +} + +/* return RA, DEC bounds of the reigon of interest */ +void define_region (CatStats *catstats, CmpCatalog *Target) { + + int NX, NY, status; + double x, y, X, Y, R, D, dX, dY, Xo, Yo; + + NX = Target[0].header.Naxis[0]; + NY = Target[0].header.Naxis[1]; + + dX = NX + NFIELD*NX; + dY = NY + NFIELD*NY; + + Xo = -0.5*NFIELD*NX; + Yo = -0.5*NFIELD*NY; + + if (ASCA) { + XY_to_RD (&R, &D, 0.5*NX, 0.5*NY, &Target[0].coords); + catstats[0].RA[0] = R - 90.0; + catstats[0].RA[1] = R + 90.0; + catstats[0].DEC[0] = MAX (-90.0, D - 90.0); + catstats[0].DEC[1] = MIN (+90.0, D + 90.0); + if (VERBOSE) fprintf (stderr, "asca region: %f - %f, %f - %f\n", + catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]); + return; + } + + catstats[0].RA[0] = 360.0; + catstats[0].RA[1] = 0.0; + catstats[0].DEC[0] = +90.0; + catstats[0].DEC[1] = -90.0; + + for (x = 0; x <= 1.0; x += 0.5) { + for (y = 0; y <= 1.0; y += 0.5) { + + X = x*dX + Xo; + Y = y*dY + Yo; + status = XY_to_RD (&R, &D, X, Y, &Target[0].coords); + if (!status) continue; + if (isinf(R) || isnan(R)) continue; + if (isinf(D) || isnan(D)) continue; + + catstats[0].RA[0] = MIN (catstats[0].RA[0], R); + catstats[0].RA[1] = MAX (catstats[0].RA[1], R); + catstats[0].DEC[0] = MIN (catstats[0].DEC[0], D); + catstats[0].DEC[1] = MAX (catstats[0].DEC[1], D); + } + } + + /* is a pole in the image? if so, include it... */ + status = RD_to_XY (&X, &Y, 0.0, 90.0, &Target[0].coords); + if (status) { + if (fabs(X - NX*0.5) > (NFIELD + 0.5)*NX) goto not_north; + if (fabs(Y - NY*0.5) > (NFIELD + 0.5)*NY) goto not_north; + catstats[0].DEC[1] = 90.0; + } +not_north: + + status = RD_to_XY (&X, &Y, 0.0, -90.0, &Target[0].coords); + if (status) { + if (fabs(X - NX*0.5) > (NFIELD + 0.5)*NX) goto not_south; + if (fabs(Y - NY*0.5) > (NFIELD + 0.5)*NY) goto not_south; + catstats[0].DEC[0] = -90.0; + } +not_south: + + if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", + catstats[0].RA[0], catstats[0].RA[1], catstats[0].DEC[0], catstats[0].DEC[1]); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/grid.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/grid.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/grid.c (revision 21563) @@ -0,0 +1,219 @@ +# include "gastro2.h" + +/* this needs to be user-configured */ +static int NX, NY, Nbin; +static double C0x, C1x, C0y, C1y; +static double *N, *DX, *DY, *D2; + +void grid (CmpCatalog *Target, RefCatalog *Subset, Answer *answer) { + + int i, j, n, Imin, Nmin, Nval; + double XMIN, XMAX, YMIN, YMAX; + double dX, dY; + double s, f, Fmin, Smin, *ntmp; + + RefCatalog Ref; + StarData *st, *sr; + + /* Ref is temporary in this function. Free Ref.stars before exiting */ + rotate (Subset, &Ref, answer[0].angle); + + st = Target[0].stars; + sr = Ref.stars; + + /* NFIELD represents the search box; it is also the extra padding used + to select the reference stars */ + XMIN = -0.5*NFIELD*Target[0].header.Naxis[0]; + XMAX = +0.5*NFIELD*Target[0].header.Naxis[0]; + YMIN = -0.5*NFIELD*Target[0].header.Naxis[1]; + YMAX = +0.5*NFIELD*Target[0].header.Naxis[1]; + + dX = dY = 0; + /* make two passes, with grids offset by 0.5 box for the second */ + for (n = 0; n < 2; n++) { + gridinit (XMIN, XMAX, YMIN, YMAX, Ref.N, Target[0].N); + + if (PLOTSTUFF) plot_init_gridplot (); + /* fill in grid points */ + for (i = 0; i < Target[0].N; i++) { + for (j = 0; j < Ref.N; j++) { + + dX = st[i].X - sr[j].X; + if (dX < XMIN) continue; + if (dX > XMAX) continue; + + dY = st[i].Y - sr[j].Y; + if (dY < YMIN) continue; + if (dY > YMAX) continue; + + gridbin (dX, dY); + if (PLOTSTUFF) plot_addpt_gridplot (dX, dY); + } + } + + /* use sorted N list to define Nmin cut */ + ALLOCATE (ntmp, double, Nbin); + bcopy (N, ntmp, Nbin*sizeof(double)); + sort (ntmp, Nbin); + for (i = 0; (ntmp[i] == 0) && (i < Nbin); i++); + Nval = MIN (Nbin - 1, (int)(0.75*(Nbin - i)) + i); + Nmin = ntmp[Nval]; + free (ntmp); + + /* select 'best' grid point - is this statistic good enough? */ + Fmin = 1e10; + Imin = -1; + for (i = 0; i < Nbin; i++) { + + if (N[i] < Nmin) continue; + + /* s is the variance, f is varience overweighted by number */ + s = fabs ((D2[i]/N[i]) - SQ(DX[i]/N[i]) - SQ(DY[i]/N[i])); + f = s / SQ(SQ(N[i])); + + if (f < Fmin) { + Smin = s; + Fmin = f; + Imin = i; + } + } + if (Imin == -1) { + fprintf (stderr, "ERROR: odd min value\n"); + exit (1); + } + + if ((n == 0) || (Fmin < answer[0].Chi)) { + Smin = fabs ((D2[Imin]/N[Imin]) - SQ(DX[Imin]/N[Imin]) - SQ(DY[Imin]/N[Imin])); + answer[0].Xoff = DX[Imin] / N[Imin]; + answer[0].Yoff = DY[Imin] / N[Imin]; + answer[0].dR = sqrt (Smin); + answer[0].Chi = Fmin; + answer[0].N = N[Imin]; + } + + fprintf (stderr, "angle: %6.1f, (%6.1f,%6.1f) - %6.2f : %10.8f for %d pairs\n", + answer[0].angle, answer[0].Xoff, answer[0].Yoff, answer[0].dR, answer[0].Chi, answer[0].N); + + if (PLOTSTUFF) plot_gridpts (N, Nbin); + if (PLOTSTUFF) plot_done_gridplot (dX, dY); + + XMIN -= 0.5*NGRID_PIX; + XMAX -= 0.5*NGRID_PIX; + YMIN -= 0.5*NGRID_PIX; + YMAX -= 0.5*NGRID_PIX; + gridfree (); + } + + free (Ref.stars); + +} + +void gridinit (double XMIN, double XMAX, double YMIN, double YMAX, int Nr, int Nt) { + + NX = (XMAX - XMIN) / NGRID_PIX; + NY = (YMAX - YMIN) / NGRID_PIX; + + C1x = NX / (XMAX - XMIN); + C0x = - XMIN * NX / (XMAX - XMIN); + + C1y = NY / (YMAX - YMIN); + C0y = - YMIN * NY / (YMAX - YMIN); + + Nbin = NX*NY; + + ALLOCATE (N, double, Nbin); + ALLOCATE (DX, double, Nbin); + ALLOCATE (DY, double, Nbin); + ALLOCATE (D2, double, Nbin); + + bzero (N, Nbin*sizeof(double)); + bzero (DX, Nbin*sizeof(double)); + bzero (DY, Nbin*sizeof(double)); + bzero (D2, Nbin*sizeof(double)); + +} + +int gridbin (double dX, double dY) { + + int bin, xbin, ybin; + + xbin = (int) (C0x + dX * C1x); + ybin = (int) (C0y + dY * C1y); + + bin = xbin + NX * ybin; + + if (bin < 0) return (0); + if (bin >= Nbin) return (0); + + N[bin] += 1.0; + DX[bin] += dX; + DY[bin] += dY; + D2[bin] += dX*dX + dY*dY; + + return (bin); + +} + +void gridfree () { + + free (N); + free (DX); + free (DY); + free (D2); + +} + +# if (0) + /* use sorted N list to define Nmin cut */ + ALLOCATE (ntmp, double, Nbin); + bcopy (N, ntmp, Nbin*sizeof(double)); + sort (ntmp, Nbin); + for (i = 0; (ntmp[i] == 0) && (i < Nbin); i++); + Nval = MIN (Nbin - 1, (int)(0.75*(Nbin - i)) + i); + Nmin = ntmp[Nval]; + free (ntmp); + + /* select 'best' grid point - is this statistic good enough? */ + Fmin = 1e10; + Imin = -1; + for (i = 0; i < Nbin; i++) { + + if (N[i] < Nmin) continue; + + /* s is the variance, f is varience overweighted by number */ + s = fabs ((D2[i]/N[i]) - SQ(DX[i]/N[i]) - SQ(DY[i]/N[i])); + f = s / SQ(SQ(N[i])); + + if (f < Fmin) { + Smin = s; + Fmin = f; + Imin = i; + } + } + if (Imin == -1) { + fprintf (stderr, "ERROR: odd min value\n"); + exit (1); + } +# endif + +# if (0) + /* find N sigma */ + Ns = Ns2 = 0; + for (i = 0; i < Nbin; i++) { + Ns += N[i]; + Ns2 += N[i]*N[i]; + } + Ns = Ns / Nbin; + Ns2 = sqrt (Ns2 / Nbin - Ns*Ns); + fprintf (stderr, "N sigma: %f\n", Ns2); + + Imin = 0; + Fmin = N[0] / Ns2; + for (i = 0; i < Nbin; i++) { + f = N[i] / Ns2; + if (f > Fmin) { + Fmin = f; + Imin = i; + } + } +# endif Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gstars2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gstars2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/gstars2.c (revision 21563) @@ -0,0 +1,177 @@ +# include "gastro2.h" + +void gstars (char *filename, CmpCatalog *Target) { + + int Nstars; + char line[80]; + double det; + int NX, NY, Nskip, FoundAstrom, extend, naxis; + StarData *stars; + FILE *f; + + /* open file for stars */ + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't open file to load stars\n"); + exit (1); + } + if (!gfits_fread_header (f, &Target[0].header)) { + fprintf (stderr, "ERROR: can't read image header\n"); + exit (1); + } + /* this line should not be needed */ + fseek (f, Target[0].header.size, SEEK_SET); + + NX = Target[0].header.Naxis[0]; + NY = Target[0].header.Naxis[1]; + + /* default values for coords */ + strcpy (Target[0].coords.ctype, "RA---TAN"); + Target[0].coords.pc1_1 = CCD_PC1_1; Target[0].coords.pc1_2 = CCD_PC1_2; + Target[0].coords.pc2_1 = CCD_PC2_1; Target[0].coords.pc2_2 = CCD_PC2_2; + Target[0].coords.cdelt1 = Target[0].coords.cdelt2 = ASEC_PIX / 3600.0; + Target[0].coords.Npolyterms = 0; + Target[0].coords.crpix1 = 0.5*NX; + Target[0].coords.crpix2 = 0.5*NY; + + /* attempt to get detailed astrometric information from header */ + FoundAstrom = FALSE; + if (!strcasecmp (ROUGH_ASTROMETRY, "header")) { + if (!HEADER[0]) { + FoundAstrom = GetCoords (&Target[0].coords, &Target[0].header); + } else { + Header header; + + if (!gfits_read_header (HEADER, &header)) { + fprintf (stderr, "ERROR: can't load external header\n"); + exit (1); + } + FoundAstrom = GetCoords (&Target[0].coords, &header); + gfits_free_header (&header); + } + if (!FoundAstrom) { + fprintf (stderr, "header coordinates incomplete, trying for rough coordinates\n"); + strcpy (ROUGH_ASTROMETRY, "config"); + } else { + /* make optional adustments to header values (standard problems) */ + if (FLIPX) { + Target[0].coords.pc1_1 *= -1; + Target[0].coords.crpix1 = NX - Target[0].coords.crpix1; + } + if (FLIPY) { + Target[0].coords.pc2_2 *= -1; + Target[0].coords.crpix2 = NY - Target[0].coords.crpix2; + } + } + } + + /*** abstract RA & DEC keywords, formats */ + /* get just RA & DEC from header, other terms from config file */ + if (!strcasecmp (ROUGH_ASTROMETRY, "config")) { + /* get RA & DEC from header */ + FoundAstrom = TRUE; + FoundAstrom &= gfits_scan (&Target[0].header, "RA", "%s", 1, line); + dms_to_ddd (&Target[0].coords.crval1, line); + Target[0].coords.crval1 = Target[0].coords.crval1 * 15.0; + FoundAstrom &= gfits_scan (&Target[0].header, "DEC", "%s", 1, line); + dms_to_ddd (&Target[0].coords.crval2, line); + } + + /* use RA & DEC from command line arguments */ + if (FORCE) { + Target[0].coords.crval1 = F_RA; + Target[0].coords.crval2 = F_DEC; + if (VERBOSE) fprintf (stderr, " forcing coordinates to: %9.4f %9.4f\n", Target[0].coords.crval1, Target[0].coords.crval2); + FoundAstrom = TRUE; + } + + if (VERBOSE) fprintf (stderr, "using coordinates: %9.4f %9.4f\n", Target[0].coords.crval1, Target[0].coords.crval2); + if (!FoundAstrom) { + fprintf (stderr, "ERROR: can't get any valid coordinates, fix config file?\n"); + exit (1); + } + + /* at this point, we need to correct the crval1, crval2, and ROT_ZERO values + based on the pole axis angle and the ra, dec offsets */ + +# define dcos(a) (cos((double)((a)*(RAD_DEG)))) +# define dsin(a) (sin((double)((a)*(RAD_DEG)))) + + if (POLAR_ALIGNMENT) { + + double X, Y, PD, PR, DE, RE, T1, T2, T3; + + X = Target[0].coords.crval1; + Y = Target[0].coords.crval2; + PD = POLE_DEC; PR = POLE_RA; + DE = DEC_OFFSET; RE = RA_OFFSET; + + T1 = dcos(Y-DE) * dcos(X-RE) * dsin(PD) + dsin(Y-DE) * dcos(PD); + T2 = dcos(Y-DE) * dsin(X-RE); + T3 = dsin(Y-DE) * dsin(PD) - dcos(Y-DE) * dcos(X-RE) * dcos(PD); + + Target[0].coords.crval1 = (DEG_RAD * atan2 (T2, T1)) + PR; + Target[0].coords.crval2 = (DEG_RAD * asin (T3)); + while (Target[0].coords.crval1 < 0) Target[0].coords.crval1 += 360.0; + while (Target[0].coords.crval1 > 360.0) Target[0].coords.crval1 -= 360.0; + + if (VERBOSE) fprintf (stderr, " after polar alignment: %9.4f %9.4f\n", Target[0].coords.crval1, Target[0].coords.crval2); + } + + /* get image area in deg^2 */ + det = Target[0].coords.pc1_1 * Target[0].coords.pc2_2 - Target[0].coords.pc1_2 * Target[0].coords.pc2_1; + Target[0].Area = fabs (NX*NY*Target[0].coords.cdelt1*Target[0].coords.cdelt2*det); + + /* read from FITS table or from text table */ + /* Is NAXIS == 0 a better test?? */ + extend = FALSE; + gfits_scan (&Target[0].header, "NAXIS", "%t", 1, &naxis); + if ((naxis == 0) && !TEXTMODE) { + Nskip = gfits_matrix_size (&Target[0].header); + fseek (f, Nskip, SEEK_CUR); + stars = rfits (f, &Nstars); + } else { + /* allocate space for stars */ + if (!gfits_scan (&Target[0].header, "NSTARS", "%d", 1, &Nstars)) { + fprintf (stderr, "ERROR: failed to find NSTARS\n"); + exit (1); + } + stars = rtext (f, &Nstars); + } + fclose (f); + + stars = remove_clumps (stars, &Nstars, NX, NY); + + sort_stars_mag (stars, Nstars); /* sorting by magnitude */ + Target[0].stars = stars; + Target[0].N = Nstars; + + /* limit number of stars */ + if (GASTRO_MAX_NSTARS && (GASTRO_MAX_NSTARS < Target[0].N)) { + Target[0].N = GASTRO_MAX_NSTARS; + REALLOCATE (Target[0].stars, StarData, Target[0].N); + } + if (VERBOSE) fprintf (stderr, "using %d stars from data file\n", Target[0].N); + + /* calculate luminosity function of stars */ + get_luminosity_func (stars, Target[0].N, &Target[0].lum); + +} + +/* + +load cmp file FITS header + +extract needed data from header: +- NX, NY? +- coords + +load stellar photometry + +sort stars by mag + +filter & limit numbers + +find luminosity function slope, area? + +*/ Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/lumfunc.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/lumfunc.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/lumfunc.c (revision 21563) @@ -0,0 +1,120 @@ +# include "gastro2.h" + +/* mag range -5 - 35, dmag = 0.25, Nbin = 160 */ +# define MMIN -5 +# define MMAX 35 +# define dM 0.5 +# define NMBIN 90 + +int get_luminosity_func (StarData *stars, int N, LumStats *lum) { + + int i, j, Nb, peaki, peakn; + double mbin[NMBIN]; + double bin[NMBIN], lbin[NMBIN], rbin[NMBIN]; + double C0, C1; + + lum[0].dNdM = lum[0].Mo = 0; + + bzero (mbin, NMBIN * sizeof (double)); + + /* sum histogram */ + for (i = 0; i < N; i++) { + if (stars[i].M < MMIN) continue; + if (stars[i].M > MMAX) continue; + + j = (stars[i].M - MMIN) / dM; + j = MIN (MAX (j, 0), (NMBIN - 1)); + mbin[j] ++; + } + + /* find peak bin */ + peaki = 0; + peakn = mbin[0]; + for (i = 0; i < NMBIN; i++) { + if (mbin[i] > peakn) { + peaki = i; + peakn = mbin[i]; + } + } + + /* select filled bins */ + for (Nb = i = 0; i < peaki; i++) { + if (mbin[i] > 0) { + bin[Nb] = i * dM + MMIN; + lbin[Nb] = log10 (mbin[i]); + Nb++; + } + } + + /* find max & min mag bins */ + lum[0].Mmin = bin[0]; + lum[0].Mmax = bin[Nb-1]; + + if (Nb < 4) { return (FALSE); } + fit_lum_bin (bin, lbin, Nb, &C0, &C1); + + /* find residuals */ + for (i = 0; i < Nb; i++) { + rbin[i] = C0 + C1*bin[i] - lbin[i]; + } + + /* keep inner 80% */ + sort_lum (rbin, bin, lbin, Nb); + for (j = 0, i = 0.1*Nb; i < 0.9*Nb; i++, j++) { + bin[j] = bin[i]; + lbin[j] = lbin[i]; + } + Nb = j; + + if (Nb < 4) { return (FALSE); } + fit_lum_bin (bin, lbin, Nb, &C0, &C1); + + lum[0].dNdM = C1; + lum[0].Mo = C0; + + if (VERBOSE) fprintf (stderr, "lum stats: dNdM = %f, Mo = %f, Mmin = %f, Mmax = %f\n", + lum[0].dNdM, lum[0].Mo, lum[0].Mmin, lum[0].Mmax); + + return (TRUE); + +} + + +void fit_lum_bin (double *x, double *y, int N, double *C0, double *C1) { + + int i; + double **c, **b; + + ALLOCATE (c, double *, 2); + ALLOCATE (b, double *, 2); + ALLOCATE (c[0], double, 2); + ALLOCATE (c[1], double, 2); + ALLOCATE (b[0], double, 1); + ALLOCATE (b[1], double, 1); + + /* fit x, y to line */ + c[0][0] = 0; c[0][1] = 0; + c[1][0] = 0; c[1][1] = 0; + b[0][0] = 0; b[1][0] = 0; + + for (i = 0; i < N; i++) { + c[0][0] += 1; + c[0][1] += x[i]; + c[1][0] += x[i]; + c[1][1] += x[i]*x[i]; + + b[0][0] += y[i]; + b[1][0] += y[i]*x[i]; + } + dgaussj (c, 2, b, 1); + *C0 = b[0][0]; + *C1 = b[1][0]; + + free (b[0]); + free (b[1]); + free (c[0]); + free (c[1]); + free (c); + free (b); + +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/misc.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/misc.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/misc.c (revision 21563) @@ -0,0 +1,63 @@ +# include "gastro2.h" + +# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) + +void hh_hms (double hh, int *hr, int *mn, double *sc) { + + int flag; + + flag = SIGN(hh); + hh *= flag; + hh = 24.0*(hh/24.0 - (int)(hh/24.0)); + *sc = 60.0*(60.0*hh - (int)(60.0*hh)); + *mn = 60.0*(hh - (int)hh); + *hr = (int) hh; + *hr *= flag; + +} + +void hms_format (char *line, double value) { + + int hr, mn; + double sc; + + hh_hms (value, &hr, &mn, &sc); + hr = (int) value; + if (isnan (value)) + sprintf (line, "xx:xx:xx.xx"); + else { + if (value < 0) { + sprintf (line, "-%02d:%02d:%05.2f", abs(hr), mn, sc); + } else { + sprintf (line, "%02d:%02d:%05.2f", hr, mn, sc); + } + } +} + +void area_of_region (CatStats *region) { + + double area; + + area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG)); + region[0].Area = area; +} + +double area_of_skyregion (SkyRegion *region) { + + double area; + + area = DEG_RAD*(region[0].Rmax - region[0].Rmin)*(sin(region[0].Dmax*RAD_DEG) - sin(region[0].Dmin*RAD_DEG)); + return (area); +} + +void dump_coords (CmpCatalog *Target) { + if (VERBOSE) { + fprintf (stderr, "%s\n", Target[0].coords.ctype); + fprintf (stderr, "%f %f\n", Target[0].coords.crval1, Target[0].coords.crval2); + fprintf (stderr, "%f %f\n", Target[0].coords.crpix1, Target[0].coords.crpix2); + fprintf (stderr, "%f %f\n", Target[0].coords.pc1_1, Target[0].coords.pc1_2); + fprintf (stderr, "%f %f\n", Target[0].coords.pc2_1, Target[0].coords.pc2_2); + fprintf (stderr, "%f %f\n", Target[0].coords.cdelt1, Target[0].coords.cdelt2); + fprintf (stderr, "\n"); + } +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plots.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plots.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plots.c (revision 21563) @@ -0,0 +1,440 @@ +# include "gastro2.h" + +int Nv, NV; +float *xv, *yv; +Graphdata gv; + +void plot_init_gridplot () { + + NV = 1000; + Nv = 0; + + ALLOCATE (xv, float, NV); + ALLOCATE (yv, float, NV); + + PlotReset (2); + + gv.xmin = -400; + gv.xmax = +400; + gv.ymin = -400; + gv.ymax = +400; + + gv.style = 2; + gv.ltype = 0; + gv.etype = 0; + gv.ebar = 0; + gv.size = 0.3; + gv.lweight = 0; + gv.color = 0; + gv.ptype = 2; + +} + +int plot_addpt_gridplot (double x, double y) { + + if (x < gv.xmin) return (0); + if (x > gv.xmax) return (0); + if (y < gv.ymin) return (0); + if (y > gv.ymax) return (0); + + xv[Nv] = x; + yv[Nv] = y; + Nv ++; + + if (Nv == NV) { + NV += 1000; + REALLOCATE (xv, float, NV); + REALLOCATE (yv, float, NV); + } + return (1); +} + +void plot_done_gridplot () { + + char c; + + /* send to Kapa */ + PrepPlotting (Nv, &gv, 2); + PlotVector (Nv, xv, 0, 2); + PlotVector (Nv, yv, 1, 2); + + DonePlotting (&gv, 2); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); + free (xv); + free (yv); +} + +static double Xm[2]; +static int Xg[2] = {0, 2}; + +void plot_resid_init (int version, double xmax) { + + if (version > 1) return; + if (version < 0) return; + + Xm[version] = xmax; +} + + +void plot_resid_plot (int version, float *xvect, float *yvect, int Nvect) { + + int i; + char c; + float xmin, xmax, ymin, ymax; + Graphdata graphdata; + + xmin = xmax = xvect[0]; + ymin = ymax = yvect[0]; + for (i = 0; i < Nvect; i++) { + xmax = MAX (xvect[i], xmax); + xmin = MIN (xvect[i], xmin); + ymax = MAX (yvect[i], ymax); + ymin = MIN (yvect[i], ymin); + } + + PlotReset (Xg[version]); + + graphdata.xmin = MIN (xmin, 0); + graphdata.xmax = MAX (xmax, Xm[version]); + graphdata.ymin = MIN (ymin, -10); + graphdata.ymax = MAX (ymax, +10); + + graphdata.style = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.ebar = 0; + graphdata.size = 0.5; + graphdata.lweight = 0; + graphdata.color = 0; + graphdata.ptype = 2; + + /* send to Kapa */ + PrepPlotting (Nvect, &graphdata, Xg[version]); + PlotVector (Nvect, xvect, 0, Xg[version]); + PlotVector (Nvect, yvect, 1, Xg[version]); + + DonePlotting (&graphdata, Xg[version]); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); +} + +void plot_gridpts (double *pts, int Npts) { + + char c; + int i; + float *xvect, *yvect, ymax; + Graphdata graphdata; + + ymax = 0; + ALLOCATE (xvect, float, Npts); + ALLOCATE (yvect, float, Npts); + for (i = 0; i < Npts; i++) { + xvect[i] = i; + yvect[i] = pts[i]; + ymax = MAX (yvect[i], ymax); + } + + PlotReset (0); + + graphdata.xmin = -10; + graphdata.xmax = Npts + 10; + graphdata.ymin = -1; + graphdata.ymax = ymax + 1; + + graphdata.style = 1; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.ebar = 0; + graphdata.size = 0.5; + graphdata.lweight = 0; + graphdata.color = 0; + graphdata.ptype = 0; + + /* send to Kapa */ + PrepPlotting (Npts, &graphdata, 0); + PlotVector (Npts, xvect, 0, 0); + PlotVector (Npts, yvect, 1, 0); + free (xvect); + free (yvect); + + DonePlotting (&graphdata, 0); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); +} + +static Graphdata gf; + +void plot_fullfield (CmpCatalog *Target, RefCatalog *Ref) { + + int i, Nvect; + float *xvect, *yvect, *zvect, M, dM, mRefMin, mRefMax; + char c; + + dM = Target[0].lum.Mmax - Target[0].lum.Mmin; + + PlotReset (1); + + gf.xmin = -100; + gf.xmax = Target[0].header.Naxis[0] + 100; + gf.ymin = -100; + gf.ymax = Target[0].header.Naxis[1] + 100; + + gf.style = 2; + gf.ltype = 0; + gf.etype = 0; + gf.ebar = 0; + gf.size = -1; + gf.lweight = 0; + + /* fill in vectors */ + Nvect = Target[0].N; + ALLOCATE (xvect, float, Nvect); + ALLOCATE (yvect, float, Nvect); + ALLOCATE (zvect, float, Nvect); + for (i = 0; i < Nvect; i++) { + xvect[i] = Target[0].stars[i].X; + yvect[i] = Target[0].stars[i].Y; + M = (Target[0].lum.Mmax - Target[0].stars[i].M) / dM; + zvect[i] = MIN (1.0, MAX (0.01, M)); + } + + /* send to Kapa */ + gf.ptype = 1; + gf.color = 0; + PrepPlotting (Nvect, &gf, 1); + PlotVector (Nvect, xvect, 0, 1); + PlotVector (Nvect, yvect, 1, 1); + PlotVector (Nvect, zvect, 1, 1); + free (xvect); + free (yvect); + free (zvect); + + /* fill in vectors */ + Nvect = Ref[0].N; + ALLOCATE (xvect, float, Nvect); + ALLOCATE (yvect, float, Nvect); + ALLOCATE (zvect, float, Nvect); + + mRefMin = 30.0; + mRefMax = -5.0; + for (i = 0; i < Nvect; i++) { + mRefMin = MIN (mRefMin, Ref[0].stars[i].M); + mRefMax = MAX (mRefMax, Ref[0].stars[i].M); + } + dM = mRefMax - mRefMin; + + for (i = 0; i < Nvect; i++) { + xvect[i] = Ref[0].stars[i].X; + yvect[i] = Ref[0].stars[i].Y; + M = (mRefMax - Ref[0].stars[i].M) / dM; + zvect[i] = MIN (1.0, MAX (0.01, M)); + } + + /* send to Kapa */ + gf.color = 2; + gf.ptype = 2; + PrepPlotting (Nvect, &gf, 1); + PlotVector (Nvect, xvect, 0, 1); + PlotVector (Nvect, yvect, 1, 1); + PlotVector (Nvect, zvect, 2, 1); + free (xvect); + free (yvect); + free (zvect); + + DonePlotting (&gf, 1); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); +} + +void plot_fullfield_pairs (float *x, float *y, int n) { + + char c; + + /* send to Kapa */ + gf.color = 6; + gf.ptype = 100; + gf.size = 1.0; + + PrepPlotting (n, &gf, 1); + PlotVector (n, x, 0, 1); + PlotVector (n, y, 1, 1); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); +} + +/* mag range -5 - 35, dmag = 0.25, Nbin = 160 */ +# define MMIN -5 +# define MMAX 35 +# define dM 0.5 +# define NMBIN 90 + +void plot_lumfunc (CmpCatalog *Target, RefCatalog *Ref) { + + int i, Nr, Nt; + float tbin[NMBIN], tval[NMBIN], rbin[NMBIN], rval[NMBIN]; + double ymin, ymax; + char c; + Graphdata graphdata; + + fill_lumfunc (Target[0].stars, Target[0].N, tval, tbin, &Nt); + fill_lumfunc (Ref[0].stars, Ref[0].N, rval, rbin, &Nr); + + /* construct plots of log(stars / degree square) */ + ymin = 5; ymax = -5; + for (i = 0; i < Nt; i++) { + tval[i] = tval[i] - log10 (Target[0].Area); + ymin = MIN (tval[i], ymin); + ymax = MAX (tval[i], ymax); + } + + for (i = 0; i < Nr; i++) { + rval[i] = rval[i] - log10 (Ref[0].Area); + rbin[i] = rbin[i] + Ref[0].Moff; + ymin = MIN (rval[i], ymin); + ymax = MAX (rval[i], ymax); + } + + graphdata.xmin = MMIN; + graphdata.xmax = MMAX; + graphdata.ymin = ymin - 0.1; + graphdata.ymax = ymax + 0.1; + + graphdata.style = 1; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.ebar = 0; + graphdata.color = 0; + + graphdata.lweight = 0; + graphdata.size = 0.5; + + PlotReset (0); + + PrepPlotting (Nt, &graphdata, 0); + PlotVector (Nt, tbin, 0, 0); + PlotVector (Nt, tval, 1, 0); + + graphdata.color = 2; + graphdata.style = 1; + graphdata.lweight = 0; + PrepPlotting (Nr, &graphdata, 0); + PlotVector (Nr, rbin, 0, 0); + PlotVector (Nr, rval, 1, 0); + + /* plot truncated lum func */ + for (i = 0; i < Nr; i++) { + if (rbin[i] < Target[0].lum.Mmin - 0.5) rval[i] = -1; + if (rbin[i] > Target[0].lum.Mmax + 0.5) rval[i] = -1; + } + graphdata.ltype = 1; + PrepPlotting (Nr, &graphdata, 0); + PlotVector (Nr, rbin, 0, 0); + PlotVector (Nr, rval, 1, 0); + + DonePlotting (&graphdata, 0); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); + +} + + +void plot_resid (StarData *st, StarData *sr, Coords *coords) { + + int i; + double x, y, dx, dy; + float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2; + int Npair, *idx1, *idx2; + + Npair = pair_lists (&idx1, &idx2); + ALLOCATE (xvect0, float, Npair); + ALLOCATE (yvect0, float, Npair); + ALLOCATE (xvect1, float, Npair); + ALLOCATE (yvect1, float, Npair); + + ALLOCATE (xvect2, float, 2*Npair); + ALLOCATE (yvect2, float, 2*Npair); + + for (i = 0; i < Npair; i++) { + RD_to_XY (&x, &y, sr[idx2[i]].R, sr[idx2[i]].D, coords); + + dx = st[idx1[i]].X - x; + dy = st[idx1[i]].Y - y; + + xvect0[i] = st[idx1[i]].X; + xvect1[i] = st[idx1[i]].Y; + yvect0[i] = dx; + yvect1[i] = dy; + + xvect2[2*i+0] = st[idx1[i]].X; + yvect2[2*i+0] = st[idx1[i]].Y; + xvect2[2*i+1] = x; + yvect2[2*i+1] = y; + } + + fprintf (stderr, "residuals using RD_to_XY\n"); + plot_resid_plot (0, xvect0, yvect0, Npair); + // plot_resid_plot (1, xvect1, yvect1, Npair); + plot_fullfield_pairs (xvect2, yvect2, 2*Npair); + + for (i = 0; i < Npair; i++) { + fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y); + + dx = (x - sr[idx2[i]].P)/coords[0].cdelt1; + dy = (y - sr[idx2[i]].Q)/coords[0].cdelt2; + + xvect0[i] = st[idx1[i]].X; + xvect1[i] = st[idx1[i]].Y; + yvect0[i] = dx; + yvect1[i] = dy; + + xvect2[2*i+0] = sr[idx2[i]].P; + yvect2[2*i+0] = sr[idx2[i]].Q; + xvect2[2*i+1] = x; + yvect2[2*i+1] = y; + } + + fprintf (stderr, "residuals using fit_apply\n"); + plot_resid_plot (1, xvect0, yvect0, Npair); + // plot_resid_plot (1, xvect1, yvect1, Npair); + // plot_fullfield_pairs (xvect2, yvect2, 2*Npair); + + free (xvect2); + free (yvect2); + free (xvect0); + free (yvect0); + free (xvect1); + free (yvect1); + +} + +void fill_lumfunc (StarData *stars, int N, float *lbin, float *bin, int *nb) { + + int i, j, Nb; + double mbin[NMBIN]; + + bzero (mbin, NMBIN * sizeof (double)); + + /* sum histogram */ + for (i = 0; i < N; i++) { + if (stars[i].M < MMIN) continue; + if (stars[i].M > MMAX) continue; + + j = (stars[i].M - MMIN) / dM; + j = MIN (MAX (j, 0), (NMBIN - 1)); + mbin[j] ++; + } + + /* select filled bins */ + for (Nb = i = 0; i < NMBIN; i++) { + if (mbin[i] > 0) { + bin[Nb] = i * dM + MMIN; + lbin[Nb] = log (mbin[i]) / log (10.0); + Nb++; + } + } + + *nb = Nb; +} + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plotstuff.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plotstuff.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/plotstuff.c (revision 21563) @@ -0,0 +1,166 @@ +# include "gastro2.h" +# include + +static int Xgraph[5] = {0,0,0,0,0}; +static int active; + +void XDead () { + signal (SIGPIPE, XDead); + fprintf (stderr, "kapa is dead, must restart\n"); + Xgraph[active] = -1; +} + +int open_graph (int N) { + +# ifdef ANSI +# define F_SETFL 4 +# define O_NONBLOCK 0200000 +# define AF_UNIX 1 +# define SOCK_STREAM 1 +#define ENOENT 2 /* No such file or directory */ +# endif /* ANSI */ + + int InitSocket, status, addreslen; + struct sockaddr_un Address; + char temp[100]; + char socket_name[100]; + + active = N; + sprintf (socket_name, "/tmp/kapa%d", N); + sprintf (temp, "rm -f %s", socket_name); + system (temp); + + strcpy (Address.sun_path, socket_name); + Address.sun_family = AF_UNIX; + InitSocket = socket (AF_UNIX, SOCK_STREAM, 0); + status = bind (InitSocket, (struct sockaddr *) &Address, sizeof (Address)); + status = listen (InitSocket, 1); + + sprintf (temp, "kapa %s -name %d &", socket_name, N); +# ifndef DEBUG + system (temp); +# else + fprintf (stderr, "start kapa, press return: %s\n", temp); + fscanf (stdin, "%d", &i); +# endif + + addreslen = sizeof (Address); + Xgraph[N] = accept (InitSocket, (struct sockaddr *) &Address, &addreslen); + if (Xgraph[N] < 0) { + fprintf (stderr, "error starting kapa\n"); + return (FALSE); + } + else { + return (TRUE); + } + +} + +void DonePlotting (Graphdata *graphmode, int N) { + char buffer[65], buffer2[65]; + + write (Xgraph[N], "DBOX", 4); + sprintf (buffer, "%f %f %f %f", graphmode[0].xmin, graphmode[0].xmax, graphmode[0].ymin, graphmode[0].ymax); + sprintf (buffer2, "NBYTES: %6d", (unsigned int) strlen (buffer)); + write (Xgraph[N], buffer2, 16); + write (Xgraph[N], buffer, strlen (buffer)); + + sprintf (buffer, "%s %s %s", "2222", "2222", "2222"); + sprintf (buffer2, "NBYTES: %6d", (unsigned int) strlen (buffer)); + write (Xgraph[N], buffer2, 16); + write (Xgraph[N], buffer, strlen (buffer)); +} + +void PrepPlotting (int Npts, Graphdata *graphmode, int N) { + + int i, status; + char buffer[128], buffer2[128]; + + active = N; + if (Npts < 1) return; + + /* test Xgraph[N], flush junk from pipe */ + signal (SIGPIPE, XDead); + fcntl (Xgraph[N], F_SETFL, O_NONBLOCK); + for (i = 0; (read (Xgraph[N], buffer, 64) > 0) && (i < 20); i++); + fcntl (Xgraph[N], F_SETFL, !O_NONBLOCK); + + if (Xgraph[N] < 1) if (!open_graph(N)) return; + + /* tell kapa to look for the incoming image */ + status = write (Xgraph[N], "PLOT", 4); + + /* send Xgraph[N] the plot details */ + sprintf (buffer, "%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); + sprintf (buffer2, "NBYTES: %6d", (unsigned int) strlen (buffer)); + write (Xgraph[N], buffer2, 16); + write (Xgraph[N], buffer, strlen (buffer)); + + sprintf (buffer, "%f %f %f %f", + graphmode[0].xmin, graphmode[0].xmax, + graphmode[0].ymin, graphmode[0].ymax); + sprintf (buffer2, "NBYTES: %6d", (unsigned int) strlen (buffer)); + write (Xgraph[N], buffer2, 16); + write (Xgraph[N], buffer, strlen (buffer)); + +} + +void PlotVector (int Npts, float *vect, int mode, int N) { + + int Nbytes; + + if (Npts < 1) return; + + active = N; + Nbytes = Npts * sizeof (float); + write (Xgraph[N], vect, Nbytes); + +} + +void PlotReset (int N) { + + char buffer[128]; + int i; + + /* test Xgraph[N], flush junk from pipe */ + signal (SIGPIPE, XDead); + fcntl (Xgraph[N], F_SETFL, O_NONBLOCK); + for (i = 0; (read (Xgraph[N], buffer, 64) > 0) && (i < 20); i++); + fcntl (Xgraph[N], F_SETFL, !O_NONBLOCK); + + if (Xgraph[N] < 1) if (!open_graph(N)) return; + + write (Xgraph[N], "ERAS", 4); + +} +/* include these lines to plot a pair of vectors: + + typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, color; + double lweight, size; + } Graphdata; + Graphdata graphdata; + + graphdata.xmin = -200; + graphdata.xmax = 4200; + graphdata.ymin = -500; + graphdata.ymax = 500; + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 0.5; + + PrepPlotting (N, &graphdata, n); + PlotVector (N, Y, 0, n); + PlotVector (N, dM, 1, n); + DonePlotting (&graphdata, n); + + */ Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/polyfit.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/polyfit.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/polyfit.c (revision 21563) @@ -0,0 +1,408 @@ +# include "gastro2.h" + +static int NTERM, NPOWR, NPARS, NORDER, Npts; +static double **sum, **xsum, **ysum; +static double **matrix, **vector; + +void fit_init (int order) { + + int i; + + NORDER = order; + NPOWR = NORDER + 1; + NTERM = 2*NORDER + 1; + NPARS = (NORDER + 1)*(NORDER + 2) / 2; + Npts = 0; + + /* allocate arrays for fit solution */ + ALLOCATE (sum, double *, NTERM); + ALLOCATE (xsum, double *, NTERM); + ALLOCATE (ysum, double *, NTERM); + for (i = 0; i < NTERM; i++) { + ALLOCATE (sum[i], double, NTERM); + bzero (sum[i], NTERM*sizeof(double)); + ALLOCATE (xsum[i], double, NTERM); + bzero (xsum[i], NTERM*sizeof(double)); + ALLOCATE (ysum[i], double, NTERM); + bzero (ysum[i], NTERM*sizeof(double)); + } + ALLOCATE (matrix, double *, NPARS); + ALLOCATE (vector, double *, NPARS); + for (i = 0; i < NPARS; i++) { + ALLOCATE (matrix[i], double, NPARS); + ALLOCATE (vector[i], double, 2); + bzero (vector[i], 2*sizeof(double)); + bzero (matrix[i], NPARS*sizeof(double)); + } +} + +/* */ +void fit_add (double x1, double y1, double x2, double y2, double wt) { + + int n, m; + double xterm, yterm, term; + + xterm = 1; + for (n = 0; n < NTERM; n++) { + yterm = 1; + for (m = 0; m < NTERM; m++) { + term = xterm*yterm; + if (n+m < NTERM) { + sum[n][m] += term; + } + if (n+m < NPOWR) { + xsum[n][m] += x2*term; + ysum[n][m] += y2*term; + } + yterm *= y1; + } + xterm *= x1; + } + Npts ++; + +} + +void fit_eval () { + + int i, j, n, m, M, N; + double max; + + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + vector[i][0] = xsum[n][m]; + vector[i][1] = ysum[n][m]; + } + } + j = 0; + for (M = 0; M < NPOWR; M++) { + for (N = 0; N < NPOWR - M; N++, j++) { + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + matrix[i][j] = sum[n+N][m+M]; + } + } + } + } + max = 0.0; + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + max = MAX (max, fabs(matrix[i][j])); + } + max = MAX (max, fabs(vector[i][0])); + max = MAX (max, fabs(vector[i][1])); + } + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + matrix[i][j] /= max; + } + vector[i][0] /= max; + vector[i][1] /= max; + } + dgaussj (matrix, NPARS, vector, 2); + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + xsum[n][m] = vector[i][0]; + ysum[n][m] = vector[i][1]; + } + } + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + if (VERBOSE) fprintf (stderr, "RA x^%dy^%d: %10.4g DEC x^%dy^%d: %10.4g \n", n, m, vector[i][0], n, m, vector[i][1]); + } + } +} + +void fit_norm () { + + xsum[0][0] = 0; + xsum[1][0] = 1; + xsum[0][1] = 0; + + ysum[0][0] = 0; + ysum[1][0] = 0; + ysum[0][1] = 1; +} + +/* evaluate the fit at (X,Y) to yield (x,y) */ +void fit_apply (double *x, double *y, double X, double Y) { + + int m, n; + double xterm, yterm; + double Xo, Yo; + + Xo = Yo = 0; + yterm = 1; + for (m = 0; m < NPOWR; m++) { + xterm = 1; + for (n = 0; n < NPOWR - m; n++) { + Xo += xterm*yterm*xsum[n][m]; + Yo += xterm*yterm*ysum[n][m]; + xterm *= X; + } + yterm *= Y; + } + + *x = Xo; + *y = Yo; +} + +/* evaluate the x-derivative of the fit at (X,Y) to yield (x,y) */ +void fit_apply_dx (double *x, double *y, double X, double Y) { + + int m, n; + double xterm, yterm; + double Xo, Yo; + + Xo = Yo = 0; + yterm = 1; + for (m = 0; m < NPOWR; m++) { + xterm = 1; + for (n = 1; n < NPOWR - m; n++) { + Xo += n*xterm*yterm*xsum[n][m]; + Yo += n*xterm*yterm*ysum[n][m]; + xterm *= X; + } + yterm *= Y; + } + + *x = Xo; + *y = Yo; +} + +/* evaluate the y-derivative of the fit at (X,Y) to yield (x,y) */ +void fit_apply_dy (double *x, double *y, double X, double Y) { + + int m, n; + double xterm, yterm; + double Xo, Yo; + + Xo = Yo = 0; + yterm = 1; + for (m = 1; m < NPOWR; m++) { + xterm = 1; + for (n = 0; n < NPOWR - m; n++) { + Xo += m*xterm*yterm*xsum[n][m]; + Yo += m*xterm*yterm*ysum[n][m]; + xterm *= X; + } + yterm *= Y; + } + + *x = Xo; + *y = Yo; +} + +/* measure the residual scatter in the fit */ +double fit_scat (StarData *st, StarData *sr, Coords *coords) { + + int i; + int Npair, *idx1, *idx2; + double x, y, dx, dy, dX, dY, dX2, dY2, dR; + + Npair = pair_lists (&idx1, &idx2); + + dX = dY = dX2 = dY2 = 0; + for (i = 0; i < Npair; i++) { + + /* projection this direction includes the error introduced by + the interation on the nonlinear solution */ + RD_to_XY (&x, &y, sr[idx2[i]].R, sr[idx2[i]].D, coords); + + dx = x - st[idx1[i]].X; + dy = y - st[idx1[i]].Y; + + dX += dx; + dY += dy; + dX2 += dx*dx; + dY2 += dy*dy; + } + + /* scatter is measured on the tangent plane in degrees */ + dX = dX / Npair; + dY = dY / Npair; + fprintf (stderr, "scatter: %f, %f\n", sqrt(dX2/Npair - dX*dX), sqrt(dY2/Npair - dY*dY)); + fprintf (stderr, "Npts: %d\n", Npair); + + dR = 0.5 * sqrt(fabs(dX2/Npair - dX*dX)) + 0.5 * sqrt (fabs(dY2/Npair - dY*dY)); + return (dR); +} + +/* convert fit terms to coords and polyterms */ +/**** what do we do with the value of ctype??? + can we leave it alone? ****/ +int fit_adjust (Coords *coords) { + + int i; + double a10, a01, a20, a11, a02, a30, a21, a12, a03; + double b10, b01, b20, b11, b02, b30, b21, b12, b03; + double Xo, Yo, det; + double **A, **B, Fx, Fy; + + /* start with the linear solution for Xo,Yo */ + coords[0].cdelt1 = hypot (xsum[1][0], ysum[1][0]); + coords[0].cdelt2 = hypot (xsum[0][1], ysum[0][1]); + // coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + + det = 1.0 / (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]); + Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]); + Yo = det*(xsum[0][0]*ysum[1][0] - ysum[0][0]*xsum[1][0]); + + coords[0].Npolyterms = NORDER; + + if (coords[0].Npolyterms > 1) { + /* use the linear solution as a starting guess */ + /* solve for L(Xo,Yo) = 0, M(Xo,Yo) = 0 */ + /* this is the Newton-Raphson method - it needs the high order terms to be small */ + ALLOCATE (A, double *, 2); + ALLOCATE (B, double *, 2); + ALLOCATE (A[0], double, 2); + ALLOCATE (A[1], double, 2); + ALLOCATE (B[0], double, 1); + ALLOCATE (B[1], double, 1); + + for (i = 0; i < 10; i++) { + fit_apply (&Fx, &Fy, Xo, Yo); + fit_apply_dx (&A[0][0], &A[0][1], Xo, Yo); + fit_apply_dy (&A[1][0], &A[1][1], Xo, Yo); + B[0][0] = -Fx; + B[1][0] = -Fy; + dgaussj (A, 2, B, 1); + Xo += B[0][0]; + Yo += B[1][0]; + } + free (A[0]); free (B[0]); + free (A[1]); free (B[1]); + free (A); free (B); + } + coords[0].crpix1 = Xo; + coords[0].crpix2 = Yo; + + switch (coords[0].Npolyterms) { + case 0: + case 1: + /* the linear solution can be analytically inverted */ + coords[0].pc1_1 = xsum[1][0] / coords[0].cdelt1; + coords[0].pc1_2 = xsum[0][1] / coords[0].cdelt2; + coords[0].pc2_1 = ysum[1][0] / coords[0].cdelt1; + coords[0].pc2_2 = ysum[0][1] / coords[0].cdelt2; + for (i = 0; i < 7; i++) { + coords[0].polyterms[i][0] = coords[0].polyterms[i][1] = 0.0; + } + break; + + case 2: + a10 = xsum[1][0] + 2.0*xsum[2][0]*Xo + xsum[1][1]*Yo; + a01 = xsum[0][1] + 2.0*xsum[0][2]*Yo + xsum[1][1]*Xo; + a20 = xsum[2][0]; + a11 = xsum[1][1]; + a02 = xsum[0][2]; + + b10 = ysum[1][0] + 2.0*ysum[2][0]*Xo + ysum[1][1]*Yo; + b01 = ysum[0][1] + 2.0*ysum[0][2]*Yo + ysum[1][1]*Xo; + b20 = ysum[2][0]; + b11 = ysum[1][1]; + b02 = ysum[0][2]; + + coords[0].pc1_1 = a10 / coords[0].cdelt1; + coords[0].pc1_2 = a01 / coords[0].cdelt2; + coords[0].pc2_1 = b10 / coords[0].cdelt1; + coords[0].pc2_2 = b01 / coords[0].cdelt2; + + coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1); + coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2); + coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2); + + coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1); + coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2); + coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2); + for (i = 3; i < 7; i++) { + coords[0].polyterms[i][0] = coords[0].polyterms[i][1] = 0.0; + } + break; + + case 3: + a10 = xsum[1][0] + 2*xsum[2][0]*Xo + xsum[1][1]*Yo + 3*xsum[3][0]*Xo*Xo + 2*xsum[2][1]*Xo*Yo + xsum[1][2]*Yo*Yo; + a01 = xsum[0][1] + 2*xsum[0][2]*Yo + xsum[1][1]*Xo + 3*xsum[0][3]*Yo*Yo + 2*xsum[1][2]*Xo*Yo + xsum[2][1]*Xo*Xo; + a20 = xsum[2][0] + 3*xsum[3][0]*Xo + xsum[2][1]*Yo; + a11 = xsum[1][1] + 2*xsum[2][1]*Xo + 2*xsum[1][2]*Yo; + a02 = xsum[0][2] + 3*xsum[0][3]*Yo + xsum[1][2]*Xo; + a30 = xsum[3][0]; + a21 = xsum[2][1]; + a12 = xsum[1][2]; + a03 = xsum[0][3]; + + b10 = ysum[1][0] + 2*ysum[2][0]*Xo + ysum[1][1]*Yo + 3*ysum[3][0]*Xo*Xo + 2*ysum[2][1]*Xo*Yo + ysum[1][2]*Yo*Yo; + b01 = ysum[0][1] + 2*ysum[0][2]*Yo + ysum[1][1]*Xo + 3*ysum[0][3]*Yo*Yo + 2*ysum[1][2]*Xo*Yo + ysum[2][1]*Xo*Xo; + b20 = ysum[2][0] + 3*ysum[3][0]*Xo + ysum[2][1]*Yo; + b11 = ysum[1][1] + 2*ysum[2][1]*Xo + 2*ysum[1][2]*Yo; + b02 = ysum[0][2] + 3*ysum[0][3]*Yo + ysum[1][2]*Xo; + b30 = ysum[3][0]; + b21 = ysum[2][1]; + b12 = ysum[1][2]; + b03 = ysum[0][3]; + + coords[0].pc1_1 = a10 / coords[0].cdelt1; + coords[0].pc1_2 = a01 / coords[0].cdelt2; + coords[0].pc2_1 = b10 / coords[0].cdelt1; + coords[0].pc2_2 = b01 / coords[0].cdelt2; + + coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1); + coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2); + coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2); + + coords[0].polyterms[3][0] = a30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1); + coords[0].polyterms[4][0] = a21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2); + coords[0].polyterms[5][0] = a12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1); + coords[0].polyterms[6][0] = a03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2); + + coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1); + coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2); + coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2); + + coords[0].polyterms[3][1] = b30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1); + coords[0].polyterms[4][1] = b21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2); + coords[0].polyterms[5][1] = b12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1); + coords[0].polyterms[6][1] = b03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2); + break; + + default: + fprintf (stderr, "error: invalid order %d\n", coords[0].Npolyterms); + exit (2); + } + + while (coords[0].crval1 < 0) coords[0].crval1 += 360.0; + while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0; + + /* test for valid solution */ + return (Npts); +} + +int mk_polyterm (int n, int m, int norder) { + + int i, nt, N; + + N = 0; + nt = n + m; + for (i = 2; i < nt; i++) { + N += i + 1; + } + N += m; + return (N); +} + +int mk_vector (int n, int m, int norder) { + + int i, N; + + N = 0; + for (i = 0; i < m; i++) { + N += (norder - i + 1); + } + N += n; + return (N); +} + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/remove_clumps.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/remove_clumps.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/remove_clumps.c (revision 21563) @@ -0,0 +1,94 @@ +# include "gastro2.h" + +StarData *remove_clumps (StarData *instars, int *nstars, int NX, int NY) { + + int i, j, nx, ny, Npix, nn, x, y, pix, Nstars; + double s1, s2, mean, sigma, cutoff; + int *hist, *xcld; + int nxcld, Nxcld, Nout; + + StarData *outstars; + + /* create histogram of pixels using coords */ + + Nstars = *nstars; + + Npix = 100; + nx = (int) (NX / Npix) + 1; + ny = (int) (NY / Npix) + 1; + nn = nx * ny; + + ALLOCATE (hist, int, nn); + ALLOCATE (xcld, int, nn); + bzero (hist, nn*sizeof(int)); + + for (i = 0; i < Nstars; i++) { + + x = (int) (instars[i].X / Npix); + y = (int) (instars[i].Y / Npix); + pix = x + nx * y; + + if ((pix < 0) || (pix >= nn)) { + fprintf (stderr, "! %f %f %d %d %d %d\n", instars[i].X, instars[i].Y, x, y, pix, nn); + continue; + } + + hist[pix] ++; + + } + + /* find stats on histogram */ + s1 = s2 = 0; + for (i = 0; i < nn; i++) { + s1 += hist[i]; + s2 += hist[i]*hist[i]; + } + mean = s1 / nn; + sigma = 2 + sqrt (s2 / nn - mean*mean); + cutoff = mean + 5*sigma; + + /* identify clumps to exclude */ + Nxcld = 0; + for (i = 0; i < nn; i++) { + if (hist[i] > cutoff) { + y = (int) (i / nx); + x = i - y*nx; + fprintf (stderr, "cut: %d %d %d %d\n", i, x, y, hist[i]); + xcld[Nxcld] = i; + Nxcld ++; + } + } + + /* identify stars to exclude (type = -1) */ + for (i = 0; i < Nxcld; i++) { + + nxcld = 0; + for (j = 0; j < Nstars; j++) { + + x = (int) (instars[j].X / Npix); + y = (int) (instars[j].Y / Npix); + pix = x + nx * y; + if (pix != xcld[i]) continue; + nxcld ++; + instars[j].type = -1; + } + fprintf (stderr, "exclude %d in clump %d\n", nxcld, i); + } + + ALLOCATE (outstars, StarData, Nstars); + Nout = 0; + + for (i = 0; i < Nstars; i++) { + if (instars[i].type == -1) continue; + outstars[Nout] = instars[i]; + Nout ++; + } + + REALLOCATE (outstars, StarData, Nout); + *nstars = Nout; + fprintf (stderr, "keeping %d of %d stars\n", Nout, Nstars); + return (outstars); + +} + + Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rfits.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rfits.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rfits.c (revision 21563) @@ -0,0 +1,46 @@ +# include "gastro2.h" + +StarData *rfits (FILE *f, int *nstars) { + + int i, N, Nstars; + Header theader; + FTable table; + SMPData *stars; + StarData *stardata; + + /* init & load in table data */ + table.header = &theader; + if (!gfits_fread_ftable (f, &table, "SMPFILE")) goto escape; + + stars = gfits_table_get_SMPData (&table, &Nstars, NULL); + + ALLOCATE (stardata, StarData, Nstars); + for (i = N = 0; i < Nstars; i++) { + /* hardwired dophot exclusions should eventually be encapsulated elsewhere */ + if (stars[i].dophot == 4) continue; + if (stars[i].dophot == 5) continue; + if (stars[i].dophot == 6) continue; + if (stars[i].dophot == 9) continue; + if ((MAX_ERROR > 0) && (stars[i].dM > MAX_ERROR)) continue; + stardata[N].X = stars[i].X; + stardata[N].Y = stars[i].Y; + stardata[N].M = stars[i].M; + stardata[N].dM = 1000*stars[i].dM; + stardata[N].type = stars[i].dophot; + stardata[N].R = 0; + stardata[N].D = 0; + N++; + } + if (N < 5) { + fprintf (stderr, "ERROR: too few stars for reliable solution, only %d\n", N); + exit (1); + } + + REALLOCATE (stardata, StarData, MAX(N,1)); + *nstars = N; + return (stardata); + +escape: + fprintf (stderr, "error reading file\n"); + exit (1); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rotate2.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rotate2.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rotate2.c (revision 21563) @@ -0,0 +1,37 @@ +# include "gastro2.h" + +void rotate (RefCatalog *Subset, RefCatalog *Ref, double angle) { + + int i; + double dX, dY, DX, DY, CS, SN; + double theta; + StarData *in, *out; + + Ref[0] = Subset[0]; + + ALLOCATE (Ref[0].stars, StarData, MAX (1, Ref[0].N)); + bcopy (Subset[0].stars, Ref[0].stars, Ref[0].N*sizeof(StarData)); + + if (angle == 0.0) return; + + theta = (angle*RAD_DEG); + CS = cos (theta); + SN = sin (theta); + + in = Subset[0].stars; + out = Ref[0].stars; + + for (i = 0; i < Ref[0].N; i++) { + dX = in[i].X; + dY = in[i].Y; + + DX = dX * CS - dY * SN; + DY = dX * SN + dY * CS; + + out[i].X = DX; + out[i].Y = DY; + } + +} + +/* rotate the star list by an angle ccw from x axis */ Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rtext.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rtext.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/rtext.c (revision 21563) @@ -0,0 +1,49 @@ +# include "gastro2.h" +/* by necesity hard wired */ +# define D_NSTARS 1000 +# define BYTES_STAR 66 +# define BLOCK 1000 + +StarData *rtext (FILE *f, int *nstars) { + + char *buffer; + int i, N, Nbytes, nbytes, Ninstar, Nstars, NSTARS; + double dmag, type; + StarData *stars; + + NSTARS = *nstars; + ALLOCATE (stars, StarData, MAX (NSTARS, 1)); + Nbytes = NSTARS*BYTES_STAR; + + N = Nstars = 0; + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR)); + + while ((nbytes = fread (buffer, 1, (BLOCK*BYTES_STAR), f)) != 0) { + Ninstar = nbytes / BYTES_STAR; + for (i = 0; i < Ninstar; i++, Nstars++) { + dparse (&stars[N].X, 1, &buffer[i*BYTES_STAR]); + dparse (&stars[N].Y, 2, &buffer[i*BYTES_STAR]); + dparse (&stars[N].M, 3, &buffer[i*BYTES_STAR]); + dparse (&dmag, 4, &buffer[i*BYTES_STAR]); + dparse (&type, 5, &buffer[i*BYTES_STAR]); + + /* hardwired dophot exclusions should eventually be encapsulated elsewhere */ + if ((type == 4) || (type == 6) || (type == 5) || (type == 9)) continue; + if ((MAX_ERROR > 0) && (dmag > 1000*MAX_ERROR)) continue; + stars[N].dM = 0.001 * dmag; + stars[N].type = type; + N++; + } + } + free (buffer); + + if (Nstars != NSTARS) { + fprintf (stderr, "WARNING: only read %d of %d stars\n", Nstars, NSTARS); + } + if (N < 5) { + fprintf (stderr, "ERROR: too few stars for reliable solution, only %d\n", N); + exit (1); + } + *nstars = N; + return (stars); +} Index: /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/sort.c =================================================================== --- /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/sort.c (revision 21563) +++ /branches/ohana/gastro-2-4/Ohana/src/gastro2/src/sort.c (revision 21563) @@ -0,0 +1,205 @@ +# include "gastro2.h" + +void sort_lum (double *R, double *X, double *Y, int N) { + + int l,j,ir,i; + double tX, tY, tR; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tR = R[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tR = R[ir]; + R[ir] = R[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + R[0] = tR; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && R[j] < R[j+1]) j++; + if (tR < R[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + R[i] = R[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + R[i] = tR; + } +} + +void sort_lists (double *X, double *Y, int *S, int N) { + + int l,j,ir,i; + double tX, tY, tS; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tX = X[l]; + tY = Y[l]; + tS = S[l]; + } + else { + tX = X[ir]; + X[ir] = X[0]; + tY = Y[ir]; + Y[ir] = Y[0]; + tS = S[ir]; + S[ir] = S[0]; + if (--ir == 0) { + X[0] = tX; + Y[0] = tY; + S[0] = tS; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tX < X[j]) { + X[i] = X[j]; + Y[i] = Y[j]; + S[i] = S[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tX; + Y[i] = tY; + S[i] = tS; + } +} + +void sort (double *X, int N) { + + int l,j,ir,i; + double tmp; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + l--; + tmp = X[l]; + } + else { + tmp = X[ir]; + X[ir] = X[0]; + if (--ir == 0) { + X[0] = tmp; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && X[j] < X[j+1]) j++; + if (tmp < X[j]) { + X[i] = X[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + X[i] = tmp; + } +} + +void sort_stars_mag (StarData *stars, int N) { + + int l,j,ir,i; + StarData tempstar; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + tempstar = stars[--l]; + } + else { + tempstar = stars[ir]; + stars[ir] = stars[0]; + if (--ir == 0) { + stars[0] = tempstar; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && stars[j].M < stars[j+1].M) ++j; + if (tempstar.M < stars[j].M) { + stars[i] = stars[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + stars[i] = tempstar; + } +} + +void sort_stars_X (StarData *stars, int N) { + + int l,j,ir,i; + StarData tempstar; + + if (N < 2) return; + + l = N >> 1; + ir = N - 1; + for (;;) { + if (l > 0) { + tempstar = stars[--l]; + } + else { + tempstar = stars[ir]; + stars[ir] = stars[0]; + if (--ir == 0) { + stars[0] = tempstar; + return; + } + } + i = l; + j = (l << 1) + 1; + while (j <= ir) { + if (j < ir && stars[j].X < stars[j+1].X) ++j; + if (tempstar.X < stars[j].X) { + stars[i] = stars[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + stars[i] = tempstar; + } +} + +