IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39688


Ignore:
Timestamp:
Sep 11, 2016, 2:36:01 PM (10 years ago)
Author:
eugene
Message:

make the ZPN projection general, not UKIRT specific

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20160809/Ohana/src/libdvo/src/coordops.c

    r39608 r39688  
    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;
    … …  
    157158          ctht = 0.0;
    158159        } else {
    159           T = DEG_RAD / R;
    160           stht =   T / sqrt ( 1.0 + T*T);
    161           ctht = 1.0 / sqrt ( 1.0 + T*T);
     160          // T = DEG_RAD / R; // T in 1/radians
     161          // stht =   T / sqrt ( 1.0 + T*T);
     162          // ctht = 1.0 / sqrt ( 1.0 + T*T);
     163
     164          T = RAD_DEG * R;
     165          stht = 1.0 / sqrt ( 1.0 + T*T);
     166          ctht =   T / sqrt ( 1.0 + T*T);
    162167        }
     168        break;
     169      case PROJ_SIN:
     170        // R = (180/pi) cos (theta)
     171        ctht = RAD_DEG * R;
     172        stht = sqrt (1 - ctht*ctht);
    163173        break;
    164174      case PROJ_STG:
    … …  
    167177        ctht = sqrt (1 - stht*stht);
    168178        break;
    169       case PROJ_SIN:
    170         // R = (180/pi) cos (theta)
    171         ctht = RAD_DEG * R;
    172         stht = sqrt (1 - ctht*ctht);
    173         break;
    174179      case PROJ_ARC:
    175180        // R = 90 - theta (degrees)
    … …  
    177182        stht = cos (RAD_DEG * R);
    178183        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         }
    211184
    212185      case PROJ_ZEA:
    … …  
    220193        ctht = sqrt (1 - stht*stht);
    221194        break;
     195
     196      case PROJ_ZPN:
     197
     198        // the forward projection is:
     199        // theta = atan2(stht, ctht)
     200        // gamma = (pi/2 - theta) : theta in radians
     201        // Ro = sum (P_i gamma^i)
     202        // R  = (180/pi) Ro
     203
     204        // given R, we need to find theta:
     205        // Ro = R * (pi / 180) = sum (P_i gamma^i)
     206        // solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson
     207
     208        // use Ro to get a guess for gamma and iterate
     209
     210        {
     211          double Ro = RAD_DEG * R;
     212         
     213          // find the roots of f(gamma) - Ro = 0
     214          // starting guess for gamma is (Ro - P0) / P1
     215          double gamma = (Ro - coords[0].polyterms[0][0]) / coords[0].polyterms[1][0];
     216         
     217          int iter;
     218          for (iter = 0; iter < 5; iter++) {
     219           
     220            double Rc = 0.0; // this will hold the ander
     221            double dR = 0.0;
     222            for (int i = coords[0].Npolyterms - 1; i > 1; i--) {
     223              double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
     224              Rc = (Rc + Pi)*gamma;
     225              dR = (dR + i*Pi)*gamma;
     226            }
     227            double P0 = coords[0].polyterms[0][0];
     228            double P1 = coords[0].polyterms[1][0];
     229            Rc = (Rc + P1)*gamma + P0;
     230            dR = (dR + P1);
     231
     232            double gamma_new = gamma - (Rc - Ro) / dR;
     233            gamma = gamma_new;
     234          }
     235         
     236          double theta = 90.0 - gamma * DEG_RAD ;
     237          ctht = cos (RAD_DEG * theta);
     238          stht = sin (RAD_DEG * theta);
     239          break;
     240        }
     241
    222242      default:
    223243        return (FALSE);
    … …  
    354374
    355375      case PROJ_ZPN:
     376        // the forward projection is:
     377        // theta = atan2(stht, ctht)
     378        // gamma = (pi/2 - theta) : theta in radians
     379        // Ro = sum (P_i gamma^i)
     380        // R  = (180/pi) Ro
     381
    356382        // Ro = (pi/180)(90 - theta)
    357383        // R = (180/pi)sum (P_i R^i)
     384
     385        // is ZPN defined for Npolyterms = 0 or 1?
     386
    358387        ctht = hypot(sphi, cphi);
    359388        theta = atan2 (stht, ctht);
    360389
    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];
     390        double Ro;
     391        double gamma = M_PI_2 - theta;
     392
     393        // i = 0 .. Npolyterms - 1 (1 <= Npolyterms <= 21)
     394        Ro = 0.0;
     395        for (i = coords[0].Npolyterms - 1; i > 0; i --) {
     396          double Pi = (i < 7) ? coords[0].polyterms[i][0] : coords[0].polyterms[i-7][1];
     397          Ro = (Ro + Pi)*gamma;
    382398        }
    383         Rc = DEG_RAD * Rc;
     399        Ro += coords[0].polyterms[0][0];
     400
     401        Rc = DEG_RAD * Ro;
    384402
    385403        *L = (ctht == 0.0) ? 0.0 : +Rc * sphi / ctht ;
    … …  
    609627enum {COORD_TYPE_NONE, COORD_TYPE_PC, COORD_TYPE_ROT, COORD_TYPE_CD, COORD_TYPE_LIN};
    610628
     629int GetRadialZPN (Coords *coords, Header *header);
     630
    611631int GetCoords (Coords *coords, Header *header) {
    612632 
    613   int i, status, status1, status2, itmp, Polynomial, Polyterm;
     633  int status, status1, status2, itmp, Polynomial, Polyterm;
    614634  double Lambda, rotate, rotate1, rotate2, scale;
    615635  double equinox;
    … …  
    666686      status &= gfits_scan (header, "PC002002", "%f",  1, &coords[0].pc2_2);
    667687
     688      ctype = &coords[0].ctype[4];
     689
     690      // read the ZPN coeffients PV2_i (i = 0 < 14)
     691      // ZPN is inconsistent with the other Polynomial types
     692      if (!strcmp (ctype, "-ZPN")) {
     693        GetRadialZPN (coords, header);
     694        break;
     695      }
     696
    668697      /* set NPLYTERM based on header.  if NPLYTERM is missing, it should have a
    669698         value of 0, unless the projection type is one of PLY, DIS, WRP, in which
    670699         case it should be set to 3 */
    671       ctype = &coords[0].ctype[4];
    672700      Polynomial = !strcmp (ctype, "-PLY") || !strcmp (ctype, "-DIS") || !strcmp (ctype, "-WRP");
    673701      Polyterm = gfits_scan (header, "NPLYTERM", "%d", 1, &itmp);
    … …  
    723751      coords[0].pc2_1 =  sin(rotate*RAD_DEG) / Lambda;
    724752      coords[0].pc2_2 =  cos(rotate*RAD_DEG);
     753
     754      // read the ZPN coeffients PV2_i (i = 0 < 14)
     755      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
    725756      break;
    726757
    … …  
    744775      coords[0].pc2_2 /= scale;
    745776
    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       }
     777      // read the ZPN coeffients PV2_i (i = 0 < 14)
     778      if (!strcmp (&coords[0].ctype[4], "-ZPN")) GetRadialZPN (coords, header);
    770779      break;
    771780
    … …  
    819828  }
    820829  return (status);
     830}
     831
     832int GetRadialZPN (Coords *coords, Header *header) {
     833
     834  // RA---ZPN can have up to 14 radial polynomial terms.  these are stored in
     835  // polyterms[0][0] - [6][0] for the first 7 and [0][1] - [6][1] for the rest
     836  // these terms are coeffients of a polynomial of the radial distances
     837
     838  // read the ZPN coeffients PV2_i (i = 0 < 14)
     839  int found;
     840  int Nmax = 0;
     841  for (int i = 0; i < 14; i++) {
     842    char name[64];
     843    snprintf (name, 64, "PV2_%d", i);
     844    if (i < 7) {
     845      coords[0].polyterms[i][0] = 0.0;
     846      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i][0]);
     847    } else {
     848      coords[0].polyterms[i-7][1] = 0.0;
     849      found = gfits_scan (header, name, "%f", 1, &coords[0].polyterms[i-7][1]);
     850    }
     851    // PV2_1 is implicit if not present
     852    if ((i == 1) && !found) {
     853      coords[0].polyterms[1][0] = 1.0;
     854      continue;
     855    }
     856    // set Npolyterms based on the largest coefficient found
     857    if (found) {
     858      Nmax = i;
     859    }
     860  }
     861  coords[0].Npolyterms = Nmax + 1;
     862  return TRUE;
    821863}
    822864
Note: See TracChangeset for help on using the changeset viewer.