satcoorddoc.cpp

gps 卫星位置的计算

// SatCoordDoc.cpp : implementation of the CSatCoordDoc class
//

#include "stdafx.h"
#include "SatCoord.h"

#include "SatCoordView.h"
#include "SatCoordDoc.h"
#include "coordtrans.h"



#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CSatCoordDoc

IMPLEMENT_DYNCREATE(CSatCoordDoc, CDocument)

BEGIN_MESSAGE_MAP(CSatCoordDoc, CDocument)
	//{{AFX_MSG_MAP(CSatCoordDoc)
	ON_COMMAND(ID_FILE_OPEN, OnFileOpen)
	ON_COMMAND(ID_COOR_CAL, OnCoorCal)
	ON_COMMAND(ID_COOR_NSAT, OnCoorNsat)
	ON_COMMAND(ID_COOR_CLO, OnCoorClo)
	ON_COMMAND(ID_COOR_GSAT, OnCoorGsat)
	ON_COMMAND(ID_WGS_BJ54, OnWgsBj54)
	ON_COMMAND(ID_PZ90_WGS, OnPz90Wgs)
	ON_COMMAND(ID_COOR_RSP3, OnCoorRsp3)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CSatCoordDoc construction/destruction



CSatCoordDoc::CSatCoordDoc()
{
	// TODO: add one-time construction code here
	Count	=	0;
}

CSatCoordDoc::~CSatCoordDoc()
{
}

BOOL CSatCoordDoc::OnNewDocument()
{
	if (!CDocument::OnNewDocument())
		return FALSE;

	// TODO: add reinitialization code here
	// (SDI documents will reuse this document)

	return TRUE;
}



/////////////////////////////////////////////////////////////////////////////
// CSatCoordDoc serialization

void CSatCoordDoc::Serialize(CArchive& ar)
{
	if (ar.IsStoring())
	{
		// TODO: add storing code here
	}
	else
	{
		// TODO: add loading code here
	}
}

/////////////////////////////////////////////////////////////////////////////
// CSatCoordDoc diagnostics

#ifdef _DEBUG
void CSatCoordDoc::AssertValid() const
{
	CDocument::AssertValid();
}

void CSatCoordDoc::Dump(CDumpContext& dc) const
{
	CDocument::Dump(dc);
}
#endif //_DEBUG

/////////////////////////////////////////////////////////////////////////////
// CSatCoordDoc commands

void CSatCoordDoc::OnFileOpen() 
{
	
}

void CSatCoordDoc::OnCoorCal() 
{
	// TODO: Add your command handler code here
	FILE* fp;

	// TODO: Add your command handler code here
	CString szFilter = "Data Files (*.??N)|*.??N|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if( dlg.DoModal() == IDOK )
	{
		fp =fopen(dlg.GetPathName(), "r");
		m_ron.Read_N_Header( fp );
		m_ron.Read_N_Data(fp,Data,ephemeris,m_ron.m_DIM);

	}
}

