A spacecraft with rate or hysteresis damping.

You can try different kinds of damping. One is 3-axis damping that applies a damping constant to each axis rate. You should find a constant that approximates whatever damper you have. If everything is working properly the spacecraft will align its torque rod with the Earth's magnetic field.

This simulation does not have any disturbances. Run for a small number of orbits to clearly see the hysteresis loops. They are minor loops since the Earth's field is not strong enough to take the rod to saturation. Run for several days to damp the rates close to zero.

The second type is magnetic hysteresis damping. Read the headers for RHSHysteresisDamper.m, TorqueHysteresis.m, and MagneticHysteresis.m MagneticHysteresis.m for more information.

Things to try:

     1. Try both types of damping.
     2. Determine a damping constant that produces the same results as
        hysteresis damping.
     3. Try different orbits.
     4. Try different dipoles.
     5. With no damping determine the oscillation frequency of the
        spacecraft.
     6. See what happens with a rotation rate about the dipole axis.
     7. Remove the dipole and initialize with nonzero rates instead.

Since version 2014.1

----------------------------------------------------------------------
See also BFromHHysteresis, RHSHysteresisDamper, RHSRigidBodyWithDamping,
HysteresisOutput
----------------------------------------------------------------------

Constants
Simulation parameters
Damper parameters
Run the simulation
Plot

%--------------------------------------------------------------------------
%   Copyright (c) 2014 Princeton Satellite Systems, Inc.
%   All Rights Reserved
%--------------------------------------------------------------------------

Constants

%-----------
rE                  = Constant('equatorial radius earth');
degToRad            = pi/180;
mu0                 = 4e-7*pi;

Simulation parameters

%-----------------------
d                   = struct;
nOrbits             = 30;	 % Number of orbits for the simulation
d.jD0               = JD2000; % Start date
d.mu                = Constant('mu earth'); % Earth's gravitational parameter
altitude            = 300; % Orbit altitude (km)
inclination         = 45*degToRad; % Orbit inclination (rad)
d.dampingType       = 2;        % 1 is 3-axis damping 2 is a damping function

% Spacecraft with a dipole from a torque rod or permament magnet
%---------------------------------------------------------------
d.dipole            = [1;0;0]; % in ATM^2 along the x-axis
d.uDipole           = Unit(d.dipole);
[d.inertia, d.mass]	= InertiaCubeSat( '3U' );
d.inertia           = d.inertia(1,1)*eye(3); % Make spherical to simplify

% Compute the orbit
%------------------
el                  = [rE+altitude inclination 0 0 0 0];
[r,v]               = El2RV(el); % Get the starting position and velocity vectors
q                   = QLVLH(r,v); % Quaternion
[bI,bIDot]          = BDipole(r,d.jD0,v);

Damper parameters

%-------------------
if( d.dampingType == 1 )
  d.dampingData       = 1e-6;     % Damping constant
  d.dampingFun        = [];       % Not used
  z                   = [];       % Not used
else
  d.dampingData.Br	= 0.004;      % Remanence (T)
  d.dampingData.Bs	= 0.025;     % Saturation flux density (T)
  d.dampingData.Hc	= 12;      % Coercive force (A/m)
  % Damper rod unit vectors
  d.dampingData.u     = [[0;1;0] [0;1;0] [0;1;0]...
                         [0;0;1] [0;0;1] [0;0;1]];
  % Dimensions are radius 1 mm by 95 mm
  d.dampingData.v     = pi*0.001^2*0.095*ones(1,size(d.dampingData.u,2));
  d.dampingFun        = @RHSHysteresisDamper;       % Damper function
  uECI                = QTForm(q,d.dampingData.u);
  hMag                = Dot(uECI,bI   )/mu0;
  hMagDot             = Dot(uECI,bIDot)/mu0;
  z                   = BFromHHysteresis( hMag, hMagDot, d.dampingData );
  %z = 0*z;
end

% Initial state vector
%---------------------
omega       = [0;-sqrt(d.mu/Mag(r)^3);0];
x0          = [r;v;q;omega;z']; % State vector

% Determine the simulation duration
%----------------------------------
orbPeriod   = Period(el(1));
d.end       = nOrbits*orbPeriod;

Run the simulation

%--------------------
disp('Beginning simulation...')
outf      = @(t,y,flag) HysteresisOutput(t,y,flag,d);
opts      = odeset('outputfcn',outf,'initialstep',2,'refine',4);
rhs       = @(t,x) RHSRigidBodyWithDamping( x, t, d);
[tout,y]	= ode45(rhs, [0 nOrbits*orbPeriod], x0, opts);
y = y';
disp('Done')

xP = HysteresisOutput([],[],'x');
zP = HysteresisOutput([],[],'z');
tP = HysteresisOutput([],[],'t');

Beginning simulation...
10% finished
20% finished
30% finished
40% finished
50% finished
60% finished
70% finished
80% finished
90% finished
Done

Plot