preprocess.cpp

综合利用三次法和blewitt法检测与修复周跳

			{
				if ( calc[j].token == b )
				{
					tre[j-3].Bw = comb[j].Bw - 3 * comb[ nextPnt[ k%3 ] ].Bw + 3 * comb[ nextPnt[ (k-1)%3 ] ].Bw - comb[ nextPnt[ (k-2)%3 ] ].Bw;
					tre[j-3].LP = comb[j].LP - 3 * comb[ nextPnt[ k%3 ] ].LP + 3 * comb[ nextPnt[ (k-1)%3 ] ].LP - comb[ nextPnt[ (k-2)%3 ] ].LP;
					BwDiff1[j-1] = comb[j].Bw - comb[ nextPnt[ k%3 ] ].Bw;
					BwDiff2[j-2] = comb[j].Bw - 2 * comb[ nextPnt[ k%3 ] ].Bw + comb[ nextPnt[ (k-1)%3 ] ].Bw;	
					nextPnt[ (k-2)%3 ] = j;
					k++;
				}
				else if ( calc[j].token > b )
				{
					break;
				}

				else
				{
				}
			}
			i = j - 1;
		}

		else
		{}


	}

	return HIGHERDIFFCALCOK;
}




static int HigherDiff_DelOutlier( CALCDAT calc[],		const COMBI comb[],
								  TREDIFF tre[],
		  const OBSDAT obs[],	double BwDiff1[],		double BwDiff2[],
						  const long int N,		const double TolZero,	const double TolScale )
{

	long int i, j, k, l, m;
	long int ncut = 0;
	int b;

	long int prePnt[3];
	long int nextPnt[3];

	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >= 0 )
		{
			for ( j = i + 1, k = 0; j < N, k < 3; j++ )
			{
				if ( calc[j].token == b )
					k++;
			}

			for ( ; j < N; j++ ) 
			{
				if ( calc[j].token == b && fabs( tre[ j-3 ].Bw ) > TolZero )
				{
					for ( k = j + 1, l = 0; k < N, l < 3;  k++ )
					{
						if ( calc[k].token == b )
							nextPnt[ l++ ] = k;

						else if ( calc[k].token > b )
							break;
					}

					if ( l == 0 )
					{
							calc[j].token = DEL_HIGHERDIFF_OUTLIER;
							ncut++;

							fprintf( fpOutPut, "%d\t", obs[j].epoch );

					}
					else
					{
						for ( k = j - 1, m = 2; k >= i, m > -1; k-- )
						{
							if ( calc[k].token == b )
								prePnt[ m-- ] = k;
						}

						if ( l == 1 )
						{
							calc[j].token = DEL_HIGHERDIFF_OUTLIER;
							ncut++;

							fprintf( fpOutPut, "%d\t", obs[j].epoch );
							BwDiff1[ nextPnt[0] - 1 ] = comb[ nextPnt[0] ].Bw - comb[ prePnt[2] ].Bw;
							BwDiff2[ nextPnt[0] - 2 ] = comb[ nextPnt[0] ].Bw - 2 * comb[ prePnt[2] ].Bw + comb[ prePnt[1] ].Bw;

							tre[ nextPnt[0] - 3 ].Bw = comb[ nextPnt[0] ].Bw - 3 * comb[ prePnt[2] ].Bw + 3 * comb[ prePnt[1] ].Bw - comb[ prePnt[0] ].Bw;
							tre[ nextPnt[0] - 3 ].LP = comb[ nextPnt[0] ].LP - 3 * comb[ prePnt[2] ].LP + 3 * comb[ prePnt[1] ].LP - comb[ prePnt[0] ].LP;
						}

						if ( l == 2 )
						{
							calc[j].token = DEL_HIGHERDIFF_OUTLIER;
							ncut++;
							fprintf( fpOutPut, "%d\t", obs[j].epoch );
							BwDiff1[ nextPnt[0] - 1 ] = comb[ nextPnt[0] ].Bw - comb[ prePnt[2] ].Bw;
							BwDiff2[ nextPnt[0] - 2 ] = comb[ nextPnt[0] ].Bw - 2 * comb[ prePnt[2] ].Bw + comb[ prePnt[1] ].Bw;
							BwDiff2[ nextPnt[1] - 2 ] = comb[ nextPnt[1] ].Bw - 2 * comb[ nextPnt[0] ].Bw + comb[ prePnt[2] ].Bw;

							tre[ nextPnt[0] - 3 ].Bw = comb[ nextPnt[0] ].Bw - 3 * comb[ prePnt[2] ].Bw + 3 * comb[ prePnt[1] ].Bw - comb[ prePnt[0] ].Bw;
							tre[ nextPnt[0] - 3 ].LP = comb[ nextPnt[0] ].LP - 3 * comb[ prePnt[2] ].LP + 3 * comb[ prePnt[1] ].LP - comb[ prePnt[0] ].LP;
							tre[ nextPnt[1] - 3 ].Bw = comb[ nextPnt[1] ].Bw - 3 * comb[ nextPnt[0] ].Bw + 3 * comb[ prePnt[2] ].Bw - comb[ prePnt[1] ].Bw;
							tre[ nextPnt[1] - 3 ].LP = comb[ nextPnt[1] ].LP - 3 * comb[ nextPnt[0] ].LP + 3 * comb[ prePnt[2] ].LP - comb[ prePnt[1] ].LP;
						}

						if ( l == 3 )
						{
							if ( ( fabs( tre[ nextPnt[0] - 3 ].Bw / tre[ j-3 ].Bw + 2.0 ) < TolScale ) &&
								 ( fabs( tre[ nextPnt[0] - 3 ].Bw / tre[ nextPnt[1] - 3 ].Bw + 2.0 ) < TolScale ) &&
								 ( fabs( tre[ nextPnt[2] - 3 ].Bw / tre[ nextPnt[1] - 3 ].Bw ) < TolScale ) &&
								 ( fabs( tre[ nextPnt[2] - 3 ].Bw ) < TolZero ) )
							{
								j = nextPnt[2];
							}
							else
							{
								while ( 1 )
								{
									calc[j].token = DEL_HIGHERDIFF_OUTLIER;
									ncut++;
									fprintf( fpOutPut, "%d\t", obs[j].epoch );
									BwDiff1[ nextPnt[0] - 1 ] = comb[ nextPnt[0] ].Bw - comb[ prePnt[2] ].Bw;
									BwDiff2[ nextPnt[0] - 2 ] = comb[ nextPnt[0] ].Bw - 2 * comb[ prePnt[2] ].Bw + comb[ prePnt[1] ].Bw;
									BwDiff2[ nextPnt[1] - 2 ] = comb[ nextPnt[1] ].Bw - 2 * comb[ nextPnt[0] ].Bw + comb[ prePnt[2] ].Bw;
									tre[ nextPnt[0] - 3 ].Bw = comb[ nextPnt[0] ].Bw - 3 * comb[ prePnt[2] ].Bw + 3 * comb[ prePnt[1] ].Bw - comb[ prePnt[0] ].Bw;
									tre[ nextPnt[0] - 3 ].LP = comb[ nextPnt[0] ].LP - 3 * comb[ prePnt[2] ].LP + 3 * comb[ prePnt[1] ].LP - comb[ prePnt[0] ].LP;
									tre[ nextPnt[1] - 3 ].Bw = comb[ nextPnt[1] ].Bw - 3 * comb[ nextPnt[0] ].Bw + 3 * comb[ prePnt[2] ].Bw - comb[ prePnt[1] ].Bw;
									tre[ nextPnt[1] - 3 ].LP = comb[ nextPnt[1] ].LP - 3 * comb[ nextPnt[0] ].LP + 3 * comb[ prePnt[2] ].LP - comb[ prePnt[1] ].LP;
									tre[ nextPnt[2] - 3 ].Bw = comb[ nextPnt[2] ].Bw - 3 * comb[ nextPnt[1] ].Bw + 3 * comb[ nextPnt[0] ].Bw - comb[ prePnt[2] ].Bw;
									tre[ nextPnt[2] - 3 ].LP = comb[ nextPnt[2] ].LP - 3 * comb[ nextPnt[1] ].LP + 3 * comb[ nextPnt[0] ].LP - comb[ prePnt[2] ].LP;
									j = nextPnt[0];
									if ( fabs( tre[ j-3 ].Bw ) <= TolZero )
										break;
									else
									{
										nextPnt[0] = nextPnt[1];
										nextPnt[1] = nextPnt[2];

										for ( k = nextPnt[2] + 1, l = 2; k < N, l < 3; k++ )
										{
											if ( calc[k].token == b )
											{
												nextPnt[ l++ ] = k;
												break;
											}
											else if ( calc[k].token > b )
											{
												break;
											}
										}

										if ( ( l >= 3 ) &&
											 ( fabs( tre[ nextPnt[0] - 3 ].Bw / tre[ j-3 ].Bw + 2.0 ) < TolScale ) &&
											 ( fabs( tre[ nextPnt[0] - 3 ].Bw / tre[ nextPnt[1] - 3 ].Bw + 2.0 ) < TolScale ) &&
											 ( fabs( tre[ nextPnt[2] - 3 ].Bw / tre[ nextPnt[1] - 3 ].Bw ) < TolScale ) &&
											 ( fabs( tre[ nextPnt[2] - 3 ].Bw ) < TolScale ) )
										{
											j = nextPnt[2];
											break;
										}
										
										else if ( l < 3 )
										{
											j = j - 1;
											break;
										}
									}
								}
							}
						}
					}
				}

				else if ( calc[j].token > b )
				{
					break;
				}
			}
			
			i = j - 1;
		}
	}

	fprintf( fpOutPut, "\n共探测野值 %ld 个\n剔除野值\n", ncut );
	fprintf( fpOutPut, "  历元编码\tb1\tb2\tbw\tbDiff1\tbDiff2\tbDiff3\t记号\n" );

	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >= 0 )
		{
			for ( j = i, l = 0; j < N; j++ )
			{
				if ( calc[j].token == b )
				{
					if ( l == 0 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%s\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, strNOVALUE, strNOVALUE, strNOVALUE, calc[j].token );
					if ( l == 1 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], strNOVALUE, strNOVALUE, calc[j].token );
					if ( l == 2 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], BwDiff2[j-2], strNOVALUE, calc[j].token );
					if ( l >= 3 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], BwDiff2[j-2], tre[j-3].Bw, calc[j].token );
					l++;
				}
				else if ( calc[j].token > b )
				{
					break;
				}
				else
				{
					fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%s\t%s\t%d\n", obs[j].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[j].token );
				}
			}
			i = j - 1;
		}
		else
		{
			fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%s\t%s\t%d\n", obs[i].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[i].token );
		}
	}

	fflush( fpOutPut );

	return HIGHERDIFFDELOK;
}




