From 43c321e36f6c5c45b879986b52840c788ca4eb2e Mon Sep 17 00:00:00 2001
From: AnkitBarik <ankit.barik@gmail.com>
Date: Thu, 18 Jan 2024 12:30:52 -0800
Subject: [PATCH] Saturn model + heat flux balance

---
 bin/autocompute.py     | 17 ++++++++++++++++-
 bin/radial_profiles.py | 24 ++++++++++++++++++------
 2 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/bin/autocompute.py b/bin/autocompute.py
index 4e1ec61..4294c5f 100644
--- a/bin/autocompute.py
+++ b/bin/autocompute.py
@@ -11,17 +11,32 @@
 import radial_profiles as rap
 
 
+def balance_heat_flux():
+
+    bc_val_bot = par.bci_thermal_val
+    bc_val_top = par.bco_thermal_val
+
+    tol = 1e-9
+
+    cd_eps = ut.chebco_f(rap.epsilon_h,par.N,par.ricb,ut.rcmb,tol)
+    epsInt = ch.chebint(cd_eps)
+
+    return eps0, bc_val_top, bc_val_bot
+
 def get_equilibrium_entropy():
     '''
     Function to compute equilibrium entropy gradient profile by solving
     Div(rho T kappa grad S) = Q
     '''
 
+    # Balance heat fluxes
+    eps0, bc_val_top, bc_val_bot = balance_heat_flux()
+
     tol = 1e-9
 
     r0  = ut.chebco(0, par.N, tol, par.ricb, ut.rcmb)
     r1  = ut.chebco(1, par.N, tol, par.ricb, ut.rcmb)
-    cd_eps = ut.chebco_f(rap.epsilon_h,par.N,par.ricb,ut.rcmb,tol)
+    cd_eps = ut.chebco_f(rap.epsilon_h,par.N,par.ricb,ut.rcmb,tol,args=eps0)
 
     D1 = ut.Dlam(1, par.N)
     D2 = ut.Dlam(2, par.N)
diff --git a/bin/radial_profiles.py b/bin/radial_profiles.py
index d35a755..596d662 100644
--- a/bin/radial_profiles.py
+++ b/bin/radial_profiles.py
@@ -93,7 +93,17 @@ def getEOSparams():
         coeffGrav = [  359.26949772994396, -1647.4951402273955, 3114.971302818801,
                       -3121.9995586185805, 1768.8806269264637, -552.869900057094,
                        78.85393983296116, 1.8991569550910268, -0.5216321175523105]
-
+    elif par.interior_model == 'saturn':
+        coeffTemp= [ 2746322.598433222, -12819667.849508610, 24913613.036177561,
+        -26159709.956658602, 16027297.768664116, -5718555.197998112,
+        1089237.928095523, -77880.376328818 ]
+        coeffDens= [ 2746322.598433222, -12819667.849508610, 24913613.036177561,
+        -26159709.956658602, 16027297.768664116, -5718555.197998112]
+        coeffGrav= [ 2746322.598433222, -12819667.849508610, 24913613.036177561,
+        1089237.928095523, -77880.376328818 ]
+        coeffAlpha= [ 2746322.598433222, -12819667.849508610, 24913613.036177561,
+        -26159709.956658602, 16027297.768664116, -5718555.197998112,
+        1089237.928095523, -77880.376328818 ]
     # elif par.interior_model == 'pns': # To be implemented
     #     pass
 
@@ -202,12 +212,14 @@ def kappa_rho(r):
     out = thermal_diffusivity(r)*density(r)
     return out
 
-def heat_source(r):
-    out = 1
+def heat_source(r,eps0=0):
+    out = eps0
     return out
 
-def epsilon_h(r):
-    out = r*heat_source(r)/(density(r)*temperature(r)*thermal_diffusivity(r))
+def epsilon_h(r,eps0=0):
+    out = (eps0*r*heat_source(r,eps0)/
+           (density(r)*temperature(r)*
+            thermal_diffusivity(r)))
     return out
 
 
@@ -286,4 +298,4 @@ def write_profiles():
     #     np.savetxt('profiles.dat',X,fmt='%.4f',header=' '.join(header),delimiter=' ')
 
     np.savetxt('profiles.dat',X,fmt='%.4f',header=' '.join(header),delimiter=' ')
-#------------------------------------------
\ No newline at end of file
+#------------------------------------------