IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 5, 2010, 10:04:42 PM (16 years ago)
Author:
eugene
Message:

handle errors in dvo I/O functions; better support for DIS images; merge in high-speed mods from Niall; systematic errors in astrometry (somewhat hackish); add fits I/O option to succeed reading a short file (padding as needed); various kapa improvements (frames consistent in X, PS, PNG; fix tick label signif digits; add function for image overlays; better thick box lines; PNG of the image); bDrawBuffer mods to drop static buffer and refer to passed through buffer; ResizeByImage; DefineSectionByImage; buutos to do png & jpeg; add program roc for raid over cluster; support for mosaic photcodes in dvo image plot; densify; section name [-image x y] : width based on current image; resize -by-image; threaded addstar -resort; threaded dvomerge; dvorepair; dvoverify; dvomerge continue; dvomerge from list (of regions)

Location:
trunk/Ohana/src/relastro/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/src/GetAstromError.c

    r27581 r29938  
    22
    33float GetAstromError (Measure *measure, int mode) {
    4 
     4  //BIG HACKXXXXXXXX
     5  return 0.1;
    56  PhotCode *code;
    67  float dPobs, dPsys, dPtotal, dM, AS, MS;
    7 
    88  switch (mode) {
    99    case ERROR_MODE_RA:
    … …  
    1919      abort();
    2020  }
    21 
    2221  /* the astrometric errors are not being carried yet (but should be!) */
    2322  /* we use the photometric mag error as a weighting term */
    … …  
    2827  dPsys = code[0].astromErrSys;
    2928  dM    = measure[0].dM;
    30  
    3129  dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
    3230
    … …  
    3533
    3634  dPtotal = MAX (dPtotal, MIN_ERROR);
    37 
    3835  return (dPtotal);
    3936}

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r28241 r29938  
    5454  memset (&fitPM,  0, sizeof(fitPM));
    5555  memset (&fitPAR, 0, sizeof(fitPAR));
    56  
    5756  initObjectData (catalog, Ncatalog);
    5857
    … …  
    7271  // use J2000 as a reference time
    7372  T2000 = ohana_date_to_sec ("2000/01/01");
    74 
    7573  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
    7674  Nsecfilt = GetPhotcodeNsecfilt ();
    … …  
    9795      N = 0;
    9896      m = catalog[i].average[j].measureOffset;
    99 
    10097      Tmin = Tmax = (catalog[i].measure[m].t - T2000) / (86400*365.25);
    10198      mode = FIT_MODE;
    … …  
    137134
    138135        // dX, dY : error in arcsec --
    139         dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
    140         dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
     136        // dX[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
     137        // dY[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
     138
     139        dX[N] = 0.1;
     140        dY[N] = 0.1;
    141141        dT[N] = catalog[i].measure[m].dt;
    142142
    … …  
    145145        // non-circular errors (different values for X and Y), then dR and dD will be
    146146        // incorrect: they would need to be rotated to take out the position angle
    147         dR[k] = dX[k] / 3600.0;
    148         dD[k] = dY[k] / 3600.0;
     147        dR[N] = dX[N] / 3600.0;
     148        dD[N] = dY[N] / 3600.0;
    149149
    150150        N++;
    … …  
    194194
    195195        FitPM (&fitPM, X, dX, Y, dY, T, N);
     196
    196197        if (XVERB) fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
    197198
    … …  
    244245
    245246      switch (result) {
    246         case FIT_AVERAGE:
    247           catalog[i].average[j].flags |= ID_STAR_USE_AVE;
    248           fit = fitAve;
    249           break;
    250         case FIT_PM_ONLY:
    251           catalog[i].average[j].flags |= ID_STAR_USE_PM;
    252           fit = fitPM;
    253           break;
    254         case FIT_PM_AND_PAR:
    255           catalog[i].average[j].flags |= ID_STAR_USE_PAR;
    256           fit = fitPAR;
    257           break;
    258       }
    259 
    260       if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
     247      case FIT_AVERAGE:
     248        catalog[i].average[j].flags |= ID_STAR_USE_AVE;
     249        fit = fitAve;
     250        break;
     251      case FIT_PM_ONLY:
     252        catalog[i].average[j].flags |= ID_STAR_USE_PM;
     253        fit = fitPM;
     254        break;
     255      case FIT_PM_AND_PAR:
     256        catalog[i].average[j].flags |= ID_STAR_USE_PAR;
     257        fit = fitPAR;
     258        break;
     259      }
     260      if (XVERB) fprintf (stderr, "A%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
    261261                          catalog[i].average[j].R,
    262262                          catalog[i].average[j].D,
    263263                          fit.Ro, fit.Do,
    264264                          3600*(catalog[i].average[j].R - fit.Ro),
    265                           3600*(catalog[i].average[j].D - fit.Do));
     265                          3600*(catalog[i].average[j].D - fit.Do),
     266                          fit.uR,
     267                          fit.uD);
    266268
    267269      //make sure that the fit succeeded
    268       status  = finite(fit.Ro);
     270      status = TRUE;
     271      status &= finite(fit.Ro);
    269272      status &= finite(fit.Do);
    270273      status &= finite(fit.dRo);
    … …  
    310313      catalog[i].average[j].Trange = (Trange * 86400 * 365.26);
    311314      catalog[i].average[j].Npos = fit.Nfit;
     315      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
     316                          catalog[i].average[j].R,
     317                          catalog[i].average[j].D,
     318                          fit.Ro, fit.Do,
     319                          3600*(catalog[i].average[j].R - fit.Ro),
     320                          3600*(catalog[i].average[j].D - fit.Do),
     321                          catalog[i].average[j].uR,
     322                          catalog[i].average[j].uD);
     323
    312324    }
    313325

  • trunk/Ohana/src/relastro/src/high_speed_objects.c

    r28811 r29938  
    1313  // XXX seems silly to require region here just to find the catalog bounds / center
    1414
    15   off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad;
    16   int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid;
     15  off_t i, j, m, J, ni, nj, *N1, Nslow, Ninvalid, NgroupA, NgroupB, NgroupAbad, NgroupBbad,l,i1,Noff,Nmeasure,Naverage, NAVERAGE, NMEASURE;
     16  int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid,Nmatchmeas,Nepoch,nv[2],Nmatchmeasobj;
    1717  double *X1, *Y1;
    18   double dX, dY, dR, RADIUS2, yJ;
     18  double dX, dY, dR, RADIUS2;
    1919  Coords tcoords;
     20  Catalog catalog1;
    2021
    2122  int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt;
     23  char filename[1024];
     24  char outdir[]="/data/ipp022.0/ndeacon/hispeedzy";
     25  Noff = strlen(CATDIR);
     26  sprintf (filename, "%s/%s", outdir, &catalog[0].filename[Noff]);
     27  dvo_catalog_init(&catalog1, TRUE); /*initialise new catalogue*/
     28  catalog1.filename = strcreate(filename);
     29
     30  SIGMA_LIM = 0.0;
    2231
    2332  Nsecfilt = GetPhotcodeNsecfilt();
     33  Naverage=catalog[0].Naverage;
     34  Nmeasure=catalog[0].Nmeasure;
     35  // ALLOCATE (catalog1.average, Average,Naverage);
     36  // ALLOCATE (catalog1.measure, Measure,Nmeasure);
     37  // ALLOCATE(catalog1.secfilt, SecFilt,catalog[0].Naverage*Nsecfilt);
     38  catalog1.filename = filename; // based on the input name, need to keep everything below the catdir portion
     39  catalog1.Nsecfilt  = Nsecfilt;
     40  // always load all of the data (if any exists)
     41  catalog1.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     42 
     43  catalog1.catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data
     44  catalog1.catmode   = dvo_catalog_catmode (CATMODE);      // set the default catmode from config data
     45
     46  if (!dvo_catalog_open (&catalog1, region, VERBOSE,"w")) {
     47    fprintf (stderr, "ERROR: failure to open catalog file %s\n",
     48             filename);
     49    exit (2);
     50  }
     51
     52  NAVERAGE = 1000;
     53  NMEASURE = 10000;
     54  REALLOCATE (catalog1.average, Average, NAVERAGE);
     55  REALLOCATE (catalog1.measure, Measure, NMEASURE);
     56  REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
    2457
    2558  ycode = GetPhotcodeCodebyName("y");
    … …  
    67100  NgroupBbad = 0;
    68101
    69   fprintf (stdout, "#      RA_A  DEC_A              RA_B  DEC_B      :     pmx     pmy      dR  :    z      dz       y      dy       J      dJ       H      dH       K      dK\n");
     102  fprintf (stdout, "#      RA_A  DEC_A              RA_B  DEC_B      :     pmx     pmy      dR  dt:    z      dz       y      dy       J      dJ       H      dH       K      dK\n");
    70103  //................270.2346670  20.2471540  270.2170434  20.2717396 :  -59.11   88.78   106.66 :   0.000  0.000   19.582  0.112   16.041  0.110   15.388  0.098   14.858  0.001
    71104
    … …  
    88121    foundA = FALSE;
    89122    for (j = 0; !foundA && (j < catalog[0].average[i].Nmeasure); j++, m++) {
     123      if((catalog[0].average[i].R>204.1923)&&(catalog[0].average[i].R<204.1924)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377))
     124        {
     125          printf("Hello");
     126        }
     127
    90128      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupA, NphotcodesGroupA)) {
    91129        foundA = TRUE;
    … …  
    97135    foundB = FALSE;
    98136    for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) {
     137                                                   
     138      if((catalog[0].average[i].R>204.192)&&(catalog[0].average[i].R<204.1925)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377))
     139        {
     140          printf("Hello");
     141        }
     142
    99143      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) {
    100144        foundB = TRUE;
    … …  
    114158    if (foundA && !foundB) {
    115159      // average-based tests:
     160                                                   
     161      if((catalog[0].average[i].R>204.1923)&&(catalog[0].average[i].R<204.1924)&&(catalog[0].average[i].D>11.376)&&(catalog[0].average[i].D<11.377))
     162        {
     163          printf("Hello");
     164        }
     165
    116166      valid = TRUE;
    117       valid &= ((catalog[0].average[i].flags & 0x40000) == 0);
    118       valid &= ((catalog[0].average[i].flags & 0x80000)  > 0);
    119       valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0);
    120       valid &= (catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0);
     167      valid &= ((catalog[0].average[i].flags & 0x02000000) == 0);
     168      valid &= ((catalog[0].average[i].flags & 0x08000000)  > 0);
     169      valid &= ((catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + yNsec].M)));
     170      valid &= (((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0))||(isnan(catalog[0].secfilt[i*Nsecfilt + zNsec].M)));
    121171      // XXX color restrictions are applied in the pair-wise matching below
    122172     
    123173      // measure-based tests:
    124       m = catalog[0].average[i].measureOffset;
     174      /*      m = catalog[0].average[i].measureOffset;
    125175      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    126176        if (catalog[0].measure[m].photcode == USNO_R) {
    … …  
    130180          valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M));
    131181        }
    132       }
     182        }*/
    133183
    134184      // ((objflags & 524288) == 524288) &&
    … …  
    154204
    155205      // average-based tests:
     206      if((catalog[0].average[i].R>204.192)&&(catalog[0].average[i].R<204.193)&&(catalog[0].average[i].D>11.372)&&(catalog[0].average[i].D<11.373))
     207        {
     208          printf("Hello");
     209        }
    156210      valid = TRUE;
    157       valid &= ((catalog[0].average[i].flags & 0x10000) == 0);
    158       valid &= ((catalog[0].average[i].flags & 0x20000)  > 0);
    159       valid &= (catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0);
     211      valid &= ((catalog[0].average[i].flags & 0x01000000) == 0);
     212
     213      valid &= ((catalog[0].average[i].flags & 0x04000000)  > 0);
     214
     215      valid &= ((catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0)||(isnan(catalog[0].secfilt[i*Nsecfilt + jNsec].M)));
    160216      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1);
    161217      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2);
    162218
    163       if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
     219      /*if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
    164220        valid &= TRUE;
    165221      } else {
    166222        valid &= ((catalog[0].secfilt[i*Nsecfilt + zNsec].M - catalog[0].secfilt[i*Nsecfilt + yNsec].M) > 0.2);
    167       }
     223        }*/
    168224
    169225      // measure-based tests:
    170       m = catalog[0].average[i].measureOffset;
     226      /*m = catalog[0].average[i].measureOffset;
    171227      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
    172228        if (catalog[0].measure[m].photcode == USNO_R) {
    … …  
    176232          valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M));
    177233        }
    178       }
     234        }*/
    179235
    180236      // (abs(glat)>15.0) &&
    … …  
    236292
    237293  Nmatch = 0;
     294  Nmatchmeas = 0;
    238295  RADIUS2 = SQ(RADIUS);
    239296
    … …  
    243300    ni = N1[i];
    244301    nj = N1[j];
    245 
    246302    // if (ni == 131030) {
    247303    //   fprintf (stderr, "got 2mass\n");
    … …  
    268324
    269325    // within match range; look for valid matches
    270     for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {
     326    for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {     
    271327      if (J == i) continue;  // avoid auto-matches
    272328
    … …  
    281337
    282338      // y_groupA - J_groupB
    283       yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
     339
     340      /*      yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
    284341      if (yJ < 1.25) continue;
    285       if (yJ > 5.0) continue;
     342      if (yJ > 5.0) continue;*/
    286343
    287344      /*** a match is found (just print it for the moment) ***/
     345      /*Define a vector of matched indices and set averages in new catalogue*/
     346      Nepoch=2;
     347      FIT_MODE = FIT_PM_ONLY;
     348/*nv=(int *)malloc(Nepoch*sizeof(int));*/
     349      nv[0]=ni; /*THESE SHOULD BE CHANGED FOR MULTIPLE EPOCHS AS SHOULD nv*/
     350      nv[1]=nj;
     351
     352      catalog1.average[Nmatch]=catalog[0].average[nv[0]];
     353      /*Loop over index values and set measurements in new catalogue*/
     354      Nmatchmeasobj=0;
     355      catalog1.average[Nmatch].measureOffset=Nmatchmeas;
     356      for(l=0;l<Nepoch;l++) {
     357          m = catalog[0].average[nv[l]].measureOffset;
     358          for(i1=0;i1<catalog[0].average[nv[l]].Nmeasure;i1++) {
     359              catalog1.measure[Nmatchmeas]=catalog[0].measure[m+i1];
     360              /*Set offset RA and Dec wrt correct average value*/
     361              catalog1.measure[Nmatchmeas].dR=catalog[0].measure[m+i1].dR+3600.0*(catalog[0].average[nv[0]].R-catalog[0].average[nv[l]].R);
     362              catalog1.measure[Nmatchmeas].dD=catalog[0].measure[m+i1].dD+3600.0*(catalog[0].average[nv[0]].D-catalog[0].average[nv[l]].D);
     363              catalog1.measure[Nmatchmeas].averef = Nmatch;
     364              Nmatchmeasobj++;
     365              Nmatchmeas++;
     366
     367              if (Nmatchmeas == NMEASURE - 1) {
     368                NMEASURE += 10000;
     369                REALLOCATE (catalog1.measure, Measure, NMEASURE);
     370              }
     371          }
     372      }
     373      catalog1.average[Nmatch].Nmeasure=Nmatchmeasobj;
    288374      Nmatch ++;
    289375
     376      if (Nmatch == NAVERAGE - 1) {
     377        NAVERAGE += 1000;
     378        REALLOCATE (catalog1.average, Average, NAVERAGE);
     379        REALLOCATE(catalog1.secfilt, SecFilt, NAVERAGE*Nsecfilt);
     380      }
    290381      // for the moment, just write the matches to stdout:
    291       float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
    292       float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
    293       fprintf (stdout, "%11.7f %11.7f  %11.7f %11.7f : %7.2f %7.2f  %7.2f : %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f\n",
     382      //float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
     383      //float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
     384      /*           fprintf (stdout, "%11.7f %11.7f  %11.7f %11.7f : %7.2f %7.2f  %7.2f: %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f\n",
    294385               catalog[0].average[ni].R, catalog[0].average[ni].D,
    295386               catalog[0].average[nj].R, catalog[0].average[nj].D,
    296                pmx, pmy, sqrt(dR),
     387               dX, dY, sqrt(dR),
    297388               catalog[0].secfilt[nj*Nsecfilt + zNsec].M,
    298389               catalog[0].secfilt[nj*Nsecfilt + zNsec].dM,
    … …  
    304395               catalog[0].secfilt[ni*Nsecfilt + hNsec].dM,
    305396               catalog[0].secfilt[ni*Nsecfilt + kNsec].M,
    306                catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);
     397               catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);*/
    307398    }
    308399    i++;
    309400  }
     401  catalog1.Naverage=Nmatch;
     402  catalog1.Nmeasure=Nmatchmeas;
     403  catalog1.Nsecfilt=Nsecfilt;
     404  catalog1.Nsecf_mem=Nmatch*Nsecfilt;
     405
     406  UpdateObjects(&catalog1, 1);
     407
    310408  fprintf (stderr, "found %d matches\n", Nmatch);
    311  
     409  dvo_catalog_save (&catalog1, VERBOSE);
     410  dvo_catalog_unlock (&catalog1);
     411  dvo_catalog_free (&catalog1);
    312412  free (slowMoving);
    313413  free (groupA);
Note: See TracChangeset for help on using the changeset viewer.