Changeset 1293 for trunk/psLib/src/astro

Timestamp:

Jul 23, 2004, 4:46:45 PM (22 years ago)

Author:

gusciora

Message:

...

Location:

trunk/psLib/src/astro

Files:

• psCoord.c (modified) (5 diffs)
• psCoord.h (modified) (5 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -22 +22 @@
 #include "psCoord.h"
 #include "psMemory.h"
-
+#include "psAbort.h"
+#include <math.h>
+#include <float.h>
 
 psPlane *psPlaneTransformApply(psPlane *out,
@@ -127 +129 @@
 }
 
+// XXX: Is this the correct way to calculate this?
+float cot(float x)
+{
+    return(1.0 / atan(x));
+}
+
+float arg(float x, float y)
+{
+    if (x > 0) {
+        return(atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(-0.5 * M_PI);
+    } else if ((x==0) && (y==0)) {
+        return(M_PI + atan(y/x));
+    } else if ((x==0) && (y==0)) {
+        return(-M_PI + atan(y/x));
+    }
+
+    psAbort(__func__, "Unacceptable range for (arg(%f, %f).\n", x, y);
+    return(0.0);
+}
+
+
+psPlane *psProject(const psSphere *coord,
+                   const psProjection *projection)
+{
+    float R = 0.0;
+    float alpha = 0.0;
+    psPlane *tmp= (psPlane *) psAlloc(sizeof(psPlane));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = cot(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = cos(coord->r) * (180.0 / M_PI);
+        tmp->x = R * sin(coord->d);
+        tmp->y = R * cos(coord->d);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->x = coord->d;
+        tmp->y = coord->r;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->x = coord->d;
+        tmp->y = log(tan(45.0 + (0.5 * coord->r))) * 180.0/M_PI;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        alpha = 1.0 / ((180.0 / M_PI) *
+                       sqrt(1.0 + (cos(coord->r) * cos(0.5 * coord->d) * 0.5)));
+
+        tmp->x = 2.0 * alpha * cos(coord->r) * sin(0.5 * coord->d);
+        tmp->y = alpha * sin(coord->d);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
+
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection)
+{
+    float R = 0.0;
+    float chu = 0.0;
+    float chu1 = 0.0;
+    float chu2 = 0.0;
+    psSphere *tmp= (psSphere *) psAlloc(sizeof(psSphere));
+
+    if (projection->type == PS_PROJ_TAN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = atan(180.0 / (R * M_PI));
+
+    } else if (projection->type == PS_PROJ_SIN) {
+        R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
+        tmp->d = arg(-coord->y, coord->x);
+        tmp->r = acos((R * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_CAR) {
+        tmp->d = coord->x;
+        tmp->r = coord->y;
+
+    } else if (projection->type == PS_PROJ_MER) {
+        tmp->d = coord->x;
+        tmp->r = (2.0 * atan(exp((coord->y * M_PI/180.0)))) - 180.0;
+
+    } else if (projection->type == PS_PROJ_AIT) {
+        chu1 = (coord->x * M_PI) / 720.0;
+        chu1*= chu1;
+        chu2 = (coord->y * M_PI) / 360.0;
+        chu2*= chu2;
+        chu = sqrt(1.0 - chu1 - chu2);
+        tmp->d = 2.0 * arg((2.0 * chu * chu) - 1.0,
+                           (coord->x * chu * M_PI)/360.0);
+        tmp->r = asin((coord->y * chu * M_PI) / 180.0);
+
+    } else if (projection->type == PS_PROJ_PAR) {
+        psAbort(__func__, "The projection type PS_PROJ_PAR is undefined.\n");
+
+    } else if (projection->type == PS_PROJ_GLS) {
+        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+    }
+
+    return(tmp);
+}
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,
@@ -133 +248 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
-
-
+    psAbort(__func__, "This function has not be dedfined.\n");
+    return(NULL);
+}
+
+
+// XXX: I copied the algorithm from the ADD exactly.  Arguments mode and unit
+// are ignored.
 psSphere *psSphereSetOffset(const psSphere *restrict position,
                             const psSphere *restrict offset,
@@ -142 +260 @@
                             psSphereOffsetUnit unit)
 {
-    return(NULL);
-}
+    psPlane lin;
+    psSphere *tmp;
+    psProjection proj;
+
+    if (mode == PS_LINEAR) {
+        proj.R = position->r;
+        proj.D = position->d;
+        proj.Xs = 0.0;
+        proj.Ys = 0.0;
+        proj.type = PS_PROJ_TAN;
+
+        lin.x = offset->r;
+        lin.y = offset->d;
+
+        tmp = psDeproject(&lin, &proj);
+        return(tmp);
+    } else if (mode == PS_SPHERICAL) {
+        psAbort(__func__, "Sperical offset modes are not defined.\n");
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
+    return(NULL);
+}

trunk/psLib/src/astro/psCoord.h

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-23 03:13:39 $
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-24 02:42:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27 +27 @@
 typedef struct
 {
-    double x;    ///< x position
-    double y;    ///< y position
+    double x;      ///< x position
+    double y;      ///< y position
     double xErr;   ///< Error in x position
     double yErr;   ///< Error in y position
@@ -36 +36 @@
 typedef struct
 {
-    double r;    ///< RA
-    double d;    ///< Dec
+    double r;      ///< RA
+    double d;      ///< Dec
     double rErr;   ///< Error in RA
     double dErr;   ///< Error in Dec
@@ -63 +63 @@
     double sinNPlat;   ///< sin of North Pole latitude
     double cosNPlat;   ///< cos of North Pole latitude
-    double sinZP;   ///< sin of Forst PT os Ares lon
-    double cosZP;   ///< cos of Forst PT os Ares lon
+    double sinZP;      ///< sin of Forst PT os Ares lon
+    double cosZP;      ///< cos of Forst PT os Ares lon
 }
 psSphereTransform;
 
 typedef enum {
-    PS_PROJ_TAN,   ///< Tangent projection
-    PS_PROJ_SIN,   ///< Sine projection
-    PS_PROJ_AIT,   ///< Aitoff projection
-    PS_PROJ_PAR,   ///< Par projection
-    PS_PROJ_GLS,   ///< GLS projection
-    PS_PROJ_NTYPE   ///< Number of types; must be last.
+    PS_PROJ_TAN,       ///< Tangent projection
+    PS_PROJ_SIN,       ///< Sine projection
+    PS_PROJ_AIT,       ///< Aitoff projection
+    PS_PROJ_PAR,       ///< Par projection
+    PS_PROJ_GLS,       ///< GLS projection
+    PS_PROJ_CAR,       ///< CAR projection
+    PS_PROJ_MER,       ///< MER projection
+    PS_PROJ_NTYPE      ///< Number of types; must be last.
 } psProjectionType;
 
 typedef struct
 {
-    double R;    ///< Coordinates of projection center
-    double D;    ///< Coordinates of projection center
+    double R;     ///< Coordinates of projection center
+    double D;     ///< Coordinates of projection center
     double Xs;    ///< plate-scale in X direction
     double Ys;    ///< plate-scale in Y direction
@@ -128 +130 @@
                    const psProjection *projection);
 
-psSphere *psProject(const psPlane *coord,
-                    const psProjection *projection);
+psSphere *psDeproject(const psPlane *coord,
+                      const psProjection *projection);
 
 psSphere *psSphereGetOffset(const psSphere *restrict position1,