Index: /trunk/Ohana/src/gastro/Makefile =================================================================== --- /trunk/Ohana/src/gastro/Makefile (revision 3302) +++ /trunk/Ohana/src/gastro/Makefile (revision 3303) @@ -1,8 +1,9 @@ +default: gastro +help: + @echo "make options: gastro (default)" + include ../../Configure + HOME = $(ROOT)/src/gastro -CONFIG = $(ROOT)/config - -default: gastro - BIN = $(HOME)/bin INC = $(HOME)/include @@ -17,152 +18,72 @@ INCS = -I$(INC) -I$(LINC) -I$(XINC) LIBS = -L$(LLIB) -lFITS -lohana -lm -lsocket -lnsl -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) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(EOBJ) -o $(BIN)/extract $(CCFLAGS) - -install.extract: - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - 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) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(SOBJ) -o $(BIN)/simulate.$(ARCH) $(CCFLAGS) - -install.simulate: - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - rm -f $(DESTBIN)/simulate - cp $(BIN)/simulate.$(ARCH) $(DESTBIN)/simulate +CFLAGS = $(INCS) +LFLAGS = $(LIBS) # ----------------------- -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 \ +OBJ = \ +$(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)/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 +$(SRC)/greference.$(ARCH).o \ +$(SRC)/getusno.$(ARCH).o \ +$(SRC)/plotstuff.$(ARCH).o \ +$(SRC)/gptolemy.$(ARCH).o \ +$(SRC)/g2mass.$(ARCH).o \ +$(SRC)/misc.$(ARCH).o gastro: $(BIN)/gastro.$(ARCH) -$(BIN)/gastro.$(ARCH): $(GOBJ) +$(BIN)/gastro.$(ARCH): $(OBJ) + +$(OBJ): $(INC)/gastro.h + +INSTALL = gastro + +# generic dependancy rules for binary code ######################### +.PRECIOUS: %.$(ARCH).o +.PRECIOUS: $(BIN)/%.$(ARCH) + +%.$(ARCH).o : %.c + $(CC) $(CFLAGS) -c $< -o $@ + +$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(GOBJ) -o $(BIN)/gastro.$(ARCH) $(CCFLAGS) + $(CC) $^ -o $@ $(LFLAGS) -install.gastro: $(DESTBIN)/gastro +$(DESTBIN)/%: $(BIN)/%.$(ARCH) + @if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi + rm -f $(DESTBIN)/$* + cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$* -$(DESTBIN)/gastro: $(BIN)/gastro.$(ARCH) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - 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 +$(INSTALL) $(DEVEL): % : $(BIN)/%.$(ARCH) -$(DOBJ): $(INC)/gastro.h +%.clean : + rm -f $(BIN)/$*.$(ARCH) -gscdump: $(BIN)/gscdump.$(ARCH) +%.install: + make $(DESTBIN)/$* -$(BIN)/gscdump.$(ARCH): $(DOBJ) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(DOBJ) -o $(BIN)/gscdump.$(ARCH) $(CCFLAGS) +# utilities ################################################# -install.gscdump: $(DESTBIN)/gscdump +install: + for i in $(INSTALL); do make $$i.install; done -$(DESTBIN)/gscdump: $(BIN)/gscdump.$(ARCH) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - 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) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(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) - @if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi - $(CC) $(ROBJ) -o $(BIN)/rastro $(CCFLAGS) - -# ----------------------- -install: install.gastro - -clean: +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: unk/Ohana/src/gastro/include/astro.h =================================================================== --- /trunk/Ohana/src/gastro/include/astro.h (revision 3302) +++ (revision ) @@ -1,63 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/alter_header.c =================================================================== --- /trunk/Ohana/src/gastro/src/alter_header.c (revision 3302) +++ (revision ) @@ -1,60 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/args.c =================================================================== --- /trunk/Ohana/src/gastro/src/args.c (revision 3302) +++ (revision ) @@ -1,148 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/astro.c =================================================================== --- /trunk/Ohana/src/gastro/src/astro.c (revision 3302) +++ (revision ) @@ -1,85 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/config.c =================================================================== --- /trunk/Ohana/src/gastro/src/config.c (revision 3302) +++ (revision ) @@ -1,110 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/coordops.c =================================================================== --- /trunk/Ohana/src/gastro/src/coordops.c (revision 3302) +++ (revision ) @@ -1,362 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/dastro.c =================================================================== --- /trunk/Ohana/src/gastro/src/dastro.c (revision 3302) +++ (revision ) @@ -1,86 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/delta.c =================================================================== --- /trunk/Ohana/src/gastro/src/delta.c (revision 3302) +++ (revision ) @@ -1,53 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/extract.c =================================================================== --- /trunk/Ohana/src/gastro/src/extract.c (revision 3302) +++ (revision ) @@ -1,36 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/extract_catalog.c =================================================================== --- /trunk/Ohana/src/gastro/src/extract_catalog.c (revision 3302) +++ (revision ) @@ -1,179 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/find_shift.c =================================================================== --- /trunk/Ohana/src/gastro/src/find_shift.c (revision 3302) +++ (revision ) @@ -1,87 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gaussian.c =================================================================== --- /trunk/Ohana/src/gastro/src/gaussian.c (revision 3302) +++ (revision ) @@ -1,60 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gdelta.c =================================================================== --- /trunk/Ohana/src/gastro/src/gdelta.c (revision 3302) +++ (revision ) @@ -1,53 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/get_argument.c =================================================================== --- /trunk/Ohana/src/gastro/src/get_argument.c (revision 3302) +++ (revision ) @@ -1,19 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/get_catalog_stars.c =================================================================== --- /trunk/Ohana/src/gastro/src/get_catalog_stars.c (revision 3302) +++ (revision ) @@ -1,40 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/get_stars.c =================================================================== --- /trunk/Ohana/src/gastro/src/get_stars.c (revision 3302) +++ (revision ) @@ -1,91 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gfit.test.c =================================================================== --- /trunk/Ohana/src/gastro/src/gfit.test.c (revision 3302) +++ (revision ) @@ -1,496 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gregions.c =================================================================== --- /trunk/Ohana/src/gastro/src/gregions.c (revision 3302) +++ (revision ) @@ -1,218 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gscdump.c =================================================================== --- /trunk/Ohana/src/gastro/src/gscdump.c (revision 3302) +++ (revision ) @@ -1,36 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/gstars1.c =================================================================== --- /trunk/Ohana/src/gastro/src/gstars1.c (revision 3302) +++ (revision ) @@ -1,260 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/line_fit.c =================================================================== --- /trunk/Ohana/src/gastro/src/line_fit.c (revision 3302) +++ (revision ) @@ -1,112 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/nrutil.c =================================================================== --- /trunk/Ohana/src/gastro/src/nrutil.c (revision 3302) +++ (revision ) @@ -1,208 +1,0 @@ -#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: unk/Ohana/src/gastro/src/precess.c =================================================================== --- /trunk/Ohana/src/gastro/src/precess.c (revision 3302) +++ (revision ) @@ -1,31 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/ralter_header.c =================================================================== --- /trunk/Ohana/src/gastro/src/ralter_header.c (revision 3302) +++ (revision ) @@ -1,47 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/ranges.c =================================================================== --- /trunk/Ohana/src/gastro/src/ranges.c (revision 3302) +++ (revision ) @@ -1,38 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/rargs.c =================================================================== --- /trunk/Ohana/src/gastro/src/rargs.c (revision 3302) +++ (revision ) @@ -1,82 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/rastro.c =================================================================== --- /trunk/Ohana/src/gastro/src/rastro.c (revision 3302) +++ (revision ) @@ -1,66 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/remove_argument.c =================================================================== --- /trunk/Ohana/src/gastro/src/remove_argument.c (revision 3302) +++ (revision ) @@ -1,21 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/simulate.c =================================================================== --- /trunk/Ohana/src/gastro/src/simulate.c (revision 3302) +++ (revision ) @@ -1,102 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/stats.c =================================================================== --- /trunk/Ohana/src/gastro/src/stats.c (revision 3302) +++ (revision ) @@ -1,82 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/svbksb.c =================================================================== --- /trunk/Ohana/src/gastro/src/svbksb.c (revision 3302) +++ (revision ) @@ -1,24 +1,0 @@ -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: unk/Ohana/src/gastro/src/svd.c =================================================================== --- /trunk/Ohana/src/gastro/src/svd.c (revision 3302) +++ (revision ) @@ -1,90 +1,0 @@ -# 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: unk/Ohana/src/gastro/src/svdcmp.c =================================================================== --- /trunk/Ohana/src/gastro/src/svdcmp.c (revision 3302) +++ (revision ) @@ -1,204 +1,0 @@ -#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: unk/Ohana/src/gastro/src/test.c =================================================================== --- /trunk/Ohana/src/gastro/src/test.c (revision 3302) +++ (revision ) @@ -1,181 +1,0 @@ -# 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; - } -*/