IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12319


Ignore:
Timestamp:
Mar 8, 2007, 11:38:26 AM (19 years ago)
Author:
eugene
Message:

basic object fitting works; images and pm/parallax not tested

Location:
branches/dvo-mods-2007-02/Ohana/src/relastro
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

  • branches/dvo-mods-2007-02/Ohana/src/relastro/Makefile

    r12276 r12319  
    5151$(SRC)/relastro.$(ARCH).o        \
    5252$(SRC)/reload_catalogs.$(ARCH).o \
     53$(SRC)/save_catalogs.$(ARCH).o   \
    5354$(SRC)/setExclusions.$(ARCH).o   \
    5455$(SRC)/write_coords.$(ARCH).o

  • branches/dvo-mods-2007-02/Ohana/src/relastro/include/relastro.h

    r12276 r12319  
    126126time_t TSTART, TSTOP;
    127127
    128 int FITTING_MODE;
    129 enum {FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR};
     128int FIT_MODE;
     129enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR};
     130
     131int FIT_TARGET;
     132enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS};
    130133
    131134SkyRegion UserPatch;
    … …  
    188191void          initstats           PROTO((char *mode));
    189192int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
    190 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog));
     193Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
    191194SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
    192195Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
    … …  
    279282double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
    280283double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);
     284int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
     285int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c

    r12220 r12319  
    4545  PMFit fit;
    4646
     47  initObjectData (catalog, Ncatalog);
     48
    4749  coords.crval1 = 0;
    4850  coords.crval2 = 0;
    … …  
    6062      /* calculate the average value of R,D for a single star */
    6163      if (catalog[i].average[j].code & STAR_BAD) continue; 
    62       m = catalog[i].average[j].offset;
    6364
    6465      N = 0;
     66      m = catalog[i].average[j].offset;
    6567      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
    6668        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
    … …  
    7880        N++;
    7981      }
     82
     83      // XXX This criterion needs to be better considered: adjust to match Ndof
    8084      if (N < STAR_TOOFEW) { /* too few measurements */
    8185        catalog[i].average[j].code |= ID_STAR_FEW;
    … …  
    9599
    96100      /* fit the model components as needed */
    97       switch (FITTING_MODE) {
     101      switch (FIT_MODE) {
    98102        case FIT_AVERAGE:
    99103          liststats (R, dR, N, &statsR);
    … …  
    121125          break;
    122126        default:
    123           fprintf (stderr, "invalid fitting mode %d\n", FITTING_MODE);
     127          fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__);
    124128          exit (2);
    125129      }   
     130
     131      if (0 && (j < 100)) {
     132          fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
     133                   catalog[i].average[j].R,
     134                   catalog[i].average[j].D,
     135                   fit.Ro, fit.Do,
     136                   3600*(catalog[i].average[j].R - fit.Ro),
     137                   3600*(catalog[i].average[j].D - fit.Do));
     138      }
     139
     140      // the measure fields must be updated before the average fields
     141      m = catalog[i].average[j].offset;
     142      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
     143        if (catalog[i].measure[m].flags & MEAS_BAD) continue;
     144        setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     145        setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
     146      }     
    126147
    127148      catalog[i].average[j].R   = fit.Ro;
    … …  
    139160
    140161      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
     162
    141163    }
    142164  }

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/args.c

    r12276 r12319  
    66  int N;
    77  double trange;
     8
     9  /* possible operations */
     10  FIT_TARGET = TARGET_NONE;
     11  FIT_MODE = FIT_AVERAGE;
     12  if ((N = get_argument (argc, argv, "-update-objects"))) {
     13    remove_argument (N, &argc, argv);
     14    FIT_TARGET = TARGET_OBJECTS;
     15
     16    // check for object fitting modes (not valid for images)
     17    if ((N = get_argument (argc, argv, "-pm"))) {
     18        remove_argument (N, &argc, argv);
     19        FIT_MODE = FIT_PM_ONLY;
     20    }
     21    if ((N = get_argument (argc, argv, "-pmpar"))) {
     22        remove_argument (N, &argc, argv);
     23        FIT_MODE = FIT_PM_AND_PAR;
     24    }
     25  }
     26  if ((N = get_argument (argc, argv, "-update-simple"))) {
     27    remove_argument (N, &argc, argv);
     28    FIT_TARGET = TARGET_SIMPLE;
     29  }
     30  if ((N = get_argument (argc, argv, "-update-chips"))) {
     31    remove_argument (N, &argc, argv);
     32    FIT_TARGET = TARGET_CHIPS;
     33  }
     34  if ((N = get_argument (argc, argv, "-update-mosaics"))) {
     35    remove_argument (N, &argc, argv);
     36    FIT_TARGET = TARGET_MOSAICS;
     37  }
     38  if (FIT_TARGET == TARGET_NONE) usage();
     39
     40  /* specify portion of the sky : allow default of all sky? */
     41  UserPatch.Rmin = 0;
     42  UserPatch.Rmax = 360;
     43  UserPatch.Dmin = -90;
     44  UserPatch.Dmax = +90;
     45  if ((N = get_argument (argc, argv, "-region"))) {
     46    remove_argument (N, &argc, argv);
     47    UserPatch.Rmin = atof (argv[N]);
     48    remove_argument (N, &argc, argv);
     49    UserPatch.Rmax = atof (argv[N]);
     50    remove_argument (N, &argc, argv);
     51    UserPatch.Dmin = atof (argv[N]);
     52    remove_argument (N, &argc, argv);
     53    UserPatch.Dmax = atof (argv[N]);
     54    remove_argument (N, &argc, argv);
     55  } else {
     56    usage ();
     57  }
    858
    959  /* define time */
    … …  
    3484  }
    3585
    36   /* specify portion of the sky : allow default of all sky? */
    37   UserPatch.Rmin = 0;
    38   UserPatch.Rmax = 360;
    39   UserPatch.Dmin = -90;
    40   UserPatch.Dmax = +90;
    41   if ((N = get_argument (argc, argv, "-region"))) {
    42     remove_argument (N, &argc, argv);
    43     UserPatch.Rmin = atof (argv[N]);
    44     remove_argument (N, &argc, argv);
    45     UserPatch.Rmax = atof (argv[N]);
    46     remove_argument (N, &argc, argv);
    47     UserPatch.Dmin = atof (argv[N]);
    48     remove_argument (N, &argc, argv);
    49     UserPatch.Dmax = atof (argv[N]);
    50     remove_argument (N, &argc, argv);
    51   } else {
    52     usage ();
    53   }
    54 
    5586  PHOTCODE_LIST = NULL;
    5687  if ((N = get_argument (argc, argv, "-photcode"))) {
    … …  
    79110  }
    80111
    81 # if (0)
    82   /* XXX is this still relevant?? */
    83112  strcpy (STATMODE, "CHI_INNER_WTMEAN");
    84113  if ((N = get_argument (argc, argv, "-statmode"))) {
    … …  
    87116    remove_argument (N, &argc, argv);
    88117  }
    89 # endif
    90118
    91119  RESET = FALSE;
    … …  
    153181    ImagSelect = TRUE;
    154182  }
    155 
    156   /* possible operations */
    157   DoUpdateObjects = FALSE;
    158   if ((N = get_argument (argc, argv, "-update-objects"))) {
    159     remove_argument (N, &argc, argv);
    160     DoUpdateObjects = TRUE;
    161   }
    162   DoUpdateSimple = FALSE;
    163   if ((N = get_argument (argc, argv, "-update-simple"))) {
    164     remove_argument (N, &argc, argv);
    165     DoUpdateSimple = TRUE;
    166   }
    167   DoUpdateChips = FALSE;
    168   if ((N = get_argument (argc, argv, "-update-chips"))) {
    169     remove_argument (N, &argc, argv);
    170     DoUpdateChips = TRUE;
    171   }
    172   DoUpdateMosaics = FALSE;
    173   if ((N = get_argument (argc, argv, "-update-mosaics"))) {
    174     remove_argument (N, &argc, argv);
    175     DoUpdateMosaics = TRUE;
    176   }
    177 
     183 
    178184  /* XXX drop this? */
    179185  DophotSelect = FALSE;

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/dvo_astrom_ops.c

    r12047 r12319  
    33double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) {
    44
    5     double ra;
     5  double ra;
    66
    7     /* the measure carries the instantaneous mean position at the epoch t */
    8     ra = average[0].R - measure[0].dR / 3600.0;
     7  /* the measure carries the instantaneous mean position at the epoch t */
     8  ra = average[0].R - measure[0].dR / 3600.0;
    99
    10     /* possible corrections to mean ra:
     10  /* possible corrections to mean ra:
    1111
    12     - proper-motion and parallax
    13     - abberation
    14     - precession and nutation, etc
    15     - refraction
    16     - DCR
     12  - proper-motion and parallax
     13  - abberation
     14  - precession and nutation, etc
     15  - refraction
     16  - DCR
    1717
    18     */
     18  */
    1919
    20     return (ra);
     20  return (ra);
    2121}
    2222
    2323double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) {
    2424
    25     double dec;
     25  double dec;
    2626
    27     /* the measure carries the instantaneous mean position at the epoch t */
    28     dec = average[0].D - measure[0].dD / 3600.0;
     27  /* the measure carries the instantaneous mean position at the epoch t */
     28  dec = average[0].D - measure[0].dD / 3600.0;
    2929
    30     /* possible corrections to mean ra:
     30  /* possible corrections to mean ra:
    3131
    32     - proper-motion and parallax
    33     - abberation
    34     - precession and nutation, etc
    35     - refraction
    36     - DCR
     32  - proper-motion and parallax
     33  - abberation
     34  - precession and nutation, etc
     35  - refraction
     36  - DCR
    3737
    38     */
     38  */
    3939
    40     return (dec);
     40  return (dec);
    4141}
     42
     43int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
     44
     45  /* math to get from new fitted position to new measure offset
     46     ra_obs = average[0].R - measure[0].dR / 3600.0;
     47     measure[0].dR = (ra_fit - ra_obs) * 3600.0;
     48     measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0
     49     measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR;
     50  */
     51
     52  /* the measure carries the instantaneous mean position at the epoch t */
     53  measure[0].dR += (ra_fit - average[0].R) * 3600.0;
     54
     55  /* possible corrections to mean ra:
     56
     57  - proper-motion and parallax
     58  - abberation
     59  - precession and nutation, etc
     60  - refraction
     61  - DCR
     62
     63  */
     64
     65  return (TRUE);
     66}
     67
     68int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) {
     69
     70  /* math to get from new fitted position to new measure offset
     71     dec_obs = average[0].D - measure[0].dD / 3600.0;
     72     measure[0].dD = (dec_fit - dec_obs) * 3600.0;
     73     measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0
     74     measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD;
     75  */
     76
     77  /* the measure carries the instantaneous mean position at the epoch t */
     78  measure[0].dD += (dec_fit - average[0].D) * 3600.0;
     79
     80  /* possible corrections to mean ra:
     81
     82  - proper-motion and parallax
     83  - abberation
     84  - precession and nutation, etc
     85  - refraction
     86  - DCR
     87
     88  */
     89
     90  return (TRUE);
     91}

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/load_catalogs.c

    r12048 r12319  
    11# include "relastro.h"
    22
    3 Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) {
     3Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
    44
    55  int i, Nstar;
    6   Catalog *catalog, tcatalog;
     6  Catalog *catalog, *pcatalog, tcatalog;
    77
    88  if (VERBOSE) fprintf (stderr, "loading catalog data\n");
    … …  
    1313  for (i = 0; i < skylist[0].Nregions; i++) {
    1414
     15    pcatalog = subselect ? &tcatalog : &catalog[i];
     16
    1517    // set up the basic catalog info
    16     tcatalog.filename = skylist[0].filename[i];
    17     tcatalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
    18     tcatalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
    19     tcatalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
     18    pcatalog[0].filename = skylist[0].filename[i];
     19    pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
     20    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
     21    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
    2022
    21     if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
    22       fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename);
     23    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
     24      fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
    2325      exit (1);
    2426    }
    25     if (VERBOSE && !tcatalog.Nave_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
     27    if (VERBOSE && !pcatalog[0].Nave_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
    2628
    2729    // select only the brighter stars
    28     bcatalog (&catalog[i], &tcatalog);
    29     dvo_catalog_unlock (&tcatalog);
    30     dvo_catalog_free (&tcatalog);
     30    if (subselect) {
     31      bcatalog (&catalog[i], &tcatalog);
     32      dvo_catalog_unlock (&tcatalog);
     33      dvo_catalog_free (&tcatalog);
     34    }
    3135  }
    3236
    … …  
    3741  }
    3842  if (Nstar < 2) {
    39     fprintf (stderr, "insufficient stars %d\n", Nstar);
    40     exit (0);
     43    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
    4144  }
    4245

  • branches/dvo-mods-2007-02/Ohana/src/relastro/src/relastro.c

    r12276 r12319  
    2424  skylist = load_images (&db, &UserPatch);
    2525
    26   /* load catalog data from region files */
    27   catalog = load_catalogs (skylist, &Ncatalog);
     26  /* load catalog data from region files : subselect only if we are not doing the objects */
     27  catalog = load_catalogs (skylist, &Ncatalog, (FIT_TARGET != TARGET_OBJECTS));
    2828
    2929  /* match measurements with images, mosaics */
    … …  
    4040
    4141  /* major modes */
    42   if (DoUpdateObjects) {
    43     UpdateObjects (catalog, Ncatalog);
    44   }
    45   if (DoUpdateSimple) {
    46     UpdateSimple (catalog, Ncatalog);
    47   }
    48   if (DoUpdateChips) {
    49     UpdateChips (catalog, Ncatalog);
    50   }
    51   if (DoUpdateMosaics) {
    52     UpdateMosaic (catalog, Ncatalog);
     42  switch (FIT_TARGET) {
     43    case TARGET_OBJECTS:
     44      UpdateObjects (catalog, Ncatalog);
     45      break;
     46    case TARGET_SIMPLE:
     47      UpdateSimple (catalog, Ncatalog);
     48      break;
     49
     50    case TARGET_CHIPS:
     51      UpdateChips (catalog, Ncatalog);
     52      break;
     53
     54    case TARGET_MOSAICS:
     55      UpdateMosaic (catalog, Ncatalog);
     56      break;
     57    default:
     58      fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__);
     59      exit (2);
    5360  }
    5461
    … …  
    6471  if (!UPDATE) exit (0);
    6572
    66   /* need to figure out how to update images, etc */
    67   dvo_image_update (&db, VERBOSE);
    68   dvo_image_unlock (&db);
     73  if (FIT_TARGET == TARGET_OBJECTS) {
     74    save_catalogs (catalog, Ncatalog);
     75  } else {
     76    dvo_image_update (&db, VERBOSE);
     77    dvo_image_unlock (&db);
     78  }
    6979
    7080  exit (0);
Note: See TracChangeset for help on using the changeset viewer.