Magnetic systems for the ADCS of a femtosatellite

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, )));