2
$\begingroup$

I am attempting to model GPS Satellite positions on the globe for a set of ephemerides. I have a verified set of ECI ( http://en.wikipedia.org/wiki/Earth-centered_inertial ) XYZ Coordinates and a verified set of ECEF ( http://en.wikipedia.org/wiki/ECEF ) XYZ coordinates.

For example ( m ):
ECEF X 4076514.50
ECEF Y 14673598.00
ECEF Z -21793674.00

ECI X 8004604.50
ECI Y 12956032.00
ECI Z -21793674.00

Sometimes the signs of the ECI X and Y coordinates are opposite of the signs for the ECEF X and Y coordinates.

When I plot these locations on a globe and, obviously they do not match unless you rotate the earth by some number of radians. I can see that if you rotate the earth, the ECI coordinates follow a groundtrack that is consistent with the ECEF groundtrack, but with my current code, it is offset by different varying amounts, depending on the time of day I plot it.

Here is my code:

// The radians earth rotates in 1 second
double earthRotation = 0.000072921151467;

// Do something different depending on the switch position
if ( self.shouldAnimate == true ) {

    // The seconds since we started rotating
    double newRotationMultiplier = [self.eciDate timeIntervalSinceDate:self.lastEciDate];

    // Seconds           // Rad/Sec        // Multiplier
    self.timeSinceOpenGlStarted = newRotationMultiplier * earthRotation;
    matrixToRotate = GLKMatrix4Rotate(matrixToRotate, self.timeSinceOpenGlStarted, 0.0, 0.0, 1.0);

    // Latch it so that when we stop rotating, it doesn't revert
    self.effect.transform.modelviewMatrix = matrixToRotate;

}

How do I determine the radians to rotate the earth to show an ECI positions for a time to be the same location on the rotated globe ( in terms of latitude and longitude )? That is, how do I calculate for any UTC time in a day, how much to rotate the earth about the Z ( "Up" ) axis so that the ECI coordinates agree with the ECEF coordinates and the satellites appear to be in the same position over my 3D globe?

UPDATE:

I got it to work - thanks to the answer below! Here is the code:

    // Get the NSDate components
    NSCalendar *currentCalendar = [NSCalendar currentCalendar];
    [currentCalendar setTimeZone:[NSTimeZone timeZoneWithAbbreviation:@"UTC"]];
    NSDateComponents *components = [currentCalendar components:NSCalendarUnitDay | NSCalendarUnitMonth | NSCalendarUnitYear | NSHourCalendarUnit | NSMinuteCalendarUnit | NSSecondCalendarUnit fromDate:self.eciDate];
    [components setTimeZone:[NSTimeZone timeZoneWithAbbreviation:@"UTC"]];
    NSInteger hour = [components hour];
    NSInteger minute = [components minute];
    NSInteger second = [components second];

    // Get UTC in terms of decimal hours
    float utc = hour + (float)minute/60.0f + (float)second/3600.0f; // UTC is expressed in hours as a decimal number

    //NSLog(@"UTC: %f H %f M %f S %f",utc,(float)hour,(float)minute/60.0f ,(float)second/3600.0f);

    // Get the seconds since the last vernal equinox
    double secondsSinceVernalEquinox = [self.eciDate timeIntervalSinceDate:appDelegate.lastVernalEquinox];

    // Do the math
    // http://physics.stackexchange.com/questions/98466/radians-to-rotate-earth-to-match-eci-lat-lon-with-ecef-lat-lon
    double d = secondsSinceVernalEquinox  / ( 60.0 * 60.0 * 24.0 ); // Give decimal days
    double p = 365.242187; // Length of tropical year
    double result = ( M_PI * 2 / 24 ) * ( utc - 12.0 ) + ( 2 * M_PI * d  / p );

    //NSLog(@"Radians = %f",result);

    // Seconds           // Rad/Sec        // Multiplier
    self.timeSinceOpenGlStarted = result;
    matrixToRotate = GLKMatrix4Rotate(matrixToRotate, self.timeSinceOpenGlStarted, 0.0, 0.0, 1.0);

    // Latch it so that when we stop rotating, it doesn't revert
    self.effect.transform.modelviewMatrix = matrixToRotate;
$\endgroup$
  • $\begingroup$ My first answer was not entirely accurate (I think it's about 1.8 degrees off), because $d$ should be measured from the perihelion instead of the vernal equinox, plus an extra constant. However, I found a more accurate and simpler formula; see my updated answer. $\endgroup$ – Pulsar Feb 10 '14 at 4:03
4
$\begingroup$

Disclaimer: I updated my previous answer as I found a more accurate formula.

The rotation angle $\theta$ is the angle between the vernal point and the meridian at Greenwich. This also corresponds with the Sidereal Time at Greenwich, converted to radians.

enter image description here

For a given UTC time and a given date, the corresponding Greenwich (Mean) Sidereal Time (in hours and decimal parts of an hour) can be calculated as $$ \text{GMST} = 18.697\,374\,558 + 24.065\,709\,824\,419\,d, $$ where $d$ is the interval, in UT1 days including any fraction of a day, since 2000 January 1, at 12h UT (source: US Naval Observatory). Of course, this must be reduced to the range 0h to 24h. The rotation angle is then simply $$ \theta = \frac{2\pi}{24}\text{GMST}. $$

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.