-
Notifications
You must be signed in to change notification settings - Fork 1
/
analyticSun2.m
70 lines (56 loc) · 1.98 KB
/
analyticSun2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%{
Author: Shane Kramer/Professor Edward Brown (UCCS)
Course: SPCE 5025 Fundamentals Of Astronautics
Date: 03.07.15
---------------------------------------------------
This function returns a low precision formula for the
Sun's position given an Julian Date (Universal Time).
Based off of 2013 Astonomical Almanac.
%}
function [ rsun ] = analyticSun()
% Establish constants for our model
% ----------------------------------------------------
global TDB;
d2r = pi/180; % degrees to radians conversion
twopi = 2*pi;
AU = 149597870700; % meters (IAU values, 2012)
% Days from J2000 epoch (UT)
% ----------------------------------------------------
n = TDB-2451545.0;
% Mean Longitude of sun, corrected for aberration
% ----------------------------------------------------
L = (280.460+0.9856474*n)*d2r; % radians
L = mod(L, twopi);
if(L<0)
L = L+twopi;
end
% Mean Anomaly
% ----------------------------------------------------
g = (357.528+0.9856003*n)*d2r; % radians
g = mod(g, twopi);
if(g<0)
g = g+twopi;
end
% Ecliptic Longitude and latitude
% ----------------------------------------------------
lambda = L + 1.915*d2r*sin(g) + 0.020*d2r*sin(2*g);
beta = 0.0;
% Obliquity of Ecliptic
% ----------------------------------------------------
epsilon = (23.439-0.0000004*n)*d2r;
% Right Ascension and declination (not used for rsun)
% ----------------------------------------------------
f=180/pi;
t = (tan(epsilon/2))^2;
alpha = lambda - f*t*sin(2*lambda) + (f/2)*t^2*sin(4*lambda);
delta = asin(sin(epsilon)*sin(lambda));
% Distance from Earth to Sun in AU
% ----------------------------------------------------
R = 1.00014-0.01671*cos(g)-0.00014*cos(2*g);
% ECI coordinates of sun in AU
% ----------------------------------------------------
rsun = [R*cos(lambda); R*cos(epsilon)*sin(lambda); R*sin(epsilon)*sin(lambda)];
% Convert to meters
% ----------------------------------------------------
rsun = AU*rsun;
end