IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16935


Ignore:
Timestamp:
Mar 11, 2008, 11:11:54 AM (18 years ago)
Author:
eugene
Message:

better api for coordinate transforms

Location:
trunk/Ohana/src/opihi
Files:
1 added
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/Makefile

    r12879 r16935  
    2222$(SRC)/cplot.$(ARCH).o             \
    2323$(SRC)/csystem.$(ARCH).o           \
     24$(SRC)/coord_systems.$(ARCH).o     \
    2425$(SRC)/ctimes.$(ARCH).o    \
    2526$(SRC)/cval.$(ARCH).o              \

  • trunk/Ohana/src/opihi/cmd.astro/csystem.c

    r8427 r16935  
    55  /* USAGE: csystem [C/G/E/H] [C/G/E/H] [epoch] */
    66  int i;
    7   double X, Y, Xo, xo, phi, T;
    8   double sin_x, sin_y, cos_x, cos_y;
    9   float *x, *y;
    10   struct timeval now;
    11   struct tm *local;
     7  double X, Y, x, y;
     8  float *xptr, *yptr;
    129  Vector *xvec, *yvec;
    13    
     10  CoordTransformSystem input, output;
     11  CoordTransform *transform;
     12
    1413  if (argc != 5) {
    1514    gprint (GP_ERR, "USAGE: csystems [C/G/E/H] [C/G/E/H] X Y\n");
    … …  
    1716  }
    1817
    19   Xo = xo = phi = 0;
    2018  switch (argv[1][0]) {
    21   case 'C':
    22     switch (argv[2][0]) {
    23     case 'C':
    24       gprint (GP_ERR, "same coordinate system\n");
    25       return (TRUE);
    26       break;
    27     case 'G':
    28       phi = -62.6*RAD_DEG;
    29       Xo = 282.25;
    30       xo = 33;
    31       break;
    32     case 'E':
    33       gettimeofday (&now, (struct timezone *) NULL);
    34       local = localtime (&now.tv_sec);
    35       T = local[0].tm_year / 100.0;
    36       phi = -1*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
    37       phi *= RAD_DEG;
    38       Xo = xo = 0.0;
    39       break;
    40     }
    41     break;
    42   case 'E':
    43     switch (argv[2][0]) {
    44     case 'C':
    45       gettimeofday (&now, (struct timezone *) NULL);
    46       local = localtime (&now.tv_sec);
    47       T = local[0].tm_year / 100.0;
    48       phi = 23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T;
    49       phi *= RAD_DEG;
    50       Xo = xo = 0.0;
    51       break;
    52     case 'G':
    53       gprint (GP_ERR, "error: conversions between galactic and ecliptic not implemented\n");
    54       return (FALSE);
    55       phi = -62.6*RAD_DEG;
    56       Xo = 282.25;
    57       xo = 33;
    58       break;
    59     case 'E':
    60       phi = Xo = xo = 0.0;
    61       break;
    62     }
    63     break;
    64   case 'G':
    65     switch (argv[2][0]) {
    66     case 'C':
    67       phi = 62.6*RAD_DEG;
    68       Xo = 33;
    69       xo = 282.25;
    70       break;
    71     case 'G':
    72       phi = Xo = xo = 0.0;
    73       break;
    74     case 'E':
    75       gprint (GP_ERR, "error: conversions between galactic and ecliptic not implemented\n");
    76       return (FALSE);
    77       gettimeofday (&now, (struct timezone *) NULL);
    78       local = localtime (&now.tv_sec);
    79       T = local[0].tm_year / 100.0;
    80       phi = -1*(23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T);
    81       Xo = xo = 0.0;
    82       break;
    83     }
     19    case 'C': input = COORD_CELESTIAL; break;
     20    case 'G': input = COORD_GALACTIC; break;
     21    case 'E': input = COORD_ECLIPTIC; break;
     22    default: abort();
    8423  }
    85  
    86   Xo *= RAD_DEG;
    8724
     25  switch (argv[2][0]) {
     26    case 'C': output = COORD_CELESTIAL; break;
     27    case 'G': output = COORD_GALACTIC; break;
     28    case 'E': output = COORD_ECLIPTIC; break;
     29    default: abort();
     30  }
     31
     32  transform = InitTransform (input, output);
     33  if (transform == NULL) {
     34    gprint (GP_ERR, "transform %c to %c is not yet defined\n", argv[1][0], argv[2][0]);
     35    return (FALSE);
     36  }
     37   
    8838  if (SelectScalar (argv[3], &X)) {
    89       if (!SelectScalar (argv[4], &Y)) return (FALSE);
     39    if (!SelectScalar (argv[4], &Y)) return (FALSE);
    9040     
    91       X *= RAD_DEG;
    92       Y *= RAD_DEG;
     41    ApplyTransform (&x, &y, X, Y, transform);
    9342
    94       sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);
    95       cos_y = sqrt (1 - sin_y*sin_y);
    96       sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
    97       cos_x = cos(Y)*cos(X - Xo) / cos_y;
    98      
    99       X = (DEG_RAD * atan2 (sin_x, cos_x) + xo + 360);
    100      
    101       while (X >= 360.0)
    102           X -= 360;
    103       Y = DEG_RAD * atan2 (sin_y, cos_y);
    104 
    105       gprint (GP_LOG, "%10.6f %10.6f\n", X, Y);
    106       return (TRUE);
     43    gprint (GP_LOG, "%10.6f %10.6f\n", x, y);
     44    return (TRUE);
    10745  }
    10846
    … …  
    11654  }
    11755 
    118   x = xvec[0].elements;
    119   y = yvec[0].elements;
     56  xptr = xvec[0].elements;
     57  yptr = yvec[0].elements;
    12058
    121   for (i = 0; i < xvec[0].Nelements; i++, x++, y++) {
    122     X = *x*RAD_DEG;
    123     Y = *y*RAD_DEG;
    124 
    125     sin_y = cos(Y)*sin(X - Xo)*sin(phi) + sin(Y)*cos(phi);
    126     cos_y = sqrt (1 - sin_y*sin_y);
    127     sin_x = (cos(Y)*sin(X - Xo)*cos(phi) - sin(Y)*sin(phi)) /  cos_y;
    128     cos_x = cos(Y)*cos(X - Xo) / cos_y;
    129    
    130     X = (DEG_RAD * atan2 (sin_x, cos_x) + xo + 360);
    131    
    132     while (X >= 360.0)
    133       X -= 360;
    134     Y = DEG_RAD * atan2 (sin_y, cos_y);
    135     *x = X;
    136     *y = Y;
     59  for (i = 0; i < xvec[0].Nelements; i++, xptr++, yptr++) {
     60    // ApplyTransform takes (double *), but xptr, yptr are (float *)
     61    ApplyTransform (&x, &y, *xptr, *yptr, transform);
     62    *xptr = x;
     63    *yptr = y;
    13764  }
    13865

  • trunk/Ohana/src/opihi/include/astro.h

    r4689 r16935  
    1111void InitAstro ();
    1212
     13typedef struct {
     14  int    isIdentity;          // identity transformation
     15  double phi;                 // saved in radians
     16  double Xo;                  // saved in radians
     17  double xo;                  // saved in degrees
     18  double sin_phi_cos_Xo;      // pre-computed values
     19  double sin_phi_sin_Xo;      // pre-computed values
     20  double cos_phi_cos_Xo;      // pre-computed values
     21  double cos_phi_sin_Xo;      // pre-computed values
     22  double cos_phi;             // pre-computed values
     23  double sin_phi;             // pre-computed values
     24  double cos_Xo;              // pre-computed values
     25  double sin_Xo;              // pre-computed values
     26} CoordTransform;
     27
     28typedef enum {COORD_NONE, COORD_CELESTIAL, COORD_GALACTIC, COORD_ECLIPTIC} CoordTransformSystem;
     29
     30CoordTransform *InitTransform (CoordTransformSystem input, CoordTransformSystem output);
     31int ApplyTransform (double *x, double *y, double X, double Y, CoordTransform *transform);
     32
    1333# endif
Note: See TracChangeset for help on using the changeset viewer.