1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
#include "calendar.h"
using namespace std;
// 计算地球日心黄经
double parameter::get_longitude(vector<double> l,double t)
{
double L=l[0];
for(int i=1;i<l.size();i++)
{
L+=l[i]*pow(t,i);
}
//返回弧度制
return L+pi;
}
// 计算地球日心黄纬
double parameter::get_latitude(vector<double> b,double t)
{
double B=b[0];
for(int i=1;i<b.size();i++)
{
B+=b[i]*pow(t,i);
}
//返回弧度制
//地心黄纬=-太阳黄纬
return -B;
}
double parameter::get_R(vector<double> r,double t)
{
double R=r[0];
for(int i=1;i<r.size();i++)
{
R+=r[i]*pow(t,i);
}
return R;
}
// 转换FK5误差，返回弧度制
double parameter::delta_FK5(double L,double B,double T)
{
//L转角度制
L*=180/pi;
//计算L'，角度制
double L_p=L-1.397*T-0.00031*T*T;
//计算L'，弧度制
L_p*=pi/180;
while(L_p<0)
L_p+=2*pi;
while(L_p>2*pi)
L_p-=2*pi;
//计算delta_L，单位为角秒
double delta_L=-0.09033+0.03916*(cos(L_p)+sin(L_p))*tan(B);
//转换为弧度制
delta_L*=pi/(180*3600);
return delta_L;
}
//获取章动有关角
vector<double> parameter::get_angle(double T)
{
vector<double> ang;
//计算平距角
ang.push_back(297.85036+445267.111480*T-0.0019142*T*T+T*T*T/189474);
//计算日地平近点角
ang.push_back(357.52772+35999.050340*T-0.0001603*T*T-T*T*T/300000);
//计算月球平近点角
ang.push_back(134.96298+477198.867398*T+0.0086972*T*T+T*T*T/56250);
//计算月球纬度参数
ang.push_back(93.27191+483202.017538*T-0.0036825*T*T+T*T*T/327270);
//计算黄道与月球平轨道交点黄经
ang.push_back(125.04452-1934.136261*T+0.0020708*T*T+T*T*T/450000);
return ang;
}
//章动修正
double parameter::nution(double T)
{
vector<double> ang=get_angle(T);
fstream fr;
fr.open("data/nutation.txt",ios::in);
string str;
double sita;
double s,s1,s2;
double c,c1,c2;
double m,n,o,p,q;
s=c=0;
while(getline(fr,str))
{
sscanf(str.c_str(),"%lf%lf%lf%lf%lf%lf%lf%lf%lf",&m,&n,&o,&p,&q,&s1,&s2,&c1,&c2);
sita=m*ang[0]+n*ang[1]+o*ang[2]+p*ang[3]+q*ang[4];
//sita为角度制，转换为弧度制
sita*=pi/180;
//计算黄经章动，单位为0.0001角秒
s+=(s1+s2*T)*sin(sita);
//计算交角章动，单位为0.0001角秒
//c+=(c1+c2*T)*cos(sita);
}
//关闭文件
fr.close();
// 计算得到的s和c单位0.0001秒，转换为弧度制
s*=0.0001;
s*=pi/(180*3600);
return s;
}
//光行差修正，返回弧度制
double parameter::aberration(double R)
{
//return -20.4898/R/3600.0*pi/180.0;
//公式，单位0.0001角秒
double delta=-20.4898/R;
//转换为弧度制
delta*=pi/(180*3600);
return delta;
}
// 获取地日运行参数，L为地球日心黄经，B为地球日心黄纬，R为地日距离
vector<vector<double>> parameter::get_parameters(double t)
{
fstream fr;
fr.open("data/earth.txt",ios::in);
string s;
double l=0;
double a,b,c,num_tmp=0;
vector<vector<double>> parameters;
vector<double> tmp;
while(getline(fr,s))
{
//当分割符为/时，表明一组参数结束，将参数存入parameters中，并清空tmp
if(s[0]=='/')
{
parameters.push_back(tmp);
tmp.resize(0);
num_tmp=0;
}
//当分割符为=时，表明一组参数中的一类结束，将参数存入tmp中，并清空num_tmp
else if(s[0]=='=')
{
tmp.push_back(num_tmp);
num_tmp=0;
}
else
{
sscanf(s.c_str(),"%lf%lf%lf",&a,&b,&c);
num_tmp+=a*cos(b+c*t);
}
}
//关闭文件
fr.close();
/* 最终传参情况：
* parameters[0]为L系列，用来计算地球日心黄经
* parameters[1]为B系列，用来计算地球日心黄纬
* parameters[2]为R系列，用来计算地日距离
*/
return parameters;
}
double parameter::sun_longitude(double t)
{
//从文件中计算地日运行参数
vector<vector<double>> p=get_parameters(t);
//取出各类参数
vector<double> l=p[0],b=p[1],r=p[2];
//计算地球日心黄经纬，返回弧度制/天文单位
double L=get_longitude(l,t);
double B=get_latitude(b,t);
double R=get_R(r,t);
//以下修正需要的都是儒略世纪数而非千年数
//转FK5误差，返回弧度制
L+=delta_FK5(L,B,10*t);
//章动修正，返回弧度制
L+=nution(10*t);
//光行差修正，返回弧度制
L+=aberration(R);
//L转角度制
L*=180/pi;
while(L<0)
L+=360;
while(L>=360)
L-=360;
if(360-L<=5)
{
L-=360;
}
return L;
}
//double moon_longitude(double t)
//{
// //文件中获取地月运行参数
// vector<vector<double>> p=get_moon_parameters(t);
// vector<double> l=p[0],b=p[1],r=p[2];
// //计算地月日心黄经纬，返回弧度制/天文单位
// double L=get_moon_longitude(l,t);
// double B=get_moon_latitude(b,t);
// double R=get_moon_R(r,t);
//}
|
