From 917e427b3e7c17e56905ce6549df4a76b1918dbb Mon Sep 17 00:00:00 2001
From: Zach Langford <langfzac@sas.upenn.edu>
Date: Wed, 21 Aug 2024 14:02:42 -0400
Subject: [PATCH] Fix up default ICs; add kepler-36 as available system

---
 src/ics/defaults.jl | 54 +++++++++++++++++++++++++++++----------------
 1 file changed, 35 insertions(+), 19 deletions(-)

diff --git a/src/ics/defaults.jl b/src/ics/defaults.jl
index 878cb89..c70c8c1 100644
--- a/src/ics/defaults.jl
+++ b/src/ics/defaults.jl
@@ -1,6 +1,6 @@
 # Some default initial conditions for testing and quick usage.
 
-function get_trappist_ICs(t0=0.0, n=0)
+function get_trappist1_elements(t0=0.0, n=0)
     # Trappist-1 orbital elements from Agol et al. 2021
     elements = [
         1.0 0.0 0.0 0.0 0.0 0.0 0.0
@@ -12,28 +12,32 @@ function get_trappist_ICs(t0=0.0, n=0)
         3.5331017018761430e-5 12.353988709624156 7257.667328639113  -0.0009722026578612578  0.001276000403979281  1.5707963267948966 0.0
         1.6410627049780406e-6 18.733535095576702 7250.524231929195  -0.010402303111464135  -0.014289870200773339  1.5707963267948966 0.0      
     ]
+    return elements
+end
 
-    # Use all the bodies if user did not specify n
-    nmax = 8
-    if n == 0
-        n = size(elements)[1]
-    elseif n > nmax
-        @warn "Maximum n is $(nmax). User asked for $(n). Setting to $(nmax)."
-        n = nmax
-    end
-
-    return ElementsIC(t0, n, elements)
+function get_kepler36_elements(t0=0.0, n=0)
+    # Mean Kepler-36 orbital elements from Carter et al. 2012
+    # Eccentricity chosen at random from e < 0.04, inclination co-planar
+    # Unit conversions using UnitfulAstro
+    elements = [
+        1.0 0.0 0.0 0.0 0.0 0.0 0.0
+        1.2479484222582662e-5 13.83989 0.0                0.002 -0.004 1.5707963267948966 0.0
+        2.2659378094037732e-5 16.23855 5.062100000213832 -0.01   0.007 1.5707963267948966 0.0
+    ]
+    return elements
 end
 
 # Available names for users to specify to get initial conditions
-# The first element in each tuple should be the key for IC_FUNCTIONS
+# The first element in each tuple should be the key for ELEMENTS_FUNCTIONS
 # rest are "aliases" for the system, if any.
 const AVAILABLE_SYSTEMS = (
-    ("trappist-1", "TRAPPIST-1", "Trappist-1"),
+    ("trappist-1", "trappist 1"),
+    ("kepler-36", "kepler 36")
 )
 
-const IC_FUNCTIONS = Dict(
-    "trappist-1" => get_trappist_ICs,
+const ELEMENTS_FUNCTIONS = Dict(
+    "trappist-1" => get_trappist1_elements,
+    "kepler-36" => get_kepler36_elements
 )
 
 """Return the AVAILABLE_SYSTEMS tuple. ONLY FOR TESTS."""
@@ -56,14 +60,26 @@ end
 
 """Get the default initial conditions for a particular system"""
 function get_default_ICs(system_name, t0=0.0, n=0)
-    # Get the key for the IC function to call
+    # Get the key for the elements function to call
     system_key = ""
     for available_names in AVAILABLE_SYSTEMS
-        if system_name ∈ available_names
+        if lowercase(system_name) ∈ available_names
             system_key = first(available_names)
         end
     end
+    system_key == "" && throw(ArgumentError("$(system_name) not an available system name."))
+
+    # Get the orbital elements for the selected system
+    elements = ELEMENTS_FUNCTIONS[system_key]()
 
-    func = IC_FUNCTIONS[system_key]
-    return func(t0, n)
+    # Use all the bodies if user did not specify n
+    nmax = size(elements)[1]
+    if n == 0
+        n = nmax
+    elseif n > nmax
+        @warn "Maximum n is $(nmax). User asked for $(n). Setting to $(nmax)."
+        n = nmax
+    end
+    
+    return ElementsIC(t0, n, elements)
 end