Small Satellites

Home » Posts tagged 'Matlab' (Page 3)

Tag Archives: Matlab

Kepler’s equation, Iterative and Non-Iterative Solver Comparison

Kepler’s equation M = E-e*sin(E) M – Mean anomaly [0..pi] rad e – Eccentricity [0..1] E – Eccentric anomaly [rad]

clear all;
clc;

Non-Iterative Solver

tic;
j = 1;
for M = 0:0.01:pi
    i = 1;
    for e = 0:0.01:1
        a   =  (1 - e)*3/(4*e+0.5);
        b   =  -M/(4*e+0.5);
        y   =  (b^2/4 + a^3/27)^0.5;
        x   =  (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
        w   =  x - 0.078*x^5/(1 + e);
        E   =  M + e*(3*w - 4*w^3);
        % Two times applying 5th-order Newton correction
        for k = 1:2
            f   =  (E - e*sin(E)- M);
            fd  =  1 - e*cos(E);
            f2d =  e*sin(E);
            f3d = -e*cos(E);
            f4d =  e*sin(E);
            E = E -f/fd*(1 + f*f2d/(2*fd^2) + f^2*(3*f2d^2 - fd*f3d)/(6*fd^4)+...
                (10*fd*f2d*f3d - 15*f2d^3-fd^2*f4d)*f^3/(24*fd^6));
        end
        E_nit(i,j) = E;
        i = i + 1;
    end
    j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 0.149860 [s]

Newton Iterative Solver

tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.01:pi
    i = 1;
    for e = 0:0.01:1
        En  = M;
        Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
        while ( abs(Ens-En) > eps )
            En = Ens;
            Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
        end;
        E_it(i,j) = Ens;
        i = i + 1;
    end
    j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 0.893673 [s]
Error = log10(abs(E_it - E_nit));
e = 0:0.01:1;
M = (0:0.01:pi);
contourf (M,e,Error);
xlabel('M - Mean anomaly [rad]');
ylabel('e - Eccentricity');
title('Non-Iterative Solver,Error in E calculation on Log scale [rad] ');
colorbar;
Kepler Equation

[6] S. Mikkola, “A Cubic Approximation for Kepler’s Equation,” Celestial
 Mechanics and Dynamical Astronomy,Vol. 40, 1987, pp. 329–334.

 

Greenwich Sidereal Time

This code is a MATLAB script that can be used to calculate Greenwich Sidereal Time, Julian Day

clc;
clear all;
% Input Date: April 11,2013.  UT time 20:11:30
year  = 2013;  month = 4; day   = 11;
hour  = 20; min   = 11; sec   = 30;
UT = hour + min/60 + sec/3600;
J0 = 367*year - floor(7/4*(year + floor((month+9)/12))) ...
    + floor(275*month/9) + day + 1721013.5;

JD = J0 + UT/24;              % Julian Day
fprintf('Julian day = %6.4f [days] \n',JD);
JC = (J0 - 2451545.0)/36525;
Julian day = 2456394.3413 [days]

JC is Julian centuries between the Julian day J0 and J2000(2,451,545.0) Greenwich sidereal time at 0 hr UT can be calculated by this equation [Seidelmann,1992]

GST0 = 100.4606184 + 36000.77004*JC + 0.000387933*JC^2 - 2.583e-8*JC^3; %[deg]
GST0 = mod(GST0, 360);  % GST0 range [0..360]
fprintf('Greenwich sidereal time at 0 hr UT %6.4f [deg]\n',GST0);
Greenwich sidereal time at 0 hr UT 199.3719 [deg]

Greenwich sidereal time at any other UT time

GST = GST0 + 360.98564724*UT/24;
GST = mod(GST, 360);  % GST0 range [0..360]
fprintf('Greenwich sidereal time at UT[hours] %6.4f [deg]\n',GST);
Greenwich sidereal time at UT[hours] 143.0761 [deg]

Best Fit Solution,Straight Line Fit

%https://smallsat.wordpress.com/

clear all
clc
T = [0,20,100,100,200,400,400,800,1000,1200,1400,1600];
R = [10.96,10.72,14.1,14.85,17.9,25.4,26,40.3,47,52.7,58,63];
plot(T,R,'r*');
title('Straight Line Fit');
xlabel('Temperature (deg C)');
ylabel('Resistivity (Ohm cm)');
sy = 0;
sx = 0;
sxy = 0;
sx2 = 0;
n = 12;
for i=1:n
sy = sy +R(i);
sx = sx + T(i);
sxy = sxy + T(i)*R(i);
sx2 = sx2 + T(i)*T(i);
end
a = (sy*sx2-sx*sxy)/(n*sx2-sx^2);
b = (n*sxy-sx*sy)/(n*sx2-sx^2);
fprintf('\n\n Best Fit Coefficients')
fprintf('\n x = %6.4f T + %6.4f \n',b,a)
cm_fit =a+b*T;
hold on;
plot(T,cm_fit,'r');
legend('Data Points','Best Fit Solution','location','northwest');
 Best Fit Coefficients
 x = 0.0337 T + 11.4695

 

Published with MATLAB® 7.10

%d bloggers like this: