Index: trunk/psLib/src/astro/psEarthOrientation.c =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.c (revision 6409) +++ trunk/psLib/src/astro/psEarthOrientation.c (revision 6425) @@ -8,6 +8,6 @@ * @author Robert Daniel DeSonia, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-06 21:57:04 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-14 00:09:26 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -356,16 +356,13 @@ psSphere *sun) { + //Check for valid inputs PS_ASSERT_PTR_NON_NULL(actual, NULL); PS_ASSERT_PTR_NON_NULL(sun, NULL); psSphere *temp = psSphereAlloc(); - // calculating the apparent angle from the actual angle and the sun position - psCube* sunVector = psSphereToCube(sun); psCube* actualVector = psSphereToCube(actual); // 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 + double dotProd = (sunVector->x*actualVector->x + sunVector->y*actualVector->y + sunVector->z*actualVector->z); @@ -378,9 +375,5 @@ 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); @@ -389,8 +382,5 @@ // make sure the deflection is not greater than 1.75 arcsec double limit = SEC_TO_RAD(1.75); - printf(" deflection = %.13g\n", deflection); - //printf(" limit = %lf\n", limit); if (deflection > limit) { - // 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"); @@ -400,19 +390,9 @@ } - /* 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: - // theta = atan(r0/d)*tan(deflection), phi = thete/tan(deflection) theta = 0.0; double phi = 0.0; @@ -477,7 +457,7 @@ double t4 = t*t*t*t; - // N.B., the following formulae are from the ADD - // Mean Anomaly of the Moon + // The following formulae are from the ADD (& s5.8 of IERS Technical Note 32): double F[14]; + // Mean Anomaly of the Moon (l) F[0] = DEG_TO_RAD(134.96340251) + SEC_TO_RAD(1717915923.2178)*t + @@ -485,6 +465,5 @@ SEC_TO_RAD(0.051635)*t3 - SEC_TO_RAD(0.00024470)*t4; - - // Mean Anomaly of the Sun + // Mean Anomaly of the Sun (l') F[1] = DEG_TO_RAD(357.52910918) + SEC_TO_RAD(129596581.0481)*t - @@ -492,6 +471,5 @@ SEC_TO_RAD(0.000136)*t3 - SEC_TO_RAD(0.00001149)*t4; - - // L - Omega + // L - Omega (L = Mean Longitude of the Moon, Omega = F4) F[2] = DEG_TO_RAD(93.27209062) + SEC_TO_RAD(1739527262.8478)*t - @@ -499,6 +477,5 @@ SEC_TO_RAD(0.001037)*t3 + SEC_TO_RAD(0.00000417)*t4; - - // Mean Elongation of the Moon from the Sun + // Mean Elongation of the Moon from the Sun (D) F[3] = DEG_TO_RAD(297.85019547) + SEC_TO_RAD(1602961601.2090)*t - @@ -506,6 +483,5 @@ SEC_TO_RAD(0.006593)*t3 - SEC_TO_RAD(0.00003169)*t4; - - // Mean Longitude of the Ascending Node of the Moon + // Mean Longitude of the Ascending Node of the Moon (Omega) F[4] = DEG_TO_RAD(125.04455501) - SEC_TO_RAD(6962890.5431)*t + @@ -513,5 +489,5 @@ SEC_TO_RAD(0.007702)*t3 - SEC_TO_RAD(0.0000593)*t4; - + //F5-F13 are the mean longitudes of the planets F[5] = 4.402608842 + 2608.7903141574*t; F[6] = 3.176146697 + 1021.3285546211*t; @@ -528,5 +504,4 @@ eocInitialized = p_psEOCInit(); if(!eocInitialized) { - // XXX: Move error message. psError(PS_ERR_UNKNOWN, false, "Could not initialize EOC tables -- check data files."); @@ -534,10 +509,9 @@ } } - // calculate the polynomial portion first double X = psPolynomial1DEval(xPoly,t); double Y = psPolynomial1DEval(yPoly,t); double S = psPolynomial1DEval(sPoly,t); - + //Units from the table & coefficients are in micro-arcseconds so convert to radians. X = SEC_TO_RAD(X * 1e-6); Y = SEC_TO_RAD(Y * 1e-6); @@ -545,35 +519,17 @@ // now calculate the non-poly portion from the tables - psF64* cols[17]; + psF64* cols[17]; //Used to store all of the table information from tab5.2?.dat for (int lcv = 0; lcv < 17; lcv++) { cols[lcv] = ((psVector*)(xTable->data[lcv]))->data.F64; } + //Get the number of rows in the table and loop through the non-poly contributions. int numRows = ((psVector*)(xTable->data[0]))->n; for (int lcv = 0; lcv < numRows; lcv++) { - // arg = sum w_(i,j,k)*F_k double arg = 0.0; + //Get the argument- from the table and corresponding F-value. Convert to radians. for (int k = 0; k < 14; k++) { arg += cols[k+3][lcv]*F[k]; } - double tj = pow(t,cols[0][lcv]); - double as = cols[1][lcv]; - double ac = cols[2][lcv]; - // printf("as-x, ac-x, = %.13g, %.13g\n", as, ac); - as = SEC_TO_RAD(as) * 1e-6; - ac = SEC_TO_RAD(ac) * 1e-6; - // X += as*tj*sin(arg) + ac*tj*cos(arg); - X += (as*tj*sin(arg) + ac*cos(arg)) * tj; - } - - for (int lcv = 0; lcv < 17; lcv++) { - cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64; - } - numRows = ((psVector*)(yTable->data[0]))->n; - for (int lcv = 0; lcv < numRows; lcv++) { - // arg = sum w_(i,j,k)*F_k - double arg = 0.0; - for (int k = 0; k < 14; k++) { - arg += cols[k+3][lcv]*F[k]; - } + //The order of t for each part is specified by a column in the tab.dat files. double tj = pow(t,cols[0][lcv]); double as = cols[1][lcv]; @@ -581,15 +537,12 @@ as = SEC_TO_RAD(as) * 1e-6; ac = SEC_TO_RAD(ac) * 1e-6; - - // Y += as*tj*sin(arg) + ac*tj*cos(arg); - Y += (as*tj*sin(arg) + ac*cos(arg)) * tj; - } - + X += (as*tj*sin(arg) + ac*cos(arg)) * tj; + } + //Do the same procedure from the previous 15 lines for Y. for (int lcv = 0; lcv < 17; lcv++) { - cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64; - } - numRows = ((psVector*)(sTable->data[0]))->n; + cols[lcv] = ((psVector*)(yTable->data[lcv]))->data.F64; + } + numRows = ((psVector*)(yTable->data[0]))->n; for (int lcv = 0; lcv < numRows; lcv++) { - // arg = sum w_(i,j,k)*F_k double arg = 0.0; for (int k = 0; k < 14; k++) { @@ -601,13 +554,28 @@ as = SEC_TO_RAD(as) * 1e-6; ac = SEC_TO_RAD(ac) * 1e-6; - - // S += as*tj*sin(arg) + ac*tj*cos(arg); + Y += (as*tj*sin(arg) + ac*cos(arg)) * tj; + } + //Again for S. + for (int lcv = 0; lcv < 17; lcv++) { + cols[lcv] = ((psVector*)(sTable->data[lcv]))->data.F64; + } + numRows = ((psVector*)(sTable->data[0]))->n; + for (int lcv = 0; lcv < numRows; lcv++) { + double arg = 0.0; + for (int k = 0; k < 14; k++) { + arg += cols[k+3][lcv]*F[k]; + } + double tj = pow(t,cols[0][lcv]); + double as = cols[1][lcv]; + double ac = cols[2][lcv]; + as = SEC_TO_RAD(as) * 1e-6; + ac = SEC_TO_RAD(ac) * 1e-6; S += (as*tj*sin(arg) + ac*cos(arg)) * tj; } - // now, the tables for S actually gives S + XY/2, so let's get the real S now + //the tables for S actually gives S + XY/2, so let's get the real S now. (from ADD) S -= X*Y/2.0; - // psEarthPole* pole = psAlloc(sizeof(psEarthPole)); + //Create the output psEarthPole and set the corresponding component values. psEarthPole* pole = psEarthPoleAlloc(); pole->x = X; @@ -621,12 +589,8 @@ psTimeBulletin bulletin) { + // Check for null parameter or invalid bulletin PS_ASSERT_PTR_NON_NULL(time,NULL); PS_ASSERT_INT_WITHIN_RANGE(bulletin, PS_IERS_A, PS_IERS_B, NULL); - - psEarthPole *out = psEarthPoleAlloc(); - out->x = 0.0; - out->y = 0.0; - out->s = 0.0; - + //Convert the input time to MJD. If NAN is returned, return NULL for the function. double MJD = psTimeToMJD(time); if (MJD == NAN) { @@ -635,16 +599,11 @@ return NULL; } - /* - if (MJD < 41684.0 || MJD > 53334.0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Invalid time input. Date, %lf, is out of range (41684-53334)\n", MJD); - return out; - } - */ + //Allocate space for the psEarthPole output. + psEarthPole *out = psEarthPoleAlloc(); + // Check if EOC data loaded if(!eocInitialized) { eocInitialized = p_psEOCInit(); if(!eocInitialized) { - // XXX: Move error message. psError(PS_ERR_UNKNOWN, false, "Could not initialize EOC tables -- check data files."); @@ -653,4 +612,5 @@ } + //Load the data table and store in the corresponding, newly-allocated vectors. psF64* cols[5]; for (int colNum = 0; colNum < 5; colNum++) { @@ -658,23 +618,4 @@ } int numRows = ((psVector*)(iersTable->data[0]))->n; - /* - for (int rowNum = 0; rowNum < numRows; rowNum++) { - if ( (MJD - cols[0][rowNum]) < 1.0 ) { - if (bulletin == PS_IERS_A) { - out->x = cols[1][rowNum]; - out->y = cols[2][rowNum]; - out->x = SEC_TO_RAD(out->x) * 1e-3; - out->y = SEC_TO_RAD(out->y) * 1e-3; - rowNum = numRows; - } else { - out->x = cols[3][rowNum]; - out->y = cols[4][rowNum]; - out->x = SEC_TO_RAD(out->x) * 1e-3; - out->y = SEC_TO_RAD(out->y) * 1e-3; - rowNum = numRows; - } - } - } - */ psVector *X = psVectorAlloc(numRows, PS_TYPE_F64); psVector *Y = psVectorAlloc(numRows, PS_TYPE_F64); @@ -694,4 +635,5 @@ } + //The following uses lagrange interpolation to calculate the corrections to X & Y. double xOut = 0.0; double yOut = 0.0; @@ -724,4 +666,5 @@ } } + //Convert the values from milli-arcsecond to radian. out->x = SEC_TO_RAD(xOut) * 1e-3; out->y = SEC_TO_RAD(yOut) * 1e-3; @@ -736,7 +679,8 @@ { int i; - psSphereRot *new = (psSphereRot*)psAlloc(sizeof(psSphereRot)); - // psMemSetDeallocator(new, (psFreeFunc)sphereRotFree); - //Convert rotation matrix to quaternions + psSphereRot *new = psSphereRotAlloc(0.0, 0.0, 0.0); + new->q3 = 0.0; + + //Convert rotation matrix to quaternions. Formula directly from ADD. double diag_sum[3]; int maxi; @@ -783,4 +727,5 @@ psSphereRot* psSphereRot_CEOtoGCRS(const psEarthPole *pole) { + // Check for null parameter PS_ASSERT_PTR_NON_NULL(pole,NULL); double A[3][3]; @@ -788,17 +733,5 @@ //Setup the rotation matrix and scalar value, a, as outlined by the ADD - //XXX: Used formula for rotation matrix D from mathworld for z-axis rotation double a = 1.0 / (1.0 + sqrt(1.0 - (pole->x*pole->x + pole->y*pole->y) ) ); - /* A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s); - A[0][1] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s); - A[0][2] = pole->x*cos(pole->s) + pole->y*sin(pole->s); - A[1][0] = -(1.0 - a*pole->x*pole->x)*sin(pole->s) - a*pole->x*pole->y*cos(pole->s); - A[1][1] = a*pole->x*pole->y*sin(pole->s) + (1.0 - a*pole->y*pole->y)*cos(pole->s); - A[1][2] = -pole->x*sin(pole->s) + pole->y*cos(pole->s); - A[2][0] = -pole->x; - A[2][1] = -pole->y; - A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y); - */ - A[0][0] = (1.0 - a*pole->x*pole->x)*cos(pole->s) - a*pole->x*pole->y*sin(pole->s); A[1][0] = -a*pole->x*pole->y*cos(pole->s) + (1.0 - a*pole->y*pole->y)*sin(pole->s); @@ -811,28 +744,5 @@ A[2][2] = 1.0 - a*(pole->x*pole->x + pole->y*pole->y); - double x, y, s; - x = pole->x; - y = pole->y; - s = pole->s; - /* A[0][0] = (1.0 - a*x*x)*cos(s) + a*x*y*sin(s); - A[0][1] = (1.0 - a*x*x)*sin(s) - a*x*y*cos(s); - A[0][2] = x; - A[1][0] = -a*x*y*cos(s) - (1.0 - a*y*y)*sin(s); - A[1][1] = -a*x*y*sin(s) + (1.0 - a*y*y)*cos(s); - A[1][2] = y; - A[2][0] = -x*cos(s) + y*sin(s); - A[2][1] = -x*sin(s) - y*cos(s); - A[2][2] = 1.0 - a * (x*x + y*y); - */ - /* A[0][0] = (1.0 - a*x*x)*cos(s) + a*x*y*sin(s); - A[1][0] = (1.0 - a*x*x)*sin(s) - a*x*y*cos(s); - A[2][0] = x; - A[0][1] = -a*x*y*cos(s) - (1.0 - a*y*y)*sin(s); - A[1][1] = -a*x*y*sin(s) + (1.0 - a*y*y)*cos(s); - A[2][1] = y; - A[0][2] = -x*cos(s) + y*sin(s); - A[1][2] = -x*sin(s) - y*cos(s); - A[2][2] = 1.0 - a * (x*x + y*y); - */ + //Convert the rotation matrix to a psSphereRot (quaternions) out = rotMatrix_To_Quat(A); @@ -843,12 +753,19 @@ psEarthPole *tidalCorr) { + // Check for null parameter PS_ASSERT_PTR_NON_NULL(time,NULL); - psTime *in = psTimeAlloc(time->type); - *in = *time; + //Create a copy of the input time that can be manipulated/changed (input time is const). + psTime *in = p_psTimeCopy(time); if (in->type != PS_TIME_UT1) { in = psTimeConvert(in, PS_TIME_UT1); } + //Check if tidal corrections should be included. + //If so, make sure values are positive and in the correct range. if (tidalCorr != NULL && tidalCorr->s != 0.0) { int nsec = in->nsec + (int)(tidalCorr->s * 1e9); + if (nsec > 1e9) { + nsec += -1e9; + in->sec += 1; + } if (nsec < 0) { in->sec += -1; @@ -858,4 +775,5 @@ } } + //Calculate the Julian Date from the input time in UT1 format. double T = psTimeToJD(in); if (T == NAN) { @@ -865,4 +783,5 @@ } T += -2451545.0; + //Formula for theta comes directly from the ADD. Create output psSphereRot from theta. double theta = 2.0 * M_PI * (0.7790572732640 + 1.00273781191135448 * T); psSphereRot *out = psSphereRotAlloc(theta, 0.0, 0.0); @@ -1003,4 +922,7 @@ static double DMOD(double x, double y) { + //Internal function for calculating the remainder of a double quotient. + //used often in the algorithm for psEOC_PolarTideCorr which is the Ray model of + //Simon et. al. from its fortran implementation. double value = x - y * trunc(x/y); return value; @@ -1011,6 +933,4 @@ // Check for null parameter PS_ASSERT_PTR_NON_NULL(time, NULL); - psEarthPole *out = psEarthPoleAlloc(); - // Convert psTime to MJD double MJD = psTimeToMJD(time); @@ -1020,14 +940,12 @@ return NULL; } - - // Calculate number of Julian centuries since 2000 - double RJD = MJD; - - //Formula comes from fortran reference - //DMOD in fortran ref. = double remainder -> x - y * trunc(x/y) + psEarthPole *out = psEarthPoleAlloc(); + + //Formula comes from fortran reference of the Ray model of Simon et. al. double T, L, LPRIME, CAPF, CAPD, OMEGA, THETA, CORX, CORY, CORZ; double ARG1, ARG2, ARG3, ARG4, ARG5, ARG6, ARG7, ARG8; double T2, T3, T4; - T = (RJD - 51544.5) / 36525.0; + // Calculate number of Julian centuries since 2000 + T = (MJD - 51544.5) / 36525.0; T2 = T*T; T3 = T*T*T; @@ -1086,5 +1004,4 @@ CORX = SEC_TO_RAD(CORX); CORY = SEC_TO_RAD(CORY); - // CORZ = SEC_TO_RAD(CORZ); out->x = CORX; @@ -1099,10 +1016,4 @@ // Check for null parameter PS_ASSERT_PTR_NON_NULL(time, NULL); - /* if (time->type == PS_TIME_UT1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "Invalid time input. Time cannot be of type UT1.\n"); - return NULL; - } - */ // Convert psTime to MJD double MJD = psTimeToMJD(time); @@ -1120,4 +1031,5 @@ double t4 = t*t*t*t; + // The following formulae are from the ADD (& s5.8 of IERS Technical Note 32): double F[5]; // Mean Anomaly of the Moon @@ -1135,5 +1047,5 @@ SEC_TO_RAD(0.00001149)*t4; - // L − Omega + // L − Omega (L = Mean Longitude of the Moon, Omega = F4) F[2] = DEG_TO_RAD(93.27209062) + SEC_TO_RAD(1739527262.8478)*t - @@ -1252,28 +1164,4 @@ double X = 0.0; double Y = 0.0; - // double arg = 0.0; - //This is from eqn 131 in the ADD - Note: pj_tj isn't included the first time. - //XXX: The xp_sin, yp_cos, etc. may need to be multiplied by pow(t,i) here? adding now... - // double tj = 0.0; - - // calculate the polynomial portion first - the pj * t^j (poly coeff's) - // X = psPolynomial1DEval(xPoly,t); - // Y = psPolynomial1DEval(yPoly,t); - // X = SEC_TO_RAD(X * 1e-6); - // Y = SEC_TO_RAD(Y * 1e-6); - /* for (int i = 0; i < 10; i++) { - double arg = 0.0; - double as = 0.0; - double ac = 0.0; - arg = w_l[i]*F[0] + w_l_p[i]*F[1] + w_F[i]*F[2] + w_D[i]*F[3] + w_Omega[i]*F[4]; - tj = 1.0; - as = xp_sin[i] * 1e-6; - ac = xp_cos[i] * 1e-6; - X += (as*tj*sin(arg) + ac*cos(arg)) * tj; - as = yp_sin[i] * 1e-6; - ac = yp_cos[i] * 1e-6; - Y += (as*tj*sin(arg) + ac*cos(arg)) * tj; - } - */ //Implementation adapted from PM_GRAVI in interp.f from hpiers.obspm.fr/eop-pc/models/interp.f @@ -1306,9 +1194,6 @@ ag = RAD_TO_SEC(ag); ag = DMOD(ag, 2.0*M_PI); - // ag = SEC_TO_RAD(ag); X += xp_sin[j] * SEC_TO_RAD(sin(ag)) + xp_cos[j] * SEC_TO_RAD(cos(ag)); Y += yp_sin[j] * SEC_TO_RAD(sin(ag)) + yp_cos[j] * SEC_TO_RAD(cos(ag)); - // X += xp_sin[j] * sin(ag) + xp_cos[j] * cos(ag); - // Y += yp_sin[j] * sin(ag) + yp_cos[j] * cos(ag); } @@ -1316,4 +1201,5 @@ pole->x = X; pole->y = Y; + //The value of s is simply: s = -4.7e-5 * t as specified by the ADD and IERS 32. pole->s = -SEC_TO_RAD(4.7e-5) * t; @@ -1323,4 +1209,5 @@ psSphereRot* psSphereRot_ITRStoTEO(const psEarthPole* motion) { + // Check for null parameter PS_ASSERT_PTR_NON_NULL(motion,NULL); double A[3][3]; @@ -1333,79 +1220,6 @@ s = -s; - - //Setup Rotation Matrix for transformation (x,y,z rotation) - //XXX: May need to be (z,y,x as in Mathworld?) - /* A[0][0] = cos(motion->x)*cos(-motion->s); - A[0][1] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s); - A[0][2] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s); - A[1][0] = -cos(motion->x)*sin(-motion->s); - A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s); - A[1][2] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s); - A[2][0] = sin(motion->x); - A[2][1] = -sin(motion->y)*cos(motion->x); - A[2][2] = cos(motion->y)*cos(motion->x); - */ - /* - A[0][0] = cos(motion->x)*cos(-motion->s); - A[1][0] = sin(motion->y)*sin(motion->x)*cos(-motion->s) + cos(motion->y)*sin(-motion->s); - A[2][0] = -cos(motion->y)*sin(motion->x)*cos(-motion->s) + sin(motion->y)*sin(-motion->s); - A[0][1] = -cos(motion->x)*sin(-motion->s); - A[1][1] = -sin(motion->y)*sin(motion->x)*sin(-motion->s) + cos(motion->y)*cos(-motion->s); - A[2][1] = cos(motion->y)*sin(motion->x)*sin(-motion->s) + sin(motion->y)*cos(-motion->s); - A[0][2] = sin(motion->x); - A[1][2] = -sin(motion->y)*cos(motion->x); - A[2][2] = cos(motion->y)*cos(motion->x); - */ - /* - A[0][0] = cos(s)*cos(x); - A[0][1] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y); - A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); - A[1][0] = -sin(s)*cos(x); - A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); - A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); - A[2][0] = sin(x); - A[2][1] = -cos(x)*sin(y); - A[2][2] = cos(x)*cos(y); - */ - /* - A[0][0] = cos(s)*cos(x); - A[1][0] = sin(s)*sin(y) + cos(s)*sin(x)*sin(y); - A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); - A[0][1] = -sin(s)*cos(x); - A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); - A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); - A[0][2] = sin(x); - A[1][2] = -cos(x)*sin(y); - A[2][2] = cos(x)*cos(y); - */ - - /* - psSphereRot r,p,t; - - r.q0=sin(y/2.0); - r.q1=0; - r.q2=0; - r.q3=cos(y/2.0); - - p.q0=0; - p.q1=sin(x/2.0); - p.q2=0; - p.q3=cos(x/2.0); - - t.q0=0; - t.q1=0; - t.q2=sin(s/2.0); - t.q3=cos(s/2.0); - - // calculate t*s*r. - psSphereRot* temp = psSphereRotCombine(NULL,&t,&p); - out = psSphereRotCombine(NULL, temp, &r); - psFree(temp); - - return out; - */ - - //Newest trial - mult. y matrix * x * z - + //Setup Rotation matrix. Rotation is constructed by rotation about the X-axis by y, + // about the Y-axis by x, and about the Z-axis by s' (where s' = -s). A[0][0] = cos(x)*cos(s); A[1][0] = cos(x)*sin(s); @@ -1418,38 +1232,4 @@ A[2][2] = cos(x)*cos(y); - /* - A[0][0] = cos(x)*cos(s); - A[0][1] = cos(x)*sin(s); - A[0][2] = -sin(x); - A[1][0] = sin(x)*sin(y)*cos(s) - cos(y)*sin(s); - A[1][1] = sin(x)*sin(y)*sin(s) + cos(y)*cos(s); - A[1][2] = cos(x)*sin(y); - A[2][0] = sin(x)*cos(y)*cos(s) + sin(y)*sin(s); - A[2][1] = sin(x)*cos(y)*sin(s) - sin(y)*cos(s); - A[2][2] = cos(x)*cos(y); - */ - //New trial - mult z * y * x - /* A[0][0] = cos(s)*cos(x); - A[1][0] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y); - A[2][0] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); - A[0][1] = -sin(s)*cos(x); - A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); - A[2][1] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); - A[0][2] = sin(x); - A[1][2] = -cos(x)*sin(y); - A[2][2] = cos(x)*cos(y); - */ - /* - A[0][0] = cos(s)*cos(x); - A[0][1] = sin(s)*cos(y) + cos(s)*sin(x)*sin(y); - A[0][2] = sin(s)*sin(y) - cos(s)*sin(x)*cos(y); - A[1][0] = -sin(s)*cos(x); - A[1][1] = cos(s)*cos(y) - sin(s)*sin(x)*sin(y); - A[1][2] = cos(s)*sin(y) + sin(s)*sin(x)*cos(y); - A[2][0] = sin(x); - A[2][1] = -cos(x)*sin(y); - A[2][2] = cos(x)*cos(y); - */ - //Convert rotation matrix to quaternions out = rotMatrix_To_Quat(A); @@ -1484,4 +1264,7 @@ psTime *from = NULL; psTime *to = NULL; + if (fromMJD > toMJD) { + psWarning("From time is later than to time in psSpherePrecess.\n"); + } if (mode == PS_PRECESS_ROUGH) { @@ -1525,5 +1308,6 @@ to = p_psTimeCopy(toTime); } - + //Calculate the earthpoles and quaternions corresponding to each time (from, to). + //Combine the quaternions to produce the output psSphereRot. psEarthPole *fromEP = psEOC_PrecessionModel(from); psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); @@ -1557,10 +1341,11 @@ to = p_psTimeCopy(toTime); } - + //Calculate the earthpoles and quaternions corresponding to each time (from, to). + //Add in the precession corrections from IERS bulletin A. + //Combine the quaternions to produce the output psSphereRot. psEarthPole *fromEP = psEOC_PrecessionModel(from); psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_A); fromEP->x += fromCorr->x; fromEP->y += fromCorr->y; - fromEP->s += fromCorr->s; psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); psEarthPole *toEP = psEOC_PrecessionModel(to); @@ -1568,5 +1353,4 @@ toEP->x += toCorr->x; toEP->y += toCorr->y; - toEP->s += toCorr->s; psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP); psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot); @@ -1600,9 +1384,11 @@ } + //Calculate the earthpoles and quaternions corresponding to each time (from, to). + //Add in the precession corrections from IERS bulletin B. + //Combine the quaternions to produce the output psSphereRot. psEarthPole *fromEP = psEOC_PrecessionModel(from); psEarthPole *fromCorr = psEOC_PrecessionCorr(from, PS_IERS_B); fromEP->x += fromCorr->x; fromEP->y += fromCorr->y; - fromEP->s += fromCorr->s; psSphereRot *fromRot = psSphereRot_CEOtoGCRS(fromEP); psEarthPole *toEP = psEOC_PrecessionModel(to); @@ -1610,5 +1396,4 @@ toEP->x += toCorr->x; toEP->y += toCorr->y; - toEP->s += toCorr->s; psSphereRot *toRot = psSphereRot_CEOtoGCRS(toEP); psSphereRot *fromConj = psSphereRotConjugate(NULL, fromRot); @@ -1625,8 +1410,3 @@ return out; } - // Apply transform to coordinates - // psSphere *out = psSphereRotApply(NULL, tmpST, coords); - // if (out->r < -0.0001) { - // out->r += 2.0 * M_PI; - // } -} +} Index: trunk/psLib/src/astro/psEarthOrientation.h =================================================================== --- trunk/psLib/src/astro/psEarthOrientation.h (revision 6409) +++ trunk/psLib/src/astro/psEarthOrientation.h (revision 6425) @@ -9,6 +9,6 @@ * @author Robert Daniel DeSonia, MHPCC * -* @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-01-23 20:04:31 $ +* @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-02-14 00:09:27 $ * * Copyright 2005 Maui High Performance Computing Center, University of Hawaii @@ -131,5 +131,8 @@ /** Calculates the rotation of the Earth about the CIP. * - * @return psSphereRot*: sphereical rotation + * If tidalCorr is non-NULL, use the S-component to provide tidal corrections to the + * UT1 time. If tidalCorr is NULL, no corrections are made & this step is skipped. + * + * @return psSphereRot*: spherical rotation of the Earth about the CIP. */ psSphereRot *psSphereRot_TEOtoCEO( @@ -187,5 +190,4 @@ */ psSphereRot* psSpherePrecess( - // psSphere *coords, ///< coordinates (modified in-place) const psTime *fromTime, ///< equinox of coords input const psTime *toTime, ///< equinox of coords output Index: trunk/psLib/src/astro/psSphereOps.h =================================================================== --- trunk/psLib/src/astro/psSphereOps.h (revision 6409) +++ trunk/psLib/src/astro/psSphereOps.h (revision 6425) @@ -8,6 +8,6 @@ * @author David Robbins, MHPCC * - * @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-23 23:52:15 $ + * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-14 00:09:27 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -20,5 +20,4 @@ /// @{ -//#include "psCoord.h" #include "psTime.h"