// satpos.cpp : Defines the entry point for the console application.
//
// PI = 4*atan(1.0)
// rou = (1*180/PI)*60*60
#include "stdafx.h"
#include <stdio.h>
#include "math.h"
#define PI 3.1415926
#include "Matrix.h"
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#ifndef _NO_NAMESPACE
using namespace std;
using namespace math;
#define STD std
#else
#define STD
#endif
#ifndef _NO_TEMPLATE
typedef matrix<double> Matrix;
#else
typedef matrix Matrix;
#endif
#define DELTETEPRN15
#ifndef _NO_EXCEPTION
# define TRYBEGIN() try {
# define CATCHERROR() } catch (const STD::exception& e) { \
cerr << "Error: " << e.what() << endl; }
#else
# define TRYBEGIN()
# define CATCHERROR()
#endif
struct ORBITPARAMETER
{
double da; //the semimajor axis
double de; //the eccentricity
double domega; //the right ascension of the ascending node
double dw; //the argument of perigee
double di; //the inclination
double dto;
};
typedef struct tagATIME
{
unsigned int wYear;
unsigned char byMonth;
unsigned char byDay;
unsigned char byHour;
unsigned char byMinute;
double dSecond;
unsigned char byDayOfWeek;
} ATIME;
typedef ATIME* PTIME;
typedef struct tagJULIANDAY
{
long lDay;
long lSecond;
double dFraction;
} JULIANDAY;
typedef JULIANDAY* PJULIANDAY;
typedef struct tagGPSTIME
{
long lWeek;
long lSecond;
double dFraction;
} GPSTIME;
typedef GPSTIME* PGPSTIME;
//年月日－儒略日
void tmTimeToJulianDay (PTIME pTime, PJULIANDAY pJulianDay)
{
int m;
int y;
double ut;
ut = pTime->byHour + pTime->byMinute/60.0 + pTime->dSecond/3600.0;
if (pTime->byMonth <= 2) {
y = pTime->wYear - 1;
m = pTime->byMonth + 12;
}
else {
y = pTime->wYear;
m = pTime->byMonth;
}
pJulianDay->lDay = (long)(365.25*y) + (long)(30.6001*(m+1)) + pTime->byDay
+ (long)(ut/24 + 1720981.5);
pJulianDay->lSecond = ((pTime->byHour+12)%24)*3600L + pTime->byMinute*60L
+ (long)pTime->dSecond;
pJulianDay->dFraction = pTime->dSecond - (long)pTime->dSecond;
}
//年月日－GPS时间
void tmTimeToGPSTime (PTIME pTime, PGPSTIME pGPSTime)
{
JULIANDAY jd;
double dTemp;
long lTemp;
tmTimeToJulianDay (pTime, &jd);
dTemp = jd.lDay + (double)(jd.lSecond + jd.dFraction) / 86400.0;
lTemp = (long)(dTemp-2444244.5);
pGPSTime->lWeek = (long)((lTemp)/7);
pGPSTime->lSecond = (long)(lTemp%7*86400L) + pTime->byHour * 3600L
+ pTime->byMinute * 60L + (long)pTime->dSecond;
pGPSTime->dFraction = pTime->dSecond - (long)pTime->dSecond;
}
/*角度换弧度*/
double degree_rad ( float n )
{
double degree;
double minute;
double second;
double rad;
degree = (int)(n);
minute = (int)((n-degree)*100);
second = (n-degree-minute/100.0)*10000;
rad = (degree+minute/60.0+second/3600.0)*PI/180.0;
return(rad);
}
int main(int argc, char* argv[])
{
FILE* fresult;
fresult = fopen("result.txt","wt");
struct ORBITPARAMETER orbit;
double dti;
double dM; //平近点角
double du; //地球引力常数
double dn; //平均角速度
double dtemp;
double dE; //偏近点角
double dEold;
double dr; //卫星到地心的距离
double dsita;//真近点角
double dSG; //格林尼治真恒星时
ATIME timeti;//原子时
GPSTIME GPSTimeti;//GPS时间
ATIME timeto;
GPSTIME GPSTimet0;
timeto.wYear = 2005;
timeto.byMonth = 5;
timeto.byDay = 1;
timeto.byHour = 8;
timeto.byMinute = 0;
timeto.dSecond = 0;
timeti.wYear = 2005;
timeti.byMonth = 5;
timeti.byDay = 1;
timeti.byHour = 8;
timeti.byMinute = 30;
timeti.dSecond = 0;
tmTimeToGPSTime (&timeto,&GPSTimet0);
tmTimeToGPSTime (&timeti,&GPSTimeti);
orbit.dto = GPSTimet0.lSecond;
dti = GPSTimeti.lSecond;
orbit.da = 2e7;
orbit.de = 0.003;
orbit.domega = degree_rad(90.00);
orbit.dw = degree_rad(60.00);
orbit.di = degree_rad(55.00);
dSG = degree_rad(30.00);
du = sqrt(3.986005e14); //
dtemp = sqrt(orbit.da);
dn = du/(pow(dtemp,3));
dM = dn*(dti-orbit.dto);
//计算平近点角M
dE = dM;
do
{
dEold = dE;
dE = dM+orbit.de*sin(dEold);
} while(fabs(dEold-dE)>0.00001);
//计算卫星的轨道半径
dr = orbit.da*(1-orbit.de*cos(dE));
//计算真近点角sita
dsita = atan((sqrt((1+orbit.de)/(1-orbit.de))*tan(dE/2)))*2;
Matrix R1(3,3);
R1.Null(); //清零
R1.Unit(); //单位化
Matrix R2(3,3);
R2.Null();
R2.Unit();
Matrix R3(3,3);
R3.Null();
R3.Unit();
Matrix R(3,3);
R.Null();
Matrix XYZ(3,1);
XYZ.Null();
Matrix xyz(3,1);
xyz.Null();
//轨道坐标系到空间直角坐标系转换的旋转矩阵R
R1(0,0) = cos(orbit.domega-dSG);
R1(0,1) = -sin(orbit.domega-dSG);
R1(1,0) = sin(orbit.domega-dSG);
R1(1,1) = cos(orbit.domega-dSG);
R2(1,1) = cos(orbit.di);
R2(1,2) = -sin(orbit.di);
R2(2,1) = sin(orbit.di);
R2(2,2) = cos(orbit.di);
R3(0,0) = cos(orbit.dw);
R3(0,1) = -sin(orbit.dw);
R3(1,0) = sin(orbit.dw);
R3(1,1) = cos(orbit.dw);
xyz(0,0) = dr*cos(dsita);
xyz(1,0) = dr*sin(dsita);
R = R1*R2*R3;
XYZ = R*xyz;
fprintf(fresult,"平均角速度:\n%lf\n",dn);
fprintf(fresult,"平近点角:\n%lf\n",dM);
fprintf(fresult,"偏近点角E:\n%lf\n",dE);
fprintf(fresult,"卫星到地心距离r:\n%lf\n",dr);
fprintf(fresult,"真近点角sita:\n%lf\n",dsita);
fprintf(fresult,"卫星在轨道直角坐标系中的位置为:\n%lf\n%lf\n%lf\n",xyz(0,0),xyz(1,0),xyz(2,0));
fprintf(fresult,"卫星在天球坐标系中的位置为:\n%lf\n%lf\n%lf\n",XYZ(0,0),XYZ(1,0),XYZ(2,0));
fclose(fresult);
return 0;
}
