Write a simple code with good commenting practice:

Question:

At this stage in the design of Sentinel-6B, you are ready to consider the attitude determination and
control system. The spacecraft is planned to have a nominal altitude of 1336km, but you want to
investigate the major disturbance torques at a selection of different altitudes in case there is a
change in the mission later. The altitudes you select for consideration (after consulting the rest of
the design team) are 1000km, 1336km, and 2000km. You decide to make the assumption that
throughout the spacecraft’s orbit, it will be under the influence of disturbance torques from two
main sources: gravity (TG), and aerodynamic drag (TA). For this initial estimation you can assume
that the satellite has a mass moment of inertia along the Y axis (IYY) of 25 kg.m−2 , a mass moment
of inertia along the Z axis (IZZ) of 40 kg.m−2 . You can also assume a coefficient of drag (Cd) of 2.2.
The cross-sectional area should be consistent with the assumptions made in the thermal tutorial and
you should make a reasonable assumption for the distance from the centre of mass to the centre of
gravity. The air density (ρ) at the different altitude can be assumed to be 3 × 10−15kg.m−3 at
1000km, 8.5 × 10−16kg.m−3 at 1336km and 7 × 10−17kg.m−3 at 2000km.

Note that you will need to calculate the orbital velocity yourself for the different altitudes. Your goal
is to investigate the maximum torques and momentum storage requirements of the spacecraft at
each of these altitudes and, based on this, suggest a suitable attitude control subsystem. Tasks To
investigate your attitude control requirements you will need to calculate the disturbance torques the
spacecraft experiences over one orbit due to gravity and drag. For the purposes of this exercise, we
will assume that the torque due to the gravity gradient will vary sinusoidally over one orbital period.
For simplicity we will also assume that the disturbance torques due to aerodynamic drag is constant
over the orbit.

Your task is to calculate the magnitude of the torques that will be acting on the satellite at each
point throughout its orbit. You will then need to integrate these over time to get the momentum
storage (H) that is required over one orbit; this has the SI unit Nms (make sure you don’t confuse H,
which here represents the required momentum storage, with h, which is used to represent altitude).
Finally, you will need to plot a graph of torque vs. time, with the area that you integrated over
clearly displayed. Note that when calculating the torques you will need to save the torque
calculations as a function, and call the function in the main script. While it is possible to complete
this exercise without using functions, this exercise is based on testing your ability to use functions.
As such failing to properly utilise them will affect your final marks. In reality, the ability to create and
use functions is a valuable one, that can save you a lot of time when coding. Calling a function
multiple times in a script is a lot easier than repetitively typing out the same section of code.
Functions also make it easier to go back and upgrade the code at a later date. This task will require
you to calculate algebraic equations, perform an integration, utilise functions and plot your results in
a graph.

Sentinel-6B is a cube satellite with total surface area (A_surf) = 20.1 m^2 with mass of 1000 kg.

1. This first step is to define the various inputs and physical constants that you will need in your
code. To make things easier, a list of the parameters that you will need to create is provided below.
Note, that now you have three different altitudes. Hence, you will have to create an array with those
altitudes. Property given in the following format

Symbol – Property – Unit – value (value stated for only some properties, rest can be found above –
or have to be calculated)
 R_e - Radius of Earth – meters [m] – 6371 km
 Mu - Standard gravitational parameter of Earth - m3 s −2 – 386400 m^3* s^ −2
 hTorque - Array of altitudes - m
 IZZ - Mass moment of inertia about the z-axis - (kgm^−2)- 40
 IYY - Mass moment of inertia about the y-axis – (kgm^−2) - 25
 A_cross - Cross sectional area of the satellite - m2 – 3.35 m^2
 Delta_d - Distance from centre of mass to centre of aerodynamic pressure - m
 Rho_ADCS - Array of density – kgm^−3
 Cd - Coefficient of drag – no units. – 2.2

2. Before you calculate the torque due to gravity, you will have to create a for loop that will run
as many times as the number of elements in the altitude array, as the gravity torque depends
on the radius of the orbit. [HINT: Check what the ‘length’ function does using the help provided
by MATLAB. This way, you will be able to run your code no matter how many elements you have
in the altitude array.]

3. The next step is to calculate the maximum torque due to the gravity gradient. As previously
discussed, you will need to do this using functions. We will start with the maximum torque due
to gravity gradient.

OrbitalRadius (R_orb) = R_e + hTorque(i)

theta = pi/2; % Theta value in radians

gravityTorque = (3 * Mu / (2 *R_orb^3)) * abs(IZZ - IYY) * sin(theta);