static int HigherDiff_CorrCycleslip( CALCDAT calc[],		COMBI comb[],
									 TREDIFF tre[],
									 const OBSDAT obs[],	double BwDiff1[],	double BwDiff2[],
									 const long int N,		const double TolZero )
{

	fprintf( fpOutPut, "检测上部数据序列，依据三差值进行周跳探测\n  周跳发生历元编码\tbw跳值\tb1跳值\tb2跳值\n" );

	long int i, j, k, l;
	int b;

	long int nextPnt[2];
	long int prePnt[3];

	double delta0, delta1, delta2, deltaLP;

	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >= 0 )
		{
			for ( j = i + 1, k = 0; j < N, k < 3; j++ )
			{
				if ( calc[j].token == b )
					k++;

				else if ( calc[j].token > b )//**
					break;//**
			}

			if ( k >= 3 )
			{
				for ( j--; j < N; j++ )
				{
					if ( calc[j].token == b && fabs( tre[ j-3 ].Bw ) > TolZero )
					{
						for ( k = j - 1, l = 2; k >= i, l > -1; k-- )
						{
							if ( calc[k].token == b )
								prePnt[ l-- ] = k;
						}

						delta0 = tre[ j-3 ].Bw;
						deltaLP = tre[ j-3 ].LP;
						
						for ( k = j + 1, l = 0; k < N, l < 2; k++  )
						{
							if ( calc[k].token == b )
							{
								delta0 += ( double )( l * 2 - 1 ) * tre[ k-3 ].Bw;
								deltaLP += ( double )( l * 2 - 1 ) * tre[ k-3 ].LP;
								nextPnt[ l++ ] = k;
							}

						
						}

						delta0 /= 4.0;													//delta0 = RoundOff( delta0 );
						deltaLP /= 4.0;
						delta2 = ( lambda1 * delta0 - deltaLP ) / lambdai;				//delta2 = RoundOff( delta2 );
						delta1 = delta0 + delta2;


						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\n", obs[j].epoch, delta0, delta1, delta2 );


						for ( k = j; k < N; k++ )
						{
							comb[k].Bw -= delta0;
							comb[k].LP -= deltaLP;
							calc[k].cb1 -= delta1;
							calc[k].cb2 -= delta2;
						}

						BwDiff1[ j-1 ] = comb[j].Bw - comb[ prePnt[2] ].Bw;
						BwDiff2[ j-2 ] = comb[j].Bw - 2 * comb[ prePnt[2] ].Bw + comb[ prePnt[1] ].Bw;
						BwDiff2[ nextPnt[0] - 2 ] = comb[ nextPnt[0] ].Bw - 2 * comb[j].Bw + comb[ prePnt[2] ].Bw;



						tre[ j-3 ].Bw = comb[ j ].Bw - 3 * comb[ prePnt[2] ].Bw + 3 * comb[ prePnt[1] ].Bw - comb[ prePnt[0] ].Bw;
						tre[ nextPnt[0] - 3 ].Bw = comb[ nextPnt[0] ].Bw - 3 * comb[j].Bw + 3 * comb[ prePnt[2] ].Bw - comb[ prePnt[1] ].Bw;
						tre[ nextPnt[1] - 3 ].Bw = comb[ nextPnt[1] ].Bw - 3 * comb[ nextPnt[0] ].Bw + 3 * comb[j].Bw - comb[ prePnt[2] ].Bw;



						j = nextPnt[1];
					}
					
					else if ( calc[j].token > b )
					{
						break;
					}
				}
			}

			i = j -1;
		}
	}


	fprintf( fpOutPut, "修复周跳\n  历元编码\tb1\tb2\tbw\tbDiff1\tbDiff2\tbDiff3\t记号\n" );
	
	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >= 0 )
		{
			for ( j = i, l = 0; j < N; j++ )
			{
				if ( calc[j].token == b )
				{
					if ( l == 0 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%s\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, strNOVALUE, strNOVALUE, strNOVALUE, calc[j].token );
					if ( l == 1 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], strNOVALUE, strNOVALUE, calc[j].token );
					if ( l == 2 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%s\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], BwDiff2[j-2], strNOVALUE, calc[j].token );
					if ( l >= 3 )
						fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%d\n", obs[j].epoch, calc[j].cb1, calc[j].cb2, comb[j].Bw, BwDiff1[j-1], BwDiff2[j-2], tre[j-3].Bw, calc[j].token );
					l++;
				}
				else if ( calc[j].token > b )
				{
					break;
				}
				else
				{
					fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%s\t%s\t%d\n", obs[j].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[j].token );
				}
			}
			i = j - 1;
		}
		else
		{
			fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%s\t%st%d\n", obs[i].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[i].token );
		}
	}

	fflush( fpOutPut );


	return HIGHERDIFFCORRCSOK;
}




