在wgs84 坐标系转到J2000 坐标系 主要 涉及到坐标的相互转换。一般给定的WGS坐标为 给定时刻的 B、L、H。 具体转换如下
step1 WGS 84 转换到协议地球坐标系。
function XYZ_m = wgs842ECEF(BLH_deg_m)
%大地经纬度高程BLH与地固坐标系的转换 BLH-->XG
%ECEF, Earth Centered Earth Fixed;地固坐标系 参考平面：平赤道面，即过地心，并且与地心与CIO点连线垂直的平面。
% x轴为参考平面与格林尼治平面交线，z轴为地心指向CIO点。
%参数
%B L H分别为大地纬度、大地经度和大地高程 输入参数为N*3矩阵
% XYZ 为地固坐标系的 x y z，输出参数为N*3矩阵
%ae = 6378137; %2000坐标系场半轴
%ee = 0.0818191910428; %第一偏心率
ae = 6378137; %wgs84坐标系场半轴
ee = 0.00669437999013 ; %第一偏心率
B_deg=BLH_deg_m(:,1);L_deg=BLH_deg_m(:,2);H_m=BLH_deg_m(:,3);
N_m=ae*sqrt(1-ee^2*sind(B_deg).^2);
XYZ_m(:,1)=(N_m+H_m).*cosd(B_deg).*cosd(L_deg);
XYZ_m(:,2)=(N_m+H_m).*cosd(B_deg).*sind(L_deg);
XYZ_m(:,3)=(N_m*(1-ee^2)+H_m).*sind(B_deg);
end
step 2 协议地球坐标系 转换为瞬时地球坐标系
这主要涉及查询 给定时刻的 xp,yp , 可以查询EOP文件获得，下载地址如下http://celestrak.com/spacedata/
earthFixedXYZ=ordinateSingleRotate(‘x’,yp)*ordinateSingleRotate(‘y’,xp)*earthFixedXYZ;
其中 ordinateSingleRotate 为 坐标转换 函数 ，该函数第一个参数为给定的坐标轴，第二个参数为 旋转的角度 ，第三个参数为 旋转的角度的度量单位。 后续要多次调用该函数。
function [R] = ordinateSingleRotate( axis,angle_deg,angleType)
%坐标轴旋转
% axis 表示围绕旋转的坐标轴
% '1'或者'x' 表示围绕 x轴逆时针旋转
% '2'或者'y' 表示围绕 y轴逆时针旋转
% '3'或者'z' 表示围绕 z轴逆时针旋转
% angle_deg 表示旋转的角度
% angleType 表示角度类型，可以为 'rad'和’deg‘，默认deg
if nargin()<3
angleType='deg';
end
switch lower(angleType)
case 'deg'
case 'rad'
angle_deg=angle_deg*180/pi;
otherwise
error(message('ordinateRotate:unknownRotation:angleType=', angleType));
end
switch axis
case {'x','X',1,'1'}
R =[1 0 0 ;...
0 cosd(angle_deg) sind(angle_deg);...
0 -sind(angle_deg) cosd(angle_deg) ];
case {'y','Y',2,'2'}
R =[cosd(angle_deg) 0 -sind(angle_deg);...
0 1 0;...
sind(angle_deg) 0 cosd(angle_deg) ];
case {'z','Z',3,'3'}
R =[cosd(angle_deg) sind(angle_deg) 0 ;...
-sind(angle_deg) cosd(angle_deg) 0;...
0 0 1 ];
otherwise
error(message('ordinateRotate:unknownRotation:axis=',axis));
end
end
step 3 瞬时地球坐标系 转换为 瞬时真天球坐标系
xyz= ordinateSingleRotate(‘z’,-gst_deg)*earthFixedXYZ;
该步骤主要涉及到格林尼治恒星时角的计算，关于格林尼治恒星时角 计算方法 很多，下面给出 一个较为精确的计算方法。其中dut1 和dat参数在EOP文件中有。EOP文件下载地址如下CelesTrak: EOP and Space Weather Data
function [gst_deg,JDTDB ] = utc2gst(UTC,dUT1,dAT)
%将utc时间转换为格林尼治恒星时
%参数
%utc 格式为y m d 其中d的数值为小数，将h m s 按（sec/3600+min/60+h)/24转换成小数，并累加到day 上
%dUT1 为ut1-utc 的差，数值不超过正负1秒，查iers可获得数值
%UNTITLED 计算当地恒星时,返回值以时秒为单位
% UTC世界协调时 协调地球时与原子时的时间差，采用闰秒的办法进行协调
%dAT 润秒数
%TAI国际原子时间 该时间最准 TAI=UTC+dAT; %国际原子时间
%TT 地球时 TT = TAI+32.184 至2014年的时候;
%TDT 地球动力学时间
% ET 历书时间
%地球时=地球动力学时间=历书时间
J2000=2451545;
if nargin()<2
dUT1=0;
end
if nargin()<3
dAT=37.0;
end
JDutc=YMD2JD(UTC(1),UTC(2),UTC(3));
JDUT1 = JDutc+dUT1/86400 ; %UT1 为世界时，世界时由于自转的不均匀，因此与UTC时间有dUT1的差别，dUT1在各国卫星授时信号中会以0.1秒的精度给出IERS经过处理后，会以1e-5的精度给出。
dT=32.184+dAT-dUT1;
JDTT=JDUT1+dT/86400;
% JDTT= YMD2JD(UTC(1),UTC(2),UTC(3))+dUT1/86400;
%首先计算TDB,TDB是质心动力学时，是太阳月球行星等天体星历表中的时间尺度
T=(JDTT-J2000)/36525;
temp= 0.001657*sin(628.3076*T+6.2401) ...
+0.000022*sin(575.3385*T+4.2970) ...
+0.000014*sin(1256.6152*T+6.1969) ...
+0.000005*sin(606.9777*T+4.0212)...
+0.000005*sin(52.9691*T+0.4444)...
+0.000002*sin(21.3299*T+5.5431)...
+0.00001*T*sin(628.3076*T+4.2490);
JDTDB=JDTT+temp/86400;
%下面计算UT2,UT2是在UT1的基础上修正地球周期性季节变化后得到的世界时间
% %根据UT1计算Tb,Tb为以贝塞尔年为单位，从历元B2000.0起算ff
% B2000=2451544.033;
% Tb=(YMD2JD(UT1(1),UT1(2),UT1(3))-B2000)/365.2422;
% dTs=0.022*sin(2*pi*Tb)-0.012*cos(2*pi*Tb)-0.006*sin(4*pi*Tb)+0.007*cos(4*pi*Tb);
% dTs=dTs/86400;
% UT2=UT1+[0 0 dTs];
%下面计算格林尼治平恒星时GMST; 单位为度
T=(JDTDB-J2000)/36525;
T2=T*T;T3=T2*T;T4=T3*T;T5=T4*T;
Du=JDUT1-J2000;
thita = 0.7790572732640+1.00273781191135448*Du; %自J2000起地球转过的圈数；
GMST_deg=(-0.0000000368*T5-0.000029956*T4-0.00000044*T3+1.3915817*T2+4612.156534*T+0.014506) /3600+(thita-floor(thita))*360 ;
[epthilongA_deg,dertaPthi_deg] = nutationInLongitudeCaculate(JDTDB);
%计算赤经章动引起的误差eclipticObliquitygama，单位为角秒
%计算月球平近点角 太阳平近点角 月球纬度辐角 日月平角距(月球平黄经-太阳平黄经） 月球升交点平黄经 单位为角秒
F_deg =(0.00000417*T4-0.001037*T3-12.7512*T2+1739527262.8478*T+335779.526232)/3600;
F_deg=mod(F_deg,360);
D_deg = (-0.00003169*T4+0.006593*T3-6.3706*T2+1602961601.2090*T+1072260.70369)/3600 ;
D_deg=mod(D_deg,360);
Omiga_deg =(-0.00005939*T4+0.007702*T3+7.4722*T2-6962890.5431*T+450160.398036 )/3600 ;
Omiga_deg=mod(Omiga_deg,360);
epthilongGama_deg=dertaPthi_deg*cosd(epthilongA_deg)...
+(0.00264096*sind( Omiga_deg )...
+0.00006352*sind(2*Omiga_deg)...
+0.00001175*sind(2*F_deg-2*D_deg+3*Omiga_deg)...
+0.00001121*sind(2*F_deg-2*D_deg+Omiga_deg)...
-0.00000455*sind(2*F_deg-2*D_deg+2*Omiga_deg)...
+0.00000202*sind(2*F_deg+3*Omiga_deg)...
+0.00000198*sind(2*F_deg+Omiga_deg)...
-0.00000172*sind(3*Omiga_deg)...
-0.00000087*T*sind(Omiga_deg))/3600;
gst_deg=epthilongGama_deg+GMST_deg;
end
该函数要调用nutationInLongitudeCaculate函数来计算章动角、黄赤交角等 ，该函数在 步骤4中给出；
同时 YMD2JD函数为 年月日转换为儒略日，具体说明 见公元纪年法（儒略历-格里高历）转儒略日年积日计算公式
代码如下：
function [jd] =YMD2JD(y,m,d) %函数公历年月日转儒略日
if size(y,2)==3 %如果输入为年月日数组
d=y(3);
m=y(2);
y=y(1);
end
if m<3
m=m+12;
y=y-1;
end
B=0;
if y>1582||(y==1582&&m>10)||(y==1582&&m==10&&d>=15)
B=2-floor(y/100)+floor(y/400);
end
jd=floor(365.25*(y+4712))+floor(30.6*(m+1))+d-63.5+B;
end
step 4 瞬时真天球坐标系 转到瞬时平天球 坐标系
xyz=ordinateSingleRotate(‘x’,-epthilongA_deg)*ordinateSingleRotate(‘z’,dertaPthi_deg)*ordinateSingleRotate(‘x’,dertaEpthilong_deg+epthilongA_deg)*xyz;
其中 epthilongA_deg 为平黄赤交角， dertaPthi_deg为黄经章动， dertaEpthilong_deg为交角章动， epthilong_deg为瞬时黄赤交角
这些角度的计算方法为
function [epthilongA_deg,dertaPthi_deg,dertaEpthilong_deg,epthilong_deg] = nutationInLongitudeCaculate(JD,accuracy)
%计算平黄赤交角，黄经章动、和交角章动、瞬时黄赤交角。
% 参数
% epthilongA_deg平黄赤交角
% dertaPthi_deg黄经章动
% dertaEpthilong_deg和交角章动
% epthilong_deg瞬时黄赤交角
%JD 儒略日
% accuracy 表示计算的精度要求 数值为’normal‘ 和’high‘ 或者用'N'和’H'表示一般精度和高精度 。默认为高精度计算
if nargin()<2
accuracy='h';
end
T=(JD-2451545)/36525;
T2=T*T;
T3=T2*T;
T4=T3*T;
T5=T4*T;
%计算月球平近点角lunarMeanAnomaly_deg (l_deg) 太阳平近点角SolarMeanAnomaly_deg(solarl_deg)
%月球纬度辐角lunarLatitudeAngle_deg(F_deg) 日月平角距diffLunarSolarElestialLongitude_deg(D_deg月球平黄经-太阳平黄经）
%月球升交点平黄经SolarAscendingNodeElestialLongitude_deg ( Omiga_deg)
l_deg = (-0.00024470*T4+0.051635*T3+31.8792*T2+1717915923.2178*T+485868.249036)/3600;
l_deg=mod(l_deg,360);
solarl_deg =(-0.00001149*T4+0.000136*T3-0.5532*T2+129596581.0481*T+1287104.79305)/3600;
solarl_deg=mod(solarl_deg,360);
F_deg =(0.00000417*T4-0.001037*T3-12.7512*T2+1739527262.8478*T+335779.526232)/3600;
F_deg=mod(F_deg,360);
D_deg = (-0.00003169*T4+0.006593*T3-6.3706*T2+1602961601.2090*T+1072260.70369)/3600;
D_deg=mod(D_deg,360);
Omiga_deg =(-0.00005939*T4+0.007702*T3+7.4722*T2-6962890.5431*T+450160.398036 )/3600;
Omiga_deg=mod(Omiga_deg,360);
basicAngle_deg=[l_deg solarl_deg F_deg D_deg Omiga_deg];
epthilongA_deg=(-0.0000000434*T5-0.000000576*T4+0.00200340*T3-0.0001831*T2-46.836769*T+84381.406)/3600;
epthilongA_deg=epthilongA_deg-floor(epthilongA_deg/360)*360;
switch lower(accuracy)
case {'h','high'} %IAU2000模型有77项
elestialLonNutation= elestialLonNutationCaculate();
dertaPthi_deg =-3.75e-08;
dertaEpthilong_deg =0.388e-3/3600;
for i = 1:77 %计算日月周期项
sumAngle_deg=0;
for j=1:5
sumAngle_deg=sumAngle_deg+elestialLonNutation(i,j)*basicAngle_deg(j);
end
sumAngle_deg=sumAngle_deg-floor(sumAngle_deg/360)*360;
dertaPthi_deg=dertaPthi_deg+((elestialLonNutation(i,6)+elestialLonNutation(i,7)*T)...
*sind(sumAngle_deg)+elestialLonNutation(i,8)*cosd(sumAngle_deg))*1e-7/3600;
dertaEpthilong_deg=dertaEpthilong_deg+((elestialLonNutation(i,9)+elestialLonNutation(i,10)*T)...
*cosd(sumAngle_deg)+elestialLonNutation(i,11)*sind(sumAngle_deg))*1e-7/3600;
end
case{'n','l','normal','low'}
Omiga_deg=basicAngle_deg(5);
F_deg=basicAngle_deg(3);
D_deg=basicAngle_deg(4);
solarl_deg=basicAngle_deg(2);
dertaPthi_deg=((-17.1996+0.01742*T)*sind(Omiga_deg)...
+(0.2062+0.00002*T)*sind(2*Omiga_deg)...
-(1.3187+0.00016*T)*sind(2*F_deg-2*D_deg+2*Omiga_deg)...
+(0.1426-0.00034*T)*sind(solarl_deg)...
-(0.2274+0.00002*T)*sind(2*F_deg+2*Omiga_deg))/3600;
dertaEpthilong_deg=((9.2025+0.00089*T)*cosd(Omiga_deg)...
-(0.0895-0.00005*T)*cosd(2*Omiga_deg)...
+(0.5736-0.00031*T)*cosd(2*F_deg-2*D_deg+2*Omiga_deg)...
+(0.0054-0.00001*T)*cosd(solarl_deg)...
+(0.0977-0.00005*T)*cosd(2*F_deg+2*Omiga_deg) )/3600;
otherwise
error(message(':unknownaccuracy:accuracy=',accuracy));
end
epthilong_deg=epthilongA_deg+dertaEpthilong_deg;
end
其中 elestialLonNutationCaculate() 函数 返回一个 77行11列 的 二维数组，函数 如下
function elestialLonNutation=elestialLonNutationCaculate()
elestialLonNutation = ...
[0 0 0 0 1 -1.72064161E+8 -174666 33386 9.2052331E+7 9086 15377;
0 0 2 -2 2 -1.3170906E+7 -1675 -13696 5.730336E+6 -3015 -4587;
0 0 2 0 2 -2.276413E+6 -234 2796 978459 -485 1374;
0 0 0 0 2 2.074554E+6 207 -698 -897492 470 -291;
0 1 0 0 0 1.475877E+6 -3633 11817 73871 -184 -1924;
0 1 2 -2 2 -516821 1226 -524 224386 -677 -174;
1 0 0 0 0 711159 73 -872 -6750 0 358;
0 0 2 0 1 -387298 -367 380 200728 18 318;
1 0 2 0 2 -301461 -36 816 129025 -63 367;
0 -1 2 -2 2 215829 -494 111 -95929 299 132;
0 0 2 -2 1 128227 137 181 -68982 -9 39;
-1 0 2 0 2 123457 11 19 -53311 32 -4;
-1 0 0 2 0 156994 10 -168 -1235 0 82;
1 0 0 0 1 63110 63 27 -33228 0 -9;
-1 0 0 0 1 -57976 -63 -189 31429 0 -75;
-1 0 2 2 2 -59641 -11 149 25543 -11 66;
1 0 2 0 1 -51613 -42 129 26366 0 78;
-2 0 2 0 1 45893 50 31 -24236 -10 20;
0 0 0 2 0 63384 11 -150 -1220 0 29;
0 0 2 2 2 -38571 -1 158 16452 -11 68;
0 -2 2 -2 2 32481 0 0 -13870 0 0;
-2 0 0 2 0 -47722 0 -18 477 0 -25;
2 0 2 0 2 -31046 -1 131 13238 -11 59;
1 0 2 -2 2 28593 0 -1 -12338 10 -3;
-1 0 2 0 1 20441 21 10 -10758 0 -3;
2 0 0 0 0 29243 0 -74 -609 0 13;
0 0 2 0 0 25887 0 -66 -550 0 11;
0 1 0 0 1 -14053 -25 79 8551 -2 -45;
-1 0 0 2 1 15164 10 11 -8001 0 -1;
0 2 2 -2 2 -15794 72 -16 6850 -42 -5;
0 0 -2 2 0 21783 0 13 -167 0 13;
1 0 0 -2 1 -12873 -10 -37 6953 0 -14;
0 -1 0 0 1 -12654 11 63 6415 0 26;
-1 0 2 2 1 -10204 0 25 5222 0 15;
0 2 0 0 0 16707 -85 -10 168 -1 10;
1 0 2 2 2 -7691 0 44 3268 0 19;
-2 0 2 0 0 -11024 0 -14 104 0 2;
0 1 2 0 2 7566 -21 -11 -3250 0 -5;
0 0 2 2 1 -6637 -11 25 3353 0 14;
0 -1 2 0 2 -7141 21 8 3070 0 4;
0 0 0 2 1 -6302 -11 2 3272 0 4;
1 0 2 -2 1 5800 10 2 -3045 0 -1;
2 0 2 -2 2 6443 0 -7 -2768 0 -4;
-2 0 0 2 1 -5774 -11 -15 3041 0 -5;
2 0 2 0 1 -5350 0 21 2695 0 12;
0 -1 2 -2 1 -4752 -11 -3 2719 0 -3;
0 0 0 -2 1 -4940 -11 -21 2720 0 -9;
-1 -1 0 2 0 7350 0 -8 -51 0 4;
2 0 0 -2 1 4065 0 6 -2206 0 1;
1 0 0 2 0 6579 0 -24 -199 0 2;
0 1 2 -2 1 3579 0 5 -1900 0 1;
1 -1 0 0 0 4725 0 -6 -41 0 3;
-2 0 2 0 2 -3075 0 -2 1313 0 -1;
3 0 2 0 2 -2904 0 15 1233 0 7;
0 -1 0 2 0 4348 0 -10 -81 0 2;
1 -1 2 0 2 -2878 0 8 1232 0 4;
0 0 0 1 0 -4230 0 5 -20 0 -2;
-1 -1 2 2 2 -2819 0 7 1207 0 3;
-1 0 2 0 0 -4056 0 5 40 0 -2;
0 -1 2 2 2 -2647 0 11 1129 0 5;
-2 0 0 0 1 -2294 0 -10 1266 0 -4;
1 1 2 0 2 2481 0 -7 -1062 0 -3;
2 0 0 0 1 2179 0 -2 -1129 0 -2;
-1 1 0 1 0 3276 0 1 -9 0 0;
1 1 0 0 0 -3389 0 5 35 0 -2;
1 0 2 0 0 3339 0 -13 -107 0 1;
-1 0 2 -2 1 -1987 0 -6 1073 0 -2;
1 0 0 0 2 -1981 0 0 854 0 0;
-1 0 0 1 0 4026 0 -353 -553 0 -139;
0 0 2 1 2 1660 0 -5 -710 0 -2;
-1 0 2 4 2 -1521 0 9 647 0 4;
-1 1 0 1 1 1314 0 0 -700 0 0;
0 -2 2 -2 1 -1283 0 0 672 0 0;
1 0 2 2 1 -1331 0 8 663 0 4;
-2 0 2 2 2 1383 0 -2 -594 0 -2;
-1 0 0 0 2 1405 0 4 -610 0 2;
1 1 2 -2 2 1290 0 0 -556 0 0];
end
step5 瞬时平天球坐标系转换为协议天球坐标系（J2000）
xyz=ordinateSingleRotate(‘Z’,zetaA)*ordinateSingleRotate(‘y’,-thitaA)*ordinateSingleRotate(‘z’,zA)*xyz;
该步主要涉及到 岁差角的计算。下面给出计算方法。
function [zetaA,thitaA,zA] = precessionAngle(JDTDB)
%zA、thitaA、zetaA为赤道岁差角,计算赤道岁差角
T=(JDTDB-2451545)/36525;
T2=T*T;
T3=T2*T;
%
% zetaA=(2306.218*T+0.30188*T2+0.017998*T3)/3600;
% zA=(2306.218*T+1.09468*T2+0.018203)/3600;
% thitaA=(2004.3109*T-0.42665*T2-0.041833*T3)/3600;
T4=T3*T;
T5=T4*T;
%
zetaA=(-0.0000003173*T5-0.000005971*T4+0.01801828*T3+0.2988499*T2+2306.083227*T+2.650545)/3600;
thitaA=(-0.0000001274*T5-0.000007089*T4-0.04182264*T3-0.4294934*T2+2004.191903*T)/3600;
zA=(-0.0000002904*T5-0.000028596*T4+0.01826837*T3+1.0927348*T2+2306.077181*T-2.650545)/3600;
end
函数中 JDTDB 为给定时刻 的地球动力学时对应的儒略日，其计算方法由步骤三中的函数给出。
参考链接：
[1]https://blog.csdn.net/hit5067/article/details/116894616
6517
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
