OceanLoading.cpp

Go to the documentation of this file.

#pragma ident "$Id: OceanLoading.cpp 1346 2008-08-05 14:59:55Z architest $"

//============================================================================
//
//  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, 2008
//
//============================================================================


#include "OceanLoading.hpp"


namespace gpstk
{


      /* Returns the effect of ocean tides loading (meters) at the given
       * station and epoch, in the Up-East-North (UEN) reference frame.
       *
       * @param name  Station name (case is NOT relevant).
       * @param time  Epoch to look up
       *
       * @return a Triple with the ocean tidas loading effect, in meters and
       * in the UEN reference frame.
       *
       * @throw InvalidRequest If the request can not be completed for any
       * reason, this is thrown. The text may have additional information
       * about the reason the request failed.
       */
   Triple OceanLoading::getOceanLoading( const string& name,
                                         const DayTime& t )
      throw(InvalidRequest)
   {

      const int NUM_COMPONENTS = 3;
      const int NUM_HARMONICS = 11;

      Matrix<double> harmonics(6,11,0.0);

         // Get harmonics data from file
      harmonics = blqData.getTideHarmonics(name);

      Vector<double> arguments(11,0.0);

         // Compute arguments
      arguments = getArg(t);

      Triple oLoading;

      for(int i=0; i<NUM_COMPONENTS;  i++)
      {

         double temp(0.0);
         for(int k=0; k<NUM_HARMONICS; k++)
         {

            temp += harmonics(i,k) * 
                    std::cos( arguments(k) - harmonics( (i+3),k)*DEG_TO_RAD );

         }

            // This Triple is in Up, West, South reference frame
         oLoading[i] = temp;

      }  // End of 'for(int i=0; i<NUM_COMPONENTS;  i++)'

         // Let's change Triple to Up, East, North [UEN] reference frame
      oLoading[1] = -oLoading[1];
      oLoading[2] = -oLoading[2];

      return oLoading;

   }  // End of method 'OceanLoading::getOceanLoading()'



      /* Sets the name of BLQ file containing ocean tides harmonics data.
       *
       * @param name      Name of BLQ tides harmonics data file.
       */
   OceanLoading& OceanLoading::setFilename(const string& name)
   {

      fileData = name;

      blqData.open(fileData);

      return (*this);

   }  // End of meters 'OceanLoading::setFilename()'



      /* Compute the value of the corresponding astronomical arguments,
       * in radians. This routine is based on IERS routine ARG.f.
       *
       * @param time      Epoch of interest
       *
       * @return A Vector<double> of 11 elements with the corresponding
       * astronomical arguments to be used in ocean loading model.
       */
   Vector<double> OceanLoading::getArg(const DayTime& time)
   {

      const int NUM_HARMONICS = 11;

         // Let's store some important values
      Vector<double> sig(NUM_HARMONICS,0.0);
      sig(0) = 1.40519e-4;
      sig(1) = 1.45444e-4;
      sig(2) = 1.37880e-4;
      sig(3) = 1.45842e-4;
      sig(4) = 0.72921e-4;
      sig(5) = 0.67598e-4;
      sig(6) = 0.72523e-4;
      sig(7) = 0.64959e-4;
      sig(8) = 0.053234e-4;
      sig(9) = 0.026392e-4;
      sig(10)= 0.003982e-4;

      Matrix<double> angfac(4,NUM_HARMONICS,0.0);
      angfac(0,0) =  2.0; angfac(1,0) = -2.0;
         angfac(2,0) =  0.0; angfac(3,0) =  0.0;
      angfac(0,1) =  0.0; angfac(1,1) =  0.0;
         angfac(2,1) =  0.0; angfac(3,1) =  0.0;
      angfac(0,2) =  2.0; angfac(1,2) = -3.0;
         angfac(2,2) =  1.0; angfac(3,2) =  0.0;
      angfac(0,3) =  2.0; angfac(1,3) =  0.0;
         angfac(2,3) =  0.0; angfac(3,3) =  0.0;
      angfac(0,4) =  1.0; angfac(1,4) =  0.0;
         angfac(2,4) =  0.0; angfac(3,4) =  0.25;
      angfac(0,5) =  1.0; angfac(1,5) = -2.0;
         angfac(2,5) =  0.0; angfac(3,5) = -0.25;
      angfac(0,6) = -1.0; angfac(1,6) =  0.0;
         angfac(2,6) =  0.0; angfac(3,6) = -0.25;
      angfac(0,7) =  1.0; angfac(1,7) = -3.0;
         angfac(2,7) =  1.0; angfac(3,7) = -0.25;
      angfac(0,8) =  0.0; angfac(1,8) =  2.0;
         angfac(2,8) =  0.0; angfac(3,8) =  0.0;
      angfac(0,9) =  0.0; angfac(1,9) =  1.0;
         angfac(2,9) = -1.0; angfac(3,9) =  0.0;
      angfac(0,10)=  2.0; angfac(1,10)=  0.0;
         angfac(2,10)=  0.0; angfac(3,10)=  0.0;


      Vector<double> arguments(NUM_HARMONICS,0.0);

         // Get day of year
      short year(time.year());

         // Fractional part of day, in seconds
      double fday(time.DOYsecond());

         // Compute time
      double d(time.DOY()+365.0*(year-1975.0)+floor((year-1973.0)/4.0));
      double t((27392.500528+1.000000035*d)/36525.0);

         // Mean longitude of Sun at beginning of day
      double H0((279.69668+(36000.768930485+3.03e-4*t)*t)*DEG_TO_RAD);

         // Mean longitude of Moon at beginning of day
      double S0( (((1.9e-6*t - 0.001133)*t +
                    481267.88314137)*t + 270.434358)*DEG_TO_RAD );

         // Mean longitude of lunar perigee at beginning of day
      double P0( (((-1.2e-5*t - 0.010325)*t +
                    4069.0340329577)*t + 334.329653)*DEG_TO_RAD );

      for(int k=0; k<NUM_HARMONICS; k++)
      {

         double temp( sig(k)*fday + angfac(0,k)*H0 + 
                      angfac(1,k)*S0 + angfac(2,k)*P0 + angfac(3,k)*TWO_PI );

         arguments(k) = fmod(temp,TWO_PI);

         if (arguments(k) < 0.0)
         {
            arguments(k) = arguments(k) + TWO_PI;
         }

      }  // End of 'for(int k=0; k<NUM_HARMONICS; k++)'

      return arguments;

   }  // End of method 'OceanLoading::getArg()'



}  // End of namespace gpstk

---