1 //==============================================================================
2 //
3 //  This file is part of GPSTk, the GPS Toolkit.
4 //
5 //  The GPSTk is free software; you can redistribute it and/or modify
6 //  it under the terms of the GNU Lesser General Public License as published
7 //  by the Free Software Foundation; either version 3.0 of the License, or
8 //  any later version.
9 //
10 //  The GPSTk is distributed in the hope that it will be useful,
11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 //  GNU Lesser General Public License for more details.
14 //
15 //  You should have received a copy of the GNU Lesser General Public
16 //  License along with GPSTk; if not, write to the Free Software Foundation,
17 //  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 //  This software was developed by Applied Research Laboratories at the
20 //  University of Texas at Austin.
21 //  Copyright 2004-2020, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24 
25 //==============================================================================
26 //
27 //  This software was developed by Applied Research Laboratories at the
28 //  University of Texas at Austin, under contract to an agency or agencies
29 //  within the U.S. Department of Defense. The U.S. Government retains all
30 //  rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 //  Pursuant to DoD Directive 523024
33 //
34 //  DISTRIBUTION STATEMENT A: This software has been approved for public
35 //                            release, distribution is unlimited.
36 //
37 //==============================================================================
38 
39 /**
40  * @file SolidTides.cpp
41  * Computes the effect of solid Earth tides on a given position and epoch.
42  */
43 
44 #include "SolidTides.hpp"
45 
46 
47 namespace gpstk
48 {
49 
50 
51       // Love numbers
52    const double SolidTides::H_LOVE(0.609), SolidTides::L_LOVE(0.0852);
53 
54       // Phase lag. Assumed as zero here because no phase lag has been
55       // detected so far.
56       // const double SolidTides::PH_LAG(0.0);
57 
58 
59       /* Returns the effect of solid Earth tides (meters) at the given
60        * position and epoch, in the Up-East-Down (UEN) reference frame.
61        *
62        * @param[in] t Epoch to look up
63        * @param[in] p Position of interest
64        *
65        * @return a Triple with the solid tidal effect, in meters and in
66        * the UED reference frame.
67        *
68        * @throw InvalidRequest If the request can not be completed for any
69        * reason, this is thrown. The text may have additional information
70        * as to why the request failed.
71        */
getSolidTide(const CommonTime & t,const Position & p) const72    Triple SolidTides::getSolidTide(const CommonTime& t,
73                                    const Position& p) const
74    {
75 
76          // We will store here the results
77       Triple res;
78 
79          // Objects to compute Sun and Moon positions
80       SunPosition  sunPosition;
81       MoonPosition moonPosition;
82 
83 
84       try
85       {
86 
87             // Variables to hold Sun and Moon positions
88          Triple sunPos(sunPosition.getPosition(t));
89          Triple moonPos(moonPosition.getPosition(t));
90 
91 
92             // Compute the factors for the Sun
93          double rpRs( p.X()*sunPos.theArray[0] +
94                       p.Y()*sunPos.theArray[1] +
95                       p.Z()*sunPos.theArray[2]);
96 
97          double Rs2(sunPos.theArray[0]*sunPos.theArray[0] +
98                     sunPos.theArray[1]*sunPos.theArray[1] +
99                     sunPos.theArray[2]*sunPos.theArray[2]);
100 
101          double rp2( p.X()*p.X() + p.Y()*p.Y() + p.Z()*p.Z() );
102 
103          double xy2p( p.X()*p.X() + p.Y()*p.Y() );
104          double sqxy2p( std::sqrt(xy2p) );
105 
106          double sqRs2(std::sqrt(Rs2));
107 
108          double fac_s( 3.0*MU_SUN*rp2/(sqRs2*sqRs2*sqRs2*sqRs2*sqRs2) );
109 
110          double g1sun( fac_s*(rpRs*rpRs/2.0 - rp2*Rs2/6.0) );
111 
112          double g2sun( fac_s * rpRs * (sunPos.theArray[1]*p.X() -
113                        sunPos.theArray[0]*p.Y()) * std::sqrt(rp2)/sqxy2p );
114 
115          double g3sun( fac_s * rpRs * ( sqxy2p* sunPos.theArray[2] -
116                        p.Z()/sqxy2p * (p.X()*sunPos.theArray[0] +
117                        p.Y()*sunPos.theArray[1]) ) );
118 
119 
120             // Compute the factors for the Moon
121          double rpRm( p.X()*moonPos.theArray[0] +
122                       p.Y()*moonPos.theArray[1] +
123                       p.Z()*moonPos.theArray[2]);
124 
125          double Rm2(moonPos.theArray[0]*moonPos.theArray[0] +
126                     moonPos.theArray[1]*moonPos.theArray[1] +
127                     moonPos.theArray[2]*moonPos.theArray[2]);
128 
129          double sqRm2(std::sqrt(Rm2));
130 
131          double fac_m( 3.0*MU_MOON*rp2/(sqRm2*sqRm2*sqRm2*sqRm2*sqRm2) );
132 
133          double g1moon( fac_m*(rpRm*rpRm/2.0 - rp2*Rm2/6.0) );
134 
135          double g2moon( fac_m * rpRm * (moonPos.theArray[1]*p.X() -
136                         moonPos.theArray[0]*p.Y()) * std::sqrt(rp2)/sqxy2p );
137 
138          double g3moon( fac_m * rpRm * ( sqxy2p* moonPos.theArray[2] -
139                         p.Z()/sqxy2p * (p.X()*moonPos.theArray[0] +
140                         p.Y()*moonPos.theArray[1]) ) );
141 
142             // Effects due to the Sun
143          double delta_sun1(H_LOVE*g1sun);
144          double delta_sun2(L_LOVE*g2sun);
145          double delta_sun3(L_LOVE*g3sun);
146 
147             // Effects due to the Moon
148          double delta_moon1(H_LOVE*g1moon);
149          double delta_moon2(L_LOVE*g2moon);
150          double delta_moon3(L_LOVE*g3moon);
151 
152             // Combined effect
153          res.theArray[0] = delta_sun1 + delta_moon1;
154          res.theArray[1] = delta_sun2 + delta_moon2;
155          res.theArray[2] = delta_sun3 + delta_moon3;
156 
157       } // End of try block
158       catch(InvalidRequest& ir)
159       {
160          GPSTK_RETHROW(ir);
161       }
162 
163       return res;
164 
165    } // End SolidTides::getSolidTide
166 
167 
168 } // end namespace gpstk
169