Examples
========

SGTO Lifetime Analysis
----------------------

This example uses DRAMA-OSCAR and DRAMA-ARES to perform the lifetime analyses
for Super-Geostationary Transit Orbit (SGTO).

First, import the needed modules, including drama.

.. code-block:: python

    from pprint import pprint
    from datetime import datetime
    from drama import ares, oscar
    from drama.monte_carlo import Gaussian


We start with OSCAR. Which parameters can I change?

.. code-block:: python

    pprint(oscar.get_basic_config())


Let's start with a simple run, changing a few parameters.

.. code-block:: python

    pprint(oscar.run(semiMajorAxis=6678.0, epoch=datetime(2024, 3, 27)))


Luckily, we already have a :download:`project<example_projects/sgto_oscar.dpz>`
ready for the SGTO lifetime analysis. Let's check it out.

.. code-block:: python

    project = 'sgto_oscar.dpz'
    pprint(oscar.get_basic_config(project=project))


Now we are ready for our parametric analysis. We want to see how the lifetime
changes varying the initial Right Ascension of Ascending Node (RAAN) and the
initial inclination.

.. code-block:: python

    # inclination values span from 10 to 20
    inc  = [i for i in range(10,25,5)]
    pprint(inc)

    # raan values span from 100 to 130
    raan = [i for i in range(100,140,10)]
    pprint(raan)

    output = oscar.run(inc=inc, raan=raan, project=project)
    pprint(output)


We can plot the lifetime results in order to better understand the trend.

.. code-block:: python

    variables = ['inclination', 'rightAscensionOfTheAscendingNode']
    target    = 'lifetime'
    oscar.plot_heatmap(output['results'], variables=variables, target=target,
      plotfile_name=target+'_heatmap', x_label=variables[0]+' (deg)',
      y_label=variables[1]+' (deg)')


The heat map shows that the lifetime can sensibly change for similar initial RAAN
and inclination, and in general can be really difficult to be precise on the
lifetime initial conditions for the SGTO orbits, since these are Highly
Elliptical Orbits (HEO).

What if our initial RAAN and inclination are on the border between two different
lifetime areas of the heat map? In this case, We can use a Monte Carlo approach.

.. code-block:: python

    raan   = Gaussian(mu=105, sigma=1.0)
    inc    = Gaussian(mu=17, sigma=0.5)
    n_runs = 10

    output = oscar.run_monte_carlo(rightAscensionOfTheAscendingNode=raan,
                                   inclination=inc,
                                   n=n_runs,
                                   project=project)
    pprint(output)


And we want to see which result value is the  most recurring. A histogram plot
can help us with this.

.. code-block:: python

    target = 'lifetime'
    oscar.plot_histogram(output['results'], target=target,
      plotfile_name=target+'_histogram', x_label=target+' (years)')


That's a short lifetime.

We still want to check the annual collision probability for our SGTO. We can do
this easily using the DRAMA-ARES tool. We are still not sure regarding the
precise initial RAAN and inclination, but we want to understand if this
probability is low enough: let's say lower than 1.0e-6.

Also this time we have a :download:`project<example_projects/sgto_ares.dpz>`
ready for the analysis.

.. code-block:: python

    project = 'sgto_ares.dpz'
    pprint(ares.get_basic_config(project=project))


The best way is to use a Monte Carlo approach with the Wilson confidence level
(95%).

.. code-block:: python

    config = {
      'rightAscensionOfTheAscendingNode': Gaussian(mu=115, sigma=1.0),
      'inclination' : Gaussian(mu=17, sigma=0.5),
    }
    target = 'annual_collision_p'
    conf_l = 95.0
    b_crit = 1.0e-6

    output = ares.run_monte_carlo_with_wilson_confidence([config],
        project=project, confidence_level=conf_l, binary_criteria=b_crit)
    pprint(output['stopping reason'])


Let's have a look at the results: how many runs did we need to reach the
confidence level?

.. code-block:: python

    annual_collisions = [[v for k,v in res.items() if k=='annual_collision_p']
      for res in output['results']]
    pprint(len(annual_collisions))
    pprint(annual_collisions)