void CSatCoordDoc::OnCoorNsat() 
{
	// TODO: Add your command handler code here

	FILE* fp;
	FILE* fp_checkEpheinfo1;
	FILE* fp_svpos;
	fp_checkEpheinfo1=fopen("N_CheckEphinfo1.txt","w");
	fp_svpos=fopen("N_svpos.txt","w");

	// TODO: Add your command handler code here
	CString szFilter = "Data Files (*.??N)|*.??N|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if( dlg.DoModal() == IDOK )
	{
		fp =fopen(dlg.GetPathName(), "r");
		m_ron.Read_N_Header( fp );
		m_ron.Read_N_Data(fp,Data,ephemeris,m_ron.m_DIM);
	}
	
			
		svposx = new double[m_ron.m_DIM];
	
		svposy = new double[m_ron.m_DIM];

		svposz = new double[m_ron.m_DIM];

		Count = m_ron.m_DIM;
	for(int i=0;i<m_ron.m_DIM;i++)
	{
		ephemeris.PRN  = int(Data[i*36+ 0]);//unsigned short  PRN/Satellite PRN number
		if(Data[i*36+ 1]>=50)         //int year;
		{
			ephemeris.year = 1900+int(Data[i*36+ 1]);
		}
		else if(int(Data[i*36+ 1])>=0
			         &&int(Data[i*36+ 1])<50)
		{
			ephemeris.year = 2000+int(Data[i*36+ 1]);
		}
		ephemeris.month = int(Data[i*36+ 2]);//int month;
		ephemeris.day   = int(Data[i*36+ 3]);//int day;
		ephemeris.hour  = int(Data[i*36+ 4]);//int hour;
		ephemeris.minute= int(Data[i*36+ 5]);//int minute;
		ephemeris.second= int(Data[i*36+ 6]);//double second;
		ephemeris.a0   = Data[i*36+ 7];//double af0;//SV_clock_bias;
		ephemeris.a1   = Data[i*36+ 8];//double af1;//SV_clock_drift;//(sec/sec)
		ephemeris.a2   = Data[i*36+ 9];//double af2;//SV_clock_drift_rate;//(sec/sec2)

		fprintf(fp_checkEpheinfo1,"%2d %d %2d %2d %2d %2d %f %20.11e %20.11e %20.11e\n",
					ephemeris.PRN,ephemeris.year,ephemeris.month,ephemeris.day,ephemeris.hour,
					ephemeris.minute,ephemeris.second,ephemeris.a0,ephemeris.a1,ephemeris.a2);

		//ORBIT 1:
		ephemeris.IODE    = Data[i*36+ 10];//double IODE;//Issue of Data, Ephemeris
		ephemeris.Crs     = Data[i*36+ 11];//double Crs;//(meters)
		ephemeris.Delta_n = Data[i*36+ 12];//double Delta_n;//(radians/sec)
	    ephemeris.M0      = Data[i*36+ 13];//double M0;//(radians)
		
		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.IODE,ephemeris.Crs,ephemeris.Delta_n,ephemeris.M0);
		
		//ORBIT 2:
		ephemeris.Cuc       = Data[i*36+ 14];//double Cuc;//(radians) 
		ephemeris.e         = Data[i*36+ 15];//double e;//Eccentricity
		ephemeris.Cus       = Data[i*36+ 16];//double Cus;//(radians)
		ephemeris.sqrt_of_A = Data[i*36+ 17];//double sqrt_of_A;//(sqrt(m))
		
		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.Cuc,ephemeris.e,ephemeris.Cus,ephemeris.sqrt_of_A);
	
		//ORBIT 3:
		ephemeris.Toe   = Data[i*36+ 18];//double Toe;//Time of Ephemeris,(sec of GPS week)
		ephemeris.Cic   = Data[i*36+ 19];//double Cic;//(radians)
		ephemeris.OMEGA0 = Data[i*36+ 20];//double OMEGA;//(radians)
		ephemeris.CIS   = Data[i*36+ 21];//double CIS;//(radians)

		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.Toe,ephemeris.Cic,ephemeris.OMEGA0,ephemeris.CIS);

		//ORBIT 4:
		ephemeris.i0        = Data[i*36+ 22];//double i0;//(radians)	
		ephemeris.Crc       = Data[i*36+ 23];//double Crc;//(meters)
		ephemeris.omega     = Data[i*36+ 24];//double omega;//(radians)
		ephemeris.OMEGA_DOT = Data[i*36+ 25];//double OMEGA_DOT;//(radians/sec)
	
		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.i0,ephemeris.Crc,ephemeris.omega,ephemeris.OMEGA_DOT);
	
		//ORBIT 5:
		ephemeris.I_DOT                 = Data[i*36+ 26];//double IDOT;//(radians/sec)
		ephemeris.Codes_on_L2			= Data[i*36+ 27];//int Codes_on_L2_channel;
		ephemeris.GPS_week				= Data[i*36+ 28];//int GPS_Week_go_with_TOE;
		ephemeris.L2_pdata_flag			= Data[i*36+ 29];//int L2_P_data_flag;

		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.I_DOT,ephemeris.Codes_on_L2,
					ephemeris.GPS_week,ephemeris.L2_pdata_flag);
	
		//ORBIT 6:
		ephemeris.SV_accuracy = Data[i*36+ 30];//double SV_accuracy;//(meters)
		ephemeris.SV_health   = Data[i*36+ 31];//int SV_health;
		ephemeris.TGD         = Data[i*36+ 32];//double TGD;//(seconds)
		ephemeris.IODC        = Data[i*36+ 33];//double IODC;//Issue of Data, Clock
	
		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e %20.11e %20.11e\n",
			        ephemeris.SV_accuracy,ephemeris.SV_health,ephemeris.TGD,ephemeris.IODC);

		//ORBIT 7:
		ephemeris.Transtime = Data[i*36+ 34];
	        //double Transmission_time_of_message;//(sec of GPS week, derived e.g.from Z-count in Hand Over Word (HOW)
		ephemeris.Fit_Interval = Data[i*36+ 35];
	
		fprintf(fp_checkEpheinfo1,"  \t%20.11e %20.11e\n",
			        ephemeris.Transtime,ephemeris.Fit_Interval);
	
		double obs_week, obs_gpstime;
		obs_week = m_ron.GetGPSTime( ephemeris.year,
			            ephemeris.month, ephemeris.day,ephemeris.hour,
						ephemeris.minute,ephemeris.second,&obs_gpstime);
		m_ron.Cal_Sate_Coord(ephemeris,obs_gpstime, svpos);

        svposx[i] = svpos[0];
		svposy[i] = svpos[1];
		svposz[i] = svpos[2];
		
		fprintf(fp_svpos,"%25f%25f%25f\n",svpos[0],svpos[1],svpos[2]);
		
	}

	fclose(fp_checkEpheinfo1);												  ///
	fclose(fp_svpos);
	//fp_resultFileData=fopen("N_SateCoordCheck.txt","w");/////////////////
	//将所有导航星历文件的数据春处在 N_DataResultCheck.txt 中////////////
	//kaiguan=0;													  ///
	//m_Read.SaveAllNaviData(fp_resultFileData,Data,m_Read.m_DIM,0);	  ///
	//fclose(fp_resultFileData);///////////////////////////////////////////

}

