📄 correctmodel.cpp
字号:
// CorrectModel.cpp: implementation of the CCorrectModel class.
//
//////////////////////////////////////////////////////////////////////
#include "stdafx.h"
#include "compute.h"
#include "CorrectModel.h"
#include "fstream.h"
#include "stdio.h"
#include "ReadJPL_EPH.h"
#include "iomanip.h" //Just for set precision
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CCorrectModel::CCorrectModel()
{
PI=3.14159265;
ET=2452277.5;
DayRest=0.2;
}
CCorrectModel::~CCorrectModel()
{
}
double Factorial(int N)
{
if(N<=1) return 1;
else
return Factorial(N-1)*N;
}
// Coordinate transform
void XYZ_FeLamRe(double X,double Y,double Z,double *Fe,double *Lam,double *Re)
{
double Rxy;
double PI=3.1415926;//add
*Re=sqrt(X*X+Y*Y+Z*Z);
Rxy=sqrt(X*X+Y*Y);
if (Rxy!= 0)
*Fe=atan(Z/Rxy);
else if(Z > 0)
*Fe=PI/2.0;
else
*Fe=-PI/2.0;
if(X!=0)
{
*Lam=atan(Y/X);
if (X < 0)
*Lam=*Lam+PI;
if (Y < 0)
*Lam=*Lam+2.0*PI;
}
else
{
if (Y > 0)
*Lam=PI/2.0;
else
*Lam=-PI/2.0;
}
}
void JUL2YMD_M(double Jul,int &Year,int &Month,int &Day,int &Hour,int &Minute,double &Second)
{
int a,b,c,d,e;
double H;
a = int(Jul + 0.5 + 1.0e-16);
b = a + 1537;
c = int((b - 122.1)/365.25 + 1.0e-16);
d = int(365.25*c + 1.0e-16);
e = int((b-d)/30.6001 + 1.0e-16);
Month = e - 1 - 12 * int(e/14.0 + 1.0e-16);
Day = b -d - int(30.6001 * e);
Year = c - 4715 -int((7+Month)/10.0 + 1.0e-16);
H = (Jul + 0.5 - int(Jul + 0.5 +1.0e-16))*24.0;
Hour = int(H + 1.0e-16);
Minute = int((H - Hour)*60.0+1.0e-16);
Second = ((H - Hour)*60.0 - Minute)*60.0;
}
void CCorrectModel::GetOceanTide(double DeltCoef[])
{
double t,Jul;
double F1,F2,F3,F4,F5; // F1=l,F2=l',F3=F;F4=D,F5=Omga;
double D1,D2,D3,D4,D5,D6;
double Tu,Theta_g;
Jul=ET; //add I think it is wrong,I will correct it in the future!
//Jul=Ymd2JDT();
//set precision
setiosflags( ios::left );
setprecision(16);
t = (ET - 2451545.0)/36525.0;
F1 = 134.96340251+(1717915923.2178*t+31.8792*t*t+0.051635*t*t*t-0.00024470*t*t*t*t)/3600.0;
F2 = 357.52910918+129596581.0481/3600.0*t+(-0.5532*t*t+0.000136*t*t*t-0.00001149*t*t*t*t)/3600.0;
F3 = 93.27209062+1739527262.8478/3600.0*t+(-12.7512*t*t-0.001037*t*t*t+0.00000417*t*t*t*t)/3600.0;
F4 = 297.85019547+1602961601.2090/3600.0*t+(-6.3706*t*t+0.006593*t*t*t-0.00003169*t*t*t*t)/3600.0;
F5 = 125.04455501-6962890.2665/3600.0*t+(7.4722*t*t+0.007702*t*t*t-0.00005939*t*t*t*t) /3600.0;
Tu = int(Jul - 2451545.0)/36525.0; //% Tu
Theta_g= (6+(41+50.54841/60.0)/60.0+(8640184.812866*Tu+0.093104*Tu*Tu-6.210e-6*Tu*Tu*Tu)/3600)*PI/12; //% in rad
// F1=l,F2=l',F3=F;F4=D,F5=Omga;
F1=F1*PI/180.0; F2=F2*PI/180.0; // in rad
F3=F3*PI/180.0; F4=F4*PI/180.0; // in rad
F5=F5*PI/180.0; // in rad
D2=F3+F5;
D3=D2-F4;
D4=D2-F1;
D5=-F5;
D6=D2-F4-F2;
D1=Theta_g+PI-D2;
/////////////////////////////////////////////////////////////////////////////
// read ocean tide file
CStdioFile OTideFile;
CString strTemp;
char Path[MAX_PATH];
CString FullPath;
GetCurrentDirectory(MAX_PATH,Path);
FullPath=Path;
FullPath+="\\Model\\OTIDES";
if(!OTideFile.Open(FullPath,CFile::modeRead))
{
AfxMessageBox("无法打开文件:"+FullPath);
return ;
}
// read the header of ocean tide file
int i;
int N1,N2,N3,N4,N5,N6;
int M_Max,N_Max;
double C_Add,S_Add,C_Minus,S_Minus; // the coeffinents for ocean model
for(i=0;i<7;i++)
{
OTideFile.ReadString(strTemp);
}
OTideType Tide;
////////////////////////////////////////////////////
while(OTideFile.ReadString(strTemp))
{
if(strTemp.IsEmpty()) break;
N1=atoi(strTemp.Mid(13,1));
N2=atoi(strTemp.Mid(14,1));
N3=atoi(strTemp.Mid(15,1));
N4=atoi(strTemp.Mid(17,1));
N5=atoi(strTemp.Mid(18,1));
N6=atoi(strTemp.Mid(19,1));
N_Max=atoi(strTemp.Mid(24,2));
M_Max=atoi(strTemp.Mid(26,2));
C_Add=atof(strTemp.Mid(31,21));
C_Minus=atof(strTemp.Mid(53,21));
S_Add=atof(strTemp.Mid(75,21));
S_Minus=atof(strTemp.Mid(97,21));
Tide.M=M_Max; Tide.N=N_Max;
Tide.C1=C_Add; Tide.C2=C_Minus;
Tide.S1=S_Add; Tide.S2=S_Minus;
Tide.Theths=N1*D1+(N2-5)*D2+(N3-5)*D3+(N4-5)*D4+(N5-5)*D5+(N6-5)*D6;
m_OTide.Add(Tide);
}
OTideFile.Close();
//////////////////////////////////////////////////////////////////////////////////
// 计算海潮摄动
int Num=m_OTide.GetUpperBound()+1;
double C20,C21,C22,S21,S22;
double F20,F21,F22;
C20=C21=C22=S21=S22=0.0;
F20=GetFnm(2,0)/100.0;
F21=GetFnm(2,1)/100.0;
F22=GetFnm(2,2)/100.0;
for(i=0;i<Num;i++)
{
Tide=m_OTide.GetAt(i);
if(Tide.N==2 && Tide.M==0)
C20+=(Tide.C1+Tide.C2)*cos(Tide.Theths)+(Tide.S1+Tide.S2)*sin(Tide.Theths);
if(Tide.N==2 && Tide.M==1)
{
C21+=(Tide.C1+Tide.C2)*cos(Tide.Theths)+(Tide.S1+Tide.S2)*sin(Tide.Theths);
S21+=(Tide.S1-Tide.S2)*cos(Tide.Theths)-(Tide.C1-Tide.C2)*sin(Tide.Theths);
}
if(Tide.N==2 && Tide.M==2)
{
C22+=(Tide.C1+Tide.C2)*cos(Tide.Theths)+(Tide.S1+Tide.S2)*sin(Tide.Theths);
S22+=(Tide.S1-Tide.S2)*cos(Tide.Theths)-(Tide.C1-Tide.C2)*sin(Tide.Theths);
}
}
C20=C20*F20; C21=C21*F21; C22=C22*F22;
S21=F21*S21; S22=S22*F22;
DeltCoef[0]=C20; DeltCoef[1]=C21; DeltCoef[2]=C22;
DeltCoef[3]=S21; DeltCoef[4]=S22;
}
double CCorrectModel::GetKn(int N)
{
double Kn;
switch(N)
{
case 2:
Kn=-0.3075;
break;
case 3:
Kn=-0.195;
break;
case 4:
Kn=-0.132;
break;
case 5:
Kn=-0.1032;
break;
case 6:
Kn=-0.0892;
break;
default:
break;
}
return Kn;
}
double CCorrectModel::GetFnm(int N, int M)
{
double Fnm,Kn,Tem,Tem1,Tem2;
double g0,Pw,G;
g0=9.798261; // m/s^2
Pw=1025.0; // kg/m^3
G=6.673e-11; // m^3/kg.s^2
/////////////////////////////////////////////////////////////
Tem=4.0*PI*G*Pw/g0;
Tem1=Factorial(N+M);
Tem2=Factorial(N-M);
Kn=GetKn(N);
if(M==0)
Fnm=Tem*sqrt(Tem1/Tem2/(2*N+1))*((1+Kn)/(2*N+1));
else
Fnm=Tem*sqrt(Tem1/Tem2/(2*N+1)/2.0)*((1+Kn)/(2*N+1));
return Fnm;
}
void CCorrectModel::GetSoildTide(double DeltCoef[])
{
// DeltCoef=[C20,C21,C22,S21,S22]
double GM,GMm,GMs;
double k2,Re;
double DC20,DC21,DC22,DS21,DS22;
double Fe1,Lam1,Fe2,Lam2,R1,R2;
double Theta1,Theta2;
double r1,r2;
double Tem1,Tem2;
double Pm20,Pm21,Pm22; // legern function about mooon
double Ps20,Ps21,Ps22; // legern function about sun
double A_GM; // A_GM stand for A^3/GM
// r1 is the distance between moon and earth
// r2 is the distance between sun and earth
// initial the GM GMm and GMs
GM=3.986005e+14; // the unit is m^3/s^2
GMm=4902.793; // the unit is km^3/s^2
GMs=132712438000.0; // the unit is km^3/s^2
Re=6378137.0; // the unit is m
k2=0.3;
XYZ_FeLamRe(m_Xm,m_Ym,m_Zm,&Fe1,&Lam1,&r1); // Transfor the X,Y,Z to Fe1,Lam1,R1 from the position of moon
XYZ_FeLamRe(m_Xs,m_Ys,m_Zs,&Fe2,&Lam2,&r2); // Transfor the X,Y,Z to Fe2,Lam2,R2 from the position of sun
Theta1=PI/2-Fe1;
Theta2=PI/2-Fe2;
//r1=sqrt(m_Xm*m_Xm+m_Ym*m_Ym+m_Zm*m_Zm);
//r2=sqrt(m_Xs*m_Xs+m_Ys*m_Ys+m_Zs*m_Zs);
A_GM=Re*Re*Re/GM;
Pm20=sqrt(5)*(1.5*cos(Theta1)*cos(Theta1)-0.5);
Pm21=sqrt(15)*cos(Theta1)*sin(Theta1);
Pm22=sqrt(15)*sin(Theta1)*sin(Theta1)/2.0;
Ps20=sqrt(5)*(1.5*cos(Theta2)*cos(Theta2)-0.5);
Ps21=sqrt(15)*cos(Theta2)*sin(Theta2);
Ps22=sqrt(15)*sin(Theta2)*sin(Theta2)/2.0;
Tem1=GMm*Pm20/r1/r1/r1;
Tem2=GMs*Ps20/r2/r2/r2;
DC20=k2*A_GM*(Tem1+Tem2)/sqrt(5.0);
Tem1=GMm*Pm21/r1/r1/r1;
Tem1=Tem1*cos(Lam1);
Tem2=GMs*Ps21/r2/r2/r2;
Tem2=Tem2*cos(Lam2);
DC21=k2*A_GM*(Tem1+Tem2)/sqrt(15.0);
Tem1=GMm*Pm21/r1/r1/r1;
Tem1=Tem1*sin(Lam1);
Tem2=GMs*Ps21/r2/r2/r2;
Tem2=Tem2*sin(Lam2);
DS21=k2*A_GM*(Tem1+Tem2)/sqrt(15.0);
Tem1=GMm*Pm22/r1/r1/r1;
Tem1=Tem1*cos(2*Lam1);
Tem2=GMs*Ps22/r2/r2/r2;
Tem2=Tem2*cos(2*Lam2);
DC22=k2*A_GM*(Tem1+Tem2)/sqrt(60.0);
Tem1=GMm*Pm22/r1/r1/r1;
Tem1=Tem1*sin(2*Lam1);
Tem2=GMs*Ps22/r2/r2/r2;
Tem2=Tem2*sin(2*Lam2);
DS22=k2*A_GM*(Tem1+Tem2)/sqrt(60.0);
DeltCoef[0]=DC20; DeltCoef[1]=DC21;DeltCoef[2]=DC22;
DeltCoef[3]=DS21;DeltCoef[4]=DS22;
}
void CCorrectModel::GetPoleTide(double DeltCoef[])
{
ReadPole(); // 读入极移参数
double DC21,DS21;
DC21=-1.348e-9*(Xp+0.0112*Yp);
DS21=1.348e-9*(Yp-0.0112*Xp);
DeltCoef[0]=DC21;
DeltCoef[1]=DS21;
}
void CCorrectModel::ReadPole()
{
BOOL m_bFind;
int month1,day1;
int month2,day2;
//double xPol,yPol,UT1_UTC,UT1_UT1R,D,ddPsi,ddEps;
double xPol1,yPol1,UT1_UTC1,UT1_UT1R1,D1,ddPsi1,ddEps1;
double xPol2,yPol2,UT1_UTC2,UT1_UT1R2,D2,ddPsi2,ddEps2;
double MJD0,MJD1,MJD2;
MJD0 = ET - 2400000.5;
m_bFind=FALSE;
CStdioFile PoleFile;
CString strTemp;
// CString Name=_T("Bulletin.txt");
char Path[MAX_PATH];
CString FullPath;
GetCurrentDirectory(MAX_PATH,Path);
FullPath=Path;
FullPath+="\\Model\\Bulletin.txt";
if(!PoleFile.Open(FullPath,CFile::modeRead))
{
AfxMessageBox("无法打开文件:"+FullPath);
return ;
}
while (PoleFile.ReadString(strTemp))
{
if(strTemp.IsEmpty()) break;
int k2=sscanf(strTemp,"%d%d%lf%lf%lf%lf%lf%lf%lf%lf",
&month2,&day2,&MJD2,&xPol2,&yPol2,&UT1_UTC2,&UT1_UT1R2,&D2,&ddPsi2,&ddEps2);
if(k2!=10) {AfxMessageBox("数据格式不对!");break;}
if(MJD2<MJD0)
{
month1=month2; day1=day2; MJD1=MJD2; xPol1=xPol2;yPol1=yPol2;
UT1_UTC1=UT1_UTC2; UT1_UT1R1=UT1_UT1R2; D1=D2;ddPsi1=ddPsi2;
ddEps1=ddEps2;
}
else
{
m_bFind=TRUE;
break;
}
}
if(m_bFind)
{
Xp = ((xPol2- xPol1)*DayRest +xPol1); // in second
Yp = ((yPol2 - yPol1)*DayRest + yPol1); // in second
}
else
{
AfxMessageBox("没有发现对应的时间!");return;
}
PoleFile.Close();
}
void CCorrectModel::SetJUL(double JD)
{
ET=JD;
JUL2YMD_M(ET,m_year,m_month,m_day,m_hour,m_minute,m_second);
DayRest=(m_hour+m_minute/60.0+m_second/3600.0)/24.0;
CalSunPos(); // Get positon of Sun from JPLEPH
CalMoonPos();// Get positon of Moon from JPLEPH
}
void CCorrectModel::CalMoonPos()
{
int ntarg,ncent;
double RRD[6];
m_ReadJPL.InitJPL_EPH();
ntarg=10;
ncent=3;
m_ReadJPL.pleph(ET,ntarg,ncent,RRD);
m_Xm=RRD[0]; m_Ym=RRD[1]; m_Zm=RRD[2];
m_ReadJPL.CloseJPL_EPH();
}
void CCorrectModel::CalSunPos()
{
int ntarg,ncent;
double RRD[6];
m_ReadJPL.InitJPL_EPH();
ntarg=11;
ncent=3;
m_ReadJPL.pleph(ET,ntarg,ncent,RRD);
m_Xs=RRD[0]; m_Ys=RRD[1]; m_Zs=RRD[2];
m_ReadJPL.CloseJPL_EPH();
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -