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

Timestamp:

Jan 17, 2006, 2:41:29 PM (21 years ago)

Author:

drobbin

Message:

Updated Table parser, Fixed misc EOC, worked primarily on GravityDeflection

File:

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

trunk/psLib/src/astro/psEarthOrientation.c

@@ -8 +8 @@
  *  @author Robert Daniel DeSonia, MHPCC
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-12 20:40:12 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-18 00:41:29 $
  *
  *  Copyright 2005 Maui High Performance Computing Center, University of Hawaii
@@ -39 +39 @@
 #define PS_AU 149597890000.0 /* meters */
 // Modified Julian Day 1/1/2000 00:00:00
-#define MJD_2000  51544.0
+#define MJD_2000  51544.5
 // Days in Julian century
 #define JULIAN_CENTURY 36525.0
@@ -299 +299 @@
     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);
     if (directionVector == NULL)
         printf("directionVector is null\n");
     psCube* actualVector = psSphereToCube(actual);
-    if (directionVector == NULL)
+    if (actualVector == NULL)
         printf("actualVector is null\n");
     //    mu = acos(directionVector->x*actualVector->x +
@@ -311 +309 @@
           directionVector->y*actualVector->y +
           directionVector->z*actualVector->z);
-
-    //    rp = apparent - mu * direction;
-    //    rp->r = actual->r - mu * direction->r;
-    //    rp->d = actual->d - mu * direction->d;
     rp->x = actualVector->x - mu * directionVector->x;
     rp->y = actualVector->y - mu * directionVector->y;
     rp->z = actualVector->z - mu * directionVector->z;
 
-    mu_p = mu + speed * ((mu * mu - 1.0) / (1.0 - speed * mu));
+    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)
@@ -345 +339 @@
     //XXX: Must be a sign error somewhere above?  Magnitude of change is correct but wrong way...
     //This will fix the problem but is somewhat of a hack.
-    r_p->x = actualVector->x - (r_p->x - actualVector->x);
-    r_p->y = actualVector->y - (r_p->y - actualVector->y);
-    r_p->z = actualVector->z - (r_p->z - actualVector->z);
+    //    r_p->x = actualVector->x - (r_p->x - actualVector->x);
+    //    r_p->y = actualVector->y - (r_p->y - actualVector->y);
+    //    r_p->z = actualVector->z - (r_p->z - actualVector->z);
 
     psSphere *r_pSphere = psCubeToSphere(r_p);
@@ -370 +364 @@
     PS_ASSERT_PTR_NON_NULL(sun, NULL);
 
+    psSphere *temp = psSphereAlloc();
     // calculating the apparent angle from the actual angle and the sun position
 
@@ -378 +373 @@
     // N.B., assuming the psSphereToCube function returns a unit vector.
     //    double theta = acos(sunVector->x*actualVector->x +
-    double theta = (sunVector->x*actualVector->x +
-                    sunVector->y*actualVector->y +
-                    sunVector->z*actualVector->z);
-    printf(" Theta = %lf\n", theta);
+    double dotProd = (sunVector->x*actualVector->x +
+                      sunVector->y*actualVector->y + sunVector->z*actualVector->z);
+    double theta, sunMag, actMag;
+    sunMag = sqrt(sunVector->x*sunVector->x + sunVector->y*sunVector->y +
+                  sunVector->z*sunVector->z);
+    actMag = sqrt(actualVector->x*actualVector->x + actualVector->y*actualVector->y +
+                  actualVector->z*actualVector->z);
+    dotProd = dotProd / (sunMag * actMag);
+    theta = acos(dotProd);
+
+    //    theta = acos( cos(sun->d)*cos(actual->d)*cos(sun->r-actual->r) + sin(sun->r)*sin(actual->r) );
+
+    printf(" Theta = %.13g\n", theta);
     //    theta = acos(-theta);
     //    printf("\n Theta = %lf\n", theta);
     double r0 = PS_AU * tan(theta);
+    printf(" r0 = %.19e\n", r0);
     double deflection = 4.0*PS_G*PS_M/(PS_C0*PS_C0*r0);
 
@@ -390 +395 @@
     double limit = SEC_TO_RAD(1.75);
     printf(" deflection = %.13g\n", deflection);
-    printf(" limit = %lf\n", limit);
+    //printf(" limit = %lf\n", limit);
     if (deflection > limit) {
         //       deflection = limit;
@@ -400 +405 @@
     }
 
-    if (apparent == NULL) {
-        apparent = psSphereAlloc();
-    } else {
-        apparent->r = 0.0;
-        apparent->d = 0.0;
-        apparent->rErr = 0.0;
-        apparent->dErr = 0.0;
-    }
-
+    /*    if (apparent == NULL) {
+            apparent = psSphereAlloc();
+        } else {
+            apparent->r = 0.0;
+            apparent->d = 0.0;
+            apparent->rErr = 0.0;
+            apparent->dErr = 0.0;
+        }
+    */
+    if (apparent != NULL) {
+        psFree(apparent);
+    }
     // bend the actual vector away from the sun vector by deflection angle.
     // XXX: Not sure how to do this.  Dave thinks the formula should be:
@@ -414 +422 @@
     theta = 0.0;
     double phi = 0.0;
-    //    deflection = SEC_TO_RAD(deflection);
+    //    deflection = SEC_TO_RAD(deflection) * 1e6;
+    //    deflection *= M_PI * 1e-2;
     theta = atan(r0/PS_AU) * tan(deflection);
-    printf(" Theta = %.13g\n", theta);
-    printf(" deflection = %.13g\n", deflection);
     //    phi = sqrt( deflection*deflection - theta*theta );
-    phi = deflection * cos(asin(theta/deflection)) * 3e-2;
+    //    phi = deflection * cos(asin(theta/deflection));
+
+    //    phi = sqrt(theta*theta - deflection*deflection);
+    //    phi = deflection * cos(asin(theta/deflection)) * 3e-2;
     //    phi = cos(asin(theta/deflection));
     //    phi = asin(theta/deflection);
-    apparent->r = theta;
-    apparent->d = phi;
+
+    //    apparent->r = theta;
+    //    apparent->d = phi;
+    /*
+        actualVector->x += actualVector->x*deflection;
+        actualVector->y += actualVector->y*deflection;
+        actualVector->z += actualVector->z*deflection;
+        apparent = psCubeToSphere(actualVector);
+    */
+    theta = tan(sun->r - actual->r) * deflection;
+    phi = tan(sun->d - actual->d) * deflection;
+
+    printf(" Theta = %.13g\n", theta);
+    printf(" phi = %.13g\n", phi);
+
+    temp->r = theta;
+    temp->d = phi;
+    apparent = psSphereSetOffset(actual, temp, PS_SPHERICAL, PS_RADIAN);
 
     psFree(actualVector);
     psFree(sunVector);
+    psFree(temp);
+
     return apparent;
 }
@@ -807 +835 @@
 }
 
-psSphereRot* psSphereRot_TEOtoCEO(const psTime *time)
+psSphereRot* psSphereRot_TEOtoCEO(const psTime *time,
+                                  psEarthPole *tidalCorr)
 {
     PS_ASSERT_PTR_NON_NULL(time,NULL);
@@ -815 +844 @@
         in = psTimeConvert(in, PS_TIME_UT1);
     }
-    //    double T = (double)(in->sec) + (double)(in->nsec / 1e9);
+    if (tidalCorr != NULL && tidalCorr->s != 0.0) {
+        int nsec = in->nsec + (int)(tidalCorr->s * 1e9);
+        if (nsec < 0.0) {
+            in->sec += -1;
+            in->nsec = (int)(1e9) - nsec;
+        } else {
+            in->nsec = nsec;
+        }
+    }
     double T = psTimeToJD(in);
     T += -2451545.0;
     double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T);
     psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0);
-    //    psSphereRot *out = psSphereRotInvert(theta, 0.0, 0.0);
 
     psFree(in);
@@ -919 +955 @@
     out->x = SEC_TO_RAD(xOut);
     out->y = SEC_TO_RAD(yOut);
+    //    psEarthPole *polarTideCorr = psEOC_PolarTideCorr(time);
+    //    out->x += polarTideCorr->x;
+    //    out->y += polarTideCorr->y;
+    //    psFree(polarTideCorr);
+
     //    out->s = SEC_TO_RAD(sOut);
 
@@ -1034 +1075 @@
     CORX = SEC_TO_RAD(CORX);
     CORY = SEC_TO_RAD(CORY);
-    CORZ = SEC_TO_RAD(CORZ);
+    //    CORZ = SEC_TO_RAD(CORZ);
 
     out->x = CORX;