AstronomicalFunctions.cpp

Go to the documentation of this file.

#pragma ident "$Id: AstronomicalFunctions.cpp 2944 2011-10-27 08:04:41Z yanweignss $"

//============================================================================
//
//  This file is part of GPSTk, the GPS Toolkit.
//
//  The GPSTk is free software; you can redistribute it and/or modify
//  it under the terms of the GNU Lesser General Public License as published
//  by the Free Software Foundation; either version 2.1 of the License, or
//  any later version.
//
//  The GPSTk is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU Lesser General Public License for more details.
//
//  You should have received a copy of the GNU Lesser General Public
//  License along with GPSTk; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  
//  Dagoberto Salazar - gAGE ( http://www.gage.es ). 2007
//
//============================================================================


#include "AstronomicalFunctions.hpp"


namespace gpstk
{

      /* Function to change from CIS to CTS(ECEF) coordinate system
       * (coordinates in meters)
       * @param posCis    Coordinates in CIS system (in meters).
       * @param t         Epoch
       *
       * @return Triple in CTS(ECEF) coordinate system.
       */
   Triple CIS2CTS(const Triple posCIS,
                  const DayTime& t)
   {

         // Angle of Earth rotation, in radians
      double ts( UTC2SID(t)*TWO_PI/24.0 );

      Triple res;

      res.theArray[0] = ::cos(ts)*posCIS.theArray[0] +
                        ::sin(ts)*posCIS.theArray[1];

      res.theArray[1] = -::sin(ts)*posCIS.theArray[0] +
                         ::cos(ts)*posCIS.theArray[1];

      res.theArray[2] = posCIS.theArray[2];

      return res;
   } // End CIS2CTS()


      /* Function to convert from UTC to sidereal time
       * @param t         Epoch
       *
       * @return sidereal time in hours.
       */
   double UTC2SID(const DayTime& t)
   {

      double y(t.year()-1.0);
      double m(13.0);
      double d(t.DOY());

         // Hours of day (decimal)
      double h(t.secOfDay()/3600.0);

         // Fraction of day
      double frofday (t.secOfDay()/86400.0);

         // Compute Julian Day, including decimals
      double jd(t.JD());

         // Temporal value, in centuries
      double tt( (jd - 2451545.0)/36525.0 );

      double sid( 24110.54841 + tt*( (8640184.812866) + 
                  tt*( (0.093104) - (6.2e-6*tt)) ) );

      sid = sid/3600.0 + h;
      sid = ::fmod(sid,24.0);

      if (sid < 0.0)
      {
         sid+=24.0;
      }

      return sid;
   }

} // end namespace gpstk

---