Changeset 1296 for trunk/psLib/src/astronomy/psCoord.c

Timestamp:

Jul 26, 2004, 7:28:05 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/astronomy/psCoord.c (modified) (12 diffs)

trunk/psLib/src/astronomy/psCoord.c

@@ -10 +10 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-24 02:42:59 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-07-27 05:28:05 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -26 +26 @@
 #include <float.h>
 
+// This is the only function in this file which I understand.
 psPlane *psPlaneTransformApply(psPlane *out,
                                const psPlaneTransform *transform,
@@ -72 +73 @@
 
 
-psSphereTransform *psSphereTransformAlloc(double NPlon,
-        double NPlat,
-        double ZP)
+psSphereTransform *psSphereTransformAlloc(double NPlat,
+        double Xo,
+        double xo)
 {
     psSphereTransform *tmp = (psSphereTransform *) psAlloc(sizeof(psSphereTransform));
 
-    tmp->sinNPlon = sin(NPlon);
-    tmp->cosNPlon = cos(NPlon);
-    tmp->sinNPlat = sin(NPlat);
-    tmp->cosNPlat = cos(NPlat);
-    tmp->sinZP = sin(ZP);
-    tmp->cosZP = cos(ZP);
+    tmp->sinPhi = sin(NPlat);
+    tmp->cosPhi = cos(NPlat);
+    tmp->Xo = Xo;
+    tmp->xo = xo;
 
     return(tmp);
 }
 
-
+// I understand this one too.
 void p_psSphereTransformFree(psSphereTransform *trans)
 {
@@ -95 +94 @@
 
 
-
+// XXX: I have no idea how this works.  This algorithm comes from an email
+// from Gene.  I assume (x,y) corresponds to (r,d) in the sphere coordinates.
+
+// XXX: In Gene's email, there are different variables with similar names
+// (y and Y) and in the code, there's cos_y, cos_Y, and cos(y).  Verify that
+// there are no typo's.
 
 psSphere *psSphereTransformApply(psSphere *out,
@@ -101 +105 @@
                                  const psSphere *coord)
 {
-    return(NULL);
-}
-
-
+    double sinY = 0.0;
+    ;
+    double cosY = 0.0;
+    ;
+    double sinX = 0.0;
+    ;
+    double cosX = 0.0;
+    ;
+    double x = 0.0;
+    double y = 0.0;
+    double dx = 0.0;
+
+    if (out == NULL) {
+        out = (psSphere *) psAlloc(sizeof(psSphere));
+    }
+
+    x = coord->r;
+    y = coord->d;
+    dx    = x - transform->xo;
+    sinY = cos(y)*sin(dx)*transform->sinPhi + sin(y)*transform->cosPhi;
+    cosY = sqrt(1.0 - sinY*sinY);
+    sinX = (cos(y)*sin(dx)*transform->cosPhi - sin(y)*transform->sinPhi) /
+           cos(y);
+    cosX = cos(y)*cos(dx) / cos(y);
+
+    out->r = atan2(sinX, cosX) + transform->Xo;
+    out->d = atan2(sinY, cosY);
+
+    return(out);
+}
+
+// XXX: Add a date argument to determine the correct year.
 psSphereTransform *psSphereTransformICRStoEcliptic(void)
 {
-    // sla_EQECL();
-    return(NULL);
-}
-
+    double year = 2004.0;
+    double T = year / 100.0;
+    double phi = -23.452294 + 0.013013*T + 0.000001639*T*T - 0.000000503*T*T*T;
+    double Xo = 0.0;
+    double xo = 0.0;
+
+    return(psSphereTransformAlloc(phi, Xo, xo));
+}
+
+// XXX: Add a date argument to determine the correct year.
 psSphereTransform *psSphereTransformEcliptictoICRS(void)
 {
-    // sla_ECLEQ();
-    return(NULL);
+    double year = 2004.0;
+    double T = year / 100.0;
+    double phi = +23.452294 - 0.013013*T - 0.000001639*T*T + 0.000000503*T*T*T;
+    double Xo = 0.0;
+    double xo = 0.0;
+
+    return(psSphereTransformAlloc(phi, Xo, xo));
 }
 
 psSphereTransform *psSphereTransformICRStoGalatic(void)
 {
-    // sla_EQGAL();
-    return(NULL);
+    return(psSphereTransformAlloc(62.6, 282.25, 33.0));
 }
 
 psSphereTransform *psSphereTransformGalatictoICRS(void)
 {
-    // sla_GALEQL();
-    return(NULL);
+    return(psSphereTransformAlloc(-62.6, 33.0, 282.25));
 }
 
@@ -135 +176 @@
 }
 
