vtecmap.cpp

来自「根据GPS观测数据」· C++ 代码 · 共 480 行 · 第 1/2 页

      vtec.push_back(data[k].VTEC);
      xtmp.push_back(range * d);       // this is probably slow
      ytmp.push_back(range * sin(bear));
      // sigma = RSS(measurement error, term ~ range = decorrelation)
      d = data[k].VTECerror * data[k].VTECerror
         + range * range * (Decorrelation/1000) * (Decorrelation/1000);
      sigma.push_back(SQRT(d));

   }  // end loop over all data

   d = ChiSqPlane(vtec,xtmp,ytmp,sigma);
   if(d < 0) {
      //if(d < -0.5) output warning: negative TEC set to 0
      d = 0.0;
   }
   
   gridpt.value = d;
}

//------------------------------------------------------------------------------------
// Chi squared calculation ... TD
// z is VTEC
double VTECMap::ChiSqPlane(vector<double>& z, vector<double>& x, vector<double>& y,
   vector<double>& sigma)
{
   int k;

   double invSigma,XinvSigma,YinvSigma,ZinvSigma;
   double s,sx,sy,sz,sxx,sxy,syy,sxz,syz;

   invSigma = 1.0/(sigma[0]*sigma[0]);
   s = invSigma;
   ZinvSigma = (z[0]-ave) * invSigma;
   sz = ZinvSigma;
   if(fittype == Linear) {
      XinvSigma = x[0] * invSigma;
      sx = XinvSigma;
      YinvSigma = y[0] * invSigma;
      sy = YinvSigma;
      sxx = x[0] * XinvSigma;
      sxy = x[0] * YinvSigma;
      syy = y[0] * YinvSigma;
      sxz = x[0] * ZinvSigma;
      syz = y[0] * ZinvSigma;
   }

   double invS2,XinvS2,YinvS2,ZinvS2;
   for(k=1; k<z.size(); k++) {
      // constant and linear
      invS2 = 1.0 / (sigma[k]*sigma[k]);
      s += invS2;
      ZinvS2 = (z[k]-ave) * invS2;
      sz += ZinvS2;
      // linear only
      if(fittype == Linear) {
         XinvS2 = x[k] * invS2;
         YinvS2 = y[k] * invS2;
         sx += XinvS2;
         sy += YinvS2;
         sxx += x[k] * XinvS2;
         sxy += x[k] * YinvS2;
         syy += y[k] * YinvS2;
         sxz += x[k] * ZinvS2;
         syz += y[k] * ZinvS2;
      }
   }

   double value;
   if(fittype == Linear) {
      double delta = sxy*(s*sxy-2*sx*sy) + sxx*sy*sy + syy*(sx*sx-s*sxx);
      value = ( sxz*(sx*syy-sxy*sy) + syz*(sxx*sy-sx*sxy)
            + sz*(sxy*sxy-sxx*syy) )/delta;
   }
   else
      value = (sz/s);
   return (ave + value);
}

//------------------------------------------------------------------------------------
void VTECMap::OutputMap(ostream& os, bool format)
{
   int i,j,k;
   os << fixed << setprecision(3);
   for(j=0; j<NumLat; j++) {
      for(i=0; i<NumLon; i++) {
         k = i * NumLat + j;
         if(format) os << grid[k].LLR[0] << " " << grid[k].LLR[1];
         os << " " << grid[k].value;
         if(format) os << endl;
      }
      os << endl;
   }
}

//------------------------------------------------------------------------------------
void MUFMap::ComputeMap(DayTime& epoch, vector<ObsData>& data)
{
   int i,k;
   double lvect1,lvect2,tmp,cosin;;
   GridData center,reflect;
   Position MUFearth;

   for(k=0; k<NumLat*NumLon; k++) {
         // Comment in the original code is:
         // "convert the lat/lon from the MUF grid
         // to XYZ positions on the surface of Earth"
         // then code uses grid[k].XYZ where MUFearth is here...
      MUFearth = grid[k].LLR;
      MUFearth[2] = MUFearth.radiusEarth();
      MUFearth.transformTo(Position::Cartesian);

      center.XYZ = (MUFearth + RefStation.xyz)*0.5;
      center.LLR = center.XYZ;
      center.LLR.transformTo(Position::Geocentric);

      reflect = center;
      reflect.LLR[2] = reflect.LLR.radiusEarth() + IonoHeight;

      ComputeGridValue(reflect, data);

      reflect.XYZ = reflect.LLR;
      reflect.XYZ.transformTo(Position::Cartesian);

      lvect1 = lvect2 = 0.0;
      for(i=0; i<3; i++) {
         tmp = MUFearth[i] - reflect.XYZ[i];
         lvect1 += tmp*tmp;
         tmp = reflect.XYZ[i] - center.XYZ[i];
         lvect2 += tmp*tmp;
      }
      cosin = SQRT(lvect2/lvect1);
      grid[k].value =
         VTECtoF0F2(0,reflect.value,epoch,reflect.LLR.longitude()) / cosin;
   }
}

