Index: /branches/elixir/Ohana/src/gastro/Makefile =================================================================== --- /branches/elixir/Ohana/src/gastro/Makefile (revision 2434) +++ /branches/elixir/Ohana/src/gastro/Makefile (revision 2434) @@ -0,0 +1,158 @@ +include ../../Configure +HOME = $(ROOT)/src/gastro +CONFIG = $(ROOT)/config + +default: 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) -I$(LINC)/$(ARCH) +LIBS = -L$(LLIB) -lFITS -lsocket -lnsl -lohana -lm +CFLAGS = -o $*.$(ARCH).o $(INCS) +CCFLAGS = $(INCS) $(LIBS) + +OBJS = \ +$(SRC)/get_stars.$(ARCH).o $(SRC)/delta.$(ARCH).o \ +$(SRC)/line_fit.$(ARCH).o $(SRC)/ranges.$(ARCH).o \ +$(SRC)/find_shift.$(ARCH).o + +# +# ----------------------- +EOBJ = \ +$(SRC)/extract.$(ARCH).o $(SRC)/extract_catalog.$(ARCH).o \ +$(SRC)/precess.$(ARCH).o $(SRC)/stats.$(ARCH).o + + +$(EOBJ): $(INC)/astro.h + +extract: $(BIN)/extract + +$(BIN)/extract: $(EOBJ) + $(CC) $(EOBJ) -o $(BIN)/extract $(CCFLAGS) + +install.extract: + rm -f $(DESTBIN)/extract + cp $(BIN)/extract $(DESTBIN)/ +# ----------------------- +SOBJ = \ +$(SRC)/simulate.$(ARCH).o $(SRC)/gregions2.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o $(SRC)/coordops.$(ARCH).o \ +$(SRC)/gproject.$(ARCH).o $(SRC)/gaussian.$(ARCH).o + +$(SOBJ): $(INC)/astro.h + +simulate: $(BIN)/simulate.$(ARCH) + +$(BIN)/simulate: $(SOBJ) + $(CC) $(SOBJ) -o $(BIN)/simulate.$(ARCH) $(CCFLAGS) + +install.simulate: + rm -f $(DESTBIN)/simulate + cp $(BIN)/simulate.$(ARCH) $(DESTBIN)/simulate + +# ----------------------- +GOBJ = \ +$(SRC)/gastro.$(ARCH).o $(SRC)/gstars.$(ARCH).o \ +$(SRC)/gargs.$(ARCH).o $(SRC)/gregions2.$(ARCH).o \ +$(SRC)/gcatalog.$(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)/gaussj.$(ARCH).o \ +$(SRC)/get_region_coords.$(ARCH).o \ +$(SRC)/greference.$(ARCH).o $(SRC)/getusno.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o $(SRC)/gptolemy.$(ARCH).o \ +$(SRC)/g2mass.$(ARCH).o $(SRC)/misc.$(ARCH).o + +$(GOBJ): $(INC)/gastro.h + +gastro: $(BIN)/gastro.$(ARCH) + +$(BIN)/gastro.$(ARCH): $(GOBJ) + $(CC) $(GOBJ) -o $(BIN)/gastro.$(ARCH) $(CCFLAGS) + +install.gastro: $(DESTBIN)/gastro + +$(DESTBIN)/gastro: $(BIN)/gastro.$(ARCH) + rm -f $(DESTBIN)/gastro + cp $(BIN)/gastro.$(ARCH) $(DESTBIN)/gastro +# ----------------------- +DOBJ = \ +$(SRC)/gscdump.$(ARCH).o $(SRC)/gstars.$(ARCH).o \ +$(SRC)/gargs.$(ARCH).o $(SRC)/gregions2.$(ARCH).o \ +$(SRC)/gcatalog.$(ARCH).o $(SRC)/coordops.$(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)/config.$(ARCH).o $(SRC)/ConfigInit.$(ARCH).o \ +$(SRC)/sort.$(ARCH).o + +$(DOBJ): $(INC)/gastro.h + +gscdump: $(BIN)/gscdump.$(ARCH) + +$(BIN)/gscdump.$(ARCH): $(DOBJ) + $(CC) $(DOBJ) -o $(BIN)/gscdump.$(ARCH) $(CCFLAGS) + +install.gscdump: $(DESTBIN)/gscdump + +$(DESTBIN)/gscdump: $(BIN)/gscdump.$(ARCH) + rm -f $(DESTBIN)/gscdump + cp $(BIN)/gscdump.$(ARCH) $(DESTBIN)/gscdump +# ----------------------- +AOBJ = $(OBJS) \ +$(SRC)/stats.$(ARCH).o $(SRC)/astro.$(ARCH).o \ +$(SRC)/get_catalog_stars.$(ARCH).o $(SRC)/extract_catalog.$(ARCH).o \ +$(SRC)/precess.$(ARCH).o $(SRC)/alter_header.$(ARCH).o \ +$(SRC)/args.$(ARCH).o + +$(AOBJ): $(INC)/astro.h + +astro: $(BIN)/astro + +$(BIN)/astro: $(AOBJ) + $(CC) $(AOBJ) -o $(BIN)/astro $(CCFLAGS) + +# ----------------------- +ROBJ = $(OBJS) \ +$(SRC)/rastro.$(ARCH).o $(SRC)/ralter_header.$(ARCH).o \ +$(SRC)/rargs.$(ARCH).o $(SRC)/get_argument.$(ARCH).o \ +$(SRC)/remove_argument.$(ARCH).o + +$(ROBJ): $(INC)/astro.h + +rastro: $(BIN)/rastro + +$(BIN)/rastro: $(ROBJ) + $(CC) $(ROBJ) -o $(BIN)/rastro $(CCFLAGS) + +# ----------------------- +install: install.gastro + +clean: + rm -f $(BIN)/*.$(ARCH) + rm -f `find . -name "*.o"` + rm -f `find . -name "*~"` + rm -f `find . -name "#*"` + +dist: clean + rm -f $(BIN)/astro + rm -f $(BIN)/rastro + + +.SUFFIXES: .$(ARCH).o + +.c.$(ARCH).o: + $(CC) $(CFLAGS) -c $< + + Index: /branches/elixir/Ohana/src/gastro/include/astro.h =================================================================== --- /branches/elixir/Ohana/src/gastro/include/astro.h (revision 2434) +++ /branches/elixir/Ohana/src/gastro/include/astro.h (revision 2434) @@ -0,0 +1,63 @@ +# include + +# define DEFAULT_RADIUS 8.0 /* rough search radius, pixels */ +# define DEFAULT_NPIX 100 +# define CCD_NX 2048 /* CCD width -- not used, determined from header */ +# define CCD_NY 4096 /* CCD height -- not used, determined from header */ +# define CCD_PC1_1 1 /* CCD orientation */ +# define CCD_PC2_2 -1 +# define CCD_PC1_2 0 +# define CCD_PC2_1 0 +# define CCD_DELT_X 1 +# define CCD_DELT_Y -1 +# define NFIELD 1.5 /* search region */ +# define MMIN 9.0 /* very saturated */ + +typedef struct { + double X; + double Y; + double mag; +} Stars; + +typedef struct { + char ctype[16]; + double crval1, crpix1, cdelt1; + double crval2, crpix2, cdelt2; + double pc1_1, pc1_2; + double pc2_1, pc2_2; +} Coords; + +static char GSCFILE[] = "./GSCregions.tbl"; +static char GSC_DIR[] = "./gsc"; + +void args PROTO((int *, char **)); +void rargs PROTO((int *, char **)); +void alter_header PROTO((char *, char **, int, double, double, double, double, double, double, double, double, int)); +void ralter_header PROTO((char *, char *, double, double, double, double, double, double, double, double, int)); +void delta PROTO((Stars *, Stars *, int, int, double *, double *, double *, double *)); +void find_shift PROTO((Stars *, Stars *, int, int, double, double, int *, double *, double *, double *)); +Stars *get_stars PROTO((char *, int *)); +int line_fit PROTO((Stars *, Stars *, int, int, double, double, double, double, + double *, double *, double *, double *, double *, double *, double *, double *)); +void ranges PROTO((Stars *, Stars *, int, int, double *, double *, double *, double *, double, double)); +void stats PROTO((char *, double *, double *, double *, double *, double)); +double dms_to_deg PROTO((char *, int)); +void extract_catalog PROTO((char *, double **, double **, double **, int *)); +Stars *get_catalog_stars PROTO((double *, double *, double *, int, double, double, double, double, int *)); +void precess PROTO((double *, double *, double, double)); + +extern double hypot PROTO((double, double)); +/* this seems to be a problem: is not included from math.h with the -ansi flag */ + +double RADIUS; +double ASEC_PIX; +int FLIP; +int ROTATE; +int N2NUMBER; +int NPIX; +int FORCE; +int DUMP; +double F_RA; +double F_DEC; +double NIM; + Index: /branches/elixir/Ohana/src/gastro/include/gastro.h =================================================================== --- /branches/elixir/Ohana/src/gastro/include/gastro.h (revision 2434) +++ /branches/elixir/Ohana/src/gastro/include/gastro.h (revision 2434) @@ -0,0 +1,140 @@ +# 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 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], GSC_DIR[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; + +typedef struct { + char name[10]; + double RA[2], DEC[2]; +} CatalogRegion; + +/* this seems to be a problem: is not included from math.h with the -ansi flag */ +extern double hypot PROTO((double, double)); + +char **find_dec_bands PROTO((CatalogRegion area, char **names, int *Nnames)); +char **gregions PROTO((CatStats *catstats, int *nnames)); +char **load_ra_blocks PROTO((int Ndec, CatalogRegion area, char **names, int *Ntnames)); +SStars *get_catalog_stars PROTO((double *, double *, double *, int, double, double, double, double, int *)); +SStars *get_stars PROTO((char *, int *)); +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 ahelp 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 args PROTO((int *, char **)); +void define_region PROTO((CatStats *catstats, Coords *coords, int NX, int NY)); +void delta PROTO((SStars *, SStars *, int, int, double *, double *, double *, double *)); +int dms_to_ddd PROTO((double *Value, char *string)); +void extract_catalog PROTO((char *, double **, double **, double **, int *)); +void find_shift PROTO((SStars *, SStars *, int, int, double, double, int *, double *, double *, double *)); +int g2mass PROTO((char *fullpath, SStars **stars, int *Nstars)); +void gargs PROTO((int *argc, char **argv, Coords *coords)); +int gaussj PROTO((double **a, int n, double **b, int m)); +int gcatalog PROTO((char *filename, SStars **stars, int *Nstars)); +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)); +int gptolemy PROTO((char *fullpath, SStars **stars, int *Nstars)); +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)); +int parse_GSC_line PROTO((CatalogRegion *region, char *line)); +void precess PROTO((double *, double *, double, double)); +void ralter_header PROTO((char *, char *, double, double, double, double, double, double, double, double, int)); +void ranges PROTO((SStars *, SStars *, int, int, double *, double *, double *, double *, double, double)); +void rargs PROTO((int *, char **)); +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)); +int str_to_radec PROTO((double *ra, double *dec, char *str1, char *str2)); +void svd PROTO((double ** matrix, int N, double **vector, int M)); +void svdcmp PROTO((double **a, int m, int n, double **w, double **v)); Index: /branches/elixir/Ohana/src/gastro/src/ConfigInit.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/ConfigInit.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/ConfigInit.c (revision 2434) @@ -0,0 +1,56 @@ +# 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); + 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_ERROR", "%lf", 0, &MAX_ERROR); + 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, "MMIN", "%lf", 0, &MMIN); + 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, "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, "ASTRO_REFCAT", "%s", 0, REFCAT); + ScanConfig (config, "CATDIR", "%s", 0, CATDIR); + ScanConfig (config, "ROUGH_ASTROMETRY", "%s", 0, ROUGH_ASTROMETRY); + ScanConfig (config, "PHOTCODE_FILE", "%s", 0, PhotCodeFile); + + 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: /branches/elixir/Ohana/src/gastro/src/alter_header.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/alter_header.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/alter_header.c (revision 2434) @@ -0,0 +1,60 @@ +# include "astro.h" + +void alter_header (head, argv, argc, X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY, N) +char head[], **argv; +int argc; +double X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY; +int N; +{ + + Header header; + char tmp[200], line[500], *c; + int i; + + fits_read_header (head, &header); + + fits_modify (&header, "RA_O", "%lf", 1, X_O); + fits_modify (&header, "DEC_O", "%lf", 1, Y_O); + if (ROTATE) { + fits_modify (&header, "RA_X", "%le", 1, X_Y); + fits_modify (&header, "RA_Y", "%le", 1, X_X); + fits_modify (&header, "DEC_X", "%le", 1, Y_Y); + fits_modify (&header, "DEC_Y", "%le", 1, Y_X); + } + else { + fits_modify (&header, "RA_X", "%le", 1, X_X); + fits_modify (&header, "RA_Y", "%le", 1, X_Y); + fits_modify (&header, "DEC_X", "%le", 1, Y_X); + fits_modify (&header, "DEC_Y", "%le", 1, Y_Y); + } + fits_modify (&header, "dRA", "%lf", 1, dX); + fits_modify (&header, "dDEC", "%lf", 1, dY); + fits_modify (&header, "Nastro", "%d", 1, N); + for (i = 3; i < argc; i++) { + if ((c = strrchr (argv[i], '/')) != NULL) + fits_modify (&header, "REF", "%s", i-2, c+1); + else + fits_modify (&header, "REF", "%s", i-2, argv[i]); + } + + fits_modify (&header, "RA_O", "%C", 1, "astro coeff"); + fits_modify (&header, "RA_X", "%C", 1, "astro coeff"); + fits_modify (&header, "RA_Y", "%C", 1, "astro coeff"); + fits_modify (&header, "dRA", "%C", 1, "error (pixels)"); + fits_modify (&header, "DEC_O", "%C", 1, "astro coeff"); + fits_modify (&header, "DEC_X", "%C", 1, "astro coeff"); + fits_modify (&header, "DEC_Y", "%C", 1, "astro coeff"); + fits_modify (&header, "dDEC", "%C", 1, "error (pixels)"); + fits_modify (&header, "Nastro", "%C", 1, "Number of stars used"); + for (i = 3; i < argc; i++) + fits_modify (&header, "REF", "%C", i-2, "astro ref file"); + + strcpy (tmp, head); + strcat (tmp, "~"); + sprintf (line, "mv %s %s\0", head, tmp); + system (line); + + fits_write_header (head, &header); + fits_free_header (&header); + +} Index: /branches/elixir/Ohana/src/gastro/src/args.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/args.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/args.c (revision 2434) @@ -0,0 +1,148 @@ +# include "astro.h" +# define NARGS 3 /* minimum is: astro ASEC_PIX catalog */ + +void ahelp () +{ + + fprintf (stderr, "astro -- an astrometry routine\n"); + + fprintf (stderr, " USAGE: astro , [,...] \n"); + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -N -- number of input stars used (default: 100) \n"); + fprintf (stderr, " -f RA,DEC -- force image center to RA,DEC \n"); + fprintf (stderr, " -R radius -- default radius of search (default: %f asec) \n", DEFAULT_RADIUS); + fprintf (stderr, " -G N -- grid search spacing (default: %d ) \n", DEFAULT_NPIX); + fprintf (stderr, " -F n -- quadrant for search (default: all) \n"); + fprintf (stderr, " -rot -- rotated axes (default: no) \n"); + fprintf (stderr, " (see below)\n"); + fprintf (stderr, "\n"); + fprintf (stderr, " -dump -- dumps the catalog stars to stdout\n"); + fprintf (stderr, " -help -- print this list\n"); + fprintf (stderr, "\n"); + fprintf (stderr, "\n"); + fprintf (stderr, " Effect of -F and -rot:\n"); + fprintf (stderr, " \n"); + fprintf (stderr, " without -rot: with -rot: \n"); + fprintf (stderr, " ------------- ---------- \n"); + fprintf (stderr, " 2 n n 1 2 e e 1 \n"); + fprintf (stderr, " | | | | \n"); + fprintf (stderr, " e - - e n - - n \n"); + fprintf (stderr, " \n"); + fprintf (stderr, " e - - e n - - n \n"); + fprintf (stderr, " | | | | \n"); + fprintf (stderr, " 3 n n 4 3 e e 4 \n"); + fprintf (stderr, "\n"); + exit (0); + +} + +void args (argc, argv) +int *argc; +char **argv; +{ + + int N; + char line[500]; + + if (N = get_argument (*argc, argv, "-dump")) { + remove_argument (N, argc, argv); + DUMP = TRUE; + } + else + DUMP = FALSE; + + if (get_argument (*argc, argv, "-help") || + get_argument (*argc, argv, "-h")) { + ahelp (); + } + + /** optional arguments **/ + if (N = get_argument (*argc, argv, "-N")) { + remove_argument (N, argc, argv); + N2NUMBER = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + N2NUMBER = 100; + } + + if (N = get_argument (*argc, argv, "-NIM")) { + remove_argument (N, argc, argv); + NIM = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + NIM = 2; + } + + if (N = get_argument (*argc, argv, "-F")) { + remove_argument (N, argc, argv); + FLIP = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + FLIP = 0; + } + + if (N = get_argument (*argc, argv, "-rot")) { + remove_argument (N, argc, argv); + ROTATE = TRUE; + } + else { + ROTATE = FALSE; + } + + if (N = get_argument (*argc, argv, "-f")) { + 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); + } + else { + FORCE = FALSE; + } + + if (N = get_argument (*argc, argv, "-R")) { + remove_argument (N, argc, argv); + RADIUS = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + RADIUS = DEFAULT_RADIUS; + } + + if (N = get_argument (*argc, argv, "-G")) { + remove_argument (N, argc, argv); + NPIX = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + NPIX = DEFAULT_NPIX; + } + + if (*argc < NARGS) { + fprintf (stderr, "USAGE: astro , [,...] \n"); + fprintf (stderr, " enter filenames, followed by EOF\n"); + exit (0); + } + + ASEC_PIX = atof (argv[1]); + remove_argument (1, argc, argv); + if (ASEC_PIX < 0.01) { + fprintf (stderr, "%f arcsec / pixel, is this really true?\n", + ASEC_PIX); + fscanf (stdin, "%s", line); + if (strcasecmp(line, "y") || strcasecmp(line, "yes")) + exit (0); + } + if (ASEC_PIX > 5) { + fprintf (stderr, "WARNING: your image may have too large\n"); + fprintf (stderr, " a field to be well modeled with\n"); + fprintf (stderr, " a linear astrometric solution\n"); + } + + +} + Index: /branches/elixir/Ohana/src/gastro/src/astro.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/astro.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/astro.c (revision 2434) @@ -0,0 +1,85 @@ +# include "astro.h" + +main (argc, argv) +int argc; +char **argv; +{ + + int i, N1, N2, Ncat, Nlast, Nused; + double dx, dy, Rx, Ry; + double RA, DEC, dRA, dDEC; + double *ra, *dec, *mag; + double X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY; + char *c, *ref, filename[200], head[200], obj[200], line[200]; + Stars *stars1, *stars2; + + args (&argc, argv); + + Nlast = Ncat = 0; + for (i = 1; i < argc; i++) { + extract_catalog (argv[i], &ra, &dec, &mag, &Ncat); + fprintf (stderr, "catalog %s, %d stars (%d total)\n", argv[i], Ncat-Nlast, Ncat); + Nlast = Ncat; + } + + if (DUMP) { + for (i = 0; i < Ncat; i++) { + fprintf (stdout, "%d %f %f %f\n", i, ra[i], dec[i], mag[i]); + } + exit (0); + } + + fprintf (stderr, "enter filenames, followed by EOF\n"); + while (fscanf (stdin, "%s", filename) != EOF) { + /**** See comment 1 below ****/ + strcpy (head, filename); + if ((c = strrchr(head, '.')) != NULL) + strcpy(strrchr(head, '.'), ".head"); + else + strcat(head, ".head"); + strcpy (obj, filename); + if ((c = strrchr(obj, '.')) != NULL) + strcpy(strrchr(obj, '.'), ".obj_out"); + else + strcat(obj, ".obj_out"); + + stats (head, &RA, &DEC, &dRA, &dDEC, ASEC_PIX); + if (FORCE) + stars1 = get_catalog_stars (ra, dec, mag, Ncat, F_RA, F_DEC, dRA, dDEC, &N1); + else { + stars1 = get_catalog_stars (ra, dec, mag, Ncat, RA, DEC, dRA, dDEC, &N1); + } + N2 = N2NUMBER; + stars2 = get_stars (obj, &N2); + if (N2 >= 3) { + if (FORCE) + Rx = ASEC_PIX / (3600 * cos (F_DEC*RAD_DEG)); + else + Rx = ASEC_PIX / (3600 * cos (DEC*RAD_DEG)); + Ry = ASEC_PIX / 3600; + delta (stars1, stars2, N1, N2, &dx, &dy, &Rx, &Ry); + fprintf (stderr, "USING: %f %f %f %f\n", dx, dy, Rx, Ry); + Nused = line_fit (stars1, stars2, N1, N2, dx, dy, Rx, Ry, + &X_O, &X_X, &X_Y, &dX, &Y_O, &Y_X, &Y_Y, &dY); + alter_header (head, argv, argc, X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY, Nused); + } + else { + fprintf (stderr, "too few stars in %s, skipped\n", filename); + } + free (stars2); + free (stars1); + } +} + + + + + /* comment 1: + assign the assumed names: + filename.head, filename.obj_out + if filename is of the form xxxx.xxx, cut off + extention before attaching .head, .obj_out + this entire problem can be removed by making the + Dophot output file contain the header info from + the image, but until this is the Dophot output + format... */ Index: /branches/elixir/Ohana/src/gastro/src/config.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/config.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/config.c (revision 2434) @@ -0,0 +1,110 @@ +# include +# define D_NBYTES 4096 + +char *LoadConfigFile (filename) +char *filename; +{ + + FILE *f; + int i, done, Nbytes, NBYTES, size; + char *config; + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "couldn't find %s\n", filename); + return ((char *) NULL); + } + + NBYTES = D_NBYTES; + ALLOCATE (config, char, NBYTES); + + size = 0; + done = FALSE; + for (i = 0; !done; i++) { + Nbytes = fread (&config[i*D_NBYTES], sizeof(char), D_NBYTES, f); + size += Nbytes; + if (Nbytes < D_NBYTES) + done = TRUE; + else { + NBYTES += D_NBYTES; + REALLOCATE (config, char, NBYTES); + } + } + + config[size] = '\n'; + config[size+1] = 0; + REALLOCATE (config, char, size + 2); + + fclose (f); + + return (config); +} + +int ScanConfig /* we expect one more field: the pointer to the value requested */ +# ifndef ANSI +(config, field, mode, va_alist) + char *config; char *field, *mode; va_dcl +# else +(char config[], + char field[], + char mode[],...) +# endif +{ + + int i, j; + char *p, *p2, tmp[256]; + va_list argp; + double value; + +# ifndef ANSI + va_start (argp); +# else + va_start (argp, N); +# endif + + /* find the correct line with field */ + p2 = config; + for (i = 0; ; i++) { + if (!strncmp (field, p2, strlen(field))) { + p = p2 + strlen (field); + break; + } + else { + p2 = strchr (p2, '\n'); + if (p2 == (char *) NULL) { + fprintf (stderr, "entry %s not found in config file\n", field); + return (FALSE); + } + else p2++; + } + } + if (!strcmp (mode, "%s")) { + p2 = strchr (p, '\n'); + if (p2 == (char *) NULL) + p2 = config + strlen(config); + bcopy (p, tmp, (p2-p)); + tmp[(p2-p)] = 0; + stripwhite (tmp); + p2 = va_arg (argp, char *); + strcpy (p2, tmp); + } + else { + + /* try to get a numerical value from the field */ + value = strtod (p, &p2); + if ((*p2 == 'd') || (*p2 == 'D')) + value *= pow (10.0, atof (p2 + 1)); + + if (!strcmp (mode, "%d")) *va_arg (argp, int *) = value; + if (!strcmp (mode, "%u")) *va_arg (argp, unsigned *) = value; + if (!strcmp (mode, "%ld")) *va_arg (argp, long *) = value; + if (!strcmp (mode, "%hd")) *va_arg (argp, short *) = value; + if (!strcmp (mode, "%f")) *va_arg (argp, float *) = value; + if (!strcmp (mode, "%lf")) *va_arg (argp, double *) = value; + + } + + va_end (argp); + return (TRUE); + +} Index: /branches/elixir/Ohana/src/gastro/src/coordops.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/coordops.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/coordops.c (revision 2434) @@ -0,0 +1,362 @@ +# include "gastro.h" + +XY_to_RD (ra, dec, x, y, coords) +double *ra, *dec; +double x, y; +Coords coords[]; +{ + + double L, M, X, Y, T, Z; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, sdp, cdp; + + *ra = 0; + *dec = 0; + stht = ctht = 1; + + /** convert pixel coordinates to cartesian system **/ + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + if (!strcmp(&coords[0].ctype[4], "-PLY")) { + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]); + Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]); + } + if (coords[0].Npolyterms > 2) { + X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]); + Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]); + } + } + + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + /* in FITS ref, L,M = x, y alpha, delta = phi, theta */ + + /**** Zenithal Projections ****/ + if (!strcmp(&coords[0].ctype[4], "-PLY") || + !strcmp(&coords[0].ctype[4], "-TAN") || + !strcmp(&coords[0].ctype[4], "-SIN") || + !strcmp(&coords[0].ctype[4], "-ZEA") || + !strcmp(&coords[0].ctype[0], "MM")) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } + else { + sphi = L / R; + cphi = -M / R; + } + + if (!strcmp(&coords[0].ctype[4], "-PLY") || + !strcmp(&coords[0].ctype[4], "-TAN")) { + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } + else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || + !strcmp(&coords[0].ctype[0], "MM")) { + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + } + if (!strcmp(&coords[0].ctype[4], "-ZEA")) { + stht = 1 - 0.5*SQ(R*RAD_DEG); + ctht = sqrt (1 - stht*stht); + } + + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + } + + /**** Locally Cartesian Projections ****/ + if (!strcmp(&coords[0].ctype[4], "-LIN") || + !strcmp(&coords[0].ctype[0], "GENE")) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + } + + /**** Other Conventional Projections ****/ + if (!strcmp(&coords[0].ctype[4], "-AIT")) { + Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M)); + alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 1.0); + delta = DEG_RAD * asin (RAD_DEG*M*Z); + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } + if (!strcmp(&coords[0].ctype[4], "-GLS")) { + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / cos (RAD_DEG * M); + delta = M; + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } + if (!strcmp(&coords[0].ctype[4], "-PAR")) { + /* L,M in degrees, alpha,delta in degrees */ + alpha = L / (1.0 - SQ(2.0*M/180)); + delta = 3 * DEG_RAD * asin (M/180.0); + *ra = alpha + coords[0].crval1; + *dec = delta + coords[0].crval2; + } +} + +RD_to_XY (x, y, ra, dec, coords) +double *x, *y; +double ra, dec; +Coords coords[]; +{ + + double phi, theta; + double tmp_d; + double X, Y, sphi, cphi, stht; + double salp, calp, sdel, cdel, sdp, cdp; + double P, CP, A, Rc; + int status; + + if (!strcmp(&coords[0].ctype[4], "-PLY")) { + /* fprintf (stderr, "approximate to polynomial TAN plane fit\n"); */ + } + + X = Y = 1; + status = TRUE; + *x = 0; + *y = 0; + + /**** Locally Cartesian Projections ****/ + if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) { + X = (ra - coords[0].crval1); + Y = (dec - coords[0].crval2); + } + + /**** Zenithal Projections ****/ + if (!strcmp(&coords[0].ctype[4], "-PLY") || + !strcmp(&coords[0].ctype[4], "-TAN") || + !strcmp(&coords[0].ctype[4], "-ZEA") || + !strcmp(&coords[0].ctype[4], "-SIN") || + !strcmp(&coords[0].ctype[0], "MM")) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; /* sin(theta) */ + sphi = cdel*salp; /* = cos(theta)*sin(phi) */ + cphi = cdel*sdp*calp - sdel*cdp; /* = cos(theta)*cos(phi) */ + if (stht < 0) status = FALSE; + + if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) { + X = DEG_RAD * sphi / stht; + Y = -DEG_RAD * cphi / stht; + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + X = DEG_RAD * sphi; + Y = -DEG_RAD * cphi; + } + if (!strcmp(&coords[0].ctype[4], "-ZEA")) { + Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht); + X = Rc * sphi; + Y = -Rc * cphi; + status = TRUE; + } + } + + /**** Other Standard Projections ****/ + if (!strcmp(&coords[0].ctype[4], "-AIT")) { + phi = RAD_DEG*(ra - coords[0].crval1); + theta = RAD_DEG*(dec - coords[0].crval2); + P = 1.0 + cos (theta) * cos (0.5*phi); + if (P != 0.0) { + A = DEG_RAD * sqrt (2.0 / P); + X = 2.0 * A * cos (theta) * sin (0.5*phi); + Y = A * sin (theta); + } else { + X = 0.0; + Y = 0.0; + } + } + if (!strcmp(&coords[0].ctype[4], "-GLS")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + X = phi * cos(RAD_DEG * theta); + Y = theta; + } + if (!strcmp(&coords[0].ctype[4], "-PAR")) { + phi = ra - coords[0].crval1; + theta = dec - coords[0].crval2; + X = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1); + Y = 180.0 * sin (RAD_DEG*theta/3.0); + } + + tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1; + *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2; + + return (status); + +} + +fRD_to_XY (x, y, ra, dec, coords) +float *x, *y; +double ra, dec; +Coords coords[]; +{ + + int status; + double tmpx, tmpy; + + status = RD_to_XY (&tmpx, &tmpy, ra, dec, coords); + *x = tmpx; + *y = tmpy; + + return (status); + +} + +fXY_to_RD (ra, dec, x, y, coords) +float *ra, *dec; +double x, y; +Coords coords[]; +{ + + int status; + double tmpr, tmpd; + + status = XY_to_RD (&tmpr, &tmpd, x, y, coords); + *ra = tmpr; + *dec = tmpd; + + return (status); + +} + +int GetCoords (coords, header) +Coords coords[]; +Header header[]; +{ + + int status; + double rotate; + + rotate = 0.0; + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + + status = TRUE; + if (fits_scan (header, "CTYPE2", "%s", 1, coords[0].ctype)) { + status = fits_scan (header, "CRVAL1", "%lf", 1, &coords[0].crval1); + status &= fits_scan (header, "CRPIX1", "%lf", 1, &coords[0].crpix1); + status &= fits_scan (header, "CRVAL2", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "CRPIX2", "%lf", 1, &coords[0].crpix2); + if (fits_scan (header, "CDELT1", "%lf", 1, &coords[0].cdelt1)) { + status &= fits_scan (header, "CDELT2", "%lf", 1, &coords[0].cdelt2); + if (fits_scan (header, "CROTA2", "%lf", 1, &rotate)) { + coords[0].pc1_1 = cos(rotate*RAD_DEG); + coords[0].pc1_2 = -sin(rotate*RAD_DEG); + coords[0].pc2_1 = sin(rotate*RAD_DEG); + coords[0].pc2_2 = cos(rotate*RAD_DEG); + } + if (fits_scan (header, "PC001001", "%lf", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "PC001002", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "PC002001", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "PC002002", "%lf", 1, &coords[0].pc2_2); + } + if (!strcmp (coords[0].ctype, "RA---PLY")) { + if (!fits_scan (header, "NPLYTERM", "%d", 1, &coords[0].Npolyterms)) { + coords[0].Npolyterms = 3; + } + switch (coords[0].Npolyterms) { + case 3: + status &= fits_scan (&header, "PCA1X3Y0", "%le", 1, &coords.polyterms[3][0]); + status &= fits_scan (&header, "PCA1X2Y1", "%le", 1, &coords.polyterms[4][0]); + status &= fits_scan (&header, "PCA1X1Y2", "%le", 1, &coords.polyterms[5][0]); + status &= fits_scan (&header, "PCA1X0Y3", "%le", 1, &coords.polyterms[6][0]); + status &= fits_scan (&header, "PCA2X3Y0", "%le", 1, &coords.polyterms[3][1]); + status &= fits_scan (&header, "PCA2X2Y1", "%le", 1, &coords.polyterms[4][1]); + status &= fits_scan (&header, "PCA2X1Y2", "%le", 1, &coords.polyterms[5][1]); + status &= fits_scan (&header, "PCA2X0Y3", "%le", 1, &coords.polyterms[6][1]); + case 2: + status &= fits_scan (&header, "PCA1X2Y0", "%le", 1, &coords.polyterms[0][0]); + status &= fits_scan (&header, "PCA1X1Y1", "%le", 1, &coords.polyterms[1][0]); + status &= fits_scan (&header, "PCA1X0Y2", "%le", 1, &coords.polyterms[2][0]); + status &= fits_scan (&header, "PCA2X2Y0", "%le", 1, &coords.polyterms[0][1]); + status &= fits_scan (&header, "PCA2X1Y1", "%le", 1, &coords.polyterms[1][1]); + status &= fits_scan (&header, "PCA2X0Y2", "%le", 1, &coords.polyterms[2][1]); + case 0: + case 1: + } + } + else { + if (fits_scan (header, "CD1_1", "%lf", 1, &coords[0].pc1_1)) { + status &= fits_scan (header, "CD1_2", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "CD2_1", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "CD2_2", "%lf", 1, &coords[0].pc2_2); + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + } + else { + status = FALSE; + } + } + } + else { + if (fits_scan (header, "RA_O", "%lf", 1, &coords[0].crval1)) { + status = fits_scan (header, "RA_X", "%lf", 1, &coords[0].pc1_1); + status &= fits_scan (header, "RA_Y", "%lf", 1, &coords[0].pc1_2); + status &= fits_scan (header, "DEC_O", "%lf", 1, &coords[0].crval2); + status &= fits_scan (header, "DEC_X", "%lf", 1, &coords[0].pc2_1); + status &= fits_scan (header, "DEC_Y", "%lf", 1, &coords[0].pc2_2); + coords[0].crpix1 = coords[0].crpix2 = 0.0; + coords[0].cdelt1 = coords[0].cdelt2 = 1.0; + strcpy (coords[0].ctype, "GENE"); + } + } + if (!status) { + fprintf (stderr, "error getting all elements for coordinate mode %d\n", coords[0].ctype); + coords[0].crval1 = coords[0].crpix1 = coords[0].cdelt1 = 0.0; + coords[0].crval2 = coords[0].crpix2 = coords[0].cdelt2 = 0.0; + coords[0].pc1_1 = coords[0].pc2_2 = 1.0; + coords[0].pc2_1 = coords[0].pc1_2 = 0.0; + strcpy (coords[0].ctype, "NONE"); + } + return (status); +} + + + /* -PLY projection is an extrapolation of the -TAN projection. + In addition to the usual linear terms of CRPIXi, PC00i00j, there are + higher order polynomial terms (up to 3rd order): + Axis 1 terms: + PCA1X2Y0 = coords.polyterm[0][0] = x^2 + PCA1X1Y1 = coords.polyterm[1][0] = xy + PCA1X0Y2 = coords.polyterm[2][0] = y^2 + PCA1X3Y0 = coords.polyterm[3][0] = x^3 + PCA1X2Y1 = coords.polyterm[4][0] = x^2 y + PCA1X1Y2 = coords.polyterm[5][0] = x y^2 + PCA1X0Y3 = coords.polyterm[6][0] = y^2 + Axis 2 terms: + PCA2X2Y0 = coords.polyterm[0][1] = x^2 + PCA2X1Y1 = coords.polyterm[1][1] = xy + PCA2X0Y2 = coords.polyterm[2][1] = y^2 + PCA2X3Y0 = coords.polyterm[3][1] = x^3 + PCA2X2Y1 = coords.polyterm[4][1] = x^2 y + PCA2X1Y2 = coords.polyterm[5][1] = x y^2 + PCA2X0Y3 = coords.polyterm[6][1] = y^2 + */ Index: /branches/elixir/Ohana/src/gastro/src/dastro.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/dastro.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/dastro.c (revision 2434) @@ -0,0 +1,86 @@ +# include "astro.h" + +main (argc, argv) +int argc; +char **argv; +{ + + int status; + Header header; + char filename[1024], reffile[1024]; + double X_O, X_X, X_Y, Y_O, Y_X, Y_Y; + double RA_O, RA_X, RA_Y, DEC_O, DEC_X, DEC_Y; + double ra_o, ra_x, ra_y, dec_o, dec_x, dec_y; + double dRA, dDEC, dX, dY, dra, ddec; + + while (fscanf (stdin, "%s", filename) != EOF) { + + status = fits_read_header (filename, &header); + if (!status) { + fprintf (stderr, "error opening file %s", filename); + continue; + } + status = fits_scan (&header, "rREF", "%s", 1, reffile); + status &= fits_scan (&header, "X_O", "%lf", 1, &X_O); + status &= fits_scan (&header, "X_X", "%lf", 1, &X_X); + status &= fits_scan (&header, "X_Y", "%lf", 1, &X_Y); + status &= fits_scan (&header, "Y_O", "%lf", 1, &Y_O); + status &= fits_scan (&header, "Y_X", "%lf", 1, &Y_X); + status &= fits_scan (&header, "Y_Y", "%lf", 1, &Y_Y); + status &= fits_scan (&header, "dX", "%lf", 1, &dX); + status &= fits_scan (&header, "dY", "%lf", 1, &dY); + if (!status) { + fprintf (stderr, "file missing rastro info: %s", filename); + fits_free_header (&header); + continue; + } + + status = fits_read_header (reffilename, &refhead); + if (!status) { + fprintf (stderr, "error opening file %s", reffilename); + continue; + } + status &= fits_scan (&refhead, "RA_O", "%lf", 1, &RA_O); + status &= fits_scan (&refhead, "RA_X", "%lf", 1, &RA_X); + status &= fits_scan (&refhead, "RA_Y", "%lf", 1, &RA_Y); + status &= fits_scan (&refhead, "DEC_O", "%lf", 1, &DEC_O); + status &= fits_scan (&refhead, "DEC_X", "%lf", 1, &DEC_X); + status &= fits_scan (&refhead, "DEC_Y", "%lf", 1, &DEC_Y); + status &= fits_scan (&refhead, "dRA", "%lf", 1, &dRA); + status &= fits_scan (&refhead, "dDEC", "%lf", 1, &dDEC); + if (!status) { + fprintf (stderr, "file missing rastro info: %s", reffile); + fits_free_header (&header); + fits_free_header (&refhead); + continue; + } + + + ra_o = RA_X*X_O + RA_Y*Y_O + RA_O; + dec_o = DEC_X*X_O + DEC_Y*Y_O + DEC_O; + + ra_x = RA_X*X_X + RA_Y*Y_X; + ra_y = RA_X*X_Y + RA_Y*Y_Y; + dec_x = DEC_X*X_X + DEC_Y*Y_X; + dec_y = DEC_X*X_Y + DEC_Y*Y_Y; + + dra = sqrt(dRA*dRA + SQ(3600*dX*RA_X) + SQ(3600*dY*RA_Y)); + ddec = sqrt(dDEC*dDEC + SQ(3600*dX*DEC_X) + SQ(3600*dY*DEC_Y)); + + status = fits_modify (&header, "rREF", "%s", 1, reffile); + status &= fits_modify (&header, "RA_O", "%lf", 1, RA_O); + status &= fits_modify (&header, "RA_X", "%lf", 1, RA_X); + status &= fits_modify (&header, "RA_Y", "%lf", 1, RA_Y); + status &= fits_modify (&header, "DEC_O", "%lf", 1, DEC_O); + status &= fits_modify (&header, "DEC_X", "%lf", 1, DEC_X); + status &= fits_modify (&header, "DEC_Y", "%lf", 1, DEC_Y); + status &= fits_modify (&header, "dRA", "%lf", 1, dRA); + status &= fits_modify (&header, "dDEC", "%lf", 1, dDEC); + + fits_free_header (&header); + fits_free_header (&refhead); + + } + +} + Index: /branches/elixir/Ohana/src/gastro/src/delta.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/delta.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/delta.c (revision 2434) @@ -0,0 +1,53 @@ +# include "astro.h" +# define FIND_SHIFT(A,B) \ + find_shift (stars1, stars2, N1, N2, (A), (B), &N, &S, &dX, &dY); \ + if (S < s) { \ + s = S; \ + *Rx = (A); \ + *Ry = (B); \ + *dx = dX; \ + *dy = dY; \ + n = N; \ + } + +void delta (stars1, stars2, N1, N2, dx, dy, Rx, Ry) +Stars stars1[], stars2[]; +int N1, N2; +double *dx, *dy, *Rx, *Ry; +{ + + int N, n; + double dX, dY, RX, RY, s, S; + + n = 0; + s = 10000; + RX = *Rx; + RY = *Ry; + + if (!FLIP || (FLIP == 1)) { +/* FIND_SHIFT( 0.9*RX, 0.9*RY); */ + FIND_SHIFT( 1.0*RX, 1.0*RY); +/* FIND_SHIFT( 1.1*RX, 1.1*RY); */ + } + + if (!FLIP || (FLIP == 2)) { +/* FIND_SHIFT(-0.9*RX, 0.9*RY); */ + FIND_SHIFT(-1.0*RX, 1.0*RY); +/* FIND_SHIFT(-1.1*RX, 1.1*RY); */ + } + + if (!FLIP || (FLIP == 3)) { +/* FIND_SHIFT(-0.9*RX, -0.9*RY); */ + FIND_SHIFT(-1.0*RX, -1.0*RY); +/* FIND_SHIFT(-1.1*RX, -1.1*RY); */ + } + + if (!FLIP || (FLIP == 4)) { +/* FIND_SHIFT( 0.9*RX, -0.9*RY); */ + FIND_SHIFT( 1.0*RX, -1.0*RY); +/* FIND_SHIFT( 1.1*RX, -1.1*RY); */ + } + + fprintf (stderr, "USING: %f %f %f %f\n", *dx, *dy, *Rx, *Ry); + +} Index: /branches/elixir/Ohana/src/gastro/src/delta.number =================================================================== --- /branches/elixir/Ohana/src/gastro/src/delta.number (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/delta.number (revision 2434) @@ -0,0 +1,41 @@ +# include "astro.h" +# define FIND_SHIFT(A,B) \ + find_shift (stars1, stars2, N1, N2, (A), (B), &N, &dX, &dY); \ + if (N > n) { \ + *Rx = (A); \ + *Ry = (B); \ + *dx = dX; \ + *dy = dY; \ + n = N; \ + } + +void delta (stars1, stars2, N1, N2, dx, dy, Rx, Ry) +Stars stars1[], stars2[]; +int N1, N2; +double *dx, *dy, *Rx, *Ry; +{ + + int N, n; + double dX, dY, RX, RY; + + n = 0; + RX = *Rx; + RY = *Ry; + FIND_SHIFT( 0.9*RX, 0.9*RY); + FIND_SHIFT(-0.9*RX, 0.9*RY); + FIND_SHIFT( 0.9*RX, -0.9*RY); + FIND_SHIFT(-0.9*RX, -0.9*RY); + + FIND_SHIFT( 1.0*RX, 1.0*RY); + FIND_SHIFT(-1.0*RX, 1.0*RY); + FIND_SHIFT( 1.0*RX, -1.0*RY); + FIND_SHIFT(-1.0*RX, -1.0*RY); + + FIND_SHIFT( 1.1*RX, 1.1*RY); + FIND_SHIFT(-1.1*RX, 1.1*RY); + FIND_SHIFT( 1.1*RX, -1.1*RY); + FIND_SHIFT(-1.1*RX, -1.1*RY); + + fprintf (stderr, "USING: %f %f %f %f\n", *dx, *dy, *Rx, *Ry); + +} Index: /branches/elixir/Ohana/src/gastro/src/extract.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/extract.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/extract.c (revision 2434) @@ -0,0 +1,36 @@ +# include "astro.h" +# define MMIN 14.5 + +main (argc, argv) +int argc; +char **argv; +{ + + int i, N1, N2, Ncat, Nlast, Nused; + double dx, dy, Rx, Ry; + double RA, DEC, dRA, dDEC; + double *ra, *dec, *mag; + double X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY; + char *c, *ref, filename[200], head[200], obj[200], line[200]; + Stars *stars1, *stars2; + + RA = 500.0; DEC = 500.0; + + if (argc < 2) { + fprintf (stderr, "USAGE: extract (gscfile) \n"); + exit (0); + } + + Ncat = 0.0; + extract_catalog (argv[1], &ra, &dec, &mag, &Ncat); + + for (i = 0; i < Ncat; i++) { + if ((hypot (ra[i] - RA, dec[i] - DEC) > 0.0006) && (mag[i] < MMIN)) { + fprintf (stdout, "%8.4f %8.4f %4.1f\n", ra[i], dec[i], mag[i]); + } + RA = ra[i]; DEC = dec[i]; + } + +} + + Index: /branches/elixir/Ohana/src/gastro/src/extract_catalog.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/extract_catalog.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/extract_catalog.c (revision 2434) @@ -0,0 +1,179 @@ +# include "astro.h" +# define NBYTES 160000 + +void extract_catalog (filename, ra, dec, mag, Nstars) +char *filename; +double **ra; +double **dec; +double **mag; +int *Nstars; +{ + + Header header; + Table table; + char **temp1; + char origin[300]; + int i; + double *temp2, R, D, minM, maxM; + double minR, maxR, minD, maxD; + + if (*Nstars == 0) { + ALLOCATE (*ra, double, 1); + ALLOCATE (*dec, double, 1); + ALLOCATE (*mag, double, 1); + } + + fits_read_header (filename, &header); + fits_scan (&header, "ORIGIN", "%s", 1, origin); + fits_free_header (&header); + + minM = 100; + maxM = 0; + minR = 390; + maxR = -10; + minD = 100; + maxD = -100; + if (!strcmp (origin, "ST ScI")) { /* Assuming it is guide star catalog.. */ + fits_read_table (filename, &table); + fits_table_column (&table, "RA_DEG", "%lf", &temp2); + REALLOCATE (*ra, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*ra)[*Nstars+i] = temp2[i]; + minR = MIN (temp2[i], minR); + maxR = MAX (temp2[i], maxR); + } + free (temp2); + fits_table_column (&table, "DEC_DEG", "%lf", &temp2); + REALLOCATE (*dec, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*dec)[*Nstars+i] = temp2[i]; + minD = MIN (temp2[i], minD); + maxD = MAX (temp2[i], maxD); + } + free (temp2); + fits_table_column (&table, "MAG", "%lf", &temp2); + REALLOCATE (*mag, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*mag)[*Nstars+i] = temp2[i]; + minM = MIN (temp2[i], minM); + maxM = MAX (temp2[i], maxM); + } + *Nstars += table.Naxis[1]; + free (temp2); + fprintf (stderr, "mag range: %f -- %f\n", minM, maxM); + fprintf (stderr, "RA range: %f -- %f, \n", minR, maxR); + fprintf (stderr, "DEC range: %f -- %f, \n", minD, maxD); + return; + } + + if (!strncmp (origin, "U.Minn. PDP11/60", strlen("U.Minn. PDP11/60"))) { /* Assuming it is Berkhuijsen etal catalog.. */ + fits_read_table (filename, &table); + fits_table_column (&table, "RA(1950)", "%s", &temp1); + REALLOCATE (*ra, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*ra)[*Nstars+i] = 15*dms_to_deg (temp1[i], ':'); + minR = MIN ((*ra)[*Nstars+i], minR); + maxR = MAX ((*ra)[*Nstars+i], maxR); + free (temp1[i]); + } + free (temp1); + fits_table_column (&table, "DEC(1950)", "%s", &temp1); + REALLOCATE (*dec, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*dec)[*Nstars+i] = dms_to_deg (temp1[i], ':'); + minD = MIN ((*dec)[*Nstars+i], minD); + maxD = MAX ((*dec)[*Nstars+i], maxD); + free (temp1[i]); + } + free (temp1); + fits_table_column (&table, "V-MAG", "%lf", &temp2); + REALLOCATE (*mag, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*mag)[*Nstars+i] = temp2[i]; + minM = MIN (temp2[i], minM); + maxM = MAX (temp2[i], maxM); + } + free (temp2); + for (i = *Nstars; i < *Nstars + table.Naxis[1]; i++) { + if (i < 20) + fprintf (stderr, "%lf %lf -> ", (*ra)[i], (*dec)[i]); + /* see Magnier et al. 1993 A&A 272 695 */ + R = 0.99999848*(*ra)[i] - 0.00152332*(*dec)[i] + 0.0614972; + D = 0.99999848*(*dec)[i] + 0.0009958* (*ra)[i] - 0.0093577; + precess (&R, &D, 1949.999790, 2000.0); /* B1950.0 -> J2000.0 */ + (*ra)[i] = R; + (*dec)[i] = D; + if (i < 20) + fprintf (stderr, "%lf %lf\n", (*ra)[i], (*dec)[i]); + } + *Nstars += table.Naxis[1]; + fprintf (stderr, "mag range: %f -- %f\n", minM, maxM); + return; + } + + if (!strncmp (origin, "MIT-CSR", strlen("MIT-CSR"))) { /* Assuming it is Magnier etal catalog.. */ + fits_read_table (filename, &table); + fits_table_column (&table, "RA", "%lf", &temp2); + REALLOCATE (*ra, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*ra)[*Nstars+i] = temp2[i]; + minR = MIN ((*ra)[*Nstars+i], minR); + maxR = MAX ((*ra)[*Nstars+i], maxR); + } + free (temp2); + fits_table_column (&table, "DEC", "%lf", &temp2); + REALLOCATE (*dec, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*dec)[*Nstars+i] = temp2[i]; + minD = MIN ((*dec)[*Nstars+i], minD); + maxD = MAX ((*dec)[*Nstars+i], maxD); + } + free (temp2); + fits_table_column (&table, "V-MAG", "%lf", &temp2); + REALLOCATE (*mag, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*mag)[*Nstars+i] = temp2[i]; + minM = MIN (temp2[i], minM); + maxM = MAX (temp2[i], maxM); + } + free (temp2); + *Nstars += table.Naxis[1]; + fprintf (stderr, "mag range: %f -- %f\n", minM, maxM); + return; + } + + if (!strncmp (origin, "RELPHOT", strlen("RELPHOT"))) { /* Assuming it is from relphot.c */ + fits_read_table (filename, &table); + fits_table_column (&table, "RA", "%lf", &temp2); + REALLOCATE (*ra, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*ra)[*Nstars+i] = temp2[i]; + minR = MIN ((*ra)[*Nstars+i], minR); + maxR = MAX ((*ra)[*Nstars+i], maxR); + } + free (temp2); + fits_table_column (&table, "DEC", "%lf", &temp2); + REALLOCATE (*dec, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*dec)[*Nstars+i] = temp2[i]; + minD = MIN ((*dec)[*Nstars+i], minD); + maxD = MAX ((*dec)[*Nstars+i], maxD); + } + free (temp2); + fits_table_column (&table, "MAG", "%lf", &temp2); + REALLOCATE (*mag, double, *Nstars + table.Naxis[1]); + for (i = 0; i < table.Naxis[1]; i++) { + (*mag)[*Nstars+i] = temp2[i]; + minM = MIN (temp2[i], minM); + maxM = MAX (temp2[i], maxM); + } + free (temp2); + *Nstars += table.Naxis[1]; + fprintf (stderr, "mag range: %f -- %f\n", minM, maxM); + return; + } + + fprintf (stderr, "unknown catalog format, ORIGIN: %s\n", origin); + exit (0); + +} Index: /branches/elixir/Ohana/src/gastro/src/find_shift.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/find_shift.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/find_shift.c (revision 2434) @@ -0,0 +1,87 @@ +# include "astro.h" +# define TIME {gettimeofday (&now, (struct timeval *) NULL); \ +fprintf (stderr, "start time: %f\n", now.tv_sec+now.tv_usec*1e-6);} + +void find_shift (stars1, stars2, N1, N2, Rx, Ry, n, S, dx, dy) +Stars stars1[], stars2[]; +int N1, N2; +double Rx, Ry; +double *dx, *dy; +double *S; +int *n; +{ + + double *DX, *DX2, *DY, *DY2, *N; + double Xrange, Yrange, Xzero, Yzero; + double dX, dY, sx, sy, s, dN, mean, sigma; + double *buf; + int i, j, dxi, dyi; + struct timeval now; + double fx1, fx2, fy1, fy2; + int k; + + *S = 1000000; + *n = *dx = *dy = 0; + NPIX = 1 + 30 * sqrt(N1*N2) / MIN (N1, N2); + mean = 2*N1*N2 / (NPIX*NPIX); + sigma = sqrt (mean); + ALLOCATE (N, double, NPIX*NPIX); + ALLOCATE (DX, double, NPIX*NPIX); + ALLOCATE (DY, double, NPIX*NPIX); + ALLOCATE (DX2, double, NPIX*NPIX); + ALLOCATE (DY2, double, NPIX*NPIX); + bzero (N, NPIX*NPIX*sizeof(double)); + bzero (DX, NPIX*NPIX*sizeof(double)); + bzero (DY, NPIX*NPIX*sizeof(double)); + bzero (DX2, NPIX*NPIX*sizeof(double)); + bzero (DY2, NPIX*NPIX*sizeof(double)); + ranges (stars1, stars2, N1, N2, &Xrange, &Yrange, &Xzero, &Yzero, Rx, Ry); + fx1 = (NPIX - 1) / Xrange; + fx2 = -fx1 * Xzero; + fy1 = (NPIX - 1) / Yrange; + fy2 = -fy1 * Yzero; + for (i = 0; i < N1; i++) { + for (j = 0; j < N2; j++) { + dX = stars1[i].X - stars2[j].X; + dY = stars1[i].Y - stars2[j].Y; + dxi = fx1*dX + fx2; + dyi = fy1*dY + fy2; + k = dxi*NPIX + dyi; + N[k] += 1.0; + DX[k] += dX; + DY[k] += dY; + DX2[k] += dX*dX; + DY2[k] += dY*dY; + } + } + + for (i = 0; i < NPIX; i++) { + for (j = 0; j < NPIX; j++) { + if (N[i*NPIX + j] < 4) + continue; + if (N[i*NPIX + j] - mean < 5*sigma) + continue; + + sx = DX2[i*NPIX + j] / (N[i*NPIX + j]) - SQ(DX[i*NPIX + j]) / SQ(N[i*NPIX + j]); + sy = DY2[i*NPIX + j] / (N[i*NPIX + j]) - SQ(DY[i*NPIX + j]) / SQ(N[i*NPIX + j]); + /* sx, sy are the positional scatter in the bin */ + s = sqrt (sx + sy + 1e-11) / SQ(N[i*NPIX + j] - 3); + if (N[i*NPIX + j] > mean + 10*sigma) { + if (s < *S) { + *S = s; + *n = N[i*NPIX + j]; + *dx = DX[i*NPIX + j] / (double) N[i*NPIX + j]; + *dy = DY[i*NPIX + j] / (double) N[i*NPIX + j]; + } + /* fprintf (stderr, "point: %f %f %f %f\n", s, N[i*NPIX + j], DX[i*NPIX + j] / (double) N[i*NPIX + j], DY[i*NPIX + j] / (double) N[i*NPIX + j]); */ + } + } + } + + free (N); + free (DX); + free (DY); + free (DX2); + free (DY2); + +} Index: /branches/elixir/Ohana/src/gastro/src/g2mass.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/g2mass.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/g2mass.c (revision 2434) @@ -0,0 +1,45 @@ +# include "gastro.h" + +/* filename has path to GSC, strip off path, use filename with CATDIR path */ + +int g2mass (char *fullpath, SStars **stars, int *Nstars) { + + int nstar, NSTARS; + float ra, dec, mag; + FILE *f; + + f = fopen (fullpath, "r"); + if (f == (FILE *) NULL) { + fprintf (stderr, "can't load catalog data\n"); + return (FALSE); + } + + /* for now, 2Mass data is just stored in a plain file with + ra, dec, mag */ + + nstar = *Nstars; + NSTARS = *Nstars + 1000; + if (*Nstars == 0) { + ALLOCATE (stars[0], SStars, NSTARS); + } else { + REALLOCATE (stars[0], SStars, NSTARS); + } + + while (fscanf (f, "%f %f %f", &ra, &dec, &mag) != EOF) { + stars[0][nstar].X = ra; + stars[0][nstar].Y = dec; + stars[0][nstar].mag = mag; + nstar ++; + if (nstar == NSTARS) { + NSTARS += 1000; + REALLOCATE (stars[0], SStars, NSTARS); + } + } + + fclose (f); + + *Nstars = nstar; + + return (TRUE); + +} Index: /branches/elixir/Ohana/src/gastro/src/gargs.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gargs.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gargs.c (revision 2434) @@ -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: /branches/elixir/Ohana/src/gastro/src/gastro.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gastro.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gastro.c (revision 2434) @@ -0,0 +1,84 @@ +# 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); + + 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: /branches/elixir/Ohana/src/gastro/src/gaussian.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gaussian.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gaussian.c (revision 2434) @@ -0,0 +1,60 @@ +# include +double drand48(); +double gaussian(); +double rnd_gauss(); + +# define NGAUSS 2048 +double gaussint[NGAUSS]; + +double rnd_gauss (double mean, double sigma) +{ + + double y; + + y = gaussint[(int)(NGAUSS*drand48())]; + + return (y*sigma + mean); + +} + +gauss_init () +{ + + int i; + long A, B; + double val, x, dx, dx1, dx2, dx3, df; + double mean, sigma; + + A = time(NULL); + for (B = 0; A == time(NULL); B++); + srand48(B); + + val = 0; + dx = 0.001; + dx1 = dx / 3.0; + dx2 = 2.0*dx/3.0; + dx3 = dx; + mean = 0.0; + sigma = 1.0; + + for (i = 0, x = -7.0; (i < NGAUSS) && (x < 7.0); x += dx) { + df = (3.0*gaussian(x , mean, sigma) + + 9.0*gaussian(x+dx1, mean, sigma) + + 9.0*gaussian(x+dx2, mean, sigma) + + 3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0); + val += df; + if (val > (i / (double) NGAUSS)) { + gaussint[i] = x + dx / 2.0; + i++; + } + } +} + +double gaussian (double x, double mean, double sigma) +{ + + double f, X; + f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma)); + return (f); + +} Index: /branches/elixir/Ohana/src/gastro/src/gaussj.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gaussj.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gaussj.c (revision 2434) @@ -0,0 +1,75 @@ +# 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; + + 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: /branches/elixir/Ohana/src/gastro/src/gcatalog.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gcatalog.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gcatalog.c (revision 2434) @@ -0,0 +1,43 @@ +# include "gastro.h" +# define BYTES_STAR 23 +# define BLOCK 1000 +# define DNSTARS 1000 + +int gcatalog (char *filename, SStars **stars, int *Nstars) { + + int i, NSTARS, nstar, Nbytes, nbytes; + char *buffer; + FILE *f; + + if (*Nstars == 0) ALLOCATE (*stars, SStars, 1); + nstar = *Nstars; + NSTARS = nstar + DNSTARS; + + f = fopen (filename, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find catalog file %s\n", filename); + return (FALSE); + } + + ALLOCATE (buffer, char, (BLOCK*BYTES_STAR)); + Nbytes = BLOCK*BYTES_STAR; + REALLOCATE (*stars, SStars, NSTARS); + 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 == NSTARS - 1) { + NSTARS += DNSTARS; + REALLOCATE (*stars, SStars, NSTARS); + } + } + } + free (buffer); + REALLOCATE (*stars, SStars, nstar); + *Nstars = nstar; + + return (TRUE); + +} Index: /branches/elixir/Ohana/src/gastro/src/gcenter.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gcenter.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gcenter.c (revision 2434) @@ -0,0 +1,253 @@ +# 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; + + if (PLOTSTUFF) { + Nvect = 0; + 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; + } + + + 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; + 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: /branches/elixir/Ohana/src/gastro/src/gdelta.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gdelta.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gdelta.c (revision 2434) @@ -0,0 +1,53 @@ +# include "astro.h" +# define FIND_SHIFT(A,B) \ + find_shift (stars1, stars2, N1, N2, (A), (B), &N, &S, &dX, &dY); \ + if (S < s) { \ + s = S; \ + *Rx = (A); \ + *Ry = (B); \ + *dx = dX; \ + *dy = dY; \ + n = N; \ + } + +void delta (stars1, stars2, N1, N2, dx, dy, Rx, Ry) +Stars stars1[], stars2[]; +int N1, N2; +double *dx, *dy, *Rx, *Ry; +{ + + int N, n; + double dX, dY, RX, RY, s, S; + + n = 0; + s = 10000; + RX = *Rx; + RY = *Ry; + + if (!FLIP || (FLIP == 1)) { +/* FIND_SHIFT( 0.9*RX, 0.9*RY); */ + FIND_SHIFT( 1.0*RX, 1.0*RY); +/* FIND_SHIFT( 1.1*RX, 1.1*RY); */ + } + + if (!FLIP || (FLIP == 2)) { +/* FIND_SHIFT(-0.9*RX, 0.9*RY); */ + FIND_SHIFT(-1.0*RX, 1.0*RY); +/* FIND_SHIFT(-1.1*RX, 1.1*RY); */ + } + + if (!FLIP || (FLIP == 3)) { +/* FIND_SHIFT(-0.9*RX, -0.9*RY); */ + FIND_SHIFT(-1.0*RX, -1.0*RY); +/* FIND_SHIFT(-1.1*RX, -1.1*RY); */ + } + + if (!FLIP || (FLIP == 4)) { +/* FIND_SHIFT( 0.9*RX, -0.9*RY); */ + FIND_SHIFT( 1.0*RX, -1.0*RY); +/* FIND_SHIFT( 1.1*RX, -1.1*RY); */ + } + + fprintf (stderr, "USING: %f %f %f %f\n", *dx, *dy, *Rx, *Ry); + +} Index: /branches/elixir/Ohana/src/gastro/src/get_argument.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/get_argument.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/get_argument.c (revision 2434) @@ -0,0 +1,19 @@ +# include "astro.h" + +int get_argument (argc, argv, arg) +int argc; +char **argv; +char *arg; +{ + + int i; + + for (i = 0; i < argc; i++) { + if (!strcmp(argv[i], arg)) + return (i); + } + + return ((int)NULL); +} + + Index: /branches/elixir/Ohana/src/gastro/src/get_catalog_stars.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/get_catalog_stars.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/get_catalog_stars.c (revision 2434) @@ -0,0 +1,40 @@ +# include "astro.h" +# define D_NSTARS 1000 + +Stars *get_catalog_stars (ra, dec, mag, Ncat, RA, DEC, dRA, dDEC, Nstars) +double *ra; +double *dec; +double *mag; +int Ncat; +double RA, DEC, dRA, dDEC; +int *Nstars; +{ + + int i, j, NSTARS; + Stars *stars; + + NSTARS = D_NSTARS; + ALLOCATE (stars, Stars, NSTARS); + fprintf (stderr,"%f %f %f %f\n", RA, DEC, dRA, dDEC); + + j = 0; + for (i = 0; i < Ncat; i++) { + if ((ra[i] > RA - NIM*dRA) && (ra[i] < RA + NIM*dRA) && + (dec[i] > DEC - NIM*dDEC) && (dec[i] < DEC + NIM*dDEC)) { + stars[j].X = ra[i]; + stars[j].Y = dec[i]; + stars[j].mag = mag[i]; + j++; + if (j == NSTARS - 1) { + NSTARS += D_NSTARS; + REALLOCATE (stars, Stars, NSTARS); + } + } + } + + *Nstars = j; + fprintf (stderr, "using %d stars from catalog\n", j); + + return (stars); + +} Index: /branches/elixir/Ohana/src/gastro/src/get_region_coords.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/get_region_coords.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/get_region_coords.c (revision 2434) @@ -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: /branches/elixir/Ohana/src/gastro/src/get_stars.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/get_stars.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/get_stars.c (revision 2434) @@ -0,0 +1,91 @@ +# include "astro.h" +# define D_NSTARS 1000 +# define NFIELD 4 + +Stars *get_stars (file, Nstars) +char file[]; +int *Nstars; +{ + + FILE *f; + int n, i, type, NSTARS; + Stars *stars; + double tmp; + + f = fopen (file, "r"); + if (f == NULL) { + fprintf (stderr, "could not find object file %s\n", file); + exit(0); + } + + NSTARS = D_NSTARS; + ALLOCATE (stars, Stars, NSTARS); + i = 0; + while ((n = fscanf (f, "%*f%d%lf%lf%lf%*f%*f%*f%*f%*f%*f%*f%*f%*f%*f", + &type, &stars[i].X, + &stars[i].Y, &stars[i].mag)) == NFIELD) { + if (ROTATE) { + tmp = stars[i].X; + stars[i].X = stars[i].Y; + stars[i].Y = tmp; + } + i++; + if (type != 1) + i--; + if (i == NSTARS - 1) { + NSTARS += D_NSTARS; + REALLOCATE (stars, Stars, NSTARS); + } + } + fclose(f); + if ((n != EOF) && (i < *Nstars)) + fprintf (stderr, "error reading in file %s\n", file); + fprintf (stderr, " %d stars read from %s, ", i, file); + NSTARS = i; + + sort_stars (stars, NSTARS); + + NSTARS = *Nstars = MIN (NSTARS, *Nstars); + REALLOCATE (stars, Stars, NSTARS); + fprintf (stderr, "using %5d\n", NSTARS); + + return (stars); +} + +sort_stars (stars, N) +Stars *stars; +int N; +{ + int l,j,ir,i; + int temp; + Stars 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; + } +} + + Index: /branches/elixir/Ohana/src/gastro/src/getusno.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/getusno.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/getusno.c (revision 2434) @@ -0,0 +1,131 @@ +# include "gastro.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}; + +USNOdata *getusno (USNOstats *usnostats, CatStats *catstats, int *Nusno) { + + 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; + USNOdata *usno; + + 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; + + 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 = 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: /branches/elixir/Ohana/src/gastro/src/gfit.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gfit.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gfit.c (revision 2434) @@ -0,0 +1,457 @@ +# 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; + + if (PLOTSTUFF) { + Nvect2 = 0; + NVECT2 = MAX(N1, N2); + Nvect = 0; + 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); + + /* 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; + + gaussj (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: /branches/elixir/Ohana/src/gastro/src/gfit.test.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gfit.test.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gfit.test.c (revision 2434) @@ -0,0 +1,496 @@ +# 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 (stars1, stars2, N1, N2, coords, NX, NY, Radius, DR, Nmatch, mode) + SStars stars1[], stars2[]; + int N1, N2, NX, NY, *Nmatch, mode; + double *Radius, *DR; + Coords *coords; +{ + + int i, j, iteration, Niter, last, halt, first_j, *Nextra, extras; + int *tmpN1, *tmpN2; + double dx, dy, Rx, Ry, dr; + 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, dX_Y, dY_X; + double RA, DEC, radius, radius2, scale, fratio; + double *tmpX1, *tmpX2, *tmpY1, *tmpY2; + double tX1, tX2, tY1, tY2; + double Rprev, Sprev, dSS, dS; + int *match; + double *dmatch; + + char c; + typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, color; + double lweight, size; + } Graphdata; + Graphdata graphdata; + float *xvect, *yvect; + float *xvect2, *yvect2; + int Nvect, NVECT; + int Nvect2, NVECT2; + + if (PLOTSTUFF) { + Nvect2 = 0; + NVECT2 = MAX(N1, N2); + Nvect = 0; + NVECT = MAX(N1, N2); + ALLOCATE (xvect, float, NVECT); + ALLOCATE (yvect, float, NVECT); + ALLOCATE (xvect2, float, NVECT2); + ALLOCATE (yvect2, float, NVECT2); + } + + + /* match arrays, init to -1 (no match) */ + ALLOCATE (match, int, N2); + ALLOCATE (dmatch, double, N2); + for (i = 0; i < N2; i++) match[i] = -1; + + /* 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); + + /* 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; + } + + /* 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, keep only if closest */ + if (DD < radius2) { + if (match[j] == -1) { + match[j] = i; + dmatch[j] = DD; + } + if (DD < dmatch[j]) { + match[j] = i; + dmatch[j] = DD; + } + } + } + j = first_j; + i++; + } + + /* measure statistics for matched stars */ + for (j = 0; j < N2; j++) { + if (match[j] == -1) continue; + i = match[j]; + + tX1 = tmpX1[i]; + tX2 = tmpX2[j]; + Dx = tX1 - tX2; + tY1 = tmpY1[i]; + tY2 = tmpY2[j]; + Dy = tY1 - tY2; + DD = Dx*Dx + Dy*Dy; + + 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); + } + } + + 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 Sx2, Sy2, Sxy, SXx, SXy, SYx, SYy; + 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; + + gaussj (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: /branches/elixir/Ohana/src/gastro/src/gfitpoly.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gfitpoly.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gfitpoly.c (revision 2434) @@ -0,0 +1,387 @@ +# 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 */ + 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); */ + gaussj (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: /branches/elixir/Ohana/src/gastro/src/gheader.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gheader.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gheader.c (revision 2434) @@ -0,0 +1,135 @@ +# 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 (!fits_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); + } + fits_modify (&header, "PHOTCODE", "%s", 1, PHOTCODE); + } + + if (Nmatch < 2) { + 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, Nmatch); + fits_modify (&header, "NASTRO", "%C", 1, "number of stars used for astrometry"); + + if (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, coords.cdelt1); + fits_modify (&header, "CDELT2", "%le", 1, coords.cdelt2); + fits_modify (&header, "CRVAL1", "%lf", 1, coords.crval1); + fits_modify (&header, "CRVAL2", "%lf", 1, coords.crval2); + fits_modify (&header, "CRPIX1", "%lf", 1, coords.crpix1); + fits_modify (&header, "CRPIX2", "%lf", 1, coords.crpix2); + fits_modify (&header, "PC001001", "%le", 1, coords.pc1_1); + fits_modify (&header, "PC001002", "%le", 1, coords.pc1_2); + fits_modify (&header, "PC002001", "%le", 1, coords.pc2_1); + fits_modify (&header, "PC002002", "%le", 1, coords.pc2_2); + fits_modify (&header, "NPLYTERM", "%d", 1, coords.Npolyterms); + if (coords.Npolyterms > 1) { + /* RA Terms */ + fits_modify (&header, "PCA1X2Y0", "%le", 1, coords.polyterms[0][0]); /* polyterms[0]); */ + fits_modify (&header, "PCA1X1Y1", "%le", 1, coords.polyterms[1][0]); /* polyterms[1]); */ + fits_modify (&header, "PCA1X0Y2", "%le", 1, coords.polyterms[2][0]); /* polyterms[2]); */ + + if (coords.Npolyterms > 2) { + fits_modify (&header, "PCA1X3Y0", "%le", 1, coords.polyterms[3][0]); /* polyterms[3]); */ + fits_modify (&header, "PCA1X2Y1", "%le", 1, coords.polyterms[4][0]); /* polyterms[4]); */ + fits_modify (&header, "PCA1X1Y2", "%le", 1, coords.polyterms[5][0]); /* polyterms[5]); */ + fits_modify (&header, "PCA1X0Y3", "%le", 1, coords.polyterms[6][0]); /* polyterms[6]); */ + } + /* Dec Terms */ + fits_modify (&header, "PCA2X2Y0", "%le", 1, coords.polyterms[0][1]); /* polyterms[7]); */ + fits_modify (&header, "PCA2X1Y1", "%le", 1, coords.polyterms[1][1]); /* polyterms[8]); */ + fits_modify (&header, "PCA2X0Y2", "%le", 1, coords.polyterms[2][1]); /* polyterms[9]); */ + + if (coords.Npolyterms > 2) { + fits_modify (&header, "PCA2X3Y0", "%le", 1, coords.polyterms[3][1]); /* polyterms[10]); */ + fits_modify (&header, "PCA2X2Y1", "%le", 1, coords.polyterms[4][1]); /* polyterms[11]); */ + fits_modify (&header, "PCA2X1Y2", "%le", 1, coords.polyterms[5][1]); /* polyterms[12]); */ + fits_modify (&header, "PCA2X0Y3", "%le", 1, coords.polyterms[6][1]); /* polyterms[13]); */ + } + } + 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*Nmatch)); + 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: /branches/elixir/Ohana/src/gastro/src/gproject.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gproject.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gproject.c (revision 2434) @@ -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: /branches/elixir/Ohana/src/gastro/src/gptolemy.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gptolemy.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gptolemy.c (revision 2434) @@ -0,0 +1,82 @@ +# include "gastro.h" + +/* filename has path to GSC, strip off path, use filename with CATDIR path */ + +int gptolemy (char *fullpath, SStars **stars, int *Nstars) { + + int i, Nsec, NSTARS, nstar; + char filename[256], *fileroot, *filepath, *subdir; + Catalog catalog; + + filepath = pathname (fullpath); + fileroot = filebasename (fullpath); + subdir = filebasename (filepath); + sprintf (filename, "%s/%s/%s", CATDIR, subdir, fileroot); + + catalog.filename = filename; + switch (lock_catalog (&catalog, LCK_SOFT)) { + case 0: + case 2: + fprintf (stderr, "can't load catalog data\n"); + return (FALSE); + case 1: + break; + } + if (!load_catalog (&catalog, LOAD_AVES | LOAD_MEAS, FALSE)) { + fprintf (stderr, "can't load catalog data\n"); + return (FALSE); + } + unlock_catalog (&catalog); + + nstar = *Nstars; + NSTARS = *Nstars + catalog.Naverage; + if (*Nstars == 0) { + ALLOCATE (stars[0], SStars, NSTARS); + } else { + REALLOCATE (stars[0], SStars, NSTARS); + } + + Nsec = catalog.Nsecfilt; + + /* select all entries, ignore magnitudes */ + for (i = 0; i < catalog.Naverage; i++, nstar++) { + stars[0][nstar].X = catalog.average[i].R; + stars[0][nstar].Y = catalog.average[i].D; + stars[0][nstar].mag = 0.001*catalog.measure[catalog.average[i].offset].M; + } + + /* select only entries with valid magnitudes */ +# if (0) + for (i = 0; i < catalog.Naverage; i++, nstar++) { + stars[0][nstar].X = catalog.average[i].R; + stars[0][nstar].Y = catalog.average[i].D; + GotIt = FALSE; + if (catalog.average[i].M != NO_MAG) { + stars[0][nstar].mag = 0.001*catalog.average[i].M; + GotIt = TRUE; + } + for (j = 0; !GotIt && (j < Nsec); j++) { + if (catalog.secfilt[i*Nsec + j].M != NO_MAG) { + stars[0][nstar].mag = 0.001*catalog.secfilt[i*Nsec + j].M; + GotIt = TRUE; + } + } + if (!GotIt) nstar --; /* skip this star */ + } +# endif + + free (catalog.average); + + *Nstars = nstar; + + return (TRUE); + +} + +/* this function is still only occasionally used. there is some trouble with the + issue of magnitudes. do we depend on the existence of an 'average.M' value, + or not? + + should this be selectable? + +*/ Index: /branches/elixir/Ohana/src/gastro/src/granges.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/granges.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/granges.c (revision 2434) @@ -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: /branches/elixir/Ohana/src/gastro/src/greference.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/greference.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/greference.c (revision 2434) @@ -0,0 +1,123 @@ +# include "gastro.h" + +# define USNO 1 +# define GSC 0 + +int greference (SStars **cat, int *Ncat, Coords *coords, int NX, int NY) { + + char **regions; + int i, Nregions, Nusno; + CatStats catstats; + USNOstats usnostats; + USNOdata *usno; + + if (VERBOSE) fprintf (stderr, "\nloading astrometric reference data from %s\n", REFCAT); + + *Ncat = 0; + + define_region (&catstats, coords, NX, NY); + + /* 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]) { + CatStats substats; + USNOdata *usno1, *usno2; + int Nusno1, Nusno2; + + substats = catstats; + substats.RA[0] = 0.0; + usno1 = getusno (&usnostats, &substats, &Nusno1); + + substats = catstats; + substats.RA[1] = 360.0; + usno2 = getusno (&usnostats, &substats, &Nusno2); + + ALLOCATE (*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 { + usno = getusno (&usnostats, &catstats, &Nusno); + + ALLOCATE (*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")) { + regions = gregions (&catstats, &Nregions); + for (i = 0; i < Nregions; i++) { + if (!gcatalog (regions[i], cat, Ncat)) + return (FALSE); + } + if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", *Ncat); + } + + if (!strcmp (REFCAT, "PTOLEMY")) { + regions = gregions (&catstats, &Nregions); + for (i = 0; i < Nregions; i++) { + if (!gptolemy (regions[i], cat, Ncat)) continue; + } + if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", *Ncat); + } + + if (!strcmp (REFCAT, "2MASS")) { + /* this needs to be written to handle the complete 2mass catalog. + at the moment, I only have a small subset */ + if (!g2mass ("/data/milo/eugene/elixir/refs/2mass/2mass.dat", cat, Ncat)) exit (1); + if (VERBOSE) fprintf (stderr, "%d stars from 2MASS\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: /branches/elixir/Ohana/src/gastro/src/gregions.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gregions.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gregions.c (revision 2434) @@ -0,0 +1,218 @@ +# include "gastro.h" +# define DEBUG 0 +double BigDecBounds[] = {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}; + +double Areas[] = {2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0, + 2700.0, 2665.4505, 2563.567, 2399.5344, 2181.6343, 1920.7299, 1629.5942, 1322.1251, 1012.5, 714.32555, 439.83856, 199.21043, 0.0}; + +char *DecSections[] = {"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"}; + +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"}; + +char *disk[] = {"disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", + "disk 1", "disk 1", "disk 1", "disk 1", "disk 1", "weirdness", + "disk 1", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", + "disk 2", "disk 2", "disk 2", "disk 2", "disk 2", "weirdness"}; + +int NBigRASections [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3, 48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 3}; + +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}; + +char **gregions (CatStats *catstats, int *Nregions) { + + char buffer[28800], temp[50], file[50], line[50]; + char *p, prefix[1024]; + char **regions; + FILE *f; + double dBigRA, RA, DEC, tmp, dr, dd, dec, ra, Ra; + double BigRABounds[2], RA0, RA1, DEC0, DEC1; + int i, NBigDec, NBigRA, NBig, VMS, N; + int NRA, NDEC, NLINES, done, NBLOCKS, nregion, NREGION; + + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE); + exit (0); + } + + /* create space to save names */ + nregion = 0; + NREGION = 50; + ALLOCATE (regions, char *, NREGION); + + /* smallest region files are 2.0 deg tall */ + /* uses arbitrary distance (3.0 deg) for width of image */ + dd = dr = 0.5; + for (Ra = catstats[0].RA[0]; Ra < catstats[0].RA[1] + dr; Ra += dr) { + for (dec = catstats[0].DEC[0]; dec < catstats[0].DEC[1] + dd; dec += dd) { + ra = Ra; + while (ra < 0) { ra += 360.0; } + while (ra >= 360.0) { ra -= 360.0; } + + /* use the pole regions, if near pole */ + if (dec > 86.25) { + sprintf (file, "%s/n8230/pole.cpt\0", GSC_DIR); + for (i = 0; i < nregion; i++) { + if (!strcmp (regions[i], file)) { + goto skip1; + } + } + regions[nregion] = strcreate (file); + nregion ++; + if (nregion == NREGION) { + NREGION += 50; + REALLOCATE (regions, char *, NREGION); + } + skip1: + continue; + } + + /* find large DEC region (directory) */ + NBigDec = -1; + for (i = 0; i < 12; i++) { + if ((dec >= BigDecBounds[i]) && (dec < BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + if (NBigDec < 0) { + for (i = 13; i < 24; i++) { + if ((dec < BigDecBounds[i]) && (dec >= BigDecBounds[i+1])) { + NBigDec = i; + break; + } + } + } + if (NBigDec < 0) { + fprintf (stderr, "ERROR: Dec out of range: %f\n", dec); + exit (0); + } + + /* count lines before section */ + NLINES = 0; + for (i = 0; i < NBigDec; i++) { + NLINES += NDecLines[i]; + } + fseek (f, 5*2880 + 48*NLINES, SEEK_SET); + + /* should be in this section. if not, there is a problem counting... */ + done = FALSE; + fread (buffer, 48*NDecLines[NBigDec], 1, f); + for (i = 0; !done && (i < NDecLines[NBigDec]); i++) { + strncpy (temp, &buffer[i*48], 48); + temp[49] = 0; + hms_to_deg (&RA0, &RA1, &DEC0, &DEC1, &temp[7]); + if (RA1 < RA0) RA1 += 360.0; + if ((dec >= 0) && (dec >= DEC0) && (dec < DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + if ((dec < 0) && (dec < DEC0) && (dec >= DEC1) && (ra >= RA0) && (ra < RA1)) { + done = TRUE; + } + } + if (!done) { + fprintf (stderr, "ERROR: can't find coords in region file: %f %f\n", ra, dec); + exit (0); + } + temp[5] = 0; + sprintf (file, "%s/%s/%s.cpt\0", GSC_DIR, Dec2Sections[NBigDec],&temp[1]); + for (i = 0; i < nregion; i++) { + if (!strcmp (regions[i], file)) { + goto skip2; + } + } + regions[nregion] = strcreate (file); + nregion ++; + if (nregion == NREGION) { + NREGION += 50; + REALLOCATE (regions, char *, NREGION); + } + skip2: + } + } + fclose (f); + + *Nregions = nregion; + return (regions); + +} + + + +/**********/ +int hms_to_deg (h0, h1, d0, d1, string) + char *string; + double *h0, *h1, *d0, *d1; +{ + + char *c; + int i, flag_d0, flag_d1, flag_h0, flag_h1; + double tmp; + + *d0 = *h0 = *d1 = *h1 = 0; + + flag_h0 = dparse (h0, 1, string); + flag_h1 = dparse (h1, 4, string); + flag_d0 = dparse (d0, 7, string); + flag_d1 = dparse (d1, 9, string); + *h0 *= flag_h0; + *h1 *= flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + dparse (&tmp, 2, string); + *h0 += tmp/60.0; + dparse (&tmp, 3, string); + *h0 += tmp/3600.0; + + dparse (&tmp, 5, string); + *h1 += tmp/60.0; + dparse (&tmp, 6, string); + *h1 += tmp/3600.0; + + dparse (&tmp, 8, string); + *d0 += tmp/60.0; + + dparse (&tmp, 10, string); + *d1 += tmp/60.0; + + *h0 *= 15*flag_h0; + *h1 *= 15*flag_h1; + *d0 *= flag_d0; + *d1 *= flag_d1; + + return (TRUE); +} + + +/* + if (buffer[i*48 + 19] == ' ') continue; + strncpy (temp, &buffer[i*48 + 19], 20); temp[20] = 0; + hms_to_deg (&RA, &DEC, temp, ' ', 3); + + + dBigRA = 360.0 / (int)(0.5 + 48*cos(3.1415927*0.5*(BigDecBounds[NBigDec] + BigDecBounds[NBigDec + 1])/180.0)); + + NBigRA = ra / dBigRA; + + NBig = NBigRA; + for (i = 0; (i < 12) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + for (i = 13; (i < 24) && (i < NBigDec); i++) { + NBig += NBigRASections[i]; + } + + + fprintf (stderr, "%d %d %d %d %f %f -> %f %f %f\n", NBigDec, NBigRA, NBigRASections[NBigDec], NBig, ra, dec, BigDecBounds[NBigDec], BigDecBounds[NBigDec + 1], dBigRA); +*/ Index: /branches/elixir/Ohana/src/gastro/src/gregions2.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gregions2.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gregions2.c (revision 2434) @@ -0,0 +1,242 @@ +# include "gastro.h" + +# 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 RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */ +char **gregions (CatStats *catstats, int *nnames) { + + CatalogRegion area; + double dtmp; + int i, Nnames; + char **names, file[256]; + + Nnames = 0; + ALLOCATE (names, char *, 1); + + /* force RA[i] to be in range 0 - 360 */ + while (catstats[0].RA[0] > 360.0) { catstats[0].RA[0] -= 360.0; } + while (catstats[0].RA[0] < 0.0) { catstats[0].RA[0] += 360.0; } + + while (catstats[0].RA[1] > 360.0) { catstats[0].RA[1] -= 360.0; } + while (catstats[0].RA[1] < 0.0) { catstats[0].RA[1] += 360.0; } + + /* force DEC[0] < DEC[1] */ + if (catstats[0].DEC[0] > catstats[0].DEC[1]) { + dtmp = catstats[0].DEC[0]; + catstats[0].DEC[0] = catstats[0].DEC[1]; + catstats[0].DEC[1] = dtmp; + } + + /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!) */ + if (catstats[0].DEC[1] > 86.25) { + sprintf (file, "n8230/pole"); + names[Nnames] = strcreate (file); + Nnames ++; + } + + /* if RA[0] > RA[1], split in two passes */ + if (catstats[0].RA[0] > catstats[0].RA[1]) { + area.RA[0] = 0.0; + area.RA[1] = catstats[0].RA[1]; + area.DEC[0] = catstats[0].DEC[0]; + area.DEC[1] = catstats[0].DEC[1]; + names = find_dec_bands (area, names, &Nnames); + + area.RA[0] = catstats[0].RA[0]; + area.RA[1] = 360.0; + area.DEC[0] = catstats[0].DEC[0]; + area.DEC[1] = catstats[0].DEC[1]; + names = find_dec_bands (area, names, &Nnames); + + } else { + + area.RA[0] = catstats[0].RA[0]; + area.RA[1] = catstats[0].RA[1]; + area.DEC[0] = catstats[0].DEC[0]; + area.DEC[1] = catstats[0].DEC[1]; + names = find_dec_bands (area, names, &Nnames); + } + + /* add the full path to the short names */ + for (i = 0; i < Nnames; i++) { + sprintf (file, "%s/%s.cpt", GSC_DIR, names[i]); + free (names[i]); + names[i] = strcreate (file); + } + + *nnames = Nnames; + return (names); + +} + +/* 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 */ + +/* must pass in an allocted pointer names, and Nnames contains the number so far */ +char **find_dec_bands (CatalogRegion area, char **names, int *Nnames) { + + int i; + + /* for NORTH region (DecBand[i] < DecBand[i+1]) */ + for (i = 0; i < NDECBANDS; i++) { + if ((area.DEC[1] > DecBand[i]) || (area.DEC[0] < DecBand[i+1])) { + names = load_ra_blocks (i, area, names, Nnames); + } + } + + /* for SOUTH region (DecBand[i] > DecBand[i+1]) */ + for (i = NDECBANDS+1; i < 2*NDECBANDS; i++) { + if ((area.DEC[0] < DecBand[i]) || (area.DEC[1] > DecBand[i+1])) { + names = load_ra_blocks (i, area, names, Nnames); + } + } + + return (names); +} + + +/* 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] */ +char **load_ra_blocks (int Ndec, CatalogRegion area, char **names, int *Ntnames) { + + CatalogRegion region; + int i, Nlines, Nnames; + char *buffer, file[256]; + FILE *f; + + /* open file */ + f = fopen (GSCFILE, "r"); + if (f == NULL) { + fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE); + exit (0); + } + + /* count lines before Dec Band */ + Nlines = 0; + for (i = 0; i < Ndec; i++) Nlines += NDecLines[i]; + + /* seek to start of Dec Band */ + fseek (f, 5*2880 + 48*Nlines, SEEK_SET); + + /* load data into buffer */ + Nlines = NDecLines[Ndec]; + ALLOCATE (buffer, char, 48*Nlines); + fread (buffer, 48*Nlines, 1, f); + fclose (f); + + /* must pass in a pointer which is already allocated */ + Nnames = *Ntnames; + REALLOCATE (names, char *, MAX (Nnames + Nlines, 1)); + + for (i = 0; i < Nlines; i++) { + parse_GSC_line (®ion, &buffer[i*48]); + if (area.RA[0] > region.RA[1]) continue; + if (area.RA[1] < region.RA[0]) continue; + if (area.DEC[0] > region.DEC[1]) continue; + if (area.DEC[1] < region.DEC[0]) continue; + sprintf (file, "%s/%s", Dec2Sections[Ndec], region.name); + names[Nnames] = strcreate (file); + Nnames ++; + } + + free (buffer); + REALLOCATE (names, char *, MAX (Nnames, 1)); + *Ntnames = Nnames; + return (names); + +} + +/* parse the GSC line, placing values in 'region'. + force RA[0] < RA[1], DEC[0] < DEC[1] */ +int parse_GSC_line (CatalogRegion *region, char *line) { + + double v0, s, dtmp; + char tmp[10]; + + /* copy name (4 chars) */ + bzero (region[0].name, 10); + strncpy (region[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; + region[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; + region[0].RA[1] = v0 * 15.0; + + /* this only happens at the 0,360 boundary */ + if (region[0].RA[1] < region[0].RA[0]) region[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; + region[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; + region[0].DEC[1] = s * v0; + + if (region[0].DEC[0] > region[0].DEC[1]) { + dtmp = region[0].DEC[1]; + region[0].DEC[1] = region[0].DEC[0]; + region[0].DEC[0] = dtmp; + } + + return (TRUE); + +} + Index: /branches/elixir/Ohana/src/gastro/src/gscdump.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gscdump.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gscdump.c (revision 2434) @@ -0,0 +1,36 @@ +# include "gastro.h" + +main (argc, argv) +int argc; +char **argv; +{ + + int i, N1, N2, Ncat, Nregions, NX, NY, Nmatch, status; + char **regions, **gregions(); + Stars *catalog, *stars1, *stars2, *gstars(); + struct timeval now, then; + Coords coords, dcoords; + double dNdM, dR, TotalArea; + float val; + + ConfigInit (); + if (argc < 4) { + fprintf (stderr, "USAGE gscdump RA DEC Radius\n"); + exit (0); + } + + coords.crval1 = atof (argv[1]); + coords.crval2 = atof (argv[2]); + NX = NY = 0; + + regions = gregions (&coords, &Nregions, &TotalArea); + + for (Ncat = i = 0; i < Nregions; i++) { + status = gcatalog (regions[i], &catalog, &Ncat, 0, 0); + } + + for (i = 0; i < Ncat; i++) { + fprintf (stdout, "%f %f %f\n", catalog[i].X, catalog[i].Y, catalog[i].mag); + } +} + Index: /branches/elixir/Ohana/src/gastro/src/gstars.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gstars.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gstars.c (revision 2434) @@ -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 (!fits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (1)\n", file); + exit(0); + } + /* get complete info from header */ + fits_scan (&header, "NAXIS1", "%d", 1, NX); + fits_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 { + fits_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 (!fits_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 (!fits_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) { + fits_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; + fits_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; + fits_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: /branches/elixir/Ohana/src/gastro/src/gstars1.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/gstars1.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/gstars1.c (revision 2434) @@ -0,0 +1,260 @@ +# 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 31 +# define BLOCK 1000 +# include +# include + +SStars *gstars (file, NSTARS, coords, NX, NY, dNdM) +char file[]; +int *NSTARS, *NX, *NY; +double *dNdM; +Coords *coords; +{ + + char line[64], *p, side[64]; + Header header, theader; + FILE *f; + int i, j, sign, Ninstar, nstar, rnumber, N, Nstars, nbytes, Nbytes, status, Nccd; + SStars *stars; + char *buffer; + double X, Y, T1, T2, T3, type, csign; + double PD, PR, DE, RE; + double ra, dec; + + /* read in image header, open image data region */ + if (!fits_read_header (file, &header)) { + fprintf (stderr, "ERROR: can't find image file %s (1)\n", file); + exit(0); + } + /* get complete info from header */ + fits_scan (&header, "NAXIS1", "%d", 1, NX); + fits_scan (&header, "NAXIS2", "%d", 1, NY); + + if (!strcasecmp (ROUGH_ASTROMETRY, "header")) { + if (!HEADER[0]) { + GetCoords (coords, &header); + } else { + fits_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 partial info from header */ + 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); + /* RA */ + fits_scan (&header, "RA", "%s", 1, line); + coords[0].crval1 = atof(&line[0]); + p = strchr (line, ':'); + coords[0].crval1 += atof((char *) (p+1)) / 60.0; + p = strchr (p+1, ':'); + coords[0].crval1 += atof((char *) (p+1)) / 3600.0; + coords[0].crval1 = coords[0].crval1 * 15.0; + /* DEC */ + fits_scan (&header, "DEC", "%s", 1, line); + coords[0].crval2 = atof(&line[0]); + sign = 1; + if (line[0] == '-') { + sign = -1; + coords[0].crval2 = fabs(coords[0].crval2 ); + } + p = strchr (line, ':'); + coords[0].crval2 += atof((char *) (p+1)) / 60.0; + p = strchr (p+1, ':'); + coords[0].crval2 += atof((char *) (p+1)) / 3600.0; + coords[0].crval2 *= sign; + } + if (VERBOSE) fprintf (stderr, "coordinates from header: %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) { + fits_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); + } + + 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); + } + + Nstars = 0; + fits_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 (&type, 5, &buffer[j*BYTES_STAR]); + if ((type == 4) || (type == 6) || (type == 8)) continue; + N++; + } + } + free (header.buffer); + free (buffer); + + 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); +} + +sort_stars (stars, N) +SStars *stars; +int N; +{ + int l,j,ir,i; + int temp; + 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; + fits_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: /branches/elixir/Ohana/src/gastro/src/line_fit.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/line_fit.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/line_fit.c (revision 2434) @@ -0,0 +1,112 @@ +# include "astro.h" + +int line_fit (stars1, stars2, N1, N2, dx, dy, Rx, Ry, + X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY) +Stars stars1[], stars2[]; +int N1, N2; +double dx, dy, Rx, Ry; +double *X_O, *X_X, *X_Y, *dX; +double *Y_O, *Y_X, *Y_Y, *dY; +{ + + int i, j, iteration; + double Xmin, Ymin, Xmax, Ymax; + double x, y, x2, y2, xy, N; + double X, Y, X2, Y2, Xy, Yx; + double Sx2, Sy2, Sxy, SXx, SXy, SYx, SYy; + double Dx, Dy, DD, d2X, d2Y, factor, dX_Y, dY_X; + FILE *f; + + factor = 1.0; + for (iteration = 0; iteration < 10; iteration ++) { + x = y = x2 = y2 = xy = X = Y = X2 = Xy = Yx = Y2 = N = 0; + for (i = 0; i < N1; i++) { + for (j = 0; j < N2; j++) { + if (iteration == 0) { + /* first iteration, more lax, use 2*RADIUS */ + Dx = 0.5*(stars1[i].X - Rx*stars2[j].X - dx) / Rx; + Dy = 0.5*(stars1[i].Y - Ry*stars2[j].Y - dy) / Ry; + } + else { + Dx = (stars1[i].X - (*X_O) - (*X_X)*stars2[j].X - (*X_Y)*stars2[j].Y) / Rx; + Dy = (stars1[i].Y - (*Y_O) - (*Y_X)*stars2[j].X - (*Y_Y)*stars2[j].Y) / Ry; + } + DD = factor * hypot (Dx, Dy); + if (DD < RADIUS) { + x += stars2[j].X; + y += stars2[j].Y; + x2 += stars2[j].X*stars2[j].X; + y2 += stars2[j].Y*stars2[j].Y; + xy += stars2[j].X*stars2[j].Y; + X += stars1[i].X; + Y += stars1[i].Y; + X2 += stars1[i].X*stars2[j].X; + Xy += stars1[i].X*stars2[j].Y; + Yx += stars1[i].Y*stars2[j].X; + Y2 += stars1[i].Y*stars2[j].Y; + N += 1.0; + } + } + } + + fprintf (stderr, "%d matched stars, iteration %d\n", (int) N, iteration + 1); + if (N < 3) { + fprintf (stderr, "too few stars, skipping\n"); + *X_O = *X_X = *X_Y = *Y_O = *Y_X = *Y_Y = 0; + return (N); + } + Sx2 = x2 - x*x/N; + Sy2 = y2 - y*y/N; + Sxy = xy - x*y/N; + SXx = X2 - X*x/N; + SXy = Xy - X*y/N; + SYx = Yx - Y*x/N; + SYy = Y2 - Y*y/N; + + *X_X = (SXx*Sy2 - SXy*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + *X_Y = (SXy*Sx2 - SXx*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + *X_O = X/N - (*X_X)*x/N - (*X_Y)*y/N; + + *Y_X = (SYx*Sy2 - SYy*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + *Y_Y = (SYy*Sx2 - SYx*Sxy) / (Sx2*Sy2 - Sxy*Sxy); + *Y_O = Y/N - (*Y_X)*x/N - (*Y_Y)*y/N; + /* make small course correction, keeping CCD rectilinear: */ +/* dX_Y = (*X_Y + *Y_X*(*Y_Y / *X_X)); */ +/* dY_X = 0.5*(*Y_X + *X_Y*(*X_Y / *Y_Y)); */ +/* *X_Y -= dX_Y; */ +/* *Y_X -= dY_X; */ + factor = factor * 1.1; + } + factor = factor / 1.1; + + d2X = d2Y = N = 0; + Ymin = Xmin = 10000; + Xmax = Ymax = 0; + for (i = 0; i < N1; i++) { + for (j = 0; j < N2; j++) { + Dx = (stars1[i].X - (*X_O) - (*X_X)*stars2[j].X - (*X_Y)*stars2[j].Y) / Rx; + Dy = (stars1[i].Y - (*Y_O) - (*Y_X)*stars2[j].X - (*Y_Y)*stars2[j].Y) / Ry; + DD = factor * hypot (Dx, Dy); + if (DD < RADIUS) { + d2X += Dx*Dx; + d2Y += Dy*Dy; + N += 1.0; + Xmin = MIN (Xmin, stars2[j].X); + Xmax = MAX (Xmax, stars2[j].X); + Ymin = MIN (Ymin, stars2[j].Y); + Ymax = MAX (Ymax, stars2[j].Y); + } + } + } + fprintf (stderr, "NPOINTS: %f\n", N); + fprintf (stderr, "Xrange: %f -- %f\n", Xmin, Xmax); + fprintf (stderr, "Yrange: %f -- %f\n", Ymin, Ymax); + (*dX) = sqrt(d2X / (N - 3.0)); + (*dY) = sqrt(d2Y / (N - 3.0)); + + fprintf (stderr, "%f %f %f %f\n", *X_O, *X_X, *X_Y, *dX); + fprintf (stderr, "%f %f %f %f\n", *Y_O, *Y_X, *Y_Y, *dY); + + return ((int)N); + +} Index: /branches/elixir/Ohana/src/gastro/src/misc.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/misc.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/misc.c (revision 2434) @@ -0,0 +1,129 @@ +# 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); + } + } +} + +/***** 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; +} + Index: /branches/elixir/Ohana/src/gastro/src/nrutil.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/nrutil.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/nrutil.c (revision 2434) @@ -0,0 +1,208 @@ +#include +#include + +void nrerror(error_text) +char error_text[]; +{ + void exit(); + + fprintf(stderr,"Numerical Recipes run-time error...\n"); + fprintf(stderr,"%s\n",error_text); + fprintf(stderr,"...now exiting to system...\n"); + exit(1); +} + + + +float *vector(nl,nh) +int nl,nh; +{ + float *v; + + v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float)); + if (!v) nrerror("allocation failure in vector()"); + return v-nl; +} + +int *ivector(nl,nh) +int nl,nh; +{ + int *v; + + v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int)); + if (!v) nrerror("allocation failure in ivector()"); + return v-nl; +} + +double *dvector(nl,nh) +int nl,nh; +{ + double *v; + + v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double)); + if (!v) nrerror("allocation failure in dvector()"); + return v-nl; +} + + + +float **matrix(nrl,nrh,ncl,nch) +int nrl,nrh,ncl,nch; +{ + int i; + float **m; + + m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*)); + if (!m) nrerror("allocation failure 1 in matrix()"); + m -= nrl; + + for(i=nrl;i<=nrh;i++) { + m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float)); + if (!m[i]) nrerror("allocation failure 2 in matrix()"); + m[i] -= ncl; + } + return m; +} + +double **dmatrix(nrl,nrh,ncl,nch) +int nrl,nrh,ncl,nch; +{ + int i; + double **m; + + m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*)); + if (!m) nrerror("allocation failure 1 in dmatrix()"); + m -= nrl; + + for(i=nrl;i<=nrh;i++) { + m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double)); + if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); + m[i] -= ncl; + } + return m; +} + +int **imatrix(nrl,nrh,ncl,nch) +int nrl,nrh,ncl,nch; +{ + int i,**m; + + m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*)); + if (!m) nrerror("allocation failure 1 in imatrix()"); + m -= nrl; + + for(i=nrl;i<=nrh;i++) { + m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int)); + if (!m[i]) nrerror("allocation failure 2 in imatrix()"); + m[i] -= ncl; + } + return m; +} + + + +float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl) +float **a; +int oldrl,oldrh,oldcl,oldch,newrl,newcl; +{ + int i,j; + float **m; + + m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*)); + if (!m) nrerror("allocation failure in submatrix()"); + m -= newrl; + + for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl; + + return m; +} + + + +void free_vector(v,nl,nh) +float *v; +int nl,nh; +{ + free((char*) (v+nl)); +} + +void free_ivector(v,nl,nh) +int *v,nl,nh; +{ + free((char*) (v+nl)); +} + +void free_dvector(v,nl,nh) +double *v; +int nl,nh; +{ + free((char*) (v+nl)); +} + + + +void free_matrix(m,nrl,nrh,ncl,nch) +float **m; +int nrl,nrh,ncl,nch; +{ + int i; + + for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); + free((char*) (m+nrl)); +} + +void free_dmatrix(m,nrl,nrh,ncl,nch) +double **m; +int nrl,nrh,ncl,nch; +{ + int i; + + for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); + free((char*) (m+nrl)); +} + +void free_imatrix(m,nrl,nrh,ncl,nch) +int **m; +int nrl,nrh,ncl,nch; +{ + int i; + + for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl)); + free((char*) (m+nrl)); +} + + + +void free_submatrix(b,nrl,nrh,ncl,nch) +float **b; +int nrl,nrh,ncl,nch; +{ + free((char*) (b+nrl)); +} + + + +float **convert_matrix(a,nrl,nrh,ncl,nch) +float *a; +int nrl,nrh,ncl,nch; +{ + int i,j,nrow,ncol; + float **m; + + nrow=nrh-nrl+1; + ncol=nch-ncl+1; + m = (float **) malloc((unsigned) (nrow)*sizeof(float*)); + if (!m) nrerror("allocation failure in convert_matrix()"); + m -= nrl; + for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl; + return m; +} + + + +void free_convert_matrix(b,nrl,nrh,ncl,nch) +float **b; +int nrl,nrh,ncl,nch; +{ + free((char*) (b+nrl)); +} Index: /branches/elixir/Ohana/src/gastro/src/plotstuff.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/plotstuff.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/plotstuff.c (revision 2434) @@ -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 ", 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 ", 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 ", 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 ", strlen (buffer)); + write (Xgraph[N], buffer2, 16); + write (Xgraph[N], buffer, strlen (buffer)); + +} + +void PlotVector (int Npts, float *vect, int mode, int N) { + + int Nbytes; + + if (Npts < 1) return; + + active = N; + Nbytes = Npts * sizeof (float); + write (Xgraph[N], vect, Nbytes); + +} + +void PlotReset (int N) { + + char buffer[128]; + int i; + + /* test Xgraph[N], flush junk from pipe */ + signal (SIGPIPE, XDead); + fcntl (Xgraph[N], F_SETFL, O_NONBLOCK); + for (i = 0; (read (Xgraph[N], buffer, 64) > 0) && (i < 20); i++); + fcntl (Xgraph[N], F_SETFL, !O_NONBLOCK); + + if (Xgraph[N] < 1) if (!open_graph(N)) return; + + write (Xgraph[N], "ERAS", 4); + +} +/* include these lines to plot a pair of vectors: + + typedef struct { + double xmin, xmax, ymin, ymax; + int style, ptype, ltype, etype, color; + double lweight, size; + } Graphdata; + Graphdata graphdata; + + graphdata.xmin = -200; + graphdata.xmax = 4200; + graphdata.ymin = -500; + graphdata.ymax = 500; + graphdata.style = 2; + graphdata.ptype = 2; + graphdata.ltype = 0; + graphdata.etype = 0; + graphdata.color = 0; + graphdata.lweight = 0; + graphdata.size = 0.5; + + PrepPlotting (N, &graphdata, n); + PlotVector (N, Y, 0, n); + PlotVector (N, dM, 1, n); + DonePlotting (&graphdata, n); + + */ Index: /branches/elixir/Ohana/src/gastro/src/precess.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/precess.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/precess.c (revision 2434) @@ -0,0 +1,31 @@ +# include "astro.h" + +void precess (Ra, Dec, in_epoch, out_epoch) +double *Ra, *Dec; +double in_epoch, out_epoch; +{ + + double T, zeta, theta, z, A, D; + double SD, CD, SA, CA; + + T = (out_epoch - in_epoch) / 100.0; + + zeta = RAD_DEG*(0.6406161*T + 0.0000839*T*T + 0.0000050*T*T*T); + theta = RAD_DEG*(0.5567530*T - 0.0001185*T*T - 0.0000116*T*T*T); + z = 0.6406161*T + 0.0003041*T*T + 0.0000051*T*T*T; + + A = *Ra; + D = *Dec; + + SD = cos(RAD_DEG*A + zeta)*sin(theta)*cos(RAD_DEG*D) + cos(theta)*sin(RAD_DEG*D); + CD = sqrt (1 - SD*SD); + SA = sin(RAD_DEG*A + zeta)*cos(RAD_DEG*D)/CD; + CA = (cos(RAD_DEG*A + zeta)*cos(theta)*cos(RAD_DEG*D) - sin(theta)*sin(RAD_DEG*D))/CD; + + *Dec = asin(SD) / RAD_DEG; + *Ra = atan2(SA, CA) / RAD_DEG + z; + + if (*Ra < 0) + *Ra += 360; +} + Index: /branches/elixir/Ohana/src/gastro/src/ralter_header.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/ralter_header.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/ralter_header.c (revision 2434) @@ -0,0 +1,47 @@ +# include "astro.h" + +void ralter_header (head, ref, X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY, N) +char head[], ref[]; +double X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY; +int N; +{ + + Header header; + char tmp[200], line[500], *c; + + fits_read_header (head, &header); + + fits_modify (&header, "X_O", "%lf", 1, X_O); + fits_modify (&header, "X_X", "%lf", 1, X_X); + fits_modify (&header, "X_Y", "%lf", 1, X_Y); + fits_modify (&header, "dX", "%lf", 1, dX); + fits_modify (&header, "Y_O", "%lf", 1, Y_O); + fits_modify (&header, "Y_X", "%lf", 1, Y_X); + fits_modify (&header, "Y_Y", "%lf", 1, Y_Y); + fits_modify (&header, "dY", "%lf", 1, dY); + if ((c = strrchr (ref, '/')) != NULL) + fits_modify (&header, "rREF", "%s", 1, c+1); + else + fits_modify (&header, "rREF", "%s", 1, ref); + fits_modify (&header, "rNSTARs", "%d", 1, N); + + fits_modify (&header, "X_O", "%C", 1, "rastro coeff"); + fits_modify (&header, "X_X", "%C", 1, "rastro coeff"); + fits_modify (&header, "X_Y", "%C", 1, "rastro coeff"); + fits_modify (&header, "dX", "%C", 1, "error (pixels)"); + fits_modify (&header, "Y_O", "%C", 1, "rastro coeff"); + fits_modify (&header, "Y_X", "%C", 1, "rastro coeff"); + fits_modify (&header, "Y_Y", "%C", 1, "rastro coeff"); + fits_modify (&header, "dY", "%C", 1, "error (pixels)"); + fits_modify (&header, "rREF", "%C", 1, "rastro ref file"); + fits_modify (&header, "rNSTARs", "%C", 1, "number of stars used"); + + strcpy (tmp, head); + strcat (tmp, "~"); + sprintf (line, "mv %s %s\0", head, tmp); + system (line); + + fits_write_header (head, &header); + fits_free_header (&header); + +} Index: /branches/elixir/Ohana/src/gastro/src/ranges.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/ranges.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/ranges.c (revision 2434) @@ -0,0 +1,38 @@ +# include "astro.h" + +void ranges (stars1, stars2, N1, N2, Xrange, Yrange, Xzero, Yzero, Rx, Ry) +Stars stars1[], stars2[]; +int N1, N2; +double *Xrange, *Yrange, *Xzero, *Yzero; +double Rx, Ry; +{ + + int i; + double maxX1, minX1, maxY1, minY1; + double maxX2, minX2, maxY2, minY2; + + 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 = Rx*stars2[0].X; + maxY2 = minY2 = Ry*stars2[0].Y; + for (i = 0; i < N2; i++) { + maxX2 = MAX (maxX2, (Rx*stars2[i].X)); + minX2 = MIN (minX2, (Rx*stars2[i].X)); + maxY2 = MAX (maxY2, (Ry*stars2[i].Y)); + minY2 = MIN (minY2, (Ry*stars2[i].Y)); + } + + *Xzero = minX1 - maxX2; + *Yzero = minY1 - maxY2; + *Xrange = ((maxX1 - minX2) - (minX1 - maxX2)); + *Yrange = ((maxY1 - minY2) - (minY1 - maxY2)); + +} Index: /branches/elixir/Ohana/src/gastro/src/rargs.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/rargs.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/rargs.c (revision 2434) @@ -0,0 +1,82 @@ +# include "astro.h" +# define NARGS 2 + +void rhelp () +{ + + fprintf (stderr, "rastro -- a relative astrometry routine\n"); + + fprintf (stderr, " USAGE: rastro \n"); + fprintf (stderr, " optional flags:\n"); + fprintf (stderr, " -N -- number of input stars used (default: 100) \n"); + fprintf (stderr, " -F n -- quadrant for search (default: all) \n"); + fprintf (stderr, " -R radius -- default radius of search (default: %d asec) \n", DEFAULT_RADIUS); + fprintf (stderr, " -G N -- grid search spacing (default: %d ) \n", DEFAULT_NPIX); + fprintf (stderr, "\n"); + fprintf (stderr, " -help -- print this list\n"); + exit (0); + +} + + + + +void rargs (argc, argv) +int *argc; +char **argv; +{ + + int N; + char line[500]; + + if (N = get_argument (*argc, argv, "-help")) { + remove_argument (N, argc, argv); + rhelp (); + } + + /** optional arguments **/ + if (N = get_argument (*argc, argv, "-N")) { + remove_argument (N, argc, argv); + N2NUMBER = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + N2NUMBER = 100; + } + + if (N = get_argument (*argc, argv, "-F")) { + remove_argument (N, argc, argv); + FLIP = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + FLIP = 0; + } + + if (N = get_argument (*argc, argv, "-R")) { + remove_argument (N, argc, argv); + RADIUS = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + RADIUS = DEFAULT_RADIUS; + } + + if (N = get_argument (*argc, argv, "-G")) { + remove_argument (N, argc, argv); + NPIX = atof (argv[N]); + remove_argument (N, argc, argv); + } + else { + NPIX = DEFAULT_NPIX; + } + + if (*argc != NARGS) { + fprintf (stderr, "USAGE: rastro [-N NN -F mode -R radius] \n"); + fprintf (stderr, " enter filenames, followed by EOF\n"); + exit (0); + } + + +} + Index: /branches/elixir/Ohana/src/gastro/src/rastro.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/rastro.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/rastro.c (revision 2434) @@ -0,0 +1,66 @@ +# include "astro.h" + +main (argc, argv) +int argc; +char **argv; +{ + + int N, N1, N2; + double dx, dy, Rx, Ry; + double X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY; + char *c, *ref, filename[200], head[200], obj[200]; + Stars *stars1, *stars2; + + rargs (&argc, argv); + + ref = argv[1]; + N1 = 5000; + stars1 = get_stars (ref, &N1); + + fprintf (stderr, "enter filenames, followed by EOF\n"); + while (fscanf (stdin, "%s", filename) != EOF) { + /**** See comment 1 below ****/ + strcpy (head, filename); + if ((c = strrchr(head, '.')) != NULL) + strcpy(strrchr(head, '.'), ".head"); + else + strcat(head, ".head"); + strcpy (obj, filename); + if ((c = strrchr(obj, '.')) != NULL) + strcpy(strrchr(obj, '.'), ".obj_out"); + else + strcat(obj, ".obj_out"); + + N2 = N2NUMBER; + stars2 = get_stars (obj, &N2); + if (N2 >= 3) { + Rx = 1.0; + Ry = 1.0; + + delta (stars1, stars2, N1, N2, &dx, &dy, &Rx, &Ry); + N = line_fit (stars1, stars2, N1, N2, dx, dy, Rx, Ry, + &X_O, &X_X, &X_Y, &dX, &Y_O, &Y_X, &Y_Y, &dY); + ralter_header (head, ref, X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY, N); + } + else { + fprintf (stderr, "too few stars in %s, skipped\n", filename); + } + free (stars2); + } + free (stars1); +} + + + + + + + /* comment 1: + assign the assumed names: + filename.head, filename.obj_out + if filename is of the form xxxx.xxx, cut off + extention before attaching .head, .obj_out + this entire problem can be removed by making the + Dophot output file contain the header info from + the image, but until this is the Dophot output + format... */ Index: /branches/elixir/Ohana/src/gastro/src/remove_argument.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/remove_argument.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/remove_argument.c (revision 2434) @@ -0,0 +1,21 @@ +# include "astro.h" + +int remove_argument (N, argc, argv) +int N; +int *argc; +char **argv; +{ + + int i; + + if ((N != (int)NULL) && (N != 0)) { + (*argc)--; + for (i = N; i < *argc; i++) { + argv[i] = argv[i+1]; + } + } + + return (N); + +} + Index: /branches/elixir/Ohana/src/gastro/src/rotate.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/rotate.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/rotate.c (revision 2434) @@ -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: /branches/elixir/Ohana/src/gastro/src/simulate.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/simulate.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/simulate.c (revision 2434) @@ -0,0 +1,102 @@ +# include "astro.h" +# define TIME {gettimeofday (&now, (struct timeval *) NULL); \ +fprintf (stderr, "start time: %f\n", now.tv_sec+now.tv_usec*1e-6);} +double rnd_gauss(); +# define T970101 852076800 +# define DRIFT_RATE 33.168 /* msec per row */ +# define MAG_LIMIT 20.100 + +main (argc, argv) +int argc; +char **argv; +{ + + int i, N1, N2, NX, NY, Ncat, Nregions, nstar; + char **regions, **gregions(); + Stars *catalog, *stars1, *stars2, *stars3; + struct timeval now; + Coords coords; + Header header; + FILE *f; + double X, Y, M, Rsigma, Msigma, Moffset; + + if (argc < 7) { + fprintf (stderr, "USAGE: simulate scale ra dec Moffset Msigma Rsigma template output\n"); + exit (0); + } + + Moffset = atof (argv[4]); + Msigma = atof (argv[5]); + Rsigma = atof (argv[6]); + + /* load in template file, determine size of "image" */ + f = fopen (argv[7], "r"); + if (f == NULL) { + fprintf (stderr, "could not find template file %s\n", argv[7]); + exit(0); + } + ALLOCATE (header.buffer, char, 8640); + fread (header.buffer, 1, 8640, f); + header.size = 8640; + fits_scan (&header, "NAXIS1", "%d", 1, &NX); + fits_scan (&header, "NAXIS2", "%d", 1, &NY); + gettimeofday (&now, (struct timeval *) NULL); + fprintf (stderr, "now: %d.%06d\n", now.tv_sec, now.tv_usec); + fits_modify (&header, "TZERO", "%d", 1, now.tv_sec - T970101); + fits_modify (&header, "TRATE", "%lf", 1, DRIFT_RATE); + fits_modify (&header, "FLIMIT", "%lf", 1, MAG_LIMIT); + fclose (f); + + /* set up coordinates for image */ + strcpy (coords.ctype, "RA---TAN"); + coords.crpix1 = 0.5*NX; + coords.crpix2 = 0.5*NY; + coords.pc1_1 = 1; coords.pc1_2 = 0; + coords.pc2_1 = 0; coords.pc2_2 = 1; + coords.cdelt1 = coords.cdelt2 = atof (argv[1]) / 3600.0; + coords.crval1 = atof (argv[2]); + coords.crval2 = atof (argv[3]); + fits_modify (&header, "RA", "%lf", 1, coords.crval1); + fits_modify (&header, "DEC", "%lf", 1, coords.crval2); + + /* find appropriate catalogs, load files */ + regions = gregions (coords, &Nregions); + for (N2 = i = 0; i < Nregions; i++) { + gcatalog (regions[i], &catalog, &Ncat); + } + gproject (catalog, &stars2, Ncat, &N2, &coords, NX, NY); + free (catalog); + + /* select catalog star within image region */ + ALLOCATE (stars3, Stars, N2); + for (nstar = i = 0; i < N2; i++) { + if ((stars2[i].X > 0) && (stars2[i].X < NX) && + (stars2[i].Y > 0) && (stars2[i].Y < NY)) { + stars3[nstar] = stars2[i]; + nstar++; + } + } + free (stars2); + fits_modify (&header, "NSTARS", "%d", 1, nstar); + + /* write header */ + f = fopen (argv[8], "w"); + if (f == NULL) { + fprintf (stderr, "could not open output file %s\n", argv[8]); + exit(0); + } + fwrite (header.buffer, 1, 8640, f); + + /* write data */ + gauss_init (); + for (i = 0; i < nstar; i++) { + X = rnd_gauss (stars3[i].X, Rsigma); + Y = rnd_gauss (stars3[i].Y, Rsigma); + M = rnd_gauss (stars3[i].mag, Msigma) + Moffset; + fprintf (f, "%6.1f %6.1f %6.3f %03.0f 1 3.2\n", + X, Y, M, Msigma*1000.0); + } + fclose (f); + +} + Index: /branches/elixir/Ohana/src/gastro/src/sort.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/sort.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/sort.c (revision 2434) @@ -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; + } +} Index: /branches/elixir/Ohana/src/gastro/src/stats.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/stats.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/stats.c (revision 2434) @@ -0,0 +1,82 @@ +# include "astro.h" + +void stats (head, RA, DEC, dRA, dDEC, ASEC_PIX) +char head[]; +double *RA, *DEC; +double *dRA, *dDEC; +double ASEC_PIX; +{ + + int j, status; + Header header; + char ra[100], dec[100]; + double epoch; + + status = fits_read_header (head, &header); + + if (!status) { + fprintf (stderr, "file %s not found\n", head); + exit (0); + } + + if (!FORCE) { + fits_scan (&header, "RA", "%s", 1, ra); + *RA = 15*dms_to_deg (ra, ':'); + + fits_scan (&header, "DEC", "%s", 1, dec); + *DEC = dms_to_deg (dec, ':'); + + /* FITS uses either EPOCH or EQUINOX. EPOCH is supposed to be + outmoded, but you can't tell,right? fits_scan will not alter the + value if it fails, so 2000.0 will act as a default, and EQUINOX + will override EPOCH if both exist */ + epoch = 2000.0; + fits_scan (&header, "EPOCH", "%lf", 1, &epoch); + fits_scan (&header, "EQUINOX", "%lf", 1, &epoch); + + fprintf (stderr, "%s -> %lf %s -> %lf, %lf\n", ra, *RA, dec, *DEC, epoch); + + precess (RA, DEC, epoch, 2000.0); /* I want only 2000.0 here! */ + + fprintf (stderr, "%s -> %lf %s -> %lf, %lf\n", ra, *RA, dec, *DEC, 2000.0); + } + + *dRA = header.Naxis[0] * ASEC_PIX / 3600.0; + *dDEC = header.Naxis[1] * ASEC_PIX / 3600.0; +} + +double dms_to_deg (string, sep) +char *string; +int sep; +{ + + char C[20]; + int c, d; + double h, flag; + + d = h = 0; + if (strchr (string, sep) == NULL) + return (HUGE_VAL); + c = strlen(string) - strlen(strchr(string, sep)); + strncpy (C, string, c); + C[c] = 0; + h = atof(C); + if (strchr(C, '-') != NULL) + flag = -1.0; + else + flag = 1.0; + h *= flag; + + d += c + 1; + c = strlen(&string[d]) - strlen(strchr(&string[d], sep)); + strncpy (C, &string[d], c); + C[c] = 0; + h += atof(C)/60.0; + + d += c + 1; + strcpy (C, &string[d]); + h += atof(C)/3600.; + + h *= flag; + return (h); +} Index: /branches/elixir/Ohana/src/gastro/src/svbksb.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/svbksb.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/svbksb.c (revision 2434) @@ -0,0 +1,24 @@ +void svbksb(u,w,v,m,n,b,x) +double **u,w[],**v,b[],x[]; +int m,n; +{ + int jj,j,i; + double s,*tmp,*vector(); + void free_vector(); + + tmp=vector(1,n); + for (j=1;j<=n;j++) { + s=0.0; + if (w[j]) { + for (i=1;i<=m;i++) s += u[i][j]*b[i]; + s /= w[j]; + } + tmp[j]=s; + } + for (j=1;j<=n;j++) { + s=0.0; + for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj]; + x[j]=s; + } + free_vector(tmp,1,n); +} Index: /branches/elixir/Ohana/src/gastro/src/svd.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/svd.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/svd.c (revision 2434) @@ -0,0 +1,90 @@ +# include + +/* + matrix[0..N-1][0..N-1] is a square matrix to be inverted, + vector[0..N-1][0..M-1] are several vectors to find solns' +*/ + +void svd (double ** matrix, int N, double **vector, int M) { + + int i, j, k; + double **u, *b, *w, **v, *x; + double wmax, wmin, term; + + ALLOCATE (u, double *, N + 1); + /* press routine wants u[1..N][1..N] */ + for (i = 0; i < N; i++) { + ALLOCATE (u[i+1], double, N + 1); + } + for (i = 0; i < N; i++) { + for (j = 0; j < N; j++) { + u[i+1][j+1] = matrix[i][j]; + } + } + ALLOCATE (v, double *, N + 1); + for (i = 0; i < N; i++) { + ALLOCATE (v[i+1], double, N + 1); + bzero (v[i+1], sizeof(double)*(N+1)); + } + ALLOCATE (w, double, N + 1); + bzero (w, sizeof(double)*(N+1)); + ALLOCATE (b, double, N + 1); + bzero (b, sizeof(double)*(N+1)); + ALLOCATE (x, double, N + 1); + bzero (x, sizeof(double)*(N+1)); + + fprintf (stderr, "sending matrix to svdcmp\n"); + svdcmp (u, N, N, w, v); + fprintf (stderr, "done with svdcmp\n"); + wmax = 0.0; + for (j = 1; j <= N; j++) if (w[j] > wmax) wmax = w[j]; + for (j = 1; j <= N; j++) fprintf (stderr, "%10.3e ", w[j]); + fprintf (stderr, "\n"); + wmin = wmax * 1e-22; + for (j = 1; j <= N; j++) if (w[j] < wmin) w[j] = 0.0; + for (j = 1; j <= N; j++) fprintf (stderr, "%10.3e ", w[j]); + fprintf (stderr, "\n"); + + for (i = 0; i < M; i++) { + fprintf (stderr, "preping backsub\n"); + for (j = 0; j < N; j++) { + b[j+1] = vector[j][i]; + } + + fprintf (stderr, "running backsub\n"); + for (j = 1; j <= N; j++) { + fprintf (stderr, "%10.3e ", b[j]); + } + fprintf (stderr, "\n"); + svbksb (u, w, v, N, N, b, x); + fprintf (stderr, "done with backsub\n"); + for (j = 0; j < N; j++) { + vector[j][i] = x[j+1]; + } + } + + /* replace matrix with inverse matrix */ + /* transpose u, invert w */ + for (i = 1; i <= N; i++) { + for (j = i+1; j <= N; j++) { + SWAP (u[i][j], u[j][i]); + } + if (w[i] != 0) { + w[i] = 1.0/w[i]; + } + } + + for (i = 1; i <= N; i++) { + for (j = 1; j <= N; j++) { + term = 0; + for (k = 1; k <= N; k++) { + term += matrix[i-1][j-1] * v[i][k]*w[k]*u[k][j]; + } + fprintf (stderr, "%6.4f ", term); + } + fprintf (stderr, "\n"); + } + + +} + Index: /branches/elixir/Ohana/src/gastro/src/svdcmp.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/svdcmp.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/svdcmp.c (revision 2434) @@ -0,0 +1,204 @@ +#include + +static double at,bt,ct; +#define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \ +(ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0)) + +static double maxarg1,maxarg2; +#define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\ + (maxarg1) : (maxarg2)) +#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) + +void svdcmp (double **a, int m, int n, double **w, double **v) { + + int flag,i,its,j,jj,k,l,nm; + double c,f,h,s,x,y,z; + double anorm=0.0,g=0.0,scale=0.0; + double *rv1,*dvector(); + void nrerror(),free_dvector(); + + if (m < n) nrerror("SVDCMP: You must augment A with extra zero rows"); + rv1=dvector(1,n); + for (i=1;i<=n;i++) { + l=i+1; + rv1[i]=scale*g; + g=s=scale=0.0; + if (i <= m) { + for (k=i;k<=m;k++) scale += fabs(a[k][i]); + if (scale) { + for (k=i;k<=m;k++) { + a[k][i] /= scale; + s += a[k][i]*a[k][i]; + } + f=a[i][i]; + g = -SIGN(sqrt(s),f); + h=f*g-s; + a[i][i]=f-g; + if (i != n) { + for (j=l;j<=n;j++) { + for (s=0.0,k=i;k<=m;k++) s += a[k][i]*a[k][j]; + f=s/h; + for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; + } + } + for (k=i;k<=m;k++) a[k][i] *= scale; + } + } + w[i]=scale*g; + g=s=scale=0.0; + if (i <= m && i != n) { + for (k=l;k<=n;k++) scale += fabs(a[i][k]); + if (scale) { + for (k=l;k<=n;k++) { + a[i][k] /= scale; + s += a[i][k]*a[i][k]; + } + f=a[i][l]; + g = -SIGN(sqrt(s),f); + h=f*g-s; + a[i][l]=f-g; + for (k=l;k<=n;k++) rv1[k]=a[i][k]/h; + if (i != m) { + for (j=l;j<=m;j++) { + for (s=0.0,k=l;k<=n;k++) s += a[j][k]*a[i][k]; + for (k=l;k<=n;k++) a[j][k] += s*rv1[k]; + } + } + for (k=l;k<=n;k++) a[i][k] *= scale; + } + } + anorm=MAX(anorm,(fabs(w[i])+fabs(rv1[i]))); + } + for (i=n;i>=1;i--) { + if (i < n) { + if (g) { + for (j=l;j<=n;j++) + v[j][i]=(a[i][j]/a[i][l])/g; + for (j=l;j<=n;j++) { + for (s=0.0,k=l;k<=n;k++) s += a[i][k]*v[k][j]; + for (k=l;k<=n;k++) v[k][j] += s*v[k][i]; + } + } + for (j=l;j<=n;j++) v[i][j]=v[j][i]=0.0; + } + v[i][i]=1.0; + g=rv1[i]; + l=i; + } + for (i=n;i>=1;i--) { + l=i+1; + g=w[i]; + if (i < n) + for (j=l;j<=n;j++) a[i][j]=0.0; + if (g) { + g=1.0/g; + if (i != n) { + for (j=l;j<=n;j++) { + for (s=0.0,k=l;k<=m;k++) s += a[k][i]*a[k][j]; + f=(s/a[i][i])*g; + for (k=i;k<=m;k++) a[k][j] += f*a[k][i]; + } + } + for (j=i;j<=m;j++) a[j][i] *= g; + } else { + for (j=i;j<=m;j++) a[j][i]=0.0; + } + ++a[i][i]; + } + for (k=n;k>=1;k--) { + for (its=1;its<=30;its++) { + flag=1; + for (l=k;l>=1;l--) { + nm=l-1; + if (fabs(rv1[l])+anorm == anorm) { + flag=0; + break; + } + if (fabs(w[nm])+anorm == anorm) break; + } + if (flag) { + c=0.0; + s=1.0; + for (i=l;i<=k;i++) { + f=s*rv1[i]; + if (fabs(f)+anorm != anorm) { + g=w[i]; + h=PYTHAG(f,g); + w[i]=h; + h=1.0/h; + c=g*h; + s=(-f*h); + for (j=1;j<=m;j++) { + y=a[j][nm]; + z=a[j][i]; + a[j][nm]=y*c+z*s; + a[j][i]=z*c-y*s; + } + } + } + } + z=w[k]; + if (l == k) { + if (z < 0.0) { + w[k] = -z; + for (j=1;j<=n;j++) v[j][k]=(-v[j][k]); + } + break; + } + if (its == 30) nrerror("No convergence in 30 SVDCMP iterations"); + x=w[l]; + nm=k-1; + y=w[nm]; + g=rv1[nm]; + h=rv1[k]; + f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); + g=PYTHAG(f,1.0); + f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; + c=s=1.0; + for (j=l;j<=nm;j++) { + i=j+1; + g=rv1[i]; + y=w[i]; + h=s*g; + g=c*g; + z=PYTHAG(f,h); + rv1[j]=z; + c=f/z; + s=h/z; + f=x*c+g*s; + g=g*c-x*s; + h=y*s; + y=y*c; + for (jj=1;jj<=n;jj++) { + x=v[jj][j]; + z=v[jj][i]; + v[jj][j]=x*c+z*s; + v[jj][i]=z*c-x*s; + } + z=PYTHAG(f,h); + w[j]=z; + if (z) { + z=1.0/z; + c=f*z; + s=h*z; + } + f=(c*g)+(s*y); + x=(c*y)-(s*g); + for (jj=1;jj<=m;jj++) { + y=a[jj][j]; + z=a[jj][i]; + a[jj][j]=y*c+z*s; + a[jj][i]=z*c-y*s; + } + } + rv1[l]=0.0; + rv1[k]=f; + w[k]=x; + } + } + free_dvector(rv1,1,n); +} + +#undef SIGN +#undef MAX +#undef PYTHAG Index: /branches/elixir/Ohana/src/gastro/src/test.c =================================================================== --- /branches/elixir/Ohana/src/gastro/src/test.c (revision 2434) +++ /branches/elixir/Ohana/src/gastro/src/test.c (revision 2434) @@ -0,0 +1,181 @@ +# include "../include/astro.h" + +main (argc, argv) +int argc; +char **argv; +{ + + Coords coords; + + coords.crpix1 = 1024; + coords.crpix2 = 512; + coords.crval1 = 10.0; + coords.crval2 = 89.0; + + coords.cdelt1 = 2.6; + coords.cdelt2 = 2.6; + + coords.pc1_1 = 1.0; coords.pc1_2 = 1.0; + coords.pc2_1 = 1.0; coords.pc2_2 = 1.0; + + + while (fscanf (stdin, "%lf%lf", &ra, &dec)!=EOF) { + RD_to_XY (&x, &y, ra, dec, &coords); + XY_to_RD (&Ra, &Dec, x, y, &coords); + fprintf (stderr, "%f %f --> %f %f --> %f %f\n", ra, dec, x, y, Ra, Dec); + } + +} + +XY_to_RD (ra, dec, x, y, coords) +double *ra, *dec; +double x, y; +Coords coords[]; +{ + + double L, M, X, Y, T; + double R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, cdel, sdp, cdp; + + *ra = 0; + *dec = 0; + + X = coords[0].cdelt1*(x - coords[0].crpix1); + Y = coords[0].cdelt2*(y - coords[0].crpix2); + L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2); + M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2); + + fprintf (stderr, "XY LM: %f %f %f %f\n", X, Y, L, M); + + if (!strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + R = hypot (L,M); + if ((L == 0) && (M == 0)) { + sphi = 0; + cphi = 1; + } + else { + sphi = L / R; + cphi = -M / R; + } + + if (!strcmp(&coords[0].ctype[4], "-TAN") ) { + if (R == 0) { + stht = 1.0; + ctht = 0.0; + } + else { + T = DEG_RAD / R; + stht = T / sqrt ( 1.0 + T*T); + ctht = 1.0 / sqrt ( 1.0 + T*T); + } + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + ctht = RAD_DEG * R; + stht = sqrt (1 - ctht*ctht); + } + + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (calp, salp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + coords[0].crval1; + *dec = DEG_RAD*delta; + } + if (!strcmp(&coords[0].ctype[4], "-LIN") || !strcmp(&coords[0].ctype[0], "GENE")) { + *ra = L + coords[0].crval1; + *dec = M + coords[0].crval2; + } + +} + + +RD_to_XY (x, y, ra, dec, coords) +double *x, *y; +double ra, dec; +Coords coords[]; +{ + + double tmp_d; + double X, Y, R, sphi, cphi, stht, ctht; + double alpha, delta, salp, calp, sdel, cdel, sdp, cdp; + + *x = 0; + *y = 0; + /* + if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) { + X = (ra - coords[0].crval1); + Y = (dec - coords[0].crval2); + } + */ + + if (!strcmp(&coords[0].ctype[4], "-TAN") || !strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + sdp = sin(RAD_DEG*coords[0].crval2); + cdp = cos(RAD_DEG*coords[0].crval2); + salp = sin(RAD_DEG*(ra - coords[0].crval1)); + calp = cos(RAD_DEG*(ra - coords[0].crval1)); + sdel = sin(RAD_DEG*dec); + cdel = cos(RAD_DEG*dec); + + stht = sdel*sdp + cdel*cdp*calp; + sphi = cdel*salp; + cphi = cdel*sdp*calp - sdel*cdp; + + if (!strcmp(&coords[0].ctype[4], "-TAN") ) { + X = DEG_RAD * sphi / stht; + Y = -DEG_RAD * cphi / stht; + } + if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) { + X = DEG_RAD * sphi; + Y = -DEG_RAD * cphi; + } + } + + tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1); + *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1; + *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2; + +} + + +/*** without sin/cos relationship simplifications ** + if (!strcmp(&ctype[4], "-TAN") || !strcmp(&ctype[4], "-SIN") || !strcmp(&ctype[0], "MM")) { + X = (x - CRPIX1)*CDELT1; + Y = (y - CRPIX2)*CDELT2; + + R = hypot (X,Y); + if ((X == 0) && (tmp_Y == 0)) + phi = 0; + else + phi = atan2 (X,-Y); *** phi in radians ** + } + sphi = sin(phi); + cphi = cos(phi); + + if (!strcmp(&ctype[4], "-TAN") ) { + if (R == 0) + tht = 0.0; + else + tht = atan (DEG_RAD / R); *** R in degrees, tht in radians ** + } + if (!strcmp(&ctype[4], "-SIN") || !strcmp(&ctype[0], "MM")) { + tht = acos (RAD_DEG * R); *** R in degrees, tht in radians ** + } + + sdp = sin(RAD_DEG*CRVAL2); + cdp = cos(RAD_DEG*CRVAL2); + + sdel = stht*sdp - ctht*cphi*cdp; + salp = ctht*sphi; + calp = stht*cdp + ctht*cphi*sdp; + alpha = atan2 (salp, calp); + delta = asin (sdel); + + *ra = DEG_RAD*alpha + CRVAL1; + *dec = DEG_RAD*delta; + } +*/