IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3323


Ignore:
Timestamp:
Feb 24, 2005, 4:40:00 PM (21 years ago)
Author:
eugene
Message:

working version of mosastro

Location:
trunk/Ohana/src/mosastro
Files:
13 added
18 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/mosastro/Makefile

    r3294 r3323  
     1default: mosastro
     2help:
     3        @echo "make options: mosastro (default) mkstandards mkobs warptest"
     4
    15include ../../Configure
     6
    27HOME    =       $(ROOT)/src/mosastro.v3
    3 PROGRAM =       mosastro
    4 default: $(PROGRAM)
    5 
    68BIN     =       $(HOME)/bin
    79INC     =       $(HOME)/include
    … …  
    3234$(SRC)/sort.$(ARCH).o \
    3335$(SRC)/ConfigInit.$(ARCH).o \
     36$(SRC)/getstone.$(ARCH).o \
    3437$(SRC)/get2mass.$(ARCH).o \
    3538$(SRC)/getgsc.$(ARCH).o \
    … …  
    4245$(SRC)/field.$(ARCH).o \
    4346$(SRC)/chips.$(ARCH).o \
     47$(SRC)/GetGradients.$(ARCH).o \
     48$(SRC)/FitGradients.$(ARCH).o \
     49$(SRC)/FitChips.$(ARCH).o \
     50$(SRC)/fitpoly.$(ARCH).o \
     51$(SRC)/output.$(ARCH).o \
     52$(SRC)/wstars.$(ARCH).o \
     53$(SRC)/gaussj.$(ARCH).o \
     54$(SRC)/mkpolyterm.$(ARCH).o \
     55$(SRC)/mkmosaic.$(ARCH).o \
     56$(SRC)/GetScatter.$(ARCH).o \
     57$(SRC)/clip.$(ARCH).o \
    4458$(SRC)/parse_time.$(ARCH).o
    4559
    … …  
    5064$(SRC)/mkstandards.$(ARCH).o \
    5165$(SRC)/mkheader.$(ARCH).o \
     66$(SRC)/random.$(ARCH).o \
    5267$(SRC)/wstars.$(ARCH).o
    5368
    … …  
    5671
    5772OBS = \
     73$(SRC)/args_obs.$(ARCH).o \
    5874$(SRC)/mkobs.$(ARCH).o \
    5975$(SRC)/mkheader.$(ARCH).o \
    … …  
    6783$(SRC)/getgsc.$(ARCH).o \
    6884$(SRC)/getptolemy.$(ARCH).o \
     85$(SRC)/getstone.$(ARCH).o \
    6986$(SRC)/gptolemy.$(ARCH).o \
    7087$(SRC)/gregions.$(ARCH).o \
    7188$(SRC)/gcatalog.$(ARCH).o \
    7289$(SRC)/getusno.$(ARCH).o \
     90$(SRC)/mkpolyterm.$(ARCH).o \
     91$(SRC)/random.$(ARCH).o \
    7392$(SRC)/wstars.$(ARCH).o
    7493
    … …  
    7897warptest: $(BIN)/warptest.(ARCH)
    7998$(BIN)/warptest.(ARCH) : $(SRC)/warptest.$(ARCH).o
     99
     100$(MOS) : $(INC)/mosastro.h
     101$(OBS) : $(INC)/mosastro.h
     102$(STD) : $(INC)/mosastro.h
     103
     104INSTALL = mosastro mkstandards mkobs
    80105
    81106# dependancy rules for binary code #########################

  • trunk/Ohana/src/mosastro/include/mosastro.h

    r3294 r3323  
    88  double X, Y;  /* Chip Coords   - pixels  */
    99  double Mag, dMag;
    10   int type;
     10  int mask;
    1111} StarData;
    1212
    … …  
    2020  char *file;
    2121  int Nstars;
     22  int NX, NY;
    2223  StarData *stars;
     24  Header header;
     25  char *buffer;
     26  int Nbuffer;
    2327
    2428  int Nmatch;
    … …  
    6064} Stars;
    6165
     66int  ChipOrder;
    6267int  Nchip;
    6368Chip *chip;
    6469Field field;
    6570double Year;  /** carried for precession - probably put this in chip data **/
    66 double RADIUS; /** raw / ref matching radius (units?) **/
     71double RADIUS; /** raw / ref matching radius (pixels on Focal Plane) **/
     72double SIGMA;
     73double SIGMA_LIM;
    6774
     75char StoneRegions[256];
    6876char CDROM[256];
    6977char REFCAT[256];
    7078char CATDIR[256];
     79char ExptimeKeyword[256];
    7180char DateKeyword[256];
    7281char DateMode[256];
    … …  
    7685char TWO_MASS_DIR[256];
    7786int VERBOSE;
     87int NO_CHIPS;
    7888
    7989char GSCFILE[256];
    8090char GSC_DIR[256];
     91
     92char *FIELD;
     93char *CHIPS;
     94char *OUTPUT;
     95char *DUMP;
     96char *FOCAL_PLANE;
    8197
    8298StarData *getusno (CatStats *catstats, int *Nstars);
    … …  
    103119int fake_field (double RA, double DEC);
    104120Gradients *GetGradients ();
     121void FitChips (int Norder);
     122void FitChipLinear (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
     123void FitChipResid (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
     124double GetScatter ();
     125Header *mkmosaic (int Nx, int Ny, int Nstars, Coords *coords);

  • trunk/Ohana/src/mosastro/src/ConfigInit.c

    r3294 r3323  
    1717  if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
    1818
    19   GetConfig (config, "USNO_CDROM",        "%s", 0, CDROM);
    20   GetConfig (config, "DATE-KEYWORD",      "%s", 0, DateKeyword);
    21   GetConfig (config, "DATE-MODE",         "%s", 0, DateMode);
    22   GetConfig (config, "UT-KEYWORD",        "%s", 0, UTKeyword);
    23   GetConfig (config, "MJD-KEYWORD",       "%s", 0, MJDKeyword);
    24   GetConfig (config, "JD-KEYWORD",        "%s", 0, JDKeyword);
    25   ScanConfig (config, "2MASS_DIR",         "%s",  0, TWO_MASS_DIR);
     19  GetConfig (config, "USNO_CDROM",        "%s",  0, CDROM);
     20  GetConfig (config, "DATE-KEYWORD",      "%s",  0, DateKeyword);
     21  GetConfig (config, "DATE-MODE",         "%s",  0, DateMode);
     22  GetConfig (config, "UT-KEYWORD",        "%s",  0, UTKeyword);
     23  GetConfig (config, "MJD-KEYWORD",       "%s",  0, MJDKeyword);
     24  GetConfig (config, "JD-KEYWORD",        "%s",  0, JDKeyword);
     25  ScanConfig (config, "2MASS_DIR",        "%s",  0, TWO_MASS_DIR);
    2626
    2727  GetConfig (config, "ASTRO_REFCAT",      "%s",  0, REFCAT);
    28   ScanConfig (config, "CATDIR",            "%s",  0, CATDIR);
    29   ScanConfig (config, "GSCFILE",           "%s",  0, GSCFILE);
    30   ScanConfig (config, "GSCDIR",            "%s",  0, GSC_DIR);
     28  ScanConfig (config, "CATDIR",           "%s",  0, CATDIR);
     29  ScanConfig (config, "GSCFILE",          "%s",  0, GSCFILE);
     30  ScanConfig (config, "GSCDIR",           "%s",  0, GSC_DIR);
     31  ScanConfig (config, "STONE_DIR",        "%s",  0, StoneRegions);
     32
     33  /* abstracted header keywords */
     34  ScanConfig (config, "DATE-KEYWORD",     "%s",  0, DateKeyword);
     35  ScanConfig (config, "DATE-MODE",        "%s",  0, DateMode);
     36  ScanConfig (config, "UT-KEYWORD",       "%s",  0, UTKeyword);
     37  ScanConfig (config, "MJD-KEYWORD",      "%s",  0, MJDKeyword);
     38  ScanConfig (config, "JD-KEYWORD",       "%s",  0, JDKeyword);
     39  ScanConfig (config, "EXPTIME-KEYWORD",  "%s",  0, ExptimeKeyword);
    3140
    3241  GetConfig (config, "RADIUS",            "%lf", 0, &RADIUS);
     42  GetConfig (config, "SIGMA_LIM",         "%lf", 0, &SIGMA_LIM);
    3343
    3444  free (config);

  • trunk/Ohana/src/mosastro/src/LoadStars.c

    r3297 r3323  
    44int LoadStars (int Nfile, char **file) {
    55
    6   Header header;
     6  Header *header;
    77  int i, j, NCHIP, itmp, Nbytes, nbytes;
    88  FILE *f;
    9   char *buffer;
    109  time_t tsval;
    1110  struct tm *tmval;
    … …  
    2019
    2120    /* load header */
    22     if (!fits_read_header (chip[Nchip].file, &header)) {
     21    if (!fits_read_header (chip[Nchip].file, &chip[Nchip].header)) {
    2322      fprintf (stderr, "ERROR: can't read header for %s\n", file[i]);
    2423      free (chip[Nchip].file);
    … …  
    3130      continue;
    3231    }
    33     fseek (f, header.size, SEEK_SET);
     32    fseek (f, chip[Nchip].header.size, SEEK_SET);
     33
     34    /* get image dimensions */
     35    fits_scan (&chip[Nchip].header, "NAXIS1",   "%d", 1, &chip[Nchip].NX);
     36    fits_scan (&chip[Nchip].header, "NAXIS2",   "%d", 1, &chip[Nchip].NY);
    3437
    3538    /* get astrometry information */
    36     if (!GetCoords (&chip[Nchip].coords, &header)) {
     39    if (!GetCoords (&chip[Nchip].coords, &chip[Nchip].header)) {
    3740      fprintf (stderr, "ERROR: no astrometric solution in header\n");
    3841      free (chip[Nchip].file);
    … …  
    4346
    4447    itmp = 0;
    45     fits_scan (&header, "NASTRO",   "%d", 1, &itmp);
     48    fits_scan (&chip[Nchip].header, "NASTRO",   "%d", 1, &itmp);
    4649    if (itmp == 0) {
    4750      fprintf (stderr, "ERROR: bad astrometric solution in header %s\n", file[i]);
    … …  
    5154
    5255    /* get time info */
    53     tsval = parse_time (&header);
     56    tsval = parse_time (&chip[Nchip].header);
    5457    tmval = gmtime (&tsval);
    5558    Year = tmval[0].tm_year;
    5659
    5760    /* find expected number of stars */
    58     fits_scan (&header, "NSTARS", "%d", 1, &chip[Nchip].Nstars);
     61    fits_scan (&chip[Nchip].header, "NSTARS", "%d", 1, &chip[Nchip].Nstars);
    5962    if (chip[Nchip].Nstars == 0) {
    6063      fprintf (stderr, "ERROR: can't get NSTARS from header for %s\n", file[i]);
    … …  
    6568   
    6669    Nbytes = BYTES_STAR * chip[Nchip].Nstars;
    67     ALLOCATE (buffer, char, Nbytes + 1);
    68     nbytes = Fread (buffer, 1, Nbytes, f, "char");
     70    ALLOCATE (chip[Nchip].buffer, char, Nbytes + 1);
     71    nbytes = Fread (chip[Nchip].buffer, 1, Nbytes, f, "char");
    6972    if (nbytes != Nbytes) { exit (1); }
     73    chip[Nchip].Nbuffer = Nbytes;
    7074
    71     /** need to save all input lines so we can write on output? (gastro does not) **/
    7275    for (j = 0; j < chip[Nchip].Nstars; j++) {
    7376      bzero (&chip[Nchip].stars[j], sizeof(StarData));
    74       dparse (&chip[Nchip].stars[j].X,   1, &buffer[j*BYTES_STAR]);
    75       dparse (&chip[Nchip].stars[j].Y,   2, &buffer[j*BYTES_STAR]);
    76       dparse (&chip[Nchip].stars[j].Mag, 3, &buffer[j*BYTES_STAR]);
     77      dparse (&chip[Nchip].stars[j].X,    1, &chip[Nchip].buffer[j*BYTES_STAR]);
     78      dparse (&chip[Nchip].stars[j].Y,    2, &chip[Nchip].buffer[j*BYTES_STAR]);
     79      dparse (&chip[Nchip].stars[j].Mag,  3, &chip[Nchip].buffer[j*BYTES_STAR]);
     80      dparse (&chip[Nchip].stars[j].dMag, 4, &chip[Nchip].buffer[j*BYTES_STAR]);
     81      chip[Nchip].stars[j].dMag *= 0.001;  /* millimag errors stored in file */
    7782    }
    78     /** free this or keep with chip data? **/
    79     free (buffer);
    8083
    8184    if (VERBOSE) fprintf (stderr, "loaded %d stars from %s\n", chip[Nchip].Nstars, chip[Nchip].file);

  • trunk/Ohana/src/mosastro/src/args.c

    r3294 r3323  
    2525  }
    2626
    27   /*
    28   DUMP_REF = FALSE;
    29   if ((N = get_argument (*argc, argv, "-dumpref"))) {
    30     DUMP_REF = TRUE;
     27  DUMP = NULL;
     28  if ((N = get_argument (*argc, argv, "-dump"))) {
    3129    remove_argument (N, argc, argv);
    32   }
    33 
    34   DUMP_RAW = FALSE;
    35   if ((N = get_argument (*argc, argv, "-dumpraw"))) {
    36     DUMP_RAW = TRUE;
     30    DUMP = strcreate(argv[N]);
    3731    remove_argument (N, argc, argv);
    3832  }
    … …  
    5145    remove_argument (N, argc, argv);
    5246  }
    53   */
    5447
     48  /** currently unused **/
    5549  field.Norder = 0;
    5650  if ((N = get_argument (*argc, argv, "-order"))) {
    … …  
    5953    remove_argument (N, argc, argv);
    6054  }
     55
     56  /** currently unused **/
     57  ChipOrder = 1;
     58  if ((N = get_argument (*argc, argv, "-chiporder"))) {
     59    remove_argument (N, argc, argv);
     60    ChipOrder = atoi (argv[N]);
     61    remove_argument (N, argc, argv);
     62  }
    6163}

  • trunk/Ohana/src/mosastro/src/chips.c

    r3294 r3323  
    55
    66  int i, j;
    7   double R, D, P, Q, L, M;
     7  double R, D, P, Q, L, M, Scale, Det;
    88
    9 # if (0)
    109  if (CHIPS != (char *) NULL) {
    1110    load_chips (CHIPS);
    1211    return (1);
    1312  }
    14 # endif
    1513
    1614  for (i = 0; i < Nchip; i++) {
    1715
    1816    /* bore site center guess */
    19     strcpy (chip[i].map.ctype, "RA---WRP");
     17    strcpy (chip[i].map.ctype, "DEC--WRP");
    2018    chip[i].map.crpix1 = 0.0;
    2119    chip[i].map.crpix2 = 0.0;
    … …  
    3028    chip[i].map.crval2 = M;
    3129
    32     /** allow guess at field rotation?? add correct handling for negative cdelt1,2 **/
    33     chip[i].map.pc1_1  = chip[i].coords.pc1_1;
    34     chip[i].map.pc2_2  = chip[i].coords.pc2_2;
    35     chip[i].map.pc1_2  = chip[i].coords.pc1_2;
    36     chip[i].map.pc2_1  = chip[i].coords.pc2_1;
     30    /** we preserve the rotation and parity of coords.pc_ij, but renormalize to unity scale **/
     31    Det = chip[i].coords.pc1_1*chip[i].coords.pc2_2 - chip[i].coords.pc1_2*chip[i].coords.pc2_1;
     32    Scale = 1.0 / sqrt(fabs(chip[i].coords.cdelt1*chip[i].coords.cdelt2*Det));
     33
     34    /** test for NaN Scale **/
     35
     36    chip[i].map.pc1_1  = Scale * chip[i].coords.pc1_1 * chip[i].coords.cdelt1;
     37    chip[i].map.pc2_2  = Scale * chip[i].coords.pc2_2 * chip[i].coords.cdelt2;
     38    chip[i].map.pc1_2  = Scale * chip[i].coords.pc1_2 * chip[i].coords.cdelt2;
     39    chip[i].map.pc2_1  = Scale * chip[i].coords.pc2_1 * chip[i].coords.cdelt1;
    3740
    3841    chip[i].map.Npolyterms = 1;
    3942    for (j = 0; j < 7; j++) {
    40       chip[i].map.polyterms[0][j] = 0;
    41       chip[i].map.polyterms[1][j] = 0;
     43      chip[i].map.polyterms[j][0] = 0;
     44      chip[i].map.polyterms[j][1] = 0;
    4245    }
    4346  }
    4447  return (1);
    4548}
     49
     50int load_chips (char *filename) {
     51
     52  fprintf (stderr, "not ready yet\n");
     53  exit (1);
     54}

  • trunk/Ohana/src/mosastro/src/dump.c

    r3294 r3323  
    66
    77  for (i = 0; i < Nstars; i++) {
    8     fprintf (f, "%4d  %10.6f %10.6f  %8.2f %8.2f  %8.2f %8.2f  %7.2f %7.2f\n",
     8    fprintf (f, "%4d  %10.6f %10.6f  %8.2f %8.2f  %8.2f %8.2f  %7.2f %7.2f  %7.2f %7.2f\n",
    99             i,
    1010             stars[i].R, stars[i].D,
    … …  
    1616  return (1);
    1717}
     18
     19int dump_grads (Gradients *grad, char *filename) {
     20
     21  int i;
     22  FILE *f;
     23
     24  fprintf (stderr, "printing to file %s\n", filename);
     25
     26  f = fopen (filename, "w");
     27
     28  for (i = 0; i < grad[0].Npts; i++) {
     29    fprintf (f, "%4d  %10.6f %10.6f  %10.6f %10.6f   %10.6f %10.6f\n",
     30             i,
     31             grad[0].dPdL[i], grad[0].dPdM[i],
     32             grad[0].dQdL[i], grad[0].dQdM[i],
     33             grad[0].Lo[i], grad[0].Mo[i]);
     34  }
     35  fclose (f);
     36  return (1);
     37}
     38
     39int dump_match () {
     40
     41  int i;
     42  FILE *f, *g;
     43  f = fopen ("raw.dat", "w");
     44  g = fopen ("ref.dat", "w");
     45  for (i = 0; i < Nchip; i++) {
     46    dump_stars (f, chip[i].raw, chip[i].Nmatch);
     47    dump_stars (g, chip[i].ref, chip[i].Nmatch);
     48  }
     49  fclose (f);
     50  fclose (g);
     51  exit (1);
     52}
     53
     54int dump_rawstars () {
     55
     56  int i;
     57  FILE *f;
     58  f = fopen ("stars.dat", "w");
     59  for (i = 0; i < Nchip; i++) {
     60    dump_stars (f, chip[i].stars, chip[i].Nstars);
     61  }
     62  fclose (f);
     63  exit (1);
     64}
     65
     66int dump_refcat (StarData *refcat, int Nrefcat) {
     67
     68  int i;
     69  FILE *f;
     70  f = fopen ("refcat.dat", "w");
     71  dump_stars (f, refcat, Nrefcat);
     72  fclose (f);
     73  exit (1);
     74}

  • trunk/Ohana/src/mosastro/src/fakefield.c

    r3297 r3323  
    22
    33/* determine an initial guess to field parameters from data */
    4 int fake_field (double RA, double DEC) {
    5 
    6   int i;
     4int fake_field_center (double RA, double DEC) {
    75
    86  field.Rmin = RA  - 0.5;
    … …  
    1210
    1311  /* bore site center guess */
    14   strcpy (field.project.ctype, "RA---TAN");
    1512  field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
    1613  field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
     14 
     15  return (1);
     16}
     17
     18/* set default field parameters (overridden in args_obs) */
     19int fake_field_defaults () {
     20
     21  int i;
     22
     23  /* bore site center guess */
     24  strcpy (field.project.ctype, "DEC--TAN");
    1725  field.project.crpix1 = 0;
    1826  field.project.crpix2 = 0;
    1927 
    20   /* measure average plate scale - would be better using parabolic min... */
     28  /* set TP plate scale */
    2129  field.project.cdelt1 = 1.0/3600.0;
    2230  field.project.cdelt2 = 1.0/3600.0;
    … …  
    2937  field.project.Npolyterms = 1;
    3038
    31   /* bore site center guess */
    32   strcpy (field.distort.ctype, "RA---WRP");
     39  /** distort only has power in polyterms **/
     40  strcpy (field.distort.ctype, "DEC--WRP");
    3341  field.distort.crval1 = 0.0;
    3442  field.distort.crval2 = 0.0;
    … …  
    3846  field.distort.cdelt2 = 1.0;
    3947 
    40   /** allow guess at field rotation?? **/
    4148  field.distort.pc1_1  = 1;
    4249  field.distort.pc2_2  = 1;
    … …  
    4451  field.distort.pc2_1  = 0;
    4552
    46   /* allow 2nd and 3rd order? */
    4753  /* how do we handle renormalization? (fixed at 1000 pixels??) */
    4854  field.distort.Npolyterms = 3;
    4955  for (i = 0; i < 7; i++) {
    50     field.distort.polyterms[0][i] = 0;
    51     field.distort.polyterms[1][i] = 0;
     56    field.distort.polyterms[i][0] = 0;
     57    field.distort.polyterms[i][1] = 0;
    5258  }
    53 
    54   /* megacam distortion : 50 arcsec @ 1800 arcsec */
    55   field.distort.polyterms[0][3] = 1e-8;
    56   field.distort.polyterms[1][6] = 1e-8;
    57 
    58   /* tiny distortion    : 0.5 arcsec @ 1800 arcsec */
    59   field.distort.polyterms[0][3] = 1e-10;
    60   field.distort.polyterms[1][6] = 1e-10;
    61 
    62   /* minor distortion    : 5 arcsec @ 1800 arcsec */
    63   field.distort.polyterms[0][3] = 1e-9;
    64   field.distort.polyterms[1][6] = 1e-9;
    6559
    6660  return (1);

  • trunk/Ohana/src/mosastro/src/field.c

    r3294 r3323  
    77  double PS;
    88
    9 # if (0)
    109  if (FIELD != (char *) NULL) {
    1110    load_field (FIELD);
    12     return (1);
     11    return (TRUE);
    1312  }
    14 # endif
    1513
    1614  /* bore site center guess */
    17   strcpy (field.project.ctype, "RA---TAN");
     15  strcpy (field.project.ctype, "DEC--TAN");
    1816  field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
    1917  field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
    … …  
    4038
    4139  /* bore site center guess */
    42   strcpy (field.distort.ctype, "RA---WRP");
     40  strcpy (field.distort.ctype, "DEC--WRP");
    4341  field.distort.crval1 = 0.0;
    4442  field.distort.crval2 = 0.0;
    … …  
    5856  field.distort.Npolyterms = 1;
    5957  for (i = 0; i < 7; i++) {
    60     field.distort.polyterms[0][i] = 0;
    61     field.distort.polyterms[1][i] = 0;
     58    field.distort.polyterms[i][0] = 0;
     59    field.distort.polyterms[i][1] = 0;
    6260  }
    63   return (1);
     61  return (TRUE);
    6462}
    6563
    … …  
    8785
    8886}
     87
     88void field_combine () {
     89
     90  int i;
     91
     92  strcpy (field.project.ctype, "DEC--WRP");
     93
     94  field.project.Npolyterms = field.distort.Npolyterms;
     95  for (i = 0; i < 7; i++) {
     96    field.project.polyterms[i][0] = field.distort.polyterms[i][0];
     97    field.project.polyterms[i][1] = field.distort.polyterms[i][1];
     98  }
     99}
     100
     101int load_field (char *filename) {
     102
     103  int i;
     104  Coords coords;
     105  Header header;
     106
     107  /* load header */
     108  if (!fits_read_header (filename, &header)) {
     109    fprintf (stderr, "ERROR: can't read header for %s\n", filename);
     110    exit (1);
     111  }
     112  /* get astrometry information */
     113  if (!GetCoords (&coords, &header)) {
     114    fprintf (stderr, "ERROR: no astrometric solution in field %s\n", filename);
     115    exit (1);
     116  }
     117 
     118  /* bore site center guess */
     119  strcpy (field.project.ctype, "DEC--TAN");
     120  field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
     121  field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
     122  field.project.crpix1 = 0;
     123  field.project.crpix2 = 0;
     124 
     125  /* measure average plate scale - would be better using parabolic min... */
     126  field.project.cdelt1 = coords.cdelt1;
     127  field.project.cdelt2 = coords.cdelt2;
     128
     129  /** allow guess at field rotation?? **/
     130  field.project.pc1_1  = coords.pc1_1;
     131  field.project.pc2_2  = coords.pc2_2;
     132  field.project.pc1_2  = coords.pc1_2;
     133  field.project.pc2_1  = coords.pc2_1;
     134  field.project.Npolyterms = 1;
     135
     136  /* bore site center guess */
     137  strcpy (field.distort.ctype, "DEC--WRP");
     138  field.distort.crval1 = 0.0;
     139  field.distort.crval2 = 0.0;
     140  field.distort.crpix1 = 0.0;
     141  field.distort.crpix2 = 0.0;
     142  field.distort.cdelt1 = 1.0;
     143  field.distort.cdelt2 = 1.0;
     144 
     145  /** allow guess at field rotation?? **/
     146  field.distort.pc1_1  = 1;
     147  field.distort.pc2_2  = 1;
     148  field.distort.pc1_2  = 0;
     149  field.distort.pc2_1  = 0;
     150
     151  /* allow 2nd and 3rd order? */
     152  /* how do we handle renormalization? (fixed at 1000 pixels??) */
     153  field.distort.Npolyterms = coords.Npolyterms;
     154  for (i = 0; i < 7; i++) {
     155    field.distort.polyterms[i][0] = coords.polyterms[i][0];
     156    field.distort.polyterms[i][1] = coords.polyterms[i][1];
     157  }
     158  return (TRUE);
     159}
     160

  • trunk/Ohana/src/mosastro/src/getusno.c

    r3297 r3323  
    2424
    2525  RA0  = catstats[0].RA[0];
    26   DEC0 = catstats[0].RA[1];
    27   RA1 = catstats[0].DEC[0];
     26  RA1 = catstats[0].RA[1];
     27  DEC0 = catstats[0].DEC[0];
    2828  DEC1 = catstats[0].DEC[1];
    2929

  • trunk/Ohana/src/mosastro/src/greference.c

    r3294 r3323  
    1616
    1717  if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", catstats.RA[0], catstats.RA[1], catstats.DEC[0], catstats.DEC[1]);
     18
     19  /* get stars from the Stone et al catalog for the given region */
     20  if (!strcmp (REFCAT, "STONE")) {
     21    stars = getstone (&catstats, &Nstars);
     22  }
    1823
    1924  /* get stars from the USNO catalog for the given region */

  • trunk/Ohana/src/mosastro/src/match.c

    r3294 r3323  
    66  int i, j, k, N, K, Ntotal;
    77  int Nmatch, NMATCH;
    8   double dp, dq, radius;
     8  double dp, dq, radius, Radius;
    99  double *p, *q, *P, *Q;
    1010  int *U, *u;
     11
     12  /* requested radius is in arcsec ; internally radius is in pixels */
     13  Radius = RADIUS / (3600.0 * field.project.cdelt2);
    1114
    1215  /* sort the REFCAT data by P in tangent plane */
    … …  
    4245    /* find star matches in the tangent plane coord system */
    4346    for (j = k = 0; (j < chip[i].Nstars) && (k < Nrefcat); ) {
     47      /* skip anything with dMag too large */
     48      if (SIGMA_LIM > 0.0) {
     49        if (chip[i].stars[u[j]].dMag > SIGMA_LIM) {
     50          j++;
     51          continue;
     52        }
     53      }
     54       
    4455      dp = p[j] - P[k];
    4556       
    46       if (dp <= -2*RADIUS) {
     57      if (dp <= -2*Radius) {
    4758        j++;
    4859        continue;
    4960      }
    50       if (dp >= 2*RADIUS) {
     61      if (dp >= 2*Radius) {
    5162        k++;
    5263        continue;
    … …  
    5465
    5566      K = k;
    56       for (; (dp > -2*RADIUS) && (k < Nrefcat); k++) {
     67      for (; (dp > -2*Radius) && (k < Nrefcat); k++) {
    5768        dp = p[j] - P[k];
    5869        dq = q[j] - Q[k];
    5970        radius = hypot (dp, dq);
    60         if (radius < RADIUS) {
     71        if (radius < Radius) {
    6172          chip[i].raw[Nmatch] = chip[i].stars[u[j]];
    6273          chip[i].ref[Nmatch] = refcat[U[k]];

  • trunk/Ohana/src/mosastro/src/mkheader.c

    r3294 r3323  
    77  ALLOCATE (header, Header, 1);
    88
     9  header[0].simple = TRUE;
    910  header[0].bitpix = -32;
    1011  header[0].unsign = FALSE;
    … …  
    2223  fits_modify (header, "UTC-OBS",  "%s", 1, "14:27:45.30");
    2324  fits_modify (header, "ZERO_PT", "%lf", 1, 25.0);
     25  fits_modify (header, "EXPTIME", "%lf", 1, 2.0);
    2426
    25   if (coords[0].Npolyterms > 1) {
    26     fits_modify (header, "CTYPE1",   "%s",  1, "RA---PLY");
    27     fits_modify (header, "CTYPE2",   "%s",  1, "DEC--PLY");
    28   } else {
    29     fits_modify (header, "CTYPE1",   "%s",  1, "RA---TAN");
    30     fits_modify (header, "CTYPE2",   "%s",  1, "DEC--TAN");
    31   }   
     27  PutCoords (coords, header);
     28
    3229  fits_modify (header, "NASTRO",   "%d", 1, 1);
    33 
    34   fits_modify (header, "CDELT1",   "%le", 1, coords[0].cdelt1);
    35   fits_modify (header, "CDELT2",   "%le", 1, coords[0].cdelt2);
    36   fits_modify (header, "CRVAL1",   "%lf", 1, coords[0].crval1);
    37   fits_modify (header, "CRVAL2",   "%lf", 1, coords[0].crval2); 
    38   fits_modify (header, "CRPIX1",   "%lf", 1, coords[0].crpix1);
    39   fits_modify (header, "CRPIX2",   "%lf", 1, coords[0].crpix2);
    40   fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
    41   fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
    42   fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
    43   fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
    44   fits_modify (header, "NPLYTERM", "%d", 1, coords[0].Npolyterms);
    45   if (coords[0].Npolyterms > 1) {
    46     fits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]);
    47     fits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]);
    48     fits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]);
    49 
    50     fits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]);
    51     fits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]);
    52     fits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]);
    53   }
    54      
    55   if (coords[0].Npolyterms > 2) {
    56     fits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]);
    57     fits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]);
    58     fits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]);
    59     fits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]);
    60 
    61     fits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]);
    62     fits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]);
    63     fits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]);
    64     fits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]);
    65   }
    6630
    6731  return (header);

  • trunk/Ohana/src/mosastro/src/mkobs.c

    r3297 r3323  
    44
    55  int i, j, N, Nrefcat, Nchip, Nstars, NSTARS;
     6  double Ro, Do, Po, Qo, Lo, Mo;
     7  double dX, dY, Dist;
    68  double RA, DEC;
    7   double Ro, Do, Po, Qo, Lo, Mo;
    8   double dX, dY;
    99  char filename[64];
    1010  Coords map, coords;
    … …  
    1414
    1515  ConfigInit (&argc, argv);
    16 
    17   if (argc != 4) {
    18     fprintf (stderr, "USAGE: mkobs (RA) (DEC) (output)\n");
    19     exit (2);
    20   }
     16  args (&argc, argv);
     17  init_random ();
    2118
    2219  RA  = atof(argv[1]);
    2320  DEC = atof(argv[2]);
    24   dX = dY = 3600*0.1;
     21  OUTPUT = argv[3];
     22  fake_field_center (RA, DEC);
    2523
    26   fake_field (RA, DEC);
     24  /* chip size in pixels (output is 5x5 grid of chips */
     25  dX = 0.1 / field.project.cdelt1;
     26  dY = 0.1 / field.project.cdelt2;
    2727
    2828  refcat = greference (&Nrefcat);
    2929  project_refcat (refcat, Nrefcat);
     30  if (FOCAL_PLANE != NULL) {
     31    FILE *f;
     32    f = fopen (FOCAL_PLANE, "w");
     33    dump_stars (f, refcat, Nrefcat);
     34    fclose (f);
     35  }
     36  if (NO_CHIPS) exit (0);
    3037
    3138  Nchip = 0;
    … …  
    3441
    3542      /* find chip center in TP, FP and Sky */
    36       Po = 0.2*i*3600.0;
    37       Qo = 0.2*j*3600.0;
     43      Po = 2*dX*i;
     44      Qo = 2*dY*j;
    3845      RD_to_XY (&Lo, &Mo, Po, Qo, &field.distort);
    3946      XY_to_RD (&Ro, &Do, Po, Qo, &field.project);
    4047
    4148      /* FP-Chip terms */
    42       strcpy (map.ctype, "RA---WRP");
     49      strcpy (map.ctype, "DEC--WRP");
    4350      map.crval1 = Lo;
    4451      map.crval2 = Mo;
    … …  
    5360      map.pc1_2  = 0;
    5461      map.pc2_1  = 0;
    55       map.Npolyterms = 1;
     62
     63      map.Npolyterms = 2;
     64      for (N = 0; N < 7; N++) {
     65        map.polyterms[N][0] = 0;
     66        map.polyterms[N][1] = 0;
     67      }
     68      map.polyterms[2][1] = 1e-5;
    5669
    5770      /* project catalog stars to chip */
    … …  
    6982        if (refcat[N].Y >= 2*dY) continue;
    7083
    71         stars[Nstars].X = refcat[N].X;
    72         stars[Nstars].Y = refcat[N].Y;
     84        /* add random noise - not gaussian noise */
     85        stars[Nstars].X = refcat[N].X + SIGMA*(drand48() - 0.5);
     86        stars[Nstars].Y = refcat[N].Y + SIGMA*(drand48() - 0.5);
    7387        stars[Nstars].M = 16.0;
    7488        stars[Nstars].dM = 0.02;
    … …  
    88102      }
    89103
     104
    90105      /* Chip-Sky terms */
    91       strcpy (coords.ctype, "RA---WRP");
     106      strcpy (coords.ctype, "DEC--WRP");
    92107      coords.crval1 = Ro;
    93108      coords.crval2 = Do;
    94109      coords.crpix1 = dX;
    95       coords.crpix2 = dX;
     110      coords.crpix2 = dY;
    96111 
    97       coords.cdelt1 = 1/3600.0;
    98       coords.cdelt2 = 1/3600.0;
     112      {
     113        double dP, dQ, dL, dM;
     114        double Mx, Lx, scale;
     115
     116        dP = 10;
     117        dQ = 10;
     118        RD_to_XY (&Lx, &Mx, Po+dP, Qo+dQ, &field.distort);
     119        dL = Lx - Lo;
     120        dM = Mx - Mo;
     121        scale = hypot (10.0, 10.0) / hypot (dL, dM);
     122        coords.cdelt1 = scale/3600.0;
     123        coords.cdelt2 = scale/3600.0;
     124      }
    99125
    100126      coords.pc1_1  = 1;
    … …  
    105131
    106132      header = mkheader (2*dX, 2*dY, Nstars, &coords);
    107       fits_modify (header, "C_L0", "%lf", 1, Lo);
    108       fits_modify (header, "C_M0", "%lf", 1, Mo);
    109133
    110       sprintf (filename, "%s.%02d.fits", argv[3], Nchip);
     134      sprintf (filename, "%s.%02d.fits", OUTPUT, Nchip);
    111135      wstars (filename, stars, Nstars, header);
    112136      fits_free_header (header);

  • trunk/Ohana/src/mosastro/src/mkstandards.c

    r3297 r3323  
    11# include "mosastro.h"
    2 # define DX 100.0
    3 # define DY 100.0
     2# define DX 30.0
     3# define DY 30.0
     4
     5extern double drand48();
    46
    57/* build a grid of reference stars */
    68int main (int argc, char **argv) {
    79
    8   int Nstars, NSTARS;
    9   double x, y, dX, dY;
     10  int i, N, Nstars, NSTARS, Random;
     11  double x, y, m, dX, dY;
    1012  Header *header;
    1113  Coords coords;
    1214  Stars *stars;
    1315
     16  Random = FALSE;
     17  if ((N = get_argument (argc, argv, "-random"))) {
     18    Random = TRUE;
     19    remove_argument (N, &argc, argv);
     20    Nstars = atoi (argv[N]);
     21    remove_argument (N, &argc, argv);
     22  }
     23
    1424  if (argc != 5) {
    15     fprintf (stderr, "USAGE: mkstandards (RA) (DEC) (half-width) (output)\n");
     25    fprintf (stderr, "USAGE: mkstandards (RA) (DEC) (half-width) (output) [-random Nstars]\n");
    1626    exit (2);
    1727  }
    … …  
    2131
    2232  /* bore site center guess */
    23   strcpy (coords.ctype, "RA---TAN");
     33  strcpy (coords.ctype, "DEC--TAN");
    2434  coords.crval1 = atof(argv[1]);
    2535  coords.crval2 = atof(argv[2]);
    … …  
    3646  coords.Npolyterms = 1;
    3747
    38   Nstars = 0;
    39   NSTARS = 1000;
    40   ALLOCATE (stars, Stars, NSTARS);
     48  if (Random) {
     49   
     50    init_random ();
    4151
    42   for (x = 0; x < 2*dX; x += DX) {
    43     for (y = 0; y < 2*dY; y += DY) {
    44       stars[Nstars].X = x;
    45       stars[Nstars].Y = y;
    46       stars[Nstars].M = 16.0;
    47       stars[Nstars].dM = 0.02;
    48       stars[Nstars].dophot = 1;
    49       stars[Nstars].sky = 1.0;
    50       stars[Nstars].Mgal = 16.0;
    51       stars[Nstars].Map = 16.0;
    52       stars[Nstars].fx = 1.0;
    53       stars[Nstars].fy = 1.0;
    54       stars[Nstars].df = 0.0;
     52    ALLOCATE (stars, Stars, Nstars);
     53    for (i = 0; i < Nstars; i++) {
     54      x = 2*dX*drand48();
     55      y = 2*dY*drand48();
     56      m = 14.0 + 4.0*drand48();
     57      stars[i].X = x;
     58      stars[i].Y = y;
     59      stars[i].M = m;
     60      stars[i].dM = 0.02;
     61    }
     62  } else {
    5563
    56       Nstars ++;
    57       if (Nstars >= NSTARS) {
    58         NSTARS += 1000;
    59         REALLOCATE (stars, Stars, NSTARS);
     64    Nstars = 0;
     65    NSTARS = 1000;
     66    ALLOCATE (stars, Stars, NSTARS);
     67
     68    for (x = 0; x < 2*dX; x += DX) {
     69      for (y = 0; y < 2*dY; y += DY) {
     70        stars[Nstars].X = x;
     71        stars[Nstars].Y = y;
     72        stars[Nstars].M = 16.0;
     73
     74        Nstars ++;
     75        if (Nstars >= NSTARS) {
     76          NSTARS += 1000;
     77          REALLOCATE (stars, Stars, NSTARS);
     78        }
    6079      }
    6180    }
     81  }
     82
     83  for (i = 0; i < Nstars; i++) {
     84    stars[i].dophot = 1;
     85    stars[i].sky = 1.0;
     86    stars[i].Mgal = 16.0;
     87    stars[i].Map = 16.0;
     88    stars[i].fx = 1.0;
     89    stars[i].fy = 1.0;
     90    stars[i].df = 0.0;
    6291  }
    6392

  • trunk/Ohana/src/mosastro/src/mosastro.c

    r3297 r3323  
    44int main (int argc, char **argv) {
    55
    6   int Nrefcat;
     6  int Nrefcat, Nastro;
     7  double Cerror;
    78  glob_t pglob;
    89  StarData *refcat;
    … …  
    1213  args (&argc, argv);
    1314
    14   if (argc != 3) {
    15     fprintf (stderr, "USAGE: mosastro (inglob) (ext)\n");
     15  if (argc != 4) {
     16    fprintf (stderr, "USAGE: mosastro (inglob) (ext) (mosaic.phu)\n");
    1617    exit (2);
    1718  }
    … …  
    3031  /* use field model to get TP & FP coords */
    3132  project_stars ();
     33  if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars();
    3234
    3335  refcat = greference (&Nrefcat);
    3436  project_refcat (refcat, Nrefcat);
     37  if ((DUMP != NULL) && !strcmp (DUMP, "refcat")) dump_refcat(refcat, Nrefcat);
    3538
    3639  match (refcat, Nrefcat);
    37 
    38   if (1) {
    39     int i;
    40     FILE *f, *g;
    41     f = fopen ("ref.dat", "w");
    42     g = fopen ("raw.dat", "w");
    43     for (i = 0; i < Nchip; i++) {
    44       dump_stars (f, chip[i].raw, chip[i].Nmatch);
    45       dump_stars (g, chip[i].ref, chip[i].Nmatch);
    46     }
    47     fclose (f);
    48     fclose (g);
    49     exit (1);
    50   }
    51 
    5240  deproject_raw ();
    5341  project_ref ();
     42  Cerror   = GetScatter (&Nastro);
     43  fprintf (stderr, "scatter (raw)  : %f for %d stars\n", Cerror, Nastro);
     44  if ((DUMP != NULL) && !strcmp (DUMP, "rawmatch")) dump_match();
    5445
    5546  grad = GetGradients ();
    5647  FitGradients (grad);
     48  Cerror   = GetScatter (&Nastro);
     49  fprintf (stderr, "scatter (grad) : %f for %d stars\n", Cerror, Nastro);
     50  if ((DUMP != NULL) && !strcmp (DUMP, "fitgrads")) dump_match();
    5751
    58   deproject_raw ();
    59   project_ref ();
    60 
    61   FitChips (1);
    62 
    63   deproject_raw ();
    64   project_ref ();
     52  FitChips (ChipOrder);
     53  Cerror   = GetScatter (&Nastro);
     54  fprintf (stderr, "scatter (chip) : %f for %d stars\n", Cerror, Nastro);
     55  if ((DUMP != NULL) && !strcmp (DUMP, "fitchips_unclip")) dump_match();
    6556
    6657  ClipOnFP (3.0);
     58  Cerror   = GetScatter (&Nastro);
     59  fprintf (stderr, "scatter (clip) : %f for %d stars\n", Cerror, Nastro);
    6760
     61  FitChips (ChipOrder);
     62  Cerror   = GetScatter (&Nastro);
     63  fprintf (stderr, "scatter (chip) : %f for %d stars\n", Cerror, Nastro);
     64  if ((DUMP != NULL) && !strcmp (DUMP, "fitchips")) dump_match();
     65
     66  output (argv[2], argv[3]);
     67  exit (0);
     68
     69/*
    6870  FitField (3);
     71  FitChips (3);
     72*/
    6973
    70   deproject_raw ();
    71   project_ref ();
    72 
    73   FitChips (3);
    74 
    75   deproject_raw ();
    76   project_ref ();
    77 
    78   wstars ();
    79   exit (0);
    8074}
    81 
    82 /** bore site terms : field.coords.crval1, field.coords.crval2, field.coords.pcN_M **/
    83 
    84 
    85   /* on the first pass, we only have the rough R,D for each input star */
    86   /* on the second pass, we have the reference R,D for each input star */
    87 
    88   /* pt[0].R,D here are from the chip coords */
    89   /* pt[0].R,D here are USNO coords */

  • trunk/Ohana/src/mosastro/src/warptest.c

    r3295 r3323  
    3838  coords.Npolyterms = 3;
    3939  for (i = 0; i < 7; i++) {
    40     coords.polyterms[0][i] = 0;
    41     coords.polyterms[1][i] = 0;
     40    coords.polyterms[i][0] = 0;
     41    coords.polyterms[i][1] = 0;
    4242  }
    4343  coords.polyterms[3][0] = 1e-10;

  • trunk/Ohana/src/mosastro/src/wstars.c

    r3294 r3323  
    11# include "mosastro.h"
    22# define NCHAR 66 /* 65 char EXCLUDING return */
     3
     4void wchip (char *filename, Chip *data) {
     5 
     6  FILE *g;
     7
     8  g = fopen (filename, "w");
     9  if (g == (FILE *) NULL) {
     10    fprintf (stderr, "ERROR: can't create output file %s\n", filename);
     11    exit (1);
     12  }
     13
     14  fwrite (data[0].header.buffer, 1, data[0].header.size, g);
     15  fwrite (data[0].buffer, 1, data[0].Nbuffer, g);
     16  fclose (g);
     17
     18}
    319
    420void wstars (char *filename, Stars *stars, int Nstars, Header *header) {
Note: See TracChangeset for help on using the changeset viewer.