static int TurboEdit( const OBSDAT obs[],		CALCDAT calc[],
					  COMBI comb[],				const long int N,
					  const int TEMinArc,
					  const double MultiFactor,	const double TolStep )
{
	int value;


	fprintf( fpOutPut, "\n解算第三步，使用TurboEdit方法剔除野值、探测和修复周跳\n" );


	if ( ( value = TurboEdit_Detect( calc, comb, N,

									 obs,

									 TEMinArc, MultiFactor, TolStep ) ) != TURBOEDITDETECTOK )
		return value;

	if ( ( value = TurboEdit_PolyFit( obs, calc, comb, N ) ) != TURBOEDITPOLYFITOK )
		return value;

	if ( ( value = TurboEdit_Repair( calc, comb,

									 obs,

									 N ) ) != TURBOEDITREPAIROK )
		return value;
	return value;

	
}




static int TurboEdit_Detect( CALCDAT calc[],	const COMBI comb[],
							 const long int N,

							 const OBSDAT obs[],

							 const int TEMinArc,	const double MultiFactor,	const double TolStep )
{



	fprintf( fpOutPut, "野值（太短弧段历元）出现的历元编码\n" );



	long int i, j, k, l;
	long int nextPnt;

	long int i1st/* = NOPNT*/;
	long int i2nd/* = NOPNT*/;
	int b;
	int jb;

	long int n;
	long int m;
	long int ncut = 0;

	double mean/* = 0.0*/;
	double sigma/* = 0.0*/;

	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >=0 )
		{
			i1st = NOPNT;
			i2nd = NOPNT;
			jb = 0;

			n = 0;
			m = 0;

			mean = 0.0;
			sigma = 0.0;
			
			for ( j = i; j < N; j++ )
			{
				if ( calc[j].token == b )
				{
					if ( i1st < 0 )
					{
						i1st = j;