MASTER THESIS
Magne ic Sys ems o he ADCS o a
Fem osa elli e
Manuel Cachón Vigil
SUPERVISED BY
Jo di Gu ié ez Cabello
Uni e si a Poli ècnica de Ca alunya
Mas e in Ae ospace Science & Technology
Feb ua y 2023
Magne ic sys ems o he ADCS o a
em osa elli e
BY
Manuel Cachón Vigil
DIPLOMA THESIS FOR DEGREE
Mas e in Ae ospace Science and Technology
AT
Uni e si a Poli ècnica de Ca alunya
SUPERVISED BY:
Jo di Gu ié ez Cabello
Depa men o Physics
ABSTRACT
This p ojec aims o selec he mos sui able ADCS o a em osa elli e in ended o
measu e he a mosphe e a low o bi al al i udes (unde 500 km); o do so, we will
employ a MATLAB code ha models he sys ems pe o mance in e ms o o al
o a ion speed, which is a ge ed o be less han a deg ee pe second a e a week.
To ha end, wo sys ems a e s udied and compa ed. Those a e hys e esis ma e ial
and a se o magne o que s.
In o de o model hei ope a ion, he o bi is p opaga ed using wo models (SGP4
and SGP8). The s a ing posi ion o hem is gi en by means o a TLE, and i is he
same o all o hem. Then, wi h he p opaga ed new loca ions o he sa elli e, he
local magne ic ield o each one o hem is calcula ed by means o wo widely-used
nume ical models (WMM and IGRF). This gene a es ou da a se s o he magne ic
ield ob ained wi h ou sligh ly di e en modelling echniques ha p o ide a way o
ge eliable esul s.
These ou da a se s a e coupled wi h a simula ion o he sa elli e's a i ude ha s a s
wi h a a he la ge ini ial o a ion speed (20 deg ees pe second in each axis o he
body e e ence ame; hese condi ions a e he s anda d wo s case o sa elli e
ejec ion om he launche ). Tha is modi ied by he ac ua ion o one o he p oposed
a i ude sys ems on he sa elli e as a unc ion o he local magne ic ield.
A pe manen magne would ne e s abilise he o a ion, and so i is dismissed as no
use ul. The hys e esis ma e ial cancels he o a ion in less han5 hou s, educing he
ime ha akes o achie e s abili y he lowe he o bi is. Rega ding he
magne o que s, hey achie e s abili y in unde 1 hou , and in mos o he cases in
unde 30 minu es.
The da a ob ained clea ly shows ha he as es de umbling me hod a e he
magne o que s as he s udy o a COTS no op imised model achie es a s abilisa ion
ime is an o de o magni ude as e han he nex bes me hod, is complian wi h he
equi emen o being less han 200 g and ha ing a peak powe consump ion o less
han 250mW, and a e de umbling less han 10 mW.

To my amily and Nadia.
Table o Con en s
CHAPTER 1 INTRODUCTION ................................................................................... 1
1.1.
A mosphe e ................................................................................................................................... 1
1.2.
How o measu e he a mosphe e ................................................................................................ 3
1.3.
ADCS .............................................................................................................................................. 3
1.4.
Modelling ools used .................................................................................................................... 5
1.4.1.
Two-Line Elemen s (TLE) .................................................................................................. 5
1.4.2.
Simpli ied pe u ba ion models .......................................................................................... 5
1.4.3.
MATLAB............................................................................................................................. 5
1.4.4.
Qua e nions ....................................................................................................................... 6
CHAPTER 2 PRELIMINARY WORK, ASSUMPTIONS AND REQUIREMENTS ....... 9
2.1.
Scope o he p ojec ..................................................................................................................... 9
2.1.1.
Ou o scope opics ............................................................................................................ 9
2.2.
Requi emen s ..............................................................................................................................10
2.3.
Assump ions ...............................................................................................................................11
CHAPTER 3 ORBIT AND ATTITUDE MODELLING ............................................... 13
3.1.
Code gene al desc ip ion ...........................................................................................................13
3.2.
Inpu s ...........................................................................................................................................14
3.3.
Ini ial pa ame e s ........................................................................................................................14
3.4.
Common Loop ............................................................................................................................15
3.5.
Hys e esis Me hod ......................................................................................................................15
3.6.
Magne o que s ............................................................................................................................16
3.7.
End o he loop ............................................................................................................................17
3.8.
Plo ...............................................................................................................................................17
CHAPTER 4 POTENTIAL ATTITUDE SYSTEMS ................................................... 19
4.1.
Hys e esis Rods ..........................................................................................................................19
4.1.1.
Desc ip ion o he sys em ................................................................................................19
4.1.2.
Simula ion de ails .............................................................................................................20
4.2.
Magne o que s ............................................................................................................................21
4.2.1.
Desc ip ion o he sys em ................................................................................................21
4.2.2.
Simula ion de ails .............................................................................................................21
CHAPTER 5 CODE VERIFICATION ........................................................................ 23
5.1.
P opaga o e i ica ion ..............................................................................................................23
5.2.
Magne ic model e i ica ion ......................................................................................................26
5.3.
Ro a ion e i ica ion ...................................................................................................................28

In oduc ion
1.1. A mosphe e
The Spanish Language
Royal Academy de ines he a mosphe e as “ he gaseous
laye ha su ounds he Ea h and o he cel
s uc u e a e di e en
in each case, as he dis ance o he pa en s a , and size o he
body –
among o he ac o s
a mosphe e is composed mainly o ni ogen (78%
amoun s o o he gases like
Ne e heless,
in e e y celes ial body
p essu e,
and he e o e densi y
dec
ease o p essu e in an iso he mal a mosphe e has a loga i hmic ela ion wi h he
al i ude inc ease, ollowing
1.1
whe e 
is he molecula mass,
al i ude, 
is he ideal gas cons an and
heigh s.
This o mula can gi e us a g oss app oxima ion o he alues, bu
a mosphe es a e
no iso he mal, and in ac his a ia ion
allows us o s a i y i in o se e al laye s as
Ea h [1].
A s a ic s anda d model called In e na ional S anda d A mosphe e (ISA) o how
p essu e, emp
e a u e, densi y and iscosi y changes wi h al i ude
he 1970s
. I includes ables o alues and he o mulas ha app oxima e he
a mosphe ic p ope ies, and
condi ions. Ne e heless,
his model is only eliable o he lowe and dense laye s o
he a mosphe e, and p edic ing i s p ope
Figu e 1.1
A mosphe ic d ag posi i e eedback lo
Chap e 1
INTRODUCTION
Royal Academy de ines he a mosphe e as “ he gaseous
laye ha su ounds he Ea h and o he cel
es ial bodies”. I s composi ion and
in each case, as he dis ance o he pa en s a , and size o he
among o he ac o s
–
in luence i s e olu ion and beha iou . Fo he Ea h
a mosphe e is composed mainly o ni ogen (78%
) and oxygen (21%) wi h ace
amoun s o o he gases like
a gon o ca bon dioxide (CO
2
).
in e e y celes ial body
he a mosphe e
expe iences a dec ease in
and he e o e densi y
,
as a unc ion o he heigh om he su ace. This
ease o p essu e in an iso he mal a mosphe e has a loga i hmic ela ion wi h he
1.1
∆pe
∆

is he molecula mass,
 is he accele a ion o g a i y,ℎ
is he ideal gas cons an and
 is he a e age
empe a u e be ween bo h
This o mula can gi e us a g oss app oxima ion o he alues, bu
no iso he mal, and in ac his a ia ion
o empe a u e wi h heigh
allows us o s a i y i in o se e al laye s as
shown in Table 1.1
o he case o he
A s a ic s anda d model called In e na ional S anda d A mosphe e (ISA) o how
e a u e, densi y and iscosi y changes wi h al i ude
was
. I includes ables o alues and he o mulas ha app oxima e he
a mosphe ic p ope ies, and
i
is i al o a ia ion, as i hea ily elies on a mosphe ic
his model is only eliable o he lowe and dense laye s o
he a mosphe e, and p edic ing i s p ope
ies in space
is a beyond i s capabili ies.
A mosphe ic d ag posi i e eedback lo
op
. ME s ands o Mechanical
Ene gy.
1
Royal Academy de ines he a mosphe e as “ he gaseous
es ial bodies”. I s composi ion and
in each case, as he dis ance o he pa en s a , and size o he
in luence i s e olu ion and beha iou . Fo he Ea h
, he
) and oxygen (21%) wi h ace
expe iences a dec ease in
as a unc ion o he heigh om he su ace. This
ease o p essu e in an iso he mal a mosphe e has a loga i hmic ela ion wi h he
1.1
is he a ia ion in
empe a u e be ween bo h
This o mula can gi e us a g oss app oxima ion o he alues, bu
eal
o empe a u e wi h heigh
o he case o he
A s a ic s anda d model called In e na ional S anda d A mosphe e (ISA) o how
was
de eloped in
. I includes ables o alues and he o mulas ha app oxima e he
is i al o a ia ion, as i hea ily elies on a mosphe ic
his model is only eliable o he lowe and dense laye s o
is a beyond i s capabili ies.
. ME s ands o Mechanical
2 Magne ic sys ems o he ADCS o a Fem osa elli e
We conside ha space begins a he on Ká mán line, es ablished a 100 km by he
Fédé a ion Aé onau ique In e na ionale (FAI), and mos o he egula o y agencies
like he UN accep his alue o some hing close o i as he bounda y o space. And
al hough mos sa elli es' o bi s will always emain well abo e his bounda y, he e
emains s ill enough gas in LEO o a ec hei dynamical beha iou . This means ha
any spacec a inside he he mosphe e ( he laye be ween 80 km and abou 500-600
km abo e sea le el) will expe ience equen collisions wi h pa icles inducing a
posi i e eedback loop (Figu e 1.1) ha ul ima ely causes he spacec a o bu n on
e-en y.
Table 1.1 Laye s o he a mosphe e
Laye Heigh s Cha ac e is ics
T oposphe e 0-12 km
Dec ease in empe a u e wi h heigh d i en by
su ace hea ing.
80% o he a mosphe e mass.
Mos o he a mosphe ic wea he occu s he e.
Only laye a ailable o p opelle d i en ai c a .
S a osphe e 12-50 km
Inc ease in empe a u e wi h heigh d i en by UV
adia ion.
S able condi ions, lack o u bulence induce
wea he .
Con ains he ozone laye .
Highes laye accessible by je powe ed ai c a .
Mesosphe e 50-80 km
Dec ease in empe a u e wi h heigh .
Coldes place on Ea h –85ºC
Mos me eo s bu n up he e.
Only accessible by sounding ocke s and ocke
p opelled ai c a .
The mosphe e
80-700 km
Inc ease in empe a u e wi h heigh .
Heigh a ies conside ably due o sola ac i i y.
Comple ely cloudless.
Con ains mos o he ionosphe e (whe e au o as
a e p oduced).
Molecules a el 1 km on a e age be ween
collisions wi h o he molecules.
Exosphe e 700-10,000 km
Mos spacec a o bi he e.
Doesn’ beha e like a collisional gas.
A oms so a apa ha can a el hund eds o km
be o e colliding wi h o he ones.
The e o e, he d ag caused by he emaining a mosphe ic gases has implica ions all
h oughou he mission as, o example, by limi ing he ope a ional li e o es ic ing
he o bi s a ailable o ensu e ha his decay akes longe han he expec ed
sa elli e’s ope a ional li e. I can also impac he a chi ec u e o he spacec a ,
demanding some kind o p opulsion o mi iga e his e ec , and hence equi ing o
de o e some o he limi ed mass budge o p opellan .
In oduc ion 3
Fu he mo e, he a mosphe ic pa ame e s a hose al i udes a y wi h sola cycles,
geomagne ic ac i i y, Ea h’s posi ion a ound he Sun, and Moon’s posi ion a ound
he Ea h, among o he ac o s. So, measu ing i s p ope ies is e y impo an o
accu a ely design space missions
1.2. How o measu e he a mosphe e
The easies way o measu e he a mosphe e a ha al i ude is by using a sa elli e o
known cha ac e is ics, and obse e how he o bi ge s modi ied by he a mosphe ic
e ec s in compa ison wi h he p edic ed o bi wi hou his d ag o ce. The e o e, his
p ojec is ocused on a mission ha will y o accomplish ha wi h he ollowing base
cha ac e is ics:
• Sphe ical em osa elli e
• No mo e han 400g
• 100 mm in diame e
• Low and as decaying o bi (ini ial al i udes in he ange 200 – 500 km)
• Nea diagonal ine ia enso .
To measu e he a mosphe ic densi y, we plan o use a e y sensible and p ecise
accele ome e o di ec ly de e mine he d ag o ce expe ienced by he sa elli e. To
p o ide accu a e posi ion and ime in o ma ion o he accele ome e ’s measu emen s,
we will employ a GNSS ecei e (GPS, GLONASS, and/o GALILEO). Fu he mo e,
hese loca ion and ime da a will p o ide a backup me hod o de e mine along-o bi
a e age densi ies i he al i ude is oo la ge o he accu acy o he accele ome e , o
in case o i s mal unc ion.
Ne e heless, he GNSS an enna has a limi a ion in o a ion a e o ensu e ha i is
able o always ecei e he signal o a leas 4 sa elli es and ind i s posi ion.
The e o e, he em osa elli e equi es some ADCS sys em in o de o educe he
angula speed o a manageable one.
1.3. ADCS
The e a e se e al me hods ha can be used o con ol he a i ude, each one wi h i s
p os and cons o ou sys em. A gene al desc ip ion o hem is shown inTable 1.2, as
well as a sho discussion abou hei sui abili y o ou p ojec .
The analysis shows ha mos o he lis ed ac ua o s ypes a e no sui able o he
mission, lea ing he magne ic me hods as he only iable op ions o be pu sued. O
hose, we can di e en ia e h ee ypes: pe manen magne s, hys e esis ods and
magne o que s. Thei gene al desc ip ion as well as hei p os and cons a e shown in
Table 1.3.¡E o ! No se encuen a el o igen de la e e encia.
4 Magne ic sys ems o he ADCS o a Fem osa elli e
Table 1.2 ADCS Me hods.
Me hod Desc ip ion Analysis
Th us e s
P opulsion sys em ha by being
i ed impa s a o que in a desi ed
axis, he e o e con olling he
a i ude.
Requi es p opellan which adds
mass o he spacec a and can
po en ially dis up he o bi .
Momen um
wheels
Wheels ha spin up o down o
change he angula momen um,
and he e o e con ol he a i ude
along he o a ional axis.
Mass and powe budge
cons ain s will make his
un easible.
Con ol
momen
gy os
Simila o he momen um wheels,
bu allow also change he axis in
which hey o a e, con olling in
his way he angula momen um
o he spacec a .
Simila o momen um wheels bu
adds e en mo e complexi y, hey
a e sui able o bigge spacec a .
The mos p ecise a i ude con ol
me hod in exis ence.
Sola
p essu e
Sola ligh impa s a small o ce
on he spacec a ; by con olling
he su ace exposed o ligh he
sa elli e can be s abilised.
Is no a signi ican o ce in he
LEO en i onmen in which he
sa elli e will be loca ed. On he
o he hand, he sa elli e's
sphe ical shape will nega e his
e ec .
Ae odynamic
A mosphe ic gas impa s a
p essu e on he spacec a ; by
con olling he c oss-sec ion o
he sa elli e, i can be s abilised
simila ly o an ai c a .
The sa elli e's sphe ical shape
and homogeneous mass
dis ibu ion will nega e his e ec .
Magne ic
ac ua o s
The sa elli e can gene a e a
magne ic ield o in e ac wi h he
Ea h's one, he e o e gene a ing
a o que ha can be used o
ADCS.
Ligh weigh , low powe
consump ion and ideal o low
o bi s.
Table 1.3 Types o magne ic ADCS.
Me hod Desc ip ion Analysis
Pe manen
magne s
The pe manen magne o ces
he sa elli e o align i sel wi h
he Ea h's magne ic ield.
The e is li le o no damping
e ec and he oscilla ion
gene a ed is simila o a
pendulum.
Hys e esis ods
A ha d e omagne ic ma e ial
ha can dissipa e ene gy while
o a ing h ough he Ea h's
magne ic ield.
Reduces oscilla ions o he
sa elli e bu is unable o selec
he o ien a ion; usually, i is
combined wi h a pe manen
magne o se said p e e ed
o ien a ion.
Magne o que s
A cu en un h ough me allic
coils gene a es a magne ic ield
ha in e ac s wi h he Ea h's,
adjus ing he posi ion.
P ecise and ac i e o ien a ion a
he expense o powe consump-
ion, mass and complexi y
added o he sa elli e.
In oduc ion 5
As seen in Table 1.3 pe manen magne s gi e no s abiliza ion, ins ead hey o ce he
spacec a o oscilla e a ound he magne ic ield lines, he e o e i alone is no
sui able o ou pu poses. Wi h ega d o he hys e esis ods and magne o que s, i s
main di e ence is ha he i s one is a e y simple passi e sys em ha does no
equi e any powe , while he second is ac i e bu a he cos o added complexi y and
powe equi ed.
Hys e esis ods equi e o be as slende as hey can o ge he mos damping
possible, bu he pa ame e ha mainly con ols hei e ec is he olume o he
ma e ial used.
The magne o que s eng h is di ec ly ela ed wi h he numbe o loops p esen in he
loop and how big hose a e, bu he longe he cable is, he bigge he esis ance, and
he lowe i s e ec .
1.4. Modelling ools used
In o de o pe o m he simula ion o he sa elli e’s a i ude, we equi e a se o ools
ha we desc ibe in he nex sec ions.
1.4.1. Two-Line Elemen s (TLE)
Two-Line Elemen Se is a o ma o encoding and dis ibu ing o bi al pa ame e s o
Ea h o bi ing objec s de eloped by NORAD in he 1960s and 70s. I s o ma o wo
lines o 69 cha ac e s is due o i s o iginal pu pose: o be used wi h punched ca ds.
Fo many yea s, USAF and la e USSF ha e kep ack and assigned TLEs o all
de ec able o bi ing objec s a ound he Ea h, and mos o hem a e published online
making i a de ac o s anda d o his ype o in o ma ion.
1.4.2. Simpli ied pe u ba ion models
Those a e a se o i e ma hema ical models used o p opaga e he o bi al s a e
ec o s accoun ing o he main e ec s o pe u ba ions caused by he Ea h and
some phenomena ound in he LEO en i onmen . O hose nume ical models, he
Simpli ied Gene al Pe u ba ions (SGP) models a e he ones used o nea Ea h
o bi s. In his s udy we will use he SGP4 [2] and i s e ised and imp o ed e sion
SGP8 [3].
The p incipal ad an age o hese models is a signi ican educ ion in he
compu a ional bu den associa ed wi h o bi p opaga ion.
1.4.3. MATLAB
Abb e ia ion o "MAT ix LABo a o y", is a p og amming and nume ical compu ing
pla o m widely used by enginee s and scien is s o analyse da a, de elop algo i hms,
and c ea e models[4].Ini ially eleased in 1984, i has been con inually upda ed and
e ined, and i s capabili ies can be enla ged and imp o ed by using Toolboxes and
Simulink, i s g aphical p og amming en i onmen .

6 Magne ic sys ems o he ADCS o a Fem osa elli e
1.4.3.1. Ae ospace and mapping oolboxes
These se s o p e-p og ammed unc ions a e no included in he basic MATLAB
e sion. Ae ospace Toolbox p o ides s anda ds-based ools and unc ions o
analysing he mo ion, mission, and en i onmen o ae ospace ehicles [5] while
Mapping Toolbox™ p o ides algo i hms and unc ions o ans o ming geog aphic
da a and c ea ing map displays [6].
1.4.3.2. Vec o T ans o ma ions
In o de o pe o m he ans o ma ion be ween ce ain coo dina e ames, bo h
oolboxes ha e speci ic commands al eady implemen ed ha help us in making
hose calcula ions. Those will be:
• "ned2ece " om he mapping oolbox, ha allows us o change ec o s om
No h-Eas -Down e e ence ame o Ea h Cen ed Ea h Fixed.
• "ece 2eci" om he ae ospace oolbox, allows us o ans o m he esul
om he p e ious command in o Ea h-Cen ed-Ine ial ame o e e ence, he
same ha SGP models p o ide.
• "eci2lla" also om he ae ospace oolbox will change ECI ec o s in o
La i ude-Longi ude-Al i ude equi ed o use he Magne ic models also included
inside he oolbox.
1.4.3.3. In e na ional Geomagne ic Re e ence Field (IGRF)
IGRF is a ma hema ical model ha desc ibes Ea h's magne ic ield and i s secula
a ia ion p oduced by IAGA since 1965 and upda ed e e y 5 yea s [7]. The cu en
las upda ed e sion, and he one used o his p ojec , is he 13 h (IGRF-13)
eleased in Decembe 2019. I is included in MATLAB’s Ae ospace oolbox.
1.4.3.4. Wo ld Magne ic Model (WMM)
The WMM is also a ma hema ical model ha desc ibes Ea h's magne ic ield and i s
secula a ia ion, bu in his case p oduced by USA and UK go e nmen agencies
and i is he s anda d o NATO [8]. I is also upda ed e e y 5 yea s and is included in
MATLAB’s Ae ospace oolbox. The las a ailable elease is he WMM-2020
published also in Decembe 2019.
1.4.4. Qua e nions
Qua e nions we e i s ly desc ibed by I ish ma hema ician William Rowan Hamil on in
1843. They a e ma hema ical objec s ha ex end he eal numbe sys em simila o
he complex numbe s (in ac , qua e nions a e hype complex numbe s) [9]. They
consis o wo pa s, a scala and a ec o (also called he imagina y pa ). They a e
o en ep esen ed in he o ms shown in 1.2
 



+ 

 + 

 + 

 1.2
In oduc ion 7
whe e 

, 

, 

and 

a e eal numbe s and ,  and  he imagina y uni s ha
sa is y equa ions 1.3






−1 1.3
Qua e nions ha e been used in a a ie y o ields such as physics, enginee ing, and
compu e science. In physics and enginee ing, qua e nions a e used o ep esen
o a ions in h ee-dimensional space, and as a esul , hey a e o en used in
ae ospace, obo ics and compu e g aphics. This o a ion can be pe o med by
ollowing 1.4.
!
 
∗
1.4
whe e is he ec o ha is o be o a ed (w i en as a qua e nion wi h 0 as scala
pa and as ec o pa ),
!
is he o a ed ec o ,  is he qua e nion and 
∗
is he
conjuga e o said qua e nion. The conjuga e o he qua e nion in exp ession 1.2 is
1.5.
#
∗
≡%#
&
−#
'
'

(#
&
−#
)
*−#
+
,−#
-
. 1.5
Then, he o a ion o equa ion 1.4 can be exp essed as show by 1.6.
!
/0, +2

3∧ 5+2∧3∧ 56 1.6
Addi ionally, he qua e nion used o o a e a ound a known axis 7 an angle 8 is gi en
by 1.7

⎝
⎜
⎜
⎛
cos?82
@A
7

BC?82
@A
7

BC?82
@A
7

BC?82
@A
⎠
⎟
⎟
⎞
1.7
This las exp ession is he usual ie be ween he axis-angle and he qua e nion
ep esen a ions o o a ions.
8 Magne ic sys ems o he ADCS o a Fem osa elli e
P elimina y wo k assump ions and equi emen s 9
Chap e 2
P elimina y wo k, assump ions and equi emen s
2.1. Scope o he p ojec
This p ojec a emp s o model he pe o mance o he ADCS o a em osa elli e
whose goal is o measu e he a mosphe ic densi y in e y low Ea h o bi (VLEO,
hus wi h al i udes below 250 km). The a i ude con ol will enable con inuous
connec ion wi h GNSS sys ems by educing he em osa elli e’s o a ion a e and
he e o e allowing densi y measu emen s o ha e accu a e ime s amps and
geog aphical loca ion in o ma ion.
To his end, wo ypes o magne ic ac ua o s will be independen ly s udied: hys e esis
ods and magne o que s. The o bi will be modelled om a syn he ic TLE gene a ed
acco ding o he equi emen s exp essed in sec ion 2.2. In o de o accomplish his
ask wi hou an eno mous compu a ional bu den, wo SGP p opaga o models will be
used: SGP4 and SGP8. The o bi al da a p o ided by hose will be used o calcula e
he local magne ic ield using wo models, WMM and IGRF. The local magne ic ield
will in e ac wi h he ac ua o s inside he sa elli e, he e o e modi ying he o a ion a e
o he spacec a .
This wo k will be pe o med in MATLAB, and he code will be explained in Chap e 3
by means o a low diag am and a s ep by s ep explana ion o he p ocesses
in ol ed. We will pa icula ise ou analysis o each ype o sys em in Chap e 4.
Then a e i ica ion p ocess will be unde aken using educed sec ions o he code, o
assess he di e ence be ween he wo p opaga o s in sec ion 5.1, he di e ence
be ween magne ic models in sec ion 5.2, and he magne ic ield o a ions will be
checked in sec ion 5.3. Chap e 6 will be de o ed o analyse he esul ing o a ion
speeds wi h each one o he me hods o h ee di e en o bi s di e ing on hei
al i ude. In sec ion 6.3, he me hod ha achie es he as es s abilisa ion ime will be
selec ed o be he inal one, p o ided ha all he pa ame e s ela ed o i a e
complian wi h he ones gi en by he equi emen s; i his is no he case, he nex
as es me hod will be selec ed.
Finally, we d aw ou conclusions and iden i y possible alleys o u u e esea ch and
imp o emen s o he sys em.
2.1.1. Ou o scope opics
The code ha is he main pa o his p ojec will no include a g aphic in e ace as a
way o inpu he in o ma ion equi ed o he simula ion; hose da a will be yped in
he sec ion alloca ed o ha pu pose inside he code and acco dingly labelled.
16 Magne ic sys ems o he ADCS o a Fem osa elli e
• Upda e angula eloci y: wi h ha a ia ion o kine ic ene gy, he new one can
be calcula ed and he e o e he new angula eloci y can also be calcula ed
applying 3.1 backwa ds as shown in 3.3, whe e ] is he axis o o a ion (a uni
ec o ) and^
P
he momen o ine ia calcula ed a ound ha axis.
H
'
'
'

Z
'
'

_
+3G.`∆G.5
JZ
3.3
• Ve i y s abiliza ion: check i angula eloci y is below he maximum accep able
alue gi en in he equi emen , and sa e he ime a which his happens. I all
ou models ha e done so, s op he loop om i e a ing any u he .
3.6. Magne o que s
In he case in which he ADCS is based on a magne o que , he loop will pe o m he
ollowing s eps:
• Find μ: as he magne ic ield in BRF changes om one i e a ion o he nex , o
ind his alue i s ly we will calcula e he a ge accele a ion in o de o s op
he o a ion in a gi en ∆a (3.4):
b =
`c
'
'
'

∆R
3.4
Then we will calcula e he o que equi ed wi h 3.5.

'

= ^ b 3.5
And hen he dipola momen is gi en by3.6.
d'
'

=‖d'
'

‖d 3.6
whe e d is he uni ec o in he di ec ion in which d'
'

is poin ed a , and ‖d'
'

‖ is
gi en by equa ion 3.7.
‖g‖=
Xh
'

×j
'

X
Xk
'
'

Xjl
3.7
m
'
'

is gi en by3.8
m
'
'

=
?n
∙j
'

Aj
'

jl
− gp 3.8
• Calcula e he Magne ic To que: by doing he c oss p oduc be ween he
dipola momen o he sa elli e and he local magne ic ield (3.9)

O bi and a i ude modelling 17
I
'
'

*
d'
'

*
×Y
'
'

*
3.9
• Calcula e angula accele a ion: The angula accele a ion can be ound sol ing
he linea sys em gene a ed by he ine ia enso and he o que (3.10).
I
'
'

*
Jq
'
'

*
3.10
• Calcula e angle u ned: using he equa ion o he uni o mly accele a ed
ci cula mo ion, he angula accele a ion and he angula o a ion a he
beginning o he i e a ion (3.11).
'
'

*s)
H
'
'
'

*
∆ +
)
+
q
'
'

*
∆
+
3.11
• Ro a e sa elli e main axis: using a qua e nion o ep esen he o a ion
depic ed by he p e ious s ep, being he axis he o a ion ec o in uni o m.
• Upda e angula eloci y: using he uni o mly accele a ed ci cula mo ion, he
p e ious angula o a ion, he angula accele a ion and he ime inc emen
(3.12).
H
'
'
'

*s)
= H
'
'
'

*
+ q
'
'

*
∆ 3.12
3.7. End o he loop
• Sa e da a: s o e he o a ion speed and ime da a in o a ays ha can be used
a e wa ds ou side he loop.
• Upda e he coun e s: add one o he cons an s ha i e a e he di e en
elemen s inside he loop o s a he new i e a ion. I also p in s on he sc een
he cu en pe cen age o he o al i e a ions scheduled al eady pe o med.
3.8. Plo
The da a s o ed in he a ays is plo ed o show he pe o mance o he me hod; his
includes angula speeds in he h ee axes, as well as he module o he angula
speed ec o o each one o he 4 scena ios. I also p in s he ime when s abiliza ion
has occu ed in he hys e esis me hod.
18 Magne ic sys ems o he ADCS o a Fem osa elli e
Po en ial A i ude Sys ems 19
Chap e 4
Po en ial A i ude Sys ems
4.1. Hys e esis Rods
4.1.1. Desc ip ion o he sys em
When a magne ic ield om an ex e nal sou ce is applied o a e omagne ic ma e ial,
i s a omic s uc u e ( he a omic magne ic dipoles) aligns wi h i , pa o i emaining in
ha aligned s a e a e he ield i s emo ed. To comple ely demagne ize he ma e ial,
hea (causing a omic ib a ions ha des oy he alignmen ) o a magne ic ield in he
opposi e di ec ion is equi ed.
Because he ela ionship be ween ield s eng h H and magne iza ion m is no linea
in such ma e ials (Figu e 4.1) and causes some in e nal ic ion, i can be used o
slow he o a ion o he sa elli e ollowing equa ion 3.2.
Figu e 4.1 Theo e ical model o magne iza ion agains magne ic ield. Figu e aken
om [11].
In ou case, he hys e esis ma e ial is aken o be AlNiCo wi h a magne ic coe ci i y
H
c
o 17.27 A/m and a densi y o 8.25 g/cm
3
; we will assume ha i is dis ibu ed in 2
ods o 5 cm in leng h, p o iding a o al olume o 5 cm
3
, he e o e amoun ing o a
o al mass o 41.25g.
20
I should be no ed ha he sa u a ion magne ic
much la ge han he
local
sa elli e ope a es. Due o he lack o eliable da a o model he eal pa ial
magne isa ion cu e,
i is assumed o be linea
whe e K
'
'

L
and K
'
'

