Context Navigation

pppParlist.cpp@ 9170

Last change on this file since 9170 was 8993, checked in by stuerze, 5 years ago
minor changes in PPP
Property svn:keywords set to `Author Date Id Rev URL;svn:eol-style=native` Property svn:mime-type set to `text/plain`
File size: 15.9 KB

Rev	Line
[7237]	1	/* -------------------------------------------------------------------------
	2	* BKG NTRIP Client
	3	* -------------------------------------------------------------------------
	4	*
	5	* Class: t_pppParlist
	6	*
	7	* Purpose: List of estimated parameters
	8	*
	9	* Author: L. Mervart
	10	*
	11	* Created: 29-Jul-2014
	12	*
	13	* Changes:
	14	*
	15	* -----------------------------------------------------------------------*/
	16
	17	#include <cmath>
	18	#include <iostream>
	19	#include <sstream>
	20	#include <iomanip>
	21	#include <algorithm>
	22	#include <newmatio.h>
	23
	24	#include "pppParlist.h"
	25	#include "pppSatObs.h"
	26	#include "pppStation.h"
	27	#include "bncutils.h"
	28	#include "bncconst.h"
	29	#include "pppClient.h"
	30
	31	using namespace BNC_PPP;
	32	using namespace std;
	33
	34	// Constructor
	35	////////////////////////////////////////////////////////////////////////////
	36	t_pppParam::t_pppParam(e_type type, const t_prn& prn, t_lc::type tLC,
	37	const vector<t_pppSatObs> obsVector) {
	38
	39	_type = type;
	40	_prn = prn;
	41	_tLC = tLC;
	42	_x0 = 0.0;
	43	_indexOld = -1;
	44	_indexNew = -1;
	45	_noise = 0.0;
	46	_ambInfo = 0;
	47
	48	switch (_type) {
	49	case crdX:
	50	_epoSpec = false;
	51	_sigma0 = OPT->_aprSigCrd[0];
	52	_noise = OPT->_noiseCrd[0];
	53	break;
	54	case crdY:
	55	_epoSpec = false;
	56	_sigma0 = OPT->_aprSigCrd[1];
	57	_noise = OPT->_noiseCrd[1];
	58	break;
	59	case crdZ:
	60	_epoSpec = false;
	61	_sigma0 = OPT->_aprSigCrd[2];
	62	_noise = OPT->_noiseCrd[2];
	63	break;
	64	case clkR:
	65	_epoSpec = true;
	66	_sigma0 = OPT->_noiseClk;
	67	break;
	68	case amb:
	69	_epoSpec = false;
	70	_sigma0 = OPT->_aprSigAmb;
	71	_ambInfo = new t_ambInfo();
	72	if (obsVector) {
	73	for (unsigned ii = 0; ii < obsVector->size(); ii++) {
	74	const t_pppSatObs* obs = obsVector->at(ii);
	75	if (obs->prn() == _prn) {
	76	double offGG = 0;
	77	if (_prn.system() == 'R' && tLC != t_lc::MW) {
	78	offGG = PPP_CLIENT->offGG();
	79	}
[8993]	80	double offGB = 0;
	81	if (_prn.system() == 'C' && tLC != t_lc::MW) {
	82	offGB = PPP_CLIENT->offGB();
	83	}
	84	_x0 = floor((obs->obsValue(tLC) - offGG - offGB - obs->cmpValue(tLC)) / obs->lambda(tLC) + 0.5);
[7237]	85	break;
	86	}
	87	}
	88	}
	89	break;
	90	case offGG:
	91	_epoSpec = true;
	92	_sigma0 = 1000.0;
	93	_x0 = PPP_CLIENT->offGG();
	94	break;
[8993]	95	case offGB:
	96	_epoSpec = true;
	97	_sigma0 = 1000.0;
	98	_x0 = PPP_CLIENT->offGB();
	99	break;
[7237]	100	case trp:
	101	_epoSpec = false;
	102	_sigma0 = OPT->_aprSigTrp;
	103	_noise = OPT->_noiseTrp;
	104	break;
[8905]	105	case ion:
	106	_epoSpec = false;
	107	_sigma0 = OPT->_aprSigIon;
	108	_noise = OPT->_noiseIon;
	109	break;
	110	case cBias1:
	111	case cBias2:
	112	_epoSpec = false;
	113	_sigma0 = OPT->_aprSigCodeBias;
	114	_noise = OPT->_noiseCodeBias;
	115	break;
	116	case pBias1:
	117	case pBias2:
	118	_epoSpec = false;
	119	_sigma0 = OPT->_aprSigPhaseBias;
	120	_noise = OPT->_noisePhaseBias;
	121	break;
[7237]	122	}
	123	}
	124
	125	// Destructor
	126	////////////////////////////////////////////////////////////////////////////
	127	t_pppParam::~t_pppParam() {
	128	delete _ambInfo;
	129	}
	130
	131	//
	132	////////////////////////////////////////////////////////////////////////////
	133	double t_pppParam::partial(const bncTime& /* epoTime /, const t_pppSatObs obs,
[8956]	134	const t_lc::type& tLC, const t_prn refPrn) const {
[7237]	135
	136	// Special Case - Melbourne-Wuebbena
	137	// ---------------------------------
	138	if (tLC == t_lc::MW && _type != amb) {
	139	return 0.0;
	140	}
	141
	142	const t_pppStation* sta = PPP_CLIENT->staRover();
[8905]	143	ColumnVector rhoV = sta->xyzApr() - obs->xc().Rows(1,3);
[7237]	144
[8905]	145	map<t_frequency::type, double> codeCoeff;
	146	map<t_frequency::type, double> phaseCoeff;
	147	map<t_frequency::type, double> ionoCoeff;
	148	obs->lcCoeff(tLC, codeCoeff, phaseCoeff, ionoCoeff);
	149
[7237]	150	switch (_type) {
	151	case crdX:
[8961]	152	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[8905]	153	return (sta->xyzApr()[0] - obs->xc()[0]) / rhoV.NormFrobenius();
[7237]	154	case crdY:
[8961]	155	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[8905]	156	return (sta->xyzApr()[1] - obs->xc()[1]) / rhoV.NormFrobenius();
[7237]	157	case crdZ:
[8961]	158	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[8905]	159	return (sta->xyzApr()[2] - obs->xc()[2]) / rhoV.NormFrobenius();
[7237]	160	case clkR:
[8961]	161	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[7237]	162	return 1.0;
	163	case offGG:
[8961]	164	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[7237]	165	return (obs->prn().system() == 'R') ? 1.0 : 0.0;
[8993]	166	case offGB:
	167	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
	168	return (obs->prn().system() == 'C') ? 1.0 : 0.0;
[7237]	169	case amb:
[8961]	170	if (tLC == t_lc::GIM \|\| tLC == t_lc::Tz0) {return 0.0;}
[8905]	171	else if ((OPT->_obsModelType == OPT->IF) \|\|
	172	(OPT->_obsModelType == OPT->PPPRTK) \|\|
	173	(OPT->_obsModelType == OPT->UncombPPP) \|\|
	174	(OPT->_obsModelType == OPT->DCMcodeBias && !obs->isReference()) \|\|
	175	(OPT->_obsModelType == OPT->DCMphaseBias && !obs->isReference()) ) {
	176	if (obs->prn() == _prn) {
	177	if (tLC == _tLC) {
	178	return (obs->lambda(tLC));
	179	}
	180	else if (tLC == t_lc::lIF && _tLC == t_lc::MW) {
	181	return obs->lambda(t_lc::lIF) * obs->lambda(t_lc::MW) / obs->lambda(t_lc::l2);
	182	}
	183	else {
	184	if (_tLC == t_lc::l1) {
	185	return obs->lambda(t_lc::l1) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
	186	}
	187	else if (_tLC == t_lc::l2) {
	188	return obs->lambda(t_lc::l2) * phaseCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
	189	}
	190	}
	191	}
	192	}
	193	break;
	194	case trp:
[8961]	195	if (tLC == t_lc::GIM) {
	196	return 0.0;
	197	}
	198	else if (tLC == t_lc::Tz0) {
	199	return 1.0;
	200	}
	201	else {
	202	return 1.0 / sin(obs->eleSat());
	203	}
[8905]	204	case ion:
[7237]	205	if (obs->prn() == _prn) {
[8905]	206	if (tLC == t_lc::c1) {
	207	return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c1)];
[7237]	208	}
[8905]	209	else if (tLC == t_lc::c2) {
	210	return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::c2)];
[7237]	211	}
[8905]	212	else if (tLC == t_lc::l1) {
	213	return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l1)];
[7237]	214	}
[8905]	215	else if (tLC == t_lc::l2) {
	216	return ionoCoeff[t_lc::toFreq(obs->prn().system(),t_lc::l2)];
	217	}
	218	else if (tLC == t_lc::GIM) {
	219	return -1.0;
	220	}
[7237]	221	}
[8956]	222	if (tLC == t_lc::GIM && _prn == refPrn) {
[8961]	223	return 1.0;
[8905]	224	}
	225	break;
	226	case cBias1:
	227	if (tLC == t_lc::c1) {
[8961]	228	return 1.0;
[8905]	229	}
	230	else {
	231	return 0.0;
	232	}
	233	break;
	234	case cBias2:
	235	if (tLC == t_lc::c2) {
[8961]	236	return 1.0;
[8905]	237	}
	238	else {
	239	return 0.0;
	240	}
	241	break;
	242	case pBias1:
	243	if (tLC == t_lc::l1) {
[8961]	244	return 1.0;
[8905]	245	}
	246	else {
	247	return 0.0;
	248	}
	249	break;
	250	case pBias2:
	251	if (tLC == t_lc::l2) {
[8961]	252	return 1.0;
[8905]	253	}
	254	else {
	255	return 0.0;
	256	}
[7237]	257	}
	258	return 0.0;
	259	}
	260
	261	//
	262	////////////////////////////////////////////////////////////////////////////
	263	string t_pppParam::toString() const {
	264	stringstream ss;
	265	switch (_type) {
	266	case crdX:
	267	ss << "CRD_X";
	268	break;
	269	case crdY:
	270	ss << "CRD_Y";
	271	break;
	272	case crdZ:
	273	ss << "CRD_Z";
	274	break;
	275	case clkR:
[8905]	276	ss << "REC_CLK ";
[7237]	277	break;
	278	case offGG:
	279	ss << "OGG ";
	280	break;
[8993]	281	case offGB:
	282	ss << "OGB ";
	283	break;
[7237]	284	case trp:
	285	ss << "TRP ";
	286	break;
[8905]	287	case amb:
	288	ss << "AMB " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
	289	break;
	290	case ion:
	291	ss << "ION " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << _prn.toString();
	292	break;
	293	case cBias1:
	294	case cBias2:
	295	case pBias1:
	296	case pBias2:
	297	ss << "BIAS " << left << setw(3) << t_lc::toString(_tLC) << right << ' ' << "REC";
	298	break;
[7237]	299	}
	300	return ss.str();
	301	}
	302
	303	// Constructor
	304	////////////////////////////////////////////////////////////////////////////
	305	t_pppParlist::t_pppParlist() {
	306	}
	307
	308	// Destructor
	309	////////////////////////////////////////////////////////////////////////////
	310	t_pppParlist::~t_pppParlist() {
	311	for (unsigned ii = 0; ii < _params.size(); ii++) {
	312	delete _params[ii];
	313	}
	314	}
	315
	316	//
	317	////////////////////////////////////////////////////////////////////////////
[8905]	318	t_irc t_pppParlist::set(const bncTime& epoTime,
[8956]	319	const std::vector<t_pppSatObs*>& obsVector,
	320	const QMap<char, t_pppRefSat*>& refSatMap) {
[7237]	321
	322	// Remove some Parameters
	323	// ----------------------
	324	vector<t_pppParam*>::iterator it = _params.begin();
	325	while (it != _params.end()) {
	326	t_pppParam* par = *it;
	327
	328	bool remove = false;
	329
	330	if (par->epoSpec()) {
	331	remove = true;
	332	}
	333
[8956]	334	else if (par->type() == t_pppParam::crdX \|\|
[7988]	335	par->type() == t_pppParam::crdY \|\|
	336	par->type() == t_pppParam::crdZ) {
[7973]	337	if (par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) {
[7237]	338	remove = true;
	339	}
	340	}
[8956]	341	else if (par->type() == t_pppParam::amb) {
	342	char system = par->prn().system();
	343	t_prn refPrn = t_prn();
	344	if (OPT->_refSatRequired) {
	345	refPrn = (refSatMap[system])->prn();
	346	}
	347	if ((par->lastObsTime().valid() && (epoTime - par->lastObsTime() > 60.0)) \|\|
	348	( refPrn == par->prn())) {
	349	remove = true;
	350	}
	351	}
[7237]	352
	353	if (remove) {
	354	delete par;
	355	it = _params.erase(it);
	356	}
	357	else {
	358	++it;
	359	}
	360	}
	361
	362	// Check whether parameters have observations
	363	// ------------------------------------------
	364	for (unsigned ii = 0; ii < _params.size(); ii++) {
	365	t_pppParam* par = _params[ii];
	366	if (par->prn() == 0) {
	367	par->setLastObsTime(epoTime);
	368	if (par->firstObsTime().undef()) {
	369	par->setFirstObsTime(epoTime);
	370	}
	371	}
	372	else {
	373	for (unsigned jj = 0; jj < obsVector.size(); jj++) {
	374	const t_pppSatObs* satObs = obsVector[jj];
	375	if (satObs->prn() == par->prn()) {
	376	par->setLastObsTime(epoTime);
	377	if (par->firstObsTime().undef()) {
	378	par->setFirstObsTime(epoTime);
	379	}
	380	break;
	381	}
	382	}
	383	}
	384	}
	385
	386	// Required Set of Parameters
	387	// --------------------------
	388	vector<t_pppParam*> required;
	389
	390	// Coordinates
	391	// -----------
	392	required.push_back(new t_pppParam(t_pppParam::crdX, t_prn(), t_lc::dummy));
	393	required.push_back(new t_pppParam(t_pppParam::crdY, t_prn(), t_lc::dummy));
	394	required.push_back(new t_pppParam(t_pppParam::crdZ, t_prn(), t_lc::dummy));
	395
	396	// Receiver Clock
	397	// --------------
	398	required.push_back(new t_pppParam(t_pppParam::clkR, t_prn(), t_lc::dummy));
	399
[8905]	400	// GPS-GLONASS Clock Offset
[7237]	401	// ------------------------
	402	if (OPT->useSystem('R')) {
	403	required.push_back(new t_pppParam(t_pppParam::offGG, t_prn(), t_lc::dummy));
	404	}
	405
[8993]	406	// GPS-BDS Clock Offset
	407	// ------------------------
	408	if (OPT->useSystem('C')) {
	409	required.push_back(new t_pppParam(t_pppParam::offGB, t_prn(), t_lc::dummy));
	410	}
	411
	412
[7237]	413	// Troposphere
	414	// -----------
	415	if (OPT->estTrp()) {
	416	required.push_back(new t_pppParam(t_pppParam::trp, t_prn(), t_lc::dummy));
	417	}
	418
[8905]	419	// Ionosphere
	420	// ----------
	421	if (OPT->_obsModelType == OPT->UncombPPP \|\|
	422	OPT->_obsModelType == OPT->DCMcodeBias \|\|
	423	OPT->_obsModelType == OPT->DCMphaseBias ) {
	424	for (unsigned jj = 0; jj < obsVector.size(); jj++) {
	425	const t_pppSatObs* satObs = obsVector[jj];
	426	required.push_back(new t_pppParam(t_pppParam::ion, satObs->prn(), t_lc::dummy));
	427	}
	428	}
	429
[7237]	430	// Ambiguities
	431	// -----------
	432	for (unsigned jj = 0; jj < obsVector.size(); jj++) {
	433	const t_pppSatObs* satObs = obsVector[jj];
[8905]	434	if ((OPT->_obsModelType == OPT->IF) \|\|
	435	(OPT->_obsModelType == OPT->PPPRTK) \|\|
	436	(OPT->_obsModelType == OPT->UncombPPP) \|\|
	437	(OPT->_obsModelType == OPT->DCMcodeBias && !satObs->isReference()) \|\|
	438	(OPT->_obsModelType == OPT->DCMphaseBias && !satObs->isReference()) ) {
	439	const vector<t_lc::type>& ambLCs = OPT->ambLCs(satObs->prn().system());
	440	for (unsigned ii = 0; ii < ambLCs.size(); ii++) {
	441	required.push_back(new t_pppParam(t_pppParam::amb, satObs->prn(), ambLCs[ii], &obsVector));
	442	}
[7237]	443	}
	444	}
	445
[8905]	446	// Receiver Code Biases
	447	// --------------------
	448	if (OPT->_obsModelType == OPT->DCMcodeBias) {
	449	required.push_back(new t_pppParam(t_pppParam::cBias1, t_prn(), t_lc::c1));
	450	required.push_back(new t_pppParam(t_pppParam::cBias2, t_prn(), t_lc::c2));
	451	}
	452
	453	// Receiver Phase Biases
	454	// ---------------------
	455	if ((OPT->_obsModelType == OPT->DCMphaseBias) \|\|
	456	(OPT->_obsModelType == OPT->PPPRTK) ) {
	457	required.push_back(new t_pppParam(t_pppParam::pBias1, t_prn(), t_lc::l1));
	458	required.push_back(new t_pppParam(t_pppParam::pBias2, t_prn(), t_lc::l2));
	459	}
	460
[7237]	461	// Check if all required parameters are present
	462	// --------------------------------------------
	463	for (unsigned ii = 0; ii < required.size(); ii++) {
	464	t_pppParam* parReq = required[ii];
	465
	466	bool found = false;
	467	for (unsigned jj = 0; jj < _params.size(); jj++) {
	468	t_pppParam* parOld = _params[jj];
	469	if (parOld->isEqual(parReq)) {
	470	found = true;
	471	break;
	472	}
	473	}
	474	if (found) {
	475	delete parReq;
	476	}
	477	else {
	478	_params.push_back(parReq);
	479	}
	480	}
	481
	482	// Set Parameter Indices
	483	// ---------------------
	484	sort(_params.begin(), _params.end(), t_pppParam::sortFunction);
	485
	486	for (unsigned ii = 0; ii < _params.size(); ii++) {
	487	t_pppParam* par = _params[ii];
	488	par->setIndex(ii);
	489	for (unsigned jj = 0; jj < obsVector.size(); jj++) {
	490	const t_pppSatObs* satObs = obsVector[jj];
	491	if (satObs->prn() == par->prn()) {
	492	par->setAmbEleSat(satObs->eleSat());
	493	par->stepAmbNumEpo();
	494	}
	495	}
	496	}
	497
	498	return success;
	499	}
	500
	501	//
	502	////////////////////////////////////////////////////////////////////////////
	503	void t_pppParlist::printResult(const bncTime& epoTime, const SymmetricMatrix& QQ,
	504	const ColumnVector& xx) const {
	505
	506	string epoTimeStr = string(epoTime);
	507	const t_pppStation* sta = PPP_CLIENT->staRover();
	508
	509	LOG << endl;
	510
	511	t_pppParam* parX = 0;
	512	t_pppParam* parY = 0;
	513	t_pppParam* parZ = 0;
	514	for (unsigned ii = 0; ii < _params.size(); ii++) {
	515	t_pppParam* par = _params[ii];
	516	if (par->type() == t_pppParam::crdX) {
	517	parX = par;
	518	}
	519	else if (par->type() == t_pppParam::crdY) {
	520	parY = par;
	521	}
	522	else if (par->type() == t_pppParam::crdZ) {
	523	parZ = par;
	524	}
	525	else {
	526	int ind = par->indexNew();
	527	double apr = (par->type() == t_pppParam::trp) ?
	528	t_tropo::delay_saast(sta->xyzApr(), M_PI/2.0) : par->x0();
	529	LOG << epoTimeStr << ' ' << par->toString() << ' '
	530	<< setw(10) << setprecision(4) << apr << ' '
	531	<< showpos << setw(10) << setprecision(4) << xx[ind] << noshowpos << " +- "
	532	<< setw(8) << setprecision(4) << sqrt(QQ[ind][ind]);
	533	if (par->type() == t_pppParam::amb) {
	534	LOG << " el = " << setw(6) << setprecision(2) << par->ambEleSat() * 180.0 / M_PI
	535	<< " epo = " << setw(4) << par->ambNumEpo();
	536	}
	537	LOG << endl;
	538	}
	539	}
	540
	541	if (parX && parY && parZ) {
	542
	543	ColumnVector xyz(3);
	544	xyz[0] = xx[parX->indexNew()];
	545	xyz[1] = xx[parY->indexNew()];
	546	xyz[2] = xx[parZ->indexNew()];
	547
	548	ColumnVector neu(3);
	549	xyz2neu(sta->ellApr().data(), xyz.data(), neu.data());
	550
	551	SymmetricMatrix QQxyz = QQ.SymSubMatrix(1,3);
	552
	553	SymmetricMatrix QQneu(3);
	554	covariXYZ_NEU(QQxyz, sta->ellApr().data(), QQneu);
	555
	556	LOG << epoTimeStr << ' ' << sta->name()
	557	<< " X = " << setprecision(4) << sta->xyzApr()[0] + xyz[0] << " +- "
	558	<< setprecision(4) << sqrt(QQxyz[0][0])
	559
	560	<< " Y = " << setprecision(4) << sta->xyzApr()[1] + xyz[1] << " +- "
	561	<< setprecision(4) << sqrt(QQxyz[1][1])
	562
	563	<< " Z = " << setprecision(4) << sta->xyzApr()[2] + xyz[2] << " +- "
	564	<< setprecision(4) << sqrt(QQxyz[2][2])
	565
	566	<< " dN = " << setprecision(4) << neu[0] << " +- "
	567	<< setprecision(4) << sqrt(QQneu[0][0])
	568
	569	<< " dE = " << setprecision(4) << neu[1] << " +- "
	570	<< setprecision(4) << sqrt(QQneu[1][1])
	571
	572	<< " dU = " << setprecision(4) << neu[2] << " +- "
	573	<< setprecision(4) << sqrt(QQneu[2][2])
	574
	575	<< endl;
	576	}
	577	}
	578

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format