Changeset 5450 for trunk/psLib/src/astro/psEarthOrientation.c

Timestamp:

Oct 27, 2005, 4:25:22 PM (21 years ago)

Author:

drobbin

Message:

Implemented and tested psAberration, but results have not been verified.

File:

• trunk/psLib/src/astro/psEarthOrientation.c (modified) (3 diffs)

trunk/psLib/src/astro/psEarthOrientation.c

@@ -9 +9 @@
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-10-27 20:38:18 $
+*  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-10-28 02:25:22 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -19 +19 @@
 
 #include "psEarthOrientation.h"
-#include "psTime.h"
+//#include "psTime.h"
 #include "psArray.h"
 #include "psPolynomial.h"
@@ -163 +163 @@
                        double speed)
 {
-
-
-    return NULL;
+    PS_ASSERT_PTR_NON_NULL(apparent, NULL);
+    PS_ASSERT_PTR_NON_NULL(direction, NULL);
+    if (fabs(speed) < DBL_EPSILON) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Aberration speed should not be equal to 0.\n");
+        return NULL;
+    }
+
+    if (actual == NULL) {
+        actual = psSphereAlloc();
+    } else {
+        actual->r = 0.0;
+        actual->d = 0.0;
+        actual->rErr = 0.0;
+        actual->dErr = 0.0;
+    }
+    psSphere *rp = psSphereAlloc();
+    psSphere *r_p = psSphereAlloc();
+    double mu = 0.0;
+    double mu_p = 0.0;
+    double a = 0.0;
+
+    //mu = apparent * direction;  I believe this should be:
+    //    mu = apparent->r * direction->r + apparent->d * direction->d;
+    psCube* directionVector = psSphereToCube(direction);
+    psCube* apparentVector = psSphereToCube(apparent);
+    mu = acos(directionVector->x*apparentVector->x +
+              directionVector->y*apparentVector->y +
+              directionVector->z*apparentVector->z);
+
+    //rp = apparent - mu * direction;
+    rp->r = apparent->r - mu * direction->r;
+    rp->d = apparent->d - mu * direction->d;
+
+    mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu));
+
+    //Not sure if this is right.  ADD gets a scalar from division of rp (a vector)
+    psCube* rpVector = psSphereToCube(rp);
+    a = sqrt( (1.0 - mu_p * mu_p) /
+              (sqrt(rpVector->x*rpVector->x
+                    + rpVector->y*rpVector->y
+                    + rpVector->z*rpVector->z)) );
+
+    //r_p = mu_p * direction + a * rp;
+    r_p->r = mu_p * direction->r + a * rp->r;
+    r_p->d = mu_p * direction->d + a * rp->d;
+
+    *actual = *r_p;
+    /*
+        psSphereRot *rot = NULL;
+        double cosR = cos(direction->r);
+        double cosD = cos(direction->d);
+        double sinR = sin(direction->r);
+        double sinD = sin(direction->d);
+        rot = psSphereRotQuat(cosR*cosD, sinR*cosD, sinD, speed);
+     
+        actual = psSphereRotApply(actual, rot, apparent);
+    */
+    psFree(rp);
+    psFree(r_p);
+    psFree(directionVector);
+    psFree(apparentVector);
+    psFree(rpVector);
+    return actual;
 }