clc;clear all;close all
%Posición del Sol en función del día del año --> Ascensión recta del Sol
%cálculo del día juliano
A=2017; %año
Me=9; %mes
D=1:1:30;%día
n=length(D);
for k=1:n
d(k)=367*A-floor((7*(A+((Me+9)/12)))/4)+floor(275*Me/9+D(k))-730530;
%Elementos orbitales del Sol
w(k)=282.9404+4.70935*10^-5*d(k);%longitud del perihelio
a=1; %distancia Tierra Sol en AU
format long
ecalc(k)=0.016709-1.151*10^-9*d(k);
e(k)=(round(ecalc(k)*10^7))/(10^7);
M(k)=(round((356.047+0.9856002585*d(k))*10^4))/(10^4); %anomalía media
if M(k)<0
Mpos(k)=round((M(k)+(floor(-M(k)/360)+1)*360)*10^7)/10^7;
M(k)=Mpos(k);
end
%Inclinación del plano de la eclíptica
oblecl(k)=23.4393-3.563*10^-7*d(k);
%Anomalía excéntrica
E(k)=M(k)+(180/pi)*e(k)*sin(M(k)*pi/180)*(1+e(k)*cos(M(k)*pi/180));
%Coordenadas rectangulares del Sol en el plano de la eclíptica. El Sol se
%mueve siempre en el plano de la eclíptica
x(k)=cos(E(k)*pi/180)-e(k);
y(k)=sin(E(k)*pi/180)*sqrt(1-e(k)*e(k));
%Convertimos a distancia y anomalía verdadera
r(k)=round((sqrt(x(k)^2+y(k)^2))*10^6)/(10^6); %distancia
tetha(k)=round(((atan2(y(k),x(k)))*180/pi)*10^4)/(10^4); %anomalía verdadera
%Longitud
lon(k)=tetha(k)+w(k);
if lon(k)>360
loncorr(k)=lon(k)-360;
lon(k)=loncorr(k);
end
%Coordenadas rectangulares del Sol
x(k)=r(k)*cos(lon(k)*pi/180);
y(k)=r(k)*sin(lon(k)*pi/180);
z=0;
%Convertir a RA y declinación
RA(k)=(atan2(y(k),x(k)))*180/pi
%decl=atan2(z,sqrt(x^2+y^2))
%ventana de lanzamiento --> escoger RAAN perpendicular a RA, sumarle 90 deg
RAAN(k)=RA(k)+90;
end
e-REdING. Biblioteca de la Escuela Superior de Ingenieros de Sevilla.