void CSatCoordDoc::OnCoorClo() 
{
	FILE* fp;

	// TODO: Add your command handler code here
	CString szFilter = "Data Files (*.??G)|*.??G|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if( dlg.DoModal() == IDOK )
	{
		fp =fopen(dlg.GetPathName(), "r");
		m_ron.Read_G_Header( fp );
		m_ron.Read_G_Data(fp,Data,gephemeris,m_ron.m_DIM);

	}
	
}

void CSatCoordDoc::OnCoorGsat() 
{
	// TODO: Add your command handler code here
	
}



void CSatCoordDoc::OnWgsBj54() 
{
	// TODO: Add your command handler code here
    FILE* fp_temp;
	FILE* fp_read;

	fp_temp=fopen("WGS_BJ54_Trans.txt","w");

	int	i=0,j;
	char	szBuf[256];
//	double	svpos[3];
	CString strFilename;
	CString szFilter = "Data Files (*svpos.txt)|*svpos.txt|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if (dlg.DoModal() == IDOK)
	{
		strFilename = dlg.GetPathName();
//		return;
	}
    COORDINATECARTESIAN cartesian1;
	COORDINATECARTESIAN cartesian2;
	fp_read=fopen(strFilename,"r");
	fgetpos(fp_read,&pos);
	
	while (!feof (fp_read)) 
	{
		fgets (szBuf, MAX_BUF, fp_read);
		if(feof (fp_read))
			break;
		i++;
	}
//	matrix	a(k,36);
	fsetpos(fp_read,&pos);

	svposx = new double[i];
	svposy = new double[i];
	svposz = new double[i];


	for(j = 0;j < i;j++)
	{	
		fgets(szBuf,128,fp_read);
		sscanf(szBuf, "%lf %lf %lf",&svpos[0],&svpos[1],&svpos[2]);
		cartesian2.dX = svpos[0];
		cartesian2.dY = svpos[1];
		cartesian2.dZ = svpos[2];
		m_trans.WGStoBJ54(cartesian2, &cartesian1);
		svposx[j] = cartesian1.dX;
		svposy[j] = cartesian1.dY;
		svposz[j] = cartesian1.dZ;
		fprintf(fp_temp,"%25lf%25lf%25lf\n",cartesian1.dX,cartesian1.dY,cartesian1.dZ);

	}

	fclose(fp_temp);
	fclose(fp_read);
}

