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 /// @file OceanLoadTides.cpp 40 41 //------------------------------------------------------------------------------------ 42 // system includes 43 #include <iostream> 44 #include <fstream> 45 #include <algorithm> 46 // GPSTk 47 #include "StringUtils.hpp" 48 #include "MiscMath.hpp" 49 #include "GNSSconstants.hpp" 50 // geomatics 51 #include "OceanLoadTides.hpp" 52 #include "RobustStats.hpp" // for QSort 53 #include "CubicSpline.hpp" 54 //#include "logstream.hpp" // TEMP 55 56 using namespace std; 57 58 namespace gpstk { 59 60 using namespace StringUtils; 61 62 // Number of standard (Schwiderski) tides read from BLQ file 63 const int OceanLoadTides::NSTD=11; 64 // Number of derived tides computed by deriveTides() 65 const int OceanLoadTides::NDER=342; 66 67 //--------------------------------------------------------------------------------- 68 // Open and read the given file, containing ocean loading coefficients, and 69 // initialize this object for the sites names in the input list that match a 70 // name in the file (case sensitive). Return the number of successfully 71 // initialized site names, and remove those sites from the input list. 72 // Ocean loading files can be obtained from the web. For example all the ITRF 73 // sites are found at ftp://maia.usno.navy.mil/conventions/chapter7/olls25.blq 74 // Also, at http://www.oso.chalmers.se/~loading one may submit site label and 75 // position for one or more sites, and the resulting ocean loading file will be 76 // emailed. 77 // @param sites vector<string> On input contains site labels found in the 78 // file, on output contains only sites that were NOT found. 79 // If sites is empty, all sites are loaded. 80 // @param filename string Input ocean loading file name. 81 // @return the number of sites successfully initialized. 82 // @throw if the file could not be opened. initializeSites(vector<string> & sites,string filename)83 int OceanLoadTides::initializeSites(vector<string>& sites, string filename) 84 { 85 try { 86 bool allsites = false; 87 if(sites.size() == 0) allsites = true; // return 0; 88 int i,n; 89 90 ifstream infile(filename.c_str()); 91 if(!infile || !infile.is_open()) { 92 Exception e("File " + filename + " could not be opened."); 93 GPSTK_THROW(e); 94 } 95 96 n = 0; // number of successes 97 bool looking=true; // true if looking for a site name 98 double lat,lon; 99 vector<double> coeff; 100 string site; 101 while(1) { // read the file 102 int count; 103 string line,word; 104 105 // get the next line 106 getline(infile,line); 107 stripTrailing(line,'\r'); 108 109 // process line 110 if(!line.empty()) { 111 word = firstWord(line); 112 //LOG(VERBOSE) << "Word is " << word << " and line is " << line; 113 114 if(word == "$$") { // NB ignore header - assume column order, etc. 115 // pick out the lat/lon 116 if(!looking) { 117 while(!line.empty()) { 118 word = stripFirstWord(line); 119 if(word == string("lon/lat:")) { 120 lon = asDouble(stripFirstWord(line)); 121 lat = asDouble(stripFirstWord(line)); 122 break; 123 } 124 } 125 } 126 } 127 // TD should test be line length <= 21 ? ... what if site name = number 128 //else if(looking && !isDecimalString(word)) { 129 else if(looking && line.length() <= 21) { 130 // site name 131 site = line; 132 stripTrailing(site,string("\n")); 133 stripTrailing(site,string("\r")); 134 stripTrailing(site); 135 stripLeading(site); 136 //LOG(VERBOSE) << "Found site " << site; 137 if(allsites) { 138 //LOG(VERBOSE) << "Push back " << site; 139 looking = false; 140 sites.push_back(site); 141 } 142 else for(i=0; i<sites.size(); i++) { 143 //LOG(VERBOSE) << "Compare " << sites[i]; 144 if(site == sites[i]) { 145 looking = false; 146 break; 147 } 148 } 149 if(!looking) { // found a site 150 count = 0; 151 coeff.clear(); 152 lat = lon = 0.0; 153 } 154 } 155 else if(!looking) { // not comment and not looking - must be data 156 if(numWords(line) != 11) { 157 Exception e("File " + filename + " is corrupted for site " + site 158 + " - offending line follows\n" + line); 159 GPSTK_THROW(e); 160 } 161 //LOG(VERBOSE) << "Push back line " << line; 162 for(i=0; i<11; i++) 163 coeff.push_back( 164 asDouble(stripFirstWord(line))); 165 count++; 166 if(count == 6) { // success 167 ostringstream oss; 168 oss << fixed; 169 for(i=0; i<coeff.size(); i++) { 170 if(i < 33) oss << " " << setprecision(5) << setw(7) << coeff[i]; 171 else oss << " " << setprecision(1) << setw(7) << coeff[i]; 172 if((i+1)%11 == 0) oss << "\n"; 173 } 174 //LOG(VERBOSE) << " Found site " << site << " with coefficients:"; 175 //LOG(VERBOSE) << oss.str(); 176 177 // update coeff map 178 coefficientMap[site] = coeff; 179 n++; 180 // update position map 181 coeff.clear(); 182 coeff.push_back(lat); 183 coeff.push_back(lon); 184 positionMap[site] = coeff; 185 186 // erase a vector element 187 if(!allsites) { 188 vector<string>::iterator pos; 189 pos = find(sites.begin(),sites.end(),site); 190 if(pos != sites.end()) sites.erase(pos); 191 } 192 looking = true; 193 } 194 } 195 196 } // end if line not empty 197 198 if(infile.eof() || !infile.good()) break; 199 200 } // end loop over lines in the file 201 202 return n; 203 } 204 catch(Exception& e) { GPSTK_RETHROW(e); } 205 catch(exception& e) { Exception E("std except: "+string(e.what())); GPSTK_THROW(E); } 206 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 207 } 208 209 //--------------------------------------------------------------------------------- 210 // Compute the site displacement vector at the given time for the given site. 211 // The site must have been successfully initialized; if not an exception is 212 // thrown. 213 // @param site string Input name of the site; must be the same as previously 214 // successfully passed to initializeSites(). 215 // @param t EphTime Input time of interest. 216 // @return Triple containing the North, East and Up components of the site 217 // displacement in meters. 218 // @throw if the site has not been initialized. computeDisplacement11(string site,EphTime time)219 Triple OceanLoadTides::computeDisplacement11(string site, EphTime time) 220 { 221 try { 222 if(!isValid(site)) { 223 Exception e("Site " + site + " has not been initialized."); 224 GPSTK_THROW(e); 225 } 226 227 // get the coefficients for this site 228 vector<double> coeff = coefficientMap[site]; 229 230 // get the astronomical arguments in radians 231 double angles[11]; 232 //inline this SchwiderskiArg(int(t.year())-1900, t.DOY(), t.secOfDay(), angles); 233 { 234 double fday(time.secOfDay()); 235 long jday(static_cast<long>(time.lMJD() + MJD_JDAY + fday/SEC_PER_DAY)); 236 int iyear,imm,iday; 237 convertJDtoCalendar(jday,iyear,imm,iday); 238 iyear -= 1900; 239 240 // ordering is: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa 241 // which are : { semi-diurnal }{ diurnal }{long-period} 242 static const double speed[11] = { 243 1.40519E-4, 1.45444E-4, 1.37880E-4, 1.45842E-4, 244 0.72921E-4, 0.67598E-4, 0.72523E-4, 0.64959E-4, 245 0.053234E-4, 0.026392E-4, 0.003982E-4 }; 246 static const double angfac[44] = 247 { 248 // sun 249 2.0, 0.0, 2.0, 2.0, // 4 : M2, S2, N2, K2 250 1.0, 1.0, -1.0, 1.0, // 8 : K1, O1, P1, Q1 251 0.0, 0.0, 2.0, // 11 : Mf, Mm, Ssa 252 // moon 253 -2.0, 0.0, -3.0, 0.0, // 15 : M2, S2, N2, K2 254 0.0, -2.0, 0.0, -3.0, // 19 : K1, O1, P1, Q1 255 2.0, 1.0, 0.0, // 22 : Mf, Mm, Ssa 256 // lunar perigee 257 0.0, 0.0, 1.0, 0.0, // 26 : M2, S2, N2, K2 258 0.0, 0.0, 0.0, 1.0, // 30 : K1, O1, P1, Q1 259 0.0, -1.0, 0.0, // 33 : Mf, Mm, Ssa 260 // two pi 261 0.0, 0.0, 0.0, 0.0, // 37 : M2, S2, N2, K2 262 0.25,-0.25,-0.25,-0.25, // 41 : K1, O1, P1, Q1 263 0.0, 0.0, 0.0 // 44 : Mf, Mm, Ssa 264 }; 265 266 int icapd = iday + 365*(iyear-75)+((iyear-73)/4); 267 268 //double capt = (27392.500528+1.000000035*double(icapd))/36525.0; 269 double capt = 0.74996579132101300 + 2.73785088295687885e-5 * double(icapd); 270 271 // mean longitude of sun at beginning of day 272 double H0 = 279.69668+(36000.768930485+0.000303*capt)*capt; 273 274 // mean longitude of moon at beginning of day 275 double S0 = ((0.0000019*capt-0.001133)*capt+481267.88314137)*capt+270.434358; 276 277 // mean longitude of lunar perigee at beginning of day 278 double P0 = ((-0.000012*capt-0.010325)*capt+4069.0340329577)*capt+334.329653; 279 280 // convert to radians 281 //static const double dtr = 0.0174532925199; 282 H0 *= DEG_TO_RAD; 283 S0 *= DEG_TO_RAD; 284 P0 *= DEG_TO_RAD; 285 286 //LOG(INFO) << "Schwiderski " << iday << " " << fixed << setprecision(5) 287 //<< setw(11) << fday << " " << icapd << " " << capt 288 //<< " " << H0 << " " << S0 << " " << P0; 289 290 static const double twopi = 6.28318530718; 291 for(int k=0; k<11; k++) { 292 angles[k] = speed[k]*fday + angfac[k]*H0 293 + angfac[11+k]*S0 294 + angfac[22+k]*P0 295 + angfac[33+k]*twopi; 296 angles[k] = ::fmod(angles[k],twopi); 297 if(angles[k] < 0.0) angles[k] += twopi; 298 } 299 } // end SchwiderskiArg() 300 301 // compute the radial, west and south components 302 // coefficients are stored by rows: radial, west, south; first amp, then phase 303 // column order same as in SchwiderskiArg() [ as in the file ] 304 Triple dc; 305 for(int i=0; i<3; i++) { // components 306 dc[i] = 0.0; 307 for(int j=0; j<11; j++) // tidal modes 308 dc[i] += coeff[i*11+j]*::cos(angles[j]-coeff[33+i*11+j]*DEG_TO_RAD); 309 } 310 311 // convert radial,west,south to north,east,up 312 double temp=dc[0]; 313 dc[0] = -dc[2]; // N = -S 314 dc[1] = -dc[1]; // E = -W 315 dc[2] = temp; // U = rad 316 317 return dc; 318 } 319 catch(Exception& e) { GPSTK_RETHROW(e); } 320 catch(exception& e) { 321 Exception E("std except: " + string(e.what())); 322 GPSTK_THROW(E); 323 } 324 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 325 } 326 327 //--------------------------------------------------------------------------------- 328 // Compute the site displacement vector at the given time for the given site. 329 // The site must have been successfully initialized; if not an exception is 330 // thrown. 331 // @param site string Input name of the site; must be the same as previously 332 // successfully passed to initializeSites(). 333 // @param t EphTime Input time of interest. 334 // @return Triple containing the North, East and Up components of the site 335 // displacement in meters. 336 // @throw if the site has not been initialized, if the time system is unknown, 337 // if there is corruption in the static arrays, or . computeDisplacement(string site,EphTime time)338 Triple OceanLoadTides::computeDisplacement(string site, EphTime time) 339 { 340 try { 341 ostringstream oss; // TEMP 342 int i; 343 344 if(!isValid(site)) { 345 Exception e("Site " + site + " has not been initialized."); 346 GPSTK_THROW(e); 347 } 348 349 // get the coefficients for this site 350 vector<double> coeff = coefficientMap[site]; 351 352 // Cartwright-Tayler numbers of Scherneck tides 353 // ordering is: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa 354 355 // standard 11 Scherneck tides: 356 static const NVector SchInd[] = { 357 { 2, 0, 0, 0, 0, 0 }, // M2 358 { 2, 2,-2, 0, 0, 0 }, // S2 359 { 2,-1, 0, 1, 0, 0 }, // N2 360 { 2, 2, 0, 0, 0, 0 }, // K2 361 { 1, 1, 0, 0, 0, 0 }, // K1 362 { 1,-1, 0, 0, 0, 0 }, // O1 363 { 1, 1,-2, 0, 0, 0 }, // P1 364 { 1,-2, 0, 1, 0, 0 }, // Q1 365 { 0, 2, 0, 0, 0, 0 }, // Mf 366 { 0, 1, 0,-1, 0, 0 }, // Mm 367 { 0, 0, 2, 0, 0, 0 }, // Ssa 368 }; 369 370 // NB there must be 11 std tides in SchInd[] 371 if((int)(sizeof(SchInd) / sizeof(NVector)) != NSTD) { 372 Exception e("Static SchInd array is corrupted"); 373 GPSTK_THROW(e); 374 } 375 376 // compute time argument 377 EphTime ttag(time); 378 ttag.convertSystemTo(TimeSystem::UTC); 379 double dayfr(ttag.secOfDay()/86400.0); 380 ttag.convertSystemTo(TimeSystem::TT); 381 // T = EarthOrientation::CoordTransTime() 382 double T((ttag.dMJD() - 51544.5)/36525.0); 383 384 // get the Delauney arguments and frequencies at t 385 double Del[5], freqDel[5]; // degrees and cycles/day 386 Del[0] = 134.9634025100 + // EarthOrientation::L() 387 T*(477198.8675605000 + 388 T*( 0.0088553333 + 389 T*( 0.0000143431 + 390 T*( -0.0000000680)))); 391 Del[1] = 357.5291091806 + // EarthOrientation::Lp() 392 T*( 35999.0502911389 + 393 T*( -0.0001536667 + 394 T*( 0.0000000378 + 395 T*( -0.0000000032)))); 396 Del[2] = 93.2720906200 + // EarthOrientation::F() 397 T*(483202.0174577222 + 398 T*( -0.0035420000 + 399 T*( -0.0000002881 + 400 T*( 0.0000000012)))); 401 Del[3] = 297.8501954694 + // EarthOrientation::D() 402 T*(445267.1114469445 + 403 T*( -0.0017696111 + 404 T*( 0.0000018314 + 405 T*( -0.0000000088)))); 406 Del[4] = 125.0445550100 + // EarthOrientation::Omega2003() 407 T*( -1934.1362619722 + 408 T*( 0.0020756111 + 409 T*( 0.0000021394 + 410 T*( -0.0000000165)))); 411 for(i=0; i<5; i++) Del[i] = ::fmod(Del[i],360.0); 412 freqDel[0] = 0.0362916471 + 0.0000000013*T; 413 freqDel[1] = 0.0027377786; 414 freqDel[2] = 0.0367481951 - 0.0000000005*T; 415 freqDel[3] = 0.0338631920 - 0.0000000003*T; 416 freqDel[4] = -0.0001470938 + 0.0000000003*T; 417 418 // convert to Doodson (Darwin) variables 419 double Dood[6], freqDood[6]; 420 Dood[0] = 360.0*dayfr - Del[3]; 421 Dood[1] = Del[2] + Del[4]; 422 Dood[2] = Dood[1] - Del[3]; 423 Dood[3] = Dood[1] - Del[0]; 424 Dood[4] = -Del[4]; 425 Dood[5] = Dood[2] - Del[1]; 426 for(i=0; i<6; i++) Dood[i] = ::fmod(Dood[i],360.0); 427 428 freqDood[0] = 1.0 - freqDel[3]; 429 freqDood[1] = freqDel[2] + freqDel[4]; 430 freqDood[2] = freqDood[1] - freqDel[3]; 431 freqDood[3] = freqDood[1] - freqDel[0]; 432 freqDood[4] = -freqDel[4]; 433 freqDood[5] = freqDood[2] - freqDel[1]; 434 435 // find amplitudes and phases for vertical, west and south components, 436 // for all 342 derived tides, from standard tides 437 double amp[NSTD],phs[NSTD]; 438 double ampS[NDER],ampW[NDER],ampU[NDER]; // south,west,up component amp.s 439 double phsS[NDER],phsW[NDER],phsU[NDER]; // south,west,up component phs.s 440 double freq[NDER]; // frequencies (same for S,W,U) 441 442 // vertical 443 int nder; // number returned, may be < NDER 444 for(i=0; i<NSTD; i++) { 445 amp[i] = coeff[i]; 446 phs[i] = -coeff[33+i]; 447 } 448 //oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 1 "; 449 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i]; 450 //LOG(INFO) << oss.str(); 451 //oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 1 "; 452 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i]; 453 //LOG(INFO) << oss.str(); 454 //LOG(INFO) << "TEST T,DAYFR,DELTA" << fixed << setprecision(15) 455 // << setw(25) << T << setw(25) << dayfr << setw(25) << ttag.secOfDay(); 456 //LOG(INFO) << "TEST Delauneys " << fixed << setprecision(15) 457 // << setw(22) << Del[0] << setw(22) << Del[1] << setw(22) << Del[2] 458 // << setw(22) << Del[3] << setw(22) << Del[4]; 459 //LOG(INFO) << "TEST Del freqs " << fixed << setprecision(15) 460 // << setw(22) << freqDel[0] << setw(22) << freqDel[1] << setw(22) 461 // << freqDel[2] << setw(22) << freqDel[3] << setw(22) << freqDel[4]; 462 //LOG(INFO) << "TEST Doods " << fixed << setprecision(15) 463 // << setw(22) << Dood[0] << setw(22) << Dood[1] << setw(22) 464 // << Dood[2] << setw(22) << Dood[3] << setw(22) << Dood[4] 465 // << setw(22) << Dood[5]; 466 //LOG(INFO) << "TEST Dood freqs" << fixed << setprecision(15) 467 // << setw(22) << freqDood[0] << setw(22) << freqDood[1] << setw(22) 468 // << freqDood[2] << setw(22) << freqDood[3] << setw(22) << freqDood[4] 469 // << setw(22) << freqDood[5]; 470 nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampU, phsU, freq, NSTD); 471 //LOG(INFO) << "Vertical returned " << nder << " derived tides"; 472 473 // west 474 for(i=0; i<NSTD; i++) { 475 amp[i] = coeff[11+i]; 476 phs[i] = -coeff[44+i]; 477 } 478 //oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 2 "; 479 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i]; 480 //LOG(INFO) << oss.str(); 481 //oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 2 "; 482 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i]; 483 //LOG(INFO) << oss.str(); 484 nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampW, phsW, freq, NSTD); 485 //LOG(INFO) << "West returned " << nder << " derived tides"; 486 487 // south 488 for(i=0; i<NSTD; i++) { 489 amp[i] = coeff[22+i]; 490 phs[i] = -coeff[55+i]; 491 } 492 //oss.str(""); oss << fixed << setprecision(5) << "TEST Amp 3 "; 493 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << amp[i]; 494 //LOG(INFO) << oss.str(); 495 //oss.str(""); oss << fixed << setprecision(1) << "TEST Phs 3 "; 496 //for(i=0; i<NSTD; i++) oss << " " << setw(8) << phs[i]; 497 //LOG(INFO) << oss.str(); 498 nder = deriveTides(SchInd, amp, phs, Dood, freqDood, ampS, phsS, freq, NSTD); 499 //LOG(INFO) << "TEST First 40 South amp, phase"; 500 //for(i=0; i<40; i++) LOG(INFO) << "TEST " << setw(2) << i+1 << fixed 501 // << setprecision(15) << setw(22) << ampS[i] << setw(22) << phsS[i]; 502 503 // sum up 504 Triple dc(0.0,0.0,0.0); // U S W 505 for(i=0; i<nder; i++) { 506 dc[0] += ampU[i] * ::cos(phsU[i]*DEG_TO_RAD); 507 //LOG(INFO)<<"TEST LOOP U " << setw(3) << i+1 << fixed << setprecision(15) 508 // << setw(22) << ampU[i]*::cos(phsU[i]*DEG_TO_RAD) << setw(22) << dc[0]; 509 } 510 //LOG(INFO) << "TEST RECURS result U " << fixed << setprecision(15) 511 // << setw(22) << dc[0]; 512 513 for(i=0; i<nder; i++) { 514 dc[1] += ampS[i] * ::cos(phsS[i]*DEG_TO_RAD); 515 //LOG(INFO)<<"TEST LOOP S " << setw(3) << i+1 << fixed << setprecision(15) 516 // << setw(22) << ampS[i]*::cos(phsS[i]*DEG_TO_RAD) << setw(22) << dc[1]; 517 } 518 //LOG(INFO) << "TEST RECURS result S " << fixed << setprecision(15) 519 // << setw(22) << dc[1]; 520 521 for(i=0; i<nder; i++) { 522 dc[2] += ampW[i] * ::cos(phsW[i]*DEG_TO_RAD); 523 //LOG(INFO)<<"TEST LOOP W " << setw(3) << i+1 << fixed << setprecision(15) 524 // << setw(22) << ampW[i]*::cos(phsW[i]*DEG_TO_RAD) << setw(22) << dc[2]; 525 } 526 //LOG(INFO) << "TEST RECURS result W " << fixed << setprecision(15) 527 // << setw(22) << dc[2]; 528 529 //LOG(INFO) << "TEST " << fixed << setprecision(6) 530 // << " " << dc[0] << " " << dc[1] << " " << dc[2]; 531 532 // convert vertical,south,west to north,east,up 533 double temp=dc[0]; 534 dc[0] = -dc[1]; // N = -S 535 dc[1] = -dc[2]; // E = -W 536 dc[2] = temp; // U = U 537 538 return dc; 539 } 540 catch(Exception& e) { GPSTK_RETHROW(e); } 541 catch(exception& e) { 542 Exception E("std except: " + string(e.what())); 543 GPSTK_THROW(E); 544 } 545 catch(...) { Exception e("Unknown exception"); GPSTK_THROW(e); } 546 547 } // end Triple OceanLoadTides::computeDisplacement 548 549 //--------------------------------------------------------------------------------- deriveTides(const NVector SchInd[],const double amp[],const double phs[],const double Dood[],const double freqDood[],double ampDer[],double phsDer[],double freqDer[],const int Nin)550 int OceanLoadTides::deriveTides(const NVector SchInd[], 551 const double amp[], const double phs[], 552 const double Dood[], const double freqDood[], 553 double ampDer[], double phsDer[], double freqDer[], 554 const int Nin) 555 { 556 // indexes for std tides: M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, Ssa 557 static const int stdindex[] = { 558 0, 1, 2, 3,109, 110, 111, 112, 263, 264, 265 }; 559 560 static const double DerAmp[] = { 561 .632208, .294107, .121046, .079915, .023818,-.023589, .022994, 562 .019333,-.017871, .017192, .016018, .004671,-.004662,-.004519, 563 .004470, .004467, .002589,-.002455,-.002172, .001972, .001947, 564 .001914,-.001898, .001802, .001304, .001170, .001130, .001061, 565 -.001022,-.001017, .001014, .000901,-.000857, .000855, .000855, 566 .000772, .000741, .000741,-.000721, .000698, .000658, .000654, 567 -.000653, .000633, .000626,-.000598, .000590, .000544, .000479, 568 -.000464, .000413,-.000390, .000373, .000366, .000366,-.000360, 569 -.000355, .000354, .000329, .000328, .000319, .000302, .000279, 570 -.000274,-.000272, .000248,-.000225, .000224,-.000223,-.000216, 571 .000211, .000209, .000194, .000185,-.000174,-.000171, .000159, 572 .000131, .000127, .000120, .000118, .000117, .000108, .000107, 573 .000105,-.000102, .000102, .000099,-.000096, .000095,-.000089, 574 -.000085,-.000084,-.000081,-.000077,-.000072,-.000067, .000066, 575 .000064, .000063, .000063, .000063, .000062, .000062,-.000060, 576 .000056, .000053, .000051, .000050, .368645,-.262232,-.121995, 577 -.050208, .050031,-.049470, .020620, .020613, .011279,-.009530, 578 -.009469,-.008012, .007414,-.007300, .007227,-.007131,-.006644, 579 .005249, .004137, .004087, .003944, .003943, .003420, .003418, 580 .002885, .002884, .002160,-.001936, .001934,-.001798, .001690, 581 .001689, .001516, .001514,-.001511, .001383, .001372, .001371, 582 -.001253,-.001075, .001020, .000901, .000865,-.000794, .000788, 583 .000782,-.000747,-.000745, .000670,-.000603,-.000597, .000542, 584 .000542,-.000541,-.000469,-.000440, .000438, .000422, .000410, 585 -.000374,-.000365, .000345, .000335,-.000321,-.000319, .000307, 586 .000291, .000290,-.000289, .000286, .000275, .000271, .000263, 587 -.000245, .000225, .000225, .000221,-.000202,-.000200,-.000199, 588 .000192, .000183, .000183, .000183,-.000170, .000169, .000168, 589 .000162, .000149,-.000147,-.000141, .000138, .000136, .000136, 590 .000127, .000127,-.000126,-.000121,-.000121, .000117,-.000116, 591 -.000114,-.000114,-.000114, .000114, .000113, .000109, .000108, 592 .000106,-.000106,-.000106, .000105, .000104,-.000103,-.000100, 593 -.000100,-.000100, .000099,-.000098, .000093, .000093, .000090, 594 -.000088, .000083,-.000083,-.000082,-.000081,-.000079,-.000077, 595 -.000075,-.000075,-.000075, .000071, .000071,-.000071, .000068, 596 .000068, .000065, .000065, .000064, .000064, .000064,-.000064, 597 -.000060, .000056, .000056, .000053, .000053, .000053,-.000053, 598 .000053, .000053, .000052, .000050,-.066607,-.035184,-.030988, 599 .027929,-.027616,-.012753,-.006728,-.005837,-.005286,-.004921, 600 -.002884,-.002583,-.002422, .002310, .002283,-.002037, .001883, 601 -.001811,-.001687,-.001004,-.000925,-.000844, .000766, .000766, 602 -.000700,-.000495,-.000492, .000491, .000483, .000437,-.000416, 603 -.000384, .000374,-.000312,-.000288,-.000273, .000259, .000245, 604 -.000232, .000229,-.000216, .000206,-.000204,-.000202, .000200, 605 .000195,-.000190, .000187, .000180,-.000179, .000170, .000153, 606 -.000137,-.000119,-.000119,-.000112,-.000110,-.000110, .000107, 607 -.000095,-.000095,-.000091,-.000090,-.000081,-.000079,-.000079, 608 .000077,-.000073, .000069,-.000067,-.000066, .000065, .000064, 609 -.000062, .000060, .000059,-.000056, .000055,-.000051 }; 610 611 static const NVector DerInd[] = { 612 { 2, 0, 0, 0, 0, 0 }, { 2, 2,-2, 0, 0, 0 }, { 2,-1, 0, 1, 0, 0 },//M2,S2,N2 613 { 2, 2, 0, 0, 0, 0 }, { 2, 2, 0, 0, 1, 0 }, { 2, 0, 0, 0,-1, 0 },//K2,x,x 614 { 2,-1, 2,-1, 0, 0 }, { 2,-2, 2, 0, 0, 0 }, { 2, 1, 0,-1, 0, 0 }, 615 { 2, 2,-3, 0, 0, 1 }, { 2,-2, 0, 2, 0, 0 }, { 2,-3, 2, 1, 0, 0 }, 616 { 2, 1,-2, 1, 0, 0 }, { 2,-1, 0, 1,-1, 0 }, { 2, 3, 0,-1, 0, 0 }, 617 { 2, 1, 0, 1, 0, 0 }, { 2, 2, 0, 0, 2, 0 }, { 2, 2,-1, 0, 0,-1 }, 618 { 2, 0,-1, 0, 0, 1 }, { 2, 1, 0, 1, 1, 0 }, { 2, 3, 0,-1, 1, 0 }, 619 { 2, 0, 1, 0, 0,-1 }, { 2, 0,-2, 2, 0, 0 }, { 2,-3, 0, 3, 0, 0 }, 620 { 2,-2, 3, 0, 0,-1 }, { 2, 4, 0, 0, 0, 0 }, { 2,-1, 1, 1, 0,-1 }, 621 { 2,-1, 3,-1, 0,-1 }, { 2, 2, 0, 0,-1, 0 }, { 2,-1,-1, 1, 0, 1 }, 622 { 2, 4, 0, 0, 1, 0 }, { 2,-3, 4,-1, 0, 0 }, { 2,-1, 2,-1,-1, 0 }, 623 { 2, 3,-2, 1, 0, 0 }, { 2, 1, 2,-1, 0, 0 }, { 2,-4, 2, 2, 0, 0 }, 624 { 2, 4,-2, 0, 0, 0 }, { 2, 0, 2, 0, 0, 0 }, { 2,-2, 2, 0,-1, 0 }, 625 { 2, 2,-4, 0, 0, 2 }, { 2, 2,-2, 0,-1, 0 }, { 2, 1, 0,-1,-1, 0 }, 626 { 2,-1, 1, 0, 0, 0 }, { 2, 2,-1, 0, 0, 1 }, { 2, 2, 1, 0, 0,-1 }, 627 { 2,-2, 0, 2,-1, 0 }, { 2,-2, 4,-2, 0, 0 }, { 2, 2, 2, 0, 0, 0 }, 628 { 2,-4, 4, 0, 0, 0 }, { 2,-1, 0,-1,-2, 0 }, { 2, 1, 2,-1, 1, 0 }, 629 { 2,-1,-2, 3, 0, 0 }, { 2, 3,-2, 1, 1, 0 }, { 2, 4, 0,-2, 0, 0 }, 630 { 2, 0, 0, 2, 0, 0 }, { 2, 0, 2,-2, 0, 0 }, { 2, 0, 2, 0, 1, 0 }, 631 { 2,-3, 3, 1, 0,-1 }, { 2, 0, 0, 0,-2, 0 }, { 2, 4, 0, 0, 2, 0 }, 632 { 2, 4,-2, 0, 1, 0 }, { 2, 0, 0, 0, 0, 2 }, { 2, 1, 0, 1, 2, 0 }, 633 { 2, 0,-2, 0,-2, 0 }, { 2,-2, 1, 0, 0, 1 }, { 2,-2, 1, 2, 0,-1 }, 634 { 2,-1, 1,-1, 0, 1 }, { 2, 5, 0,-1, 0, 0 }, { 2, 1,-3, 1, 0, 1 }, 635 { 2,-2,-1, 2, 0, 1 }, { 2, 3, 0,-1, 2, 0 }, { 2, 1,-2, 1,-1, 0 }, 636 { 2, 5, 0,-1, 1, 0 }, { 2,-4, 0, 4, 0, 0 }, { 2,-3, 2, 1,-1, 0 }, 637 { 2,-2, 1, 1, 0, 0 }, { 2, 4, 0,-2, 1, 0 }, { 2, 0, 0, 2, 1, 0 }, 638 { 2,-5, 4, 1, 0, 0 }, { 2, 0, 2, 0, 2, 0 }, { 2,-1, 2, 1, 0, 0 }, 639 { 2, 5,-2,-1, 0, 0 }, { 2, 1,-1, 0, 0, 0 }, { 2, 2,-2, 0, 0, 2 }, 640 { 2,-5, 2, 3, 0, 0 }, { 2,-1,-2, 1,-2, 0 }, { 2,-3, 5,-1, 0,-1 }, 641 { 2,-1, 0, 0, 0, 1 }, { 2,-2, 0, 0,-2, 0 }, { 2, 0,-1, 1, 0, 0 }, 642 { 2,-3, 1, 1, 0, 1 }, { 2, 3, 0,-1,-1, 0 }, { 2, 1, 0, 1,-1, 0 }, 643 { 2,-1, 2, 1, 1, 0 }, { 2, 0,-3, 2, 0, 1 }, { 2, 1,-1,-1, 0, 1 }, 644 { 2,-3, 0, 3,-1, 0 }, { 2, 0,-2, 2,-1, 0 }, { 2,-4, 3, 2, 0,-1 }, 645 { 2,-1, 0, 1,-2, 0 }, { 2, 5, 0,-1, 2, 0 }, { 2,-4, 5, 0, 0,-1 }, 646 { 2,-2, 4, 0, 0,-2 }, { 2,-1, 0, 1, 0, 2 }, { 2,-2,-2, 4, 0, 0 }, 647 { 2, 3,-2,-1,-1, 0 }, { 2,-2, 5,-2, 0,-1 }, { 2, 0,-1, 0,-1, 1 }, 648 { 2, 5,-2,-1, 1, 0 }, { 1, 1, 0, 0, 0, 0 }, { 1,-1, 0, 0, 0, 0 },//x,K1,O1 649 { 1, 1,-2, 0, 0, 0 }, { 1,-2, 0, 1, 0, 0 }, { 1, 1, 0, 0, 1, 0 },//P1,Q1,x 650 { 1,-1, 0, 0,-1, 0 }, { 1, 2, 0,-1, 0, 0 }, { 1, 0, 0, 1, 0, 0 }, 651 { 1, 3, 0, 0, 0, 0 }, { 1,-2, 2,-1, 0, 0 }, { 1,-2, 0, 1,-1, 0 }, 652 { 1,-3, 2, 0, 0, 0 }, { 1, 0, 0,-1, 0, 0 }, { 1, 1, 0, 0,-1, 0 }, 653 { 1, 3, 0, 0, 1, 0 }, { 1, 1,-3, 0, 0, 1 }, { 1,-3, 0, 2, 0, 0 }, 654 { 1, 1, 2, 0, 0, 0 }, { 1, 0, 0, 1, 1, 0 }, { 1, 2, 0,-1, 1, 0 }, 655 { 1, 0, 2,-1, 0, 0 }, { 1, 2,-2, 1, 0, 0 }, { 1, 3,-2, 0, 0, 0 }, 656 { 1,-1, 2, 0, 0, 0 }, { 1, 1, 1, 0, 0,-1 }, { 1, 1,-1, 0, 0, 1 }, 657 { 1, 4, 0,-1, 0, 0 }, { 1,-4, 2, 1, 0, 0 }, { 1, 0,-2, 1, 0, 0 }, 658 { 1,-2, 2,-1,-1, 0 }, { 1, 3, 0,-2, 0, 0 }, { 1,-1, 0, 2, 0, 0 }, 659 { 1,-1, 0, 0,-2, 0 }, { 1, 3, 0, 0, 2, 0 }, { 1,-3, 2, 0,-1, 0 }, 660 { 1, 4, 0,-1, 1, 0 }, { 1, 0, 0,-1,-1, 0 }, { 1, 1,-2, 0,-1, 0 }, 661 { 1,-3, 0, 2,-1, 0 }, { 1, 1, 0, 0, 2, 0 }, { 1, 1,-1, 0, 0,-1 }, 662 { 1,-1,-1, 0, 0, 1 }, { 1, 0, 2,-1, 1, 0 }, { 1,-1, 1, 0, 0,-1 }, 663 { 1,-1,-2, 2, 0, 0 }, { 1, 2,-2, 1, 1, 0 }, { 1,-4, 0, 3, 0, 0 }, 664 { 1,-1, 2, 0, 1, 0 }, { 1, 3,-2, 0, 1, 0 }, { 1, 2, 0,-1,-1, 0 }, 665 { 1, 0, 0, 1,-1, 0 }, { 1,-2, 2, 1, 0, 0 }, { 1, 4,-2,-1, 0, 0 }, 666 { 1,-3, 3, 0, 0,-1 }, { 1,-2, 1, 1, 0,-1 }, { 1,-2, 3,-1, 0,-1 }, 667 { 1, 0,-2, 1,-1, 0 }, { 1,-2,-1, 1, 0, 1 }, { 1, 4,-2, 1, 0, 0 }, 668 { 1,-4, 4,-1, 0, 0 }, { 1,-4, 2, 1,-1, 0 }, { 1, 5,-2, 0, 0, 0 }, 669 { 1, 3, 0,-2, 1, 0 }, { 1,-5, 2, 2, 0, 0 }, { 1, 2, 0, 1, 0, 0 }, 670 { 1, 1, 3, 0, 0,-1 }, { 1,-2, 0, 1,-2, 0 }, { 1, 4, 0,-1, 2, 0 }, 671 { 1, 1,-4, 0, 0, 2 }, { 1, 5, 0,-2, 0, 0 }, { 1,-1, 0, 2, 1, 0 }, 672 { 1,-2, 1, 0, 0, 0 }, { 1, 4,-2, 1, 1, 0 }, { 1,-3, 4,-2, 0, 0 }, 673 { 1,-1, 3, 0, 0,-1 }, { 1, 3,-3, 0, 0, 1 }, { 1, 5,-2, 0, 1, 0 }, 674 { 1, 1, 2, 0, 1, 0 }, { 1, 2, 0, 1, 1, 0 }, { 1,-5, 4, 0, 0, 0 }, 675 { 1,-2, 0,-1,-2, 0 }, { 1, 5, 0,-2, 1, 0 }, { 1, 1, 2,-2, 0, 0 }, 676 { 1, 1,-2, 2, 0, 0 }, { 1,-2, 2, 1, 1, 0 }, { 1, 0, 3,-1, 0,-1 }, 677 { 1, 2,-3, 1, 0, 1 }, { 1,-2,-2, 3, 0, 0 }, { 1,-1, 2,-2, 0, 0 }, 678 { 1,-4, 3, 1, 0,-1 }, { 1,-4, 0, 3,-1, 0 }, { 1,-1,-2, 2,-1, 0 }, 679 { 1,-2, 0, 3, 0, 0 }, { 1, 4, 0,-3, 0, 0 }, { 1, 0, 1, 1, 0,-1 }, 680 { 1, 2,-1,-1, 0, 1 }, { 1, 2,-2, 1,-1, 0 }, { 1, 0, 0,-1,-2, 0 }, 681 { 1, 2, 0, 1, 2, 0 }, { 1, 2,-2,-1,-1, 0 }, { 1, 0, 0, 1, 2, 0 }, 682 { 1, 0, 1, 0, 0, 0 }, { 1, 2,-1, 0, 0, 0 }, { 1, 0, 2,-1,-1, 0 }, 683 { 1,-1,-2, 0,-2, 0 }, { 1,-3, 1, 0, 0, 1 }, { 1, 3,-2, 0,-1, 0 }, 684 { 1,-1,-1, 0,-1, 1 }, { 1, 4,-2,-1, 1, 0 }, { 1, 2, 1,-1, 0,-1 }, 685 { 1, 0,-1, 1, 0, 1 }, { 1,-2, 4,-1, 0, 0 }, { 1, 4,-4, 1, 0, 0 }, 686 { 1,-3, 1, 2, 0,-1 }, { 1,-3, 3, 0,-1,-1 }, { 1, 1, 2, 0, 2, 0 }, 687 { 1, 1,-2, 0,-2, 0 }, { 1, 3, 0, 0, 3, 0 }, { 1,-1, 2, 0,-1, 0 }, 688 { 1,-2, 1,-1, 0, 1 }, { 1, 0,-3, 1, 0, 1 }, { 1,-3,-1, 2, 0, 1 }, 689 { 1, 2, 0,-1, 2, 0 }, { 1, 6,-2,-1, 0, 0 }, { 1, 2, 2,-1, 0, 0 }, 690 { 1,-1, 1, 0,-1,-1 }, { 1,-2, 3,-1,-1,-1 }, { 1,-1, 0, 0, 0, 2 }, 691 { 1,-5, 0, 4, 0, 0 }, { 1, 1, 0, 0, 0,-2 }, { 1,-2, 1, 1,-1,-1 }, 692 { 1, 1,-1, 0, 1, 1 }, { 1, 1, 2, 0, 0,-2 }, { 1,-3, 1, 1, 0, 0 }, 693 { 1,-4, 4,-1,-1, 0 }, { 1, 1, 0,-2,-1, 0 }, { 1,-2,-1, 1,-1, 1 }, 694 { 1,-3, 2, 2, 0, 0 }, { 1, 5,-2,-2, 0, 0 }, { 1, 3,-4, 2, 0, 0 }, 695 { 1, 1,-2, 0, 0, 2 }, { 1,-1, 4,-2, 0, 0 }, { 1, 2, 2,-1, 1, 0 }, 696 { 1,-5, 2, 2,-1, 0 }, { 1, 1,-3, 0,-1, 1 }, { 1, 1, 1, 0, 1,-1 }, 697 { 1, 6,-2,-1, 1, 0 }, { 1,-2, 2,-1,-2, 0 }, { 1, 4,-2, 1, 2, 0 }, 698 { 1,-6, 4, 1, 0, 0 }, { 1, 5,-4, 0, 0, 0 }, { 1,-3, 4, 0, 0, 0 }, 699 { 1, 1, 2,-2, 1, 0 }, { 1,-2, 1, 0,-1, 0 }, { 0, 2, 0, 0, 0, 0 },//x,x,Mf 700 { 0, 1, 0,-1, 0, 0 }, { 0, 0, 2, 0, 0, 0 }, { 0, 0, 0, 0, 1, 0 },//Mm,SSa 701 { 0, 2, 0, 0, 1, 0 }, { 0, 3, 0,-1, 0, 0 }, { 0, 1,-2, 1, 0, 0 }, 702 { 0, 2,-2, 0, 0, 0 }, { 0, 3, 0,-1, 1, 0 }, { 0, 0, 1, 0, 0,-1 }, 703 { 0, 2, 0,-2, 0, 0 }, { 0, 2, 0, 0, 2, 0 }, { 0, 3,-2, 1, 0, 0 }, 704 { 0, 1, 0,-1,-1, 0 }, { 0, 1, 0,-1, 1, 0 }, { 0, 4,-2, 0, 0, 0 }, 705 { 0, 1, 0, 1, 0, 0 }, { 0, 0, 3, 0, 0,-1 }, { 0, 4, 0,-2, 0, 0 }, 706 { 0, 3,-2, 1, 1, 0 }, { 0, 3,-2,-1, 0, 0 }, { 0, 4,-2, 0, 1, 0 }, 707 { 0, 0, 2, 0, 1, 0 }, { 0, 1, 0, 1, 1, 0 }, { 0, 4, 0,-2, 1, 0 }, 708 { 0, 3, 0,-1, 2, 0 }, { 0, 5,-2,-1, 0, 0 }, { 0, 1, 2,-1, 0, 0 }, 709 { 0, 1,-2, 1,-1, 0 }, { 0, 1,-2, 1, 1, 0 }, { 0, 2,-2, 0,-1, 0 }, 710 { 0, 2,-3, 0, 0, 1 }, { 0, 2,-2, 0, 1, 0 }, { 0, 0, 2,-2, 0, 0 }, 711 { 0, 1,-3, 1, 0, 1 }, { 0, 0, 0, 0, 2, 0 }, { 0, 0, 1, 0, 0, 1 }, 712 { 0, 1, 2,-1, 1, 0 }, { 0, 3, 0,-3, 0, 0 }, { 0, 2, 1, 0, 0,-1 }, 713 { 0, 1,-1,-1, 0, 1 }, { 0, 1, 0, 1, 2, 0 }, { 0, 5,-2,-1, 1, 0 }, 714 { 0, 2,-1, 0, 0, 1 }, { 0, 2, 2,-2, 0, 0 }, { 0, 1,-1, 0, 0, 0 }, 715 { 0, 5, 0,-3, 0, 0 }, { 0, 2, 0,-2, 1, 0 }, { 0, 1, 1,-1, 0,-1 }, 716 { 0, 3,-4, 1, 0, 0 }, { 0, 0, 2, 0, 2, 0 }, { 0, 2, 0,-2,-1, 0 }, 717 { 0, 4,-3, 0, 0, 1 }, { 0, 3,-1,-1, 0, 1 }, { 0, 0, 2, 0, 0,-2 }, 718 { 0, 3,-3, 1, 0, 1 }, { 0, 2,-4, 2, 0, 0 }, { 0, 4,-2,-2, 0, 0 }, 719 { 0, 3, 1,-1, 0,-1 }, { 0, 5,-4, 1, 0, 0 }, { 0, 3,-2,-1,-1, 0 }, 720 { 0, 3,-2, 1, 2, 0 }, { 0, 4,-4, 0, 0, 0 }, { 0, 6,-2,-2, 0, 0 }, 721 { 0, 5, 0,-3, 1, 0 }, { 0, 4,-2, 0, 2, 0 }, { 0, 2, 2,-2, 1, 0 }, 722 { 0, 0, 4, 0, 0,-2 }, { 0, 3,-1, 0, 0, 0 }, { 0, 3,-3,-1, 0, 1 }, 723 { 0, 4, 0,-2, 2, 0 }, { 0, 1,-2,-1,-1, 0 }, { 0, 2,-1, 0, 0,-1 }, 724 { 0, 4,-4, 2, 0, 0 }, { 0, 2, 1, 0, 1,-1 }, { 0, 3,-2,-1, 1, 0 }, 725 { 0, 4,-3, 0, 1, 1 }, { 0, 2, 0, 0, 3, 0 }, { 0, 6,-4, 0, 0, 0 }, 726 }; 727 728 if((int)(sizeof(DerAmp) / sizeof(double)) != NDER 729 || (int)(sizeof(DerInd) / sizeof(NVector)) != NDER) { 730 Exception e("Static arrays are corrupted"); 731 GPSTK_THROW(e); 732 } 733 734 int i,j,k,kk; 735 static const double dtr(0.01745329252); 736 737 // get amplitude, phase and frequency for each of the standard tides 738 double RealAmp[NSTD], ImagAmp[NSTD], Freq[NSTD]; 739 double phsrad, freq, phas; 740 for(i=0; i<Nin; i++) { // Nin is NSTD 741 // first find the index for this tide 742 j = stdindex[i]; 743 744 // amplitudes 745 phsrad = phs[i] * dtr; //DEG_TO_RAD; // phase in radians 746 RealAmp[i] = amp[i] * ::cos(phsrad) / ::fabs(DerAmp[j]); 747 ImagAmp[i] = amp[i] * ::sin(phsrad) / ::fabs(DerAmp[j]); 748 //LOG(INFO) << "TEST " << setw(2) << i+1 << fixed << setprecision(15) 749 // << setw(19) << amp[i] << setw(19) << phsrad << setw(19) << DerAmp[j]; 750 751 // phase and freq 752 freq = phas = 0.0; 753 for(k=0; k<6; k++) { 754 freq += DerInd[j].n[k] * freqDood[k]; 755 // not used phas += DerInd[j].n[k] * Dood[k]; 756 } 757 Freq[i] = freq; 758 759 // make 0 <= phas < 360 -- why? 760 // not used phas = ::fmod(phas,360.0); 761 // not used if(phas < 0.0) phas += 360.0; 762 763 //LOG(INFO) << "Dood " << setw(2) << i << " at index " << setw(3) << j 764 // << " (" << DerInd[j].n[0] << "," << setw(2) << DerInd[j].n[1] 765 // << "," << setw(2) << DerInd[j].n[2] << "," << setw(2) << DerInd[j].n[3] 766 // << "," << setw(2) << DerInd[j].n[4] << "," << setw(2) << DerInd[j].n[5] 767 // << ") rA iA F " << fixed << setprecision(10) << setw(13) 768 // << RealAmp[i] << " " << setw(13) << ImagAmp[i] 769 // << " " << setw(12) << Freq[i]; 770 } 771 772 // sort the frequency, and keep the key 773 int key[NSTD]; 774 for(i=0; i<Nin; ++i) key[i] = i; 775 QSort(Freq, key, Nin); 776 777 // use key to sort amplitudes 778 double tmpR[NSTD],tmpI[NSTD]; 779 for(i=0; i<Nin; ++i) { 780 tmpR[i] = RealAmp[i]; 781 tmpI[i] = ImagAmp[i]; 782 } 783 for(i=0; i<Nin; ++i) { 784 RealAmp[i] = tmpR[key[i]]; 785 ImagAmp[i] = tmpI[key[i]]; 786 } 787 788 // count the shells 789 int nl(0),nm(0),nh(0); 790 //LOG(INFO) << "TEST Sorted reamp, imamp, freq\n"; 791 for(i=0; i<Nin; i++) { // Nin is NSTD 792 //LOG(INFO) << "Sorted Dood " << setw(2) << key[i] << " rA iA F P " 793 // << fixed << setprecision(10) << setw(13) << RealAmp[i] 794 // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i]; 795 //LOG(INFO) << "TEST " << setw(2) << i+1 << fixed << setprecision(15) 796 // << setw(19)<< RealAmp[i] << setw(19) << ImagAmp[i] << setw(19) << Freq[i]; 797 if( Freq[i] < 0.5) nl++; 798 else if(Freq[i] < 1.5) nm++; 799 else if(Freq[i] < 2.5) nh++; 800 // so freq cannot be >= 2.5?? 801 } 802 //LOG(INFO) << "Shells contain " << nl << " " << nm << " " << nh; 803 804 // split arrays into vector<double> for each shell 805 vector<double> Flow,Rlow,Ilow,Fmed,Rmed,Imed,Fhi,Rhi,Ihi; 806 for(i=0; i<nl; i++) { 807 Flow.push_back(Freq[i]); 808 Rlow.push_back(RealAmp[i]); 809 Ilow.push_back(ImagAmp[i]); 810 //LOG(INFO) << "Low shell Dood " << setw(2) << key[i] << " rA iA F " 811 // << fixed << setprecision(10) << setw(13) << RealAmp[i] 812 // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i]; 813 } 814 for(i=nl; i<nl+nm; i++) { 815 Fmed.push_back(Freq[i]); 816 Rmed.push_back(RealAmp[i]); 817 Imed.push_back(ImagAmp[i]); 818 //LOG(INFO) << "Med shell Dood " << setw(2) << key[i] << " rA iA F " 819 // << fixed << setprecision(10) << setw(13) << RealAmp[i] 820 // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i]; 821 } 822 for(i=nl+nm; i<nl+nm+nh; i++) { 823 Fhi.push_back(Freq[i]); 824 Rhi.push_back(RealAmp[i]); 825 Ihi.push_back(ImagAmp[i]); 826 //LOG(INFO) << "Hi shell Dood " << setw(2) << key[i] << " rA iA F " 827 // << fixed << setprecision(10) << setw(13) << RealAmp[i] 828 // << " " << setw(13) << ImagAmp[i] << " " << setw(12) << Freq[i]; 829 } 830 831 // find splines of amp vs frequency in each shell 832 CubicSpline<double> csRlow, csIlow, csRmed, csImed, csRhi, csIhi; 833 if(nl > 0) { 834 csRlow.Initialize(Flow, Rlow); 835 csIlow.Initialize(Flow, Ilow); 836 } 837 csRmed.Initialize(Fmed, Rmed); 838 csImed.Initialize(Fmed, Imed); 839 csRhi.Initialize(Fhi, Rhi); 840 csIhi.Initialize(Fhi, Ihi); 841 842 // evaluate splines at each of the NDER waves; not all will contribute 843 int nout(0); 844 for(j=0; j<NDER; j++) { // loop over 342 derived tides 845 // this is why nout may be < NDER 846 if(DerInd[j].n[0] == 0 && nl == 0) continue; 847 848 // get phase and freq for this tide 849 freqDer[nout] = phsDer[nout] = 0.0; 850 for(k=0; k<6; k++) { 851 freqDer[nout] += DerInd[j].n[k] * freqDood[k]; 852 phsDer[nout] += DerInd[j].n[k] * Dood[k]; 853 } 854 phsDer[nout] = ::fmod(phsDer[nout],360.0); 855 if(phsDer[nout] < 0.0) phsDer[nout] += 360.0; 856 857 //LOG(INFO) << "TEST TDFRPH " 858 // << setw(3) << j+1 << fixed << setprecision(15) 859 // << setw(22) << freqDer[nout] << setw(22) << phsDer[nout] 860 // << setw(3) << DerInd[j].n[0] << setw(3) << DerInd[j].n[1] 861 // << setw(3) << DerInd[j].n[2] << setw(3) << DerInd[j].n[3] 862 // << setw(3) << DerInd[j].n[4] << setw(3) << DerInd[j].n[5]; 863 864 if( DerInd[j].n[0] == 0) phsDer[nout] += 180.0; 865 else if(DerInd[j].n[0] == 1) phsDer[nout] += 90.0; 866 867 // get amplitudes at freq 868 freq = freqDer[nout]; 869 double ramp,iamp; 870 if( DerInd[j].n[0] == 0) { 871 if(csRlow.testLimits(freq,ramp)) ramp = csRlow.Evaluate(freq); 872 if(csIlow.testLimits(freq,iamp)) iamp = csIlow.Evaluate(freq); 873 } 874 else if(DerInd[j].n[0] == 1) { 875 if(csRmed.testLimits(freq,ramp)) ramp = csRmed.Evaluate(freq); 876 if(csImed.testLimits(freq,iamp)) iamp = csImed.Evaluate(freq); 877 } 878 else if(DerInd[j].n[0] == 2) { 879 if(csRhi.testLimits(freq,ramp)) ramp = csRhi.Evaluate(freq); 880 if(csIhi.testLimits(freq,iamp)) iamp = csIhi.Evaluate(freq); 881 } 882 883 ampDer[nout] = DerAmp[j] * RSS(ramp,iamp); 884 phsDer[nout] += ::atan2(iamp,ramp)/dtr; //*RAD_TO_DEG; // TEMP 885 if(phsDer[nout] > 180.0) phsDer[nout] -= 360.0; 886 887 //LOG(INFO) << "TEST RE AM " << setw(3) << j+1 << fixed 888 // << setprecision(15) << setw(22) << ramp << setw(22) << iamp 889 // << setw(22) << ampDer[nout] << setw(22) << phsDer[nout]; 890 891 nout++; 892 } 893 894 return nout; 895 896 } // end int OceanLoadTides::deriveTides() 897 898 } // end namespace gpstk 899 //------------------------------------------------------------------------------------ 900 //------------------------------------------------------------------------------------ 901