计算行星位置matlab程序

function a=Getposition(time,n) %n=1表示地球 n=2 表示火星 n=3表示ceres %Ceres轨道参数 ac= 2.7654089; ec=0.21432115; Mc=80.39176*0.01745329 ; ic=10.58663*0.01745329; CapitalOmegac= 0.0790083*0.01745329; Omegac=10.58663*0.01745329; Tc=1679.667171431635; %地球参数 ae = 1; ee = 0.01671022; Me = -2.48074*0.01745329; ie = 0; CapitalOmegae = 0; Omegae = 102.94719*0.01745329; Te = 365.2564; %火星参数 am = 1.52366231; em = 0.09341233; Mm = -4.55343205*0.01745329; im = 1.85061*0.01745329; CapitalOmegam = 49.57854*0.01745329; Omegam = 286.4623*0.01745329; Tm = 686.98; % a=5.20336301; % e=0.04839266; % M=19.65053*0.01745329; % i=1.30530*0.01745329; % w=100.55615*0.01745329; % q=274.1977*0.01745329; jD=JD(time(1),time(2),time(3)); orbitparam=[am,em,Mm,im,CapitalOmegam,Omegam,Tm]; orbitparae=[ae,ee,Me,ie,CapitalOmegae,Omegae,Te]; orbitparac=[ac,ec,Mc,ic,CapitalOmegac,Omegac,Tc]; if n==1 a=ElipCartHelioCent(jD,orbitparae); elseif n==2 a=ElipCartHelioCent(jD,orbitparam); elseif n==3 a=ElipCartHelioCent(jD,orbitparac); else a=[]; end end function time=ElipCartHelioCent(jd1,orbitpara) % 轨道参数等a半长轴，e椭率，M平近点角，i倾角，\[CapitalOmega]升交点黄经， % \[Omega]\近日点幅角，后面加1表示角度制，不加1表示弧度制，填入参数历元year0、month0、day0 % 以及你所需要的日期year1、month1、\day1、hour1 % a=5.20336301; % e=0.04839266; % M1=19.65053; % i1=1.30530; % CapitalOmega1=100.55615; % Omega1=274.1977; a=orbitpara(1); e=orbitpara(2); M=orbitpara(3); i=orbitpara(4); CapitalOmega=orbitpara(5); Omega=orbitpara(6); T1=orbitpara(7); year0=2000; month0=1; day0=1; % T1=4.6*365.2564; mass=1898.6; massS = 1988910; if T1 > 0 T=T1; else T= 2*pi*(massS/(massS + mass))^0.5*a^1.5/0.0172021; end jd0 = JD(year0, month0, day0); ETinverse = ET(i, CapitalOmega, Omega); t = jd1 - jd0 + M/2/pi*T; t1 = t - T*floor(t/T + 0.5); %*计算位置角\[CurlyPhi]的值,得到行星运动平面坐标系下坐标, %经欧拉矩阵变换后，变换到日心黄道坐标系的坐标，并求出日心黄经黄纬*) CurlyPhi= fzero(@(x) F(x),[-pi,pi]); r = a*(1 - e^2)/(1 + e*cos(CurlyPhi)); Coordinate1 = [r*cos(CurlyPhi), r*sin(CurlyPhi), 0]; time= (ETinverse*Coordinate1')'; function f=F(fai) f=atan(sqrt((1 - e)/(1 + e))*tan(fai/2))/pi - e*(1 - e^2)^0.5*sin(fai)/(1 + e*cos(fai))/2/pi - t1/T; end function et=ET(ii, CapitalOmega2, Omega2) E11 = cos(CapitalOmega2)*cos(Omega2)- cos(ii)*sin(CapitalOmega2)*sin(Omega2); E12 = -cos(CapitalOmega2)*sin(Omega2)-cos(ii)*sin(CapitalOmega2)*cos(Omega2); E13 = sin(ii)*cos(CapitalOmega2); E21 = sin(CapitalOmega2)*cos(Omega2) + cos(ii)*cos(CapitalOmega2)*sin(Omega2); E22 = -sin(CapitalOmega2)*sin(Omega2) +cos(ii)*cos(CapitalOmega2)*cos(Omega2); E23 = -sin(ii)*cos(CapitalOmega2); E31= sin(ii)*sin(Omega2); E32 = sin(ii)*cos(Omega2); E33 = cos(ii); et=[E11, E12, E13; E21, E22, E23; E31, E32, E33]; end end %Ceres轨道参数 ac= 2.7654089; ec=0.21432115; Mc=80.39176*0.01745329 ; ic=10.58663*0.01745329; CapitalOmegac= 0.0790083*0.01745329; Omegac=10.58663*0.01745329; Tc=1679.667171431635; %地球参数 ae = 1; ee = 0.01671022; Me = -2.48074*0.01745329; ie = 0; CapitalOmegae = 0; Omegae = 102.94719*0.01745329; Te = 365.2564; yeare = 2000; monthe = 1; daye = 1; %火星参数 am = 1.52366231; em = 0.09341233; Mm = -4.55343205*0.01745329; im = 1.85061*0.01745329; CapitalOmegam = 49.57854*0.01745329; Omegam = 286.4623*0.01745329; yearm = 2000; monthm = 1; daym = 1; Tm = 686.98; massm = 0.64185; a=5.20336301; e=0.04839266; M=19.65053*0.01745329; i=1.30530*0.01745329; w=100.55615*0.01745329; q=274.1977*0.01745329; orbitparam=[am,em,Mm,im,CapitalOmegam,Omegam,Tm]; orbitparae=[ae,ee,Me,ie,CapitalOmegae,Omegae,Te]; orbitparac=[ac,ec,Mc,ic,CapitalOmegac,Omegac,Tc]; for jD=2450000:2450009 ElipCartHelioCent(jD,orbitparac) end function jd=JD(yeaR, montH, daY) %日期转化为儒略日 a1 = floor((14 - montH)/12); y = yeaR + 4800 - a1; m = montH + 12*a1 - 3; jd=daY + floor((153*m + 2)/5) + 365*y + floor(y/4) - 32083; end