{"id":853,"date":"2018-07-15T00:23:58","date_gmt":"2018-07-14T22:23:58","guid":{"rendered":"https:\/\/aerospaceresearch.net\/?p=853"},"modified":"2018-07-22T14:14:08","modified_gmt":"2018-07-22T12:14:08","slug":"gsoc2018orbitdeterminatorarya-week-7-8-implementing-cowell-propagator","status":"publish","type":"post","link":"https:\/\/aerospaceresearch.net\/?p=853","title":{"rendered":"[GSoC2018|OrbitDeterminator|Arya] Week #7-8: Implementing Cowell propagator"},"content":{"rendered":"<p>Last time when we tried to use the SGP4 propagator in our Kalman Filter, we ran into problems. This was because you can&#8217;t stop the algorithm at a certain stage and continue from there. Suppose you want to propagate from <img decoding=\"async\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_1\" \/> to <img decoding=\"async\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_2\" name=\"equationview\" \/>, you can directly propagate from <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_1\" name=\"equationview\" \/> to <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_2\" name=\"equationview\" \/> but you can&#8217;t propagate it from <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_1\" name=\"equationview\" \/> to <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_3\" name=\"equationview\" \/> to <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?t_2\" name=\"equationview\" \/> . So I set out to make our own propagator taking into account two major perturbations &#8211; oblateness of Earth and atmospheric drag. Instead of making an analytical propagator like SGP4, I went with a numerical one because that can be stopped and started at any point. So the first challenge was making a numerical integrator.<\/p>\n<h3>Numerical Integration<\/h3>\n<p>Numerical integration are used to integrate a differential equation from one point to another point. I studied and tested various integrators and finally settled on Runge-Kutta 4th order integrator (RK4).<\/p>\n<p>Advantages of RK4:<\/p>\n<ul>\n<li>Good accuracy for small to medium step sizes<\/li>\n<li>Low computation cost<\/li>\n<li>Easy to implement<\/li>\n<\/ul>\n<p>Let the initial value problem be <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cdot%20y%20%3D%20f%28t%2Cy%29\" name=\"equationview\" \/> and <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?y%28t_0%29%20%3D%20y_0\" name=\"equationview\" \/> .<\/p>\n<p>Choose a step size <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?h%20%3E%200\" name=\"equationview\" \/> and define<\/p>\n<p><img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cbegin%7Balign*%7D%20y_%7Bn&amp;plus;1%7D%20%26%3D%20y_n%20&amp;plus;%20%5Cfrac%7B1%7D%7B6%7D%28k_1&amp;plus;2k_2&amp;plus;2k_3&amp;plus;k_4%29%20%5C%5C%20t_%7Bn&amp;plus;1%7D%20%26%3D%20t_n%20&amp;plus;%20h%20%5Cend%7Balign*%7D\" name=\"equationview\" \/><\/p>\n<p>for\u00a0 <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?n%20%3D%200%2C1%2C2%2C3%2C%5Cdots\" name=\"equationview\" \/> where,<\/p>\n<p><img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cbegin%7Balign*%7D%20k_1%20%26%3D%20h%5C%2Cf%5Cleft%28%20t_n%2Cy_n%20%5Cright%29%20%5C%5C%20k_2%20%26%3D%20h%5C%2Cf%5Cleft%28%20t_n%20&amp;plus;%20%5Cfrac%7Bh%7D%7B2%7D%2C%20y_n%20&amp;plus;%20%5Cfrac%7Bk_1%7D%7B2%7D%20%5Cright%29%20%5C%5C%20k_3%20%26%3D%20h%5C%2Cf%5Cleft%28%20t_n%20&amp;plus;%20%5Cfrac%7Bh%7D%7B2%7D%2C%20y_n%20&amp;plus;%20%5Cfrac%7Bk_2%7D%7B2%7D%20%5Cright%29%20%5C%5C%20k_4%20%26%3D%20h%5C%2Cf%5Cleft%28%20t_n%20&amp;plus;%20h%2C%20y_n%20&amp;plus;%20k_3%20%5Cright%29%20%5C%5C%20%5Cend%7Balign*%7D\" name=\"equationview\" \/><\/p>\n<p>Note that these equations work equally well whether y is a scalar or a vector. In our case, the differential equation is <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cmathbf%7B%5Cddot%20r%7D%3D-%5Cfrac%7B%5Cmu%7D%7Br%5E3%7D%5Cvec%7B%5Cmathbf%20r%7D%20&amp;plus;%20%5Cvec%7B%5Cmathbf%20a%7D_p\" name=\"equationview\" \/>. It is not in a form suitable to be used directly in RK4. To use it in RK4 we have to vectorize the equation. First, let&#8217;s define the state of the system <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cmathbf%20s\" name=\"equationview\" \/> as <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Bx%2C%20y%2C%20z%2C%20%5Cdot%20x%2C%20%5Cdot%20y%2C%20%5Cdot%20z%5D%5ET\" name=\"equationview\" \/> . Let <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?r%20%3D%20%5Csqrt%7Bx%5E2&amp;plus;y%5E2&amp;plus;z%5E2%7D\" name=\"equationview\" \/> . Then,<\/p>\n<p><img decoding=\"async\" class=\"aligncenter\" src=\"http:\/\/quicklatex.com\/cache3\/1a\/ql_5929c6ae4452caf9c17fd49252aa7e1a_l3.png\" \/><\/p>\n<p>This form can be directly used in the RK4 integrator.<\/p>\n<p><!--more--><\/p>\n<h3>Modelling oblateness of the Earth<\/h3>\n<p>Spherical harmonics is used to model the gravity field around the Earth. The method is similar to how Taylor series works. Just like Taylor series, there are various coefficients for each degree of a term in the expansion. The most important coefficient is J<sub>2<\/sub>. The next non-zero coefficient J<sub>4<\/sub> is 1000 times smaller than J<sub>2<\/sub>. Considering only J<sub>2<\/sub>, the gravitational potential around Earth is:<\/p>\n<p><img decoding=\"async\" id=\"equationview\" class=\"aligncenter\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?U%28r%2C%5Cphi%29%20%3D%20-%5Cfrac%7B%5Cmu%7D%7Br%7D&amp;plus;%5Cfrac%7B1%7D%7B2%7D%20%5Cfrac%7B%5Cmu%7D%7Br%5E3%7D%20J_2%20R_%7B%5Coplus%7D%5E2%283%5Ccos%5E2%5Cphi%20-1%29\" name=\"equationview\" \/><\/p>\n<p>where,<br \/>\n<img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?r\" name=\"equationview\" \/> is the radius from the centre<br \/>\n<img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cphi\" name=\"equationview\" \/> is the latitude of the point with respect to the equator<br \/>\n<img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?J_2\" name=\"equationview\" \/> is a dimensionless constant equal to <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?0.001083\" name=\"equationview\" \/><br \/>\n<img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?R_%7B%5Coplus%7D\" name=\"equationview\" \/> is the equatorial radius of the Earth, equal to <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?6378.1%5C%2C%20%5Ctext%7Bkm%7D\" name=\"equationview\" \/><\/p>\n<p>To get the acceleration in xyz axes, find the gradient of U and transform coordinates to cartesian. Doing this, we get:<\/p>\n<p><img decoding=\"async\" id=\"equationview\" class=\"aligncenter\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cbegin%7Balign*%7D%20a_x%20%26%3D%20%5Cfrac%7B3%7D%7B2%7D%5Cfrac%7BJ_2%20R_%7B%5Coplus%7D%5E2%7D%7Br%5E5%7D%5Cmu%20x%5Cleft%285%5Cfrac%7Bz%5E2%7D%7Br%5E2%7D-1%5Cright%20%29%20%5C%5C%20a_y%20%26%3D%20%5Cfrac%7B3%7D%7B2%7D%5Cfrac%7BJ_2%20R_%7B%5Coplus%7D%5E2%7D%7Br%5E5%7D%5Cmu%20y%5Cleft%285%5Cfrac%7Bz%5E2%7D%7Br%5E2%7D-1%5Cright%20%29%20%5C%5C%20a_z%20%26%3D%20%5Cfrac%7B3%7D%7B2%7D%5Cfrac%7BJ_2%20R_%7B%5Coplus%7D%5E2%7D%7Br%5E5%7D%5Cmu%20z%5Cleft%285%5Cfrac%7Bz%5E2%7D%7Br%5E2%7D-3%5Cright%20%29%20%5C%5C%20%5Cend%7Balign*%7D\" name=\"equationview\" \/><\/p>\n<p>These equations implement J<sub>2<\/sub> perturbation. The main effect of J<sub>2<\/sub> perturbation is precession of the nodes and the argument of perigee.<\/p>\n<h3>Adding drag to the propagator<\/h3>\n<p>The simplest model for drag force is <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cvec%7B%5Cmathbf%20a%7D%20%3D%20-%5Cfrac%7B1%7D%7B2%7D%5Cfrac%7BC_DA%7D%7Bm%7D%5Crho%5C%2Cv_%7Brel%7D%5Cvec%7B%5Cmathbf%20v%7D_%7Brel%7D\" name=\"equationview\" \/> . We have to determine the air density, relative velocity and drag coefficients. For air density I used the model given on <a href=\"https:\/\/www.spaceacademy.net.au\/watch\/debris\/atmosmod.htm\">this<\/a> website and adjusted the constants for the current solar cycle. The final air density model is:<\/p>\n<p><img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Crho%20%3D%200.6%5Cexp%5Cleft%28-%5Cfrac%7B1%7D%7B915%7D%28h-175%29%2829.4-0.012h%29%5Cright%29%20%5Ctext%7Bkg%7D%5C%2C%20%7B%5Ctext%20km%7D%5E%7B-3%7D\" name=\"equationview\" \/><\/p>\n<p>where h is the altitude of the satellite above the surface of the Earth.<\/p>\n<p>To find the relative velocity, one must account for the rotation of the Earth. <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cvec%7B%5Cmathbf%20v%7D_%7Bsurface%7D%20%3D%20%5Cvec%7B%5Cmathbf%20%5Comega%7D%20%5Ctimes%20%5Cvec%7B%5Cmathbf%20r%7D\" name=\"equationview\" \/> where <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cvec%7B%5Cmathbf%20%5Comega%7D%20%3D%207.29115%20%5Ctimes%2010%5E%7B-5%7D%20%5C%20%5Chat%7B%5Cmathbf%20z%7D%20%5C%20%5Ctext%7Brad\/s%7D\" name=\"equationview\" \/> .<\/p>\n<p><img decoding=\"async\" id=\"equationview\" class=\"aligncenter\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cvec%7B%5Cmathbf%20v%7D_%7Brel%7D%20%3D%20%5Cvec%7B%5Cmathbf%20v%7D%20-%20%5Cvec%7B%5Cmathbf%20v%7D_%7Bsurface%7D\" name=\"equationview\" \/><\/p>\n<p>After putting together all this, I adjusted the drag coefficient so that satellites decay at approximately the current rate of decay. This finishes the implementation of drag perturbation.<\/p>\n<p>Combining both the perturbations we get: <img decoding=\"async\" id=\"equationview\" title=\"This is the rendered form of the equation. You can not edit this directly. Right click will give you the option to save the image, and in most browsers you can drag the image onto your desktop or another program.\" src=\"https:\/\/latex.codecogs.com\/svg.latex?%5Cddot%7B%5Cmathbf%20r%7D%20%3D%20-%5Cfrac%7B%5Cmu%7D%7Br%5E3%7D%5Cvec%7B%5Cmathbf%20r%7D%20&amp;plus;%20%5Cvec%7B%5Cmathbf%20a%7D_%7B%5Ctext%7Boblateness%7D%7D%20&amp;plus;%20%5Cvec%7B%5Cmathbf%20a%7D_%7B%5Ctext%7Bdrag%7D%7D\" name=\"equationview\" \/><br \/>\nPutting it in the RK4 integrator completes the propagator.<\/p>\n<figure id=\"attachment_873\" aria-describedby=\"caption-attachment-873\" style=\"width: 1050px\" class=\"wp-caption alignnone\"><img decoding=\"async\" loading=\"lazy\" class=\"size-full wp-image-873\" src=\"https:\/\/aerospaceresearch.net\/wp-content\/uploads\/2018\/07\/newplot.png\" alt=\"\" width=\"1050\" height=\"700\" \/><figcaption id=\"caption-attachment-873\" class=\"wp-caption-text\">Propagated orbit of the ISS.<\/figcaption><\/figure>\n<h3>Problems faced<\/h3>\n<p>I tried the algorithm on a recent TLE of the ISS (shown below).<\/p>\n<pre>ISS (ZARYA)             \r\n1 25544U 98067A   18190.66798432  .00001104  00000-0  24004-4 0  9994\r\n2 25544  51.6418 266.6543 0003456 290.0933 212.4518 15.54021918122018<\/pre>\n<p>According to this, it makes 15.54021 revolutions a day. This translates to a time period of 5559.8 secs. However, according to other simulators I found online, the actual period is around 5556 secs. A difference of 4 secs might not sound like much, but over the course of a day it adds up. The <a href=\"https:\/\/en.wikipedia.org\/wiki\/International_Space_Station\">Wikipedia page of the ISS<\/a> lists the time period to be 92.49 mins (5549.4 secs) but also states that it makes 15.54 orbits in a day (5559.8 secs). This is a contradiction. I searched for this issue and tried to correct it, but failed. Right now, the only method is to adjust the semi-major axis by a few kilometres to fix the time period. I will try to correct it later.<\/p>\n<h3>Future Work<\/h3>\n<ul>\n<li>Implementing the Kalman Filter<\/li>\n<li>Testing against real world satellites<\/li>\n<\/ul>\n<h3>References<\/h3>\n<ul>\n<li><a href=\"https:\/\/en.wikipedia.org\/wiki\/Runge%E2%80%93Kutta_methods\">https:\/\/en.wikipedia.org\/wiki\/Runge%E2%80%93Kutta_methods<\/a><\/li>\n<li><a href=\"https:\/\/vtechworks.lib.vt.edu\/bitstream\/handle\/10919\/49102\/Keil_EM_T_2014.pdf?sequence=1\">https:\/\/vtechworks.lib.vt.edu\/bitstream\/handle\/10919\/49102\/Keil_EM_T_2014.pdf?sequence=1<\/a><\/li>\n<li><a href=\"https:\/\/ocw.mit.edu\/courses\/earth-atmospheric-and-planetary-sciences\/12-201-essentials-of-geophysics-fall-2004\/lecture-notes\/ch2.pdf\">https:\/\/ocw.mit.edu\/courses\/earth-atmospheric-and-planetary-sciences\/12-201-essentials-of-geophysics-fall-2004\/lecture-notes\/ch2.pdf<\/a><\/li>\n<li><a href=\"https:\/\/www.spaceacademy.net.au\/watch\/debris\/atmosmod.htm\">https:\/\/www.spaceacademy.net.au\/watch\/debris\/atmosmod.htm<\/a><\/li>\n<\/ul>\n","protected":false},"excerpt":{"rendered":"<p>Last time when we tried to use the SGP4 propagator in our Kalman Filter, we ran into problems. This was because you can&#8217;t stop the algorithm at a certain stage and continue from there. Suppose you want to propagate from to , you can directly propagate from to but you can&#8217;t propagate it from to &hellip; <a href=\"https:\/\/aerospaceresearch.net\/?p=853\" class=\"more-link\"><span class=\"screen-reader-text\">\u201e[GSoC2018|OrbitDeterminator|Arya] Week #7-8: Implementing Cowell propagator\u201c<\/span> weiterlesen<\/a><\/p>\n","protected":false},"author":10,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":[],"categories":[1],"tags":[],"_links":{"self":[{"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/posts\/853"}],"collection":[{"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/users\/10"}],"replies":[{"embeddable":true,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=853"}],"version-history":[{"count":13,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/posts\/853\/revisions"}],"predecessor-version":[{"id":898,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=\/wp\/v2\/posts\/853\/revisions\/898"}],"wp:attachment":[{"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=853"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=853"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/aerospaceresearch.net\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=853"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}