+// This is some kind of arc tan function.
 float arg(float x, float y)
 {
@@ -153 +195 @@
 }
 
-
+// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
+// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
 psPlane *psProject(const psSphere *coord,
                    const psProjection *projection)
@@ -190 +233 @@
 
     } else if (projection->type == PS_PROJ_GLS) {
-        psAbort(__func__, "The projection type PS_PROJ_GLG is undefined.\n");
+        psAbort(__func__, "The projection type PS_PROJ_GLS is undefined.\n");
     }
 
@@ -196 +239 @@
 }
 
+// XXX: Waiting for the definition of the PS_PROJ_PAR projection.
+// XXX: Waiting for the definition of the PS_PROJ_GLS projection.
 psSphere *psDeproject(const psPlane *coord,
                       const psProjection *projection)
@@ -248 +293 @@
                             psSphereOffsetUnit unit)
 {
-    psAbort(__func__, "This function has not be dedfined.\n");
+    //    psPlane lin;
+    psSphere *tmp;
+    //    psProjection proj;
+    double tmpR = 0.0;
+    double tmpD = 0.0;
+
+    if (mode == PS_LINEAR) {
+        // XXX: I have no idea how to construct this.
+        return(NULL);
+    } else if (mode == PS_SPHERICAL) {
+        tmpR = position2->r - position1->r;
+        tmpD = position2->d - position1->d;
+
+        if (unit == PS_ARCSEC) {
+            tmpR = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
+            tmpD = (tmpR * 180.0 * 60.0 * 60.0) / M_PI;
+        } else if (unit == PS_ARCMIN) {
+            tmpR = (tmpR * 180.0 * 60.0) / M_PI;
+            tmpD = (tmpR * 180.0 * 60.0) / M_PI;
+        } else if (unit == PS_DEGREE) {
+            tmpR = (tmpR * 180.0) / M_PI;
+            tmpD = (tmpR * 180.0) / M_PI;
+        } else if (unit == PS_RADIAN) {}
+        else {
+            psAbort(__func__, "Unknow offset unit: 0x%x\n", unit);
+        }
+
+        tmp = (psSphere *) psAlloc(sizeof(psSphere));
+        tmp->r = tmpR;
+        tmp->d = tmpD;
+        tmp->rErr = 0.0;
+        tmp->dErr = 0.0;
+        // XXX: Do we need to wrap these to an acceptable range?
+        return(tmp);
+    }
+    psAbort(__func__, "Unrecognized offset mode\n");
     return(NULL);
 }
 
 
-// XXX: I copied the algorithm from the ADD exactly.  Arguments mode and unit
-// are ignored.
+// XXX: I copied the algorithm from the ADD exactly.
 psSphere *psSphereSetOffset(const psSphere *restrict position,
                             const psSphere *restrict offset,
@@ -263 +342 @@
     psSphere *tmp;
     psProjection proj;
+    double tmpR = 0.0;
+    double tmpD = 0.0;
 
     if (mode == PS_LINEAR) {
         proj.R = position->r;
         proj.D = position->d;
-        proj.Xs = 0.0;
-        proj.Ys = 0.0;
+        proj.Xs = 1.0;
+        proj.Ys = 1.0;
         proj.type = PS_PROJ_TAN;
 
@@ -276 +357 @@
         tmp = psDeproject(&lin, &proj);
         return(tmp);
+
     } else if (mode == PS_SPHERICAL) {
-        psAbort(__func__, "Sperical offset modes are not defined.\n");
+        if (unit == PS_ARCSEC) {
+            tmpR = (M_PI * offset->r) / (180.0 * 60.0 * 60.0);
+            tmpD = (M_PI * offset->d) / (180.0 * 60.0 * 60.0);
+        } else if (unit == PS_ARCMIN) {
+            tmpR = (M_PI * offset->r) / (180.0 * 60.0);
+            tmpD = (M_PI * offset->d) / (180.0 * 60.0);
+        } else if (unit == PS_DEGREE) {
+            tmpR = (M_PI * offset->r) / (180.0);
+            tmpD = (M_PI * offset->d) / (180.0);
+        } else if (unit == PS_RADIAN) {
+            tmpR = offset->r;
+            tmpD = offset->d;
+        } else {
+            psAbort(__func__, "Unknow offset unit: 0x%x\n", unit);
+        }
+
+        tmp = (psSphere *) psAlloc(sizeof(psSphere));
+        tmp->r = position->r + tmpR;
+        tmp->r = position->d + tmpD;
+        tmp->rErr = 0.0;
+        tmp->dErr = 0.0;
+
+        // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI)
+        // and (0 to 2*PI).
+        return(tmp);
     }
     psAbort(__func__, "Unrecognized offset mode\n");