IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3296


Ignore:
Timestamp:
Feb 19, 2005, 4:20:04 PM (21 years ago)
Author:
eugene
Message:

added WRP and DIS models, adjusted PLY terms

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2413 r3296  
    44int XY_to_RD (double *ra, double *dec, double x, double y, Coords *coords) {
    55
     6  int Polynomial, Zenithal, Cartesian, PseudoCyl;
     7  char *type;
    68  double L, M, X, Y, T, Z;
    79  double R, sphi, cphi, stht, ctht;
    … …  
    1113  *dec = 0;
    1214  stht = ctht = 1;
     15  type = &coords[0].ctype[4];
     16 
     17  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
     18  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
     19  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
     20  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
     21  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
    1322
    1423  /** convert pixel coordinates to cartesian system **/
    1524  X = coords[0].cdelt1*(x - coords[0].crpix1);
    1625  Y = coords[0].cdelt2*(y - coords[0].crpix2);
    17   if (!strcmp(&coords[0].ctype[4], "-PLY")) {
    18     if (coords[0].Npolyterms > 2) {
    19       X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
    20       Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]);
    21     }
    22     if (coords[0].Npolyterms > 2) {
    23       X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]);
    24       Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]);
    25     }
    26   }
    2726
    2827  L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
    2928  M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
    30   /* in FITS ref, L,M = x, y  alpha, delta = phi, theta */
    31 
     29
     30  /** extra polynomial terms **/
     31  if (Polynomial) {
     32    if (coords[0].Npolyterms > 1) {
     33      L += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
     34      M += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
     35    }
     36    if (coords[0].Npolyterms > 2) {
     37      L += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
     38      M += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
     39    }
     40  }
     41
     42  /**** Locally Cartesian Projections ****/
     43  if (Cartesian) {
     44    *ra  = L + coords[0].crval1;
     45    *dec = M + coords[0].crval2;
     46    return (TRUE);
     47  }
     48 
    3249  /**** Zenithal Projections ****/
    33   if (!strcmp(&coords[0].ctype[4], "-PLY") ||
    34       !strcmp(&coords[0].ctype[4], "-TAN") ||
    35       !strcmp(&coords[0].ctype[4], "-SIN") ||
    36       !strcmp(&coords[0].ctype[4], "-ZEA") ||
    37       !strcmp(&coords[0].ctype[0], "MM")) {
     50  if (Zenithal) {
    3851    R = hypot (L,M);
    3952    if ((L == 0) && (M == 0)) {
    4053      sphi = 0;
    4154      cphi = 1;
    42     }
    43     else {
     55    } else {
    4456      sphi =  L / R;
    4557      cphi = -M / R;
    4658    }
    4759
    48     if (!strcmp(&coords[0].ctype[4], "-PLY") ||
    49         !strcmp(&coords[0].ctype[4], "-TAN")) {
     60    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
    5061      if (R == 0) {
    5162        stht = 1.0;
    5263        ctht = 0.0;
    53       }
    54       else {
     64      } else {
    5565        T = DEG_RAD / R;
    5666        stht =   T / sqrt ( 1.0 + T*T);
    … …  
    5868      }
    5969    }
    60     if (!strcmp(&coords[0].ctype[4], "-SIN") ||
    61         !strcmp(&coords[0].ctype[0], "MM")) {
     70    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
    6271      ctht = RAD_DEG * R;
    6372      stht = sqrt (1 - ctht*ctht);
    6473    }
    65     if (!strcmp(&coords[0].ctype[4], "-ZEA")) {
     74    if (!strcmp(type, "-ZEA")) {
    6675      stht = 1 - 0.5*SQ(R*RAD_DEG);
    6776      ctht = sqrt (1 - stht*stht);
    … …  
    7988    *ra  = DEG_RAD*alpha + coords[0].crval1;
    8089    *dec = DEG_RAD*delta;
    81   }
    82  
    83   /**** Locally Cartesian Projections ****/
    84   if (!strcmp(&coords[0].ctype[4], "-LIN") ||
    85       !strcmp(&coords[0].ctype[0], "GENE")) {
    86     *ra  = L + coords[0].crval1;
    87     *dec = M + coords[0].crval2;
     90    return (TRUE);
    8891  }
    8992 
    9093  /**** Other Conventional Projections ****/
    91   if (!strcmp(&coords[0].ctype[4], "-AIT")) {
    92     Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
    93     alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 1.0);
    94     delta = DEG_RAD * asin (RAD_DEG*M*Z);
    95     *ra  = alpha + coords[0].crval1;
    96     *dec = delta + coords[0].crval2;
    97   }
    98   if (!strcmp(&coords[0].ctype[4], "-GLS")) {
    99     /* L,M in degrees, alpha,delta in degrees */
    100     alpha = L / cos (RAD_DEG * M);
    101     delta = M;
    102     *ra  = alpha + coords[0].crval1;
    103     *dec = delta + coords[0].crval2;
    104   }
    105   if (!strcmp(&coords[0].ctype[4], "-PAR")) {
    106     /* L,M in degrees, alpha,delta in degrees */
    107     alpha = L / (1.0 - SQ(2.0*M/180));
    108     delta = 3 * DEG_RAD * asin (M/180.0);
    109     *ra  = alpha + coords[0].crval1;
    110     *dec = delta + coords[0].crval2;
    111   }
    112   return (TRUE);
     94  if (PseudoCyl) {
     95    if (!strcmp(type, "-AIT")) {
     96      Z = sqrt (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
     97      alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*SQ(Z) - 1.0);
     98      delta = DEG_RAD * asin (RAD_DEG*M*Z);
     99      *ra  = alpha + coords[0].crval1;
     100      *dec = delta + coords[0].crval2;
     101    }
     102    if (!strcmp(type, "-GLS")) {
     103      /* L,M in degrees, alpha,delta in degrees */
     104      alpha = L / cos (RAD_DEG * M);
     105      delta = M;
     106      *ra  = alpha + coords[0].crval1;
     107      *dec = delta + coords[0].crval2;
     108    }
     109    if (!strcmp(type, "-PAR")) {
     110      /* L,M in degrees, alpha,delta in degrees */
     111      alpha = L / (1.0 - SQ(2.0*M/180));
     112      delta = 3 * DEG_RAD * asin (M/180.0);
     113      *ra  = alpha + coords[0].crval1;
     114      *dec = delta + coords[0].crval2;
     115    }
     116    return (TRUE);
     117  }
     118  return (FALSE);
    113119}
    114120
    115121int RD_to_XY (double *x, double *y, double ra, double dec, Coords *coords) {
    116122
     123  char *type;
     124  int i, status, Polynomial, Zenithal, Cartesian, PseudoCyl;
    117125  double phi, theta;
    118   double tmp_d;
    119   double X, Y, sphi, cphi, stht;
     126  double determ;
     127  double X, Y, L, M, Lo, Mo, dL, dM;
     128  double sphi, cphi, stht;
    120129  double salp, calp, sdel, cdel, sdp, cdp;
    121130  double P, A, Rc;
    122   int status;
    123 
    124   if (!strcmp(&coords[0].ctype[4], "-PLY")) {
    125     /* fprintf (stderr, "approximate to polynomial TAN plane fit\n");  */
    126   }
    127 
    128   X = Y = 1;
     131
    129132  status = TRUE;
    130133  *x = 0;
    131134  *y = 0;
     135  type = &coords[0].ctype[4];
     136  L = M = 0;
     137
     138  Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP");
     139  Zenithal   = !strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-SIN") || !strcmp(type, "-ZEA") || !strcmp(&coords[0].ctype[0], "MM");
     140  Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
     141  PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
     142  if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
    132143
    133144  /**** Locally Cartesian Projections ****/
    134   if (!strcmp(&coords[0].ctype[0], "GENE") || !strcmp(&coords[0].ctype[4], "-LIN")) {
    135     X = (ra  - coords[0].crval1);
    136     Y = (dec - coords[0].crval2);
     145  if (Cartesian) {
     146    L = (ra  - coords[0].crval1);
     147    M = (dec - coords[0].crval2);
    137148  }
    138149 
    139150  /**** Zenithal Projections ****/
    140   if (!strcmp(&coords[0].ctype[4], "-PLY") ||
    141       !strcmp(&coords[0].ctype[4], "-TAN") ||
    142       !strcmp(&coords[0].ctype[4], "-ZEA") ||
    143       !strcmp(&coords[0].ctype[4], "-SIN") ||
    144       !strcmp(&coords[0].ctype[0], "MM")) {
     151  if (Zenithal)  {
    145152    sdp  = sin(RAD_DEG*coords[0].crval2);
    146153    cdp  = cos(RAD_DEG*coords[0].crval2);
    … …  
    155162    if (stht < 0) status = FALSE;
    156163
    157     if (!strcmp(&coords[0].ctype[4], "-PLY") || !strcmp(&coords[0].ctype[4], "-TAN") ) {
    158       X =  DEG_RAD * sphi / stht;
    159       Y = -DEG_RAD * cphi / stht;
    160     }
    161     if (!strcmp(&coords[0].ctype[4], "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
    162       X =  DEG_RAD * sphi;
    163       Y = -DEG_RAD * cphi;
    164     }
    165     if (!strcmp(&coords[0].ctype[4], "-ZEA")) {
     164    if (!strcmp(type, "-PLY") || !strcmp(type, "-TAN") || !strcmp(type, "-DIS")) {
     165      L =  DEG_RAD * sphi / stht;
     166      M = -DEG_RAD * cphi / stht;
     167    }
     168    if (!strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM")) {
     169      L =  DEG_RAD * sphi;
     170      M = -DEG_RAD * cphi;
     171    }
     172    if (!strcmp(type, "-ZEA")) {
    166173      Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
    167       X =  Rc * sphi;
    168       Y = -Rc * cphi;
     174      L =  Rc * sphi;
     175      M = -Rc * cphi;
    169176      status = TRUE;
    170177    }
    … …  
    172179
    173180  /**** Other Standard Projections ****/
    174   if (!strcmp(&coords[0].ctype[4], "-AIT")) {
    175     phi = RAD_DEG*(ra - coords[0].crval1);
    176     theta = RAD_DEG*(dec - coords[0].crval2);
    177     P = 1.0 + cos (theta) * cos (0.5*phi);
    178     if (P != 0.0) {
    179       A =  DEG_RAD * sqrt (2.0 / P);
    180       X =  2.0 * A * cos (theta) * sin (0.5*phi);
    181       Y =  A * sin (theta);
    182     } else {
    183       X =  0.0;
    184       Y =  0.0;
    185     }   
    186   }
    187   if (!strcmp(&coords[0].ctype[4], "-GLS")) {
    188     phi = ra - coords[0].crval1;
    189     theta = dec - coords[0].crval2;
    190     X = phi * cos(RAD_DEG * theta);
    191     Y = theta;
    192   }
    193   if (!strcmp(&coords[0].ctype[4], "-PAR")) {
    194     phi = ra - coords[0].crval1;
    195     theta = dec - coords[0].crval2;
    196     X = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
    197     Y = 180.0 * sin (RAD_DEG*theta/3.0);
    198   }
    199    
    200   tmp_d = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
    201   *x = tmp_d * (coords[0].pc2_2*X - coords[0].pc1_2*Y) / coords[0].cdelt1 + coords[0].crpix1;
    202   *y = tmp_d * (coords[0].pc1_1*Y - coords[0].pc2_1*X) / coords[0].cdelt2 + coords[0].crpix2;
     181  if (PseudoCyl) {
     182    if (!strcmp(type, "-AIT")) {
     183      phi = RAD_DEG*(ra - coords[0].crval1);
     184      theta = RAD_DEG*(dec - coords[0].crval2);
     185      P = 1.0 + cos (theta) * cos (0.5*phi);
     186      if (P != 0.0) {
     187        A =  DEG_RAD * sqrt (2.0 / P);
     188        L =  2.0 * A * cos (theta) * sin (0.5*phi);
     189        M =  A * sin (theta);
     190      } else {
     191        L =  0.0;
     192        M =  0.0;
     193      }
     194    }
     195    if (!strcmp(type, "-GLS")) {
     196      phi = ra - coords[0].crval1;
     197      theta = dec - coords[0].crval2;
     198      L = phi * cos(RAD_DEG * theta);
     199      M = theta;
     200    }
     201    if (!strcmp(type, "-PAR")) {
     202      phi = ra - coords[0].crval1;
     203      theta = dec - coords[0].crval2;
     204      L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
     205      M = 180.0 * sin (RAD_DEG*theta/3.0);
     206    }
     207  }
     208
     209  /* convert L,M to X,Y */
     210  determ = 1.0 / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
     211  X = determ * (coords[0].pc2_2*L - coords[0].pc1_2*M);
     212  Y = determ * (coords[0].pc1_1*M - coords[0].pc2_1*L);
     213
     214  /** extra polynomial terms **/
     215  if (Polynomial) {
     216    for (i = 0; i < 3; i++) {
     217      Lo = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
     218      Mo = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
     219      if (coords[0].Npolyterms > 1) {
     220        Lo += X*X*coords[0].polyterms[0][0] + X*Y*coords[0].polyterms[1][0] + Y*Y*coords[0].polyterms[2][0];
     221        Mo += X*X*coords[0].polyterms[0][1] + X*Y*coords[0].polyterms[1][1] + Y*Y*coords[0].polyterms[2][1];
     222      }
     223      if (coords[0].Npolyterms > 2) {
     224        Lo += X*X*X*coords[0].polyterms[3][0] + X*X*Y*coords[0].polyterms[4][0] + X*Y*Y*coords[0].polyterms[5][0] + Y*Y*Y*coords[0].polyterms[6][0];
     225        Mo += X*X*X*coords[0].polyterms[3][1] + X*X*Y*coords[0].polyterms[4][1] + X*Y*Y*coords[0].polyterms[5][1] + Y*Y*Y*coords[0].polyterms[6][1];
     226      }
     227      dL = (L - Lo);
     228      dM = (M - Mo);
     229
     230      X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM);
     231      Y += determ * (coords[0].pc1_1*dM - coords[0].pc2_1*dL);
     232    }
     233  }
     234  /* check for correct size (iterate?) */
     235
     236  *x = X / coords[0].cdelt1 + coords[0].crpix1;
     237  *y = Y / coords[0].cdelt2 + coords[0].crpix2;
    203238
    204239  return (status);
    … …  
    391426   PCA2X0Y3 = coords.polyterm[6][1] = y^2                                               
    392427*/
     428
     429# if (0)
     430
     431  /** convert pixel coordinates to cartesian system **/
     432  X = coords[0].cdelt1*(x - coords[0].crpix1);
     433  Y = coords[0].cdelt2*(y - coords[0].crpix2);
     434  if (Polynomi) {
     435    if (coords[0].Npolyterms > 2) {
     436      X += coords[0].cdelt1*(x*x*coords[0].polyterms[0][0] + x*y*coords[0].polyterms[1][0] + y*y*coords[0].polyterms[2][0]);
     437      Y += coords[0].cdelt2*(x*x*coords[0].polyterms[0][1] + x*y*coords[0].polyterms[1][1] + y*y*coords[0].polyterms[2][1]);
     438    }
     439    if (coords[0].Npolyterms > 2) {
     440      X += coords[0].cdelt1*(x*x*x*coords[0].polyterms[3][0] + x*x*y*coords[0].polyterms[4][0] + x*y*y*coords[0].polyterms[5][0] + y*y*y*coords[0].polyterms[6][0]);
     441      Y += coords[0].cdelt2*(x*x*x*coords[0].polyterms[3][1] + x*x*y*coords[0].polyterms[4][1] + x*y*y*coords[0].polyterms[5][1] + y*y*y*coords[0].polyterms[6][1]);
     442    }
     443  }
     444
     445  L = (X*coords[0].pc1_1 + Y*coords[0].pc1_2);
     446  M = (X*coords[0].pc2_1 + Y*coords[0].pc2_2);
     447/** this code is the old method used for higher order terms.  they
     448    are essentially 6th order, with weird coupled terms.
     449    I don't think any real data used these terms, but they should
     450    be re-calculated, I would think
     451**/
     452
     453# endif
Note: See TracChangeset for help on using the changeset viewer.