OceanLoading.cpp

Go to the documentation of this file.

//==============================================================================
//
//  This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
//
//  The GNSSTk 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 3.0 of the License, or
//  any later version.
//
//  The GNSSTk 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 GNSSTk; if not, write to the Free Software Foundation,
//  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin.
//  Copyright 2004-2022, The Board of Regents of The University of Texas System
//
//==============================================================================
 
//==============================================================================
//
//  This software was developed by Applied Research Laboratories at the
//  University of Texas at Austin, under contract to an agency or agencies
//  within the U.S. Department of Defense. The U.S. Government retains all
//  rights to use, duplicate, distribute, disclose, or release this software.
//
//  Pursuant to DoD Directive 523024
//
//  DISTRIBUTION STATEMENT A: This software has been approved for public
//                            release, distribution is unlimited.
//
//==============================================================================
 
#include "OceanLoading.hpp"
#include "YDSTime.hpp"
 
using namespace std;
 
namespace gnsstk
{
 
 
      /* 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 CommonTime& t )
   {
 
      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 CommonTime& 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
      int year(static_cast<YDSTime>(time).year);
 
         // Fractional part of day, in seconds
      double fday(static_cast<YDSTime>(time).sod);
 
         // Compute time
      double d(static_cast<YDSTime>(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 gnsstk