The transformation from wgs84 coordinate system to J2000 coordinate system mainly involves the mutual transformation of coordinates. Generally, the given WGS coordinates are t and BLH at a given time. Specific conversion Ru Xia
step1 WGS 84 is converted to the protocol earth coordinate system.
function XYZ = wgs842xyz(BLH)
%UNTITLED2 Geodetic coordinate system to North Tiandong coordinate system
% Detailed description here thita When it is a local star
ae = 6378137; %wgs84 Field half axis of coordinate system
ee = 0.00669437999013 ; %First eccentricity
%ee=0.00669437999; %84 system
thitaG = pi/2; %Angle between Greenwich surface and vernal equinox surface
%Normalized dimension(Earth)
% Horizon coordinates
thita = thitaG+BLH(2);
temp =ae/sqrt(1-ee^2*sin(BLH(1)));
temp2 = (temp+BLH(3))*cos(BLH(1));
x = temp2*cos(thita) ; %east
y = temp2*sin(thita) ; %north
z = (temp*(1-ee^2)+BLH(3))*sin(BLH(1)) ;%day
CL2I =[-sin(thita) -cos(thita)*sin(BLH(1)) cos(thita)*cos(BLH(1));
cos(thita) -sin(thita)*sin(BLH(1)) sin(thita)*cos(BLH(1));
0 cos(BLH(1)) sin(BLH(1)) ]; %Transformation matrix of horizon to earth inertial system
XYZ=CL2I * [ x y z]' ;% X Point to the spring equinox for the center of the earth, Y 90 degrees from the earth's center to the equatorial plane, Z Point north for the center of the earth
endThe step 2 # protocol converts the earth coordinate system to the instantaneous coordinate system
This mainly involves querying XP and YP at a given time, which can be obtained by querying the IERS bulletin. If the calculation requirements are not high, this type of coordinate transformation can be ignored.
earthFixedXYZ=ordinateSingleRotate('x',yp)*ordinateSingleRotate('y',xp)*earthFixedXYZ;
Where ordinateSingleRotate is the coordinate conversion function. The first parameter of the function is the given coordinate axis, the second parameter is the rotation angle, and the third parameter is the measurement unit of the rotation angle. This function will be called many times later.
function [R] = ordinateSingleRotate( axis,angle_deg,angleType)
%Axis rotation
% axis Represents the axis of rotation
% '1'perhaps'x' Indicates around x The shaft rotates counterclockwise
% '2'perhaps'y' Indicates around y The shaft rotates counterclockwise
% '3'perhaps'z' Indicates around x Counterclockwise rotation axis
% angle_deg Represents the angle of rotation
% angleType Represents the angle type, which can be 'rad'and'deg'，default 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
endstep 3: transform the instantaneous earth coordinate system into the instantaneous celestial coordinate system
xyz= ordinateSingleRotate('z',-gst_deg)*earthFixedXYZ;
This step mainly involves the calculation of Greenwich star time angle. There are many calculation methods for Greenwich star time angle. A more accurate calculation method is given below.
function [gst_deg,JDTDB ] = utc2gst(UTC,dUT1,dAT)
%take utc Time is converted to Greenwich star time
%parameter
%utc Format is y m d among d The value of is decimal, and h m s Press( sec/3600+min/60+h)/24 Convert to decimal and add up day upper
%dUT1 by ut1-utc If the value does not exceed plus or minus 1 second, check iers Available values
%UNTITLED When calculating local stars,The return value is in hours and seconds
% UTC Universal coordinated time coordinates the time difference between earth time and atomic time by using leap seconds
%dAT Run seconds
%TAI International Atomic Time is the most accurate time TAI=UTC+dAT; %International Atomic Time
%TT Earth time TT = TAI+32.184 By 2014;
%TDT Geodynamic time
% ET Almanac time
%Earth time=Geodynamic time=Almanac time
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 It is universal time. Because the rotation of universal time is uneven, it is different from UTC Time has dUT1 The difference between, dUT1 Countries will send 0 signals in the satellite.1 The accuracy of seconds is given IERS After processing, it will be 1 e-5 The accuracy is given.
dT=32.184+dAT-dUT1;
JDTT=JDUT1+dT/86400;
% JDTT= YMD2JD(UTC(1),UTC(2),UTC(3))+dUT1/86400;
%First calculate TDB,TDB When it is centroid dynamics, it is the time scale in the ephemeris of celestial bodies such as the sun, moon and planets
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;
%Calculate below UT2,UT2 Yes UT1 Based on the correction of the earth's periodic seasonal changes, the world time is obtained
% %according to UT1 calculation Tb,Tb In Bessel years, from calendar year B2000.0 Start 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];
%The rotation angle of the earth is calculated below thita，thita In circles,
%Let's calculate Greenwich Mean stellar time GMST; The unit is seconds, using Zhang Hongbo's textbook method
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; %since J2000 The number of turns the earth has made;
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);
%Calculate the error caused by the nutation of the right ascension eclipticObliquitygama，Unit: arcsec
%Calculate the lunar perigee angle, solar perigee angle, lunar latitude, radiation angle, Solar Lunar perigee distance(Moon flat yellow meridian-Sun flat yellow meridian) lunar ascending node flat yellow meridian, unit: arcsec
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
This function calls the nutationinlongitudecculate function to calculate the nutation intersection, yellow red cross angle, etc. this function is given in step 4;
At the same time, YMD2JD function converts year, month and day into Julian day. See for details https://blog.csdn.net/qq_24172609/article/details/112244135?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522162115452716780357279082%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=162115452716780357279082&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~first_rank_v2~rank_v29-3-112244135.first_rank_v2_pc_rank_v29&utm_term=%E5%85%AC%E5%85%83%E7%BA%AA%E5%B9%B4%E6%B3%95&spm=1018.2226.3001.4187
The code is as follows
function [jd] =YMD2JD(y,m,d) %Function Gregorian calendar month day to Julian day
if size(y,2)==3 %If the input is an array of year, month and day
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;
endstep 4 switch from the instantaneous celestial coordinate system to the instantaneous celestial coordinate system
xyz=ordinateSingleRotate('x',-epthilongA_deg)*ordinateSingleRotate('z',dertaPthi_deg)*ordinateSingleRotate('x',dertaEpthilong_deg)*xyz;
Of which, epthilongA_deg is the intersection angle of flat yellow and red, dertaPthi_deg is the nutation of the Yellow Sutra, derta epthilong_deg is angular nutation, epthilong_deg is the instantaneous yellow red cross angle
The calculation method of these angles is
function [epthilongA_deg,dertaPthi_deg,dertaEpthilong_deg,epthilong_deg] = nutationInLongitudeCaculate(JD,accuracy)
%Calculate the horizontal yellow red cross angle, yellow meridians nutation, sum angle nutation and instantaneous yellow red cross angle.
% parameter
% epthilongA_deg Flat yellow red angle
% dertaPthi_deg Huang Jing nutation
% dertaEpthilong_deg Angular nutation
% epthilong_deg Instantaneous yellow red cross angle
%JD Julian day
% accuracy Indicates the required accuracy of the calculation. The value is'normal' and'high' Or use'N'and'H'Indicates general precision and high precision. The default is high-precision calculation
if nargin()<2
accuracy='h';
end
T=(JD-2451545)/36525;
T2=T*T;
T3=T2*T;
T4=T3*T;
T5=T4*T;
%Calculation of lunar perigee angle lunarMeanAnomaly_deg (l_deg) Solar perigee angle SolarMeanAnomaly_deg(solarl_deg)
%Lunar latitude and radian lunarLatitudeAngle_deg(F_deg) Horizontal angular distance between sun and moon diffLunarSolarElestialLongitude_deg(D_deg Moon flat yellow meridian-Taiyangpinghuang meridian)
%Lunar ascending node flat yellow meridian 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 There are 77 models
elestialLonNutation= elestialLonNutationCaculate();
dertaPthi_deg =-3.75e-08;
dertaEpthilong_deg =0.388e-3/3600;
for i = 1:77 %Calculate daily and monthly cycle items
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
The # electiallonnutationcalculate() function returns a two-dimensional array of 77 rows and 11 columns. The function is as follows
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: convert the instantaneous celestial coordinate system into the protocol celestial coordinate system (J2000)
xyz=ordinateSingleRotate('Z',zetaA)*ordinateSingleRotate('y',-thitaA)*ordinateSingleRotate('z',zA)*xyz;
This step mainly involves the calculation of precession angle. The calculation method is given below.
function [zetaA,thitaA,zA] = precessionAngle(JDTDB) %zA,thitaA,zetaA Equatorial precession angle,Calculate equatorial precession angle T=(JD-2451545)/36525; T2=T*T; T3=T2*T; %Zhang Yuxiang satellite orbit measurement technology and method % 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; %Zhang Hongbo 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
In the function, JDTDB is the Julian day corresponding to the geodynamic time at a given time, and its calculation method is given by the functions of step three.