Simulate semi-major axis change in a dawn-dusk sun-synchronous orbit.

Uses a simple specular model for the solar pressure force.
Integration tolerance must be tight enough to produce circular
precessing orbit when sail is edge-on.
Multiple sail angle laws are implemented in the code. Scroll down to
select a method or add your own.

Default sail sizes are from 40 to 60 m square. The total spacecraft
mass is on the order of 500 kg.

If one axis of the sail is aligned with LVLH, which alleviates gravity
gradient torques, the semi-major axis change is on the order of 40 km
in 15 days. The sun-synch orbit is closely maintained.
The oscillations in semi-major axis change are due to the inclusion of
the J2 gravity term, which is necessary for a sun-synch orbit.

Since version 7.
------------------------------------------------------------------------
See also Inertias, Cone, Eul2Q, QLVLH, QMult, QTForm, NPlot, Plot2D,
TimeLabl, Cross, Dot, RK45, Unit, Date2JD, LoadGravityModel, El2RV,
PltOrbit, RV2El, Eclipse, SunV1, Accel, UToConeClock, delta
------------------------------------------------------------------------

Contents

%--------------------------------------------------------------------------
%   Copyright (c) 2006 Princeton Satellite Systems, Inc.
%   All rights reserved.
%--------------------------------------------------------------------------

clear d

Set up the physical parameters

%-------------------------------
Asail         = 2000; % 3600 1600]; % Sail area in m2
arealDensity  = 0.025; % Sail subsystem areal density in g/m2
mSail         = Asail*arealDensity;   % Sail subsystem mass in kg
Isail         = Inertias(mSail(1),[60 60],'plate'); % Inertia
mSC           = 450; % mass of spacecraft, kg

Constants

%----------
m2km          = 0.001;
cLight        = 3e8;  % speed of light, m/s
flux          = 1367; % Solar flux at 1 AU, W/m2
dW            = 360/(365*86400); % deg/s
rE            = 6378; % km
gravityModel  = LoadGEM( 1 );

nM            = 36; % Tesseral harmonics
nN            = 36; % Zonal harmonics

Epoch

%------
jD0           = Date2JD([2010 3 15, 16 0 0]);
[uS, rS]      = SunV1( jD0 );

Orbit geometry

%---------------
sma   = 7978.1;
meanm = sqrt( gravityModel.mu/sma^3 );
J2    = 0.001082;
inc   = acos( -dW*pi/180/ (3*meanm*J2*rE^2/(2*sma^2)) );
RA    = atan2(uS(1),uS(2));
el0   = [sma;inc;RA;0;0;0];
PltOrbit( el0, jD0, [0;0;1],'earth' )

% ensure one or other is vector, not both
if length(Asail) > 1 && length(mSC) > 1
  error('Please make either Asail or mSC a vector, but not both.')
end

ans = 
  Figure (PlotPSS) with properties:

      Number: 1
        Name: 'earth Orbit'
       Color: [0.94 0.94 0.94]
    Position: [560 528 560 420]
       Units: 'pixels'

  Use GET to show all properties

Total mass

%-----------
mass      = mSail + mSC;
mSim      = length(mass);

Number of sim steps and duration

%---------------------------------
nDays     = 2;
dT        = 150;
nSim      = ceil(nDays*86400/dT);

Plotting array

%---------------
nPlot     = zeros(mSim,nSim);
aPlot     = zeros(mSim,nSim);
iPlot     = aPlot;
WPlot     = aPlot;
ePlot     = aPlot;
accPlot   = aPlot;
conePlot  = aPlot;
clckPlot  = aPlot;
sunAPlot  = aPlot;
eclPlot   = zeros(mSim,nSim);

Run the simulation(s)