4. There are multiple ways to generate functions, but since we want the final output to be a
single Live Script, we will generate a function within the script. These must go at the very end of
the script or they will not work. The first line of a function should be in the format:

‘function [outputArg1,outputArg2] = untitled(inputArg1,inputArg2)’

Here ‘outputArg1,outputArg2’ represent two placeholder outputs, untitled represents a

placeholder function name, and ‘inputArg1, inputArg2’ represent two placeholder inputs. We
want a function that will take in the IZZ, IY Y , Radius of orbit and µ values, and output a TG
value. With this in mind, you need to replace the inputs with the variables you created for the
Mass moments of inertia, the radius of orbit and the standard gravitational parameter of Earth,
and replace the output with a variable name for the torque due to gravity. Also make sure to
change the function name to something that clearly explains what it does.

5. As the radius of orbit is one of the inputs of the function and because it changes every time
the loop runs, you will have to calculate the radius of orbit before calling the function. Don’t
forget to identify which of the several elements of the altitude array you are using for each run.

6. Within the actual function, you will need to include the algebraic formula. The formula will
use the IZZ, IYY , Radius of orbit and µ values as inputs, but you will still need to create variables
for the constants. Within the function, you will need to define two variables, one for the radius
of Earth, and one for µ, Earth’s gravitational constant. Once these are defined, you can input the
equation for the gravitation torque (TG)

(gravityTorque = (3 * Mu / (2 * R_orb^3)) * abs(IZZ - IYY) * sin(theta);).

You will also need to assume a value for the roll angle of the spacecraft θ=pi/2; since you are
doing an initial sizing it would be sensible to consider a worst case value for this. It is important
that the variable name you give for TG is the same as the name given as the output in the first
line of the function. Finally, make sure that the last line of the function is the word end. This will
signal to MATLAB that this is a complete function that it can call.

7. You should now be familiar enough with creating functions to generate some for the torque
due to drag. You will do this using equation

V = sqrt(Mu / orbitalRadius)

TD = ½* Rho_ADCS*Cd*A_cross*V^2*(delta_d) – equation for torque due to drag

8. Now the functions have been created you will need to return to the main script. You will need
to use these functions to calculate the values for TG and TD. This will involve calling the
functions like you would any other MATLAB function, using the format
‘outputname=functionname(inputnames)’. You should now have values for all of the torque
values. As the functions were called inside a for loop, don’t forget add the counter of the loop.

9. Next you will need to calculate the orbital period (τ ). You should be familiar with this
equation by now which is:

tau = 2*pi*sqrt(R_orb^3/mu)

10. Now you will need an equation to calculate the total torque over time. As previously
mentioned, the TG value will vary over time, with the value you calculated previously
representing the maximum value. You will need to take this into account when coding your
equation. You should create a temporary variable (we’ll refer to it as t), and create a function of
this variable that reflects the change in torque over time. You will need to use the following
equation.

f(t) = cos [(t/τ)* 2π]* TG + TD

There is a specific format for inputting equations like this into MATLAB. For example, if you were
to include the formula y=t+1, you could input this to MATLAB as ‘y=@(t) t+1’ Note that you
would not need to have defined t before this line, as the ‘@(t)’ section of the code counts as
defining the variable.

Using this format, you will need to input the above equation, and save f(t) under a suitable
name.

11. Next you will need to calculate the momentum storage that is required over one orbit ‘(H)’.
You will need to do this every time the loop runs. To get this, you will need to evaluate the
numerical integral of f(t) over the orbital period. To do this, you can use the ‘integral’ function.
Once this is done, set your code to output the value of the angular momentum over one orbit
‘H’. There are a variety of functions that you can use for this, such as ‘disp’ or ‘fprintf’, but make
sure that your output is understandable and user friendly.

12. Finally, you will need to plot the results to a graph every time the loop runs. Try using
‘figure(counter)’. This will create as many figures as the maximum value of your counter. The
function f(t) represents the y-axis values, but you will need to create a vector that will represent
the x axis values. Using the easy method of vector creation that we have previously 4 covered
‘(minvalue:stepsize:maxvalue)’, create a vector from zero to the orbital period. It’s up to you to
pick a step size, you should try and make it as large as possible, to save processing time, but not
so large that it skews the results. We will refer to this vector as x.

13. Ordinarily, plotting x against f(x) would require the ‘plot’ function. However, as we would
like you to indicate the area that we have integrated over, we advise you to use the ‘area’
function. Once this is done, label the x and y axes of the graph.