Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/Makefile =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/Makefile (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/Makefile (revision 21724) @@ -0,0 +1,91 @@ +default: gastro +help: + @echo "make options: gastro (default)" + +include ../../Configure + +HOME = $(ROOT)/src/gastro +BIN = $(HOME)/bin +INC = $(HOME)/include +SRC = $(HOME)/src +MAN = $(HOME)/doc +DESTBIN = $(LBIN) +DESTLIB = $(LLIB) +DESTINC = $(LINC) +DESTMAN = $(LMAN) + +# +INCS = -I$(INC) -I$(LINC) -I$(XINC) +LIBS = -L$(LLIB) -ldvo -lFITS -lohana -lm -lsocket -lnsl +CFLAGS = $(INCS) +LFLAGS = $(LIBS) + +# ----------------------- +OBJ = \ +$(SRC)/gastro.$(ARCH).o \ +$(SRC)/gstars.$(ARCH).o \ +$(SRC)/gargs.$(ARCH).o \ +$(SRC)/gcenter.$(ARCH).o \ +$(SRC)/granges.$(ARCH).o \ +$(SRC)/gfit.$(ARCH).o \ +$(SRC)/rotate.$(ARCH).o \ +$(SRC)/gproject.$(ARCH).o \ +$(SRC)/gheader.$(ARCH).o \ +$(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o \ +$(SRC)/gfitpoly.$(ARCH).o \ +$(SRC)/get_region_coords.$(ARCH).o \ +$(SRC)/greference.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/getgsc.$(ARCH).o \ +$(SRC)/getptolemy.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/misc.$(ARCH).o + +#$(SRC)/gregions2.$(ARCH).o \ +#$(SRC)/gcatalog.$(ARCH).o \ +#$(SRC)/gptolemy.$(ARCH).o \ +#$(SRC)/g2mass.$(ARCH).o \ + +gastro: $(BIN)/gastro.$(ARCH) + +$(BIN)/gastro.$(ARCH): $(OBJ) + +$(OBJ): $(INC)/gastro.h + +INSTALL = gastro + +# generic dependancy rules for binary code ######################### +.PRECIOUS: %.$(ARCH).o +.PRECIOUS: $(BIN)/%.$(ARCH) + +%.$(ARCH).o : %.c + $(CC) $(CFLAGS) -c $< -o $@ + +$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o + @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi + $(CC) $^ -o $@ $(LFLAGS) + +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* + +$(INSTALL) $(DEVEL): % : $(BIN)/%.$(ARCH) + +%.clean : + rm -f $(BIN)/$*.$(ARCH) + +%.install: + make $(DESTBIN)/$* + +# utilities ################################################# + +install: + for i in $(INSTALL); do make $$i.install; done + +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/bin/.cvsignore =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/bin/.cvsignore (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/bin/.cvsignore (revision 21724) @@ -0,0 +1,2 @@ +*.linux *.lin64 *.sol *.sun *.sid *.hp *.irix +*.linrh Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/doc/ChangeLog.txt =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/doc/ChangeLog.txt (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/doc/ChangeLog.txt (revision 21724) @@ -0,0 +1,13 @@ + + gastro-1-4: + * converted to gfits APIs (forces libfits 1.6) + * converted to new DVO APIs (forces libdvo 1.3) + * removed unused code + * better error-checking for ConfigInit + + gastro-1-3: + minor fixes to use updates to average.d + + gastro-1-2: + added support for dvo table modes / formats (libdvo v1.0) + minor changes to sync with libfits modes (libfits v1.4) Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/include/gastro.h =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/include/gastro.h (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/include/gastro.h (revision 21724) @@ -0,0 +1,120 @@ +# include +# include + +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 MMIN; +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 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 FORCE; +double F_RA; +double F_DEC; + +char GSCFILE[256], GSCDIR[256], LONEOS_REGION_FILE[256]; +double ASEC_PIX; +char ROUGH_ASTROMETRY[64]; + +typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, color; + double lweight, size; +} Graphdata; + +/* simple structure to carry around data on an array of stars */ +typedef struct { + double X; + double Y; + double mag; +} SStars; + +typedef struct { + Coords coords; + float *X, *Y; + int *N; + double RA[2], DEC[2]; + double Area, density, spacing; +} CatStats; + +typedef struct { + double R, D; + double r, b; +} USNOdata; + +typedef struct { + int *match; + float *X, *Y; + int *N; +} USNOstats; + +/* this seems to be a problem: is not included from math.h with the -ansi flag */ +extern double hypot PROTO((double, double)); + +SStars *getptolemy PROTO((CatStats *catstats, int *NSTARS)); +SStars *getgsc PROTO((CatStats *catstats, int *NSTARS)); +USNOdata *getusno PROTO((USNOstats *usnostats, CatStats *catstats, int *Nusno)); +SStars *gstars PROTO((char *file, int *NSTARS, Coords *coords, int *NX, int *NY, double *dNdM)); +void ConfigInit PROTO((int *argc, char **argv)); +void DonePlotting PROTO((Graphdata *graphmode, int N)); +void PlotReset PROTO((int N)); +void PlotVector PROTO((int Npts, float *vect, int mode, int N)); +void PrepPlotting PROTO((int Npts, Graphdata *graphmode, int N)); +void XDead PROTO(()); +void alter_header PROTO((char *, char **, int, double, double, double, double, double, double, double, double, int)); +void area_of_region PROTO((CatStats *region)); +void define_region PROTO((CatStats *catstats, Coords *coords, int NX, int NY)); +void find_shift PROTO((SStars *, SStars *, int, int, double, double, int *, double *, double *, double *)); +void gargs PROTO((int *argc, char **argv, Coords *coords)); +int gaussj PROTO((double **a, int n, double **b, int m)); +int gcenter PROTO((SStars *stars1in, SStars *stars2, int N1, int N2, Coords *coords, int NX, int NY, double *dR)); +int get_region_coords PROTO((double *ra, double *dec, int rnumber, char *side)); +int gfit PROTO((SStars *stars1, SStars *stars2, int N1, int N2, Coords *coords, int NX, int NY, double *Radius, double *DR, int *Nmatch, int mode)); +void gfitpoly PROTO((SStars *stars1, SStars *stars2, int N1, int N2, Coords *coords, double *Radius, double *DR, int *Nmatch)); +void gheader PROTO((char *file, Coords coords, double dR, int Nmatch)); +void gproject PROTO((SStars *catalog, SStars **stars, int Ncat, int *Nstars, Coords *coords, int NX, int NY, double dNdM, int N1)); +void granges PROTO((SStars *stars1, SStars *stars2, int N1, int N2, int NPIX, double *gx, double *gy, double *gx0, double *gy0)); +int greference PROTO((SStars **cat, int *Ncat, Coords *coords, int NX, int NY)); +void hh_hms PROTO((double hh, int *hr, int *mn, double *sc)); +void hms_format PROTO((char *line, double value)); +int line_fit PROTO((SStars *, SStars *, int, int, double, double, double, double, double *, double *, double *, double *, double *, double *, double *, double *)); +int mk_polyterm PROTO((int n, int m, int norder)); +int mk_vector PROTO((int n, int m, int norder)); +int open_graph PROTO((int N)); +void precess PROTO((double *, double *, double, double)); +void ranges PROTO((SStars *, SStars *, int, int, double *, double *, double *, double *, double, double)); +void rotate PROTO((SStars *stars, int Nstars, double angle, int Xo, int Yo)); +void sort_stars PROTO((SStars *stars, int N)); +void sort_lists PROTO((double *X, double *Y, int *S, int N)); +void stats PROTO((char *, double *, double *, double *, double *, double)); Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/ConfigInit.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/ConfigInit.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/ConfigInit.c (revision 21724) @@ -0,0 +1,55 @@ +# include "gastro.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); // max allowed offset in gcenter + ScanConfig (config, "MIN_MATCHES", "%d", 0, &MIN_MATCHES); // min allowed fitted stars + ScanConfig (config, "DEFAULT_RADIUS", "%lf", 0, &DEFAULT_RADIUS); // starting radius for matched fit + ScanConfig (config, "MINIMUM_RADIUS", "%lf", 0, &MINIMUM_RADIUS); // min allowed radius for matched fit + ScanConfig (config, "MAX_ERROR", "%lf", 0, &MAX_ERROR); // max allowed error for valid solution + ScanConfig (config, "MAX_NONLINEAR", "%lf", 0, &MAX_NONLINEAR); // max allowed shear |(11*22 - 12*21) / (L1*L2)| + ScanConfig (config, "CCD_PC1_1", "%lf", 0, &CCD_PC1_1); // guess if WCS is missing + ScanConfig (config, "CCD_PC2_2", "%lf", 0, &CCD_PC2_2); // guess if WCS is missing + ScanConfig (config, "CCD_PC1_2", "%lf", 0, &CCD_PC1_2); // guess if WCS is missing + ScanConfig (config, "CCD_PC2_1", "%lf", 0, &CCD_PC2_1); // guess if WCS is missing + ScanConfig (config, "ASEC_PIX", "%lf", 0, &ASEC_PIX); // guess if WCS is missing + ScanConfig (config, "NFIELD", "%lf", 0, &NFIELD); // search region in field units + ScanConfig (config, "NPOLYTERMS", "%d", 0, &NPOLYTERMS); // fit order + ScanConfig (config, "ROT_ZERO", "%lf", 0, &ROT_ZERO); // rotation search region + ScanConfig (config, "dROT", "%lf", 0, &dROT); // rotation search region + ScanConfig (config, "NROT", "%d", 0, &NROT); // rotation search region + ScanConfig (config, "POLAR_ALIGNMENT", "%d", 0, &POLAR_ALIGNMENT); // apply polar alignment correction + ScanConfig (config, "POLAR_AXIS_RA", "%lf", 0, &POLE_RA); // true coords of pole (should be HA, not RA) + ScanConfig (config, "POLAR_AXIS_DEC", "%lf", 0, &POLE_DEC); // true coords of pole + ScanConfig (config, "RA_OFFSET", "%lf", 0, &RA_OFFSET); // ?? not well defined (should be euler angle) + ScanConfig (config, "DEC_OFFSET", "%lf", 0, &DEC_OFFSET); // ?? not well defined (should be euler angle) + ScanConfig (config, "LONEOS_REGIONS", "%s", 0, LONEOS_REGION_FILE); // table of LONEOS regions to fix guess astrometry + + ScanConfig (config, "GSCFILE", "%s", 0, GSCFILE); // location of sky table + ScanConfig (config, "GSCDIR", "%s", 0, GSCDIR); // location of HST GSC data + ScanConfig (config, "USNO_CDROM", "%s", 0, CDROM); // location of USNO A data (USNO_A_DIR in gastro2) + ScanConfig (config, "ASTRO_REFCAT", "%s", 0, REFCAT); // which astrometry catalog to use + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); // location of ptolemy-format ref data + ScanConfig (config, "ROUGH_ASTROMETRY", "%s", 0, ROUGH_ASTROMETRY); // where to get initial guess (header, config) + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); // location of photcode table to convert supplied photcode + + 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/ohana/gastro-1-4/Ohana/src/gastro/src/gargs.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gargs.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gargs.c (revision 21724) @@ -0,0 +1,125 @@ +# include "gastro.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 gargs (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/ohana/gastro-1-4/Ohana/src/gastro/src/gastro.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gastro.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gastro.c (revision 21724) @@ -0,0 +1,85 @@ +# include "gastro.h" +# define ABORT \ +{ Nmatch = 1; dR = 0; \ + gheader (argv[1], coords, dR, Nmatch); \ + exit (0); } + +int main (int argc, char **argv) { + + int N1, N2, Ncat, NX, NY, Nmatch, Nminterms; + SStars *catalog, *stars1, *stars2; + struct timeval now, then; + Coords coords; + double dNdM, dR, Radius; + + gettimeofday (&then, (void *) NULL); + + ConfigInit (&argc, argv); + gargs (&argc, argv, &coords); + + /* load stars from image (*.cmp file) */ + N1 = NMAX_STARS; /* we only want a small number of stars */ + stars1 = gstars (argv[1], &N1, &coords, &NX, &NY, &dNdM); + + /* load stars from reference catalogs */ + greference (&catalog, &Ncat, &coords, NX, NY); + + /* get rough alignment with reference stars */ + gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY, dNdM, N1); + gcenter (stars1, stars2, N1, N2, &coords, NX, NY, &Radius); + free (stars2); + + /* reload reference, get good astrometry */ + if (!greference (&catalog, &Ncat, &coords, NX, NY)) ABORT; + /* NFIELD = 1.2; */ + gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY, dNdM, N1); + if (!gfit (stars1, stars2, N1, N2, &coords, NX, NY, &Radius, &dR, &Nmatch, 1)) ABORT; + free (stars2); + + gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY, dNdM, N1); + gfitpoly (stars1, stars2, N1, N2, &coords, &Radius, &dR, &Nmatch); + + if (dR > MAX_ERROR) { + fprintf (stderr, "ERROR: bad solution! %f %f (%d stars)\n", dR, (dR / sqrt(1.0*Nmatch)), Nmatch); + ABORT; + } + fprintf (stderr, "good solution: %f %f (%d stars)\n", dR, (dR / sqrt(1.0*Nmatch)), Nmatch); + + Nminterms = MIN_MATCHES; + switch (NPOLYTERMS) { + case 0: + case 1: + /* we don't really allow zero order fits */ + Nminterms = MIN_MATCHES; + break; + case 2: + Nminterms = 20; + break; + case 3: + Nminterms = 40; + break; + } + + if (Nmatch <= Nminterms) { + fprintf (stderr, "ERROR: too few stars for reliable solution, only %d\n", Nmatch); + ABORT; + } + + gheader (argv[1], coords, dR, Nmatch); + + if (VERBOSE) { + gettimeofday (&now, (void *) NULL); + fprintf (stderr, "%s: elapsed time = %.2f sec\n", argv[1], + (now.tv_sec - then.tv_sec) + 1e-6*(now.tv_usec - then.tv_usec)); + } + fprintf (stderr, "SUCCESS\n"); + exit (0); +} + +/* + free (stars2); + + DEFAULT_RADIUS = DEFAULT_RADIUS / 4; + gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY, dNdM, N1); + if (!gfit (stars1, stars2, N1, N2, &coords, NX, NY, &Radius, &dR, &Nmatch, 0)) ABORT; +*/ Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gcenter.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gcenter.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gcenter.c (revision 21724) @@ -0,0 +1,255 @@ +# include "gastro.h" +# define NSIGMA 2.0 + +/* stars1.X,Y and stars2.X,Y are in image pixels + stars2 (ref catalog) is approx, based on guess for scale */ + +int gcenter (SStars *stars1in, SStars *stars2, int N1, int N2, Coords *coords, int NX, int NY, double *dR) { + + double mean, sigma, gx, gy, gx0, gy0, n, SearchRadius; + int NPIX, minN, Nmin, Nmin0; + int i, j, k; + double *N, *DX, *DY, *D2, *tX1, *tY1, *tX2, *tY2; + double rot, Rot, Smin, Smin0, s, Fmin, Fmin0, f; + double Xmin, Xmin0, Ymin, Ymin0; + double RA, DEC, RAo, DECo; + double dX, dY, refX, refY; + double cs, sn; + double dR1, dD1, dR2, dD2, d1, d2; + double *sx1, *sy1, *sx2, *sy2; + SStars *stars1; + char c; + Graphdata graphdata; + float *xvect, *yvect; + int Nvect, NVECT; + + SearchRadius = NX * SEARCH_RADIUS; + xvect = yvect = NULL; + Nvect = 0; + + if (PLOTSTUFF) { + NVECT = N1*N2; + ALLOCATE (xvect, float, NVECT); + ALLOCATE (yvect, float, NVECT); + graphdata.xmin = 0; + graphdata.xmax = 2000; + graphdata.ymin = 0; + graphdata.ymax = 4000; + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 0.5; + } + + /* NPIX = MAX (300, sqrt (20*(N1*N2) / MIN (N1, N2))); */ + /* NPIX = 10 * sqrt (N2); */ + NPIX = 300; + mean = 2 * N1*N2 / (NPIX*NPIX); + sigma = sqrt (mean); + minN = MAX (6, mean + NSIGMA*sigma); + fprintf (stderr, "N1: %d, N2: %d, minN: %d\n", N1, N2, minN); + + ALLOCATE (N, double, NPIX*NPIX); + ALLOCATE (DX, double, NPIX*NPIX); + ALLOCATE (DY, double, NPIX*NPIX); + ALLOCATE (D2, double, NPIX*NPIX); + + ALLOCATE (stars1, SStars, N1); + for (i = 0; i < N1; i++) { + stars1[i] = stars1in[i]; + } + + if (PLOTSTUFF) { + for (i = 0; i < N1; i++) { + xvect[Nvect] = stars1in[i].X; + yvect[Nvect] = stars1in[i].Y; + Nvect ++; + } + } + if (PLOTSTUFF) { + PlotReset (0); + PrepPlotting (Nvect, &graphdata, 0); + PlotVector (Nvect, xvect, 0, 0); + PlotVector (Nvect, yvect, 1, 0); + DonePlotting (&graphdata, 0); + usleep (300000); + fprintf (stderr, "plotting %d points\n", Nvect); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); + Nvect = 0; + } + + Xmin = Ymin = 0; + Fmin = Smin = 1000000.0; + rot = Rot = ROT_ZERO-dROT*NROT; + rotate (stars1, N1, Rot, (int)coords[0].crpix1, (int)coords[0].crpix2); + + graphdata.xmin = -200; + graphdata.xmax = 200; + graphdata.ymin = -200; + graphdata.ymax = 200; + + ALLOCATE (tX1, double, N1); + ALLOCATE (tY1, double, N1); + ALLOCATE (tX2, double, N2); + ALLOCATE (tY2, double, N2); + for (i = 0; i < N2; i++) { + tX2[i] = stars2[i].X; + tY2[i] = stars2[i].Y; + } + for (Rot = ROT_ZERO-dROT*NROT; Rot <= ROT_ZERO+dROT*NROT; Rot += dROT) { + granges (stars1, stars2, N1, N2, NPIX, &gx, &gy, &gx0, &gy0); + bzero (N, NPIX*NPIX*sizeof(double)); + bzero (DX, NPIX*NPIX*sizeof(double)); + bzero (DY, NPIX*NPIX*sizeof(double)); + bzero (D2, NPIX*NPIX*sizeof(double)); + for (i = 0; i < N1; i++) { + tX1[i] = stars1[i].X; + tY1[i] = stars1[i].Y; + } + sx1 = tX1; sy1 = tY1; + for (i = 0; i < N1; i++, sx1++, sy1++) { + sx2 = tX2; sy2 = tY2; + for (j = 0; j < N2; j++, sx2++, sy2++) { + dX = *sx1 - *sx2; + dY = *sy1 - *sy2; + if (hypot (dX, dY) > SearchRadius) continue; + if (PLOTSTUFF) { + xvect[Nvect] = stars1[i].X - stars2[j].X; + yvect[Nvect] = stars1[i].Y - stars2[j].Y; + Nvect ++; + } + k = NPIX*(int)(gx*dX+gx0) + (int)(gy*dY+gy0); + N[k] += 1.0; + DX[k] += dX; + DY[k] += dY; + D2[k] += dX*dX + dY*dY; + } + } + + Fmin0 = 1000.0; + Smin0 = Nmin0 = Xmin0 = Ymin0 = 0; + for (k = 0; k < NPIX*NPIX; k++) { + n = N[k]; /* 1*/ + if (n < minN) + continue; + s = D2[k] - (SQ(DX[k]) + SQ(DY[k])) / n; + f = s / (pow(n,4.0)); /* = 12 ops */ + + if (f < Fmin0) { + Fmin0 = f; + Smin0 = s; + Nmin0 = n; + Xmin0 = DX[k] / n; + Ymin0 = DY[k] / n; + } + } + if (VERBOSE) fprintf (stderr, "best offset: %7.1f %7.1f at %.1f deg (%f %f %d)\n", Xmin0, Ymin0, Rot, Fmin0, sqrt(Smin0), Nmin0); + if (Fmin0 < Fmin) { + Fmin = Fmin0; + Smin = Smin0; + Nmin = Nmin0; + Xmin = Xmin0; + Ymin = Ymin0; + rot = Rot; + } + if (PLOTSTUFF) { + PlotReset (0); + PrepPlotting (Nvect, &graphdata, 0); + PlotVector (Nvect, xvect, 0, 0); + PlotVector (Nvect, yvect, 1, 0); + DonePlotting (&graphdata, 0); + usleep (300000); + fprintf (stderr, "plotting %d points\n", Nvect); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); + Nvect = 0; + } + + rotate (stars1, N1, dROT, (int)coords[0].crpix1, (int)coords[0].crpix2); + } + rotate (stars1, N1, -ROT_ZERO-dROT*NROT, (int)coords[0].crpix1, (int)coords[0].crpix2); + + free (N); + free (DX); + free (DY); + free (D2); + + /* dx = dy = 10 pix */ + refX = coords[0].crpix1 - Xmin; + refY = coords[0].crpix2 - Ymin; + XY_to_RD (&RA, &DEC, refX, refY, coords); + XY_to_RD (&RAo, &DECo, (refX + 10), refY, coords); + dR1 = (RAo - RA)*cos(DEC*RAD_DEG); + dD1 = (DECo - DEC); + XY_to_RD (&RAo, &DECo, refX, (refY + 10), coords); + dR2 = (RAo - RA)*cos(DEC*RAD_DEG); + dD2 = (DECo - DEC); + d1 = coords[0].cdelt1; d2 = coords[0].cdelt2; + cs = cos(RAD_DEG*rot); sn = sin(RAD_DEG*rot); + + coords[0].pc1_1 = cs*dR1 / (10*d1) + sn*dR2 / (10*d1); coords[0].pc1_2 = cs*dR2 / (10*d1) - sn*dR1 / (10*d1); + coords[0].pc2_1 = cs*dD1 / (10*d2) + sn*dD2 / (10*d2); coords[0].pc2_2 = cs*dD2 / (10*d2) - sn*dD1 / (10*d2); + coords[0].crval1 = RA; + coords[0].crval2 = DEC; + + /* diameter of 1 pixel box */ + *dR = 1*sqrt(Smin); + + fprintf (stderr, "%f x %f, %f\n", 1/gx, 1/gy, *dR); + if (VERBOSE) fprintf (stderr, "using: %7.1f %7.1f at %.1f deg\n", Xmin, Ymin, rot); + +# if (0) + if (VERBOSE) { + fprintf (stderr, "%s\n", coords[0].ctype); + fprintf (stderr, "%f %f\n", coords[0].crval1, coords[0].crval2); + fprintf (stderr, "%f %f\n", coords[0].crpix1, coords[0].crpix2); + fprintf (stderr, "%f %f\n", coords[0].pc1_1, coords[0].pc1_2); + fprintf (stderr, "%f %f\n", coords[0].pc2_1, coords[0].pc2_2); + fprintf (stderr, "%f %f\n", coords[0].cdelt1, coords[0].cdelt2); + } +# endif + + DEFAULT_RADIUS = MAX ((5*NX / NPIX), DEFAULT_RADIUS); + return (TRUE); +} + + /* we now have Xmin, Ymin, rot, get coords in unrotate stars1 frame of correct crref + RAo = coords[0].crval1; + DECo = coords[0].crval2; + + cs = cos(RAD_DEG*rot); sn = sin(RAD_DEG*rot); + refX = coords[0].crpix1; + refY = coords[0].crpix2; + + coords[0].crpix1 = Xmin*cs + Ymin*sn + refX; + coords[0].crpix2 = -Xmin*sn + Ymin*cs + refY; + + xx = coords[0].pc1_1; xy = coords[0].pc1_2; + yx = coords[0].pc2_1; yy = coords[0].pc2_2; + coords[0].pc1_1 = cs*xx + sn*xy; coords[0].pc1_2 = cs*xy - sn*xx; + coords[0].pc2_1 = cs*yx + sn*yy; coords[0].pc2_2 = cs*yy - sn*yx; + + XY_to_RD (&RA, &DEC, refX, refY, coords); + if (fabs(RAo - RA) > 90) { + RA = (RA > 180.0) ? (RA - 180) : (RA + 180); + DEC = (DEC > 0.0) ? (180.0 - DEC) : (-180.0 - DEC); + } + coords[0].crval1 = RA; + coords[0].crval2 = DEC; + coords[0].crpix1 = refX; + coords[0].crpix2 = refY; + + Rot = RAo - RA; + cs = cos(RAD_DEG*rot); sn = sin(RAD_DEG*rot); + xx = coords[0].pc1_1; xy = coords[0].pc1_2; + yx = coords[0].pc2_1; yy = coords[0].pc2_2; + coords[0].pc1_1 = cs*xx + sn*xy; coords[0].pc1_2 = cs*xy - sn*xx; + coords[0].pc2_1 = cs*yx + sn*yy; coords[0].pc2_2 = cs*yy - sn*yx; + + + +*/ Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/get_region_coords.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/get_region_coords.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/get_region_coords.c (revision 21724) @@ -0,0 +1,57 @@ +# include "gastro.h" +# define NBYTE_LINE 53 +# define NLINES 100 + +int get_region_coords (double *ra, double *dec, int rnumber, char *side) { + + FILE *f; + int i, j, done, found, Nbytes, Nline, num, NBYTES; + char *buffer; + double R, D; + + f = fopen (LONEOS_REGION_FILE, "r"); + if (f == NULL) { + fprintf (stderr, "couldn't find region map %s\n", LONEOS_REGION_FILE); + return (FALSE); + } + + NBYTES = NBYTE_LINE * NLINES; + ALLOCATE (buffer, char, NBYTES); + + found = done = FALSE; + for (i = 0; !done && !found; i++) { + Nbytes = fread (buffer, sizeof(char), NBYTES, f); + if (Nbytes < 1) done = TRUE; + Nline = Nbytes / NBYTE_LINE; + for (j = 0; !found && (j < Nline); j++) { + num = atof (&buffer[j*NBYTE_LINE]); + if (num == rnumber) { + found = TRUE; + sscanf (&buffer[j*NBYTE_LINE], "%*d %lf %lf", &R, &D); + fwrite (&buffer[j*NBYTE_LINE], 1, 106, stderr); + fprintf (stderr, "\n\n%f %f\n", R, D); + if (!strncasecmp (side, "east", 4)) { + R += 0.026 / cos (D); + /* if the word says "east", we need to offset by 1 chip width, + R and D are in radians here, so 0.026 is 1.5 deg in radians */ + } + R *= (180.0 / M_PI); + D *= (180.0 / M_PI); + } + } + } + + free (buffer); + fclose (f); + + if (!found) { + fprintf (stderr, "error: can't find desired region number %d\n", rnumber); + *ra = *dec = 0; + return (FALSE); + } + + *ra = R; + *dec = D; + return (TRUE); + +} Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getgsc.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getgsc.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getgsc.c (revision 21724) @@ -0,0 +1,86 @@ +# include "gastro.h" +# define BYTES_STAR 23 +# define BLOCK 1000 + +SStars *rd_gsc (char *filename, int *Nstars); + +SStars *getgsc (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Ngsc, Nstars; + SStars *gsc; + SStars *stars; + SkyList *skylist; + SkyTable *sky; + SkyRegion patch; + + patch.Rmin = catstats[0].RA[0]; + patch.Rmax = catstats[0].RA[1]; + patch.Dmin = catstats[0].DEC[0]; + patch.Dmax = catstats[0].DEC[1]; + + /* load regions from GSC table, restrict to patch */ + sky = SkyTableFromGSC (GSCFILE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (sky, GSCDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &patch); + + Nstars = 0; + ALLOCATE (stars, SStars, 1); + + for (i = 0; i < skylist[0].Nregions; i++) { + gsc = rd_gsc (skylist[0].filename[i], &Ngsc); + + Ns = Nstars; + Nstars += Ngsc; + + REALLOCATE (stars, SStars, MAX (1, Nstars)); + for (k = Ns, j = 0; j < Ngsc; k++, j++) { + stars[k].X = gsc[j].X; + stars[k].Y = gsc[j].Y; + stars[k].mag = gsc[j].mag; + } + free (gsc); + } + SkyTableFree (sky); + + *NSTARS = Nstars; + return (stars); +} + +SStars *rd_gsc (char *filename, int *Nstars) { + + SStars *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, SStars, 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].X, 1, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].Y, 2, &buffer[i*BYTES_STAR]); + dparse (&stars[nstar].mag, 3, &buffer[i*BYTES_STAR]); + nstar++; + if (nstar == NSTAR) { + NSTAR += 1000; + REALLOCATE (stars, SStars, NSTAR); + } + } + } + + free (buffer); + + *Nstars = nstar; + return (stars); +} Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getptolemy.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getptolemy.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/getptolemy.c (revision 21724) @@ -0,0 +1,60 @@ +# include "gastro.h" + +SStars *getptolemy (CatStats *catstats, int *NSTARS) { + + int i, j, k, Ns, Nstars; + Catalog catalog; + SStars *stars; + SkyList *skylist; + SkyTable *sky; + SkyRegion patch; + + patch.Rmin = catstats[0].RA[0]; + patch.Rmax = catstats[0].RA[1]; + patch.Dmin = catstats[0].DEC[0]; + patch.Dmax = catstats[0].DEC[1]; + + Nstars = 0; + ALLOCATE (stars, SStars, 1); + + /* load regions from GSC table, restrict to patch */ + sky = SkyTableLoadOptimal (CATDIR, NULL, GSCFILE, SKY_DEPTH_HST, VERBOSE); + SkyTableSetFilenames (sky, CATDIR, "cpt"); + skylist = SkyListByPatch (sky, -1, &patch); + + for (i = 0; i (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 = 5000; + ALLOCATE (usno, USNOdata, NUSNO); + nusno = 0; + + 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, "still can't open file %s, is cdrom %d in drive?\n", filename, disk); + exit (0); + } + } + 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, "RA out of range\n"); + exit (0); + } + + /* 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, "can't open file %s\n", filename); + exit (0); + } + /* 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 reading data from file %s (%d, %d, %d)\n", + filename, start[bin], number[bin], nitems); + exit (0); + } + 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)) { + usno[nusno].R = buf[0]/360000.0; + usno[nusno].D = buf[1]/360000.0 - 90.0; + usno[nusno].r = 0.1*(buf[2] - 1000*((int)(buf[2]/1000))); + usno[nusno].b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); + nusno ++; + if (nusno == NUSNO - 1) { + NUSNO += 5000; + REALLOCATE (usno, USNOdata, NUSNO); + } + } + } + free (buffer); + } + fclose (f); + } + + REALLOCATE (usno, USNOdata, MAX (nusno, 1)); + *Nusno = nusno; + return (usno); + +} + + Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfit.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfit.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfit.c (revision 21724) @@ -0,0 +1,460 @@ +# include "gastro.h" + +/* stars1.X,Y and stars2.X,Y are both in image pixels. + The conversion terms (X_O, X_X, etc) make linear + corrections in pixel coordinates */ + +int gfit (SStars *stars1, SStars *stars2, int N1, int N2, Coords *coords, int NX, int NY, double *Radius, double *DR, int *Nmatch, int mode) { + + int i, j, iteration, Niter, last, halt, first_j, *Nextra, extras; + int *tmpN1, *tmpN2; + double X_O, X_X, X_Y, dX; + double Y_O, Y_X, Y_Y, dY; + double x, y, x2, y2, xy, N, R, wt; + double X, Y, Xx, Yy, Xy, Yx; + double Dx, Dy, DD, d2X, d2Y; + double RA, DEC, radius, radius2, fratio; + double *tmpX1, *tmpX2, *tmpY1, *tmpY2; + double tX1, tX2, tY1, tY2; + double dS; + + char c; + Graphdata graphdata; + float *xvect, *yvect; + float *xvect2, *yvect2; + int Nvect, NVECT; + int Nvect2, NVECT2; + + Nvect = NVECT = Nvect2 = NVECT2 = 0; + xvect = yvect = xvect2 = yvect2 = NULL; + + if (PLOTSTUFF) { + NVECT2 = MAX(N1, N2); + NVECT = MAX(N1, N2); + ALLOCATE (xvect, float, NVECT); + ALLOCATE (yvect, float, NVECT); + ALLOCATE (xvect2, float, NVECT2); + ALLOCATE (yvect2, float, NVECT2); + } + + + /* allocate space for star coords */ + ALLOCATE (tmpX1, double, N1); + ALLOCATE (tmpY1, double, N1); + ALLOCATE (tmpN1, int, N1); + + ALLOCATE (tmpX2, double, N2); + ALLOCATE (tmpY2, double, N2); + ALLOCATE (tmpN2, int, N2); + ALLOCATE (Nextra, int, N2); + bzero (Nextra, sizeof(int) * N2); + + dX = dY = N = 0; + + /* assign and sort list */ + for (i = 0; i < N1; i++) { + tmpX1[i] = stars1[i].X; + tmpY1[i] = stars1[i].Y; + tmpN1[i] = i; + } + if (N1 > 1) sort_lists (tmpX1, tmpY1, tmpN1, N1); + + + /* choose iteration ranges */ + fratio = 1.41421; + extras = halt = last = FALSE; + radius = *Radius; + Niter = 2 + log (radius/MINIMUM_RADIUS) / log (fratio); + + /* initial values for fit coeffs */ + X_X = 1; X_Y = 0; X_O = 0; + Y_X = 0; Y_Y = 1; Y_O = 0; + + for (iteration = 0; iteration < Niter; iteration ++) { + + if (iteration >= Niter - 1) { /* next loop is the last one */ + radius *= fratio; + last = TRUE; + } + + /* setup and define */ + radius2 = radius*radius; + dX = dY = d2X = d2Y = x = y = x2 = y2 = xy = X = Y = Xx = Xy = Yx = Yy = N = R = 0; + for (i = 0; i < N2; i++) { + tmpX2[i] = (X_O) + (X_X)*stars2[i].X + (X_Y)*stars2[i].Y; + tmpY2[i] = (Y_O) + (Y_X)*stars2[i].X + (Y_Y)*stars2[i].Y; + tmpN2[i] = i; + } + + if (PLOTSTUFF) { + for (i = 0; i < N1; i++) { + xvect2[i] = tmpX1[i]; + yvect2[i] = tmpY1[i]; + } + Nvect2 = N1; + + graphdata.xmin = 0; + graphdata.xmax = 2000; + graphdata.ymin = 4000; + graphdata.ymax = 0; + graphdata.style = 2; + graphdata.ptype = 3; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 1.5; + + PlotReset (1); + PrepPlotting (Nvect2, &graphdata, 1); + PlotVector (Nvect2, xvect2, 0, 1); + PlotVector (Nvect2, yvect2, 1, 1); + DonePlotting (&graphdata, 1); + Nvect2 = 0; + + for (i = 0; i < N2; i++) { + xvect2[i] = tmpX2[i]; + yvect2[i] = tmpY2[i]; + } + Nvect2 = N2; + + graphdata.ptype = 1; + + PrepPlotting (Nvect2, &graphdata, 1); + PlotVector (Nvect2, xvect2, 0, 1); + PlotVector (Nvect2, yvect2, 1, 1); + DonePlotting (&graphdata, 1); + Nvect2 = 0; + } + + if (N2 > 1) sort_lists (tmpX2, tmpY2, tmpN2, N2); + + /* find matched stars */ + for (i = j = 0; (i < N1) && (j < N2); ) { + tX1 = tmpX1[i]; + tX2 = tmpX2[j]; + Dx = tX1 - tX2; + if (Dx <= -2.0*radius) { + i++; + continue; + } + if (Dx >= 2.0*radius) { + j++; + continue; + } + + /**** possible improvement: find only the closest match + for stars that have more than one (save DD and i for each + cat star j */ + /* in the right range */ + first_j = j; + for (; (Dx > -2.0*radius) && (j < N2); j++) { + tY1 = tmpY1[i]; + tX2 = tmpX2[j]; + tY2 = tmpY2[j]; + Dx = tX1 - tX2; + Dy = tY1 - tY2; + DD = Dx*Dx + Dy*Dy; + /* stars matched */ + if (DD < radius2) { + if (PLOTSTUFF) { + xvect[Nvect] = Dx; + yvect[Nvect] = tmpY1[i]; + Nvect ++; + if (Nvect == NVECT) { + NVECT += 100; + REALLOCATE (xvect, float, NVECT); + REALLOCATE (yvect, float, NVECT); + } + xvect2[Nvect2] = tmpX1[i]; + yvect2[Nvect2] = tmpY1[i]; + Nvect2 ++; + if (Nvect2 == NVECT2 - 1) { + NVECT2 += 100; + REALLOCATE (xvect2, float, NVECT2); + REALLOCATE (yvect2, float, NVECT2); + } + } + /* wt = sqrt(sqrt(DD)); */ + wt = DD + 1; + dX += Dx; + dY += Dy; + d2X += Dx*Dx; + d2Y += Dy*Dy; + x += tX2/wt; + y += tY2/wt; + x2 += tX2*tX2/wt; + y2 += tY2*tY2/wt; + xy += tX2*tY2/wt; + X += tX1/wt; + Y += tY1/wt; + Xx += tX1*tX2/wt; + Xy += tX1*tY2/wt; + Yx += tY1*tX2/wt; + Yy += tY1*tY2/wt; + N += 1.0; + R += 1.0/wt; + if (last && MATCHDUMP && mode) { + XY_to_RD (&RA, &DEC, stars2[tmpN2[j]].X, stars2[tmpN2[j]].Y, coords); + fprintf (stdout, "%f %f %f %f %f\n", RA, DEC, stars2[tmpN2[j]].X, stars2[tmpN2[j]].Y, stars2[tmpN2[j]].mag); + } + } + } + j = first_j; + i++; + } + + if (PLOTSTUFF) { + + graphdata.xmin = 0; + graphdata.xmax = 2000; + graphdata.ymin = 4000; + graphdata.ymax = 0; + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 1.5; + + PrepPlotting (Nvect2, &graphdata, 1); + PlotVector (Nvect2, xvect2, 0, 1); + PlotVector (Nvect2, yvect2, 1, 1); + DonePlotting (&graphdata, 1); + + graphdata.xmin = -150; + graphdata.xmax = 150; + graphdata.ymin = 0; + graphdata.ymax = 4000; + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 1.5; + + PlotReset (0); + PrepPlotting (Nvect, &graphdata, 0); + PlotVector (Nvect, xvect, 0, 0); + PlotVector (Nvect, yvect, 1, 0); + DonePlotting (&graphdata, 0); + usleep (300000); + fprintf (stderr, "plotting %d points\n", Nvect); + fprintf (stderr, "type return to continue"); + fscanf (stdin, "%c", &c); + Nvect = 0; + } + + /* + if (extras == 1) { + extras = 2; + halt = last = FALSE; + iteration -= 2; + radius *= fratio*fratio; + } + */ + + if (MATCHDUMP && mode && last) exit (0); + + if (NOMATCHDUMP && mode && last) { + for (i = 0; i < N2; i++) { + XY_to_RD (&RA, &DEC, stars2[i].X, stars2[i].Y, coords); + fprintf (stdout, "%f %f %f %f %f\n", RA, DEC, stars2[i].X, stars2[i].Y, stars2[i].mag); + } + exit (0); + } + + /* calculate the fit parameters */ + if (!last) { + double XX1, XY1, YX1, YY1, XO1, YO1; + double XX0, XY0, YX0, YY0, XO0, YO0; + double **matrix, **vector; + int NR, NC; + + if (VERBOSE) fprintf (stderr, "radius: %f, No. of matched stars: %d\n", radius, (int) N); + + if (N < 3) { + fprintf (stderr, "ERROR: too few stars\n"); + X_O = X_X = X_Y = Y_O = Y_X = Y_Y = 0; + return (FALSE); + } + + NR = 3; NC = 2; + ALLOCATE (matrix, double *, NR); + ALLOCATE (vector, double *, NR); + for (i = 0; i < NR; i++) { + ALLOCATE (matrix[i], double, NR); + ALLOCATE (vector[i], double, NC); + bzero (vector[i], NC*sizeof(double)); + bzero (matrix[i], NR*sizeof(double)); + } + + matrix[0][0] = R; + matrix[0][1] = matrix[1][0] = x; + matrix[0][2] = matrix[2][0] = y; + matrix[1][2] = matrix[2][1] = xy; + matrix[1][1] = x2; + matrix[2][2] = y2; + + vector[0][0] = X; + vector[1][0] = Xx; + vector[2][0] = Xy; + + vector[0][1] = Y; + vector[1][1] = Yx; + vector[2][1] = Yy; + + dgaussj (matrix, NR, vector, NC); + + /* + Sx2 = x2 - x*x/N; + Sy2 = y2 - y*y/N; + Sxy = xy - x*y/N; + SXx = Xx - X*x/N; + SXy = Xy - X*y/N; + SYx = Yx - Y*x/N; + SYy = Yy - Y*y/N; + */ + + XX0 = X_X; XY0 = X_Y; XO0 = X_O; + YX0 = Y_X; YY0 = Y_Y; YO0 = Y_O; + + /* fit parameters relative to rotated frame */ + /* + XX1 = (SXx*Sy2 - SXy*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + XY1 = (SXy*Sx2 - SXx*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + XO1 = X/N - (XX1)*x/N - (XY1)*y/N; + + YX1 = (SYx*Sy2 - SYy*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + YY1 = (SYy*Sx2 - SYx*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + YO1 = Y/N - (YX1)*x/N - (YY1)*y/N; + */ + + XO1 = vector[0][0]; + XX1 = vector[1][0]; + XY1 = vector[2][0]; + + YO1 = vector[0][1]; + YX1 = vector[1][1]; + YY1 = vector[2][1]; + + fprintf (stderr, "%f %f %f\n", vector[0][0], vector[1][0], vector[2][0]); + fprintf (stderr, "%f %f %f\n", vector[0][1], vector[1][1], vector[2][1]); + + /* fit parameters relative to original frame */ + X_X = XX1*XX0 + XY1*YX0; + X_Y = XX1*XY0 + XY1*YY0; + X_O = XX1*XO0 + XY1*YO0 + XO1; + + Y_X = YX1*XX0 + YY1*YX0; + Y_Y = YX1*XY0 + YY1*YY0; + Y_O = YX1*XO0 + YY1*YO0 + YO1; + + } + dX = sqrt(d2X/N - dX*dX/(N*N)); /* scatter in pixels in the X direction */ + dY = sqrt(d2Y/N - dY*dY/(N*N)); /* scatter in pixels in the Y direction */ + dS = hypot (dX, dY) / sqrt(N); + if (VERBOSE) { + fprintf (stderr, "scatter in pixels: %5.2f x %5.2f -- %5.2f %d %d %d %d\n", dX, dY, dS, halt, last, extras, iteration); + } +# if (0) + if (!halt && !extras) { + if (iteration > 0) { + dSS = (Sprev - dS) / Sprev; + } else { + dSS = 1.0; + } + if (dSS < 0.05) { /* no fractional improvement, stop at last value */ + radius = Rprev * fratio; + halt = TRUE; + /* recalculate fit parameters on first, get scatter on second, exit on third */ + } + Rprev = radius; + Sprev = dS; + } +# endif + radius /= fratio; + } + radius *= fratio; + + /* convert X_X, etc to coords */ + { + double X0, Y0, S1, S2, p11, p21, p12, p22; + double delt, A, B, C, D, dRot; + + *Nmatch = N; + + delt = 1.0 / (X_X*Y_Y - X_Y*Y_X); + X = (coords[0].crpix1 - X_O); + Y = (coords[0].crpix2 - Y_O); + X0 = delt * (X*Y_Y - Y*X_Y); + Y0 = delt * (Y*X_X - X*Y_X); + XY_to_RD (&RA, &DEC, X0, Y0, coords); + + 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; + + A = S1*p11*Y_Y - S2*p12*Y_X; B = S2*p12*X_X - S1*p11*X_Y; + C = S1*p21*Y_Y - S2*p22*Y_X; D = S2*p22*X_X - S2*p21*X_Y; + + coords[0].cdelt1 = sqrt (A*A + C*C); + coords[0].cdelt2 = sqrt (B*B + D*D); + + coords[0].pc1_1 = A / coords[0].cdelt1; + coords[0].pc1_2 = B / coords[0].cdelt2; + coords[0].pc2_1 = C / coords[0].cdelt1; + coords[0].pc2_2 = D / coords[0].cdelt2; + + coords[0].cdelt1 = coords[0].cdelt1 * delt; + coords[0].cdelt2 = coords[0].cdelt2 * delt; + if ((fabs(coords[0].cdelt1) > 2.0*ASEC_PIX/3600.0) || + (fabs(coords[0].cdelt1) < 0.5*ASEC_PIX/3600.0) || + (fabs(coords[0].cdelt2) > 2.0*ASEC_PIX/3600.0) || + (fabs(coords[0].cdelt2) < 0.5*ASEC_PIX/3600.0)) { + fprintf (stderr, "ERROR: absurd solution\n"); + return (FALSE); + } + + dRot = coords[0].crval1 - RA ; + coords[0].crval1 = RA; + coords[0].crval2 = DEC; + + A = cos (dRot*RAD_DEG); + B = sin (dRot*RAD_DEG); + + p11 = coords[0].pc1_1; p12 = coords[0].pc1_2; + p21 = coords[0].pc2_1; p22 = coords[0].pc2_2; + + coords[0].pc1_1 = p11*A - p12*B; + coords[0].pc1_2 = p11*B + p12*A; + coords[0].pc2_1 = p21*A - p22*B; + coords[0].pc2_2 = p21*B + p22*A; + + } + + while (coords[0].crval1 < 0) coords[0].crval1 += 360.0; + while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0; + + *DR = sqrt (SQ(dX*coords[0].cdelt1*3600.0) + SQ(dY*coords[0].cdelt1*3600.0)); + *Radius = radius; + + if ((mode == 1) && VERBOSE) { + fprintf (stderr, "linear astrometric solution:\n"); + fprintf (stderr, "mode: %s\n", coords[0].ctype); + fprintf (stderr, "ref value: %f %f\n", coords[0].crval1, coords[0].crval2); + fprintf (stderr, "ref pixel: %f %f\n", coords[0].crpix1, coords[0].crpix2); + fprintf (stderr, "ra terms: %f %f\n", coords[0].pc1_1, coords[0].pc1_2); + fprintf (stderr, "dec terms: %f %f\n", coords[0].pc2_1, coords[0].pc2_2); + fprintf (stderr, "plt scale: %f %f\n", coords[0].cdelt1, coords[0].cdelt2); + fprintf (stderr, "accuracy: %f %f\n", *DR, *DR / sqrt ((double)(*Nmatch))); + } + + return (TRUE); + + + +} Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfitpoly.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfitpoly.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gfitpoly.c (revision 21724) @@ -0,0 +1,389 @@ +# include "gastro.h" + +void gfitpoly (SStars *stars1, SStars *stars2, int N1, int N2, Coords *coords, double *Radius, double *DR, int *Nmatch) { + + int i, j, m, n, M, N; + int first_j, last; + int *tmpN1, *tmpN2; + int NORDER, NTERM, NPARS, NPOWR; + double **sum, **xsum, **ysum; + double **matrix, **vector; + double Dx, Dy, DD, dX, dY, d2X, d2Y; + double radius, radius2; + double *tmpX1, *tmpX2, *tmpY1, *tmpY2; + double tX1, tX2, tY1, tY2; + double xterm, yterm, term, max; + + NORDER = NPOLYTERMS; + NPOWR = NORDER + 1; + NTERM = 2*NORDER + 1; + NPARS = (NORDER + 1)*(NORDER + 2) / 2; + if (NPOLYTERMS < 2) { + coords[0].Npolyterms = 0; + { + double nominal_det, measure_det, min_det, max_det, d1, d2, diffangle, tmp; + + nominal_det = CCD_PC1_1 * CCD_PC2_2 - CCD_PC1_2 * CCD_PC2_1; + min_det = nominal_det / 1.05; + max_det = nominal_det * 1.05; + if (min_det > max_det) { + tmp = max_det; max_det = min_det; min_det = tmp; + } + + measure_det = coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1; + if ((measure_det > max_det) || (measure_det < min_det)) { + fprintf (stderr, "absurd solution, not cartesian\n"); + *DR = 1e9; + } + d1 = hypot (coords[0].pc1_2, coords[0].pc1_1); + d2 = hypot (coords[0].pc2_2, coords[0].pc2_1); + diffangle = fabs (coords[0].pc2_1*coords[0].pc1_1 + coords[0].pc1_2*coords[0].pc2_2) / (d1*d2); + if (diffangle > MAX_NONLINEAR) { + fprintf (stderr, "absurd solution, not cartesian\n"); + *DR = 1e9; + } + } + return; + } + + fprintf (stderr, "\nattempting higher order fit\n"); + + /* allocate space for star coords */ + ALLOCATE (tmpX1, double, N1); + ALLOCATE (tmpY1, double, N1); + ALLOCATE (tmpN1, int, N1); + + ALLOCATE (tmpX2, double, N2); + ALLOCATE (tmpY2, double, N2); + ALLOCATE (tmpN2, int, N2); + + /* assign and sort list */ + for (i = 0; i < N1; i++) { + tmpX1[i] = stars1[i].X; + tmpY1[i] = stars1[i].Y; + tmpN1[i] = i; + } + if (N1 > 1) sort_lists (tmpX1, tmpY1, tmpN1, N1); + for (i = 0; i < N2; i++) { + tmpX2[i] = stars2[i].X; + tmpY2[i] = stars2[i].Y; + tmpN2[i] = i; + } + if (N2 > 1) sort_lists (tmpX2, tmpY2, tmpN2, N2); + + /* choose iteration ranges */ + radius = MINIMUM_RADIUS; + radius2 = radius*radius; + + /* 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)); + } + + /* do the following loop twice. + on first pass, find the coeffs. + on next pass, find just the residuals */ + + dX = dY = d2X = d2Y = N = 0; + for (last = 0; last < 2; last++) { + if (last) { + /* assign values based on determined coeffs */ + for (i = 0; i < N1; i++) { + tmpX1[i] = tmpY1[i] = 0; + yterm = 1; + for (m = 0; m < NPOWR; m++) { + xterm = 1; + for (n = 0; n < NPOWR - m; n++) { + tmpX1[i] += xterm*yterm*xsum[n][m]; + tmpY1[i] += xterm*yterm*ysum[n][m]; + xterm *= stars1[i].X; + } + yterm *= stars1[i].Y; + } + tmpN1[i] = i; + } + if (N1 > 1) sort_lists (tmpX1, tmpY1, tmpN1, N1); + dX = dY = d2X = d2Y = N = 0; + } + /* find matched stars */ + for (i = j = 0; (i < N1) && (j < N2); ) { + tX1 = tmpX1[i]; + tX2 = tmpX2[j]; + Dx = tX1 - tX2; + if (Dx <= -2.0*radius) { + i++; + continue; + } + if (Dx >= 2.0*radius) { + j++; + continue; + } + /* in the right range */ + first_j = j; + for (; (Dx > -2.0*radius) && (j < N2); j++) { + tY1 = tmpY1[i]; + tX2 = tmpX2[j]; + tY2 = tmpY2[j]; + Dx = tX1 - tX2; + Dy = tY1 - tY2; + DD = Dx*Dx + Dy*Dy; + /* stars matched */ + if (DD < radius2) { + if (last) { /* calculate residuals */ + dX += Dx; + dY += Dy; + d2X += Dx*Dx; + d2Y += Dy*Dy; + N += 1.0; + } else { /* accumulate data for coeffs */ + 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] += tX2*term; + ysum[n][m] += tY2*term; + } + yterm *= tY1; + } + xterm *= tX1; + } + } + } + } + j = first_j; + i++; + } + + if (!last) { /* calculate polyterm coeffs */ + fprintf (stderr, "matched %.0f stars for polyterms\n", sum[0][0]); + 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; + } + /* svd (matrix, NPARS, vector, 2); */ + dgaussj (matrix, NPARS, vector, 2); + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + xsum[n][m] = vector[i][0]; + ysum[n][m] = vector[i][1]; + } + } + i = 0; + for (m = 0; m < NPOWR; m++) { + for (n = 0; n < NPOWR - m; n++, i++) { + 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]); + } + } + } else { + fprintf (stderr, "%d stars matched for residuals\n", N); + + dX = sqrt(d2X/N - dX*dX/(N*N)); /* scatter in pixels in the X direction */ + dY = sqrt(d2Y/N - dY*dY/(N*N)); /* scatter in pixels in the Y direction */ + if (VERBOSE) fprintf (stderr, "scatter in pixels: %5.2f x %5.2f -- %5.2f\n", dX, dY, hypot(dX,dY) / sqrt(N)); + *DR = sqrt (SQ(dX*coords[0].cdelt1*3600.0) + SQ(dY*coords[0].cdelt1*3600.0)); + *Nmatch = N; + + } + } + + /* convert new terms to adjustments in coords and to polyterms */ + { + 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; + 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; + + { + double nominal_det, measure_det, min_det, max_det, tmp; + + nominal_det = CCD_PC1_1 * CCD_PC2_2 - CCD_PC1_2 * CCD_PC2_1; + min_det = nominal_det / 1.05; + max_det = nominal_det * 1.05; + if (min_det > max_det) { + tmp = max_det; max_det = min_det; min_det = tmp; + } + + measure_det = coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1; + if ((measure_det > max_det) || (measure_det < min_det)) { + fprintf (stderr, "absurd solution, not cartesian: %f (%f - %f)\n", measure_det, max_det, min_det); + *DR = 1e9; + } + } +} + +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); +} + + + +/********************** + + vector vs polyterms: + + the variable 'vector' is similar to, but not exactly like polyterms. + vector[i][j] provides coeffs for all of the x,y terms, + polyterms[i][j] only provides coeffs for the terms or order > 2. + + vector[i][j] and polyterms[i][j] also use a slightly different order: + + vector[i][j] runs in this order: + + n m norder = 3 + vector[0][0] * x^0 * y^0 + vector[1][0] * x^1 * y^0 + vector[2][0] * x^2 * y^0 + vector[3][0] * x^3 * y^0 + vector[4][0] * x^0 * y^1 + vector[5][0] * x^1 * y^1 + vector[6][0] * x^2 * y^1 + vector[7][0] * x^0 * y^2 + vector[8][0] * x^1 * y^2 + vector[9][0] * x^0 * y^3 + + to generate the vector entry from n, m, norder: + N = 0; + for (i = 0; i < m; i++) { + N += (norder - i + 1); + } + N += n; + + polyterms[i][j] runs in this order: + + n m + polyterms[0][0] * x^2 * y^0 = vector[2][0] + polyterms[1][0] * x^1 * y^1 = vector[5][0] + polyterms[2][0] * x^0 * y^2 = vector[7][0] + polyterms[3][0] * x^3 * y^0 = vector[3][0] + polyterms[4][0] * x^2 * y^1 = vector[6][0] + polyterms[5][0] * x^1 * y^2 = vector[8][0] + polyterms[6][0] * x^0 * y^3 = vector[9][0] + + to generate the polyterms entry from n, m, norder: + + N = 0; + nt = n + m; + for (i = 2; i < nt; i++) { + N += i + 1; + } + N += m; + + */ Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gheader.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gheader.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gheader.c (revision 21724) @@ -0,0 +1,136 @@ +# include "gastro.h" + +void gheader (char *file, Coords coords, double dR, int Nmatch) { + + Header header; + FILE *f, *g; + int i, oldsize, nbytes, status; + char line[1024]; + + if (!gfits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (3)\n", file); + exit(0); + } + oldsize = header.size; + + if (NEWPHOTCODE) { + /* we are going to write the photcode into the header + here we are just checking that the photcode provided + is a valid code */ + if (!LoadPhotcodes (PhotCodeFile)) { + fprintf (stderr, "error loading photcodes\n"); + exit (0); + } + if (!GetPhotcodeCodebyName (PHOTCODE)) { + fprintf (stderr, "ERROR: photcode not found in photcode table\n"); + exit (0); + } + gfits_modify (&header, "PHOTCODE", "%s", 1, PHOTCODE); + } + + if (Nmatch < 2) { + gfits_modify (&header, "NASTRO", "%d", 1, 0); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + goto skipstuff; + } + + gfits_modify (&header, "NASTRO", "%d", 1, Nmatch); + gfits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + + /*** use PutCoords to update header ***/ + if (coords.Npolyterms > 1) { + gfits_modify (&header, "CTYPE1", "%s", 1, "RA---PLY"); + gfits_modify (&header, "CTYPE2", "%s", 1, "DEC--PLY"); + } else { + gfits_modify (&header, "CTYPE1", "%s", 1, "RA---TAN"); + gfits_modify (&header, "CTYPE2", "%s", 1, "DEC--TAN"); + } + gfits_modify (&header, "CDELT1", "%le", 1, coords.cdelt1); + gfits_modify (&header, "CDELT2", "%le", 1, coords.cdelt2); + gfits_modify (&header, "CRVAL1", "%lf", 1, coords.crval1); + gfits_modify (&header, "CRVAL2", "%lf", 1, coords.crval2); + gfits_modify (&header, "CRPIX1", "%lf", 1, coords.crpix1); + gfits_modify (&header, "CRPIX2", "%lf", 1, coords.crpix2); + gfits_modify (&header, "PC001001", "%le", 1, coords.pc1_1); + gfits_modify (&header, "PC001002", "%le", 1, coords.pc1_2); + gfits_modify (&header, "PC002001", "%le", 1, coords.pc2_1); + gfits_modify (&header, "PC002002", "%le", 1, coords.pc2_2); + gfits_modify (&header, "NPLYTERM", "%d", 1, coords.Npolyterms); + if (coords.Npolyterms > 1) { + /* RA Terms */ + gfits_modify (&header, "PCA1X2Y0", "%le", 1, coords.polyterms[0][0]); /* polyterms[0]); */ + gfits_modify (&header, "PCA1X1Y1", "%le", 1, coords.polyterms[1][0]); /* polyterms[1]); */ + gfits_modify (&header, "PCA1X0Y2", "%le", 1, coords.polyterms[2][0]); /* polyterms[2]); */ + + if (coords.Npolyterms > 2) { + gfits_modify (&header, "PCA1X3Y0", "%le", 1, coords.polyterms[3][0]); /* polyterms[3]); */ + gfits_modify (&header, "PCA1X2Y1", "%le", 1, coords.polyterms[4][0]); /* polyterms[4]); */ + gfits_modify (&header, "PCA1X1Y2", "%le", 1, coords.polyterms[5][0]); /* polyterms[5]); */ + gfits_modify (&header, "PCA1X0Y3", "%le", 1, coords.polyterms[6][0]); /* polyterms[6]); */ + } + /* Dec Terms */ + gfits_modify (&header, "PCA2X2Y0", "%le", 1, coords.polyterms[0][1]); /* polyterms[7]); */ + gfits_modify (&header, "PCA2X1Y1", "%le", 1, coords.polyterms[1][1]); /* polyterms[8]); */ + gfits_modify (&header, "PCA2X0Y2", "%le", 1, coords.polyterms[2][1]); /* polyterms[9]); */ + + if (coords.Npolyterms > 2) { + gfits_modify (&header, "PCA2X3Y0", "%le", 1, coords.polyterms[3][1]); /* polyterms[10]); */ + gfits_modify (&header, "PCA2X2Y1", "%le", 1, coords.polyterms[4][1]); /* polyterms[11]); */ + gfits_modify (&header, "PCA2X1Y2", "%le", 1, coords.polyterms[5][1]); /* polyterms[12]); */ + gfits_modify (&header, "PCA2X0Y3", "%le", 1, coords.polyterms[6][1]); /* polyterms[13]); */ + } + } + gfits_modify (&header, "CERROR", "%lf", 1, dR); + gfits_modify (&header, "CERROR", "%C", 1, "scatter in astrometry soln (arcsec)"); + gfits_modify (&header, "CPRECISE", "%lf", 1, dR / sqrt(1.0*Nmatch)); + gfits_modify (&header, "CPRECISE", "%C", 1, "precision of astrometry soln (arcsec)"); + gfits_modify (&header, "EQUINOX", "%lf", 1, 2000.0); + /* we force equinox to be 2000.0 for all images */ + +skipstuff: + if (header.size > oldsize) { + if (VERBOSE) fprintf (stderr, "header expanded, creating new copy\n"); + sprintf (line, "mv %s %s~", file, file); + status = system (line); + if (status) { + fprintf (stderr, "ERROR: unable to create %s~, exiting\n", file); + exit (0); + } + sprintf (line, "%s~", file); + f = fopen (line, "r"); + g = fopen (file, "w"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find image file %s (4)\n", line); + exit(0); + } + if (g == NULL) { + fprintf (stderr, "ERROR: can't open output image file %s (4)\n", file); + exit(0); + } + nbytes = fwrite (header.buffer, 1, header.size, g); + fseek (f, oldsize, SEEK_SET); + for (i = 0; (nbytes = fread (header.buffer, 1, header.size, f)) > 0; i++) { + if (nbytes != fwrite (header.buffer, 1, nbytes, g)) { + fprintf (stderr, "ERROR: failure writing output data file\n"); + exit (0); + } + } + fclose (f); + fclose (g); + } else { + f = fopen (file, "r+"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find image file %s (4)\n", file); + exit(0); + } + + fseek (f, 0, SEEK_SET); + nbytes = fwrite (header.buffer, 1, header.size, f); + + fclose (f); + } + free (header.buffer); + +} + + Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gproject.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gproject.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gproject.c (revision 21724) @@ -0,0 +1,97 @@ +# include "gastro.h" +# define MMIN 2 +# define dM 0.5 +# define NMBIN 64 + +void gproject (SStars *catalog, SStars **stars, int Ncat, int *Nstars, Coords *coords, int NX, int NY, double dNdM, int N1) { + + double mbin[NMBIN], ratio; + int i, j, NSTARS, nstar, nstar2; + double X, Y, m0, Mmin, Mmax; + int XMIN, XMAX, YMIN, YMAX; + SStars *tstars, *t2stars; + + NSTARS = Ncat; + ALLOCATE (tstars, SStars, NSTARS); + + XMIN = 0.5*(1.0 - NFIELD)*NX; + XMAX = 0.5*(1.0 + NFIELD)*NX; + YMIN = 0.5*(1.0 - NFIELD)*NY; + YMAX = 0.5*(1.0 + NFIELD)*NY; + + /* project to local coords, select stars within region */ + for (nstar = i = 0; i < Ncat; i++) { + RD_to_XY (&X, &Y, catalog[i].X, catalog[i].Y, coords); + if ((X > XMIN) && (X < XMAX) && (Y > YMIN) && (Y < YMAX)) { + tstars[nstar].X = X; + tstars[nstar].Y = Y; + tstars[nstar].mag = catalog[i].mag; + nstar++; + if (CATDUMP) { + fprintf (stdout, "%f %f %f %f %f\n", catalog[i].X, catalog[i].Y, X, Y, catalog[i].mag); + } + } + } + if (CATDUMP) { + exit (0); + } + + REALLOCATE (tstars, SStars, nstar); + if (VERBOSE) fprintf (stderr, "%d total reference stars\n", nstar); + + /* find appropriate magnitude range */ + ratio = NFIELD * NFIELD; + m0 = 0; + bzero (mbin, NMBIN * sizeof (double)); + for (i = 0; i < nstar; i++) { + if (tstars[i].mag < MMIN) { + fprintf (stderr, "%d %f %f %f\n", i, tstars[i].X, tstars[i].Y, tstars[i].mag); + } + j = (tstars[i].mag - MMIN) / dM; + j = MIN (MAX (j, 0), (NMBIN - 1)); + mbin[j] ++; + } + for (i = 0; i < NMBIN; i++) { + if (ratio * dNdM < mbin[i] / dM) { + m0 = (i - 1) * dM + MMIN; + break; + } + } + + Mmin = m0 - 1.0; + Mmax = m0 + MAX (N1/dNdM, 1) + 1.0; + if (m0 == 0) { + Mmin = 0; + Mmax = 32; + } + if (VERBOSE) fprintf (stderr, "choosing magnitude range for reference stars: "); + if (VERBOSE) fprintf (stderr, " %5.2f to %5.2f mags\n", Mmin, Mmax); + + + ALLOCATE (t2stars, SStars, nstar); + if (MAGLIMS) { + /* make cut on mags */ + for (nstar2 = i = 0; i < nstar; i++) { + if ((tstars[i].mag > Mmin) && (tstars[i].mag < Mmax)) { + t2stars[nstar2] = tstars[i]; + nstar2 ++; + } + } + } else { + bcopy (tstars, t2stars, nstar*sizeof(SStars)); + nstar2 = nstar; + } + + + free (tstars); + *stars = t2stars; + *Nstars = nstar2; + if (VERBOSE) fprintf (stderr, "%d reference stars in mag range\n\n", nstar2); + +} + + +/* in this routine, + catalog[i].X,Y are RA,DEC in dec degrees + stars[i].X.Y are projected coords in approx pixels on image + */ Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/granges.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/granges.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/granges.c (revision 21724) @@ -0,0 +1,40 @@ +# include "gastro.h" + +void granges (SStars *stars1, SStars *stars2, int N1, int N2, int NPIX, double *gx, double *gy, double *gx0, double *gy0) { + + int i; + double maxX1, minX1, maxY1, minY1; + double maxX2, minX2, maxY2, minY2; + double Xzero, Xrange, Yzero, Yrange; + + maxX1 = minX1 = stars1[0].X; + maxY1 = minY1 = stars1[0].Y; + for (i = 0; i < N1; i++) { + maxX1 = MAX (maxX1, stars1[i].X); + minX1 = MIN (minX1, stars1[i].X); + maxY1 = MAX (maxY1, stars1[i].Y); + minY1 = MIN (minY1, stars1[i].Y); + } + + maxX2 = minX2 = stars2[0].X; + maxY2 = minY2 = stars2[0].Y; + for (i = 0; i < N2; i++) { + maxX2 = MAX (maxX2, (stars2[i].X)); + minX2 = MIN (minX2, (stars2[i].X)); + maxY2 = MAX (maxY2, (stars2[i].Y)); + minY2 = MIN (minY2, (stars2[i].Y)); + } + + Xzero = minX1 - maxX2; + Yzero = minY1 - maxY2; + Xrange = ((maxX1 - minX1) + (maxX2 - minX2)); + Yrange = ((maxY1 - minY1) + (maxY2 - minY2)); + + *gx = (NPIX - 1.0) / Xrange; + *gy = (NPIX - 1.0) / Yrange; + *gx0 = (1.0 - NPIX)*Xzero/Xrange; + *gy0 = (1.0 - NPIX)*Yzero/Yrange; + + fprintf (stderr, "gx, gy: %f %f\n", *gx, *gy); + +} Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/greference.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/greference.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/greference.c (revision 21724) @@ -0,0 +1,119 @@ +# include "gastro.h" + +# define USNO 1 +# define GSC 0 + +int greference (SStars **cat, int *Ncat, Coords *coords, int NX, int NY) { + + int i; + CatStats catstats; + + if (VERBOSE) fprintf (stderr, "\nloading astrometric reference data from %s\n", REFCAT); + + *Ncat = 0; + + define_region (&catstats, coords, NX, NY); + ALLOCATE (*cat, SStars, 1); + + /* get stars from the USNO catalog for the given region */ + if (!strcmp (REFCAT, "USNO") || !strcmp (REFCAT, "BOTH")) { + while (catstats.RA[0] > 360.0) catstats.RA[0] -= 360.0; + while (catstats.RA[0] < 0.0) catstats.RA[0] += 360.0; + while (catstats.RA[1] > 360.0) catstats.RA[1] -= 360.0; + while (catstats.RA[1] < 0.0) catstats.RA[1] += 360.0; + + /* if RA crosses 0,360 boundary, do 2 passes */ + if (catstats.RA[0] > catstats.RA[1]) { + int Nusno1, Nusno2; + USNOdata *usno1, *usno2; + USNOstats usnostats; + CatStats substats; + + substats = catstats; + substats.RA[0] = 0.0; + usno1 = getusno (&usnostats, &substats, &Nusno1); + + substats = catstats; + substats.RA[1] = 360.0; + usno2 = getusno (&usnostats, &substats, &Nusno2); + + REALLOCATE (*cat, SStars, MAX (Nusno1 + Nusno2, 1)); + for (i = 0; i < Nusno1; i++) { + cat[0][i].X = usno1[i].R; + cat[0][i].Y = usno1[i].D; + cat[0][i].mag = fabs(usno1[i].r); + } + for (i = Nusno1; i < Nusno2; i++) { + cat[0][i].X = usno2[i].R; + cat[0][i].Y = usno2[i].D; + cat[0][i].mag = fabs(usno2[i].r); + } + *Ncat = Nusno1 + Nusno2; + free (usno1); + free (usno2); + } else { + int Nusno; + USNOdata *usno; + USNOstats usnostats; + + usno = getusno (&usnostats, &catstats, &Nusno); + + REALLOCATE (*cat, SStars, MAX (Nusno, 1)); + for (i = 0; i < Nusno; i++) { + cat[0][i].X = usno[i].R; + cat[0][i].Y = usno[i].D; + cat[0][i].mag = fabs(usno[i].r); + } + *Ncat = Nusno; + free (usno); + } + if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", *Ncat); + } + + if (!strcmp (REFCAT, "GSC") || !strcmp (REFCAT, "BOTH")) { + int j, Ngsc; + SStars *gsc; + + gsc = getgsc (&catstats, &Ngsc); + REALLOCATE (*cat, SStars, MAX (Ngsc + *Ncat, 1)); + for (j = *Ncat, i = 0; i < Ngsc; i++) { + cat[0][j] = gsc[i]; + } + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Ngsc); + *Ncat += Ngsc; + } + + if (!strcmp (REFCAT, "PTOLEMY")) { + free (*cat); + *cat = getptolemy (&catstats, Ncat); + if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", *Ncat); + } + return (TRUE); +} + +void define_region (CatStats *catstats, Coords *coords, int NX, int NY) { + + int i, j; + double X, Y, R, D; + + catstats[0].RA[0] = catstats[0].DEC[0] = 360.0; + catstats[0].RA[1] = catstats[0].DEC[1] = -360.0; + + for (i = -1; i < 2; i++) { + for (j = -1; j < 2; j++) { + X = 0.5*(1.0 + i*NFIELD)*NX; + Y = 0.5*(1.0 + j*NFIELD)*NY; + XY_to_RD (&R, &D, X, Y, coords); + /* coords returns a region all in same phase + while (R < 0.0) R += 360.0; + while (R >= 360.0) R -= 360.0; + */ + 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); + } + } + 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/ohana/gastro-1-4/Ohana/src/gastro/src/gstars.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gstars.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/gstars.c (revision 21724) @@ -0,0 +1,254 @@ +# include "gastro.h" +# define dcos(a) (cos((double)((a)*(RAD_DEG)))) +# define dsin(a) (sin((double)((a)*(RAD_DEG)))) + +/* by necesity hard wired */ +# define D_NSTARS 1000 +# define BYTES_STAR 66 +# define BLOCK 1000 +# include +# include + +SStars *gstars (char *file, int *NSTARS, Coords *coords, int *NX, int *NY, double *dNdM) { + + char line[64], side[64]; + Header header, theader; + FILE *f; + int j, Ninstar, nstar, rnumber, N, Nstars, nbytes, Nbytes; + SStars *stars; + char *buffer; + double X, Y, T1, T2, T3, type, csign; + double PD, PR, DE, RE; + double ra, dec, dmag; + + /* read in image header, open image data region */ + if (!gfits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (1)\n", file); + exit(0); + } + /* get complete info from header */ + gfits_scan (&header, "NAXIS1", "%d", 1, NX); + gfits_scan (&header, "NAXIS2", "%d", 1, NY); + + /* attempt to get detailed astrometric information from header */ + if (!strcasecmp (ROUGH_ASTROMETRY, "header")) { + if (!HEADER[0]) { + GetCoords (coords, &header); + } else { + gfits_read_header (HEADER, &theader); + GetCoords (coords, &theader); + } + if (!strcmp (coords[0].ctype, "NONE") || (coords[0].cdelt1 == 0) || (coords[0].cdelt2 == 0)) { + fprintf (stderr, "header coordinates incomplete, trying for rough coordinates\n"); + strcpy (ROUGH_ASTROMETRY, "config"); + } else { + if (FLIPX) { + coords[0].pc1_1 *= -1; + coords[0].crpix1 = *NX - coords[0].crpix1; + } + if (FLIPY) { + coords[0].pc2_2 *= -1; + coords[0].crpix2 = *NY - coords[0].crpix2; + } + ASEC_PIX = fabs (coords[0].cdelt1 * 3600.0); + csign = coords[0].cdelt1 / fabs (coords[0].cdelt1); + CCD_PC1_1 = coords[0].pc1_1 * csign; + CCD_PC2_1 = coords[0].pc2_1 * csign; + csign = coords[0].cdelt2 / fabs (coords[0].cdelt2); + CCD_PC1_2 = coords[0].pc1_2 * csign; + CCD_PC2_2 = coords[0].pc2_2 * csign; + if (!strcmp (&coords[0].ctype[4], "-PLY")) { + strcpy (coords[0].ctype, "DEC--TAN"); + } + } + } + + /* get just RA & DEC from header, other terms from config file */ + if (!strcasecmp (ROUGH_ASTROMETRY, "config")) { + /* default values for coords */ + strcpy (coords[0].ctype, "RA---TAN"); + coords[0].pc1_1 = CCD_PC1_1; coords[0].pc1_2 = CCD_PC1_2; + coords[0].pc2_1 = CCD_PC2_1; coords[0].pc2_2 = CCD_PC2_2; + coords[0].cdelt1 = coords[0].cdelt2 = ASEC_PIX / 3600.0; + coords[0].Npolyterms = 0; + coords[0].crpix1 = 0.5*(*NX); + coords[0].crpix2 = 0.5*(*NY); + + /* get RA & DEC from header, unless FORCE is ste */ + if (!FORCE) { + /* RA (in hours, not degrees) */ + if (!gfits_scan (&header, "RA", "%s", 1, line)) { + fprintf (stderr, "ERROR: no astrometry in header\n"); + exit (1); + } + dms_to_ddd (&coords[0].crval1, line); + coords[0].crval1 = coords[0].crval1 * 15.0; + + /* DEC */ + if (!gfits_scan (&header, "DEC", "%s", 1, line)) { + fprintf (stderr, "ERROR: no astrometry in header\n"); + exit (1); + } + dms_to_ddd (&coords[0].crval2, line); + } + } + if (VERBOSE) fprintf (stderr, "coordinates from header: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2); + + /* use RA & DEC from command line arguments */ + if (FORCE) { + coords[0].crval1 = F_RA; + coords[0].crval2 = F_DEC; + if (VERBOSE) fprintf (stderr, " forcing coordinates to: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2); + } + + /* the following two sections are LONEOS derived and may not be needed elsewhere */ + if (LONEOS_COORDS) { + gfits_scan (&header, "COMMENT", "%s", 1, line); + sscanf (line, "%*s%d%s", &rnumber, side); + if (get_region_coords (&ra, &dec, rnumber, side)) { + if (fabs(ra - coords[0].crval1) > 0.1) { + fprintf (stderr, "large offset from claimed position, using region coords %f %f -> %f %f (%d %s)\n", + coords[0].crval1, coords[0].crval2, ra, dec, rnumber, side); + coords[0].crval1 = ra; + coords[0].crval2 = dec; + } + } + } + + /* 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 */ + if (POLAR_ALIGNMENT) { + X = coords[0].crval1; + Y = coords[0].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); + + coords[0].crval1 = (DEG_RAD * atan2 (T2, T1)) + PR; + coords[0].crval2 = (DEG_RAD * asin (T3)); + while (coords[0].crval1 < 0) coords[0].crval1 += 360.0; + while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0; + + if (VERBOSE) fprintf (stderr, " after polar alignment: %9.4f %9.4f\n", coords[0].crval1, coords[0].crval2); + } + + Nstars = 0; + gfits_scan (&header, "NSTARS", "%d", 1, &Nstars); + if (Nstars == -1) { + fprintf (stderr, "ERROR: failed to find NSTARS\n"); + exit (0); + } + ALLOCATE (stars, SStars, Nstars); + Nbytes = Nstars*BYTES_STAR; + + /* re-open file for stars */ + f = fopen (file, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find image file %s (2)\n", file); + exit(0); + } + fseek (f, header.size, SEEK_SET); + + N = nstar = 0; + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR)); + + while ((nbytes = fread (buffer, 1, (BLOCK*BYTES_STAR), f)) != 0) { + Ninstar = nbytes / BYTES_STAR; + for (j = 0; j < Ninstar; j++, nstar++) { + dparse (&stars[N].X, 1, &buffer[j*BYTES_STAR]); + dparse (&stars[N].Y, 2, &buffer[j*BYTES_STAR]); + dparse (&stars[N].mag, 3, &buffer[j*BYTES_STAR]); + dparse (&dmag, 4, &buffer[j*BYTES_STAR]); + dparse (&type, 5, &buffer[j*BYTES_STAR]); + if ((type == 4) || (type == 6) || (type == 5) || (type == 9)) continue; + if (dmag > 100) continue; + N++; + } + } + free (header.buffer); + free (buffer); + fclose (f); + + if (nstar != Nstars) { + fprintf (stderr, "WARNING: only read %d of %d stars\n", nstar, Nstars); + } + if (N < 5) { + fprintf (stderr, "ERROR: too few stars for reliable solution, only %d\n", + N); + exit (0); + } + + sort_stars (stars, N); /* sorting by magnitude */ + Nstars = N; + + if (VERBOSE) fprintf (stderr, "\nread %d stars from data file", Nstars); + if (*NSTARS < Nstars) { + REALLOCATE (stars, SStars, *NSTARS); + if (VERBOSE) fprintf (stderr, ", using %d\n", *NSTARS); + } else { + *NSTARS = Nstars; + if (VERBOSE) fprintf (stderr, "\n"); + } + + *dNdM = *NSTARS / (stars[(*NSTARS-1)].mag - stars[0].mag) ; + if (VERBOSE) fprintf (stderr, "brightest star in datafile: %f mag\n", stars[0].mag); + + return (stars); +} + +void sort_stars (SStars *stars, int N) { + + int l,j,ir,i; + SStars 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].mag < stars[j+1].mag) ++j; + if (tempstar.mag < stars[j].mag) { + stars[i] = stars[j]; + j += (i=j) + 1; + } + else j = ir + 1; + } + stars[i] = tempstar; + } +} + + +# if (0) + + /*** this is tailored for LONEOS ***/ + Nccd = -1; + gfits_scan (&header, "NCCD", "%d", 1, &Nccd); + if (Nccd == -1) { /* no ccd info in header, not loneos */ + coords[0].crpix1 = 0.5*(*NX); + coords[0].crpix2 = 0.5*(*NY); + } + if (Nccd == 0) { /* chip 0 (a) *** might be wrong *** */ + coords[0].crpix1 = 0; + coords[0].crpix2 = 0.5*(*NY); + } + if (Nccd == 1) { /* chip 1 (b) */ + coords[0].crpix1 = (*NX); + coords[0].crpix2 = 0.5*(*NY); + } + +# endif Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/misc.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/misc.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/misc.c (revision 21724) @@ -0,0 +1,44 @@ +# include "gastro.h" + +# define SIGN(X) (((X) == 0) ? 0 : ((fabs((double)(X))) / (X))) + +void hh_hms (double hh, int *hr, int *mn, double *sc) { + + int flag; + + flag = SIGN(hh); + hh *= flag; + hh = 24.0*(hh/24.0 - (int)(hh/24.0)); + *sc = 60.0*(60.0*hh - (int)(60.0*hh)); + *mn = 60.0*(hh - (int)hh); + *hr = (int) hh; + *hr *= flag; + +} + +void hms_format (char *line, double value) { + + int hr, mn; + double sc; + + hh_hms (value, &hr, &mn, &sc); + hr = (int) value; + if (isnan (value)) + sprintf (line, "xx:xx:xx.xx"); + else { + if (value < 0) { + sprintf (line, "-%02d:%02d:%05.2f", abs(hr), mn, sc); + } else { + sprintf (line, "%02d:%02d:%05.2f", hr, mn, sc); + } + } +} + +void area_of_region (CatStats *region) { + + double area; + + area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG)); + region[0].Area = area; +} + Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/plotstuff.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/plotstuff.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/plotstuff.c (revision 21724) @@ -0,0 +1,166 @@ +# include "gastro.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 %f %f ", + Npts, graphmode[0].style, + graphmode[0].ptype, graphmode[0].ltype, + graphmode[0].etype, 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/ohana/gastro-1-4/Ohana/src/gastro/src/rotate.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/rotate.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/rotate.c (revision 21724) @@ -0,0 +1,44 @@ +# include "gastro.h" + +void rotate (SStars *stars, int Nstars, double angle, int Xo, int Yo) { + + int i; + double dX, dY, DX, DY, CS, SN; + double theta, theta2; + + if (angle == 0.0) { + return; + } + + theta = (angle*RAD_DEG); + + if (fabs (angle) < 10) { + theta2 = 0.5*theta*theta; + for (i = 0; i < Nstars; i++) { + dX = (stars[i].X - Xo); + dY = (stars[i].Y - Yo); + stars[i].X += -theta*dY - theta2*dX; + stars[i].Y += theta*dX - theta2*dY; + } + } else { + + CS = cos (theta); + SN = sin (theta); + + for (i = 0; i < Nstars; i++) { + dX = (stars[i].X - Xo); + dY = (stars[i].Y - Yo); + + DX = dX * CS - dY * SN; + DY = dX * SN + dY * CS; + + stars[i].X = DX + Xo; + stars[i].Y = DY + Yo; + } + } + +} + + + /* rotate the list 'stars' by an angle ccw from x axis */ + Index: /tags/ohana/gastro-1-4/Ohana/src/gastro/src/sort.c =================================================================== --- /tags/ohana/gastro-1-4/Ohana/src/gastro/src/sort.c (revision 21724) +++ /tags/ohana/gastro-1-4/Ohana/src/gastro/src/sort.c (revision 21724) @@ -0,0 +1,46 @@ + +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; + } +}