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

Timestamp:

Oct 25, 2005, 3:20:15 PM (21 years ago)

Author:

drobbin

Message:

Finished and tested psGravityDeflection function

File:

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

trunk/psLib/src/astro/psEarthOrientation.c

@@ -6 +6 @@
 *  @ingroup EarthOrientation
 *
+*  @author Dave Robbins, MHPCC
 *  @author Robert Daniel DeSonia, MHPCC
 *
-*  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-10-07 21:27:49 $
+*  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-10-26 01:20:15 $
 *
 *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -166 +167 @@
 psSphere *psGravityDeflection(psSphere *apparent, psSphere *actual, psSphere *sun)
 {
-    // calculating the apparent angle from the actual angle and the sun position
-
-    // first, calculate the angle between the sun vector and the actual vector
+    PS_ASSERT_PTR_NON_NULL(apparent, NULL);
+    PS_ASSERT_PTR_NON_NULL(sun, NULL);
+
+    // calculating the actual angle from the apparent angle and the sun position
+
+    // first, calculate the angle between the sun vector and the apparent vector
 
     // Moving to cartesian first:  XXX -- is this required?
     psCube* sunVector = psSphereToCube(sun);
-    psCube* actualVector = psSphereToCube(actual);
+    psCube* apparentVector = psSphereToCube(apparent);
 
     // use dot product to calculate the angle of separation
     // N.B., assuming the psSphereToCube function returns a unit vector.
-    double theta = acos(sunVector->x*actualVector->x +
-                        sunVector->y*actualVector->y +
-                        sunVector->z*actualVector->z);
+    double theta = acos(sunVector->x*apparentVector->x +
+                        sunVector->y*apparentVector->y +
+                        sunVector->z*apparentVector->z);
 
     double r0 = PS_AU * tan(theta);
-
     double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
 
@@ -187 +190 @@
     double limit = SEC_TO_RAD(1.75);
     if (deflection > limit) {
-        deflection = limit;
-    }
-
-    // bend the actual vector away from the sun vector by deflection angle.
+        //       deflection = limit;
+        //if deflection is greater than limit, the light rays will hit the sun
+        psWarning("Invalid positions.  Light ray will not be seen on earth.\n");
+        psFree(apparentVector);
+        psFree(sunVector);
+        return actual;
+    }
+
+    if (actual == NULL) {
+        actual = psSphereAlloc();
+    } else {
+        actual->r = 0.0;
+        actual->d = 0.0;
+        actual->rErr = 0.0;
+        actual->dErr = 0.0;
+    }
+
+    // bend the apparent vector away from the sun vector by deflection angle.
     // XXX: Not sure how to do this.  Dave thinks the formula should be:
     //      theta = atan(r0/d)*tan(deflection), phi = thete/tan(deflection)
-
-    return NULL;
+    theta = 0.0;
+    double phi = 0.0;
+    deflection = SEC_TO_RAD(deflection);
+    theta = atan(r0/PS_AU) * tan(deflection);
+    phi = sqrt( deflection*deflection - theta*theta );
+    actual->r = theta;
+    actual->d = phi;
+    psFree(apparentVector);
+    psFree(sunVector);
+    return actual;
 }