IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jan 6, 2017, 11:30:10 AM (10 years ago)
Author:
eugene
Message:

merging changes from czw dev branch (compare with r39924)

Location:
trunk
Files:
7 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/Ohana

  • trunk/Ohana/src/libdvo/include/dvo.h

    r39670 r39926  
    114114/* Measure.flags values -- these values are 32 bit (as of PS1_V1) */
    115115typedef enum {
    116   ID_MEAS_NOCAL          = 0x00000001,  // detection ignored for this analysis (photcode, time range) -- internal only
    117   ID_MEAS_POOR_PHOTOM    = 0x00000002,  // detection is photometry outlier                                               
    118   ID_MEAS_SKIP_PHOTOM    = 0x00000004,  // detection was ignored for photometry measurement                               
    119   ID_MEAS_AREA           = 0x00000008,  // detection near image edge                                                 
    120   ID_MEAS_POOR_ASTROM    = 0x00000010,  // detection is astrometry outlier                                               
    121   ID_MEAS_SKIP_ASTROM    = 0x00000020,  // detection was ignored for astrometry measurement                               
    122   ID_MEAS_USED_OBJ       = 0x00000040,  // detection was used during update objects 
    123   ID_MEAS_USED_CHIP      = 0x00000080,  // detection was used during update chips (XXX this probably does not make it into the db)
    124   ID_MEAS_BLEND_MEAS     = 0x00000100,  // detection is within radius of multiple objects
    125   ID_MEAS_BLEND_OBJ      = 0x00000200,  // multiple detections within radius of object
    126   ID_MEAS_WARP_USED      = 0x00000400,  // measurement used to find mean warp photometry
    127   ID_MEAS_UNMASKED_ASTRO = 0x00000800,  // measurement was not masked in final astrometry fit
    128   ID_MEAS_BLEND_MEAS_X   = 0x00001000,  // detection is within radius of multiple objects across catalogs                   
    129   ID_MEAS_ARTIFACT       = 0x00002000,  // detection is thought to be non-astronomical                               
    130   ID_MEAS_SYNTH_MAG      = 0x00004000,  // magnitude is synthetic
    131   ID_MEAS_PHOTOM_UBERCAL = 0x00008000,  // externally-supplied zero point from ubercal analysis
    132   ID_MEAS_STACK_PRIMARY  = 0x00010000,  // this stack measurement is in the primary skycell
    133   ID_MEAS_STACK_PHOT_SRC = 0x00020000,  // this measurement supplied the stack photometry
    134   ID_MEAS_ICRF_QSO       = 0x00040000,  // this measurement is an ICRF reference position
    135   ID_MEAS_IMAGE_EPOCH    = 0x00080000,  // this measurement is registered to the image epoch (not tied to ref catalog epoch)
    136   ID_MEAS_PHOTOM_PSF     = 0x00100000,  // this measurement is used for the mean psf mag
    137   ID_MEAS_PHOTOM_APER    = 0x00200000,  // this measurement is used for the mean ap mag
    138   ID_MEAS_PHOTOM_KRON    = 0x00400000,  // this measurement is used for the mean kron mag
    139   ID_MEAS_MASKED_PSF     = 0x01000000,  // this measurement is masked based on IRLS weights for mean psf mag
    140   ID_MEAS_MASKED_APER    = 0x02000000,  // this measurement is masked based on IRLS weights for mean ap mag
    141   ID_MEAS_MASKED_KRON    = 0x04000000,  // this measurement is masked based on IRLS weights for mean kron mag
     116  ID_MEAS_NOCAL            = 0x00000001,  // detection ignored for this analysis (photcode, time range) -- internal only
     117  ID_MEAS_POOR_PHOTOM      = 0x00000002,  // detection is photometry outlier                                             
     118  ID_MEAS_SKIP_PHOTOM      = 0x00000004,  // detection was ignored for photometry measurement                             
     119  ID_MEAS_AREA             = 0x00000008,  // detection near image edge                                               
     120  ID_MEAS_POOR_ASTROM      = 0x00000010,  // detection is astrometry outlier                                             
     121  ID_MEAS_SKIP_ASTROM      = 0x00000020,  // detection was ignored for astrometry measurement                             
     122  ID_MEAS_USED_OBJ         = 0x00000040,  // detection was used during update objects 
     123  ID_MEAS_USED_CHIP        = 0x00000080,  // detection was used during update chips (XXX this probably does not make it into the db)
     124  ID_MEAS_BLEND_MEAS       = 0x00000100,  // detection is within radius of multiple objects
     125  ID_MEAS_BLEND_OBJ        = 0x00000200,  // multiple detections within radius of object
     126  ID_MEAS_WARP_USED        = 0x00000400,  // measurement used to find mean warp photometry
     127  ID_MEAS_UNMASKED_ASTRO   = 0x00000800,  // measurement was not masked in final astrometry fit
     128  ID_MEAS_BLEND_MEAS_X     = 0x00001000,  // detection is within radius of multiple objects across catalogs                 
     129  ID_MEAS_ARTIFACT         = 0x00002000,  // detection is thought to be non-astronomical                                     
     130  ID_MEAS_SYNTH_MAG        = 0x00004000,  // magnitude is synthetic
     131  ID_MEAS_PHOTOM_UBERCAL   = 0x00008000,  // externally-supplied zero point from ubercal analysis
     132  ID_MEAS_STACK_PRIMARY    = 0x00010000,  // this stack measurement is in the primary skycell
     133  ID_MEAS_STACK_PHOT_SRC   = 0x00020000,  // this measurement supplied the stack photometry
     134  ID_MEAS_ICRF_QSO         = 0x00040000,  // this measurement is an ICRF reference position
     135  ID_MEAS_IMAGE_EPOCH      = 0x00080000,  // this measurement is registered to the image epoch (not tied to ref catalog epoch)
     136  ID_MEAS_PHOTOM_PSF       = 0x00100000,  // this measurement is used for the mean psf mag
     137  ID_MEAS_PHOTOM_APER      = 0x00200000,  // this measurement is used for the mean ap mag
     138  ID_MEAS_PHOTOM_KRON      = 0x00400000,  // this measurement is used for the mean kron mag
     139  ID_MEAS_MASKED_PSF       = 0x01000000,  // this measurement is masked based on IRLS weights for mean psf mag
     140  ID_MEAS_MASKED_APER      = 0x02000000,  // this measurement is masked based on IRLS weights for mean ap mag
     141  ID_MEAS_MASKED_KRON      = 0x04000000,  // this measurement is masked based on IRLS weights for mean kron mag
     142  ID_MEAS_OBJECT_HAS_2MASS = 0x10000000,  // measurement comes from an object with 2mass data
     143  ID_MEAS_OBJECT_HAS_GAIA  = 0x20000000,  // measurement comes from an object with gaia data
     144  ID_MEAS_OBJECT_HAS_TYCHO = 0x40000000,  // measurement comes from an object with tycho data
    142145} DVOMeasureFlags;
    143146
    … …  
    233236  ID_SECF_STACK_BESTDET = 0x00008000, // PS1 stack best measurement is a detection (not forced)
    234237  ID_SECF_STACK_PRIMDET = 0x00010000, // PS1 stack primary measurement is a detection (not forced)
     238
     239  ID_SECF_HAS_SDSS      = 0x00100000, // this photcode has SDSS photometry
     240  ID_SECF_HAS_HSC       = 0x00200000, // this photcode has HSC  photometry
     241  ID_SECF_HAS_CFH       = 0x00400000, // this photcode has CFH  photometry (mostly Megacam)
     242  ID_SECF_HAS_DES       = 0x00800000, // this photcode has DES  photometry
     243
    235244  ID_SECF_OBJ_EXT       = 0x01000000, // extended in this band
    236245

  • trunk/Ohana/src/libdvo/src/ImageMetadataSelection.c

    r37807 r39926  
    5959}
    6060
     61// note that this function is called by the dvo clients and uses the 'metadata' table
     62// generated by the dvo shell for the mosaic images only.
    6163Coords *MatchMosaicMetadata (unsigned int imageID) {
    6264
    … …  
    7375  mosaic.crval2 = image[m].crval2;
    7476
     77  // note that image->theta is calculated based on pc1_1, pc1_2 when the metadata file is
     78  // generated
    7579  mosaic.pc1_1 =  cos(RAD_DEG*image[m].theta);
    7680  mosaic.pc1_2 =  sin(RAD_DEG*image[m].theta);

  • trunk/Ohana/src/libdvo/src/coordops.c

    r39608 r39926  
    8484
    8585  /** extra polynomial terms **/
     86  // for ZPN, these are used to modify the radial distance and not the X,Y coords
    8687  if ((coords[0].Npolyterms > 1) && (proj != PROJ_ZPN)) {
    8788    X2 = X*X;
    … …  
    126127    if (proj == PROJ_WRP) {
    127128      if (!coords->mosaic) {
    128         myAbort ("missing mosaic element");
     129        // myAbort ("missing mosaic element");
     130        *ra  = L;
     131        *dec = M;
    129132        return (FALSE);
    130133      }
    … …  
    157160          ctht = 0.0;
    158161        } else {
    159           T = DEG_RAD / R;
    160           stht =   T / sqrt ( 1.0 + T*T);
    161           ctht = 1.0 / sqrt ( 1.0 + T*T);
     162          // T = DEG_RAD / R; // T in 1/radians
     163          // stht =   T / sqrt ( 1.0 + T*T);
     164          // ctht = 1.0 / sqrt ( 1.0 + T*T);
     165
     166          T = RAD_DEG * R;
     167          stht = 1.0 / sqrt ( 1.0 + T*T);
     168          ctht =   T / sqrt ( 1.0 + T*T);
    162169        }
     170        break;
     171      case PROJ_SIN:
     172        // R = (180/pi) cos (theta)
     173        ctht = RAD_DEG * R;
     174        stht = sqrt (1 - ctht*ctht);
    163175        break;
    164176      case PROJ_STG:
    … …  
    167179        ctht = sqrt (1 - stht*stht);
    168180        break;
    169       case PROJ_SIN:
    170         // R = (180/pi) cos (theta)
    171         ctht = RAD_DEG * R;
    172         stht = sqrt (1 - ctht*ctht);
    173         break;
    174181      case PROJ_ARC:
    175182        // R = 90 - theta (degrees)
    … …  
    177184        stht = cos (RAD_DEG * R);
    178185        break;
    179 
    180       case PROJ_ZPN:
    181         // R = 90 - theta (degrees)
    182         // this is wrong because we are ignoring the distortion
    183         // XXX For now, just solve for terms up to n = 3
    184 
    185         // Ro = (pi/180)(90 - theta)
    186         // R = (180/pi)sum (P_i R^i)
    187 
    188         if (UKIRT_ONLY) {
    189           double Ro = RAD_DEG * R;
    190           double P1 = coords[0].polyterms[0][1];
    191           double P3 = coords[0].polyterms[0][3];
    192 
    193           // find the roots of f(gamma) = P1 gamma + P3 gamma^3 - Ro
    194           // starting guess for gamma is Ro / P1
    195           double gamma = Ro / P1;
    196 
    197           int i;
    198           for (i = 0; i < 5; i++) {
    199             double F = P1*gamma + P3*gamma*gamma*gamma - Ro;
    200             double dFdgamma = P1 + 3.0*P3*gamma*gamma;
    201 
    202             double gamma_new = gamma - F / dFdgamma;
    203             gamma = gamma_new;
    204           }
    205          
    206           double theta = 90.0 - gamma * DEG_RAD ;
    207           ctht = cos (RAD_DEG * theta);
    208           stht = sin (RAD_DEG * theta);
    209           break;
    210         }
    211186
    212187      case PROJ_ZEA:
    … …  
    220195        ctht = sqrt (1 - stht*stht);
    221196        break;
     197
     198      case PROJ_ZPN:
     199
     200        // the forward projection is:
     201        // theta = atan2(stht, ctht)
     202        // gamma = (pi/2 - theta) : theta in radians
     203        // Ro = sum (P_i gamma^i)
     204        // R  = (180/pi) Ro
     205
     206        // given R, we need to find theta:
     207        // Ro = R * (pi / 180) = sum (P_i gamma^i)
     208        // solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson
     209
     210        // use Ro to get a guess for gamma and iterate
     211
     212        {
     213          double Ro = RAD_DEG * R;
     214         
     215          // find the roots of f(gamma) - Ro = 0
     216          // starting guess for gamma is (Ro - P0) / P1
     217          double gamma = (Ro - coords[0].polyterms[0][0]) / coords[0].polyterms[1][0];
     218         
     219          int iter;
     220          for (iter = 0; iter < 5; iter++) {
     221           
     222            double Rc = 0.0; // this will hold the ander
     223            double dR = 0.0;
     224            for (int i = coords[0].Npolyterms - 1; i > 1; i--) {
     225              double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
     226              Rc = (Rc + Pi)*gamma;
     227              dR = (dR + i*Pi)*gamma;
     228            }
     229            double P0 = coords[0].polyterms[0][0];
     230            double P1 = coords[0].polyterms[1][0];
     231            Rc = (Rc + P1)*gamma + P0;
     232            dR = (dR + P1);
     233
     234            double gamma_new = gamma - (Rc - Ro) / dR;
     235            gamma = gamma_new;
     236          }
     237         
     238          double theta = 90.0 - gamma * DEG_RAD ;
     239          ctht = cos (RAD_DEG * theta);
     240          stht = sin (RAD_DEG * theta);
     241          break;
     242        }
     243
    222244      default:
    223245        return (FALSE);
    … …  
    354376
    355377      case PROJ_ZPN:
     378        // the forward projection is:
     379        // theta = atan2(stht, ctht)
     380        // gamma = (pi/2 - theta) : theta in radians
     381        // Ro = sum (P_i gamma^i)
     382        // R  = (180/pi) Ro
     383
    356384        // Ro = (pi/180)(90 - theta)
    357385        // R = (180/pi)sum (P_i R^i)
     386
     387        // is ZPN defined for Npolyterms = 0 or 1?
     388
    358389        ctht = hypot(sphi, cphi);
    359390        theta = atan2 (stht, ctht);
    360391
    361         double Rc;
    362         if (UKIRT_ONLY) {
    363           double P1 = coords[0].polyterms[0][1];
    364           double P3 = coords[0].polyterms[0][3];
    365           double gamma = RAD_DEG * (90 - DEG_RAD * theta);
    366           Rc = P1*gamma + P3*gamma*gamma*gamma;
    367         } else {
    368           double Ro = RAD_DEG * (90 - DEG_RAD * theta);
    369 
    370           // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
    371           i = coords[0].Npolyterms - 1;
    372           Rc = 0.0;
    373           while (i > 0) {
    374             if (i < 7) {
    375               Rc = (Rc + coords[0].polyterms[0][i])*Ro;
    376             } else {
    377               Rc = (Rc + coords[0].polyterms[1][i-7])*Ro;
    378             }
    379             i --;
    380           }
    381           Rc += coords[0].polyterms[0][i];
     392        double Ro;
     393        double gamma = M_PI_2 - theta;
     394
     395        // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
     396        Ro = 0.0;
     397        for (i = coords[0].Npolyterms - 1; i > 0; i --) {
     398          double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
     399          Ro = (Ro + Pi)*gamma;
    382400        }
    383         Rc = DEG_RAD * Rc;
     401        Ro += coords[0].polyterms[0][0];
     402
     403        Rc = DEG_RAD * Ro;
    384404
    385405        *L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ;
    … …  
    609629enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN};
    610630
     631int GetRadialZPN (Coords *coords, Header *header);
     632
    611633int GetCoords (Coords *coords, Header *header) {
    612634 
    613   int i, status, status1, status2, itmp, Polynomial, Polyterm;
     635  int status, status1, status2, itmp, Polynomial, Polyterm;
    614636  double Lambda, rotate, rotate1, rotate2, scale;
    615637  double equinox;
    … …  
    666688      status &= gfits_scan (header, "PC002002", "%f",  1, &coords[0].pc2_2);
    667689
     690      ctype = &coords[0].ctype[4];
     691
     692      // read the ZPN coeffients PV2_i (i = 0 < 14)
     693      // ZPN is inconsistent with the other Polynomial types
     694      if (!strcmp (ctype, "-ZPN")) {
     695        GetRadialZPN (coords, header);
     696        break;
     697      }
     698
    668699      /* set NPLYTERM based on header.  if NPLYTERM is missing, it should have a
    669700         value of 0, unless the projection type is one of PLY, DIS, WRP, in which
    670701         case it should be set to 3 */
    671       ctype = &coords[0].ctype[4];
    672702      Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
    673703      Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp);
    … …  
    723753      coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
    724754      coords[0].pc2_2 =  cos(rotate*RAD_DEG);
     755
     756      // read the ZPN coeffients PV2_i (i = 0 < 14)
     757      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
    725758      break;
    726759
    … …  
    744777      coords[0].pc2_2 /= scale;
    745778
    746       if (!strcmp (&coords[0].ctype[4], "-ZPN")) {
    747         int found;
    748         for (i = 0; i < 14; i++) {
    749           char name[64];
    750           snprintf (name, 64, "PV2_%d", i);
    751           if (i < 7) {
    752             found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[0][i]);
    753           } else {
    754             found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[1][i-7]);
    755           }
    756           if ((i == 0) && !found) {
    757             coords[0].polyterms[0][0] = 0.0;
    758             continue;
    759           }
    760           if ((i == 1) && !found) {
    761             coords[0].polyterms[0][1] = 1.0;
    762             continue;
    763           }
    764           if (!found) {
    765             coords[0].Npolyterms = i;
    766             break;
    767           }
    768         }
    769       }
     779      // read the ZPN coeffients PV2_i (i = 0 < 14)
     780      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
    770781      break;
    771782
    … …  
    819830  }
    820831  return (status);
     832}
     833
     834int GetRadialZPN (Coords *coords, Header *header) {
     835
     836  // RA---ZPN can have up to 14 radial polynomial terms.  these are stored in
     837  // polyterms[0][0] - [6][0] for the first 7 and [0][1] - [6][1] for the rest
     838  // these terms are coeffients of a polynomial of the radial distances
     839
     840  // read the ZPN coeffients PV2_i (i = 0 < 14)
     841  int found;
     842  int Nmax = 0;
     843  for (int i = 0; i < 14; i++) {
     844    char name[64];
     845    snprintf (name, 64, "PV2_%d", i);
     846    if (i < 7) {
     847      coords[0].polyterms[i][0] = 0.0;
     848      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i][0]);
     849    } else {
     850      coords[0].polyterms[i-7][1] = 0.0;
     851      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i-7][1]);
     852    }
     853    // PV2_1 is implicit if not present
     854    if ((i == 1) && !found) {
     855      coords[0].polyterms[1][0] = 1.0;
     856      continue;
     857    }
     858    // set Npolyterms based on the largest coefficient found
     859    if (found) {
     860      Nmax = i;
     861    }
     862  }
     863  coords[0].Npolyterms = Nmax + 1;
     864  return TRUE;
    821865}
    822866

  • trunk/Ohana/src/libdvo/src/dbExtractMeasures.c

    r39670 r39926  
    427427      break;
    428428    case MEAS_TMEAN: /* OK */
     429      value.Flt = TimeValue (average[0].Tmean, TimeReference, TimeFormat);
     430      break;
     431    case MEAS_TRANGE: /* OK */
    429432      value.Flt = GetTimeRange (average[0].Trange, TimeFormat);
    430       break;
    431     case MEAS_TRANGE: /* OK */
    432       value.Flt = TimeValue (average[0].Trange, 0, TimeFormat);
    433433      break;
    434434    case MEAS_NMEAS: /* OK */
    … …  
    461461      break;
    462462    case MEAS_RA_FIT_OFFSET: /* OK */
     463      // RA_epoch_fit = RA_mean + uR*(t - Tmean)/cos(dec) + plx*parR
     464      // note that this extraction ignores parallax
    463465      dT = (measure[0].t - average[0].Tmean) / (86400*365.25);
    464       dR = dvoOffsetR (measure, average);
    465       value.Flt = average[0].uR * dT + dR;
     466      dR = dvoOffsetR (measure, average); // RA_epoch - RA_mean (** NOT local linear distance **)
     467      value.Flt = dR*cos(RAD_DEG*measure[0].D) - average[0].uR * dT;
     468      // this is the local linear distance of the measurement from the fit
    466469      break;
    467470    case MEAS_DEC_FIT_OFFSET: /* OK */
    468471      dT = (measure[0].t - average[0].Tmean) / (86400*365.25);
    469472      dD = dvoOffsetD (measure, average);
    470       value.Flt = average[0].uD * dT + dD;
     473      value.Flt = dD - average[0].uD * dT;
    471474      break;
    472475    case MEAS_RA_OFFSET_ERR: /* OK */

  • trunk/Ohana/src/libdvo/src/galaxy_model.c

    r39588 r39926  
    2828    V_sol  =  11.18; // km/sec
    2929    W_sol  =   7.61; // km/sec
     30    return TRUE;
     31  }
     32  if (!strcmp(version, "TEST-CONSTANT")) {
     33    // use for testing
     34    A_oort = +47.40; // km/sec/kpc
     35    B_oort = -47.40; // km/sec/kpc
     36    U_sol  =   0.00; // km/sec
     37    V_sol  =   0.00; // km/sec
     38    W_sol  =   0.00; // km/sec
     39    return TRUE;
     40  }
     41  if (!strcmp(version, "TEST-ZERO")) {
     42    // use for testing
     43    A_oort =   0.00; // km/sec/kpc
     44    B_oort =   0.00; // km/sec/kpc
     45    U_sol  =   0.00; // km/sec
     46    V_sol  =   0.00; // km/sec
     47    W_sol  =   0.00; // km/sec
    3048    return TRUE;
    3149  }
Note: See TracChangeset for help on using the changeset viewer.