IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 12069


Ignore:
Timestamp:
Feb 26, 2007, 5:14:14 PM (19 years ago)
Author:
eugene
Message:

defined projection and proj mode, cleaner use of cases

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/dvo-mods-2007-02/Ohana/src/libdvo/src/coordops.c

    r10510 r12069  
    1818  *dec = 0;
    1919  stht = ctht = 1;
    20   type = &coords[0].ctype[4];
    21  
    22   /* PLY is equiv to LIN with higher order terms
    23      ZPL is equiv to ZEA with higher order terms
    24      DIS is equiv to TAN with higher order terms
    25      WRP is equiv to PLY, with implied mosaic */
    26 
    27   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP") || !strcmp(type, "-ZPL");
    28   Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
    29   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
    30   Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
    31   Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
    32   Zenithal   = !strcmp(type, "-ZEA") || !strcmp(type, "-ZPL") || Zenith1 || Zenith2;
    33   if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
     20 
     21  proj = GetProjection (coords[0].ctype);
     22  mode = GetProjectionMode (proj);
     23  if (proj == PROJ_NONE) return (FALSE);
     24  if (proj == PROJ_MODE_NONE) return (FALSE);
    3425
    3526  /** convert pixel coordinates to cartesian system **/
    … …  
    5041  }
    5142
    52   /**** Locally Cartesian Projections ****/
    53   if (Cartesian) {
     43  /** Locally Cartesian Projections **/
     44  if (mode == PROJ_MODE_CARTESIAN) {
    5445    *ra  = L + coords[0].crval1;
    5546    *dec = M + coords[0].crval2;
    5647
    5748    /* mosaic astrometry : WRP is chip astrometry; apply mosaic (DIS) term */
    58     if (!strcmp(type, "-WRP")) {
     49    if (proj == PROJ_WRP) {
    5950      if (mosaic == NULL) return (FALSE);
    6051      XY_to_RD (ra, dec, L + coords[0].crval1, M + coords[0].crval2, mosaic);
    … …  
    6354  }
    6455 
    65   /**** Zenithal Projections ****/
    66   if (Zenithal) {
     56  /** Zenithal Projections **/
     57  if (mode == PROJ_MODE_ZENITHAL) {
    6758    R = hypot (L,M);
    6859    if ((L == 0) && (M == 0)) {
    … …  
    7465    }
    7566
    76     /* this is wrong : STG is not TAN - need to put in correct relationships.  but is a close approx */
    77     if (Zenith1) {
    78       if (R == 0) {
    79         stht = 1.0;
    80         ctht = 0.0;
    81       } else {
    82         T = DEG_RAD / R;
    83         stht =   T / sqrt ( 1.0 + T*T);
    84         ctht = 1.0 / sqrt ( 1.0 + T*T);
    85       }
    86     }
    87     if (Zenith2) {
    88       ctht = RAD_DEG * R;
    89       stht = sqrt (1 - ctht*ctht);
    90     }
    91     if (!strcmp(type, "-ZEA") || !strcmp(type, "-ZPL")) {
    92       if (R > 2*DEG_RAD) {
    93         *ra = L;
    94         *dec = M;
    95         return (FALSE);
    96       }
    97       stht = 1 - 0.5*SQ(R*RAD_DEG);
    98       ctht = sqrt (1 - stht*stht);
    99     }
    100 
     67    switch (proj) {
     68      case PROJ_TAN:
     69        if (R == 0) {
     70          stht = 1.0;
     71          ctht = 0.0;
     72        } else {
     73          T = DEG_RAD / R;
     74          stht =   T / sqrt ( 1.0 + T*T);
     75          ctht = 1.0 / sqrt ( 1.0 + T*T);
     76        }
     77        break;
     78      case PROJ_STG:
     79        stht = (4 - RAD_DEG*R) / (4 + RAD_DEC*R);
     80        ctht = sqrt (1 - stht*stht);
     81        break;
     82      case PROJ_SIN:
     83        ctht = RAD_DEG * R;
     84        stht = sqrt (1 - ctht*ctht);
     85        break;
     86      case PROJ_ZEA:
     87      case PROJ_ZPL:
     88        if (R > 2*DEG_RAD) {
     89          *ra = L;
     90          *dec = M;
     91          return (FALSE);
     92        }
     93        stht = 1 - 0.5*SQ(R*RAD_DEG);
     94        ctht = sqrt (1 - stht*stht);
     95        break;
     96    }
    10197    sdp  = sin(RAD_DEG*coords[0].crval2);
    10298    cdp  = cos(RAD_DEG*coords[0].crval2);
    … …  
    119115 
    120116  /**** Other Conventional Projections ****/
    121   if (PseudoCyl) {
    122     if (!strcmp(type, "-AIT")) {
     117  if (mode == PROJ_MODE_PSEUDOCLY) {
     118    switch (proj) {
     119      case PROJ_AIT:
    123120      Z2 = (1.0 - SQ(RAD_DEG*0.25*L) - SQ(RAD_DEG*0.5*M));
    124121      if (Z2 < 0) return (FALSE);
    … …  
    126123      alpha = 2.0 * DEG_RAD * atan2 (RAD_DEG*0.5*Z*L, 2.0*Z2 - 1.0);
    127124      delta = DEG_RAD * asin (RAD_DEG*M*Z);
    128       *ra  = alpha + coords[0].crval1;
    129       *dec = delta + coords[0].crval2;
    130     }
    131     if (!strcmp(type, "-GLS")) {
    132       /* L,M in degrees, alpha,delta in degrees */
    133       alpha = L / cos (RAD_DEG * M);
    134       delta = M;
    135       *ra  = alpha + coords[0].crval1;
    136       *dec = delta + coords[0].crval2;
    137     }
    138     if (!strcmp(type, "-PAR")) {
    139       /* L,M in degrees, alpha,delta in degrees */
    140       alpha = L / (1.0 - SQ(2.0*M/180));
    141       delta = 3 * DEG_RAD * asin (M/180.0);
    142       *ra  = alpha + coords[0].crval1;
    143       *dec = delta + coords[0].crval2;
    144     }
     125      break;
     126     
     127      case PROJ_GLS:
     128        /* L,M in degrees, alpha,delta in degrees */
     129        alpha = L / cos (RAD_DEG * M);
     130        delta = M;
     131        break;
     132      case PROJ_PAR:
     133        /* L,M in degrees, alpha,delta in degrees */
     134        alpha = L / (1.0 - SQ(2.0*M/180));
     135        delta = 3 * DEG_RAD * asin (M/180.0);
     136        break;
     137    }
     138    *ra  = alpha + coords[0].crval1;
     139    *dec = delta + coords[0].crval2;
    145140
    146141    /* rationalize ra range 0 - 360.0 */
    … …  
    170165  L = M = 0;
    171166
    172   Polynomial = !strcmp(type, "-PLY") || !strcmp(type, "-DIS") || !strcmp(type, "-WRP") || !strcmp(type, "-ZPL");
    173   Cartesian  = !strcmp(type, "-LIN") || !strcmp(type, "-PLY") || !strcmp(type, "-WRP") || !strcmp(&coords[0].ctype[0], "GENE");
    174   PseudoCyl  = !strcmp(type, "-AIT") || !strcmp(type, "-GLS") || !strcmp(type, "-PAR");
    175   Zenith1    = !strcmp(type, "-DIS") || !strcmp(type, "-TAN") || !strcmp(type, "-STG");
    176   Zenith2    = !strcmp(type, "-SIN") || !strcmp(&coords[0].ctype[0], "MM");
    177   Zenithal   = !strcmp(type, "-ZEA") || !strcmp(type, "-ZPL") || Zenith1 || Zenith2;
    178   if (!Zenithal && !Cartesian && !PseudoCyl) return (FALSE);
     167  proj = GetProjection (coords[0].ctype);
     168  mode = GetProjectionMode (proj);
     169  if (proj == PROJ_NONE) return (FALSE);
     170  if (proj == PROJ_MODE_NONE) return (FALSE);
    179171
    180172  /**** Locally Cartesian Projections ****/
    181   if (Cartesian) {
    182     if (!strcmp(type, "-WRP")) {
     173  if (mode == PROJ_MODE_CARTESIAN) {
     174    if (proj == PROJ_WRP) {
    183175      if (mosaic == NULL) return (FALSE);
    184176      RD_to_XY (&Lo, &Mo, ra, dec, mosaic);
    … …  
    192184 
    193185  /**** Zenithal Projections ****/
    194   if (Zenithal) {
     186  if (mode == PROJ_MODE_ZENITHAL) {
    195187    sdp  = sin(RAD_DEG*coords[0].crval2);
    196188    cdp  = cos(RAD_DEG*coords[0].crval2);
    … …  
    205197    if (stht < 0) status = FALSE;
    206198
    207     if (Zenith1) {
    208       L =  DEG_RAD * sphi / stht;
    209       M = -DEG_RAD * cphi / stht;
    210     }
    211     if (Zenith2) {
    212       L =  DEG_RAD * sphi;
    213       M = -DEG_RAD * cphi;
    214     }
    215     if (!strcmp(type, "-ZEA") || !strcmp(type, "-ZPL")) {
    216       Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
    217       L =  Rc * sphi;
    218       M = -Rc * cphi;
    219       status = TRUE;
     199    switch (proj) {
     200      case PROJ_TAN:
     201      case PROJ_DIS:
     202        L =  DEG_RAD * sphi / stht;
     203        M = -DEG_RAD * cphi / stht;
     204        break;
     205      case PROJ_SIN:
     206        L =  DEG_RAD * sphi;
     207        M = -DEG_RAD * cphi;
     208        break;
     209      case PROJ_ZEA:
     210      case PROJ_ZPL:
     211        Rc = DEG_RAD * M_SQRT2 / sqrt (1 + stht);
     212        L =  Rc * sphi;
     213        M = -Rc * cphi;
     214        status = TRUE;
     215        break;
    220216    }
    221217  }
    222218
    223219  /**** Other Standard Projections ****/
    224   if (PseudoCyl) {
    225     if (!strcmp(type, "-AIT")) {
    226       phi = RAD_DEG*(ra - coords[0].crval1);
    227       theta = RAD_DEG*(dec - coords[0].crval2);
    228       P = 1.0 + cos (theta) * cos (0.5*phi);
    229       if (P != 0.0) {
    230         A =  DEG_RAD * sqrt (2.0 / P);
    231         L =  2.0 * A * cos (theta) * sin (0.5*phi);
    232         M =  A * sin (theta);
    233       } else {
    234         L =  0.0;
    235         M =  0.0;
    236       }
    237     }
    238     if (!strcmp(type, "-GLS")) {
    239       phi = ra - coords[0].crval1;
    240       theta = dec - coords[0].crval2;
    241       L = phi * cos(RAD_DEG * theta);
    242       M = theta;
    243     }
    244     if (!strcmp(type, "-PAR")) {
    245       phi = ra - coords[0].crval1;
    246       theta = dec - coords[0].crval2;
    247       L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
    248       M = 180.0 * sin (RAD_DEG*theta/3.0);
     220  if (mode == PROJ_MODE_PSEUDOCYL) {
     221    switch (proj) {
     222      case PROJ_AIT:
     223        phi = RAD_DEG*(ra - coords[0].crval1);
     224        theta = RAD_DEG*(dec - coords[0].crval2);
     225        P = 1.0 + cos (theta) * cos (0.5*phi);
     226        if (P != 0.0) {
     227          A =  DEG_RAD * sqrt (2.0 / P);
     228          L =  2.0 * A * cos (theta) * sin (0.5*phi);
     229          M =  A * sin (theta);
     230        } else {
     231          L =  0.0;
     232          M =  0.0;
     233        }       
     234        break;
     235      case PROJ_GLS:
     236        phi = ra - coords[0].crval1;
     237        theta = dec - coords[0].crval2;
     238        L = phi * cos(RAD_DEG * theta);
     239        M = theta;
     240        break;
     241      case PROJ_PAR:
     242        phi = ra - coords[0].crval1;
     243        theta = dec - coords[0].crval2;
     244        L = phi * (2.0*cos(2*RAD_DEG*theta/3.0) - 1);
     245        M = 180.0 * sin (RAD_DEG*theta/3.0);
     246        break;
    249247    }
    250248  }
    … …  
    270268      dL = (L - Lo);
    271269      dM = (M - Mo);
    272       // fprintf (stderr, "%d: %f,%f : %f,%f : %f,%f : %f,%f\n", i, L, M, X, Y, Lo, Mo, dL, dM);
    273270
    274271      X += determ * (coords[0].pc2_2*dL - coords[0].pc1_2*dM);
    … …  
    282279
    283280  return (status);
    284  
    285281}
    286282
    … …  
    599595
    600596*/
     597
     598  /* PLY is equiv to LIN with higher order terms
     599     ZPL is equiv to ZEA with higher order terms
     600     DIS is equiv to TAN with higher order terms
     601     WRP is equiv to PLY, with implied mosaic */
     602
     603int GetProjection (char *ctype) {
     604  if (!strcmp(&ctype[4], "-ZEA")) return PROJ_ZEA;
     605  if (!strcmp(&ctype[4], "-ZPL")) return PROJ_ZPL;
     606  if (!strcmp(&ctype[4], "-ARC")) return PROJ_ARC;
     607  if (!strcmp(&ctype[4], "-STG")) return PROJ_STG;
     608  if (!strcmp(&ctype[4], "-SIN")) return PROJ_SIN;
     609  if (!strcmp(&ctype[0], "MM"))   return PROJ_SIN; // note ctype[0]
     610  if (!strcmp(&ctype[4], "-TAN")) return PROJ_TAN;
     611  if (!strcmp(&ctype[4], "-DIS")) return PROJ_DIS;
     612  if (!strcmp(&ctype[4], "-LIN")) return PROJ_LIN;
     613  if (!strcmp(&ctype[0], "GENE")) return PROJ_LIN; // note ctype[0]
     614  if (!strcmp(&ctype[4], "-PLY")) return PROJ_PLY;
     615  if (!strcmp(&ctype[4], "-WRP")) return PROJ_WRP;
     616  if (!strcmp(&ctype[4], "-AIT")) return PROJ_AIT;
     617  if (!strcmp(&ctype[4], "-GLS")) return PROJ_GLS;
     618  if (!strcmp(&ctype[4], "-PAR")) return PROJ_PAR;
     619  return PROJ_NONE;
     620}
     621 
     622int SetProjection (char *ctype, int proj) {
     623  switch (proj) {
     624    case PROJ_ZEA: strcpy(&ctype[4], "-ZEA") return TRUE;
     625    case PROJ_ZPL: strcpy(&ctype[4], "-ZPL") return TRUE;
     626    case PROJ_ARC: strcpy(&ctype[4], "-ARC") return TRUE;
     627    case PROJ_STG: strcpy(&ctype[4], "-STG") return TRUE;
     628    case PROJ_SIN: strcpy(&ctype[4], "-SIN") return TRUE;
     629    case PROJ_TAN: strcpy(&ctype[4], "-TAN") return TRUE;
     630    case PROJ_DIS: strcpy(&ctype[4], "-DIS") return TRUE;
     631    case PROJ_LIN: strcpy(&ctype[4], "-LIN") return TRUE;
     632    case PROJ_PLY: strcpy(&ctype[4], "-PLY") return TRUE;
     633    case PROJ_WRP: strcpy(&ctype[4], "-WRP") return TRUE;
     634    case PROJ_AIT: strcpy(&ctype[4], "-AIT") return TRUE;
     635    case PROJ_GLS: strcpy(&ctype[4], "-GLS") return TRUE;
     636    case PROJ_PAR: strcpy(&ctype[4], "-PAR") return TRUE;
     637  }
     638  return FALSE;
     639} 
     640
     641int GetProjectionMode (int proj) {
     642  switch (proj) {
     643    case PROJ_ZEA:
     644    case PROJ_ZPL:
     645    case PROJ_ARC:
     646    case PROJ_STG:
     647    case PROJ_SIN:
     648    case PROJ_TAN:
     649    case PROJ_DIS:
     650      return PROJ_MODE_ZENITHAL;
     651    case PROJ_LIN:
     652    case PROJ_PLY:
     653    case PROJ_WRP:
     654      return PROJ_MODE_CARTESIAN;
     655    case PROJ_AIT:
     656    case PROJ_GLS:
     657    case PROJ_PAR:
     658      return PROJ_MODE_PSEUDOCLY;
     659  }
     660  return PROJ_MODE_NONE;
     661}
     662
Note: See TracChangeset for help on using the changeset viewer.