a e he coe ci i y in he cases o sa u a ing magne ic induc ion (
and e es ial magne ic induc ion (
4.1.2. Simula ion de ails
In his case he o a ional speed is upda ed by calcula ing he kine ic ene gy o
o a ion o he sa elli e (
3.1
checked i i s new kine ic
ene gy
new o a ion a e is calcula ed wi h equa ion
nega i e
(a physically un ealis ic si ua ion)
he da e is s o ed, and once he ou models a e s
Figu e 4.2
Loop pe o med by he use o hys e esis ma e ial
Magne ic sys ems o he ADCS o a Fem osa elli e
I should be no ed ha he sa u a ion magne ic
induc ion
o his ma e ial is 1.53 T,
local
e es ial magne
ic ield a he al i udes
sa elli e ope a es. Due o he lack o eliable da a o model he eal pa ial
i is assumed o be linea
u
'
'

j
'

w

u
'
'

x
j
'

y
a e he coe ci i y in he cases o sa u a ing magne ic induc ion (
and e es ial magne ic induc ion (
O
'

z
).
In his case he o a ional speed is upda ed by calcula ing he kine ic ene gy o
3.1
), hen he ene gy dissipa ed in one cycle (
ene gy
is posi i e ( ha is, whe he {
|
−
new o a ion a e is calcula ed wi h equa ion
3.3. I
he o a ional kine ic ene gy
(a physically un ealis ic si ua ion)
, he
sa elli e is conside ed
he da e is s o ed, and once he ou models a e s
able he simula ion is s opped.
Loop pe o med by he use o hys e esis ma e ial
Magne ic sys ems o he ADCS o a Fem osa elli e
o his ma e ial is 1.53 T,
ic ield a he al i udes
in which he
sa elli e ope a es. Due o he lack o eliable da a o model he eal pa ial
4.1
a e he coe ci i y in he cases o sa u a ing magne ic induc ion (
O
'

}
)
In his case he o a ional speed is upda ed by calcula ing he kine ic ene gy o
), hen he ene gy dissipa ed in one cycle (
3.2), hen i is
∆{
|
~0) and he
he o a ional kine ic ene gy
is
sa elli e is conside ed
o be s abilized,
able he simula ion is s opped.
Loop pe o med by he use o hys e esis ma e ial
.
Po en ial A i ude Sys ems 21
4.2. Magne o que s
4.2.1. Desc ip ion o he sys em
Magne o que s a e essen ially elec omagne s, made by looping a wi e mul iple imes
wi h a known sec ion, ha gene a e a magne ic dipola momen gi en by equa ion
4.2.
d'
'

• J €
'
'

4.2
The dipola magne ic momen depends on he numbe o loops (C) o he wi e, he
in ensi y o he cu en un h ough i (^), and he a ea inside he loops (•
R
); his a ea
is conside ed as a ec o , wi h a di ec ion pe pendicula o he su ace in he
s anda d way desc ibed in classical opology. This g is gene ally es ic ed (4.3) by
he maximum powe allowed o eed he ac ua o ‚
ƒQ„
.
d'
'

…[†
• €
'
'

_
‡…[†
ˆ‰*Š‹
4.3
whe e 
Œ•Ž•
is he esis ance o he cu en low gene a ed by he wi e. In ou s udy,
we will use a comme cial magne o que om New Space Sys ems (NSS), NCTR-
M003 [11] ha has a maximum powe consump ion o less han 250 mW, weigh o
unde 30 g and a size o 72mm×15mm×13mm, alues sui able o ou mission, as
h ee o hem will be equi ed. This sys em p o ides a maximum magne ic momen o
0.29 A m
2
.
Wi h hose pa ame e s and ollowing om 4.3 we can deduce ha he ela ion
be ween d'
'

and ‡ is gi en by 4.4
d•€
_
‡
ˆ‰*Š‹
√‡ 4.4
whe e  can be ound by ollowing 4.5

d
√‡

d…[†
‘‡…[†
4.5
and he powe consumed in e e y momen is gi en by 4.6
‡∑
d*
.


4.6
wi h he index being used o iden i y he di e en magne o que s.
4.2.2. Simula ion de ails
In his case, he dipola momen changes om one i e a ion o he nex . To de e mine
i , we will i s ly calcula e he a ge dipola momen (as explained in sec ion 3.6) ha

22
would be equi ed
in o de o s op
alue is unde he limi s o he sys em
momen
is o ced o be he maximum a ailable. Then
e olu ion is modelled.
Figu e 4.3
Loop pe o med by he use o magne o que s
Magne ic sys ems o he ADCS o a Fem osa elli e
in o de o s op
comple ely he o a ion. Then
, we will
alue is unde he limi s o he sys em
; i , on he con a y,
i is oo
is o ced o be he maximum a ailable. Then
,
wi h his alue he
Loop pe o med by he use o magne o que s
Magne ic sys ems o he ADCS o a Fem osa elli e
, we will
check i he
i is oo
la ge, he dipola
wi h his alue he
o a ion a e
Loop pe o med by he use o magne o que s
.
Code e i ica ion 23
Chap e 5
Code e i ica ion
A i al pa o he modelling p ocess is o e i y he pe o mance o he p ocesses
used. In his case as se e al me hods will be compa ed, i s ly i will be in e es ing o
compa e how bo h p opaga o s used di e ge o he same ime a ia ion; in his way,
we can ha e a hindsigh o which pa o a po en ial di e ence can be a ibu able o
he o bi p opaga o s. Secondly, i will be necessa y o compa e he wo magne ic
models o se e al posi ions a ound he o bi . Las ly, he o a ion o hese ields
should be moni o ed, in o de o e i y ha he o a ions a e no modi ying he alue,
jus changing he e e ence sys em in which hey a e ep esen ed.
5.1. P opaga o e i ica ion
In o de o compa e SGP4 and SGP8, he same eal li e TLE will be used o check
he di e ence in beha iou . Fo his, only a educed sec ion o he code will be
gene a ed (included in he annexes), p opaga ing he o bi o a o al o 20 days, and
s o ing a da a poin e e y 300s.
We will use he TLE om he I idium 65 sa elli e on 16/05/1999 (Figu e 5.1), which a
ha ime o bi ed a a heigh o a ound 785 km o e he Ea h's su ace.
1 25288U 98021D 99117.03750742 -.00000015 00000+0 -12318-4 0 1471
2 25288 86.3956 112.7366 0002518 67.8845 292.2618 14.34216487 55284
Figu e 5.1 I idium 65 TLE.
Plo ing he absolu e e o o each di ec ion, we can see ha he di e ence is ne e
la ge han 0.4%, and mos o he ime emains below 0.05% Figu e 5.2.
I we plo he e o o he no m o he ec o o he same se o da a, we can see
ha he e o always emains below 0.00018% (see Figu e 5.3) which shows ha in
p ac ice he di e ence is negligible.
24
Figu e 5.2
P opaga ed ec o componen s e
Figu e 5.3
P opaga ed ec o module e o I idium 65
Magne ic sys ems o he ADCS o a Fem osa elli e
P opaga ed ec o componen s e
o I idium 65
P opaga ed ec o module e o I idium 65
Magne ic sys ems o he ADCS o a Fem osa elli e
o I idium 65
.
P opaga ed ec o module e o I idium 65
.
Code e i ica ion
To check o possible al i ude dependen issues, we selec a
al i ude sa elli e Sa elio 1(
Fi
da a se co esponds o
18/01/2023.
1 47961U 21022AF
23018.66453903
2 47961
97.5171 280.0729 0019951 121.2971 239.0217 15.13768716 99338
Plo ing again he same bo h g aphs, we can see ha due o he a mosphe ic d ag
di e ence, he e is a s eady inc ease in he e o o m one o he o he (
and Figu e 5.6).
Figu e 5.5
P opaga ed ec o componen s e o I idium 65
This is no unexpec ed, as TLEs should no be p opaga ed o such long pe iods o
ime using SGP4 o SGP8. The inaccu acies in he TLEs, and he app oxima ions in
he e y na u e o SGPs will cause a apid inc ease in
why he NORAD publishes daily TLEs o all ac i e sa elli es and he la ges pieces o
space deb is.
To check o possible al i ude dependen issues, we selec a
TLE
se
Fi
gu e 5.4), hen loca ed
a a ound 540 km o heigh . The
18/01/2023.
23018.66453903
.00029170 00000+0 16080
-
97.5171 280.0729 0019951 121.2971 239.0217 15.13768716 99338
Figu e 5.4 Sa elio 1 TLE.
Plo ing again he same bo h g aphs, we can see ha due o he a mosphe ic d ag
di e ence, he e is a s eady inc ease in he e o o m one o he o he (
P opaga ed ec o componen s e o I idium 65
This is no unexpec ed, as TLEs should no be p opaga ed o such long pe iods o
ime using SGP4 o SGP8. The inaccu acies in he TLEs, and he app oxima ions in
he e y na u e o SGPs will cause a apid inc ease in
he e o . This is he eason
why he NORAD publishes daily TLEs o all ac i e sa elli es and he la ges pieces o
25
se
om he lowe
a a ound 540 km o heigh . The
-
2 0 9999
97.5171 280.0729 0019951 121.2971 239.0217 15.13768716 99338
Plo ing again he same bo h g aphs, we can see ha due o he a mosphe ic d ag
di e ence, he e is a s eady inc ease in he e o o m one o he o he (
Figu e 5.5
P opaga ed ec o componen s e o I idium 65
.
This is no unexpec ed, as TLEs should no be p opaga ed o such long pe iods o
ime using SGP4 o SGP8. The inaccu acies in he TLEs, and he app oxima ions in
he e o . This is he eason
why he NORAD publishes daily TLEs o all ac i e sa elli es and he la ges pieces o
32
Figu e 6.4
S abilisa ion imes a 300 km al i ude wi h hys e esis
Figu e 6.5
Ro a ion speed a 300 km al i ude wi h hys e esis
• 400 km:
his simula ion uns o
s abiliza ion imes
a e
shown in Figu e 6.
7
o he ou models.
Magne ic sys ems o he ADCS o a Fem osa elli e
S abilisa ion imes a 300 km al i ude wi h hys e esis
Ro a ion speed a 300 km al i ude wi h hys e esis
his simula ion uns o
0.2 days wi h a s ep o 2
second
a e
shown in Figu e 6.6,
and he o a ion a e e olu ion is
7
.
Again, he angula eloci y e olu ion is un
Magne ic sys ems o he ADCS o a Fem osa elli e
S abilisa ion imes a 300 km al i ude wi h hys e esis
.
Ro a ion speed a 300 km al i ude wi h hys e esis
.
second
s. Now he
and he o a ion a e e olu ion is
Again, he angula eloci y e olu ion is un
dis inguishable

Magne ic ADCS Pe o mance
Simula ions
Figu e 6.6
S abilisa ion imes a 400 km al i ude wi h hys e esis
Figu e 6.7
Ro a ion speed a 400 km al i ude wi
• 500 km:
his simula ion uns o 0.
beha iou shown in
wo cases al eady
discussed
Simula ions
S abilisa ion imes a 400 km al i ude wi h hys e esis
Ro a ion speed a 400 km al i ude wi
h hys e esis
his simula ion uns o 0.
25 days wi h a ime
s ep o
Figu e 6.8 and Figu e 6.9
is comple ely analogous o he
discussed
.
33
S abilisa ion imes a 400 km al i ude wi h hys e esis
.
h hys e esis
.
s ep o
2 seconds. The
is comple ely analogous o he
34
Figu e 6.8
S abilisa ion imes a 500 km al i ude wi h hys e esis
Figu e 6.9
Ro a ion speed a 500 km al i ude wi h hys e esis
Magne ic sys ems o he ADCS o a Fem osa elli e
S abilisa ion imes a 500 km al i ude wi h hys e esis
Ro a ion speed a 500 km al i ude wi h hys e esis
Magne ic sys ems o he ADCS o a Fem osa elli e
S abilisa ion imes a 500 km al i ude wi h hys e esis
.
Ro a ion speed a 500 km al i ude wi h hys e esis
.
Magne ic ADCS Pe o mance Simula ions 35
6.2.2. Magne o que s
• 300 km: his simula ion is un o 0.012 days wi h a ime s ep o 1 second,
ob aining Figu e 6.10.
Figu e 6.10 Ro a ion speed a 300 km al i ude wi h magne o que s.
36 Magne ic sys ems o he ADCS o a Fem osa elli e
• 400 km: his simula ion is un o 0.03 days wi h a ime s ep o 1 second,
ob aining Figu e 6.11.
Figu e 6.11 Ro a ion speed a 400 km al i ude wi h magne o que s.
Magne ic ADCS Pe o mance Simula ions 37
• 500 km: his simula ion is un o 0.02 days wi h a ime s ep o 1 second
Figu e 6.12.
Figu e 6.12 Ro a ion speed a 500 km al i ude wi h magne o que s.
6.3. Analysis o esul s
The esul s o he p e ious g aphs a e summa ised in se e al ables, Table 6.1 o
he hys e esis ods and Table 6.2 o he magne o que s.
Table 6.1 Hys e esis ods esul s.
Al i ude
P opaga o
Magne ic
model S abiliza ion
Resul
300 km
SGP4 WMM
4h 25min The e is no signi ican di e ence
be ween he ou di e en models a
each heigh .
The e is a decline in o a ion speed
whose s eepness oscilla es due o
he a ia ion in magne ic along he
o bi .
IGRF
SGP8 WMM
IGRF
400 km
SGP4 WMM
4h 40min
IGRF
SGP8 WMM
IGRF

38 Magne ic sys ems o he ADCS o a Fem osa elli e
Table 6.1 Hys e esis ods esul s. (Con .)
Al i ude
P opaga o
Magne ic
model S abiliza ion
Resul
500 km
SGP4 WMM
4h 50min As nea he Ea h he magne ic ield
is s onge , he dampening e ec is
also mo e signi ican .
IGRF
SGP8 WMM
IGRF
Table 6.2 Magne o que esul s.
Al i ude
P opaga o
Magne ic
model S abiliza ion
Resul
300 km
SGP4 WMM 26 min
Fo he i s 2 s all plo s a e he
same.
IGRF 26 min
SGP8 WMM 16 min
IGRF 11 min
400 km
SGP4 WMM 11 min 25 s
All excep SGP4 WMM a e he same
o 42 s.
IGRF 37 min
SGP8 WMM 10 min 50 s
IGRF 21 min
500 km
SGP4 WMM 15 min 50 s
All a e equal o he i s 54 s.
IGRF 26 min 30 s
SGP8 WMM 14 min 50 s
IGRF 15 min
In gene al, i can be seen ha hys e esis ods ake longe han magne o que s o
s abilize he sa elli e, a esul in line wi h ou expec a ions. All he magne o que
g aphs ha e a e y s eep dec ease in o a ion speed a he beginning, ha inally
ge s he sa elli e o speeds o unde 1 deg ee pe second in less han an hou , and in
mos cases unde 20 minu es.
The powe consump ion o he magne o que s can be calcula ed as s a ed in
equa ion 4.5 and plo ed o each one o he simula ions.
Fo he 300km case, Figu e 6.13 is p oduced ( he i s peak is c opped ou o make
he es o he g aph no iceable). In his case he e is a la ge 250mW peak a he
beginning, and hen se e al smalle ones ha consume less han 25 mW.
Magne ic ADCS Pe o mance Simula ions 39
Figu e 6.13 Powe used du ing s abiliza ion in he 300 km al i ude o bi .
Fo he 400 km case we ob ain Figu e 6.14 (again, he i s peak is c opped ou o
make he es o he g aph no iceable); now he peak is smalle han in he p e ious
case, being jus a 157mW peak. Fu he on in he simula ion, he e appea se e al
smalle peaks o less han 20 mW. The es o he alues will only consume less han
5 mW.
40 Magne ic sys ems o he ADCS o a Fem osa elli e
Figu e 6.14 Powe used a he400 km al i ude o bi .
Fo he 500 km case, as seen in Figu e6.15, he ea ly peak is jus 67 mW high
(ne e heless, his i s peak is c opped ou o make he es o he g aph no iceable),
Seconda y peaks a e a smalle , wi h a maximum 25 mW peak a 48 seconds, and
se e al peaks o less han 5 mW
Magne ic ADCS Pe o mance Simula ions 41
Figu e6.15 Powe consumed 500 km al i ude.
Figu e 6.13, Figu e 6.14 and Figu e6.15 can be summa ised in Table 6.3.
Table 6.3 Powe consump ion o each simula ion.
Al i ude Bigges peak
Residual Analysis
300 km 250 mW <1 mW Has a e y la ge peak o 250 mW and o he
4 smalle ones o less han 25 mW, bu in
gene al he powe consump ion is low.
400 km 157 mW <1 mW Has a la ge peak o 157 mW and a couple
mo e o less han 10 mW, bu in gene al he
powe consump ion is low.
500 km 67 mW <1 mW Has a la ge peak o 67 mW o he one o 24
m-W and a se e al bellow 5 mW, bu in
gene al he powe consump ion is low.
Rega ding he equi emen s, he Table 6.4 analyses he compliance o bo h sys ems
o he equi emen s applicable o he subsys em as s udied in his documen .
48 Magne ic sys ems o he ADCS o a Fem osa elli e
%--------------------------------------------------------------------------
%--------------------------Ge ini ial pa ame e s--------------------------
%--------------------------------------------------------------------------
% Re e ence ec o o he ECI ame coo dina e sys em
Ine ial ame=[1,1,1]/(sq (3));
% Ob ain o bi da a o m TLE
[sa da a]= O bi alp ( name);
% Change Tenso o Ine ia in uni s om [g*mm^2] o [kg*m^2]
I=I/(10^9); % [kg*m^2]
% Time pa ame e s o he p opaga ion o sa elli e's s a e ec o
% o wa d (+) o backwa d (-) [minu es]
since = p opdays*1440; % P opaga ion ime [minu es]
s ep=s ep_sec/60; % S ep [minu es]
=0; % Fi s s ep
len= ix( since/s ep); % Numbe o i e a ions
%-----------------------Gene a e he s o age a ays------------------------
% SGP4 WMM o a ion speed s o age
omega14X(len)=0; %S o age o X alue
omega14Y(len)=0; %S o age o Y alue
omega14Z(len)=0; %S o age o Z alue
omega14M(len)=0; %S o age o module
% SGP4 IGRF o a ion speed s o age
omega24X(len)=0; %S o age o X alue
omega24Y(len)=0; %S o age o Y alue
omega24Z(len)=0; %S o age o Z alue
omega24M(len)=0; %S o age o module
% SGP8 WMM o a ion speed s o age
omega18X(len)=0; %S o age o X alue
omega18Y(len)=0; %S o age o Y alue
omega18Z(len)=0; %S o age o Z alue
omega18M(len)=0; %S o age o module
% SGP8 IGRF o a ion speed s o age
omega28X(len)=0; %S o age o X alue
omega28Y(len)=0; %S o age o Y alue
omega28Z(len)=0; %S o age o Z alue
omega28M(len)=0; %S o age o module
% Powe consumed s o age
powe c14(len)=0; % SGP4 WMM
powe c24(len)=0; % SGP4 IGRF
powe c18(len)=0; % SGP8 WMM
powe c28(len)=0; % SGP8 IGRF
% S o e ime o he i e a ion
imed(len)=0;
imeh(len)=0;
imem(len)=0;
% S abilisa ion da e s o age
s abiliza ion(4)=0;
% I e a ion coun e ini ializa ion
i e a ion=1;

Annexes 49
%--------------------------Ge he beginning yea --------------------------
i (sa da a.yea < 57)
yea 0 = sa da a.yea + 2000;
else
yea 0 = sa da a.yea + 1900;
end
%--------------------------------------------------------------------------
%-----------------------------------Loop-----------------------------------
%--------------------------------------------------------------------------
while < since
%-------------------------------Time upda e--------------------------------
% Calcula e da e o cu en i e a ion
i sa da a.doy+( /1440) < 365
doy = (sa da a.doy+( /1440));
yea =yea 0;
else
doy = (sa da a.doy+( /1440)-365);
yea =yea 0+1;
end
% Calcula e da e pa ame e s
[mon,day,h ,minu e,sec] = days2mdh(yea ,doy);
% Da e as decimal yea
dy = decyea (yea ,mon,day,h ,minu e,sec);
% Gene a e UTC ec o
u c=[yea ,mon,day,h ,minu e,sec];
%---------------------------P opaga e he o bi ----------------------------
[ eme4, eme4] = sgp4( , sa da a); % SGP4 [km, km/s]
[ eme8, eme8] = sgp8( , sa da a); % SGP8 [km, km/s]
% Con e posi ion ec o uni s om [km] o [m]
eme4= eme4*1000; % [m]
eme8= eme8*1000; % [m]
% Calcula e LLA coo dina es SGP4
lla4=eci2lla([ eme4(1), eme4(2), eme4(3)],u c);
la 4=lla4(1); lon4=lla4(2); h4=lla4(3);
% Calcula e LLA coo dina es SGP8
lla8=eci2lla([ eme8(1), eme8(2), eme8(3)],u c);
la 8=lla8(1); lon8=lla8(2); h8=lla8(3);
%------------------------Calcula e magne ic ields-------------------------
% SGP4 [nT]
[XYZ14, H14, D14, I14, F14] = w ldmagm(h4,la 4,lon4,dy,'2020'); % WMM
[XYZ24,H24,D24,I24,F24] = ig magm(h4,la 4,lon4,dy,13); % IGRF
XYZ24= anspose(XYZ24); % Adap magne ic ield ec o o ma
% SGP8 [nT]
[XYZ18, H18, D18, I18, F18] = w ldmagm(h8,la 8,lon8,dy,'2020'); % WMM
[XYZ28,H28,D28,I28,F28] = ig magm(h8,la 8,lon8,dy,13); % IGRF
XYZ28= anspose(XYZ28); % Adap magne ic ield ec o o ma
%--------------------T ans o m magne ic ields in o BRF--------------------
%...............................NED o ECEF................................
% SGP4 WMM [nT]
[ECEF14(1),ECEF14(2),ECEF14(3)]=ned2ece (XYZ14(1),XYZ14(2),XYZ14(3),la
4,lon4);
% SGP4 IGRF [nT]
[ECEF24(1),ECEF24(2),ECEF24(3)]=ned2ece (XYZ24(1),XYZ24(2),XYZ24(3),la
4,lon4);
50 Magne ic sys ems o he ADCS o a Fem osa elli e
% SGP8 WMM [nT]
[ECEF18(1),ECEF18(2),ECEF18(3)]=ned2ece (XYZ18(1),XYZ18(2),XYZ18(3),la
8,lon8);
% SGP8 IGRF [nT]
[ECEF28(1),ECEF28(2),ECEF28(3)]=ned2ece (XYZ28(1),XYZ28(2),XYZ28(3),la
8,lon8);
%...............................ECEF o ECI................................
eci14=ece 2eci(u c,ECEF14); % SGP4 WMM [nT]
eci24=ece 2eci(u c,ECEF24); % SGP4 IGRF [nT]
eci18=ece 2eci(u c,ECEF18); % SGP8 WMM [nT]
eci28=ece 2eci(u c,ECEF28); % SGP8 IGRF [nT]
%....................Qua e nion o ECI o BRF o a ion....................
% SGP4 WMM
S14=Sx14+Sy14+Sz14; % Gene a e sa elli e o ien a ion ec o
S14=S14/no m(S14); % No malise o ien a ion ec o
p 14=c oss(Ine ial ame,S14); % Calcula e o a ion axis
% Calcula e o a ion angle
i S14(1)+S14(2)+S14(3)<0
i14=deg2 ad(180)-
(asin(no m(p 14))/(no m(Ine ial ame)*no m(S14)));
else
i14=asin(no m(p 14)/(no m(Ine ial ame)*no m(S14)));
end
% Angle o e e se o a ion
an i i14=- i14;
% Calcula e qua e nion scala pa
p014=cos( i14/2);
an ip014=cos(an i i14/2);
% Calcula e qua e nion ec o pa
i no m(p 14)<0.0000001
p14= anspose((p 14)*0);
an ip14=p14;
else
p14= anspose((p 14/no m(p 14))*sin( i14/2));
an ip14= anspose((p 14/no m(p 14))*sin(an i i14/2));
end
% SGP4 IGRF
S24=Sx24+Sy24+Sz24; % Gene a e sa elli e o ien a ion ec o
S24=S24/no m(S24); % No malise o ien a ion ec o
p 24=c oss(Ine ial ame,S24); % Calcula e o a ion axis
% Calcula e o a ion angle
i S24(1)+S24(2)+S24(3)<0
i24=deg2 ad(180)-
(asin(no m(p 24))/(no m(Ine ial ame)*no m(S24)));
else
i24=asin(no m(p 24)/(no m(Ine ial ame)*no m(S24)));
end
% Angle o e e se o a ion
an i i24=- i24;
% Calcula e qua e nion scala pa
p024=cos( i24/2);
an ip024=cos(an i i24/2);
% Calcula e qua e nion ec o pa
i no m(p 24)<0.0000001
p24= anspose((p 24)*0);
an ip24=p24;
else
p24= anspose((p 24/no m(p 24))*sin( i24/2));
an ip24= anspose((p 24/no m(p 24))*sin(an i i24/2));
end
Annexes 51
% SGP8 WMM
S18=Sx18+Sy18+Sz18; % Gene a e sa elli e o ien a ion ec o
S18=S18/no m(S18); % No malise o ien a ion ec o
p 18=c oss(Ine ial ame,S18); % Calcula e o a ion axis
% Calcula e o a ion angle
i S18(1)+S18(2)+S18(3)<0
i18=deg2 ad(180)-
(asin(no m(p 18))/(no m(Ine ial ame)*no m(S18)));
else
i18=asin(no m(p 18)/(no m(Ine ial ame)*no m(S18)));
end
% Angle o e e se o a ion
an i i18=- i18;
% Calcula e qua e nion scala pa
p018=cos( i18/2);
an ip018=cos(an i i18/2);
% Calcula e qua e nion ec o pa
i no m(p 18)<0.0000001
p18= anspose((p 18)*0);
an ip18=p18;
else
p18= anspose((p 18/no m(p 18))*sin( i18/2));
an ip18= anspose((p 18/no m(p 18))*sin(an i i18/2));
end
% SGP8 IGRF
S28=Sx28+Sy28+Sz28; % Gene a e sa elli e o ien a ion ec o
S28=S28/no m(S28); % No malise o ien a ion ec o
p 28=c oss(Ine ial ame,S28); % Calcula e o a ion axis
% Calcula e o a ion angle
i S28(1)+S28(2)+S28(3)<0
i28=deg2 ad(180)-
(asin(no m(p 28))/(no m(Ine ial ame)*no m(S28)));
else
i28=asin(no m(p 28)/(no m(Ine ial ame)*no m(S28)));
end
% Angle o e e se o a ion
an i i28=- i28;
an ip028=cos(an i i28/2);
% Calcula e qua e nion scala pa
p028=cos( i28/2);
% Calcula e qua e nion ec o pa
i no m(p 28)<0.0000001
p28= anspose((p 28)*0);
an ip28=p28;
else
p28= anspose((p 28/no m(p 28))*sin( i28/2));
an ip28= anspose((p 28/no m(p 28))*sin(an i i28/2));
end
%................................ECI o BRF................................
b 14=qua _ o a ion(p014,p14,eci14); %SGP4 WMM [nT]
b 24=qua _ o a ion(p024,p24,eci24); %SGP4 IGRF [nT]
b 18=qua _ o a ion(p018,p18,eci18); %SGP8 WMM [nT]
b 28=qua _ o a ion(p028,p28,eci28); %SGP8 IGRF [nT]
% Adap uni s o m [nT] o [T]
b 14T=b 14/(10^9); %SGP4 WMM [T]
b 24T=b 24/(10^9); %SGP4 IGRF [T]
b 18T=b 18/(10^9); %SGP8 WMM [T]
b 28T=b 28/(10^9); %SGP8 IGRF [T]
52 Magne ic sys ems o he ADCS o a Fem osa elli e
%-------------------Calcula e mu alues o he i e a ion------------------
%..........................Hys e esis calcula ion..........................
i ype==1 %I i is hys e esis
% SGP4 WMM
EK 14=0.5* o speed14*I* anspose( o speed14); % Ro a ional kine ic
ene gy [J]
Del aEK14=Chys *no m(b 14T)*no m( ad2deg( o speed14)); % Kine ic
ene gy a ia ion [J]
EK 14=EK 14-Del aEK14; % Upda e kine ic ene gy [J]
i EK 14<0 % Ve i y i all has dissipa ed
o speed14= o speed14*0; % Finish he calcula ions
s abiliza ion(1)=1; % S abilisa ion o be check
else
Ine ia=do (S14*I,S14); % Momen o ine ia a he
axis o o a ion [kg*m^2]
o speed14=(S14)*sq (2*EK 14/Ine ia); % Upda e o a ion speed
end
% SGP4 IGRF
EK 24=0.5* o speed24*I* anspose( o speed24); % Ro a ional kine ic
ene gy [J]
Del aEK24=Chys *no m(b 24T)*no m( ad2deg( o speed24)); % Kine ic
ene gy a ia ion [J]
EK 24=EK 24-Del aEK24; % Upda e kine ic ene gy [J]
i EK 24<0 % Ve i y i all has dissipa ed
o speed24= o speed24*0; % Finish he calcula ions
s abiliza ion(2)=1; % S abilisa ion o be s o ed
else
Ine ia=do (S24*I,S24); % Momen o ine ia a he
axis o o a ion [kg*m^2]
o speed24=(S24)*sq (2*EK 24/Ine ia); % Upda e o a ion speed
end
% SGP8 WMM
EK 18=0.5* o speed18*I* anspose( o speed18); % Ro a ional kine ic
ene gy [J]
Del aEK18=Chys *no m(b 18T)*no m( ad2deg( o speed18)); % Kine ic
ene gy a ia ion [J]
EK 18=EK 18-Del aEK18; % Upda e kine ic ene gy [J]
i EK 18<0 % Ve i y i all has
dissipa ed
o speed18= o speed18*0; % Finish he calcula ions
s abiliza ion(3)=1; % S abilisa ion o be s o ed
else
Ine ia=do (S18*I,S18); % Momen o ine ia a he
axis o o a ion [kg*m^2]
o speed18=(S18)*sq (2*EK 18/Ine ia);% Upda e o a ion speed
end
Annexes 53
% SGP8 IGRF
EK 28=0.5* o speed28*I* anspose( o speed28); % Ro a ional kine ic
ene gy [J]
Del aEK28=Chys *no m(b 28T)*no m( ad2deg( o speed28)); % Kine ic
ene gy a ia ion [J]
EK 28=EK 28-Del aEK28; % Upda e kine ic ene gy [J]
i EK 28<0 % Ve i y i all has dissipa ed
o speed28= o speed28*0; % Finish he calcula ions
s abiliza ion(4)=1; % S abilisa ion o be s o ed
else
Ine ia=do (S28*I,S28); % Momen o ine ia a he
axis o o a ion [kg*m^2]
o speed28=(S28)*sq (2*EK 28/Ine ia); % Upda e o a ion speed
end
% Sa e s abilisa ion da es
i s abiliza ion(1)==1
s abda e(1)= /1440;
s abiliza ion(1)=s abiliza ion(1)+1;
end
i s abiliza ion(2)==1
s abda e(2)= /1440;
s abiliza ion(2)=s abiliza ion(2)+1;
end
i s abiliza ion(3)==1
s abda e(3)= /1440;
s abiliza ion(3)=s abiliza ion(3)+1;
end
i s abiliza ion(4)==1
s abda e(4)= /1440;
s abiliza ion(4)=s abiliza ion(4)+1;
end
i s abiliza ion(1)>1 && s abiliza ion(2)>1 && s abiliza ion(3)>1 &&
s abiliza ion(4)>1
= since;
end
%.........................Magne o que calcula ion.........................
elsei ype==2 % I i is magne o que
% SGP4 WMM
accele a ion14=(- o speed14)/s ep_sec; % [ ad/s^2]
To que14=I* anspose(accele a ion14); % [N*m]
muha 14=c oss(b 14T,To que14)/no m(c oss(b 14T,To que14));
k14=((do (muha 14,b 14T)*b 14T)/do (b 14T,b 14T))-muha 14;
mu14=(no m(c oss(To que14,b 14T))/(no m(k14)*do (b 14T,b 14T)));
mu14=mu14*muha 14; % [A*m^2]
% Ve i y mu is no o e he max alue
i (abs(mu14(1))+abs(mu14(2))+abs(mu14(3)))>mumax
muk=mumax/(abs(mu14(1))+abs(mu14(2))+abs(mu14(3)));
mu14(1)=mu14(1)*muk;
mu14(2)=mu14(2)*muk;
mu14(3)=mu14(3)*muk;
end

54 Magne ic sys ems o he ADCS o a Fem osa elli e
% SGP4 IGRF
accele a ion24=(- o speed24)/s ep_sec; % [ ad/s^2]
To que24=I* anspose(accele a ion24); % [N*m]
muha 24=c oss(b 24T,To que24)/no m(c oss(b 24T,To que24));
k24=((do (muha 24,b 24T)*b 24T)/do (b 24T,b 24T))-muha 24;
mu24=(no m(c oss(To que24,b 24T))/(no m(k24)*do (b 24T,b 24T)));
mu24=mu24*muha 24; % [A*m^2]
% Ve i y mu is no o e he max alue
i (abs(mu24(1))+abs(mu24(2))+abs(mu24(3)))>mumax
muk=mumax/(abs(mu24(1))+abs(mu24(2))+abs(mu24(3)));
mu24(1)=mu24(1)*muk;
mu24(2)=mu24(2)*muk;
mu24(3)=mu24(3)*muk;
end
% SGP8 WMM
accele a ion18=(- o speed18)/s ep_sec; % [ ad/s^2]
To que18=I* anspose(accele a ion18); % [N*m]
muha 18=c oss(b 18T,To que18)/no m(c oss(b 18T,To que18));
k18=((do (muha 18,b 18T)*b 18T)/do (b 18T,b 18T))-muha 18;
mu18=(no m(c oss(To que18,b 18T))/(no m(k18)*do (b 18T,b 18T)));
mu18=mu18*muha 18; % [A*m^2]
% Ve i y mu is no o e he max alue
i (abs(mu18(1))+abs(mu18(2))+abs(mu18(3)))>mumax
muk=mumax/(abs(mu18(1))+abs(mu18(2))+abs(mu18(3)));
mu18(1)=mu18(1)*muk;
mu18(2)=mu18(2)*muk;
mu18(3)=mu18(3)*muk;
end
% SGP8 IGRF
accele a ion28=(- o speed28)/s ep_sec; % [ ad/s^2]
To que28=I* anspose(accele a ion28); % [N*m]
muha 28=c oss(b 28T,To que28)/no m(c oss(b 28T,To que28));
k28=((do (muha 28,b 28T)*b 28T)/do (b 28T,b 28T))-muha 28;
mu28=(no m(c oss(To que28,b 28T))/(no m(k28)*do (b 28T,b 28T)));
mu28=mu28*muha 28; % [A*m^2]
i (abs(mu28(1))+abs(mu28(2))+abs(mu28(3)))>mumax
muk=mumax/(abs(mu28(1))+abs(mu28(2))+abs(mu28(3)));
mu28(1)=mu28(1)*muk;
mu28(2)=mu28(2)*muk;
mu28(3)=mu28(3)*muk;
end
%--------------------Calcula e o a ion o he main axis-------------------
% Calcula e o ques
T14=c oss(mu14,b 14T); %[N*m] SGP4 WMM
T24=c oss(mu24,b 24T); %[N*m] SGP4 IGRF
T18=c oss(mu18,b 18T); %[N*m] SGP8 WMM
T28=c oss(mu28,b 28T); %[N*m] SGP8 IGRF
% Calcula e accele a ion ec o s
accele a ion14= anspose(linsol e(I,T14)); %[ ad/s^2] SGP4 WMM
accele a ion24= anspose(linsol e(I,T24)); %[ ad/s^2] SGP4 IGRF
accele a ion18= anspose(linsol e(I,T18)); %[ ad/s^2] SGP8 WMM
accele a ion28= anspose(linsol e(I,T28)); %[ ad/s^2] SGP8 IGRF
Annexes 55
% Calcula e angle u ned
% SGP4 WMM [ ad]
he a14= o speed14*s ep_sec+(1/2)*accele a ion14*s ep_sec*s ep_sec;
% SGP4 IGRF [ ad]
he a24= o speed24*s ep_sec+(1/2)*accele a ion24*s ep_sec*s ep_sec;
% SGP8 WMM [ ad]
he a18= o speed18*s ep_sec+(1/2)*accele a ion18*s ep_sec*s ep_sec;
% SGP8 IGRF [ ad]
he a28= o speed28*s ep_sec+(1/2)*accele a ion28*s ep_sec*s ep_sec;
% Calcula e he o a ion axis o he i e a ion % [uni ec o ]
o axis14= o speed14/no m( o speed18); % SGP4 WMM
o axis24= o speed24/no m( o speed24); % SGP4 IGRF
o axis18= o speed18/no m( o speed18); % SGP8 WMM
o axis28= o speed28/no m( o speed28); % SGP8 IGRF
% Calcula e he s ep qua e nion
% SGP4 WMM
q014=cos(no m( he a14)/2); % Qua e nion scala pa
q14= o axis14*(sin(no m( he a14)/2)); % Qua e nion ec o pa
% SGP4 IGRF
q024=cos(no m( he a24)/2); % Qua e nion scala pa
q24= o axis24*(sin(no m( he a24)/2)); % Qua e nion ec o pa
% SGP8 WMM
q018=cos(no m( he a18)/2); % Qua e nion scala pa
q18= o axis18*(sin(no m( he a18)/2)); % Qua e nion ec o pa
% SGP8 IGRF
q028=cos(no m( he a28)/2); % Qua e nion scala pa
q28= o axis28*(sin(no m( he a28)/2)); % Qua e nion ec o pa
% T ans o m sa elli es main axis o ECI
% SGP4 WMM
Sx14=qua _ o a ion(q014,q14,Sx14); % Sa elli e X axis in ECI
Sy14=qua _ o a ion(q014,q14,Sy14); % Sa elli e Y axis in ECI
Sz14=qua _ o a ion(q014,q14,Sz14); % Sa elli e Z axis in ECI
% SGP4 IGRF
Sx24=qua _ o a ion(q024,q24,Sx24); % Sa elli e X axis in ECI
Sy24=qua _ o a ion(q024,q24,Sy24); % Sa elli e Y axis in ECI
Sz24=qua _ o a ion(q024,q24,Sz24); % Sa elli e Z axis in ECI
% SGP8 WMM
Sx18=qua _ o a ion(q018,q18,Sx18); % Sa elli e X axis in ECI
Sy18=qua _ o a ion(q018,q18,Sy18); % Sa elli e Y axis in ECI
Sz18=qua _ o a ion(q018,q18,Sz18); % Sa elli e Z axis in ECI
% SGP8 IGRF
Sx28=qua _ o a ion(q028,q28,Sx28); % Sa elli e X axis in ECI
Sy28=qua _ o a ion(q028,q28,Sy28); % Sa elli e Y axis in ECI
Sz28=qua _ o a ion(q028,q28,Sz28); % Sa elli e Z axis in ECI
% Calcula e new speed o o a ion
o speed14= o speed14+accele a ion14*s ep_sec; % [ ad/s] SGP4 WMM
o speed24= o speed24+accele a ion24*s ep_sec; % [ ad/s] SGP4 IGRF
o speed18= o speed18+accele a ion18*s ep_sec; % [ ad/s] SGP8 WMM
o speed28= o speed28+accele a ion28*s ep_sec; % [ ad/s] SGP8 IGRF
else
56 Magne ic sys ems o he ADCS o a Fem osa elli e
p in ('W ong ype o sys em, ype mus be be ween 1 and 3');
= since;
end
%--------------------------------Sa e da a---------------------------------
% Ro a ion a es
% SGP4 WMM
omega14X(i e a ion)= ad2deg( o speed14(1)); % [deg/s]
omega14Y(i e a ion)= ad2deg( o speed14(2)); % [deg/s]
omega14Z(i e a ion)= ad2deg( o speed14(3)); % [deg/s]
omega14M(i e a ion)= ad2deg(no m( o speed14)); % [deg/s]
% SGP4 IGRF
omega24X(i e a ion)= ad2deg( o speed24(1)); % [deg/s]
omega24Y(i e a ion)= ad2deg( o speed24(2)); % [deg/s]
omega24Z(i e a ion)= ad2deg( o speed24(3)); % [deg/s]
omega24M(i e a ion)= ad2deg(no m( o speed24)); % [deg/s]
% SGP8 WMM
omega18X(i e a ion)= ad2deg( o speed18(1)); % [deg/s]
omega18Y(i e a ion)= ad2deg( o speed18(2)); % [deg/s]
omega18Z(i e a ion)= ad2deg( o speed18(3)); % [deg/s]
omega18M(i e a ion)= ad2deg(no m( o speed18)); % [deg/s]
% SGP8 IGRF
omega28X(i e a ion)= ad2deg( o speed28(1)); % [deg/s]
omega28Y(i e a ion)= ad2deg( o speed28(2)); % [deg/s]
omega28Z(i e a ion)= ad2deg( o speed28(3)); % [deg/s]
omega28M(i e a ion)= ad2deg(no m( o speed28)); % [deg/s]
% S o e powe consump ion o magne o que s [mW]
i ype==2
powe c14(i e a ion)=((mu14(1)^2)+(mu14(2)^2)+(mu14(3))^2)/(pow cons
^2);
powe c24(i e a ion)=((mu24(1)^2)+(mu24(2)^2)+(mu24(3))^2)/(pow cons
^2);
powe c18(i e a ion)=((mu18(1)^2)+(mu18(2)^2)+(mu18(3))^2)/(pow cons
^2);
powe c28(i e a ion)=((mu28(1)^2)+(mu28(2)^2)+(mu28(3)^2))/(pow cons
^2);
end
% S o e ime alue
imed(i e a ion)= /1440; %[days]
imeh(i e a ion)= /60; %[hou s]
imem(i e a ion)= ; %[min]
%------------------------------Upda e coun e ------------------------------
i e a ion % P in he i e a ion numbe
i e a ion=i e a ion+1;
= +s ep;
end
Annexes 57
%--------------------------------------------------------------------------
%------------------------------Gene a e plo s------------------------------
%--------------------------------------------------------------------------
i ype==1 % Plo i alid hys e esis calcula ion
igu e
plo ( imeh,omega28M, imeh,omega18M, imeh,omega24M, imeh,omega14M,'LineWi
d h',1.1)
i le('Hys e esis 500 km')
xlabel('Time elapsed [hou s]')
ylabel('Speed o o a ion [deg/s]')
legend('SGP8 IGRF','SGP8 WMM','SGP4 IGRF','SGP4 WMM')
se (gca,'linewid h',1.5)
% P in he hys e esis s abiliza ion ime
i s abiliza ion(1)==2
p in ('Sa elli e hys e esis s abilised a e %.4 days. SGP4
WMM n',s abda e(1));
else
p in ('Sa elli e hys e esis does no s abilise a e %.4 days.
SGP4 WMM n', since);
end
i s abiliza ion(2)==2
p in ('Sa elli e hys e esis s abilised a e %.4 days. SGP4
IGRF n',s abda e(2));
else
p in ('Sa elli e hys e esis does no s abilise a e %.4 days.
SGP4 IGRF n', since);
end
i s abiliza ion(3)==2
p in ('Sa elli e hys e esis s abilised a e %.4 days. SGP8
WMM n',s abda e(3));
else
p in ('Sa elli e hys e esis does no s abilise a e %.4 days.
SGP8 WMM n', since);
end
i s abiliza ion(4)==2
p in ('Sa elli e hys e esis s abilised a e %.4 days. SGP8
IGRF n',s abda e(4));
else
p in ('Sa elli e hys e esis does no s abilise a e %.4 days.
SGP4 IGRF n', since);
end
elsei ype==2 % Plo i alid magne o que calcula ion
igu e
plo ( imem,omega28M, imem,omega18M, imem,omega24M, imem,omega14M,'LineWi
d h',1.1)
i le('Magne o que s 500 km')
xlabel('Time elapsed [min]')
ylabel('Speed o o a ion [deg/s]')
legend('SGP8 IGRF','SGP8 WMM','SGP4 IGRF','SGP4 WMM')
se (gca,'linewid h',1.5)
64 Magne ic sys ems o he ADCS o a Fem osa elli e

Annexes 65
Magne ic ield ec o o a ion e i ica ion code
%--------------------------------------------------------------------------
% Ro a ion e i ica ion subcode
%
%--------------------------------------------------------------------------
clc
clea
o ma long g
%--------------------------------------------------------------------------
%----------------------------------INPUTS----------------------------------
%--------------------------------------------------------------------------
% TLE iles names
name = 'SATELIOT_1. x ';
% Tenso o Ine ia [g*mm^2]
I=[[37971.2,-17.4,-15.7];
[-17.4,37596.6,238.8];
[-15.7,238.8,37727.4]];
% Ini ial o a ion speed
o speed=[1,1,1]*deg2 ad(20); % [ ad/s]
% P opaga ion pa ame e s
p opdays= 0.01; % p opaga ion ime [days]
s ep_sec = 1; % s ep [seconds]
% Dipola momen
mu=[0,0.2,0]; % [A*m^2]
% Sa elli e main axis a ECI ame a =0
Sx=[1,0,0]; % X axis [uni ]
Sy=[0,1,0]; % Y axis [uni ]
Sz=[0,0,1]; % Z axis [uni ]
%--------------------------------------------------------------------------
%--------------------------Ge ini ial pa ame e s--------------------------
%--------------------------------------------------------------------------
% Re e ence ec o o he ECI ame coo dina e sys em
Ine ial ame=[1,1,1]/(sq (3)); % [uni ]
% Ob ain o bi da a o m TLE
[sa da a]= O bi alp ( name);
% Change Tenso o Ine ia in uni s om [g*mm^2] o [kg*m^2]
I=I/(10^9); % [kg*m^2]
% Time pa ame e s o he p opaga ion o sa elli e's s a e ec o
% o wa d (+) o backwa d (-) [minu es]
since = p opdays*1440; % P opaga ion ime [minu es]
s ep=s ep_sec/60; % S ep [minu es]
=0; % Fi s s ep
len= ix( since/s ep)-1; % Numbe o i e a ions
66 Magne ic sys ems o he ADCS o a Fem osa elli e
%-----------------------Gene a e he s o age a ays------------------------
% Di e ence in o a ion o he magne ic ield s o age
es 1(len)=0; % S o age o e o om NED o ECEF
es 2(len)=0; % S o age o e o om ECEF o ECI
es 3(len)=0; % S o age o e o om ECI o BRF
% S o e ime o he i e a ion
ime(len)=0;
% I e a ion coun e ini ializa ion
i e a ion=1;
%--------------------------Ge he beginning yea --------------------------
i (sa da a.yea < 57)
yea 0 = sa da a.yea + 2000;
else
yea 0 = sa da a.yea + 1900;
end
%--------------------------------------------------------------------------
%-----------------------------------Loop-----------------------------------
%--------------------------------------------------------------------------
while < since
%-------------------------------Time upda e--------------------------------
% Calcula e da e o i e a ion
i sa da a.doy+( /1440) < 365
doy = (sa da a.doy+( /1440));
yea =yea 0;
else
doy = (sa da a.doy+( /1440)-365);
yea =yea 0+1;
end
% Calcula e da e pa ame e s
[mon,day,h ,minu e,sec] = days2mdh(yea ,doy);
% Da e as decimal yea
dy = decyea (yea ,mon,day,h ,minu e,sec);
% Gene a e UTC ec o
u c=[yea ,mon,day,h ,minu e,sec];
%---------------------------P opaga e he o bi ----------------------------
[ eme, eme] = sgp8( , sa da a); % [km, km/s]
% Con e posi ion ec o uni s om [km] o [m]
eme= eme*1000; % [m]
% Calcula e LLA coo dina es
lla=eci2lla([ eme(1), eme(2), eme(3)],u c);
la =lla(1); lon=lla(2); h=lla(3);
%------------------------Calcula e magne ic ields-------------------------
% SGP8 [nT]
[XYZ,H2,D2,I2,F2] = ig magm(h,la ,lon,dy,13); %IGRF
XYZ= anspose(XYZ); % Adap magne ic ield ec o o ma
Annexes 67
%--------------------T ans o m magne ic ields in o BRF--------------------
%...............................NED o ECEF................................
[ECEF(1),ECEF(2),ECEF(3)]=ned2ece (XYZ(1),XYZ(2),XYZ(3),la ,lon);%[nT]
%...............................ECEF o ECI................................
ECI=ece 2eci(u c,ECEF); % [nT]
%....................Qua e nion o ECI o BRF o a ion....................
S=Sx+Sy+Sz; % Gene a e sa elli e o ien a ion ec o
S=S/no m(S); % No malise o ien a ion ec o
p =c oss(Ine ial ame,S); % Calcula e o a ion axis
% Calcula e o a ion angle
i S(1)+S(2)+S(3)<0
i=deg2 ad(180)-(asin(no m(p ))/(no m(Ine ial ame)*no m(S)));
else
i=asin(no m(p )/(no m(Ine ial ame)*no m(S)));
end
% Angle o e e se o a ion
an i i=- i;
% Calcula e qua e nion scala pa
p0=cos( i/2);
an ip0=cos(an i i/2);
% Calcula e qua e nion ec o pa
i no m(p )<0.0000001
p= anspose((p )*0);
an ip=p;
else
p= anspose((p /no m(p ))*sin( i/2));
an ip= anspose((p /no m(p ))*sin(an i i/2));
end
%................................ECI o BRF................................
BRF=qua _ o a ion(p0,p,ECI); % [nT]
an iBRF=qua _ o a ion(an ip0,an ip,BRF);
% Adap uni s o m [nT] o [T]
BRF_ esla=BRF/(10^9); % [T]
%--------------------Calcula e o a ion o he main axis-------------------
% Calcula e o ques
T= anspose(c oss(mu,BRF_ esla)); % [N*m]
% Calcula e accele a ion ec o s
accele a ion= anspose(linsol e(I,T)); % [ ad/s^2]
% Ro a e accele a ion om BRF back o ECI
accele a ion=qua _ o a ion(an ip0,an ip,accele a ion);
% Calcula e angle u ned
he a= o speed*s ep_sec+(1/2)*accele a ion*s ep_sec*s ep_sec; % [ ad]
% Ro a ion axis o he i e a ion
o axis= o speed/no m( o speed); % [uni a y ec o ]
% Calcula e he s ep qua e nion [To be e iewed]
q0=cos(no m( he a)/2); % Qua e nion scala pa
q= o axis*(sin(no m( he a)/2)); % Qua e nion ec o pa
68 Magne ic sys ems o he ADCS o a Fem osa elli e
% T ans o m sa elli es main axis o ECI
Sx=qua _ o a ion(q0,q,Sx); % Sa elli e X axis in ECI
Sy=qua _ o a ion(q0,q,Sy); % Sa elli e Y axis in ECI
Sz=qua _ o a ion(q0,q,Sz); % Sa elli e Z axis in ECI
% Calcula e new speed o o a ion
o speed= o speed+accele a ion*s ep_sec; % [ ad/s]
%------------------------------Sa e es da a------------------------------
% S o e o a ion e o in pe cen age
% NED o ECEF
es 1(i e a ion)=abs((no m(XYZ)-no m(ECEF))/no m(XYZ))*100;
% ECEF o ECI
es 2(i e a ion)=abs((no m(ECEF)-no m(ECI))/no m(ECEF))*100;
% ECI o BRF
es 3(i e a ion)=abs((no m(ECI)-no m(BRF))/no m(ECI))*100;
% S o e ime alue
ime(i e a ion)= /1440; % [days]
%------------------------------Upda e coun e ------------------------------
i e a ion % P in he i e a ion numbe
i e a ion=i e a ion+1;
= +s ep;
end
%--------------------------------------------------------------------------
%------------------------------Gene a e plo s------------------------------
%--------------------------------------------------------------------------
% NED o ECEF
igu e
plo ( ime, es 1)
i le('E o NED o ECEF')
xlabel('Time elapsed [days]')
ylabel('% o e o ')
% ECEF o ECI
igu e
plo ( ime, es 2)
i le('E o ECEF o ECI')
xlabel('Time elapsed [days]')
ylabel('% o e o ')
% ECI o BRF
igu e
plo ( ime, es 3)
i le('E o ECI o BRF')
xlabel('Time elapsed [days]')
ylabel('% o e o ')
Annexes 69
TLE o bi al pa ame e s code
%-------------------------------------------------------------------
%-------------- Ob ain o bi al pa ame e s o m TLE -----------------
%-------------------------------------------------------------------
unc ion [sa da a] = O bi alp( name)
ge = 398600.8; % Ea h g a i a ional cons an [km^3/s^2]
TWOPI = 2*pi;
MINUTES_PER_DAY = 1440.;
MINUTES_PER_DAY_SQUARED = (MINUTES_PER_DAY * MINUTES_PER_DAY);
MINUTES_PER_DAY_CUBED = (MINUTES_PER_DAY * MINUTES_PER_DAY_SQUARED);
% Open he TLE ile and ead TLE elemen s
id = open( name, ' ');
% ead i s line
line = ge l( id);
Cnum = line(3:7); % Ca alogue Numbe (NORAD)
SC = line(8); % Secu i y Classi ica ion
ID = line(10:17); % Iden i ica ion Numbe
yea = s 2num( line(19:20)); % Yea
doy = s 2num( line(21:32)); % Day o yea
epoch = s 2num( line(19:32)); % Epoch
TD1 = s 2num( line(34:43)); % i s ime de i a i e
TD2 = s 2num( line(45:50)); % 2nd Time De i a i e
ExTD2 = s 2num( line(51:52)); % Exponen o 2nd Time
De i a i e
BS a = s 2num( line(54:59)); % Bs a /d ag Te m
ExBS a = s 2num( line(60:61)); % Exponen o Bs a /d ag Te m
BS a = BS a *1e-5*10^ExBS a ;
E ype = line(63); % Epheme is Type
Enum = s 2num( line(65:end)); % Elemen Numbe
% ead second line
line = ge l( id);
i = s 2num( line(9:16)); % O bi Inclina ion (deg ees)
aan = s 2num( line(18:25)); % Righ Ascension o Ascending
Node (deg ees)
e = s 2num(s ca ('0.', line(27:33))); % Eccen ici y
omega = s 2num( line(35:42)); % A gumen o Pe igee (deg ees)
M = s 2num( line(44:51)); % Mean Anomaly (deg ees)
no = s 2num( line(53:63)); % Mean Mo ion
a = ( ge/(no*2*pi/86400)^2 )^(1/3); % semi majo axis (km)
No = s 2num( line(65:end)); % Re olu ion Numbe a Epoch
close( id);
sa da a.epoch = epoch;
sa da a.yea = yea ;
sa da a.doy = doy;
sa da a.no ad_numbe = Cnum;
sa da a.bulle in_numbe = ID;
sa da a.classi ica ion = SC; % almos always 'U'
sa da a. e olu ion_numbe = No;
sa da a.epheme is_ ype = E ype;
sa da a.xmo = M * (pi/180);
sa da a.xnodeo = aan * (pi/180);
sa da a.omegao = omega * (pi/180);
sa da a.xincl = i * (pi/180);
sa da a.eo = e;
sa da a.xno = no * TWOPI / MINUTES_PER_DAY;
sa da a.semimajo = a;
sa da a.xnd 2o = TD1 * TWOPI / MINUTES_PER_DAY_SQUARED;
sa da a.xndd6o = TD2 * 10^ExTD2 * TWOPI / MINUTES_PER_DAY_CUBED;
sa da a.bs a = BS a ;

70 Magne ic sys ems o he ADCS o a Fem osa elli e
Annexes 71
Qua e nion o a ion
%-------------------------------------------------------------------
%------------------------Qua e nion o a ion------------------------
%-------------------------------------------------------------------
unc ion [ ]=qua _ o a ion(p0,p, )
= + (2*p0*c oss(p, )) + c oss((2*p),(c oss(p, )));