1 /// @file AtmLoadTides.cpp 2 /// Classes to handle site displacement due to atmospheric loading. 3 4 //============================================================================== 5 // 6 // This file is part of GPSTk, the GPS Toolkit. 7 // 8 // The GPSTk is free software; you can redistribute it and/or modify 9 // it under the terms of the GNU Lesser General Public License as published 10 // by the Free Software Foundation; either version 3.0 of the License, or 11 // any later version. 12 // 13 // The GPSTk is distributed in the hope that it will be useful, 14 // but WITHOUT ANY WARRANTY; without even the implied warranty of 15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 16 // GNU Lesser General Public License for more details. 17 // 18 // You should have received a copy of the GNU Lesser General Public 19 // License along with GPSTk; if not, write to the Free Software Foundation, 20 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA 21 // 22 // This software was developed by Applied Research Laboratories at the 23 // University of Texas at Austin. 24 // Copyright 2004-2020, The Board of Regents of The University of Texas System 25 // 26 //============================================================================== 27 28 //============================================================================== 29 // 30 // This software was developed by Applied Research Laboratories at the 31 // University of Texas at Austin, under contract to an agency or agencies 32 // within the U.S. Department of Defense. The U.S. Government retains all 33 // rights to use, duplicate, distribute, disclose, or release this software. 34 // 35 // Pursuant to DoD Directive 523024 36 // 37 // DISTRIBUTION STATEMENT A: This software has been approved for public 38 // release, distribution is unlimited. 39 // 40 //============================================================================== 41 42 //------------------------------------------------------------------------------------ 43 // system includes 44 #include <iostream> 45 #include <fstream> 46 #include <algorithm> 47 // GPSTk 48 #include "StringUtils.hpp" 49 #include "MiscMath.hpp" 50 #include "GNSSconstants.hpp" 51 // geomatics 52 #include "AtmLoadTides.hpp" 53 //#include "logstream.hpp" // TEMP 54 55 using namespace std; 56 57 namespace gpstk { 58 59 using namespace StringUtils; 60 61 //--------------------------------------------------------------------------------- 62 // Open and read the given file, containing atmospheric loading coefficients, and 63 // initialize this object for the sites names in the input list that match a 64 // name in the file (case sensitive). Return the number of successfully 65 // initialized site names, and remove those sites from the input list. 66 // Atmospheric loading files can be obtained by running the program grdinterp.f 67 // param sites vector<string> On input contains site labels found in the 68 // file, on output contains only sites that were NOT found. 69 // If sites is empty, all sites are loaded. 70 // param filename string Input atmospheric loading file name. 71 // return the number of sites successfully initialized. 72 // throw if the file could not be opened. initializeSites(vector<string> & sites,string filename)73 int AtmLoadTides::initializeSites(vector<string>& sites, string filename) 74 { 75 try { 76 bool allsites = false; 77 if(sites.size() == 0) allsites = true; 78 int i,nwant,nfound; 79 80 ifstream infile(filename.c_str()); 81 if(!infile || !infile.is_open()) { 82 Exception e("File " + filename + " could not be opened."); 83 GPSTK_THROW(e); 84 } 85 86 nwant = sites.size(); 87 nfound = 0; // number of successes 88 bool looking=true; // true if looking for a site name 89 double lat,lon; 90 vector<double> coeff; 91 string site; 92 while(1) { // read the file 93 int count; 94 string line,word; 95 96 // get the next line 97 getline(infile,line); 98 stripTrailing(line,'\r'); 99 stripLeading(line); 100 101 // process line 102 if(!line.empty()) { 103 word = firstWord(line); 104 //LOG(INFO) << "Word is " << word << " and line is " << line; 105 106 if(word == "$$") { // NB ignore header - assume column order, etc. 107 // pick out the lat/lon 108 while(!line.empty()) { 109 word = stripFirstWord(line); 110 if(word == string("station")) { 111 site = stripFirstWord(line); 112 stripTrailing(site,';'); 113 } 114 if(word == string("coord.(long,lat)")) { 115 lon = asDouble(stripFirstWord(line)); 116 lat = asDouble(stripFirstWord(line)); 117 //LOG(INFO) << "Found coords for site " << site << " " 118 // << fixed << setprecision(6) << lon << " " << lat; 119 break; 120 } 121 } 122 } 123 else if(looking && line.length() <= 40) { // site 124 // site name 125 site = line; 126 //LOG(INFO) << "Found site " << site; 127 if(allsites) { 128 //LOG(INFO) << "Push back " << site; 129 looking = false; 130 sites.push_back(site); 131 } 132 else for(i=0; i<sites.size(); i++) { 133 //LOG(INFO) << "Compare " << sites[i]; 134 if(site == sites[i]) { 135 looking = false; 136 break; 137 } 138 } 139 if(!looking) { // found a site 140 count = 0; 141 coeff.clear(); 142 } 143 } 144 else if(!looking) { // not comment and not looking - must be data 145 if(numWords(line) != 4) { 146 Exception e("File " + filename + " is corrupted for site " + site 147 + " - offending line follows\n" + line); 148 GPSTK_THROW(e); 149 } 150 //LOG(INFO) << "Push back line " << line; 151 for(i=0; i<4; i++) 152 coeff.push_back(asDouble(stripFirstWord(line))); 153 count++; 154 if(count == 3) { // success 155 ostringstream oss; 156 oss << fixed; 157 for(i=0; i<coeff.size(); i++) { 158 oss << " " << setprecision(4) << setw(7) << coeff[i]; 159 if((i+1)%4 == 0) oss << "\n"; 160 } 161 //LOG(INFO) << " Found site " << site << " with coefficients:"; 162 //LOG(INFO) << oss.str(); 163 164 // update coeff map 165 coefficientMap[site] = coeff; 166 nfound++; 167 // update position map 168 coeff.clear(); 169 coeff.push_back(lat); 170 coeff.push_back(lon); 171 positionMap[site] = coeff; 172 //LOG(INFO) << "pushback coords for site " << site << " " 173 // << fixed << setprecision(6) << lon << " " << lat; 174 175 // erase a vector element 176 if(!allsites) { 177 vector<string>::iterator pos; 178 pos = find(sites.begin(),sites.end(),site); 179 if(pos != sites.end()) sites.erase(pos); 180 } 181 looking = true; 182 } 183 } 184 185 } // end if line not empty 186 187 if(!allsites && nfound >= nwant) break; 188 189 if(infile.eof() || !infile.good()) break; 190 191 } // end loop over lines in the file 192 193 return nfound; 194 } 195 catch(Exception& e) { GPSTK_RETHROW(e); } 196 catch(exception& e) { 197 Exception E("std except: "+string(e.what())); 198 GPSTK_THROW(E); 199 } 200 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 201 } 202 203 //--------------------------------------------------------------------------------- 204 // Compute the site displacement vector at the given time for the given site. 205 // The site must have been successfully initialized; if not an exception is 206 // thrown. 207 // param site string Input name of the site; must be the same as previously 208 // successfully passed to initializeSites(). 209 // param t EphTime Input time of interest. 210 // param UT1mUTC Difference of UT1 and UTC, a very small correction to t. 211 // return Triple containing the North, East and Up components of the site 212 // displacement in meters. 213 // throw if the site has not been initialized, if the time system is unknown, 214 // if there is corruption in the static arrays, or . computeDisplacement(string site,EphTime time,double UT1mUTC)215 Triple AtmLoadTides::computeDisplacement(string site, EphTime time, double UT1mUTC) 216 { 217 try { 218 219 // get the coefficients for this site 220 if(coefficientMap.find(site) == coefficientMap.end()) 221 GPSTK_THROW(Exception("Site not found in atmospheric loading store")); 222 vector<double> coeff = coefficientMap[site]; 223 224 // compute time argument 225 EphTime ttag(time); 226 ttag.convertSystemTo(TimeSystem::UTC); 227 // ignoring UT1-UTC is probably fine, since this is extremely small 228 ttag += UT1mUTC; 229 double dayfr(ttag.secOfDay()/86400.0); // fraction of day 230 static const double w1(2*PI), w2(4*PI); 231 const double cos1(::cos(w1*dayfr)); 232 const double sin1(::sin(w1*dayfr)); 233 const double cos2(::cos(w2*dayfr)); 234 const double sin2(::sin(w2*dayfr)); 235 236 // NB Displacement is defined positive up, north 237 // and east directions, and is in millimeters 238 239 // total d(t) = d(1)*cos(t*w1) + d(2)*sin(t*w1) 240 // + d(3)*cos(t*w2) + d(4)*sin(t*w2) 241 // If t is fractions of a UT1 day, 242 // then w1=2pi radians/day and w2=4pi radians/day. 243 244 // Column order coss1 sins1 coss2 sins2 245 // Row (coeff) order: RADIAL NS EW 246 247 Triple dc; // dc is NEU, so must RAD,N,E -> NEU 248 dc[2] = coeff[0]*cos1 + coeff[1]*sin1 + coeff[2]*cos2 + coeff[3]*sin2; 249 dc[0] = coeff[4]*cos1 + coeff[5]*sin1 + coeff[6]*cos2 + coeff[7]*sin2; 250 dc[1] = coeff[8]*cos1 + coeff[9]*sin1 + coeff[10]*cos2 + coeff[11]*sin2; 251 // convert to meters 252 dc[0] /= 1000.0; 253 dc[1] /= 1000.0; 254 dc[2] /= 1000.0; 255 256 return dc; 257 } 258 catch(Exception& e) { GPSTK_RETHROW(e); } 259 catch(exception& e) { 260 Exception E("std except: " + string(e.what())); 261 GPSTK_THROW(E); 262 } 263 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 264 265 } // end Triple AtmLoadTides::computeDisplacement 266 267 } // end namespace gpstk 268 //------------------------------------------------------------------------------------ 269 //------------------------------------------------------------------------------------ 270