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 IonoModel.cpp 41 * Implementation of the ICD-GPS-200 Ionosphere model. 42 */ 43 44 #include <math.h> 45 #include "GNSSconstants.hpp" 46 #include "IonoModel.hpp" 47 #include "YDSTime.hpp" 48 #include "GNSSconstants.hpp" 49 50 namespace gpstk 51 { IonoModel(const double a[4],const double b[4])52 IonoModel::IonoModel(const double a[4], const double b[4]) throw() 53 { 54 setModel(a, b); 55 } 56 IonoModel(const EngAlmanac & engalm)57 IonoModel::IonoModel(const EngAlmanac& engalm) 58 throw() 59 { 60 try 61 { 62 engalm.getIon(alpha, beta); 63 valid = true; 64 } 65 catch(InvalidRequest& e) 66 { 67 valid = false; 68 } 69 } 70 71 setModel(const double a[4],const double b[4])72 void IonoModel::setModel(const double a[4], const double b[4]) throw() 73 { 74 for (int n = 0; n < 4; n++) 75 { 76 alpha[n] = a[n]; 77 beta[n] = b[n]; 78 } 79 valid = true; 80 } 81 82 getCorrection(const CommonTime & time,const Position & rxgeo,double svel,double svaz,CarrierBand band) const83 double IonoModel::getCorrection(const CommonTime& time, 84 const Position& rxgeo, 85 double svel, 86 double svaz, 87 CarrierBand band) const 88 { 89 90 if (!valid) 91 { 92 InvalidIonoModel e("Alpha and beta parameters invalid."); 93 GPSTK_THROW(e); 94 } 95 96 // all angle units are in semi-circles (radians / TWO_PI) 97 // Note: math functions (cos, sin, etc.) require arguments in 98 // radians so all semi-circles must be multiplied by TWO_PI 99 100 double azRad = svaz * DEG_TO_RAD; 101 double svE = svel / 180.0; 102 103 double phi_u = rxgeo.getGeodeticLatitude() / 180.0; 104 double lambda_u = rxgeo.getLongitude() / 180.0; 105 106 double psi = (0.0137 / (svE + 0.11)) - 0.022; 107 108 double phi_i = phi_u + psi * std::cos(azRad); 109 if (phi_i > 0.416) 110 phi_i = 0.416; 111 if (phi_i < -0.416) 112 phi_i = -0.416; 113 114 double lambda_i = lambda_u + psi * ::sin(azRad) / ::cos(phi_i*PI); 115 116 double phi_m = phi_i + 0.064 * ::cos((lambda_i - 1.617)*PI); 117 118 double iAMP = 0.0; 119 double iPER = 0.0; 120 iAMP = alpha[0]+phi_m*(alpha[1]+phi_m*(alpha[2]+phi_m*alpha[3])); 121 iPER = beta[0]+phi_m*( beta[1]+phi_m*( beta[2]+phi_m* beta[3])); 122 123 if (iAMP < 0.0) 124 iAMP = 0.0; 125 if (iPER < 72000.0) 126 iPER = 72000.0; 127 128 double t = 43200.0 * lambda_i + YDSTime(time).sod; 129 if (t >= 86400.0) 130 t -= 86400.0; 131 if (t < 0) 132 t += 86400.0; 133 134 double x = TWO_PI * (t - 50400.0) / iPER; // x is in radians 135 136 double iF = 1.0 + 16.0 * (0.53 - svE)*(0.53 - svE)*(0.53 - svE); 137 138 double t_iono = 0.0; 139 if (fabs(x) < 1.57) 140 t_iono = iF * (5.0e-9 + iAMP * (1 + x*x * (-0.5 + x*x/24.0))); 141 else 142 t_iono = iF * 5.0e-9; 143 144 // Correction factor for GPS band; see ICD-GPS-200 20.3.3.3.3.2. 145 if (band == CarrierBand::L2) 146 { 147 t_iono *= GAMMA_GPS_12; // GAMMA_GPS = (fL1 / fL2)^2 148 } 149 else if (band == CarrierBand::L5) 150 { 151 t_iono *= GAMMA_GPS_15; // GAMMA_GPS = (fL1 / fL5)^2 152 } 153 else if (band != CarrierBand::L1) 154 { 155 throw InvalidIonoModel("Invalid CarrierBand, not one of L1,L2,L5."); 156 } 157 158 double correction = t_iono * C_MPS; // return correction in [m] 159 160 return correction; 161 } 162 operator ==(const IonoModel & right) const163 bool IonoModel::operator==(const IonoModel& right) const 164 throw() 165 { 166 for (int n = 0; n < 4; n++) 167 { 168 if (alpha[n] != right.alpha[n] || beta[n] != right.beta[n]) 169 return false; 170 } 171 return true; 172 } 173 operator !=(const IonoModel & right) const174 bool IonoModel::operator!=(const IonoModel&right) const 175 throw() 176 { 177 return !(operator==(right)); 178 } 179 } 180