Updated tutorial

- Updated .gitignore for pycharm
- Updated the gradient_based_optimization tutorial with my understanding of what is going on with an explanation at a level that an engineer with a bachelors in comp sci understands
- Added/updated comments for the gradient_based_optimization tutorial
- Removed the duplicate copy of the gradient_based_optimization notebook
- Fixed a bug in the gradient_based_optimization code. It was trying to call tle load from kessler which is actually a local function
- Shamelessly added my name to the bottom of that tutorial so that I have something to point to for the time spend with leadership. Feel free to delete it if my explanation is terrible :-p
- Added comments to initl.py
- Removed unused import from newton_method.py
- Updated the README.md to reflect the website with the compiled docs

Signed-off-by: Grant Curell [email protected]

.gitignore

@@ -139,3 +139,5 @@ dmypy.json
 cython_debug/
 
 .vscode
+
+.idea/

README.md

@@ -25,8 +25,9 @@ Differentiable SGP4.
 
 * Differentiable version of SGP4
 
-## Installation, documentation, and examples
+## Documentation
 
+For documentation see [our GitHub pages site](https://esa.github.io/dSGP4/).
 
 ## Authors:
 * [Giacomo Acciarini](https://www.esa.int/gsp/ACT/team/giacomo_acciarini/)

doc/notebooks/code/gradient_descent_image_code.py

@@ -0,0 +1,46 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Defining a simple quadratic function f(x) = x^2
+def f(x):
+    return x ** 2
+
+# Derivative of the function f(x) = x^2
+def df(x):
+    return 2 * x
+
+# Initial parameter and learning rate
+x_start = 4.5
+learning_rate = 0.15
+# Gradient descent iterations
+iterations = 15
+
+# Storing the progression of x values for plotting
+x_progression = [x_start]
+y_progression = [f(x_start)]
+
+# Performing the gradient descent
+for _ in range(iterations):
+    x_gradient = df(x_progression[-1])
+    x_next = x_progression[-1] - learning_rate * x_gradient
+    x_progression.append(x_next)
+    y_progression.append(f(x_next))
+
+# Creating a range of x values for plotting the function
+x_values = np.linspace(-5, 5, 100)
+y_values = f(x_values)
+
+# Plotting the function f(x)
+plt.plot(x_values, y_values, label=r'$f(x) = x^2$')
+
+# Plotting the gradient descent progression
+plt.scatter(x_progression, y_progression, color='red', zorder=5)
+plt.plot(x_progression, y_progression, linestyle='--', color='red', label='Gradient Descent')
+
+# Adding details to the plot
+plt.title('2D Graph Illustrating Gradient Descent')
+plt.xlabel('x')
+plt.ylabel('f(x)')
+plt.legend()
+plt.grid(True)
+plt.show()

dsgp4/initl.py

@@ -4,59 +4,104 @@
 from . import util
 
 def initl(
-       xke, j2,
-       ecco,   epoch,  inclo,   no,
-       method,
-       opsmode,
-       ):
-
-     x2o3   = torch.tensor(2.0 / 3.0);
-
-     eccsq  = ecco * ecco;
-     omeosq = 1.0 - eccsq;
-     rteosq = omeosq.sqrt();
-     cosio  = inclo.cos();
-     cosio2 = cosio * cosio;
-
-     ak    = torch.pow(xke / no, x2o3);
-     d1    = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq);
-     del_  = d1 / (ak * ak);
-     adel  = ak * (1.0 - del_ * del_ - del_ *
-             (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0));
-     del_  = d1/(adel * adel);
-     no    = no / (1.0 + del_);
-
-     ao    = torch.pow(xke / no, x2o3);
-     sinio = inclo.sin();
-     po    = ao * omeosq;
-     con42 = 1.0 - 5.0 * cosio2;
-     con41 = -con42-cosio2-cosio2;
-     ainv  = 1.0 / ao;
-     posq  = po * po;
-     rp    = ao * (1.0 - ecco);
-     method = 'n';
-
-     if opsmode == 'a':
-         #  gst time
-         ts70  = epoch - 7305.0;
-         ds70 = torch.floor_divide(ts70 + 1.0e-8,1);
-         tfrac = ts70 - ds70;
-         #  find greenwich location at epoch
-         c1    = torch.tensor(1.72027916940703639e-2);
-         thgr70= torch.tensor(1.7321343856509374);
-         fk5r  = torch.tensor(5.07551419432269442e-15);
-         c1p2p = c1 + (2*numpy.pi);
-         gsto  = (thgr70 + c1*ds70 + c1p2p*tfrac + ts70*ts70*fk5r) % (2*numpy.pi)
-         if gsto < 0.0:
-             gsto = gsto + (2*numpy.pi);
-
-     else:
-        gsto = util.gstime(epoch + 2433281.5);
-
-     return (
-       no,
-       method,
-       ainv,  ao,    con41,  con42, cosio,
-       cosio2,eccsq, omeosq, posq,
-       rp,    rteosq,sinio , gsto,
-       )
+        xke, j2,
+        ecco, epoch, inclo, no,
+        opsmode,
+):
+    # Initialize a tensor for the fraction 2/3, used in Kepler's third law calculations.
+    # Kepler third law is: The ratio of the square of an object's orbital period with the cube of the semi-major axis
+    # of its orbit is the same for all objects orbiting the same primary.
+    x2o3 = torch.tensor(2.0 / 3.0)
+
+    # Calculate the square of the eccentricity to evaluate orbit shape and perturbation effects
+    # ecco = eccentricity
+    eccsq = ecco * ecco
+
+    # Compute one minus the square of eccentricity. This is used to describe how much an object's orbit differs from
+    # a perfect circle. If the eccentricity were zero, the value would be 1, IE: a perfect circle.
+    omeosq = 1.0 - eccsq
+
+    # Calculate the square root of (1 - eccsq), used in later orbital calculations. Ex: calculating orbital radius,
+    # and other orbital geometry
+    rteosq = omeosq.sqrt()
+
+    # Compute the cosine of the inclination, determines the object's orientation which we can use later to calculate
+    # how the object will move along the orbit
+    cosio = inclo.cos()
+
+    # Square the cosine of the inclination, used in perturbation calculations.
+    cosio2 = cosio * cosio
+
+    # Compute the semi-major axis (ak - half the diameter of an ellipse between its two furthest points) from Kepler's
+    # third law, essential for orbit scaling.
+    ak = torch.pow(xke / no, x2o3)
+
+    # Calculate the first part of the drag term, important for orbital decay predictions. The equation effectively
+    # calculates the drag induced by the "oblateness" of the earth. IE: the fact that it is flatter at the poles than
+    # in the middle.
+    # j2: represents the flattening of the Earth at the poles and its effect on the gravitational field.
+    d1 = 0.75 * j2 * (3.0 * cosio2 - 1.0) / (rteosq * omeosq)
+
+    # Estimate the perturbation factor to adjust the semi-major axis.
+    del_ = d1 / (ak * ak)
+
+    # Adjust the semi-major axis with the perturbation factor. This improves the orbit accuracy.
+    adel = ak * (1.0 - del_ * del_ - del_ *
+                 (1.0 / 3.0 + 134.0 * del_ * del_ / 81.0))
+
+    # Recalculate the perturbation factor with the adjusted semi-major axis.
+    del_ = d1 / (adel * adel)
+
+    # Correct the mean motion with the recalculated perturbation
+    no = no / (1.0 + del_)
+
+    # Recompute the semi-major axis
+    ao = torch.pow(xke / no, x2o3)
+
+    # Calculate the sine of the inclination, used in orientation and perturbation computations.
+    sinio = inclo.sin()
+
+    # Compute the orbital period
+    po = ao * omeosq
+
+    # Define constants related to the inclination, used in gravitational perturbation corrections.
+    con42 = 1.0 - 5.0 * cosio2
+    con41 = -con42 - cosio2 - cosio2
+
+    # Calculate the inverse of the semi-major axis, important for orbit adjustments.
+    ainv = 1.0 / ao
+
+    # Compute the square of the orbital period, used in timing and synchronization.
+    posq = po * po
+
+    # Determine the perigee radius, essential for closest approach calculations.
+    rp = ao * (1.0 - ecco)
+
+    # Set the method to 'n' - standard processing mode.
+    # TODO
+    method = 'n'
+
+    if opsmode == 'a':
+        # Perform calculations for the Greenwich Sidereal Time in alternative mode.
+        ts70 = epoch - 7305.0
+        ds70 = torch.floor_divide(ts70 + 1.0e-8, 1)
+        tfrac = ts70 - ds70
+        c1 = torch.tensor(1.72027916940703639e-2)
+        thgr70 = torch.tensor(1.7321343856509374)
+        fk5r = torch.tensor(5.07551419432269442e-15)
+        c1p2p = c1 + (2 * numpy.pi)
+        gsto = (thgr70 + c1 * ds70 + c1p2p * tfrac + ts70 * ts70 * fk5r) % (2 * numpy.pi)
+        if gsto < 0.0:
+            gsto += 2 * numpy.pi
+    else:
+        # Compute the Greenwich Sidereal Time for standard processing.
+        gsto = util.gstime(epoch + 2433281.5)
+
+    # Return a tuple of updated orbital elements and computed variables
+    return (
+        no,
+        method,
+        ainv, ao, con41, con42, cosio,
+        cosio2, eccsq, omeosq, posq,
+        rp, rteosq, sinio, gsto,
+    )

dsgp4/newton_method.py

@@ -1,6 +1,5 @@
 import numpy as np
 import torch
-import datetime
 from .sgp4 import sgp4
 from .sgp4init import sgp4init
 from . import util