IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

add the ZPN projection used by HypersuprimeCam and UKIRT : this projection includes a radial profile which needs to be added to the structure

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20160809/psLib/src/astro/psCoord.c

    r19504 r39684  
    234234    const psPlane* coords)
    235235{
     236  if (!transform) {
     237    fprintf (stderr, "problem\n");
     238  }
    236239    PS_ASSERT_PTR_NON_NULL(transform, NULL);
    237240    PS_ASSERT_PTR_NON_NULL(transform->x, NULL);
    … …  
    305308void projectionFree(psProjection *p)
    306309{
    307     // There are no dynamically allocated items
     310  if (!p) return;
     311  if (!p->radial) return;
     312  psFree (p->radial);
    308313}
    309314
    … …  
    321326    p->Ys = Ys;
    322327    p->type = type;
     328    p->radial = NULL;
    323329
    324330    psMemSetDeallocator(p, (psFreeFunc) projectionFree);
    … …  
    368374      case PS_PROJ_ZEA:
    369375      case PS_PROJ_ZPL:
     376      case PS_PROJ_ZPN:
    370377        class = PS_PROJECTION_CLASS_ZENITHAL;
    371378        break;
    … …  
    442449            out->y = -zeta * cosPhiCT;
    443450            break;
     451
     452          case PS_PROJ_ZPN: {
     453            // the forward projection is:
     454            // theta = atan2(stht, ctht)
     455            // gamma = (pi/2 - theta) : theta in radians
     456            // Ro = sum (P_i gamma^i)
     457            // R  = (180/pi) Ro
     458           
     459            // Ro = (pi/180)(90 - theta)
     460            // R = (180/pi)sum (P_i R^i)
     461           
     462            // is ZPN defined for Npolyterms = 0 or 1?
     463
     464            double cosTheta = hypot(sinPhiCT, cosPhiCT);
     465            double theta = atan2 (sinTheta, cosTheta);
     466           
     467            double gamma = M_PI_2 - theta;
     468           
     469            // i = 0 .. N - 1 (1 <= Npolyterms <= 21)
     470            double Ro = 0.0;
     471            for (int i = projection->radial->n - 1; i > 0; i--) {
     472              double Pi = projection->radial->data.F64[i];
     473              Ro = (Ro + Pi)*gamma;
     474            }
     475            Ro += projection->radial->data.F64[0];
     476           
     477            out->x = (cosTheta == 0.0) ? 0.0 : +Ro * sinPhiCT / cosTheta ;
     478            out->y = (cosTheta == 0.0) ? 0.0 : -Ro * cosPhiCT / cosTheta ;
     479            break;
     480          }
     481
    444482          default:
    445483            psAbort("invalid projection");
    … …  
    525563      case PS_PROJ_ZEA:
    526564      case PS_PROJ_ZPL:
     565      case PS_PROJ_ZPN:
    527566        class = PS_PROJECTION_CLASS_ZENITHAL;
    528567        break;
    … …  
    574613          case PS_PROJ_SIN:
    575614            // Orhtographic deprojection
     615            cosTheta = R;
    576616            sinTheta = sqrt (1 - R*R);
    577             cosTheta = R;
    578617            break;
    579618          case PS_PROJ_STG:
    580             psAbort("STG not defined");
     619            sinTheta = (4 - R) / (4 + R);
     620            cosTheta = sqrt (1 - sinTheta*sinTheta);
    581621            break;
    582622          case PS_PROJ_ZEA:
    … …  
    587627            cosTheta = sqrt (1 - PS_SQR(sinTheta));
    588628            break;
     629
     630      case PS_PROJ_ZPN:
     631
     632        // the forward projection is:
     633        // theta = atan2(stht, ctht)
     634        // gamma = (pi/2 - theta) : theta in radians
     635        // Ro = sum (P_i gamma^i)
     636        // R  = (180/pi) Ro
     637
     638        // given R, we need to find theta:
     639        // Ro = R * (pi / 180) = sum (P_i gamma^i)
     640        // solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson
     641
     642        // use Ro to get a guess for gamma and iterate
     643
     644        {
     645          // find the roots of f(gamma) - R = 0
     646          // starting guess for gamma is (R - P0) / P1
     647          double gamma = (R - projection->radial->data.F64[0]) / projection->radial->data.F64[1];
     648         
     649          for (int iter = 0; iter < 5; iter++) {
     650           
     651            double Rc = 0.0; // this will hold the ander
     652            double dR = 0.0;
     653            for (int i = projection->radial->n - 1; i > 1; i--) {
     654              double Pi = projection->radial->data.F64[i];
     655              Rc = (Rc + Pi)*gamma;
     656              dR = (dR + i*Pi)*gamma;
     657            }
     658            double P0 = projection->radial->data.F64[0];
     659            double P1 = projection->radial->data.F64[1];
     660            Rc = (Rc + P1)*gamma + P0;
     661            dR = (dR + P1);
     662
     663            double gamma_new = gamma - (Rc - R) / dR;
     664            gamma = gamma_new;
     665          }
     666         
     667          double theta = M_PI_2 - gamma ;
     668          cosTheta = cos (theta);
     669          sinTheta = sin (theta);
     670          break;
     671        }
     672       
    589673          default:
    590674            psAbort("invalid projection");
    … …  
    11831267
    11841268    switch (type) {
    1185       case PS_PROJ_LIN:
    1186         psStringAppend (&name, "%s-LIN", prefix);
    1187         return name;
    1188       case PS_PROJ_PLY:
    1189         psStringAppend (&name, "%s-PLY", prefix);
    1190         return name;
    1191       case PS_PROJ_WRP:
    1192         psStringAppend (&name, "%s-WRP", prefix);
    1193         return name;
    1194       case PS_PROJ_TAN:
    1195         psStringAppend (&name, "%s-TAN", prefix);
    1196         return name;
    1197       case PS_PROJ_TNX:
    1198         psStringAppend (&name, "%s-TNX", prefix);
    1199         return name;
    1200       case PS_PROJ_DIS:
    1201         psStringAppend (&name, "%s-DIS", prefix);
    1202         return name;
    1203       case PS_PROJ_SIN:
    1204         psStringAppend (&name, "%s-SIN", prefix);
    1205         return name;
    1206       case PS_PROJ_STG:
    1207         psStringAppend (&name, "%s-STG", prefix);
    1208         return name;
    1209       case PS_PROJ_AIT:
    1210         psStringAppend (&name, "%s-AIT", prefix);
    1211         return name;
    1212       case PS_PROJ_PAR:
    1213         psStringAppend (&name, "%s-PAR", prefix);
    1214         return name;
    1215       case PS_PROJ_GLS:
    1216         psStringAppend (&name, "%s-GLS", prefix);
    1217         return name;
    1218       case PS_PROJ_CAR:
    1219         psStringAppend (&name, "%s-CAR", prefix);
    1220         return name;
    1221       case PS_PROJ_MER:
    1222         psStringAppend (&name, "%s-MER", prefix);
    1223         return name;
     1269      case PS_PROJ_LIN: psStringAppend (&name, "%s-LIN", prefix); return name;
     1270      case PS_PROJ_PLY: psStringAppend (&name, "%s-PLY", prefix); return name;
     1271      case PS_PROJ_WRP: psStringAppend (&name, "%s-WRP", prefix); return name;
     1272      case PS_PROJ_TAN: psStringAppend (&name, "%s-TAN", prefix); return name;
     1273
     1274      case PS_PROJ_DIS: psStringAppend (&name, "%s-DIS", prefix); return name;
     1275      case PS_PROJ_SIN: psStringAppend (&name, "%s-SIN", prefix); return name;
     1276      case PS_PROJ_STG: psStringAppend (&name, "%s-STG", prefix); return name;
     1277      case PS_PROJ_TNX: psStringAppend (&name, "%s-TNX", prefix); return name;
     1278
     1279      case PS_PROJ_ZEA: psStringAppend (&name, "%s-ZEA", prefix); return name;
     1280      case PS_PROJ_ZPL: psStringAppend (&name, "%s-ZPL", prefix); return name;
     1281      case PS_PROJ_ZPN: psStringAppend (&name, "%s-ZPN", prefix); return name;
     1282      case PS_PROJ_AIT: psStringAppend (&name, "%s-AIT", prefix); return name;
     1283
     1284      case PS_PROJ_PAR: psStringAppend (&name, "%s-PAR", prefix); return name;
     1285      case PS_PROJ_GLS: psStringAppend (&name, "%s-GLS", prefix); return name;
     1286      case PS_PROJ_CAR: psStringAppend (&name, "%s-CAR", prefix); return name;
     1287      case PS_PROJ_MER: psStringAppend (&name, "%s-MER", prefix); return name;
    12241288
    12251289      default:
    … …  
    12381302    if (!strcmp (&name[4], "-WRP")) return PS_PROJ_WRP;
    12391303    if (!strcmp (&name[4], "-TAN")) return PS_PROJ_TAN;
    1240     if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
    12411304    if (!strcmp (&name[4], "-DIS")) return PS_PROJ_DIS;
    12421305    if (!strcmp (&name[4], "-SIN")) return PS_PROJ_SIN;
    12431306    if (!strcmp (&name[4], "-STG")) return PS_PROJ_STG;
     1307    if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
     1308    if (!strcmp (&name[4], "-ZEA")) return PS_PROJ_ZEA;
     1309    if (!strcmp (&name[4], "-ZPL")) return PS_PROJ_ZPL;
     1310    if (!strcmp (&name[4], "-ZPN")) return PS_PROJ_ZPN;
    12441311    if (!strcmp (&name[4], "-AIT")) return PS_PROJ_AIT;
    12451312    if (!strcmp (&name[4], "-PAR")) return PS_PROJ_PAR;
Note: See TracChangeset for help on using the changeset viewer.