Ignore:
Timestamp:
Mar 29, 2011, 7:15:07 PM (13 years ago)
Author:
mervart
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/BNC/upload/bncuploadcaster.cpp

    r3174 r3182  
    2020#include "bncversion.h"
    2121#include "bncapp.h"
     22#include "bncclockrinex.h"
     23#include "bncsp3.h"
    2224
    2325using namespace std;
    … …  
    2931                                 const QString& password,
    3032                                 const QString& crdTrafo, bool  CoM,
     33                                 const QString& rnxFileName,
     34                                 const QString& sp3FileName,
    3135                                 const QString& outFileName) {
    3236
    3337  bncSettings settings;
     38
     39  connect(this, SIGNAL(newMessage(QByteArray,bool)),
     40          ((bncApp*)qApp), SLOT(slotMessage(const QByteArray,bool)));
    3441
    3542  _mountpoint = mountpoint;
    … …  
    4350  _sOpenTrial = 0;
    4451
     52  _append = Qt::CheckState(settings.value("rnxAppend").toInt()) == Qt::Checked;
     53
    4554  if (outFileName.isEmpty()) {
    4655    _outFile   = 0;
    … …  
    5059    _outFile = new QFile(outFileName);
    5160    QIODevice::OpenMode oMode;
    52     if (Qt::CheckState(settings.value("rnxAppend").toInt()) == Qt::Checked) {
     61    if (_append) {
    5362      oMode = QIODevice::WriteOnly | QIODevice::Unbuffered | QIODevice::Append;
    5463    }
    … …  
    6170    }
    6271  }
    63   connect(this, SIGNAL(newMessage(QByteArray,bool)),
    64           ((bncApp*)qApp), SLOT(slotMessage(const QByteArray,bool)));
     72
     73  // RINEX writer
     74  // ------------
     75  if ( settings.value("rnxPath").toString().isEmpty() ) {
     76    _rnx = 0;
     77  }
     78  else {
     79    QString prep  = "BNC";
     80    QString ext   = ".clk";
     81    QString path  = settings.value("rnxPath").toString();
     82    QString intr  = settings.value("rnxIntr").toString();
     83    int     sampl = settings.value("rnxSampl").toInt();
     84    _rnx = new bncClockRinex(prep, ext, path, intr, sampl);
     85  }
     86
     87  // SP3 writer
     88  // ----------
     89  if ( settings.value("sp3Path").toString().isEmpty() ) {
     90    _sp3 = 0;
     91  }
     92  else {
     93    QString prep  = "BNC";
     94    QString ext   = ".sp3";
     95    QString path  = settings.value("sp3Path").toString();
     96    QString intr  = settings.value("sp3Intr").toString();
     97    int     sampl = settings.value("sp3Sampl").toInt();
     98    _sp3 = new bncSP3(prep, ext, path, intr, sampl);
     99  }
     100
     101  // Set Transformation Parameters
     102  // -----------------------------
     103  if      (_crdTrafo == "ETRF2000") {
     104    _dx  =    0.0541;
     105    _dy  =    0.0502;
     106    _dz  =   -0.0538;
     107    _dxr =   -0.0002;
     108    _dyr =    0.0001;
     109    _dzr =   -0.0018;
     110    _ox  =  0.000891;
     111    _oy  =  0.005390;
     112    _oz  = -0.008712;
     113    _oxr =  0.000081;
     114    _oyr =  0.000490;
     115    _ozr = -0.000792;
     116    _sc  =      0.40;
     117    _scr =      0.08;
     118    _t0  =    2000.0;
     119  }
     120  else if (_crdTrafo == "NAD83") {
     121    _dx  =    0.9963;
     122    _dy  =   -1.9024;
     123    _dz  =   -0.5210;
     124    _dxr =    0.0005;
     125    _dyr =   -0.0006;
     126    _dzr =   -0.0013;
     127    _ox  =  0.025915;
     128    _oy  =  0.009426;
     129    _oz  =  0.011599;
     130    _oxr =  0.000067;
     131    _oyr = -0.000757;
     132    _ozr = -0.000051;
     133    _sc  =      0.78;
     134    _scr =     -0.10;
     135    _t0  =    1997.0;
     136  }
     137  else if (_crdTrafo == "GDA94") {
     138    _dx  =   -0.07973;
     139    _dy  =   -0.00686;
     140    _dz  =    0.03803;
     141    _dxr =    0.00225;
     142    _dyr =   -0.00062;
     143    _dzr =   -0.00056;
     144    _ox  =  0.0000351;
     145    _oy  = -0.0021211;
     146    _oz  = -0.0021411;
     147    _oxr = -0.0014707;
     148    _oyr = -0.0011443;
     149    _ozr = -0.0011701;
     150    _sc  =      6.636;
     151    _scr =      0.294;
     152    _t0  =     1994.0;
     153  }
     154  else if (_crdTrafo == "SIRGAS2000") {
     155    _dx  =   -0.0051;
     156    _dy  =   -0.0065;
     157    _dz  =   -0.0099;
     158    _dxr =    0.0000;
     159    _dyr =    0.0000;
     160    _dzr =    0.0000;
     161    _ox  =  0.000150;
     162    _oy  =  0.000020;
     163    _oz  =  0.000021;
     164    _oxr =  0.000000;
     165    _oyr =  0.000000;
     166    _ozr =  0.000000;
     167    _sc  =     0.000;
     168    _scr =     0.000;
     169    _t0  =    0000.0;
     170  }
     171  else if (_crdTrafo == "SIRGAS95") {
     172    _dx  =    0.0077;
     173    _dy  =    0.0058;
     174    _dz  =   -0.0138;
     175    _dxr =    0.0000;
     176    _dyr =    0.0000;
     177    _dzr =    0.0000;
     178    _ox  =  0.000000;
     179    _oy  =  0.000000;
     180    _oz  = -0.000030;
     181    _oxr =  0.000000;
     182    _oyr =  0.000000;
     183    _ozr =  0.000000;
     184    _sc  =     1.570;
     185    _scr =     0.000;
     186    _t0  =    0000.0;
     187  }
     188  else if (_crdTrafo == "Custom") {
     189    _dx  = settings.value("trafo_dx").toDouble();
     190    _dy  = settings.value("trafo_dy").toDouble();
     191    _dz  = settings.value("trafo_dz").toDouble();
     192    _dxr = settings.value("trafo_dxr").toDouble();
     193    _dyr = settings.value("trafo_dyr").toDouble();
     194    _dzr = settings.value("trafo_dzr").toDouble();
     195    _ox  = settings.value("trafo_ox").toDouble();
     196    _oy  = settings.value("trafo_oy").toDouble();
     197    _oz  = settings.value("trafo_oz").toDouble();
     198    _oxr = settings.value("trafo_oxr").toDouble();
     199    _oyr = settings.value("trafo_oyr").toDouble();
     200    _ozr = settings.value("trafo_ozr").toDouble();
     201    _sc  = settings.value("trafo_sc").toDouble();
     202    _scr = settings.value("trafo_scr").toDouble();
     203    _t0  = settings.value("trafo_t0").toDouble();
     204  }
    65205}
    66206
    … …  
    150290  }
    151291}
     292
     293// Encode and Upload Clocks, Orbits, and Biases
     294////////////////////////////////////////////////////////////////////////////
     295void bncUploadCaster::uploadClockOrbitBias(const QStringList& lines,
     296                           const QMap<QString, bncEphUser::t_ephPair*>& ephMap,
     297                           int year, int month, int day,
     298                           int GPSweek, double GPSweeks) {
     299
     300  this->open();
     301 
     302  struct ClockOrbit co;
     303  memset(&co, 0, sizeof(co));
     304  co.GPSEpochTime      = (int)GPSweeks;
     305  co.GLONASSEpochTime  = (int)fmod(GPSweeks, 86400.0)
     306                       + 3 * 3600 - gnumleap(year, month, day);
     307  co.ClockDataSupplied = 1;
     308  co.OrbitDataSupplied = 1;
     309  co.SatRefDatum       = DATUM_ITRF;
     310 
     311  struct Bias bias;
     312  memset(&bias, 0, sizeof(bias));
     313  bias.GPSEpochTime      = (int)GPSweeks;
     314  bias.GLONASSEpochTime  = (int)fmod(GPSweeks, 86400.0)
     315                         + 3 * 3600 - gnumleap(year, month, day);
     316 
     317  for (int ii = 0; ii < lines.size(); ii++) {
     318 
     319    QString      prn;
     320    ColumnVector xx(14); xx = 0.0;
     321    bncEphUser::t_ephPair*   pair = 0;
     322 
     323    QTextStream in(lines[ii].toAscii());
     324    in >> prn;
     325    if ( ephMap.contains(prn) ) {
     326      in >> xx(1) >> xx(2) >> xx(3) >> xx(4) >> xx(5)
     327         >> xx(6) >> xx(7) >> xx(8) >> xx(9) >> xx(10)
     328         >> xx(11) >> xx(12) >> xx(13) >> xx(14);
     329      xx(1) *= 1e3;          // x-crd
     330      xx(2) *= 1e3;          // y-crd
     331      xx(3) *= 1e3;          // z-crd
     332      xx(4) *= 1e-6;         // clk
     333      xx(5) *= 1e-6;         // rel. corr.
     334                             // xx(6), xx(7), xx(8) ... PhaseCent - CoM
     335                             // xx(9)  ... P1-C1 DCB in meters
     336                             // xx(10) ... P1-P2 DCB in meters
     337                             // xx(11) ... dT
     338      xx(12) *= 1e3;         // x-crd at time + dT
     339      xx(13) *= 1e3;         // y-crd at time + dT
     340      xx(14) *= 1e3;         // z-crd at time + dT
     341 
     342      pair = ephMap[prn];
     343    }
     344 
     345    // Use old ephemeris if the new one is too recent
     346    // ----------------------------------------------
     347    t_eph* ep = 0;
     348    if (pair) {
     349      ep = pair->last;
     350      if (pair->prev && ep &&
     351          ep->receptDateTime().secsTo(QDateTime::currentDateTime()) < 60) {
     352        ep = pair->prev;
     353      }
     354    }
     355 
     356    if (ep != 0) {
     357      struct ClockOrbit::SatData* sd = 0;
     358      if      (prn[0] == 'G') {
     359        sd = co.Sat + co.NumberOfGPSSat;
     360        ++co.NumberOfGPSSat;
     361      }
     362      else if (prn[0] == 'R') {
     363        sd = co.Sat + CLOCKORBIT_NUMGPS + co.NumberOfGLONASSSat;
     364        ++co.NumberOfGLONASSSat;
     365      }
     366      if (sd) {
     367        QString outLine;
     368        processSatellite(ep, GPSweek, GPSweeks, prn, xx, sd, outLine);
     369        this->printAscii(outLine);
     370      }
     371 
     372      struct Bias::BiasSat* biasSat = 0;
     373      if      (prn[0] == 'G') {
     374        biasSat = bias.Sat + bias.NumberOfGPSSat;
     375        ++bias.NumberOfGPSSat;
     376      }
     377      else if (prn[0] == 'R') {
     378        biasSat = bias.Sat + CLOCKORBIT_NUMGPS + bias.NumberOfGLONASSSat;
     379        ++bias.NumberOfGLONASSSat;
     380      }
     381 
     382      // Coefficient of Ionosphere-Free LC
     383      // ---------------------------------
     384      const static double a_L1_GPS =  2.54572778;
     385      const static double a_L2_GPS = -1.54572778;
     386      const static double a_L1_Glo =  2.53125000;
     387      const static double a_L2_Glo = -1.53125000;
     388 
     389      if (biasSat) {
     390        biasSat->ID = prn.mid(1).toInt();
     391        biasSat->NumberOfCodeBiases = 3;
     392        if      (prn[0] == 'G') {
     393          biasSat->Biases[0].Type = CODETYPEGPS_L1_Z;
     394          biasSat->Biases[0].Bias = - a_L2_GPS * xx(10);
     395          biasSat->Biases[1].Type = CODETYPEGPS_L1_CA;
     396          biasSat->Biases[1].Bias = - a_L2_GPS * xx(10) + xx(9);
     397          biasSat->Biases[2].Type = CODETYPEGPS_L2_Z;
     398          biasSat->Biases[2].Bias = a_L1_GPS * xx(10);
     399        }
     400        else if (prn[0] == 'R') {
     401          biasSat->Biases[0].Type = CODETYPEGLONASS_L1_P;
     402          biasSat->Biases[0].Bias = - a_L2_Glo * xx(10);
     403          biasSat->Biases[1].Type = CODETYPEGLONASS_L1_CA;
     404          biasSat->Biases[1].Bias = - a_L2_Glo * xx(10) + xx(9);
     405          biasSat->Biases[2].Type = CODETYPEGLONASS_L2_P;
     406          biasSat->Biases[2].Bias = a_L1_Glo * xx(10);
     407        }
     408      }
     409    }
     410  }
     411 
     412  if ( this->usedSocket() &&
     413       (co.NumberOfGPSSat > 0 || co.NumberOfGLONASSSat > 0) ) {
     414    char obuffer[CLOCKORBIT_BUFFERSIZE];
     415 
     416    int len = MakeClockOrbit(&co, COTYPE_AUTO, 0, obuffer, sizeof(obuffer));
     417    if (len > 0) {
     418      this->write(obuffer, len);
     419    }
     420  }
     421 
     422  if ( this->usedSocket() &&
     423       (bias.NumberOfGPSSat > 0 || bias.NumberOfGLONASSSat > 0) ) {
     424    char obuffer[CLOCKORBIT_BUFFERSIZE];
     425    int len = MakeBias(&bias, BTYPE_AUTO, 0, obuffer, sizeof(obuffer));
     426    if (len > 0) {
     427      this->write(obuffer, len);
     428    }
     429  }
     430}
     431
     432//
     433////////////////////////////////////////////////////////////////////////////
     434void bncUploadCaster::processSatellite(t_eph* eph, int GPSweek,
     435                                       double GPSweeks, const QString& prn,
     436                                       const ColumnVector& xx,
     437                                       struct ClockOrbit::SatData* sd,
     438                                       QString& outLine) {
     439
     440  const double secPerWeek = 7.0 * 86400.0;
     441
     442  ColumnVector rsw(3);
     443  ColumnVector rsw2(3);
     444  double dClk;
     445
     446  for (int ii = 1; ii <= 2; ++ii) {
     447
     448    int    GPSweek12  = GPSweek;
     449    double GPSweeks12 = GPSweeks;
     450    if (ii == 2) {
     451      GPSweeks12 += xx(11);
     452      if (GPSweeks12 > secPerWeek) {
     453        GPSweek12  += 1;
     454        GPSweeks12 -= secPerWeek;
     455      }
     456    }
     457
     458    ColumnVector xB(4);
     459    ColumnVector vv(3);
     460
     461    eph->position(GPSweek12, GPSweeks12, xB.data(), vv.data());
     462   
     463    ColumnVector xyz;
     464    if (ii == 1) {
     465      xyz = xx.Rows(1,3);
     466    }
     467    else {
     468      xyz = xx.Rows(12,14);
     469    }
     470   
     471    // Correction Center of Mass -> Antenna Phase Center
     472    // -------------------------------------------------
     473    if (! _CoM) {
     474      xyz(1) += xx(6);
     475      xyz(2) += xx(7);
     476      xyz(3) += xx(8);
     477    }
     478   
     479    if (_crdTrafo != "IGS05") {
     480      crdTrafo(GPSweek12, xyz);
     481    }
     482   
     483    ColumnVector dx = xB.Rows(1,3) - xyz ;
     484   
     485    if (ii == 1) {
     486      XYZ_to_RSW(xB.Rows(1,3), vv, dx, rsw);
     487      dClk = (xx(4) + xx(5) - xB(4)) * 299792458.0;
     488    }
     489    else {
     490      XYZ_to_RSW(xB.Rows(1,3), vv, dx, rsw2);
     491    }
     492  }
     493
     494  if (sd) {
     495    sd->ID                    = prn.mid(1).toInt();
     496    sd->IOD                   = eph->IOD();
     497    sd->Clock.DeltaA0         = dClk;
     498    sd->Orbit.DeltaRadial     = rsw(1);
     499    sd->Orbit.DeltaAlongTrack = rsw(2);
     500    sd->Orbit.DeltaCrossTrack = rsw(3);
     501    sd->Orbit.DotDeltaRadial     = (rsw2(1) - rsw(1)) / xx(11);
     502    sd->Orbit.DotDeltaAlongTrack = (rsw2(2) - rsw(2)) / xx(11);
     503    sd->Orbit.DotDeltaCrossTrack = (rsw2(3) - rsw(3)) / xx(11);
     504  }
     505
     506  outLine.sprintf("%d %.1f %s  %3d  %10.3f  %8.3f %8.3f %8.3f\n",
     507                  GPSweek, GPSweeks, eph->prn().c_str(),
     508                  eph->IOD(), dClk, rsw(1), rsw(2), rsw(3));
     509
     510  if (_rnx) {
     511    _rnx->write(GPSweek, GPSweeks, prn, xx);
     512  }
     513  if (_sp3) {
     514    _sp3->write(GPSweek, GPSweeks, prn, xx, _append);
     515  }
     516}
     517
     518// Transform Coordinates
     519////////////////////////////////////////////////////////////////////////////
     520void bncUploadCaster::crdTrafo(int GPSWeek, ColumnVector& xyz) {
     521
     522  // Current epoch minus 2000.0 in years
     523  // ------------------------------------
     524  double dt = (GPSWeek - (1042.0+6.0/7.0)) / 365.2422 * 7.0 + 2000.0 - _t0;
     525
     526  ColumnVector dx(3);
     527
     528  dx(1) = _dx + dt * _dxr;
     529  dx(2) = _dy + dt * _dyr;
     530  dx(3) = _dz + dt * _dzr;
     531
     532  static const double arcSec = 180.0 * 3600.0 / M_PI;
     533
     534  double ox = (_ox + dt * _oxr) / arcSec;
     535  double oy = (_oy + dt * _oyr) / arcSec;
     536  double oz = (_oz + dt * _ozr) / arcSec;
     537
     538  double sc = 1.0 + _sc * 1e-9 + dt * _scr * 1e-9;
     539
     540  Matrix rMat(3,3);
     541  rMat(1,1) = 1.0;
     542  rMat(1,2) = -oz;
     543  rMat(1,3) =  oy;
     544  rMat(2,1) =  oz;
     545  rMat(2,2) = 1.0;
     546  rMat(2,3) = -ox;
     547  rMat(3,1) = -oy;
     548  rMat(3,2) =  ox;
     549  rMat(3,3) = 1.0;
     550
     551  xyz = sc * rMat * xyz + dx;
     552}
Note: See TracChangeset for help on using the changeset viewer.