// Part of BNC, a utility for retrieving decoding and
// converting GNSS data streams from NTRIP broadcasters.
//
// Copyright (C) 2007
// German Federal Agency for Cartography and Geodesy (BKG)
// http://www.bkg.bund.de
// Czech Technical University Prague, Department of Geodesy
// http://www.fsv.cvut.cz
//
// Email: euref-ip@bkg.bund.de
//
// This program is free software; you can redistribute it and/or
// modify it under the terms of the GNU General Public License
// as published by the Free Software Foundation, version 2.
//
// This program 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 General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
/* -------------------------------------------------------------------------
* BKG NTRIP Client
* -------------------------------------------------------------------------
*
* Class: t_astro, t_tides, t_tropo
*
* Purpose: Observation model
*
* Author: L. Mervart
*
* Created: 29-Jul-2014
*
* Changes:
*
* -----------------------------------------------------------------------*/
#include <cmath>
#include "pppModel.h"
using namespace BNC_PPP;
using namespace std;
Matrix t_astro::rotX(double Angle) {
const double C = cos(Angle);
const double S = sin(Angle);
Matrix UU(3, 3);
UU[0][0] = 1.0;
UU[0][1] = 0.0;
UU[0][2] = 0.0;
UU[1][0] = 0.0;
UU[1][1] = +C;
UU[1][2] = +S;
UU[2][0] = 0.0;
UU[2][1] = -S;
UU[2][2] = +C;
return UU;
}
Matrix t_astro::rotY(double Angle) {
const double C = cos(Angle);
const double S = sin(Angle);
Matrix UU(3, 3);
UU[0][0] = +C;
UU[0][1] = 0.0;
UU[0][2] = -S;
UU[1][0] = 0.0;
UU[1][1] = 1.0;
UU[1][2] = 0.0;
UU[2][0] = +S;
UU[2][1] = 0.0;
UU[2][2] = +C;
return UU;
}
Matrix t_astro::rotZ(double Angle) {
const double C = cos(Angle);
const double S = sin(Angle);
Matrix UU(3, 3);
UU[0][0] = +C;
UU[0][1] = +S;
UU[0][2] = 0.0;
UU[1][0] = -S;
UU[1][1] = +C;
UU[1][2] = 0.0;
UU[2][0] = 0.0;
UU[2][1] = 0.0;
UU[2][2] = 1.0;
return UU;
}
// Greenwich Mean Sidereal Time
///////////////////////////////////////////////////////////////////////////
double t_astro::GMST(double Mjd_UT1) {
const double Secs = 86400.0;
double Mjd_0 = floor(Mjd_UT1);
double UT1 = Secs * (Mjd_UT1 - Mjd_0);
double T_0 = (Mjd_0 - MJD_J2000) / 36525.0;
double T = (Mjd_UT1 - MJD_J2000) / 36525.0;
double gmst = 24110.54841 + 8640184.812866 * T_0 + 1.002737909350795 * UT1
+ (0.093104 - 6.2e-6 * T) * T * T;
return 2.0 * M_PI * Frac(gmst / Secs);
}
// Nutation Matrix
///////////////////////////////////////////////////////////////////////////
Matrix t_astro::NutMatrix(double Mjd_TT) {
const double T = (Mjd_TT - MJD_J2000) / 36525.0;
double ls = 2.0 * M_PI * Frac(0.993133 + 99.997306 * T);
double D = 2.0 * M_PI * Frac(0.827362 + 1236.853087 * T);
double F = 2.0 * M_PI * Frac(0.259089 + 1342.227826 * T);
double N = 2.0 * M_PI * Frac(0.347346 - 5.372447 * T);
double dpsi = (-17.200 * sin(N) - 1.319 * sin(2 * (F - D + N))
- 0.227 * sin(2 * (F + N))
+ 0.206 * sin(2 * N) + 0.143 * sin(ls)) / RHO_SEC;
double deps = (+9.203 * cos(N) + 0.574 * cos(2 * (F - D + N))
+ 0.098 * cos(2 * (F + N))
- 0.090 * cos(2 * N)) / RHO_SEC;
double eps = 0.4090928 - 2.2696E-4 * T;
return rotX(-eps - deps) * rotZ(-dpsi) * rotX(+eps);
}
// Precession Matrix
///////////////////////////////////////////////////////////////////////////
Matrix t_astro::PrecMatrix(double Mjd_1, double Mjd_2) {
const double T = (Mjd_1 - MJD_J2000) / 36525.0;
const double dT = (Mjd_2 - Mjd_1) / 36525.0;
double zeta = ((2306.2181 + (1.39656 - 0.000139 * T) * T) +
((0.30188 - 0.000344 * T) + 0.017998 * dT) * dT) * dT / RHO_SEC;
double z = zeta
+ ((0.79280 + 0.000411 * T) + 0.000205 * dT) * dT * dT / RHO_SEC;
double theta = ((2004.3109 - (0.85330 + 0.000217 * T) * T) -
((0.42665 + 0.000217 * T) + 0.041833 * dT) * dT) * dT / RHO_SEC;
return rotZ(-z) * rotY(theta) * rotZ(-zeta);
}
// Sun's position
///////////////////////////////////////////////////////////////////////////
ColumnVector t_astro::Sun(double Mjd_TT) {
const double eps = 23.43929111 / RHO_DEG;
const double T = (Mjd_TT - MJD_J2000) / 36525.0;
double M = 2.0 * M_PI * Frac(0.9931267 + 99.9973583 * T);
double L = 2.0 * M_PI * Frac(0.7859444 + M / 2.0 / M_PI +
(6892.0 * sin(M) + 72.0 * sin(2.0 * M)) / 1296.0e3);
double r = 149.619e9 - 2.499e9 * cos(M) - 0.021e9 * cos(2 * M);
ColumnVector r_Sun(3);
r_Sun << r * cos(L) << r * sin(L) << 0.0;
r_Sun = rotX(-eps) * r_Sun;
return rotZ(GMST(Mjd_TT))
* NutMatrix(Mjd_TT)
* PrecMatrix(MJD_J2000, Mjd_TT)
* r_Sun;
}
// Moon's position
///////////////////////////////////////////////////////////////////////////
ColumnVector t_astro::Moon(double Mjd_TT) {
const double eps = 23.43929111 / RHO_DEG;
const double T = (Mjd_TT - MJD_J2000) / 36525.0;
double L_0 = Frac(0.606433 + 1336.851344 * T);
double l = 2.0 * M_PI * Frac(0.374897 + 1325.552410 * T);
double lp = 2.0 * M_PI * Frac(0.993133 + 99.997361 * T);
double D = 2.0 * M_PI * Frac(0.827361 + 1236.853086 * T);
double F = 2.0 * M_PI * Frac(0.259086 + 1342.227825 * T);
double dL = +22640 * sin(l) - 4586 * sin(l - 2 * D) + 2370 * sin(2 * D)
+ 769 * sin(2 * l)
- 668 * sin(lp) - 412 * sin(2 * F) - 212 * sin(2 * l - 2 * D)
- 206 * sin(l + lp - 2 * D)
+ 192 * sin(l + 2 * D) - 165 * sin(lp - 2 * D) - 125 * sin(D)
- 110 * sin(l + lp)
+ 148 * sin(l - lp) - 55 * sin(2 * F - 2 * D);
double L = 2.0 * M_PI * Frac(L_0 + dL / 1296.0e3);
double S = F + (dL + 412 * sin(2 * F) + 541 * sin(lp)) / RHO_SEC;
double h = F - 2 * D;
double N = -526 * sin(h) + 44 * sin(l + h) - 31 * sin(-l + h)
- 23 * sin(lp + h)
+ 11 * sin(-lp + h) - 25 * sin(-2 * l + F) + 21 * sin(-l + F);
double B = (18520.0 * sin(S) + N) / RHO_SEC;
double cosB = cos(B);
double R = 385000e3 - 20905e3 * cos(l) - 3699e3 * cos(2 * D - l)
- 2956e3 * cos(2 * D)
- 570e3 * cos(2 * l) + 246e3 * cos(2 * l - 2 * D)
- 205e3 * cos(lp - 2 * D)
- 171e3 * cos(l + 2 * D) - 152e3 * cos(l + lp - 2 * D);
ColumnVector r_Moon(3);
r_Moon << R * cos(L) * cosB << R * sin(L) * cosB << R * sin(B);
r_Moon = rotX(-eps) * r_Moon;
return rotZ(GMST(Mjd_TT))
* NutMatrix(Mjd_TT)
* PrecMatrix(MJD_J2000, Mjd_TT)
* r_Moon;
}
// Tidal Correction
////////////////////////////////////////////////////////////////////////////
ColumnVector t_tides::earth(const bncTime& time, const ColumnVector& xyz) {
if (time.undef()) {
ColumnVector dX(3);
dX = 0.0;
return dX;
}
double Mjd = time.mjd() + time.daysec() / 86400.0;
if (Mjd != _lastMjd) {
_lastMjd = Mjd;
_xSun = t_astro::Sun(Mjd);
_rSun = sqrt(DotProduct(_xSun, _xSun));
_xSun /= _rSun;
_xMoon = t_astro::Moon(Mjd);
_rMoon = sqrt(DotProduct(_xMoon, _xMoon));
_xMoon /= _rMoon;
}
double rRec = sqrt(DotProduct(xyz, xyz));
ColumnVector xyzUnit = xyz / rRec;
// Love's Numbers
// --------------
const double H2 = 0.6078;
const double L2 = 0.0847;
// Tidal Displacement
// ------------------
double scSun = DotProduct(xyzUnit, _xSun);
double scMoon = DotProduct(xyzUnit, _xMoon);
double p2Sun = 3.0 * (H2 / 2.0 - L2) * scSun * scSun - H2 / 2.0;
double p2Moon = 3.0 * (H2 / 2.0 - L2) * scMoon * scMoon - H2 / 2.0;
double x2Sun = 3.0 * L2 * scSun;
double x2Moon = 3.0 * L2 * scMoon;
const double gmWGS = 398.6005e12;
const double gms = 1.3271250e20;
const double gmm = 4.9027890e12;
double facSun = gms / gmWGS *
(rRec * rRec * rRec * rRec) / (_rSun * _rSun * _rSun);
double facMoon = gmm / gmWGS *
(rRec * rRec * rRec * rRec) / (_rMoon * _rMoon * _rMoon);
ColumnVector dX = facSun * (x2Sun * _xSun + p2Sun * xyzUnit)
+ facMoon * (x2Moon * _xMoon + p2Moon * xyzUnit);
return dX;
}
// Constructor
///////////////////////////////////////////////////////////////////////////
t_tides::t_tides() {
_lastMjd = 0.0;
_rSun = 0.0;
_rMoon = 0.0;
newBlqData = 0;
}
t_tides::~t_tides() {
if (newBlqData) {
delete newBlqData;
}
QMapIterator<QString, t_blqData*> it(blqMap);
while (it.hasNext()) {
it.next();
delete it.value();
}
blqMap.clear();
}
t_irc t_tides::readBlqFile(const char* fileName) {
QFile inFile(fileName);
inFile.open(QIODevice::ReadOnly | QIODevice::Text);
QTextStream in(&inFile);
int row = 0;
QString site = QString();
while (!in.atEnd()) {
QString line = in.readLine();
// skip empty lines and comments
if (line.indexOf("$$") != -1) {
continue;
}
line = line.trimmed();
QTextStream inLine(line.toLatin1(), QIODevice::ReadOnly);
switch (row) {
case 0:
site = line;
site = site.toUpper();
newBlqData = new t_blqData;
newBlqData->amplitudes.ReSize(3, 11);
newBlqData->phases.ReSize(3, 11);
break;
case 1:
case 2:
case 3:
for (int ii = 0; ii < 11; ii++) {
inLine >> newBlqData->amplitudes[row - 1][ii];
}
break;
case 4:
case 5:
for (int ii = 0; ii < 11; ii++) {
inLine >> newBlqData->phases[row - 4][ii];
}
break;
case 6:
for (int ii = 0; ii < 11; ii++) {
inLine >> newBlqData->phases[row - 4][ii];
}
if (newBlqData && !site.isEmpty()) {
blqMap[site] = newBlqData;
site = QString();
newBlqData = 0;
}
row = -1;
break;
}
row++;
}
inFile.close();
return success;
}
ColumnVector t_tides::ocean(const bncTime& time, const ColumnVector& xyz,
const std::string& station) {
ColumnVector dX(3); dX = 0.0;
if (time.undef()) {
return dX;
}
QString stationQ = station.c_str();
if (blqMap.find(stationQ) == blqMap.end()) {
return dX;
}
t_blqData* blqSet = blqMap[stationQ]; //printBlqSet(station, blqSet);
// angular argument: see arg2.f from IERS Conventions software collection
double speed[11] = {1.40519e-4, 1.45444e-4, 1.3788e-4, 1.45842e-4, 7.2921e-5,
6.7598e-5, 7.2523e-5, 6.4959e-5, 5.3234e-6, 2.6392e-6, 3.982e-7};
double angfac[4][11];
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] = .25;
angfac[0][5] = 1.0;
angfac[1][5] =-2.0;
angfac[2][5] = 0.0;
angfac[3][5] =-.25;
angfac[0][6] =-1.0;
angfac[1][6] = 0.0;
angfac[2][6] = 0.0;
angfac[3][6] =-.25;
angfac[0][7] = 1.0;
angfac[1][7] =-3.0;
angfac[2][7] = 1.0;
angfac[3][7] =-.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;
double twopi = 6.283185307179586476925287e0;
double dtr = 0.0174532925199;
// fractional part of the day in seconds
unsigned int year, month, day;
time.civil_date(year, month, day);
int iyear = year - 2000;
QDateTime datTim = QDateTime::fromString(QString::fromStdString(time.datestr()), Qt::ISODate);
int doy = datTim.date().dayOfYear();
double fday = time.daysec();
int icapd = doy + 365 * (iyear - 75) + ((iyear - 73) / 4);
double capt = (icapd * 1.000000035 + 27392.500528) / 36525.0;
// mean longitude of the sun at the beginning of the day
double h0 = (279.69668e0 + (36000.768930485e0 + 3.03e-4 * capt) * capt) * dtr;
// mean longitude of moon at the beginning of the day
double s0 = (((1.9e-6 * capt - .001133e0) * capt + 481267.88314137e0) * capt + 270.434358e0) * dtr;
// mean longitude of lunar perigee at the beginning of the day
double p0 = (((-1.2e-5 * capt - .010325e0) * capt + 4069.0340329577e0) * capt + 334.329653e0) * dtr;
// tidal angle arguments
double angle[11];
for (int k = 0; k < 11; ++k) {
angle[k] = speed[k] * fday
+ angfac[0][k] * h0
+ angfac[1][k] * s0
+ angfac[2][k] * p0
+ angfac[3][k] * twopi;
angle[k] = fmod(angle[k], twopi);
if (angle[k] < 0.0) {
angle[k] += twopi;
}
}
// displacement by 11 constituents
ColumnVector rwsSta(3); rwsSta = 0.0; // radial, west, south
for (int rr = 0; rr < 3; rr++) {
for (int cc = 0; cc < 11; cc++) {
rwsSta[rr] += blqSet->amplitudes[rr][cc] * cos((angle[cc] - (blqSet->phases[rr][cc]/RHO_DEG)));
}
}
// neu2xyz
ColumnVector dneu(3); // neu
dneu[0] = -rwsSta[2];
dneu[1] = -rwsSta[1];
dneu[2] = rwsSta[0];
double recEll[3]; xyz2ell(xyz.data(), recEll) ;
neu2xyz(recEll, dneu.data(), dX.data());
return dX;
}
// Print
////////////////////////////////////////////////////////////////////////////
void t_tides::printAllBlqSets() const {
QMapIterator<QString, t_blqData*> it(blqMap);
while (it.hasNext()) {
it.next();
t_blqData* blq = it.value();
QString site = it.key();
cout << site.toStdString().c_str() << "\n===============\n";
for (int rr = 0; rr < 3; rr++) {
for (int cc = 0; cc < 11; cc++) {
cout << blq->amplitudes[rr][cc] << " ";
}
cout << endl;
}
for (int rr = 0; rr < 3; rr++) {
for (int cc = 0; cc < 11; cc++) {
cout << blq->phases[rr][cc] << " ";
}
cout << endl;
}
}
}
// Print
////////////////////////////////////////////////////////////////////////////
void t_tides::printBlqSet(const std::string& station, t_blqData* blq) {
cout << station << endl;
for (int rr = 0; rr < 3; rr++) {
for (int cc = 0; cc < 11; cc++) {
cout << blq->amplitudes[rr][cc] << " ";
}
cout << endl;
}
for (int rr = 0; rr < 3; rr++) {
for (int cc = 0; cc < 11; cc++) {
cout << blq->phases[rr][cc] << " ";
}
cout << endl;
}
}
// Constructor
///////////////////////////////////////////////////////////////////////////
t_windUp::t_windUp() {
for (unsigned ii = 0; ii <= t_prn::MAXPRN; ii++) {
sumWind[ii] = 0.0;
lastEtime[ii] = 0.0;
}
}
// Phase Wind-Up Correction
///////////////////////////////////////////////////////////////////////////
double t_windUp::value(const bncTime& etime, const ColumnVector& rRec,
t_prn prn, const ColumnVector& rSat, bool ssr,
double yaw, const ColumnVector& vSat) {
if (etime.mjddec() != lastEtime[prn.toInt()]) {
// Unit Vector GPS Satellite --> Receiver
// --------------------------------------
ColumnVector rho = rRec - rSat;
rho /= rho.NormFrobenius();
// GPS Satellite unit Vectors sz, sy, sx
// -------------------------------------
ColumnVector sHlp;
if (!ssr) {
sHlp = t_astro::Sun(etime.mjddec());
}
else {
ColumnVector Omega(3);
Omega[0] = 0.0;
Omega[1] = 0.0;
Omega[2] = t_CST::omega;
sHlp = vSat + crossproduct(Omega, rSat);
}
sHlp /= sHlp.NormFrobenius();
ColumnVector sz = -rSat / rSat.NormFrobenius();
ColumnVector sy = crossproduct(sz, sHlp);
ColumnVector sx = crossproduct(sy, sz);
if (ssr) {
// Yaw angle consideration
Matrix SXYZ(3, 3);
SXYZ.Column(1) = sx;
SXYZ.Column(2) = sy;
SXYZ.Column(3) = sz;
SXYZ = DotProduct(t_astro::rotZ(yaw), SXYZ);
sx = SXYZ.Column(1);
sy = SXYZ.Column(2);
sz = SXYZ.Column(3);
}
// Effective Dipole of the GPS Satellite Antenna
// ---------------------------------------------
ColumnVector dipSat = sx - rho * DotProduct(rho, sx)
- crossproduct(rho, sy);
// Receiver unit Vectors rx, ry
// ----------------------------
ColumnVector rx(3);
ColumnVector ry(3);
double recEll[3];
xyz2ell(rRec.data(), recEll);
double neu[3];
neu[0] = 1.0;
neu[1] = 0.0;
neu[2] = 0.0;
neu2xyz(recEll, neu, rx.data());
neu[0] = 0.0;
neu[1] = -1.0;
neu[2] = 0.0;
neu2xyz(recEll, neu, ry.data());
// Effective Dipole of the Receiver Antenna
// ----------------------------------------
ColumnVector dipRec = rx - rho * DotProduct(rho, rx)
+ crossproduct(rho, ry);
// Resulting Effect
// ----------------
double alpha = DotProduct(dipSat, dipRec)
/ (dipSat.NormFrobenius() * dipRec.NormFrobenius());
if (alpha > 1.0)
alpha = 1.0;
if (alpha < -1.0)
alpha = -1.0;
double dphi = acos(alpha) / 2.0 / M_PI; // in cycles
if (DotProduct(rho, crossproduct(dipSat, dipRec)) < 0.0) {
dphi = -dphi;
}
if (lastEtime[prn.toInt()] == 0.0) {
sumWind[prn.toInt()] = dphi;
}
else {
sumWind[prn.toInt()] = nint(sumWind[prn.toInt()] - dphi) + dphi;
}
lastEtime[prn.toInt()] = etime.mjddec();
}
return sumWind[prn.toInt()];
}
// Tropospheric Model (Saastamoinen)
////////////////////////////////////////////////////////////////////////////
double t_tropo::delay_saast(const ColumnVector& xyz, double Ele) {
Tracer tracer("bncModel::delay_saast");
if (xyz[0] == 0.0 && xyz[1] == 0.0 && xyz[2] == 0.0) {
return 0.0;
}
double ell[3];
xyz2ell(xyz.data(), ell);
double height = ell[2];
// Prevent pp from causing segmentation fault (Loukis)
if (height > 40000.0 ) {
return 0.000000001;
}
double pp = 1013.25 * pow(1.0 - 2.26e-5 * height, 5.225);
double TT = 18.0 - height * 0.0065 + 273.15;
double hh = 50.0 * exp(-6.396e-4 * height);
double ee = hh / 100.0
* exp(-37.2465 + 0.213166 * TT - 0.000256908 * TT * TT);
double h_km = height / 1000.0;
if (h_km < 0.0)
h_km = 0.0;
if (h_km > 5.0)
h_km = 5.0;
int ii = int(h_km + 1);
if (ii > 5)
ii = 5;
double href = ii - 1;
double bCor[6];
bCor[0] = 1.156;
bCor[1] = 1.006;
bCor[2] = 0.874;
bCor[3] = 0.757;
bCor[4] = 0.654;
bCor[5] = 0.563;
double BB = bCor[ii - 1] + (bCor[ii] - bCor[ii - 1]) * (h_km - href);
double zen = M_PI / 2.0 - Ele;
return (0.002277 / cos(zen))
* (pp + ((1255.0 / TT) + 0.05) * ee - BB * (tan(zen) * tan(zen)));
}
// Constructor
///////////////////////////////////////////////////////////////////////////
t_iono::t_iono() {
_psiPP = _phiPP = _lambdaPP = _lonS = 0.0;
}
t_iono::~t_iono() {
}
double t_iono::stec(const t_vTec* vTec, double signalPropagationTime,
const ColumnVector& rSat, const bncTime& epochTime,
const ColumnVector& xyzSta) {
// Latitude, longitude, height are defined with respect to a spherical earth model
// -------------------------------------------------------------------------------
ColumnVector geocSta(3);
if (xyz2geoc(xyzSta.data(), geocSta.data()) != success) {
return 0.0;
}
// satellite position rotated to the epoch of signal reception
// -----------------------------------------------------------
ColumnVector xyzSat(3);
double omegaZ = t_CST::omega * signalPropagationTime;
xyzSat[0] = rSat[0] * cos(omegaZ) + rSat[1] * sin(omegaZ);
xyzSat[1] = rSat[1] * cos(omegaZ) - rSat[0] * sin(omegaZ);
xyzSat[2] = rSat[2];
// elevation and azimuth with respect to a spherical earth model
// -------------------------------------------------------------
ColumnVector rhoV = xyzSat - xyzSta;
double rho = rhoV.NormFrobenius();
ColumnVector neu(3);
xyz2neu(geocSta.data(), rhoV.data(), neu.data());
double sphEle = acos(sqrt(neu[0] * neu[0] + neu[1] * neu[1]) / rho);
if (neu[2] < 0) {
sphEle *= -1.0;
}
double sphAzi = atan2(neu[1], neu[0]);
double epoch = fmod(epochTime.gpssec(), 86400.0);
double stec = 0.0;
for (unsigned ii = 0; ii < vTec->_layers.size(); ii++) {
piercePoint(vTec->_layers[ii]._height, epoch, geocSta.data(), sphEle,
sphAzi);
double vtec = vtecSingleLayerContribution(vTec->_layers[ii]);
stec += vtec * sin(sphEle + _psiPP);
}
return stec;
}
double t_iono::vtecSingleLayerContribution(const t_vTecLayer& vTecLayer) {
double vtec = 0.0;
int N = vTecLayer._C.Nrows() - 1;
int M = vTecLayer._C.Ncols() - 1;
double fac;
for (int n = 0; n <= N; n++) {
for (int m = 0; m <= min(n, M); m++) {
double pnm = associatedLegendreFunction(n, m, sin(_phiPP));
double a = factorial(n - m);
double b = factorial(n + m);
if (m == 0) {
fac = sqrt(2.0 * n + 1);
}
else {
fac = sqrt(2.0 * (2.0 * n + 1) * a / b);
}
pnm *= fac;
double Cnm_mlambda = vTecLayer._C[n][m] * cos(m * _lonS);
double Snm_mlambda = vTecLayer._S[n][m] * sin(m * _lonS);
vtec += (Snm_mlambda + Cnm_mlambda) * pnm;
}
}
if (vtec < 0.0) {
vtec = 0.0;
}
return vtec;
}
void t_iono::piercePoint(double layerHeight, double epoch,
const double* geocSta,
double sphEle, double sphAzi) {
double q = (t_CST::rgeoc + geocSta[2]) / (t_CST::rgeoc + layerHeight);
_psiPP = M_PI / 2 - sphEle - asin(q * cos(sphEle));
_phiPP = asin(
sin(geocSta[0]) * cos(_psiPP)
+ cos(geocSta[0]) * sin(_psiPP) * cos(sphAzi));
if (((geocSta[0] * 180.0 / M_PI > 0)
&& (tan(_psiPP) * cos(sphAzi) > tan(M_PI / 2 - geocSta[0])))
||
((geocSta[0] * 180.0 / M_PI < 0)
&& (-(tan(_psiPP) * cos(sphAzi)) > tan(M_PI / 2 + geocSta[0])))) {
_lambdaPP = geocSta[1] + M_PI
- asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
}
else {
_lambdaPP = geocSta[1] + asin((sin(_psiPP) * sin(sphAzi) / cos(_phiPP)));
}
_lonS = fmod((_lambdaPP + (epoch - 50400) * M_PI / 43200), 2 * M_PI);
return;
}