///////////////////////////////////////////////////////////////////////////
//                       PZ90坐标转WGS84坐标                             //
//                                                                       //
///////////////////////////////////////////////////////////////////////////
void CSatCoordDoc::OnPz90Wgs() 
{
	// TODO: Add your command handler code here
	FILE* fp_temp;
	FILE* fp_read;
	fp_temp=fopen("PZ_WGS_Trans.txt","w");
	
	char	szBuf[128];
//	double	svpos[3];
	CString strFilename;
	CString szFilter = "Data Files (*svpos.txt)|*svpos.txt|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if (dlg.DoModal() != IDOK)
	{
		return;
	}
    COORDINATECARTESIAN cartesian1;
	COORDINATECARTESIAN cartesian2;
	strFilename = dlg.GetPathName();
	fp_read=fopen(strFilename,"r");
	
	while (!feof (fp_read)) 
	{
		fgets(szBuf,128,fp_read);
		sscanf(szBuf, "%lf %lf %lf",&svpos[0],&svpos[1],&svpos[2]);
		cartesian2.dX = svpos[0];
		cartesian2.dY = svpos[1];
		cartesian2.dZ = svpos[2];
		m_trans.PZ90toWGS(cartesian2, &cartesian1);
		fprintf(fp_temp,"%25lf%25lf%25lf\n",cartesian1.dX,cartesian1.dY,cartesian1.dZ);
    }
	fclose(fp_temp);
	fclose(fp_read);
	
}
///////////////////////////////////////////////////////////////////////////
//                       读取GPS精密星历文件(*.SP3)                      //
//                                                                       //
///////////////////////////////////////////////////////////////////////////
void CSatCoordDoc::OnCoorRsp3() 
{
	// TODO: Add your command handler code here
	FILE* fp;
//	double*	S_data;
	int i,j;
	CString szFilter = "Data Files (*.SP3)|*.SP3|All Files (*.*)|*.*||";
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT, szFilter);
	if( dlg.DoModal() != IDOK )
	{
		return;
	}
	fp =fopen(dlg.GetPathName(), "r");
	m_ron.m_DIM = m_ron.Get_SP3_Dimension( fp );
	S_data.SetSize(m_ron.m_DIM * 28 * 5 + m_ron.m_DIM * 6);
//	S_data = new double [m_ron.m_DIM * 28 * 5];

	fp =fopen(dlg.GetPathName(), "r");
	m_ron.Read_SP3_Data(fp,S_data,sephemeris,m_ron.m_DIM);

	pSP3TIME	=	new	SP3TIME[m_ron.m_DIM];
	svposx = new double[m_ron.m_DIM * 28];
	svposy = new double[m_ron.m_DIM * 28];
	svposz = new double[m_ron.m_DIM * 28];
	clockbias=new double[m_ron.m_DIM * 28];
	xuhao	= new int[m_ron.m_DIM * 28];

	for(i = 0;i < m_ron.m_DIM;i++)
	{
		pSP3TIME[i].year	= (int)S_data[ i*146];
		pSP3TIME[i].month	= (int)S_data[ i*146 + 1];
		pSP3TIME[i].day		= (int)S_data[ i*146 + 2];
		pSP3TIME[i].hour	= (int)S_data[ i*146 + 3];
		pSP3TIME[i].minute	= (int)S_data[ i*146 + 4];
		pSP3TIME[i].second	= S_data[ i*146 + 5];
		for(j = 0; j < 28;j++)
		{
			xuhao[i * 28 + j]		= (int)S_data[ i*146 + j * 5 + 6];
			svposx[ i*28 + j]		= S_data[ i*146 + j * 5 + 7];
			svposy[ i*28 + j]		= S_data[ i*146 + j * 5 + 8];
			svposz[ i*28 + j]		= S_data[ i*146 + j * 5 + 9];
			clockbias[ i*28 + j]	= S_data[ i*146 + j * 5 + 10];
		}
	}
	Count = m_ron.m_DIM;
}