diff --git a/sub_chandra/init_1d.H b/sub_chandra/init_1d.H
index f876d84..635b049 100644
--- a/sub_chandra/init_1d.H
+++ b/sub_chandra/init_1d.H
@@ -462,7 +462,6 @@ AMREX_INLINE void init_1d() {
} // end loop over zones
-
mass_he = 0.0;
mass_wd = 0.0;
@@ -568,6 +567,11 @@ AMREX_INLINE void init_1d() {
amrex::Error("ERROR: mass did not converge");
}
+ std::cout << "converged" << std::endl;
+ std::cout << "central density of WD: " << rho_c << std::endl;
+ std::cout << "density at base of He layer: " << rho_he << std::endl;
+
+
std::cout << "final masses: " << std::endl;
std::cout << " mass WD: " << mass_wd / C::M_solar << std::endl;
std::cout << " mass He: " << mass_he / C::M_solar << std::endl;
diff --git a/toy_atm/README b/toy_atm/README
deleted file mode 100644
index 00eecac..0000000
--- a/toy_atm/README
+++ /dev/null
@@ -1,12 +0,0 @@
-This setup is used for toy_convect and for the XRB stuff. The init
-code takes an input file (one of the _param* files) to determine what
-is going on.
-
--- _params.xrb_mixed.hi_dens.tall.CNO
-
- This is the _params file that is used for the xrb_mixed stuff
-
--- _params.xrb_mixed.hi_dens.tall2.CNO
-
- This version differs from the above one with more buffer beneath
- the convective layer
\ No newline at end of file
diff --git a/toy_atm/README.md b/toy_atm/README.md
new file mode 100644
index 0000000..27c1563
--- /dev/null
+++ b/toy_atm/README.md
@@ -0,0 +1,131 @@
+# `toy_atm`
+
+Create a 1-d hydrostatic, atmosphere with an isothermal region
+(`T_star`) representing the NS, a hyperbolic tangent rise to a
+peak temperature (`T_base`) representing the base of an accreted
+layer, an isoentropic profile down to a lower temperature (`T_lo`),
+and then isothermal. This can serve as an initial model for a
+nova or XRB.
+
+The temperature profile is:
+
+```
+ ^
+ |
+ T_base + /\
+ | / \
+ | / . \
+ T_star +-----+ \
+ | . . \
+ | \
+ | . . \
+ T_lo + +-----------
+ | . .
+ +-----+---+---------------> r
+ | \ /
+ | delta
+ |< H_star>|
+```
+
+We take `dens_base`, the density at the base of the isentropic layer
+as input. The composition is "ash" in the lower isothermal region
+and "fuel" in the isentropic and upper isothermal regions. In the
+linear transition region, we linearly interpolate the composition.
+
+The fuel and ash compositions are specified by the `fuel?_name`,
+`fuel?_frac` and `ash?_name`, `ash?_frac` parameters (name of the species
+and mass fraction).
+
+The model is placed into HSE by the following differencing:
+
+```
+ (1/dr) [ _i -
_{i-1} ] = (1/2) [ _i + _{i-1} ] g
+```
+
+This will be iterated over in tandem with the EOS call,
+
+```
+ P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1)
+```
+
+## Parameters
+
+The following parameters can be set in the inputs file (prefixed with
+`problem.`).
+
+### Thermodynamics
+
+* `dens_base` : density at the base of the fuel layer
+
+* `T_star` : temperature in the isothermal region beneath the fuel layer
+
+* `T_base` : temperature at the base of the fuel layer
+
+* `T_lo` : lowest temperature in the fuel layer. Below this we switch from
+ being isentropic to isothermal
+
+* `H_star` : Height to base of the fuel layer
+
+* `delta` : with the transition region between `T_star` and `T_base`
+
+* `low_density_cutoff` : the density below which we don't worry about HSE.
+
+### Composition:
+
+The composition of the fuel is set by a pair of parameters, one giving
+the name of the species and the other giving the mass fraction. These
+all have the form:
+
+* `fuel?_name`, `fuel?_frac`, where `?` can be 1, 2, ... 7. At least the
+ first needs to be set.
+
+The ash composition is similarly set via:
+
+* `ash?_name`, `ash?_frac`, where `?` can be 1, 2, ... 7. At least the
+ first needs to be set.
+
+Finally:
+
+* `smallx` : the smallest allowed mass fraction
+
+### Domain
+
+The domain size and resolution are set via:
+
+* `xmin` : coordinate of bottom of the model
+
+* `xmax` : coordinate of the top of the model
+
+* `nx` : number of points.
+
+These should be chosen to ensure that the grid spacing, dx = (xmax - xmin) / nx
+corresponds to the finest resolution in the simulation you will be running.
+
+### Gravity
+
+There are 2 options for gravity: constant and a `1/r**2` profile:
+
+* `g_const` : the constant gravitational acceleration
+
+* `M_enclosed` : for `1/r**2` gravity, the mass enclosed beneath the
+ model. This will be used together with `xmin` to compute the
+ gravitational acceleration at the base of the model.
+
+* `do_invsq_grav` : do we use constant? or `1/r**2` gravity?
+
+
+## Science problems where this is used
+
+This setup is used for ``toy_convect`` and for the many of the XRB problems.
+The init code takes an input file:
+
+* `inputs.xrb_mixed.hi_dens.tall.CNO`
+
+ This is the `_params` file that is used for the `xrb_mixed` file
+
+* `inputs.xrb_mixed.hi_dens.tall2.CNO`
+
+ This version differs from the above one with more buffer beneath
+ the convective layer
+
+
diff --git a/toy_atm/init_1d.H b/toy_atm/init_1d.H
index 32ea245..a97248b 100755
--- a/toy_atm/init_1d.H
+++ b/toy_atm/init_1d.H
@@ -1,45 +1,3 @@
-// Create a 1-d hydrostatic, atmosphere with an isothermal region
-// (T_star) representing the NS, a hyperbolic tangent rise to a
-// peak temperature (T_base) representing the base of an accreted
-// layer, an isoentropic profile down to a lower temperature (T_lo),
-// and then isothermal. This can serve as an initial model for a
-// nova or XRB.
-//
-// The temperature profile is:
-//
-// ^
-// |
-// T_base + /\
-// | / \
-// | / . \
-// T_star +-----+ \
-// | . . \
-// | \
-// | . . \
-// T_lo + +-----------
-// | . .
-// +-----+---+---------------> r
-// | \ /
-// | delta
-// |< H_star>|
-//
-// We take dens_base, the density at the base of the isentropic layer
-// as input. The composition is "ash" in the lower isothermal region
-// and "fuel" in the isentropic and upper isothermal regions. In the
-// linear transition region, we linearly interpolate the composition.
-//
-// The fuel and ash compositions are specified by the fuel?_name,
-// fuel?_frac and ash?_name, ash?_frac parameters (name of the species
-// and mass fraction). Where ? = 1,2,3.
-//
-// The model is placed into HSE by the following differencing:
-//
-// (1/dr) [ _i -
_{i-1} ] = (1/2) [ _i + _{i-1} ] g
-//
-// This will be iterated over in tandem with the EOS call,
-// P(i-1) = P_eos(rho(i-1), T(i-1), X(i-1)
-//
-
#include
#include
diff --git a/toy_atm/inputs_subch_planar_Temp_base-1.75e8_dens_base-1.117e6_g_const-1.1742e9 b/toy_atm/inputs_subch_planar_Temp_base-1.75e8_dens_base-1.117e6_g_const-1.1742e9
new file mode 100644
index 0000000..7aed063
--- /dev/null
+++ b/toy_atm/inputs_subch_planar_Temp_base-1.75e8_dens_base-1.117e6_g_const-1.1742e9
@@ -0,0 +1,33 @@
+problem.model_prefix = "toy_subch"
+
+problem.nx = 500
+
+problem.dens_base = 1.117e6
+
+problem.T_star = 1.e7
+problem.T_base = 1.75e8
+problem.T_lo = 7.5e7
+
+problem.H_star = 8.e7
+problem.delta = 5.e6
+
+problem.fuel1_name = "helium-4"
+problem.fuel1_frac = 0.99
+
+problem.fuel2_name = "nitrogen-14"
+problem.fuel2_frac = 0.01
+
+problem.ash1_name = "carbon-12"
+problem.ash1_frac = 0.5
+
+problem.ash2_name = "oxygen-16"
+problem.ash2_frac = 0.5
+
+problem.xmin = 0.0
+problem.xmax = 1.e9
+
+problem.g_const = -1.1742e9
+
+problem.low_density_cutoff = 1.e-4
+
+problem.index_base_from_temp = 1
\ No newline at end of file