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 the given time. The specific conversion is as follows
step1 WGS 84 converts to the protocol earth coordinate system.
function XYZ_m = wgs842ECEF(BLH_deg_m) %Geodetic latitude and longitude elevation BLH Transformation with earth fixed coordinate system BLH-->XG %ECEF, Earth Centered Earth Fixed;The reference plane of the earth fixed coordinate system: the plane of the equator, that is, passing through the center of the earth and with the center of the earth CIO The plane perpendicular to the point line. % x The axis is the intersection of the reference plane and the Greenwich plane, z The axis is geocentric CIO Point. %parameter %B L H The input parameters are geodetic latitude, geodetic longitude and geodetic elevation respectively N*3 matrix % XYZ It is the of the earth fixed coordinate system x y z，The output parameter is N*3 matrix %ae = 6378137; %2000 Field half axis of coordinate system %ee = 0.0818191910428; %First eccentricity ae = 6378137; %wgs84 Field half axis of coordinate system ee = 0.00669437999013 ; %First eccentricity 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
The step 2 protocol converts the earth coordinate system into the instantaneous earth coordinate system
This mainly involves querying XP and YP at a given time, which can be obtained by querying the EOP file. The download address is as follows http://celestrak.com/spacedata/
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 z The shaft rotates counterclockwise
% 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
end
step 3 transformation from instantaneous earth coordinate system to 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. The dut1 and dat parameters are in the EOP file. The download address of EOP file is CelesTrak: EOP and Space Weather Data
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 In the satellite time service signals of various countries, 0 will be used.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];
%Let's calculate Greenwich Mean stellar time GMST; In degrees
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 nutation angle, yellow red cross angle, etc. this function is given in step 4;
At the same time, YMD2JD function is to convert year, month and day into Julian day. For details, see the conversion of year method (Julian calendar Gregorian calendar) to Julian day_$ {Wang Xiaojian}'s blog - CSDN blog_ Annual product day calculation formula
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;
end
step 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+epthilongA_deg)*xyz;
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 elestialonnutationcaculate() 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=(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
The time of TDB in the given time function of Julian is the corresponding time of the earth.
Reference link:
[1]https://blog.csdn.net/hit5067/article/details/116894616