//------------------------------------------------------------------------------------
// First cut at foF2 assuming constant slab thickness of 280 km and 
// TEC = 1.24e10 (foF2)^2 tau / 10^16
void F0F2Map::ComputeMap(DayTime& epoch, vector<ObsData>& data)
{
   int i,j,k;
   for(i=0; i<NumLon; i++) {
      for(j=0; j<NumLat; j++) {
         k = i * NumLat + j;
         ComputeGridValue(grid[k],data);
         grid[k].value = VTECtoF0F2(1,grid[k].value,epoch,grid[k].LLR.longitude());
      }
   }
}

//------------------------------------------------------------------------------------
double VTECMap::VTECtoF0F2(int method, double vtec, DayTime& epoch, double lon)
{
try {
   double fof2,tau,dt;
   static DayTime computeTime=DayTime::BEGINNING_OF_TIME;
   const double con[4]={0.019600827088077529, -1.549245071973630372,
      29.890989537102175433, 237.467144625490760745};

   if(method == 0) {
      tau = 280;
   }
   else if(method == 1) {
      if(epoch != computeTime) {
         computeTime = epoch;
         dt = epoch.hour()+epoch.minute()/60.;
         dt += (lon - 262.2743352)/15;
         if(dt > 24) dt -= 24;
         if(dt <  0) dt += 24;
         tau = con[0];
         for(int i=1; i<4; i++) tau = tau * dt + con[i];
      }
   }
   else {
      throw Exception("VTECtoF0F2 finds unknown method");
   }

   fof2 = SQRT(806.4 * vtec / tau);

   return fof2;
}
catch(Exception& e) {
   GPSTK_RETHROW(e);
}
}

//------------------------------------------------------------------------------------
double VTECMap::VTECError(double t, double el, double vtec)
{
   double slant = VTECErrorMultipath * VTECErrorMultipath / (1+2*t/3);
   slant += VTECErrorSat * VTECErrorSat;
   //slant += BiasError * BiasError;
   double q = Obliquity(el);
   double delta = VTECErrorCFC[0];
   for(int i=1; i<4; i++)
      delta = delta * el + VTECErrorCFC[i];
   delta *= 0.01;
   double vert = slant/(q*q);
   vert += delta * vtec * delta * vtec;
   vert = SQRT(vert);
   return vert;
}

//------------------------------------------------------------------------------------
double VTECMap::Obliquity(double el)
{
   if(el > 85) return 1.0;
   double obq = ObliqCoef[0];
   for(int i=1; i<4; i++)
      obq = obq * el + ObliqCoef[i];
   return obq;
}

//------------------------------------------------------------------------------------
void gpstk::PlaneCoefficients(double cof[3], double p1[3], double p2[3], double p3[3])
   throw(Exception)
{
   double denom=(p2[0]-p3[0])*p1[1]+(p3[0]-p1[0])*p2[1]+(p1[0]-p2[0])*p3[1];
   if(denom == 0) throw Exception("PlaneCoefficients finds singular problem");
   denom = 1.0/denom;
   cof[0] = ((p3[1]-p2[1])*p1[2]+(p1[1]-p3[1])*p2[2]+(p2[1]-p1[1])*p3[2])*denom;
   cof[1] = ((p2[0]-p3[0])*p1[2]+(p3[0]-p1[0])*p2[2]+(p1[0]-p2[0])*p3[2])*denom;
   cof[2] = ((p3[0]*p2[1]-p2[0]*p3[1])*p1[2]+(p1[0]*p3[1]-p3[0]*p1[1])*p2[2]+
             (p2[0]*p1[1]-p1[0]*p2[1])*p3[2])*denom;
}

//------------------------------------------------------------------------------------
ostream& gpstk::operator<<(ostream& os, const Station& s)
{
   os << "Station filename: " << s.filename << endl << fixed;
   os << s.xyz.printf("   %.3x m %.3y m %.3z m (ECEF)\n"
                      "   %.9a degN %.9l degE %.3r m (Geocentric)");
   return os;
}

//------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------