Index: /tags/gastro-2-2/Ohana/src/gastro2/Makefile =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/Makefile (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/Makefile (revision 6697) @@ -0,0 +1,90 @@ +include ../../Configure +HOME = $(ROOT)/src/gastro2 +PROGRAM = gastro2 + +BIN = $(HOME)/bin +INC = $(HOME)/include +SRC = $(HOME)/src +MAN = $(HOME)/doc +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +# +INCS = -I$(INC) -I$(LINC) -I$(XINC) +LIBS = -L$(LLIB) -ldvo -lFITS -lohana -lsocket -lnsl -lm +CFLAGS = -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +2MASS = \ +$(SRC)/extr2mass.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/gregions2.$(ARCH).o + +GASTRO = \ +$(SRC)/getptolemy.$(ARCH).o \ +$(SRC)/gptolemy2.$(ARCH).o \ +$(SRC)/plots.$(ARCH).o \ +$(SRC)/gheader2.$(ARCH).o \ +$(SRC)/gfit2.$(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)/gregions2.$(ARCH).o \ +$(SRC)/gaussj.$(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)/gcatalog.$(ARCH).o \ +$(SRC)/remove_clumps.$(ARCH).o \ +$(SRC)/rfits.$(ARCH).o \ +$(SRC)/rtext.$(ARCH).o + +OBJ = $(GASTRO) + +default: $(PROGRAM) + +all: $(PROGRAM) extr2mass + +$(GASTRO): $(INC)/gastro2.h + +extr2mass: $(2MASS) + $(CC) $(2MASS) -o $(BIN)/extr2mass.linux $(CCFLAGS) + +# dependancy rules for binary code ########################## +$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + +$(BIN)/$(PROGRAM).$(ARCH): $(OBJ) + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi + $(CC) $(OBJ) -o $(BIN)/$(PROGRAM).$(ARCH) $(CCFLAGS) + +install: $(DESTBIN)/$(PROGRAM) + +$(DESTBIN)/$(PROGRAM): $(BIN)/$(PROGRAM).$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$(PROGRAM) + cp $(BIN)/$(PROGRAM).$(ARCH) $(DESTBIN)/$(PROGRAM) + +# utilities ################################################# +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $< + + Index: /tags/gastro-2-2/Ohana/src/gastro2/bin/.cvsignore =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/bin/.cvsignore (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/bin/.cvsignore (revision 6697) @@ -0,0 +1,2 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix +*.linrh Index: /tags/gastro-2-2/Ohana/src/gastro2/doc/ChangeLog.txt =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/doc/ChangeLog.txt (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/doc/ChangeLog.txt (revision 6697) @@ -0,0 +1,8 @@ + + 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: /tags/gastro-2-2/Ohana/src/gastro2/include/gastro2.h =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/include/gastro2.h (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/include/gastro2.h (revision 6697) @@ -0,0 +1,217 @@ +# include +# include + +typedef struct { + double R; + double D; + double X; + double Y; + double M; + double dM; + int type; +} StarData; + +# if (0) +typedef struct { + float X; + float Y; + float M; + float dM; + float Mgal; + float Map; + float sky; + float fx; + float fy; + float df; + char dophot; + char dummy[3]; +} Stars; +# endif + +typedef struct { + double 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; + +typedef struct { + double R, D; + double M; +} GSCdata; + +/* 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 CDROM[256]; +char USNO_B_DIR[256]; +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 FORCE; +double F_RA; +double F_DEC; + +char GSCFILE[256], GSC_DIR[256], LONEOS_REGION_FILE[256]; +char TWO_MASS_DIR[256]; +double ASEC_PIX; +char ROUGH_ASTROMETRY[64]; + +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); + +GSCdata *gptolemy (char *fullpath, int *Nstars); +CatStats *gregions (CatStats *patch, int *nregion); +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); +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); +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 get2mass (CatStats *catstats, RefCatalog *Ref); +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 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); +void plot_resid_init (int version, double xmax); +void plot_resid_plot (int version, float *xvect, float *yvect, int Nvect); +double ptolemy_area (GSCdata *ptolemy, int N, CatStats *patch); +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: /tags/gastro-2-2/Ohana/src/gastro2/src/ConfigInit.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/ConfigInit.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/ConfigInit.c (revision 6697) @@ -0,0 +1,63 @@ +# 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 (0); + } + if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file); + + ScanConfig (config, "OFFSET_RADIUS", "%lf", 0, &SEARCH_RADIUS); + ScanConfig (config, "MIN_MATCHES", "%d", 0, &MIN_MATCHES); + ScanConfig (config, "DEFAULT_RADIUS", "%lf", 0, &DEFAULT_RADIUS); + ScanConfig (config, "MINIMUM_RADIUS", "%lf", 0, &MINIMUM_RADIUS); + ScanConfig (config, "MAX_NONLINEAR", "%lf", 0, &MAX_NONLINEAR); + ScanConfig (config, "MAX_PRECISE", "%lf", 0, &MIN_PRECISE); + ScanConfig (config, "CCD_PC1_1", "%lf", 0, &CCD_PC1_1); + ScanConfig (config, "CCD_PC2_2", "%lf", 0, &CCD_PC2_2); + ScanConfig (config, "CCD_PC1_2", "%lf", 0, &CCD_PC1_2); + ScanConfig (config, "CCD_PC2_1", "%lf", 0, &CCD_PC2_1); + ScanConfig (config, "ASEC_PIX", "%lf", 0, &ASEC_PIX); + ScanConfig (config, "NFIELD", "%lf", 0, &NFIELD); + ScanConfig (config, "NPOLYTERMS", "%d", 0, &NPOLYTERMS); + ScanConfig (config, "ROT_ZERO", "%lf", 0, &ROT_ZERO); + ScanConfig (config, "dROT", "%lf", 0, &dROT); + ScanConfig (config, "NROT", "%d", 0, &NROT); + ScanConfig (config, "POLAR_ALIGNMENT", "%d", 0, &POLAR_ALIGNMENT); + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); + ScanConfig (config, "GSCDIR", "%s", 0, GSC_DIR); + ScanConfig (config, "2MASS_DIR", "%s", 0, TWO_MASS_DIR); + ScanConfig (config, "POLAR_AXIS_RA", "%lf", 0, &POLE_RA); + ScanConfig (config, "POLAR_AXIS_DEC", "%lf", 0, &POLE_DEC); + ScanConfig (config, "RA_OFFSET", "%lf", 0, &RA_OFFSET); + ScanConfig (config, "DEC_OFFSET", "%lf", 0, &DEC_OFFSET); + ScanConfig (config, "LONEOS_REGIONS", "%s", 0, LONEOS_REGION_FILE); + ScanConfig (config, "USNO_CDROM", "%s", 0, CDROM); + ScanConfig (config, "USNO_B_DIR", "%s", 0, USNO_B_DIR); + ScanConfig (config, "ASTRO_REFCAT", "%s", 0, REFCAT); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "CATMODE", "%s", 0, CATMODE); + ScanConfig (config, "CATFORMAT", "%s", 0, CATFORMAT); + ScanConfig (config, "ROUGH_ASTROMETRY", "%s", 0, ROUGH_ASTROMETRY); + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + ScanConfig (config, "GASTRO_MAX_NSTARS", "%d", 0, &GASTRO_MAX_NSTARS); + ScanConfig (config, "GASTRO_MAX_MAG_ERROR", "%lf", 0, &MAX_ERROR); + + if (*CATMODE == 0) strcpy (CATMODE, "RAW"); + if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR"); + + if (strcasecmp (ROUGH_ASTROMETRY, "header") && + strcasecmp (ROUGH_ASTROMETRY, "config")) { + fprintf (stderr, "ROUGH_ASTROMETRY must be one of: header, config\n"); + exit (0); + } + free (config); + free (file); +} Index: /tags/gastro-2-2/Ohana/src/gastro2/src/extr2mass.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/extr2mass.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/extr2mass.c (revision 6697) @@ -0,0 +1,183 @@ +# include "gastro2.h" + +# define NDECLIST 24 +double DecList[] = {-90.0, -82.5, -75.0, -67.5, -60.0, -52.5, -45.0, -37.5, -30.0, -22.5, -15.0, -7.5, 0.0, + +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, +90.0}; + +typedef struct { + char name[256]; + float rmin, rmax; + int Nstar, Nbyte, Ncbyte; +} TwoMassAcc; + +int main (int argc, char **argv) { + + int i, j, k, Nregions; + TwoMassAcc acc[49]; + CatStats patch, *regions; + char line[1024], filename[256], saveline[1024], *path, *Fullpath, *Basename, *Pathname; + FILE *f, *fout; + int status, done, Nacc, Nstar, NSTAR, Nout, Nscan; + double *ra, *dec, *mag, DecMin, DecMax; + + ConfigInit (&argc, argv); + + if (argc != 5) { + fprintf (stderr, "USAGE: extr2mass (input) (acc) (outdir) (gsctbl)\n"); + exit (2); + } + + /* load accelerator file */ + f = fopen (argv[2], "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't load accelerator file %s\n", argv[2]); + return (FALSE); + } + for (i = 0; i < 49; i++) { + status = fscanf (f, "%s %f %f %d %d %d", acc[i].name, &acc[i].rmin, &acc[i].rmax, + &acc[i].Nstar, &acc[i].Nbyte, &acc[i].Ncbyte); + if (status != 6) { + fprintf (stderr, "ERROR: can't read 2Mass accelerator file\n"); + return (FALSE); + } + acc[i].rmin *= 15.0; + acc[i].rmax *= 15.0; + } + fclose (f); + + /* find file in acc file */ + Nacc = -1; + for (i = 0; (Nacc == -1) && (i < 49); i++) { + if (!strcmp (argv[1], acc[i].name)) Nacc = i; + } + if (Nacc == -1) { + fprintf (stderr, "can't find file in acc file\n"); + exit (1); + } + + fprintf (stderr, "reading from %s, %d stars\n", acc[Nacc].name, acc[Nacc].Nstar); + + /* load 2mass data file */ + /* + sprintf (line, "gzip -cd %s > %s", acc[Nacc].name, acc[Nacc].name); + system (line); + */ + + sprintf (line, "gzip -cd %s", acc[Nacc].name); + f = popen (line, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't read 2Mass data file: %s\n", line); + return (FALSE); + } + + /* all data in 2mass file are sequential in DEC. read stars while in dec range */ + + /* load data in DEC blocks */ + saveline[0] = 0; + for (i = 0; i < NDECLIST; i++) { + + Nstar = 0; + NSTAR = 1000; + ALLOCATE (ra, double, NSTAR); + ALLOCATE (dec, double, NSTAR); + ALLOCATE (mag, double, NSTAR); + + DecMin = DecList[i]; + DecMax = DecList[i+1]; + + done = FALSE; + while (!done) { + if (saveline[0]) { + strcpy (line, saveline); + saveline[0] = 0; + } else { + if (scan_line (f, line) == EOF) { goto escape; } + } + Nscan = sscanf (line, "%lf %lf %*s %*s %*s %*s %lf", &ra[Nstar], &dec[Nstar], &mag[Nstar]); + if (Nscan != 3) { + fprintf (stderr, "read wrong number of entries: %d\n", Nscan); + fprintf (stderr, "%s\n", line); + exit (1); + } + if (dec[Nstar] > DecMax) { + strcpy (saveline, line); + goto escape; + } + Nstar ++; + if (Nstar == NSTAR) { + NSTAR += 1000; + REALLOCATE (ra, double, NSTAR); + REALLOCATE (dec, double, NSTAR); + REALLOCATE (mag, double, NSTAR); + } + if (Nstar > acc[Nacc].Nstar) { + fprintf (stderr, "problem!\n %s\n", line); + exit (1); + } + } + escape: + + fprintf (stderr, "read %d stars for DEC block %f %f\n", Nstar, DecMin, DecMax); + + /* now we write the Nstar objects to the appropriate files */ + patch.RA[0] = acc[Nacc].rmin; + patch.RA[1] = acc[Nacc].rmax; + patch.DEC[0] = DecMin; + patch.DEC[1] = DecMax; + + regions = gregions (&patch, &Nregions); + + for (j = 0; j < Nregions; j++) { + /* name has complete path to gsc, strip off path - 1 */ + Fullpath = pathname (regions[j].name); + Basename = filebasename (regions[j].name); + Pathname = filebasename (Fullpath); + sprintf (filename, "%s/%s/%s", argv[3], Pathname, Basename); + + path = pathname (filename); + mkdirhier (path, S_IRWXU | S_IRWXG | S_IRWXO); + + fout = fopen (filename, "a+"); + if (fout == (FILE *) NULL) { + fprintf (stderr, "can't open output file %s\n", filename); + exit (1); + } + + Nout = 0; + for (k = 0; k < Nstar; k++) { + if (ra[k] < regions[j].RA[0]) continue; + if (ra[k] >= regions[j].RA[1]) continue; + if (dec[k] < regions[j].DEC[0]) continue; + if (dec[k] >= regions[j].DEC[1]) continue; + fprintf (fout, "%8.4f %8.4f %4.1f\n", ra[k], dec[k], mag[k]); + Nout ++; + } + if (Nout > 0) { + fprintf (stderr, "wrote %d stars to %s\n", Nout, filename); + } + + free (path); + free (Fullpath); + free (Basename); + free (Pathname); + fclose (fout); + } + + free (ra); + free (dec); + free (mag); + } + + fclose (f); + exit (0); +} + + +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; +} + Index: /tags/gastro-2-2/Ohana/src/gastro2/src/gargs.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gargs.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gargs.c (revision 6697) @@ -0,0 +1,125 @@ +# 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); + } + + LONEOS_COORDS = FALSE; + if ((N = get_argument (*argc, argv, "-loneos"))) { + LONEOS_COORDS = TRUE; + remove_argument (N, argc, argv); + } + + 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); + } + + 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); + } + + 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); + } + + 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); + } + + if (*argc != NARGS) { + fprintf (stderr, "USAGE: gastro filename\n"); + exit (0); + } + +} + Index: /tags/gastro-2-2/Ohana/src/gastro2/src/gastro2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gastro2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gastro2.c (revision 6697) @@ -0,0 +1,91 @@ +# 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 (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 < 3; i++) { + gfit (&Target, &Ref, MIN (MAX (1, NPOLYTERMS), 3)); + fprintf (stderr, "precision: %f\n", Target.answer.dR / sqrt(Target.answer.N)); + } + + 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: /tags/gastro-2-2/Ohana/src/gastro2/src/gaussj.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gaussj.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gaussj.c (revision 6697) @@ -0,0 +1,77 @@ +# include + +int gaussj (double **a, int n, double **b, int m) { + + int *indxc,*indxr,*ipiv; + int i, icol, irow, j, k, l, ll; + double big,dum,pivinv; + + irow = icol = 0; + + ALLOCATE (indxc, int, n); + ALLOCATE (indxr, int, n); + ALLOCATE (ipiv, int, n); + for (j = 0; j < n; j++) + ipiv[j] = 0; + + for (i = 0; i < n; i++) { + big = 0.0; + for (j = 0; j < n; j++) { + if (ipiv[j] != 1) { + for (k = 0; k < n; k++) { + if (ipiv[k] == 0) { + if (fabs (a[j][k]) >= big) { + big = fabs (a[j][k]); + irow = j; + icol = k; + } + } + else + if (ipiv[k] > 1) { + fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); + return (0); + } + } + } + } + ipiv[icol]++; + if (irow != icol) { + for (l = 0; l < n; l++) + SWAP (a[irow][l], a[icol][l]); + for (l = 0; l < m; l++) + SWAP (b[irow][l], b[icol][l]); + } + indxr[i] = irow; + indxc[i] = icol; + if (a[icol][icol] == 0.0) { + fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); + return (0); + } + pivinv = 1.0 / a[icol][icol]; + a[icol][icol] = 1.0; + for (l = 0; l < n; l++) + a[icol][l] *= pivinv; + for (l = 0; l < m; l++) + b[icol][l] *= pivinv; + for (ll = 0; ll < n; ll++) { + if (ll != icol) { + dum = a[ll][icol]; + a[ll][icol] = 0.0; + for (l = 0; l < n; l++) + a[ll][l] -= a[icol][l]*dum; + for (l = 0; l < m; l++) + b[ll][l] -= b[icol][l]*dum; + } + } + } + + for (l = n - 1; l >= 0; l--) { + if (indxr[l] != indxc[l]) + for (k = 0; k < n; k++) + SWAP (a[k][indxr[l]], a[k][indxc[l]]); + } + free (ipiv); + free (indxr); + free (indxc); + return (1); +} Index: /tags/gastro-2-2/Ohana/src/gastro2/src/gcatalog.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gcatalog.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gcatalog.c (revision 6697) @@ -0,0 +1,125 @@ +# include "gastro2.h" +CatStats *gregions (CatStats *patch, int *nregion); +GSCdata *gcatalog (char *filename, int *Nstars); + +int get2mass (CatStats *catstats, RefCatalog *Ref) { + + int i, j, k, Ns, Ngsc, Nregions; + GSCdata *gsc; + CatStats *regions; + + Ref[0].N = 0; + Ref[0].Area = 0; + ALLOCATE (Ref[0].stars, StarData, 1); + + set_catalog (TWO_MASS_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + + Ns = Ref[0].N; + Ref[0].N += Ngsc; + Ref[0].Area += regions[i].Area; + + 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; + } + } + + 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 2MASS\n", Ref[0].N); + return (TRUE); +} + +int getgsc (CatStats *catstats, RefCatalog *Ref) { + + int i, j, k, Ns, Ngsc, Nregions; + GSCdata *gsc; + CatStats *regions; + + Ref[0].N = 0; + Ref[0].Area = 0; + ALLOCATE (Ref[0].stars, StarData, 1); + + set_catalog (GSC_DIR); + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + gsc = gcatalog (regions[i].name, &Ngsc); + + Ns = Ref[0].N; + Ref[0].N += Ngsc; + Ref[0].Area += regions[i].Area; + + 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; + } + } + + 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 +GSCdata *gcatalog (char *filename, int *Nstars) { + + GSCdata *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, GSCdata, 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, GSCdata, NSTAR); + } + } + } + + free (buffer); + + *Nstars = nstar; + return (stars); +} + Index: /tags/gastro-2-2/Ohana/src/gastro2/src/gcenter2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gcenter2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gcenter2.c (revision 6697) @@ -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: /tags/gastro-2-2/Ohana/src/gastro2/src/getptolemy.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/getptolemy.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/getptolemy.c (revision 6697) @@ -0,0 +1,82 @@ +# include "gastro2.h" + +int getptolemy (CatStats *catstats, RefCatalog *Ref) { + + int i, j, k, Ns, Nptolemy, Nregions; + double FracArea; + GSCdata *ptolemy; + CatStats *regions; + + Ref[0].N = 0; + Ref[0].Area = 0; + ALLOCATE (Ref[0].stars, StarData, 1); + + regions = gregions (catstats, &Nregions); + + for (i = 0; i < Nregions; i++) { + ptolemy = gptolemy (regions[i].name, &Nptolemy); + if (Nptolemy == 0) continue; + + FracArea = ptolemy_area (ptolemy, Nptolemy, ®ions[i]); + + Ns = Ref[0].N; + Ref[0].N += Nptolemy; + Ref[0].Area += regions[i].Area * FracArea; + + REALLOCATE (Ref[0].stars, StarData, MAX (1, Ref[0].N)); + for (k = Ns, j = 0; j < Nptolemy; k++, j++) { + Ref[0].stars[k].R = ptolemy[j].R; + Ref[0].stars[k].D = ptolemy[j].D; + Ref[0].stars[k].M = ptolemy[j].M; + } + free (ptolemy); + } + + 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 ptolemy_area (GSCdata *ptolemy, int N, CatStats *patch) { + + int i, xb, yb, Nb; + int bin[10][10]; + double frac, Rmin, Rmax, Dmin, Dmax, dR, dD; + + Rmin = patch[0].RA[0]; + Rmax = patch[0].RA[1]; + Dmin = patch[0].DEC[0]; + Dmax = patch[0].DEC[1]; + 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 < N; i++) { + xb = MAX (MIN (0, 10 * (ptolemy[i].R - Rmin) / dR), 9); + yb = MAX (MIN (0, 10 * (ptolemy[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: /tags/gastro-2-2/Ohana/src/gastro2/src/getusno.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/getusno.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/getusno.c (revision 6697) @@ -0,0 +1,145 @@ +# 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], c; + 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 cdrom for spd %d\n", spd); + exit (0); + } + + /* load accelerator file */ + sprintf (filename, "%s/zone%04d.acc", CDROM, spd); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (1); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; + else + last = 1 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", CDROM, spd); + 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: /tags/gastro-2-2/Ohana/src/gastro2/src/getusnob.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/getusnob.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/getusnob.c (revision 6697) @@ -0,0 +1,167 @@ +# include "gastro2.h" +# define NZONE 180 + +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], c; + 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, disk; + 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); + fprintf (stderr, "reading from %s\n", filename); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk); + fprintf (stderr, "press return when ready to continue: "); + fscanf (stdin, "%c", &c); + fprintf (stderr, "trying again...\n"); + f = fopen (filename, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (1); + } + } + for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++); + Nbins = i; + fclose (f); + + first = catstats[0].RA[0] / 3.75; + if ((catstats[0].RA[1] / 3.75) == (int) (catstats[0].RA[1] / 3.75)) + last = catstats[0].RA[1] / 3.75; + else + last = 1 + catstats[0].RA[1] / 3.75; + + if ((first > Nbins) || (last > Nbins)) { + fprintf (stderr, "ERROR: RA out of range\n"); + exit (1); + } + + /* open data file */ + sprintf (filename, "%s/zone%04d.cat", CDROM, spd); + 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 */ + /**** 80 is record-length for USNO-B ****/ + offset = 80*sizeof(int)*(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 = 80*Nstars/4; /* number of integer blocks; need to use Fread for byte-swapping read */ + + /* allocate space for stars in segment */ + 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 < Nstars; i++, buf+=20) { + 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: /tags/gastro-2-2/Ohana/src/gastro2/src/gfit2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gfit2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gfit2.c (revision 6697) @@ -0,0 +1,446 @@ +# include "gastro2.h" + +static int NPAIR, Npair; +static int *idx1, *idx2; + +void gfit (CmpCatalog *Target, RefCatalog *Ref, int order) { + + int i, j, j0; + double Radius, Radius2; + double dX, dY, dR; + RefCatalog Subset; + StarData *st, *sr; + + 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.0 * 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 ++; + } + +# if (0) + pair_init (); + for (i = 0; i < Target[0].N; i++) { + pair_add (i, i); + } +# endif + + fit_init (order); + fit_norm (); + + plot_resid_init (0, (double) Target[0].header.Naxis[0]); + plot_resid_init (1, (double) Target[0].header.Naxis[1]); + if (PLOTSTUFF) plot_resid (st, sr); + + /* 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]].X, sr[idx2[i]].Y, 1.0); + } + fit_eval (); + + Target[0].answer.dR = fit_scat (st, sr); + Target[0].answer.N = fit_adjust (&Target[0].coords); + + plot_resid_init (0, (double) Target[0].header.Naxis[0]); + plot_resid_init (1, (double) Target[0].header.Naxis[1]); + if (PLOTSTUFF) plot_resid (st, sr); + free (idx1); + free (idx2); +} + +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; + + fprintf (stderr, "npts: %d\n", Npts); + + 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; + } +# if (0) + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + fprintf (stderr, "%10.4e ", matrix[i][j]); + } + fprintf (stderr, " %10.4e ", vector[i][0]); + fprintf (stderr, " %10.4e \n", vector[i][1]); + } +# endif + gaussj (matrix, NPARS, vector, 2); +# if (0) + fprintf (stderr, "\n\n"); + for (i = 0; i < NPARS; i++) { + for (j = 0; j < NPARS; j++) { + fprintf (stderr, "%11.4e ", matrix[i][j]); + } + fprintf (stderr, " %11.4e ", vector[i][0]); + fprintf (stderr, " %11.4e \n", vector[i][1]); + } +# endif + 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++) { + 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; +} + +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; +} + + +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); + } + +} + +double fit_scat (StarData *st, StarData *sr) { + + int i; + double x, y, dx, dy, dX, dY, dX2, dY2, dR; + + dX = dY = dX2 = dY2 = 0; + for (i = 0; i < Npair; i++) { + + fit_apply (&x, &y, sr[idx2[i]].X, sr[idx2[i]].Y); + + dx = x - st[idx1[i]].X; + dy = y - st[idx1[i]].Y; + + dX += dx; + dY += dy; + dX2 += dx*dx; + dY2 += dy*dy; + + } + + dX = dX / Npair; + dY = dY / Npair; + fprintf (stderr, "scatter: %f, %f\n", sqrt(dX2/Npair - dX*dX), sqrt(dY2/Npair - dY*dY)); + fprintf (stderr, "precise: %f, %f\n", sqrt((dX2/Npair - dX*dX) / Npair), sqrt((dY2/Npair - dY*dY)/Npair)); + + dR = 0.5 * sqrt(fabs(dX2/Npair - dX*dX)) + 0.5 * sqrt (fabs(dY2/Npair - dY*dY)); + return (dR); +} + +/* convert new terms to adjustments in coords and to polyterms */ +int fit_adjust (Coords *coords) { + + int i, j, N; + double S1, S2, p11, p12, p21, p22; + double a0, a1, a2, b0, b1, b2, det; + double X, Y; + int Np, Nv; + + S1 = coords[0].cdelt1; + S2 = coords[0].cdelt2; + p11 = coords[0].pc1_1; p12 = coords[0].pc1_2; + p21 = coords[0].pc2_1; p22 = coords[0].pc2_2; + + /* get the correct vector entries for the linear terms */ + N = mk_vector (0, 0, NORDER); + a0 = vector[N][0]; b0 = vector[N][1]; + N = mk_vector (1, 0, NORDER); + a1 = vector[N][0]; b1 = vector[N][1]; + N = mk_vector (0, 1, NORDER); + a2 = vector[N][0]; b2 = vector[N][1]; + + det = 1.0 / (a1*b2 - a2*b1); + + coords[0].pc1_1 = p11*a1 + p12*b1*(S2/S1); + coords[0].pc2_1 = p21*a1 + p22*b1*(S2/S1); + + coords[0].pc1_2 = p12*b2 + p11*a2*(S1/S2); + coords[0].pc2_2 = p22*b2 + p21*a2*(S1/S2); + + X = (coords[0].crpix1 - a0); + Y = (coords[0].crpix2 - b0); + coords[0].crpix1 = det*(X*b2 - Y*a2); + coords[0].crpix2 = det*(Y*a1 - X*b1); + + coords[0].Npolyterms = NORDER; + if (NORDER > 1) strcpy (coords[0].ctype, "DEC--PLY"); + + /* generate higher order terms from vector */ + for (i = 0; i < NORDER + 1; i++) { + for (j = 0; j < (NORDER - i + 1); j++) { + if (i + j < 2) continue; + Np = mk_polyterm (i, j, NORDER); + Nv = mk_vector (i, j, NORDER); + coords[0].polyterms[Np][0] = det*(vector[Nv][0]*b2 - vector[Nv][1]*a2); /* x2 y0 */ + coords[0].polyterms[Np][1] = det*(vector[Nv][1]*a1 - vector[Nv][0]*b1); /* x2 y0 */ + } + } + 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); +} + +void plot_resid (StarData *st, StarData *sr) { + + int i; + double x, y, dx, dy; + float *xvect0, *yvect0, *xvect1, *yvect1, *xvect2, *yvect2; + int Nvect; + + Nvect = Npair; + ALLOCATE (xvect0, float, Nvect); + ALLOCATE (yvect0, float, Nvect); + ALLOCATE (xvect1, float, Nvect); + ALLOCATE (yvect1, float, Nvect); + + ALLOCATE (xvect2, float, 2*Nvect); + ALLOCATE (yvect2, float, 2*Nvect); + + for (i = 0; i < Npair; i++) { + fit_apply (&x, &y, st[idx1[i]].X, st[idx1[i]].Y); + + dx = x - sr[idx2[i]].X; + dy = y - sr[idx2[i]].Y; + + xvect0[i] = sr[idx2[i]].X; + xvect1[i] = sr[idx2[i]].Y; + yvect0[i] = dx; + yvect1[i] = dy; + + xvect2[2*i+0] = sr[idx2[i]].X; + yvect2[2*i+0] = sr[idx2[i]].Y; + xvect2[2*i+1] = x; + yvect2[2*i+1] = y; + } + + plot_resid_plot (0, xvect0, yvect0, Nvect); + plot_resid_plot (1, xvect1, yvect1, Nvect); + plot_fullfield_pairs (xvect2, yvect2, 2*Nvect); + + free (xvect2); + free (yvect2); + free (xvect0); + free (yvect0); + free (xvect1); + free (yvect1); + +} Index: /tags/gastro-2-2/Ohana/src/gastro2/src/gheader2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gheader2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gheader2.c (revision 6697) @@ -0,0 +1,124 @@ +# 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 (!fits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (3)\n", file); + exit(0); + } + oldsize = header.size; + + if (Target[0].answer.N < 6) { + fits_modify (&header, "NASTRO", "%d", 1, 0); + fits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + + fits_modify (&header, "NASTRO", "%d", 1, Target[0].answer.N); + fits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + + /*** use PutCoords to update header ***/ + if (Target[0].coords.Npolyterms > 1) { + fits_modify (&header, "CTYPE1", "%s", 1, "RA---PLY"); + fits_modify (&header, "CTYPE2", "%s", 1, "DEC--PLY"); + } else { + fits_modify (&header, "CTYPE1", "%s", 1, "RA---TAN"); + fits_modify (&header, "CTYPE2", "%s", 1, "DEC--TAN"); + } + fits_modify (&header, "CDELT1", "%le", 1, Target[0].coords.cdelt1); + fits_modify (&header, "CDELT2", "%le", 1, Target[0].coords.cdelt2); + fits_modify (&header, "CRVAL1", "%lf", 1, Target[0].coords.crval1); + fits_modify (&header, "CRVAL2", "%lf", 1, Target[0].coords.crval2); + fits_modify (&header, "CRPIX1", "%lf", 1, Target[0].coords.crpix1); + fits_modify (&header, "CRPIX2", "%lf", 1, Target[0].coords.crpix2); + fits_modify (&header, "PC001001", "%le", 1, Target[0].coords.pc1_1); + fits_modify (&header, "PC001002", "%le", 1, Target[0].coords.pc1_2); + fits_modify (&header, "PC002001", "%le", 1, Target[0].coords.pc2_1); + fits_modify (&header, "PC002002", "%le", 1, Target[0].coords.pc2_2); + fits_modify (&header, "NPLYTERM", "%d", 1, Target[0].coords.Npolyterms); + if (Target[0].coords.Npolyterms > 1) { + /* RA Terms */ + fits_modify (&header, "PCA1X2Y0", "%le", 1, Target[0].coords.polyterms[0][0]); /* polyterms[0]); */ + fits_modify (&header, "PCA1X1Y1", "%le", 1, Target[0].coords.polyterms[1][0]); /* polyterms[1]); */ + fits_modify (&header, "PCA1X0Y2", "%le", 1, Target[0].coords.polyterms[2][0]); /* polyterms[2]); */ + + if (Target[0].coords.Npolyterms > 2) { + fits_modify (&header, "PCA1X3Y0", "%le", 1, Target[0].coords.polyterms[3][0]); /* polyterms[3]); */ + fits_modify (&header, "PCA1X2Y1", "%le", 1, Target[0].coords.polyterms[4][0]); /* polyterms[4]); */ + fits_modify (&header, "PCA1X1Y2", "%le", 1, Target[0].coords.polyterms[5][0]); /* polyterms[5]); */ + fits_modify (&header, "PCA1X0Y3", "%le", 1, Target[0].coords.polyterms[6][0]); /* polyterms[6]); */ + } + /* Dec Terms */ + fits_modify (&header, "PCA2X2Y0", "%le", 1, Target[0].coords.polyterms[0][1]); /* polyterms[7]); */ + fits_modify (&header, "PCA2X1Y1", "%le", 1, Target[0].coords.polyterms[1][1]); /* polyterms[8]); */ + fits_modify (&header, "PCA2X0Y2", "%le", 1, Target[0].coords.polyterms[2][1]); /* polyterms[9]); */ + + if (Target[0].coords.Npolyterms > 2) { + fits_modify (&header, "PCA2X3Y0", "%le", 1, Target[0].coords.polyterms[3][1]); /* polyterms[10]); */ + fits_modify (&header, "PCA2X2Y1", "%le", 1, Target[0].coords.polyterms[4][1]); /* polyterms[11]); */ + fits_modify (&header, "PCA2X1Y2", "%le", 1, Target[0].coords.polyterms[5][1]); /* polyterms[12]); */ + fits_modify (&header, "PCA2X0Y3", "%le", 1, Target[0].coords.polyterms[6][1]); /* polyterms[13]); */ + } + } + + dR = fabs (Target[0].answer.dR*Target[0].coords.cdelt1*3600.0); + fits_modify (&header, "CERROR", "%lf", 1, dR); + fits_modify (&header, "CERROR", "%C", 1, "scatter in astrometry soln (arcsec)"); + fits_modify (&header, "CPRECISE", "%lf", 1, dR / sqrt(1.0*Target[0].answer.N)); + fits_modify (&header, "CPRECISE", "%C", 1, "precision of astrometry soln (arcsec)"); + fits_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: /tags/gastro-2-2/Ohana/src/gastro2/src/gproject2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gproject2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gproject2.c (revision 6697) @@ -0,0 +1,78 @@ +# include "gastro2.h" + +void gproject (CmpCatalog *Target, RefCatalog *Ref, RefCatalog *Subset) { + + int i, N; + double X, Y, M, Moff; + double XMIN, XMAX, YMIN, YMAX, MMIN, MMAX; + Coords *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; + + 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 */ + MMAX = Target[0].lum.Mmax + 0.5; + MMIN = Target[0].lum.Mmin - 0.5; + + 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); + + 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; + } + + out[N] = in[i]; + out[N].X = X; + out[N].Y = Y; + 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)); + +} + +/* 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: /tags/gastro-2-2/Ohana/src/gastro2/src/gptolemy2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gptolemy2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gptolemy2.c (revision 6697) @@ -0,0 +1,56 @@ +# include "gastro2.h" + +/* filename has path to GSC, strip off path, use filename with CATDIR path */ + +GSCdata *gptolemy (char *fullpath, int *Nstars) { + + int i, nstar; + char filename[256], *fileroot, *filepath, *subdir; + Catalog catalog; + GSCdata *stars; + + filepath = pathname (fullpath); + fileroot = filebasename (fullpath); + subdir = filebasename (filepath); + sprintf (filename, "%s/%s/%s", CATDIR, subdir, fileroot); + + *Nstars = 0; + catalog.filename = filename; + switch (lock_catalog (&catalog, LCK_SOFT)) { + case 0: + case 2: + fprintf (stderr, "can't lock catalog data %s\n", filename); + return ((GSCdata *)NULL); + case 1: + break; + } + + /* CATMODE and CATFORMAT determined from catalog */ + catalog.catflags = LOAD_AVES | LOAD_MEAS; + if (!load_catalog (&catalog, FALSE)) { + fprintf (stderr, "can't load catalog data %s\n", filename); + return ((GSCdata *)NULL); + } + unlock_catalog (&catalog); + + nstar = 0; + *Nstars = catalog.Naverage; + if (catalog.Naverage == 0) return ((GSCdata *)NULL); + + ALLOCATE (stars, GSCdata, *Nstars); + + /* select all entries, ignore magnitudes */ + for (i = 0; i < catalog.Naverage; i++, nstar++) { + stars[nstar].R = catalog.average[i].R; + stars[nstar].D = catalog.average[i].D; + stars[nstar].M = catalog.measure[catalog.average[i].offset].M_PS; + } + + free (catalog.average); + + *Nstars = nstar; + REALLOCATE (stars, GSCdata, *Nstars); + + return (stars); + +} Index: /tags/gastro-2-2/Ohana/src/gastro2/src/greference2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/greference2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/greference2.c (revision 6697) @@ -0,0 +1,88 @@ +# 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 (!strcmp (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 (!strcmp (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 (!strcmp (REFCAT, "GSC")) { + getgsc (&catstats, Ref); + } + + /* get stars from 2MASS for the given region -- add PHOTCODE check? */ + if (!strcmp (REFCAT, "2MASS")) { + strcpy (CATDIR, TWO_MASS_DIR); + getptolemy (&catstats, Ref); + /* get2mass (&catstats, Ref); */ + } + + /* get stars from the DVO CATDIR for the given region */ + if (!strcmp (REFCAT, "PTOLEMY")) { + getptolemy (&catstats, Ref); + } + + if (Ref[0].N == 0) { + fprintf (stderr, "no ref objs: %s\n", REFCAT); + exit (1); + } + + Ref[0].lum.Mz = (2 + log (Ref[0].Area) / log (10.) - Ref[0].lum.Mo) / Ref[0].lum.dNdM; + Target[0].lum.Mz = (2 + log (Target[0].Area) / log (10.) - Target[0].lum.Mo) / Target[0].lum.dNdM; + Ref[0].Moff = Target[0].lum.Mz - Ref[0].lum.Mz; + fprintf (stderr, "mag offset: %f (Areas: %f vs %f)\n", Ref[0].Moff, Target[0].Area, Ref[0].Area); + + 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; + 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; + + catstats[0].RA[0] = catstats[0].DEC[0] = 360.0; + catstats[0].RA[1] = catstats[0].DEC[1] = -360.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; + XY_to_RD (&R, &D, X, Y, &Target[0].coords); + + 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); + } + } + 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: /tags/gastro-2-2/Ohana/src/gastro2/src/gregions2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gregions2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gregions2.c (revision 6697) @@ -0,0 +1,247 @@ +# include "gastro2.h" +# define MAX_NAME 256 + +static int Nregion, NREGION; +static CatStats *region; +static char *CatDir; + +/* data which defines the layout of the regions */ +# define NDECBANDS 12 +double DecBand[] = {0.0, +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, +90.0, + 0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, -52.5, -60.0, -67.5, -75.0, -82.5, -90.0}; + +int NDecLines[] = {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, 0, + 597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48, 0}; + +int NRaBlocks [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0, + 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0}; + +char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", + "n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", + "s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", + "s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"}; + +/* given patch on sky with RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */ +CatStats *gregions (CatStats *patch, int *nregion) { + + CatStats area; + double dtmp; + char file[MAX_NAME]; + + init_regions (); + + /* force RA[i] to be in range 0 - 360 , DEC[0] < DEC[1] */ + while (patch[0].RA[0] > 360.0) { patch[0].RA[0] -= 360.0; } + while (patch[0].RA[0] < 0.0) { patch[0].RA[0] += 360.0; } + while (patch[0].RA[1] > 360.0) { patch[0].RA[1] -= 360.0; } + while (patch[0].RA[1] < 0.0) { patch[0].RA[1] += 360.0; } + if (patch[0].DEC[0] > patch[0].DEC[1]) { + dtmp = patch[0].DEC[0]; + patch[0].DEC[0] = patch[0].DEC[1]; + patch[0].DEC[1] = dtmp; + } + + /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!) */ + if (patch[0].DEC[1] > 86.25) { + area.RA[0] = 0.00; + area.RA[1] = 360.00; + area.DEC[0] = 86.25; + area.DEC[1] = 90.00; + sprintf (file, "n8230/pole"); + area.name = strcreate (file); + add_to_regions (&area); + free (area.name); + } + + /* if RA[0] > RA[1], split in two passes */ + if (patch[0].RA[0] > patch[0].RA[1]) { + area = patch[0]; + area.RA[0] = 0.0; + find_dec_bands (&area); + + area = patch[0]; + area.RA[1] = 360.0; + find_dec_bands (&area); + } else { + area = patch[0]; + find_dec_bands (&area); + } + + *nregion = Nregion; + return (region); +} + +/* given an area [RA[0] - RA[1], DEC[0] - DEC[1], + where 0.0 < RA[0] < RA[1] < 360.0 + and -90.0 < (DEC[0], DEC[1]) < 90.0 + and DEC[0] < DEC[1] + search through the dec bands to identify ones which overlap the area */ + +int find_dec_bands (CatStats *area) { + + int i; + + /* for NORTH region (DecBand[i] < DecBand[i+1]) */ + for (i = 0; i < NDECBANDS; i++) { + if ((area[0].DEC[1] > DecBand[i]) || (area[0].DEC[0] < DecBand[i+1])) { + load_ra_blocks (i, area); + } + } + + /* for SOUTH region (DecBand[i] > DecBand[i+1]) */ + for (i = NDECBANDS+1; i < 2*NDECBANDS; i++) { + if ((area[0].DEC[0] < DecBand[i]) || (area[0].DEC[1] > DecBand[i+1])) { + load_ra_blocks (i, area); + } + } + return (TRUE); +} + +/* given the number of a DecBlock and a defined area, return all area names in that block */ +/* area is defined so that RA[0] < RA[1], DEC[0] < DEC[1] */ +int load_ra_blocks (int Ndec, CatStats *area) { + + CatStats tregion; + int i, Nlines, Nskip; + char *buffer, file[MAX_NAME]; + FILE *f; + + /* count lines before Dec Band, and in Dec Band */ + Nskip = 0; + for (i = 0; i < Ndec; i++) Nskip += NDecLines[i]; + Nlines = NDecLines[Ndec]; + + ALLOCATE (buffer, char, 48*Nlines); + ALLOCATE (tregion.name, char, MAX_NAME); + + /* get data from GSC file */ + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE); + exit (1); + } + fseek (f, 5*2880 + 48*Nskip, SEEK_SET); + fread (buffer, 48*Nlines, 1, f); + fclose (f); + + /* find overlapping regions */ + for (i = 0; i < Nlines; i++) { + parse_GSC_line (&tregion, &buffer[i*48]); + if (area[0].RA[0] > tregion.RA[1]) continue; + if (area[0].RA[1] < tregion.RA[0]) continue; + if (area[0].DEC[0] > tregion.DEC[1]) continue; + if (area[0].DEC[1] < tregion.DEC[0]) continue; + sprintf (file, "%s/%s", Dec2Sections[Ndec], tregion.name); + strcpy (tregion.name, file); + add_to_regions (&tregion); + } + + free (buffer); + free (tregion.name); + return (TRUE); +} + +/* parse the GSC line, placing values in 'region'. force RA[0] < RA[1], DEC[0] < DEC[1] */ +int parse_GSC_line (CatStats *tregion, char *line) { + + double v0, s, dtmp; + char tmp[10]; + + /* copy name (4 chars) */ + bzero (tregion[0].name, 10); + strncpy (tregion[0].name, &line[1], 4); + + /* get RA[0] */ + bzero (tmp, 10); + strncpy (tmp, &line[7], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[10], 2); + v0 += atof (tmp) / 60.0; + + bzero (tmp, 10); + strncpy (tmp, &line[13], 5); + v0 += atof (tmp) / 3600.0; + tregion[0].RA[0] = v0 * 15.0; + + /* get RA[1] */ + bzero (tmp, 10); + strncpy (tmp, &line[19], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[22], 2); + v0 += atof (tmp) / 60.0; + + bzero (tmp, 10); + strncpy (tmp, &line[25], 5); + v0 += atof (tmp) / 3600.0; + tregion[0].RA[1] = v0 * 15.0; + + /* this only happens at the 0,360 boundary */ + if (tregion[0].RA[1] < tregion[0].RA[0]) tregion[0].RA[1] += 360.0; + + /* get DEC[0] */ + if (line[31] == '-') s = -1; + else s = +1; + + bzero (tmp, 10); + strncpy (tmp, &line[32], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[35], 2); + v0 += atof (tmp) / 60.0; + tregion[0].DEC[0] = s * v0; + + /* get DEC[1] */ + if (line[40] == '-') s = -1; + else s = +1; + + bzero (tmp, 10); + strncpy (tmp, &line[41], 2); + v0 = atof (tmp); + + bzero (tmp, 10); + strncpy (tmp, &line[44], 2); + v0 += atof (tmp) / 60.0; + tregion[0].DEC[1] = s * v0; + + if (tregion[0].DEC[0] > tregion[0].DEC[1]) { + dtmp = tregion[0].DEC[1]; + tregion[0].DEC[1] = tregion[0].DEC[0]; + tregion[0].DEC[0] = dtmp; + } + return (TRUE); +} + +void init_regions () { + + Nregion = 0; + NREGION = 100; + + ALLOCATE (region, CatStats, NREGION); + +} + +void add_to_regions (CatStats *area) { + + region[Nregion] = area[0]; + ALLOCATE (region[Nregion].name, char, MAX_NAME); + sprintf (region[Nregion].name, "%s/%s.%s", CatDir, area[0].name, "cpt"); + + area_of_region (®ion[Nregion]); + + Nregion ++; + + if (Nregion == NREGION) { + NREGION += 100; + REALLOCATE (region, CatStats, NREGION); + } + +} + +void set_catalog (char *catdir) { + CatDir = catdir; +} Index: /tags/gastro-2-2/Ohana/src/gastro2/src/grid.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/grid.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/grid.c (revision 6697) @@ -0,0 +1,217 @@ +# include "gastro2.h" + +# define NPIX 50 +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; + + 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*NPIX; + XMAX -= 0.5*NPIX; + YMIN -= 0.5*NPIX; + YMAX -= 0.5*NPIX; + gridfree (); + } + + free (Ref.stars); + +} + +void gridinit (double XMIN, double XMAX, double YMIN, double YMAX, int Nr, int Nt) { + + NX = (XMAX - XMIN) / NPIX; + NY = (YMAX - YMIN) / NPIX; + + 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: /tags/gastro-2-2/Ohana/src/gastro2/src/gstars2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/gstars2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/gstars2.c (revision 6697) @@ -0,0 +1,176 @@ +# include "gastro2.h" + +void gstars (char *filename, CmpCatalog *Target) { + + int Nstars; + char line[80]; + double det; + int NX, NY, Nskip, FoundAstrom, extend; + 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 (!fits_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 (!fits_read_header (HEADER, &header)) { + fprintf (stderr, "ERROR: can't load external header\n"); + exit (1); + } + FoundAstrom = GetCoords (&Target[0].coords, &header); + fits_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 &= fits_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 &= fits_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 */ + extend = FALSE; + fits_scan (&Target[0].header, "EXTEND", "%t", 1, &extend); + if (extend) { + Nskip = fits_matrix_size (&Target[0].header); + fseek (f, Nskip, SEEK_CUR); + stars = rfits (f, &Nstars); + } else { + /* allocate space for stars */ + if (!fits_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: /tags/gastro-2-2/Ohana/src/gastro2/src/lumfunc.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/lumfunc.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/lumfunc.c (revision 6697) @@ -0,0 +1,110 @@ +# 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; + 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] ++; + } + + /* 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++; + } + } + + /* 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]; + } + gaussj (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: /tags/gastro-2-2/Ohana/src/gastro2/src/misc.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/misc.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/misc.c (revision 6697) @@ -0,0 +1,140 @@ +# 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); + } + } +} + +/***** convert [-]00:00:00 to 0.0000 ****/ +int dms_to_ddd (double *Value, char *string) { + + int valid, neg, status; + double tmp, value; + char *p1, *p2, *px; + + valid = FALSE; + neg = FALSE; + stripwhite (string); + p1 = string; + px = string + strlen(string); + + if (string[0] == '-') { + valid = TRUE; + neg = TRUE; + p1 = &string[1]; + } + if (string[0] == '+') { + valid = TRUE; + neg = FALSE; + p1 = &string[1]; + } + if (isdigit(string[0])) { + valid = TRUE; + p1 = &string[0]; + } + if (!valid) { return (FALSE); } + + status = 1; + tmp = strtod (p1, &p2); + if (p2 == p1) return (FALSE); /* entry not a number: +fred */ + value = tmp; + if (p2 == px) goto escape; /* entry only number: +1.0 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:fred */ + status = 2; + value += tmp / 60.0; + if (p2 == px) goto escape; /* entry only number: +1:1 */ + p1 = p2 + 1; + + tmp = strtod (p1, &p2); + if (p2 == p1) goto escape; /* entry not a number: +1:1:fred */ + value += tmp / 3600.0; + + escape: + if (neg) { + value *= -1; + } + *Value = value; + + return (status); +} + +/**********/ +int str_to_radec (double *ra, double *dec, char *str1, char *str2) { + + double Ra, Dec; + + *ra = *dec = 0; + switch (dms_to_ddd (&Ra, str1)) { + case 0: + fprintf (stderr, "syntax error in RA\n"); + return (FALSE); + case 1: + break; + case 2: + Ra = Ra * 15; + break; + } + switch (dms_to_ddd (&Dec, str2)) { + case 0: + fprintf (stderr, "syntax error in DEC\n"); + return (FALSE); + case 1: + case 2: + break; + } + *ra = Ra; + *dec = Dec; + return (TRUE); +} + +void area_of_region (CatStats *region) { + + double area; + + area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG)); + region[0].Area = area; +} + +void 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: /tags/gastro-2-2/Ohana/src/gastro2/src/plots.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/plots.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/plots.c (revision 6697) @@ -0,0 +1,360 @@ +# 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; + 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); + for (i = 0; i < Nvect; i++) { + xvect[i] = Ref[0].stars[i].X; + yvect[i] = Ref[0].stars[i].Y; + M = (Target[0].lum.Mmax - 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); + + ymin = 5; ymax = -5; + for (i = 0; i < Nt; i++) { + tval[i] = tval[i] - log (Target[0].Area) / log (10.0); + ymin = MIN (tval[i], ymin); + ymax = MAX (tval[i], ymax); + } + + for (i = 0; i < Nr; i++) { + rval[i] = rval[i] - log (Ref[0].Area) / log (10.0); + 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 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: /tags/gastro-2-2/Ohana/src/gastro2/src/plotstuff.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/plotstuff.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/plotstuff.c (revision 6697) @@ -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 &", socket_name); +# 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: /tags/gastro-2-2/Ohana/src/gastro2/src/remove_clumps.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/remove_clumps.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/remove_clumps.c (revision 6697) @@ -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: /tags/gastro-2-2/Ohana/src/gastro2/src/rfits.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/rfits.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/rfits.c (revision 6697) @@ -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 (!fits_fread_ftable (f, &table, "SMPFILE")) goto escape; + + stars = fits_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: /tags/gastro-2-2/Ohana/src/gastro2/src/rotate2.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/rotate2.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/rotate2.c (revision 6697) @@ -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: /tags/gastro-2-2/Ohana/src/gastro2/src/rtext.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/rtext.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/rtext.c (revision 6697) @@ -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: /tags/gastro-2-2/Ohana/src/gastro2/src/sort.c =================================================================== --- /tags/gastro-2-2/Ohana/src/gastro2/src/sort.c (revision 6697) +++ /tags/gastro-2-2/Ohana/src/gastro2/src/sort.c (revision 6697) @@ -0,0 +1,195 @@ +# include "gastro2.h" + +void sort_lum (double *R, double *X, double *Y, int N) { + + int l,j,ir,i; + double tX, tY, tR; + + 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; + + 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; + + 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; + + 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; + + 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; + } +} + +