IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37523


Ignore:
Timestamp:
Oct 22, 2014, 10:24:44 AM (12 years ago)
Author:
eugene
Message:

working on fakeastro; still not ready

Location:
branches/eam_branches/ipp-20140904/Ohana/src/fakeastro
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/Makefile

    r37481 r37523  
    5353$(SRC)/save_fake_stars.$(ARCH).o \
    5454$(SRC)/match_fake_stars.$(ARCH).o \
     55$(SRC)/fakeastro_2mass.$(ARCH).o \
     56$(SRC)/make_2mass_measures.$(ARCH).o \
    5557$(SRC)/remote_hosts.$(ARCH).o
    5658

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/include/fakeastro.h

    r37481 r37523  
    66# include <pthread.h>
    77
    8 typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES} FakeastroOp;
     8# define RESETTIME { gettimeofday (&startTimer, (void *) NULL); }
     9
     10typedef enum {OP_NONE, OP_GALAXY, OP_IMAGES, OP_2MASS} FakeastroOp;
    911
    1012typedef struct {
    … …  
    3032  Average  average;
    3133  Measure  measure;
     34  StarPar  starpar;
    3235  int found;
    3336} Stars;
    … …  
    9497float FAKEASTRO_RGAL; // parsec
    9598char FAKEASTRO_REF_EPOCH[80];
     99char FAKEASTRO_2MASS_EPOCH[80];
    96100
    97101float  RADIUS;
    … …  
    188192FakeAstro_Stars *make_fakeqsos (int Nstars, int ICRF);
    189193
     194int fakeastro_2mass ();
     195int make_2mass_measures (Catalog *catalog);

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/ConfigInit.c

    r37481 r37523  
    5757    strcpy (FAKEASTRO_REF_EPOCH, "2011/05/11,00:00:00"); // epoch of reference catalog
    5858  }
     59  if (!ScanConfig (config, "FAKEASTRO_2MASS_EPOCH", "%s", 0, FAKEASTRO_2MASS_EPOCH)) {
     60    strcpy (FAKEASTRO_2MASS_EPOCH, "2000/01/01,00:00:00"); // epoch of reference catalog
     61  }
    5962
    6063  /* set the default search radius */

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro.c

    r37496 r37523  
    1111    case OP_GALAXY:
    1212      /* the object analysis is a separate process iterating over catalogs */
    13       fakeastro_galaxy ();
     13      fakeastro_galaxy ();  // generate fake stars following a galaxy model
     14      /* make_fakestars() -- generates FakeAstro_Stars with density,distance and uR,uD following galaxy model
     15         make_fakeqsos()
     16         sortStars()
     17         make_subset()
     18         save_fakestars()
     19           fakestar_catalog()
     20             insert_fakestar() -- converts suppied FakeAstro_Stars and populates average/secfilt/starpar structures
     21           (fakestar_save_stars) [fakeastro_client / fakestar_load_stars]
     22       */
    1423      exit (0);
    1524
    1625    case OP_IMAGES:
    1726      fakeastro_images ();
     27      /* fakeastro_images_region()
     28           load_fake_stars()
     29           make_fake_images()
     30           make_fake_stars()
     31           fit_fake_stars()
     32           save_fake_stars()
     33             match_fake_stars()
     34       */
    1835      exit (0);
    1936
    2037    case OP_2MASS:
    2138      fakeastro_2mass ();
     39      /* make_2mass_measures()
     40       */
    2241      exit (0);
    2342

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakeastro_2mass.c

    r37496 r37523  
    88  SkyTableSetFilenames (skyTable, CATDIR, "cpt");
    99
    10   SkyList *skyList  = SkyListByPatch (skyTable, -1, &UserPatch);
     10  SkyList *skylist  = SkyListByPatch (skyTable, -1, &UserPatch);
    1111
    1212  Catalog catalog;
    … …  
    2121    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    2222    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    23     catalog.catflags  = LOAD_AVES | LOAD_SECF | LOAD_MEAS  | LOAD_STARPAR;
     23    catalog.catflags  = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR;
    2424    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    2525    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "r")) {
    … …  
    4343    dvo_catalog_free (&catalog);
    4444    if (VERBOSE) MARKTIME ("save cpt: %f sec\n", dtime); RESETTIME;
    45    
    46 
     45  }
    4746
    4847  exit (0);

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fakestar_io.c

    r37449 r37523  
    2727  gfits_define_bintable_column (&theader, "D", "RA",       "", "degree", 1.0, 0.0);
    2828  gfits_define_bintable_column (&theader, "D", "DEC",      "", "degree", 1.0, 0.0);
    29   gfits_define_bintable_column (&theader, "D", "GLON",     "", "degree", 1.0, 0.0);
    30   gfits_define_bintable_column (&theader, "D", "GLAT",     "", "degree", 1.0, 0.0);
     29  gfits_define_bintable_column (&theader, "E", "GLON",     "", "degree", 1.0, 0.0);
     30  gfits_define_bintable_column (&theader, "E", "GLAT",     "", "degree", 1.0, 0.0);
    3131  gfits_define_bintable_column (&theader, "E", "Ebv",      "", "", 1.0, 0.0);
    3232  gfits_define_bintable_column (&theader, "E", "dEbv",     "", "", 1.0, 0.0);
    … …  
    4444  double       *ra       = NULL; ALLOCATE (ra      , double      , Nstars); for (i = 0; i < Nstars; i++) ra      [i] = stars[i].R               ;
    4545  double       *dec      = NULL; ALLOCATE (dec     , double      , Nstars); for (i = 0; i < Nstars; i++) dec     [i] = stars[i].D               ;
    46   double       *glon     = NULL; ALLOCATE (glon    , double      , Nstars); for (i = 0; i < Nstars; i++) glon    [i] = stars[i].starpar.galLon  ;
    47   double       *glat     = NULL; ALLOCATE (glat    , double      , Nstars); for (i = 0; i < Nstars; i++) glat    [i] = stars[i].starpar.galLat  ;
     46  float        *glon     = NULL; ALLOCATE (glon    , float       , Nstars); for (i = 0; i < Nstars; i++) glon    [i] = stars[i].starpar.galLon  ;
     47  float        *glat     = NULL; ALLOCATE (glat    , float       , Nstars); for (i = 0; i < Nstars; i++) glat    [i] = stars[i].starpar.galLat  ;
    4848  float        *Ebv      = NULL; ALLOCATE (Ebv     , float       , Nstars); for (i = 0; i < Nstars; i++) Ebv     [i] = stars[i].starpar.Ebv     ;
    4949  float        *dEbv     = NULL; ALLOCATE (dEbv    , float       , Nstars); for (i = 0; i < Nstars; i++) dEbv    [i] = stars[i].starpar.dEbv    ;
    … …  
    135135  GET_COLUMN(ra      , "RA"        ,   double);
    136136  GET_COLUMN(dec     , "DEC"       ,   double);
    137   GET_COLUMN(glon    , "GLON"      ,   double);
    138   GET_COLUMN(glat    , "GLAT"      ,   double);
     137  GET_COLUMN(glon    , "GLON"      ,   float);
     138  GET_COLUMN(glat    , "GLAT"      ,   float);
    139139  GET_COLUMN(Ebv     , "Ebv"       ,   float);
    140140  GET_COLUMN(dEbv    , "dEbv"      ,   float);
    … …  
    159159    stars[i].R                = ra      [i];
    160160    stars[i].D                = dec     [i];
     161    stars[i].starpar.R        = ra      [i];
     162    stars[i].starpar.D        = dec     [i];
    161163    stars[i].starpar.galLon   = glon    [i];
    162164    stars[i].starpar.galLat   = glat    [i];

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/fit_fake_stars.c

    r37478 r37523  
    1212  if (Nstars > 225) order_use = 5; //  9 stars per cell
    1313  if (Nstars > 360) order_use = 6; // 10 stars per cell
     14  if (order_use < 1) {
     15    fprintf (stderr, "insufficient measurements (%d) for linear fit\n", Nstars);
     16    image[0].flags |= ID_IMAGE_ASTROM_FEW;
     17    return FALSE;
     18  }
    1419
    1520  // I have Nstars with "reference" positions stars[].Rref,Dref and
    … …  
    1924
    2025  int CHIPMAP = (image[0].coords.Npolyterms == -1);
    21   image[0].coords.Npolyterms = CHIPMAP ? -1 : order_use;
     26  image[0].coords.Npolyterms = CHIPMAP ? 1 : order_use;
    2227
    2328  CoordFit *fit = fit_init (image[0].coords.Npolyterms);

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_2mass_measures.c

    r37496 r37523  
    33// region corresponds to the catalog
    44int make_2mass_measures (Catalog *catalog) {
     5
     6  int Nsecfilt = GetPhotcodeNsecfilt ();
    57
    68  Average *average = catalog->average;
    … …  
    911  Measure *measure = catalog->measure;
    1012 
     13  time_t timeRef = ohana_date_to_sec(FAKEASTRO_REF_EPOCH);
     14  time_t tzero_2mass = ohana_date_to_sec(FAKEASTRO_2MASS_EPOCH);
     15
     16  // tzero, timeRef are in UNIX seconds, Toffset should be in years
     17  double Toffset = (tzero_2mass - timeRef) / 365.25 / 86400.0;
     18
     19  // use photcode to get zero point
     20  PhotCode *codeJ = GetPhotcodebyName ("2MASS_J");
     21  // PhotCode *codeH = GetPhotcodebyName ("2MASS_H");
     22  // PhotCode *codeK = GetPhotcodebyName ("2MASS_K");
     23
     24  float ZP = 21.0; // this is chosen to give 16.0 mag stars a flux of 100.0 counts (crude, yes)
     25  float SkyCts = 300.0; // this is chosen to make SN @ 16.0 = 5.0 (
     26
     27  int Nmeasure = catalog->Nmeasure;
     28  int NMEASURE = catalog->Nmeasure;
     29
     30  int Nsec = 4; // y-band
     31
    1132  int i;
    1233  for (i = 0; i < catalog->Naverage; i++) {
    1334
    14     int nStar = average[i].starparOffset;
     35    if (Nmeasure >= NMEASURE) {
     36      NMEASURE = Nmeasure + 1000;
     37      REALLOCATE (catalog[0].measure, Measure, NMEASURE);
     38      measure = catalog[0].measure;
     39    }
    1540
    16     // which filter?
    17     double y_PS1 = secfilt[i*Nsecfilt + 4].M  - ZP;
     41    int Nstarpar = average[i].starparOffset;
     42
     43    // make a crude JHK color for now:
     44    double J_PS1 = secfilt[i*Nsecfilt + Nsec].M - 1.0; // make all stars have y-J = 1.0
     45    if (J_PS1 > 16.0) continue; // only generate 2MASS detections for objects with J_PS1 < 16.0
     46
     47    double Minst = J_PS1 - ZP;
     48
    1849    double Counts = pow(10.0, -0.4*Minst);
    1950
    … …  
    2152
    2253    // true position from src catalog
    23     // XXX need to save the true position somewhere...
    24     double Rtru = average[i].R;
    25     double Dtru = average[i].D;
     54    double Rtru = starpar[Nstarpar].R;
     55    double Dtru = starpar[Nstarpar].D;
    2656
    2757    // observed position is scattered from true position by:
    2858    // * proper motion
    2959    // * gaussian scatter (~ seeing)
    30     double uR = starpar[nStar].uRA; // starpar are in arcsec / year
    31     double uD = starpar[nStar].uDEC;
     60    double uR = starpar[Nstarpar].uRA; // starpar are in arcsec / year
     61    double uD = starpar[Nstarpar].uDEC;
    3262   
    33     // tzero, timeRef are in UNIX seconds, Toffset should be in years
    34     double Toffset = (tzero_2mass - timeRef) / 365.25 / 86400.0;
    35 
    3663    // uR,uD in linear (arcsec / yr)
    3764    double dRpm = uR*Toffset;
    … …  
    4875    double Dobs = Dtru + dDoff;
    4976
    50     measure[nMeasure].M      = Minst + ZP;
    51     measure[nMeasure].dM     = 1.0 / SN;
     77    dvo_measure_init (&measure[Nmeasure]);
    5278
    53     // stars[Nstars].measure.dMcal      = ps1data[i].dMcal;
    54     measure[nMeasure].Sky        = sky[Nsec];
    55     measure[nMeasure].dSky       = sqrt(sky[Nsec]);
    56     measure[nMeasure].photFlags  = 0;
    57     measure[nMeasure].photFlags2 = 0;
    58     measure[nMeasure].R = Robs;
    59     measure[nMeasure].D = Dobs;
    60     measure[nMeasure].airmass = airmass (image->secz, stars[Nstars].average.R, stars[Nstars].average.D, image->sidtime, image->latitude);
    61     measure[nMeasure].az      = azimuth (15.0*image->sidtime - stars[Nstars].average.R, stars[Nstars].average.D, image->latitude);
    62     measure[nMeasure].Mcal    = image->Mcal;
    63     measure[nMeasure].t       = image->tzero + 1e-4*stars[Nstars].measure.Yccd*image->trate;  // trate is in 0.1 msec / row
    64     measure[nMeasure].dt      = Mtime;
     79    measure[Nmeasure].R = Robs;
     80    measure[Nmeasure].D = Dobs;
    6581
    66     measure[nMeasure].imageID = 0;
     82    measure[Nmeasure].M      = J_PS1;
     83    measure[Nmeasure].dM     = 1.0 / SN;
     84
     85    measure[Nmeasure].Sky        = SkyCts;
     86    measure[Nmeasure].dSky       = sqrt(SkyCts);
     87    measure[Nmeasure].photFlags  = 0;
     88    measure[Nmeasure].photFlags2 = 0;
     89    measure[Nmeasure].airmass = 1.0;
     90    measure[Nmeasure].az      = 0.0; // irrelevant
     91    measure[Nmeasure].Mcal    = 0.0;
     92    measure[Nmeasure].t       = tzero_2mass;
     93    measure[Nmeasure].dt      = 0.0;
     94    measure[Nmeasure].photcode = codeJ->code;
     95
     96    measure[Nmeasure].averef   = i;
     97    measure[Nmeasure].objID    = average[i].objID;
     98    measure[Nmeasure].catID    = average[i].catID;
     99
     100    measure[Nmeasure].imageID = 0;
    67101
    68102    // This is may optionally be replaced by the internal sequence (see FilterStars.c)
    69     measure[nMeasure].detID      = Nstars + 1;
    70     Nstars ++;
     103    measure[Nmeasure].detID      = Nmeasure;
     104    Nmeasure ++;
    71105  }
     106
     107  if (VERBOSE) fprintf (stderr, "added %d 2mass entries to %d stars\n",  (int) Nmeasure, (int) catalog->Naverage);
     108
     109  catalog->Nmeasure = Nmeasure;
     110  return TRUE;
    72111}

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fake_stars_catalog.c

    r37478 r37523  
    7171    double Dtru = average[i].D;
    7272
     73    stars[Nstars].starpar = starpar[nStar]; // keep this so we can save it on the generated stars
     74
    7375    // XXX what is the epoch for the src catalog? J2000?
    7476    stars[Nstars].Rref = Rtru; // add in any deviations we want here
    … …  
    113115    // stars[Nstars].measure.pltscale   = ps1data[i].pltscale;
    114116
    115     stars[Nstars].measure.M      = Minst + ZP;
     117    stars[Nstars].measure.M      = Minst + ZP; // XX I need to compensate for the internal zero point of 25.0
    116118    stars[Nstars].measure.dM     = 1.0 / SN;
    117119

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakeqsos.c

    r37481 r37523  
    3636    stars[i].found = FALSE;
    3737
     38    stars[i].starpar.R      = R;
     39    stars[i].starpar.D      = D;
    3840    stars[i].starpar.galLon = NAN;
    3941    stars[i].starpar.galLat = NAN;

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/make_fakestars.c

    r37478 r37523  
    9696    stars[i].found = FALSE;
    9797
     98    stars[i].starpar.R      = R;
     99    stars[i].starpar.D      = D;
    98100    stars[i].starpar.galLon = L;
    99101    stars[i].starpar.galLat = B;

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/match_fake_stars.c

    r37464 r37523  
    88  double dX, dY, dR;
    99  off_t *N1, *N2;
    10   off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
     10  off_t Nave, NAVE, Nmeas, NMEAS, Nstarpar, NSTARPAR, Nmatch;
    1111  unsigned int objID, catID;
    1212  Coords tcoords;
    … …  
    3535  Nmatch = 0;
    3636  NMEAS = Nmeas = catalog[0].Nmeasure;
     37  NSTARPAR = Nstarpar = catalog[0].Nstarpar;
    3738
    3839  // current max obj ID for this catalog
    … …  
    221222      REALLOCATE (catalog[0].measure, Measure, NMEAS);
    222223    }
     224    /* make sure there is space for next entry */
     225    if (Nstarpar >= NSTARPAR - 1) {
     226      NSTARPAR = Nstarpar + 1000;
     227      REALLOCATE (catalog[0].starpar, StarPar, NSTARPAR);
     228    }
    223229    if (Nave >= NAVE) {
    224230      NAVE = Nave + 1000;
    … …  
    231237    catalog[0].average[Nave].D             = stars[i].average.D;
    232238
     239    catalog[0].average[Nave].Nstarpar      = 1;
     240    catalog[0].average[Nave].starparOffset = Nstarpar;
    233241    catalog[0].average[Nave].Nmeasure      = 1;
    234242    catalog[0].average[Nave].measureOffset = Nmeas;
    … …  
    260268    }
    261269   
     270    // supply the starpar values from this detection
     271    dvo_starpar_init (&catalog[0].starpar[Nstarpar]);
     272    catalog[0].starpar[Nstarpar] = stars[i].starpar;
     273    Nstarpar ++;
     274
    262275    /* next[Nmeas] should always be -1 in this context (it is always the only
    263276       measurement for the star) */
    … …  
    269282  REALLOCATE (catalog[0].average, Average, Nave);
    270283  REALLOCATE (catalog[0].measure, Measure, Nmeas);
     284  REALLOCATE (catalog[0].starpar, StarPar, Nstarpar);
    271285 
    272286  catalog[0].sorted = FALSE;
    … …  
    283297  catalog[0].Naverage = Nave;
    284298  catalog[0].Nmeasure = Nmeas;
     299  catalog[0].Nstarpar = Nstarpar;
    285300  catalog[0].Nsecf_mem = Nave*Nsecfilt;
    286301  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);

  • branches/eam_branches/ipp-20140904/Ohana/src/fakeastro/src/save_fake_stars.c

    r37464 r37523  
    11# include "fakeastro.h"
    2 
    3 # define RESETTIME { gettimeofday (&startTimer, (void *) NULL); }
    42
    53int save_fake_stars (SkyTable *sky, Image *image, Stars *stars, int Nstars) {
    … …  
    2523    catalog.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
    2624    catalog.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
    27     catalog.catflags  = LOAD_AVES | LOAD_SECF;
     25    catalog.catflags  = LOAD_AVES | LOAD_SECF | LOAD_MEAS | LOAD_STARPAR;
    2826    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
    2927
    … …  
    4139
    4240    SetProtect (TRUE);
    43     catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
    4441    dvo_catalog_update (&catalog, VERBOSE);
    4542    SetProtect (FALSE);
Note: See TracChangeset for help on using the changeset viewer.