%----------------------
for j = 1:mSim

  % Assume a specular sail and compute acceleration
  %------------------------------------------------
  a0      = 2*Asail(j)*flux/cLight/mass(j)*m2km;

  % Initial conditions
  %-------------------
  [r0,v0] = El2RV( el0, [], gravityModel.mu );
  x       = [r0;v0];
  jD      = jD0;

  disp('-----------------------------')
  disp(['Simulation ' num2str(j)])
  disp('-----------------------------')
  kS = 1;

  [ ratioRealTime, tToGoMem ] = TimeGUI( nSim );

  for k = 1:nSim

    [ ratioRealTime, tToGoMem ] = TimeGUI( nSim, k, tToGoMem, ratioRealTime, dT );

    % Local variables
    %----------------
    r = x(1:3);
    v = x(4:6);

    % Control (thrust vector determination)
    %--------------------------------------

First determine sail unit vector.

    [uS, rS] = SunV1( jD );
    method = 2;
    switch method
      case 1

Method A: point to sun and check results

        u     = -uS;
        alpha = 0;
        delta = 0;

      case 2

Method B: align towards v, semi-major axis grows

Constant cone angle Clock angle rotates through 2pi each orbit The method creates an (unmodeled) gravity gradient torque

        alpha           = atan(1/sqrt(2));
        u               = cos(alpha)*-uS + sin(alpha)*Cross(uS,Cross(Unit(v),uS));
        [cone, delta]   = UToConeClock( u, r, v, uS );

      case 3
        % Method C: Align with LVLH and then yaw about z
        % Cone and clock angles oscillate
        % Orbit follows sun-synch
        qLVLH           = QLVLH( r, v );
        % Align body x axis to orbit normal
        qLVLHToBody     = Eul2Q([0;0;-pi/2]);
        q0              = QMult( qLVLH, qLVLHToBody );
        % Rotate about z
        cone            = atan(1/sqrt(2));
        q1              = QMult( q0, Eul2Q([0;0;cone]) );
        u               = QTForm( q1, [1;0;0] );
        [alpha, delta]  = UToConeClock( u, r, v, uS );
      case 4

Method D: point edge-on and check results (orbit precession)

        u               = -uS;
        alpha           = pi/2;
        delta           = 0;

    end

    % Check for eclipses (should be permanently illuminated)
    [n,eclPlot(j,k)] = Eclipse( r, rS*uS );
    accel          = a0*cos(alpha)^2;
    d.a            = accel*u;

    % Plotting
    el             = RV2El( r, v, gravityModel.mu )';
    aPlot(j,k)     = el(1);
    iPlot(j,k)     = el(2)*180/pi;
    WPlot(j,k)     = el(3)*180/pi;
    ePlot(j,k)     = el(5);
    accPlot(j,k)   = accel;
    conePlot(j,k)  = alpha;
    clckPlot(j,k)  = delta;
    sunAPlot(j,k)  = acos(Dot(uS,Unit(Cross(r,v))));

    % Integrate the trajectory one step
    [x,hLast] = RK45( 'FOrbCartH', x, dT, dT, 0.01*dT, 1e-10,...
		0, d.a, gravityModel, nN, nM, jD );
    jD = jD + dT/86400;

  end
end

TimeGUI( 'close' );

-----------------------------
Simulation 1
-----------------------------

Plot the results

tSec    = (0:(k-1))*dT;
[t, tL] = TimeLabl( tSec );

vec1    = 1:mSim;
vec2    = mSim+1:2*mSim;
vec3    = 2*mSim+1:3*mSim;

Plot2D(t,[aPlot;iPlot;ePlot],tL,{'a' 'i (deg)' 'e'},['Elements'],'lin',...
  {vec1,vec2,vec3})
W2 = WPlot(:,1) + repmat(dW*tSec,mSim,1);
Plot2D(t,[WPlot;W2],tL,'Right Ascension','Orbit Progression Against Analytical');
legend(num2str(sqrt(Asail)',3))
Plot2D(t,sunAPlot*180/pi,tL,'Sun Angle (deg)','Angle between orbit normal and the sun')
Plot2D(t,accPlot*1e6,tL,'Acceleration (mm/s^2)')
Plot2D(t,[conePlot;clckPlot],tL,{'\alpha' '\delta'},'Angles','lin',...
  {vec1,vec2})
Plot2D(t,eclPlot,tL,'Eclipse Type','Eclipse Type (0 = none, 3 = total)')

%--------------------------------------


Published with MATLAB® R2019a