diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml
new file mode 100644
index 0000000..6b9cf15
--- /dev/null
+++ b/.github/workflows/docs.yaml
@@ -0,0 +1,29 @@
+name: Documentation
+on:
+ push:
+ branches:
+ - master
+ - main
+permissions:
+ contents: write
+jobs:
+ deploy:
+ runs-on: ubuntu-latest
+ steps:
+ - uses: actions/checkout@v4
+ - name: Configure Git Credentials
+ run: |
+ git config user.name github-actions[bot]
+ git config user.email 41898282+github-actions[bot]@users.noreply.github.com
+ - uses: actions/setup-python@v5
+ with:
+ python-version: "3.x"
+ - run: echo "cache_id=$(date --utc '+%V')" >> $GITHUB_ENV
+ - uses: actions/cache@v4
+ with:
+ key: mkdocs-material-${{ env.cache_id }}
+ path: ~/.cache
+ restore-keys: |
+ mkdocs-material-
+ - run: pip install mkdocs-material markdown-include pymdown-extensions mkdocstrings mkdocstrings-python
+ - run: mkdocs gh-deploy --force --clean
\ No newline at end of file
diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md
index 8894c10..a8500d8 100644
--- a/CONTRIBUTING.md
+++ b/CONTRIBUTING.md
@@ -14,7 +14,7 @@ mkdocs serve
You can find the documentation source in the [docs](https://github.com/FormingWorlds/MORS/tree/main/docs) directory.
If you are adding new pages, make sure to update the listing in the [`mkdocs.yml`](https://github.com/FormingWorlds/MORS/blob/main/mkdocs.yml) under the `nav` entry.
-The documentation is hosted on [readthedocs](https://fwl-mors.readthedocs.io).
+The documentation is hosted on [PROTEUS Framework](https://proteus-framework.org/MORS/).
### Running tests
MORS uses [pytest](https://docs.pytest.org/en/latest/) to run the tests. You can run the tests for yourself using:
diff --git a/README.md b/README.md
index 5ccf69d..b8ee276 100644
--- a/README.md
+++ b/README.md
@@ -1,526 +1,60 @@
-
-
# MODEL FOR ROTATION OF STARS (MORS)
-**This code is distributed as a python package for the purpose of the [PROTEUS framework](https://github.com/FormingWorlds/PROTEUS), a coupled simulation tool for the long-term evolution of atmospheres and interiors of rocky planets.
-The MORS package solves specifically the stellar rotation and evolution. It is based on the [original code](https://www.aanda.org/articles/aa/pdf/2021/05/aa38407-20.pdf) and model developed by Colin P. Johnstone.**
-
-**Original Author:** Colin P. Johnstone
-
-This code solves the stellar rotation and XUV evolution model presented in Johnstone et al. (2021). The package can be used to calculate evolutionary tracks for stellar rotaton and X-ray, EUV, and Ly-alpha emission for stars with masses between 0.1 and 1.25 Msun and has additional functionality such as allowing the user to get basic stellar parameters such as stellar radius and luminosity as functions of mass and age using the stellar evolution models of Spada et al. (2013). When publishing results that were calculated using this code, both the Johnstone et al. (2020) paper and Spada et al. (2013) should be cited.
-
-**NOTE:** This version contains the fix for the error in the equation converting EUV1 to EUV2.
-
-## 1. INSTALLATION
+MORS is a Python package (distributed as `fwl-mors`) used in the [PROTEUS framework](https://github.com/FormingWorlds/PROTEUS) to model **stellar rotation** and **high-energy emission (X-ray, EUV, Ly-α)** evolution.
+It implements the model of **Johnstone et al. (2021)** and provides stellar evolution quantities based on **Spada et al. (2013)** (plus optional Baraffe tracks).
-### 1.1. Basic install
+> **Note:** This version includes the fix for the EUV1 → EUV2 conversion.
-The Forming Worlds Mors package is available on PyPI. Run the following command to install
+## Install
-```
+```bash
pip install fwl-mors
```
-### 1.2. Developer install
-
-You can alternatively download the source code from GitHub somewhere on your computer using
-
-```
-git clone git@github.com:FormingWorlds/MORS.git
-```
-
-Then run the following command inside the main directory to install the code (check the pyproject.toml file for dependencies)
-
-```
-pip install -e .
-```
-
-
-### 1.3. Stellar evolution tracks
-
-The code requires also a set of stellar evolution data, stored in the [OSF repository](https://osf.io/9u3fb/).
-
-You can use `mors download all` to download the data. This will download and extract package stellar evolution tracks data.
-By default, MORS stores the data in based on the [XDG specification](https://specifications.freedesktop.org/basedir-spec/latest/).
-You can check the location by typing `mors env` in your terminal.
-You can override the path using the `FWL_DATA` environment variable, e.g.
+### Required data: stellar evolution tracks
+Download the stellar evolution tracks (stored on OSF):
-```console
-export FWL_DATA=...
+```bash
+mors download all
+mors env
```
-Where ... should be replaced with the path to your main data directory. To make this permanent on Ubuntu, use
+By default, data follow the XDG base directory convention. You can override the data root with:
-```console
-gedit ~/.profile
+```bash
+export FWL_DATA=/path/to/data
```
-and add the export command to the bottom of the file.
-
-Alternatively, when creating a star object in your Python script, you can specify the path to a directory where evolution tracks are stored using the starEvoDir keyword
+Or set a per-script directory:
```python
import mors
-myStar = mors.StarEvo(starEvoDir=...)
-```
-
-where ... can be given as the path relative to the current directory. When this is done, no environmental variable needs to be set.
-
-## 2. EVOLUTIONARY CALCULATIONS
-
-The main way that the user is meant to interact with the code is through the Star class. The user can create an instance of the star class, specifying only the mass and star's initial (1 Myr) rotation rate using the Mstar and Omega0 keyword arguments. This can be done for a star with a mass of 1 Msun and an initial rotation of 10x the modern Sun using
-
-```python
-star = mors.Star(Mstar=1.0, Omega=10.0)
-```
-
-The user can instead set the initial rotation rate using the Prot keyword argument giving the surface rotation period in days.
-
-```python
-star = mors.Star(Mstar=1.0, Prot=2.7)
-```
-
-Alternatively, the user can specify starting values for both the core and envelope rotation rates using the OmegaEnv and OmegaCore arguments, though this is usually not recommended.
-
-If the user instead specifies an age the Age and Omega keyword argument, the code will look for the track that passes through the specified surface rotation rate at the specified age. The surface rotation rate can be specified either using Omega or OmegaEnv.
-
-```python
-star = mors.Star(Mstar=1.0, Age=100.0, Omega=50.0)
-```
-
-In this case, the code will search for a rotation track that has a surface rotation rate that is 50x the current Sun's at an age of 100 Myr. It is not recommended to do this when the age of the star is so large that the rotation rates of stars with different initial rotation rates have converged (i.e. Age should only be specified if it is low enough that there is still a large spread in rotation rates for that mass at that age). In not all cases will it actually be able to find a physically realistic rotation rate, for example if the rotation rate specified is unreaslistically large.
-
-The code will calculate evolutionary tracks for all rotation and activity quantities available. These include surface rotation rate, OmegaEnv, and X-ray luminosity, Lx. To see a list of quantities with calculated evolutionary tracks, use the PrintAvailableTracks() attribute of the Star class.
-
-```python
-star.PrintAvailableTracks()
-```
-
-This gives units for all quantities. This can also be used to get units for each of the available quantities. Tracks for each of the quantities are held in numpy arrays and can be accessed using the Tracks dictionary, which is an attribute of the Star class. For example, the user can plot evolutionary tracks for Lx using
-
-```python
-plt.plot(star.Tracks['Age'], star.Tracks['Lx'])
-```
-
-Alternatively, numpy arrays for each quantity are attributes of the Star class with the name of the quantity followed by 'Track', so the above can be replaced with
-
-```python
-plt.plot(star.AgeTrack, star.LxTrack)
-```
-
-Units are held in the dictionary Units, which is also an attribute of the Star class. For example, to see the units of Lx, the user can use
-
-```python
-print(star.Units['Lx'])
-```
-
-The value of one of these quantities at a given age can be output using the Value() attribute of the Star class giving age the Myr and a string with the name of the desired quantity. These can be given either as positional or as keyword arguments. For example, to print Lx at 150 Myr to the screen you can use
-
-```python
-print(star.Value(150.0,'Lx'))
-```
-
-or
-
-```python
-print(star.Value(Age=150.0, Quantity='Lx'))
-```
-
-More simply, the user can use the function with the name of the desired quantity, so the two above lines could be replaced with
-
-```python
-print(star.Lx(150.0))
-```
-
-Instances of the Star class can be saved using the Save (or save) functions,
-
-```python
-star.Save()
-```
-
-By default, the data will be saved into a file called 'star.pickle' but using the user can specify the name of the file using the filename argument in the call to Save(). Saved stars can be loaded using the Load() function.
-
-```python
-star = mors.Load("star.pickle")
-```
-
-## 3. MODEL DISTRIBUTION AND PERCENTILES
-
-It is possible when creating an instance of the Star class to specify the initial rotation as a percentile of the model distribution from Johnstone et al. (2020). This model distribution is composed of measured rotation distributions from several clusters with ages of ~150 Myr evolved back to 1 Myr. To get the masses and initial rotation rates of the model cluster, use the function ModelCluster().
-
-```python
-Mstar , Omega = ModelCluster()
-```
-
-If this model cluster is used in any research, please cite the papers listed in Table 1 of Johnstone et al. (2020) for the 150 Myr age bin as the original sources for the rotation measurements (note that the function ModelCluster does not return these measurements, but 1 Myr rotation rates for these stars derived from their measurments using the rotational evolution model of Johnstone et al. 2020). To evolve this cluster, use the Cluster class as described below.
-
-When creating an instance of the Star class, use the keyword argument percentile to set the initial rotation rate to a percentile of this distribution. This can either be given as a float or int between 0 and 100 or as a string, with the three options being 'slow', 'medium', and 'fast', corresponding to the 5th, 50th, and 95th percentiles of the rotation distribution. For example
-
-```python
-star = mors.Star(Mstar=1.0, percentile=5.0)
-```
-
-This creates a solar mass star with an initial rotation rate equal to the 5th percentile of the rotation distribution at 1 Myr, as determined from our model cluster. This is equivalent to
-
-```python
-star = mors.Star(Mstar=1.0, percentile='slow')
-```
-
-To calculate this percentile, the parameter dMstarPer is used and can be set if the user sets up their own parameters as discussed above. This is set by default to 0.1, meaning that all stars within 0.1 Msun of the specified Mstar will be considered.
-
-Regardless of how a star is setup, the percentile in the model distribution is calculated after the evolutionary tracks are calculated and stored in the percentile attribute of the class. It can be seen using
-
-```python
-print(star.percentile)
-```
-
-If the user wants to find out the percentile of a given rotation race for a given mass, they can use the the Percentile function, specifying the star's mass and surface rotation rate, either as a rotation velocity in OmegaSun using the Omega (or OmegaEnv) keyword argument, or as a rotation period in days using the Prot keyword argument. For example
-
-```python
-print(mors.Percentile(Mstar=1.0, Omega=10.0))
-```
-
-This will print to the screen where in the starting (1 Myr) distribution a solar mass star rotating at 10x the current Sun is, using the model distribution from Johnstone et al. (2020). Or for a star was a 1 day rotation period
-
-```python
-print(mors.Percentile(Mstar=1.0, Prot=1.0))
-```
-
-Alternatively, the user can specify the percentile and get the corresponding rotation rate
-
-```python
-print(mors.Percentile(Mstar=1.0, percentile=10.0))
-```
-
-This prints the rotation rate of the 10th percentile for solar mass stars in OmegaSun. If the user wants to use a different cluster distribution, instead of the 1 Myr model distribution used in Johnstone et al. (2020), the masses and rotation rates of the stars in this distribution can be input. The masses are input as an array of masses in Msun using the MstarDist keyword argument, and the rotation rates can be input either as an array of surface angular velocities in OmegaSun using the OmegaDist keyword argument or as an array of rotation periods in days using the ProtDist keyword argument.
-
-## 4. SETTING SIMULATION PARAMETERS
-
-In addition to just the masses and initial rotation rates, the basic behavior of the code depends on a large number of parameters all of which have default values and do not need to be changed by the user. These default values cause the code to run the evolutionary model for rotation and XUV emission described in Johnstone et al. (2020) and generally do not need to be changed by the user. However, if the user wishes, these parameters can be changed.
-
-In general, simulation parameters are held in a dictionary called 'params' within the code, and the user can specify this parameter dictionary when creating a instance of the Star class. When this is not done, the code will use the paramsDefault dictionary created in parameters.py in the source code. To use user specified set of parameters, the user must first generate a parameter dictionary using the NewParams() function.
-
-```python
-myParams = mors.NewParams()
-```
-
-This will create a dictionary that is identical to paramsDefault that the user can edit as they want. For example, to change the parameter param1 to 1.5 and param2 to 2.5, use
-
-```python
-myParams = mors.NewParams()
-myParams['param1'] = 1.5
-myParams['param2'] = 2.5
+star_evo = mors.StarEvo(starEvoDir="path/to/evolution-tracks")
```
-Alternatively, and probably preferrably, this can also be done in the initial call to NewParams using keyword arguments.
-
-```python
-myParams = mors.NewParams(param1=1.5, param2=2.5)
-```
-
-This user should then input this as a keyword argument when creating an instance of the Star class.
-
-```python
-star = mors.Star(Mstar=1.0, Omega=10.0, params=myParams)
-```
-
-To see a complete list of all parameters that should be set, the user can use the PrintParams() function.
-
-```python
-mors.PrintParams()
-```
-
-Or can simply look into the parameters.py file in the main directory where Mors is installed.
-
-The user can set the keyword parameter AgesOut when creating a new Star object to a number or a numpy array and this will cause the code to only include these ages and the starting age in the evolutionary tracks. The ages should be given in Myr. For example
-
-```python
-star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=100.0)
-```
-
-will cause the evolutionary tracks to only contain the starting age of 1 Myr and the specified age of 100 Myr. Similarly
-
-```python
-star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=np.array([100.0,200.0,300.0,400.0]))
-```
-
-will cause the evolutionary tracks to contain 1 Myr, and the specified 100, 200, 300, and 400 Myr ages. The simulations will also end at the oldest year in the array. This array should be in ascending order. This is not recommended if the user wants to extract quantities at arbitrary ages along evolutionary tracks, for example using star.Lx(Age), since these functions interpolate between the values and if there are not enough age bins in the evolutionary tracks then these results can be inaccurate. Note that having two many age bins in the AgesOut array, e.g. with AgesOut=np.linspace(1.0,5000.0,10000), will cause the calculations of the evolutionary tracks to be very slow.
-
-## 5. ROTATION AND ACTIVITY QUANTITES
-
-The code gives the user the ability to use the basic functions for calculating many activity related quantities as functions os stellar properties. For example, the user can calculate XUV luminosities from stellar mass, age, and rotation rates using the Lxuv function. For example,
-
-```python
-Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
-```
-
-This gives a dictionary holding X-ray, EUV, and Lyman-alpha luminosities for a 5000 Myr old solar mass star rotating 10 times faster than the Sun. The surface rotation rate can be specified as a rotation velocity using the Omega or OmegaEnv arguments or as a rotation period in days using the Prot argument, e.g.
-
-```python
-Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Prot=1.0)
-```
-
-This is similar to above, but for a star with a rotation period of 1 day.
-
-The dictionary returned contains the following luminosities, all in erg s-1
-
-- Lxuv - 0.517 to 92 nm
-- Lx - 0.517 to 12.4 nm (0.1 to 24 keV)
-- Leuv1 - 10 to 36 nm
-- Leuv2 - 36 to 92 nm
-- Leuv - 10 to 92 nm
-- Lly - Lymann-alpha emission line
-
-The dictionary also contains surface fluxes (e.g. Fxuv, Fx, Feuv1,...) in erg s-1 cm-2 and luminosities normalised to bolometric luminosities (e.g Rxuv, Rx, Reuv1,...).
-
-The user can also just get the X-ray luminosity as a float using Lx()
-
-```python
-Lx = mors.Lx(Mstar=1.0, Age=5000.0, Omega=10.0)
-```
-
-And similarly for the EUV and Lymann-alpha luminosities
-
-```python
-Leuv = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0)
-Lly = mors.Lly(Mstar=1.0, Age=5000.0, Omega=10.0)
-```
-
-The function Leuv() also takes the keyword argument 'band' which can be used to specify the band desired (=0 for 10-92 nm; =1 for 10-32 nm; =2 for 32-92 nm).
-
-The above function calculate an average Lx for a star given its mass, age, and rotation. In reality there is a scatter around these values that appear somewhat random, likely due to variability. This scatter can be described as a log normal probability density function. To calculate this scatter for X-ray luminosity, the user can use the XrayScatter() function.
-
-```python
-deltaLx = mors.XrayScatter(Lx)
-```
-
-Here, the user should just enter the X-ray luminosity either as a single value or a numpy array and it returns the deltaLx from the equation LxScattered = Lx + deltaLx, where Lx is the value with this scatter not included and LxScattered is the value with it included. The user can also send in Fx or Rx to the function and get deltaFx and deltaRx. Note that these values are random and will be different each time. The width of the random distribution is set by sigmaXray in the parameter file and can be set using the params keyword argument.
-
-To get scatter values for EUV, Ly-alpha, and all other parameters returned by Lxuv() described above, use XUVScatter() which takes the dictionary returned by Lxuv().
-
-```python
-Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
-deltaLxuv = mors.XUVScatter(Lxuv)
-```
-
-The return value is a dictionary containing delta values for the same quantities in Lxuv.
-
-If the user wants a more detailed set of parameters for a given star, the ExtendedQuantities() function can be used and in this case, the star's mass, age, and envelope and core rotation rates must be specified.
-
-```python
-quantities = ExtendedQuantities(Mstar=1.0, Age=5000.0, OmegaEnv=10.0, OmegaCore=10.0)
-```
-
-This returns a larger dictionary with many other quantities, including some basic stellar properties such as radius and moment of inertia, the mass loss rate in the wind, the dipole field strength, and all of the torques acting on the star to influence its rotation. A list of the available parameters can be seen with
-
-```python
-list(quantities)
-```
-
-## 6. STELLAR EVOLUTION QUANTITIES
-
-The rotation and activity evolution model requires that various basic stellar properties such as bolometric luminosity, radius, and convective turnover time, can be accessed at all ages for for all masses within the mass range considered (0.1 to 1.25 Msun). These are calculated using the evolutionary tracks fronm Spada et al. (2013) and the functions that do this are available to the user. First import the Mors package
+## Quickstart
```python
import mors
-```
-
-The stellar mass for a 1.0 Msun star with an age of 1000 Myr can be calculated using
-
-```python
-Rstar = mors.Rstar(1.0, 1000.0)
-```
-
-The functions available are
-
-- Rstar - radius (Rsun)
-- Lbol - bolometric luminosity (Lsun)
-- Teff - effective temperature (K)
-- Itotal - total moment of inertia in (g cm2)
-- Icore - core moment of inertia in (g cm2)
-- Ienv - envelope moment of inertia in (g cm2)
-- Mcore - core mass in (Msun)
-- Menv - envelope mass in (Msun)
-- Rcore - core radius in (Rsun)
-- tau - convective turnover time (days)
-- dItotaldt - rate of change of total moment of inertia (g cm2 Myr-1)
-- dIcoredt - rate of change of core moment of inertia (g cm2 Myr-1)
-- dIenvdt - rate of change of envelope moment of inertia (g cm2 Myr-1)
-- dMcoredt - rate of change of core mass (Msun Myr-1)
-- dRcoredt - rate of change of core radius (Rsun Myr^-1)
+import matplotlib.pyplot as plt
-Note that core and envelope have the definitions from the rotation model, so the 'core' is not the hydrogen burning core as it would typically be defined but everything interior to the outer convective zone, and the envelope is the outer convective zone. All of these functions take stellar mass in Msun and age in Myr. Alternatively, the user can call the function Value which takes the same two inputs and then the parameter name as a string. For example
-
-```python
-Rstar = mors.Value(1.0, 1000.0, 'Rstar')
-```
-
-which is equivalent to the above example.
-
-In all cases, multiple values for mass and age can be input using lists or numpy arrays and the functions will return numpy arrays holding the results for each input case. For example, if stellar mass is input as a 1D array or list of length 5, the result will be a length 5 numpy array. If both stellar mass and age are input as arrays, the result will be a 2-dimensional array of length len(Mstar)xlen(Age). Additionally, in the call to Value(), the string holding the parameter to calculate for can also be input as a list of strings, in which case the variable retuned will have an additional dimension with values for each of the input quantities.
-
-The first time one of these functions is called, the code loads all of the evolutionary tracks from the Spada et al. (2013) model. These will be written to the directors in the file SEmodels.pickle in the main directory of the evolutionary models and this file will be used to load the models in the future. This file can be safely deleted if the user wants since this will just cause the code to remake it next time it is needed. The code then does a linear interpolation between mass and age bins to get values at the input masses and ages. This can be time consuming, especially if an interpolation between mass bins is required, which will be the case if the input age is not an integer multiple of 0.05 Msun. The user can load a full evolutionary track for a specific stellar mass using LoadTrack(). For example, to load a track for a 1.02 Msun star, use
-
-```python
-mors.LoadTrack(1.02)
-```
+star = mors.Star(Mstar=1.0, Prot=2.7) # 1 Msun star, initial rotation period in days (at age ~1 Myr)
+print(star.Lx(150.0)) # X-ray luminosity at 150 Myr
-If a track for that specific mass is already loaded, this will do nothing.
-
-## 7. CLUSTER EVOLUTION CALCULATIONS
-
-The code allows the user to calculate the evolution of stellar clusters in addition to single stars. This is done using the Cluster class. An instance of the Cluster class can be created in much the same way as an instance of the Star class using the same input keyword arguments. The only difference is that the masses and rotation rates of the stars should be given as arrays or lists
-
-```python
-Mstar = np.array([1.0, 0.5, 0.75])
-Omega = np.array([ 10.0, 10.0, 10.0])
-cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
-```
-
-This will create a cluster composed of three stars and immediately calculate evolutionary tracks for each. To see the progress of these calculations printed to the screen, use the verbose keyword argument
-
-```python
-cluster = mors.Cluster(Mstar=Mstar, Omega=Omega, verbose=True)
-```
-
-As in the Star class, it is possible to specify the Age keyword argument. If this is not specified, the code will assume the Omega values are starting (1 Myr) rotation rates. If it is specified, the code will fit the rotation tracks such that they pass through the specified rotation rates at the speficied ages. If Age is specified as a single value (in Myr), it will be assumed that all stars have that age, and if it is specified as a numpy array then the array should have the same lengths as Mstar and Omega and the ages will be assumed individually for each star. As with the Star class, it is possible to input simulation parameters using the params keyword argument (see above).
-
-As in the Star class, an instance of the Cluster class can be saved using the Save (or save) function
-
-```python
-cluster.Save()
-```
-
-By default, the data will be saved into a file called 'cluster.pickle' but using the user can specify the name o the file using the filename argument in the call to Save(). Saved clusters can be loaded using the Load() function.
-
-```python
-cluster = mors.Load("cluster.pickle")
-```
-
-This is advisable since the evolutionary tracks for clusters can take a lot of time to calculate if there are lots of stars.
-
-The Star class instances for each star in the cluster is held in the stars list, which is an attribute of this class. To plot Lx tracks for each star for example, the user can use
-
-```python
-plt.plot(cluster.stars[0].AgeTrack, cluster.stars[0].LxTrack)
-plt.plot(cluster.stars[1].AgeTrack, cluster.stars[1].LxTrack)
-plt.plot(cluster.stars[2].AgeTrack, cluster.stars[2].LxTrack)
-```
-
-Alternatively, each star is an attribute of the class instance and has the name 'star' followed by its place in the list (starting at zero), so the above can be replaced with
-
-```python
-plt.plot(cluster.stars0.AgeTrack, cluster.stars0.LxTrack)
-plt.plot(cluster.stars1.AgeTrack, cluster.stars1.LxTrack)
-plt.plot(cluster.stars2.AgeTrack, cluster.stars2.LxTrack)
-```
-
-To get values for each star at a given age of a given parameter, the user can use the Values() function, specifying the age in Myr and a string for the quantity to retrieve. For example,
-
-```python
-Lx = cluster.Values(Age=100.0, Quantity='Lx')
-```
-
-This gives an array with Lx for each star in the cluster. Alternatively, as with the Star class, each quantity can be specified using its own function and the above is equivalent to.
-
-```python
-Lx = cluster.Lx(100.0)
-```
-
-So to plot the Lx distribution for a cluster at 100 Myr, the user can use
-
-```python
-plt.scatter(cluster.Mstar, cluster.Lx(100.0))
-```
-
-In Johnstone et al. (2020), many of the results are based on a composite cluster that is often referred to in the paper as the model cluster. This is composed of measured rotation distributions from several clusters with ages of ~150 Myr evolved back to 1 Myr.
-To get the masses and initial rotation rates of the model cluster, use the function ModelCluster().
-
-```python
-Mstar , Omega = ModelCluster()
-```
-
-If this model cluster is used in any research, please cite the papers listed in Table 1 of Johnstone et al. (2020) for the 150 Myr age bin as the original sources for the rotation measurements (note that the function ModelCluster does not return these measurements, but 1 Myr rotation rates for these stars derived from their measurments using the rotational evolution model of Johnstone et al. 2020). This cluster can then be evolved in the usual way
-
-```python
-Mstar , Omega = ModelCluster()
-cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
-```
-
-The Cluster class also has a function that allows the user to find where in the distribution a star with a given mass and surface rotation rate rate is at a given age. This uses the Percentile function discussed above applied to the rotation distribution of this cluster at the specified age.
-
-## 8. HABITABLE ZONE BOUNDARIES
-
-Using the formulae of Kopparapu et al. (2013) and the luminosities and effective temperatures from the stellar models of Spada et al. (2013), the user can calculate the orbital distances of the habitable zone boundaries as a function of stellar mass and age. Please cite these two papers if using the output of this function. The HZ boundaries are calculated using the aOrbHZ function.
-
-```python
-aOrbHZ = mors.aOrbHZ(Mstar=1.0)
-```
-
-In this case, the function returns a dictionary holding the orbital distances in AU in a dictionary. The six values in the dictionary are 'RecentVenus', 'RunawayGreenhouse', 'MoistGreenhouse', 'MaximumGreenhouse', 'EarlyMars', and 'HZ'. While the first five are obvious and correspond to the boundaries in Kopparapu et al. (2013), the final one is defined in Johnstone et al. (2020) as the average of the runaway greenhouse and moist greenhouse orbital distances. For example, to see these values
-
-```python
-print(aOrbHZ['RecentVenus'])
-print(aOrbHZ['RunawayGreenhouse'])
-print(aOrbHZ['MoistGreenhouse'])
-print(aOrbHZ['MaximumGreenhouse'])
-print(aOrbHZ['EarlyMars'])
-print(aOrbHZ['HZ'])
-```
-
-Mstar can be input as a numpy array in which case the dictionary will contain arrays with values for each mass. By default the orbital distances are calculated assuming stellar parameters at 5000 Myr. The user can also set the age in Myr using the Age keyword argument
-
-```python
-aOrbHZ = mors.aOrbHZ(Mstar=1.0, Age=1000.0)
-```
-
-## 9. ACTIVITY LIFETIMES
-
-The code contains functions that calculate how long stars remain above certain acitivty thresholds. Firstly, the function ActivityLifetime() takes an evolutionary track and returns when the value of the track goes below a certain threshold. For example
-
-```python
-AgeThreshold = mors.ActivityLifetime(Age=star.AgeTrack, Track=star.LxTrack, Threshold=1.0e28)
-```
-
-This will tell us when the star's Lx dropped below 1028 erg s-1. If the star crosses this threshold multiple times, only the final time will be returned. If it never goes below this threshold then the final age in the track will be returned and if it is never above this threshold then 0.0 will be returned. If the user wants to set a maximum age so that the code only looks for crossings of the threshold below this age then this can be done using the AgeMax keyword argument
-
-More recommended though is to use the function of the same name in the Star class, where the quantity of interest should be specified as a string.
-
-```python
-AgeThreshold = star.ActivityLifetime(Quantity='Lx', Threshold=1.0e28)
-```
-
-This does the same thing. Valid options for Quantity are 'Lx', 'Fx', 'Rx', 'FxHZ', 'Leuv1', 'Feuv1', 'Reuv1', 'Feuv1HZ', 'Leuv2', 'Feuv2', 'Reuv2', 'Feuv2HZ', 'Leuv', 'Feuv', 'Reuv', 'FeuvHZ', 'Lly', 'Fly', 'Rly', and 'FlyHZ'. It is also possible here to specify AgeMax. If the user wants to know how long the star remained saturated, Threshold can be given with the string 'sat'.
-
-```python
-AgeSaturated = star.ActivityLifetime(Quantity='Lx', Threshold='sat')
+plt.plot(star.AgeTrack, star.LxTrack)
+plt.xlabel(f"Age [{star.Units['Age']}]")
+plt.ylabel(f"Lx [{star.Units['Lx']}]")
+plt.show()
```
-Note that when doing this, the value of Quantity should not influence the return value since all XUV quantities saturate at the same time, though it still must be set to a valid option.
+## Documentation
-The Cluster class also contains this function and is called in the exact same way, returning this time the values for each star in the cluster as a numpy array.
+You can find the complete documentation [here](https://proteus-framework.org/MORS/).
-```python
-AgeThreshold = cluster.ActivityLifetime(Quantity='Lx', Threshold=1.0e28)
-```
-## 10. BARAFFE STELLAR EVOLUTION QUANTITIES
+## Citation
-The Forming Worlds MORS package also provides access to stellar evolution quantities according to the Baraffe model (Baraffe et al., 2002). The evolution track for a given star mass Mstar, between 0.01 and 1.4 Msun, can be loaded into memory with the command
+When publishing results computed with MORS, please cite:
+- **Johnstone et al. (2021)** for the rotation/XUV evolution model
+- **Spada et al. (2013)** for the stellar evolution tracks used for stellar properties
-```python
-baraffe = mors.BaraffeTrack(Mstar)
-```
-The command performs mass interpolation (if needed) and time interpolation into a fine grid of 5e4 points. For a given time, the user can then access to the luminosity, the solar constant and the stellar radius with the following functions
-
-- Stellar radius in Rsun, time in yr
-```
-Rstar = baraffe.BaraffeStellarRadius(time)
-```
-- Bolometric flux in Lsun, time in yr
-```
-Lbol = baraffe.BaraffeLuminosity(time)
-```
-- Flux scaled by the star-planet distance in W m-2, time in yr, distance in AU
-```
-Flux = baraffe.BaraffeSolarConstant(time, distance)
-```
+If you use the model cluster distribution/percentiles, also cite the rotation-measurement sources referenced in **Johnstone et al. (2020)** (Table 1, ~150 Myr bin).
diff --git a/README_old.md b/README_old.md
new file mode 100644
index 0000000..5ccf69d
--- /dev/null
+++ b/README_old.md
@@ -0,0 +1,526 @@
+
+
+# MODEL FOR ROTATION OF STARS (MORS)
+
+**This code is distributed as a python package for the purpose of the [PROTEUS framework](https://github.com/FormingWorlds/PROTEUS), a coupled simulation tool for the long-term evolution of atmospheres and interiors of rocky planets.
+The MORS package solves specifically the stellar rotation and evolution. It is based on the [original code](https://www.aanda.org/articles/aa/pdf/2021/05/aa38407-20.pdf) and model developed by Colin P. Johnstone.**
+
+**Original Author:** Colin P. Johnstone
+
+This code solves the stellar rotation and XUV evolution model presented in Johnstone et al. (2021). The package can be used to calculate evolutionary tracks for stellar rotaton and X-ray, EUV, and Ly-alpha emission for stars with masses between 0.1 and 1.25 Msun and has additional functionality such as allowing the user to get basic stellar parameters such as stellar radius and luminosity as functions of mass and age using the stellar evolution models of Spada et al. (2013). When publishing results that were calculated using this code, both the Johnstone et al. (2020) paper and Spada et al. (2013) should be cited.
+
+**NOTE:** This version contains the fix for the error in the equation converting EUV1 to EUV2.
+
+## 1. INSTALLATION
+
+### 1.1. Basic install
+
+The Forming Worlds Mors package is available on PyPI. Run the following command to install
+
+```
+pip install fwl-mors
+```
+### 1.2. Developer install
+
+You can alternatively download the source code from GitHub somewhere on your computer using
+
+```
+git clone git@github.com:FormingWorlds/MORS.git
+```
+
+Then run the following command inside the main directory to install the code (check the pyproject.toml file for dependencies)
+
+```
+pip install -e .
+```
+
+
+### 1.3. Stellar evolution tracks
+
+The code requires also a set of stellar evolution data, stored in the [OSF repository](https://osf.io/9u3fb/).
+
+You can use `mors download all` to download the data. This will download and extract package stellar evolution tracks data.
+
+By default, MORS stores the data in based on the [XDG specification](https://specifications.freedesktop.org/basedir-spec/latest/).
+You can check the location by typing `mors env` in your terminal.
+You can override the path using the `FWL_DATA` environment variable, e.g.
+
+```console
+export FWL_DATA=...
+```
+
+Where ... should be replaced with the path to your main data directory. To make this permanent on Ubuntu, use
+
+```console
+gedit ~/.profile
+```
+
+and add the export command to the bottom of the file.
+
+Alternatively, when creating a star object in your Python script, you can specify the path to a directory where evolution tracks are stored using the starEvoDir keyword
+
+```python
+import mors
+myStar = mors.StarEvo(starEvoDir=...)
+```
+
+where ... can be given as the path relative to the current directory. When this is done, no environmental variable needs to be set.
+
+## 2. EVOLUTIONARY CALCULATIONS
+
+The main way that the user is meant to interact with the code is through the Star class. The user can create an instance of the star class, specifying only the mass and star's initial (1 Myr) rotation rate using the Mstar and Omega0 keyword arguments. This can be done for a star with a mass of 1 Msun and an initial rotation of 10x the modern Sun using
+
+```python
+star = mors.Star(Mstar=1.0, Omega=10.0)
+```
+
+The user can instead set the initial rotation rate using the Prot keyword argument giving the surface rotation period in days.
+
+```python
+star = mors.Star(Mstar=1.0, Prot=2.7)
+```
+
+Alternatively, the user can specify starting values for both the core and envelope rotation rates using the OmegaEnv and OmegaCore arguments, though this is usually not recommended.
+
+If the user instead specifies an age the Age and Omega keyword argument, the code will look for the track that passes through the specified surface rotation rate at the specified age. The surface rotation rate can be specified either using Omega or OmegaEnv.
+
+```python
+star = mors.Star(Mstar=1.0, Age=100.0, Omega=50.0)
+```
+
+In this case, the code will search for a rotation track that has a surface rotation rate that is 50x the current Sun's at an age of 100 Myr. It is not recommended to do this when the age of the star is so large that the rotation rates of stars with different initial rotation rates have converged (i.e. Age should only be specified if it is low enough that there is still a large spread in rotation rates for that mass at that age). In not all cases will it actually be able to find a physically realistic rotation rate, for example if the rotation rate specified is unreaslistically large.
+
+The code will calculate evolutionary tracks for all rotation and activity quantities available. These include surface rotation rate, OmegaEnv, and X-ray luminosity, Lx. To see a list of quantities with calculated evolutionary tracks, use the PrintAvailableTracks() attribute of the Star class.
+
+```python
+star.PrintAvailableTracks()
+```
+
+This gives units for all quantities. This can also be used to get units for each of the available quantities. Tracks for each of the quantities are held in numpy arrays and can be accessed using the Tracks dictionary, which is an attribute of the Star class. For example, the user can plot evolutionary tracks for Lx using
+
+```python
+plt.plot(star.Tracks['Age'], star.Tracks['Lx'])
+```
+
+Alternatively, numpy arrays for each quantity are attributes of the Star class with the name of the quantity followed by 'Track', so the above can be replaced with
+
+```python
+plt.plot(star.AgeTrack, star.LxTrack)
+```
+
+Units are held in the dictionary Units, which is also an attribute of the Star class. For example, to see the units of Lx, the user can use
+
+```python
+print(star.Units['Lx'])
+```
+
+The value of one of these quantities at a given age can be output using the Value() attribute of the Star class giving age the Myr and a string with the name of the desired quantity. These can be given either as positional or as keyword arguments. For example, to print Lx at 150 Myr to the screen you can use
+
+```python
+print(star.Value(150.0,'Lx'))
+```
+
+or
+
+```python
+print(star.Value(Age=150.0, Quantity='Lx'))
+```
+
+More simply, the user can use the function with the name of the desired quantity, so the two above lines could be replaced with
+
+```python
+print(star.Lx(150.0))
+```
+
+Instances of the Star class can be saved using the Save (or save) functions,
+
+```python
+star.Save()
+```
+
+By default, the data will be saved into a file called 'star.pickle' but using the user can specify the name of the file using the filename argument in the call to Save(). Saved stars can be loaded using the Load() function.
+
+```python
+star = mors.Load("star.pickle")
+```
+
+## 3. MODEL DISTRIBUTION AND PERCENTILES
+
+It is possible when creating an instance of the Star class to specify the initial rotation as a percentile of the model distribution from Johnstone et al. (2020). This model distribution is composed of measured rotation distributions from several clusters with ages of ~150 Myr evolved back to 1 Myr. To get the masses and initial rotation rates of the model cluster, use the function ModelCluster().
+
+```python
+Mstar , Omega = ModelCluster()
+```
+
+If this model cluster is used in any research, please cite the papers listed in Table 1 of Johnstone et al. (2020) for the 150 Myr age bin as the original sources for the rotation measurements (note that the function ModelCluster does not return these measurements, but 1 Myr rotation rates for these stars derived from their measurments using the rotational evolution model of Johnstone et al. 2020). To evolve this cluster, use the Cluster class as described below.
+
+When creating an instance of the Star class, use the keyword argument percentile to set the initial rotation rate to a percentile of this distribution. This can either be given as a float or int between 0 and 100 or as a string, with the three options being 'slow', 'medium', and 'fast', corresponding to the 5th, 50th, and 95th percentiles of the rotation distribution. For example
+
+```python
+star = mors.Star(Mstar=1.0, percentile=5.0)
+```
+
+This creates a solar mass star with an initial rotation rate equal to the 5th percentile of the rotation distribution at 1 Myr, as determined from our model cluster. This is equivalent to
+
+```python
+star = mors.Star(Mstar=1.0, percentile='slow')
+```
+
+To calculate this percentile, the parameter dMstarPer is used and can be set if the user sets up their own parameters as discussed above. This is set by default to 0.1, meaning that all stars within 0.1 Msun of the specified Mstar will be considered.
+
+Regardless of how a star is setup, the percentile in the model distribution is calculated after the evolutionary tracks are calculated and stored in the percentile attribute of the class. It can be seen using
+
+```python
+print(star.percentile)
+```
+
+If the user wants to find out the percentile of a given rotation race for a given mass, they can use the the Percentile function, specifying the star's mass and surface rotation rate, either as a rotation velocity in OmegaSun using the Omega (or OmegaEnv) keyword argument, or as a rotation period in days using the Prot keyword argument. For example
+
+```python
+print(mors.Percentile(Mstar=1.0, Omega=10.0))
+```
+
+This will print to the screen where in the starting (1 Myr) distribution a solar mass star rotating at 10x the current Sun is, using the model distribution from Johnstone et al. (2020). Or for a star was a 1 day rotation period
+
+```python
+print(mors.Percentile(Mstar=1.0, Prot=1.0))
+```
+
+Alternatively, the user can specify the percentile and get the corresponding rotation rate
+
+```python
+print(mors.Percentile(Mstar=1.0, percentile=10.0))
+```
+
+This prints the rotation rate of the 10th percentile for solar mass stars in OmegaSun. If the user wants to use a different cluster distribution, instead of the 1 Myr model distribution used in Johnstone et al. (2020), the masses and rotation rates of the stars in this distribution can be input. The masses are input as an array of masses in Msun using the MstarDist keyword argument, and the rotation rates can be input either as an array of surface angular velocities in OmegaSun using the OmegaDist keyword argument or as an array of rotation periods in days using the ProtDist keyword argument.
+
+## 4. SETTING SIMULATION PARAMETERS
+
+In addition to just the masses and initial rotation rates, the basic behavior of the code depends on a large number of parameters all of which have default values and do not need to be changed by the user. These default values cause the code to run the evolutionary model for rotation and XUV emission described in Johnstone et al. (2020) and generally do not need to be changed by the user. However, if the user wishes, these parameters can be changed.
+
+In general, simulation parameters are held in a dictionary called 'params' within the code, and the user can specify this parameter dictionary when creating a instance of the Star class. When this is not done, the code will use the paramsDefault dictionary created in parameters.py in the source code. To use user specified set of parameters, the user must first generate a parameter dictionary using the NewParams() function.
+
+```python
+myParams = mors.NewParams()
+```
+
+This will create a dictionary that is identical to paramsDefault that the user can edit as they want. For example, to change the parameter param1 to 1.5 and param2 to 2.5, use
+
+```python
+myParams = mors.NewParams()
+myParams['param1'] = 1.5
+myParams['param2'] = 2.5
+```
+
+Alternatively, and probably preferrably, this can also be done in the initial call to NewParams using keyword arguments.
+
+```python
+myParams = mors.NewParams(param1=1.5, param2=2.5)
+```
+
+This user should then input this as a keyword argument when creating an instance of the Star class.
+
+```python
+star = mors.Star(Mstar=1.0, Omega=10.0, params=myParams)
+```
+
+To see a complete list of all parameters that should be set, the user can use the PrintParams() function.
+
+```python
+mors.PrintParams()
+```
+
+Or can simply look into the parameters.py file in the main directory where Mors is installed.
+
+The user can set the keyword parameter AgesOut when creating a new Star object to a number or a numpy array and this will cause the code to only include these ages and the starting age in the evolutionary tracks. The ages should be given in Myr. For example
+
+```python
+star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=100.0)
+```
+
+will cause the evolutionary tracks to only contain the starting age of 1 Myr and the specified age of 100 Myr. Similarly
+
+```python
+star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=np.array([100.0,200.0,300.0,400.0]))
+```
+
+will cause the evolutionary tracks to contain 1 Myr, and the specified 100, 200, 300, and 400 Myr ages. The simulations will also end at the oldest year in the array. This array should be in ascending order. This is not recommended if the user wants to extract quantities at arbitrary ages along evolutionary tracks, for example using star.Lx(Age), since these functions interpolate between the values and if there are not enough age bins in the evolutionary tracks then these results can be inaccurate. Note that having two many age bins in the AgesOut array, e.g. with AgesOut=np.linspace(1.0,5000.0,10000), will cause the calculations of the evolutionary tracks to be very slow.
+
+## 5. ROTATION AND ACTIVITY QUANTITES
+
+The code gives the user the ability to use the basic functions for calculating many activity related quantities as functions os stellar properties. For example, the user can calculate XUV luminosities from stellar mass, age, and rotation rates using the Lxuv function. For example,
+
+```python
+Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+```
+
+This gives a dictionary holding X-ray, EUV, and Lyman-alpha luminosities for a 5000 Myr old solar mass star rotating 10 times faster than the Sun. The surface rotation rate can be specified as a rotation velocity using the Omega or OmegaEnv arguments or as a rotation period in days using the Prot argument, e.g.
+
+```python
+Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Prot=1.0)
+```
+
+This is similar to above, but for a star with a rotation period of 1 day.
+
+The dictionary returned contains the following luminosities, all in erg s-1
+
+- Lxuv - 0.517 to 92 nm
+- Lx - 0.517 to 12.4 nm (0.1 to 24 keV)
+- Leuv1 - 10 to 36 nm
+- Leuv2 - 36 to 92 nm
+- Leuv - 10 to 92 nm
+- Lly - Lymann-alpha emission line
+
+The dictionary also contains surface fluxes (e.g. Fxuv, Fx, Feuv1,...) in erg s-1 cm-2 and luminosities normalised to bolometric luminosities (e.g Rxuv, Rx, Reuv1,...).
+
+The user can also just get the X-ray luminosity as a float using Lx()
+
+```python
+Lx = mors.Lx(Mstar=1.0, Age=5000.0, Omega=10.0)
+```
+
+And similarly for the EUV and Lymann-alpha luminosities
+
+```python
+Leuv = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+Lly = mors.Lly(Mstar=1.0, Age=5000.0, Omega=10.0)
+```
+
+The function Leuv() also takes the keyword argument 'band' which can be used to specify the band desired (=0 for 10-92 nm; =1 for 10-32 nm; =2 for 32-92 nm).
+
+The above function calculate an average Lx for a star given its mass, age, and rotation. In reality there is a scatter around these values that appear somewhat random, likely due to variability. This scatter can be described as a log normal probability density function. To calculate this scatter for X-ray luminosity, the user can use the XrayScatter() function.
+
+```python
+deltaLx = mors.XrayScatter(Lx)
+```
+
+Here, the user should just enter the X-ray luminosity either as a single value or a numpy array and it returns the deltaLx from the equation LxScattered = Lx + deltaLx, where Lx is the value with this scatter not included and LxScattered is the value with it included. The user can also send in Fx or Rx to the function and get deltaFx and deltaRx. Note that these values are random and will be different each time. The width of the random distribution is set by sigmaXray in the parameter file and can be set using the params keyword argument.
+
+To get scatter values for EUV, Ly-alpha, and all other parameters returned by Lxuv() described above, use XUVScatter() which takes the dictionary returned by Lxuv().
+
+```python
+Lxuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+deltaLxuv = mors.XUVScatter(Lxuv)
+```
+
+The return value is a dictionary containing delta values for the same quantities in Lxuv.
+
+If the user wants a more detailed set of parameters for a given star, the ExtendedQuantities() function can be used and in this case, the star's mass, age, and envelope and core rotation rates must be specified.
+
+```python
+quantities = ExtendedQuantities(Mstar=1.0, Age=5000.0, OmegaEnv=10.0, OmegaCore=10.0)
+```
+
+This returns a larger dictionary with many other quantities, including some basic stellar properties such as radius and moment of inertia, the mass loss rate in the wind, the dipole field strength, and all of the torques acting on the star to influence its rotation. A list of the available parameters can be seen with
+
+```python
+list(quantities)
+```
+
+## 6. STELLAR EVOLUTION QUANTITIES
+
+The rotation and activity evolution model requires that various basic stellar properties such as bolometric luminosity, radius, and convective turnover time, can be accessed at all ages for for all masses within the mass range considered (0.1 to 1.25 Msun). These are calculated using the evolutionary tracks fronm Spada et al. (2013) and the functions that do this are available to the user. First import the Mors package
+
+```python
+import mors
+```
+
+The stellar mass for a 1.0 Msun star with an age of 1000 Myr can be calculated using
+
+```python
+Rstar = mors.Rstar(1.0, 1000.0)
+```
+
+The functions available are
+
+- Rstar - radius (Rsun)
+- Lbol - bolometric luminosity (Lsun)
+- Teff - effective temperature (K)
+- Itotal - total moment of inertia in (g cm2)
+- Icore - core moment of inertia in (g cm2)
+- Ienv - envelope moment of inertia in (g cm2)
+- Mcore - core mass in (Msun)
+- Menv - envelope mass in (Msun)
+- Rcore - core radius in (Rsun)
+- tau - convective turnover time (days)
+- dItotaldt - rate of change of total moment of inertia (g cm2 Myr-1)
+- dIcoredt - rate of change of core moment of inertia (g cm2 Myr-1)
+- dIenvdt - rate of change of envelope moment of inertia (g cm2 Myr-1)
+- dMcoredt - rate of change of core mass (Msun Myr-1)
+- dRcoredt - rate of change of core radius (Rsun Myr^-1)
+
+Note that core and envelope have the definitions from the rotation model, so the 'core' is not the hydrogen burning core as it would typically be defined but everything interior to the outer convective zone, and the envelope is the outer convective zone. All of these functions take stellar mass in Msun and age in Myr. Alternatively, the user can call the function Value which takes the same two inputs and then the parameter name as a string. For example
+
+```python
+Rstar = mors.Value(1.0, 1000.0, 'Rstar')
+```
+
+which is equivalent to the above example.
+
+In all cases, multiple values for mass and age can be input using lists or numpy arrays and the functions will return numpy arrays holding the results for each input case. For example, if stellar mass is input as a 1D array or list of length 5, the result will be a length 5 numpy array. If both stellar mass and age are input as arrays, the result will be a 2-dimensional array of length len(Mstar)xlen(Age). Additionally, in the call to Value(), the string holding the parameter to calculate for can also be input as a list of strings, in which case the variable retuned will have an additional dimension with values for each of the input quantities.
+
+The first time one of these functions is called, the code loads all of the evolutionary tracks from the Spada et al. (2013) model. These will be written to the directors in the file SEmodels.pickle in the main directory of the evolutionary models and this file will be used to load the models in the future. This file can be safely deleted if the user wants since this will just cause the code to remake it next time it is needed. The code then does a linear interpolation between mass and age bins to get values at the input masses and ages. This can be time consuming, especially if an interpolation between mass bins is required, which will be the case if the input age is not an integer multiple of 0.05 Msun. The user can load a full evolutionary track for a specific stellar mass using LoadTrack(). For example, to load a track for a 1.02 Msun star, use
+
+```python
+mors.LoadTrack(1.02)
+```
+
+If a track for that specific mass is already loaded, this will do nothing.
+
+## 7. CLUSTER EVOLUTION CALCULATIONS
+
+The code allows the user to calculate the evolution of stellar clusters in addition to single stars. This is done using the Cluster class. An instance of the Cluster class can be created in much the same way as an instance of the Star class using the same input keyword arguments. The only difference is that the masses and rotation rates of the stars should be given as arrays or lists
+
+```python
+Mstar = np.array([1.0, 0.5, 0.75])
+Omega = np.array([ 10.0, 10.0, 10.0])
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
+```
+
+This will create a cluster composed of three stars and immediately calculate evolutionary tracks for each. To see the progress of these calculations printed to the screen, use the verbose keyword argument
+
+```python
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega, verbose=True)
+```
+
+As in the Star class, it is possible to specify the Age keyword argument. If this is not specified, the code will assume the Omega values are starting (1 Myr) rotation rates. If it is specified, the code will fit the rotation tracks such that they pass through the specified rotation rates at the speficied ages. If Age is specified as a single value (in Myr), it will be assumed that all stars have that age, and if it is specified as a numpy array then the array should have the same lengths as Mstar and Omega and the ages will be assumed individually for each star. As with the Star class, it is possible to input simulation parameters using the params keyword argument (see above).
+
+As in the Star class, an instance of the Cluster class can be saved using the Save (or save) function
+
+```python
+cluster.Save()
+```
+
+By default, the data will be saved into a file called 'cluster.pickle' but using the user can specify the name o the file using the filename argument in the call to Save(). Saved clusters can be loaded using the Load() function.
+
+```python
+cluster = mors.Load("cluster.pickle")
+```
+
+This is advisable since the evolutionary tracks for clusters can take a lot of time to calculate if there are lots of stars.
+
+The Star class instances for each star in the cluster is held in the stars list, which is an attribute of this class. To plot Lx tracks for each star for example, the user can use
+
+```python
+plt.plot(cluster.stars[0].AgeTrack, cluster.stars[0].LxTrack)
+plt.plot(cluster.stars[1].AgeTrack, cluster.stars[1].LxTrack)
+plt.plot(cluster.stars[2].AgeTrack, cluster.stars[2].LxTrack)
+```
+
+Alternatively, each star is an attribute of the class instance and has the name 'star' followed by its place in the list (starting at zero), so the above can be replaced with
+
+```python
+plt.plot(cluster.stars0.AgeTrack, cluster.stars0.LxTrack)
+plt.plot(cluster.stars1.AgeTrack, cluster.stars1.LxTrack)
+plt.plot(cluster.stars2.AgeTrack, cluster.stars2.LxTrack)
+```
+
+To get values for each star at a given age of a given parameter, the user can use the Values() function, specifying the age in Myr and a string for the quantity to retrieve. For example,
+
+```python
+Lx = cluster.Values(Age=100.0, Quantity='Lx')
+```
+
+This gives an array with Lx for each star in the cluster. Alternatively, as with the Star class, each quantity can be specified using its own function and the above is equivalent to.
+
+```python
+Lx = cluster.Lx(100.0)
+```
+
+So to plot the Lx distribution for a cluster at 100 Myr, the user can use
+
+```python
+plt.scatter(cluster.Mstar, cluster.Lx(100.0))
+```
+
+In Johnstone et al. (2020), many of the results are based on a composite cluster that is often referred to in the paper as the model cluster. This is composed of measured rotation distributions from several clusters with ages of ~150 Myr evolved back to 1 Myr.
+To get the masses and initial rotation rates of the model cluster, use the function ModelCluster().
+
+```python
+Mstar , Omega = ModelCluster()
+```
+
+If this model cluster is used in any research, please cite the papers listed in Table 1 of Johnstone et al. (2020) for the 150 Myr age bin as the original sources for the rotation measurements (note that the function ModelCluster does not return these measurements, but 1 Myr rotation rates for these stars derived from their measurments using the rotational evolution model of Johnstone et al. 2020). This cluster can then be evolved in the usual way
+
+```python
+Mstar , Omega = ModelCluster()
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
+```
+
+The Cluster class also has a function that allows the user to find where in the distribution a star with a given mass and surface rotation rate rate is at a given age. This uses the Percentile function discussed above applied to the rotation distribution of this cluster at the specified age.
+
+## 8. HABITABLE ZONE BOUNDARIES
+
+Using the formulae of Kopparapu et al. (2013) and the luminosities and effective temperatures from the stellar models of Spada et al. (2013), the user can calculate the orbital distances of the habitable zone boundaries as a function of stellar mass and age. Please cite these two papers if using the output of this function. The HZ boundaries are calculated using the aOrbHZ function.
+
+```python
+aOrbHZ = mors.aOrbHZ(Mstar=1.0)
+```
+
+In this case, the function returns a dictionary holding the orbital distances in AU in a dictionary. The six values in the dictionary are 'RecentVenus', 'RunawayGreenhouse', 'MoistGreenhouse', 'MaximumGreenhouse', 'EarlyMars', and 'HZ'. While the first five are obvious and correspond to the boundaries in Kopparapu et al. (2013), the final one is defined in Johnstone et al. (2020) as the average of the runaway greenhouse and moist greenhouse orbital distances. For example, to see these values
+
+```python
+print(aOrbHZ['RecentVenus'])
+print(aOrbHZ['RunawayGreenhouse'])
+print(aOrbHZ['MoistGreenhouse'])
+print(aOrbHZ['MaximumGreenhouse'])
+print(aOrbHZ['EarlyMars'])
+print(aOrbHZ['HZ'])
+```
+
+Mstar can be input as a numpy array in which case the dictionary will contain arrays with values for each mass. By default the orbital distances are calculated assuming stellar parameters at 5000 Myr. The user can also set the age in Myr using the Age keyword argument
+
+```python
+aOrbHZ = mors.aOrbHZ(Mstar=1.0, Age=1000.0)
+```
+
+## 9. ACTIVITY LIFETIMES
+
+The code contains functions that calculate how long stars remain above certain acitivty thresholds. Firstly, the function ActivityLifetime() takes an evolutionary track and returns when the value of the track goes below a certain threshold. For example
+
+```python
+AgeThreshold = mors.ActivityLifetime(Age=star.AgeTrack, Track=star.LxTrack, Threshold=1.0e28)
+```
+
+This will tell us when the star's Lx dropped below 1028 erg s-1. If the star crosses this threshold multiple times, only the final time will be returned. If it never goes below this threshold then the final age in the track will be returned and if it is never above this threshold then 0.0 will be returned. If the user wants to set a maximum age so that the code only looks for crossings of the threshold below this age then this can be done using the AgeMax keyword argument
+
+More recommended though is to use the function of the same name in the Star class, where the quantity of interest should be specified as a string.
+
+```python
+AgeThreshold = star.ActivityLifetime(Quantity='Lx', Threshold=1.0e28)
+```
+
+This does the same thing. Valid options for Quantity are 'Lx', 'Fx', 'Rx', 'FxHZ', 'Leuv1', 'Feuv1', 'Reuv1', 'Feuv1HZ', 'Leuv2', 'Feuv2', 'Reuv2', 'Feuv2HZ', 'Leuv', 'Feuv', 'Reuv', 'FeuvHZ', 'Lly', 'Fly', 'Rly', and 'FlyHZ'. It is also possible here to specify AgeMax. If the user wants to know how long the star remained saturated, Threshold can be given with the string 'sat'.
+
+```python
+AgeSaturated = star.ActivityLifetime(Quantity='Lx', Threshold='sat')
+```
+
+Note that when doing this, the value of Quantity should not influence the return value since all XUV quantities saturate at the same time, though it still must be set to a valid option.
+
+The Cluster class also contains this function and is called in the exact same way, returning this time the values for each star in the cluster as a numpy array.
+
+```python
+AgeThreshold = cluster.ActivityLifetime(Quantity='Lx', Threshold=1.0e28)
+```
+## 10. BARAFFE STELLAR EVOLUTION QUANTITIES
+
+The Forming Worlds MORS package also provides access to stellar evolution quantities according to the Baraffe model (Baraffe et al., 2002). The evolution track for a given star mass Mstar, between 0.01 and 1.4 Msun, can be loaded into memory with the command
+
+```python
+baraffe = mors.BaraffeTrack(Mstar)
+```
+The command performs mass interpolation (if needed) and time interpolation into a fine grid of 5e4 points. For a given time, the user can then access to the luminosity, the solar constant and the stellar radius with the following functions
+
+- Stellar radius in Rsun, time in yr
+```
+Rstar = baraffe.BaraffeStellarRadius(time)
+```
+- Bolometric flux in Lsun, time in yr
+```
+Lbol = baraffe.BaraffeLuminosity(time)
+```
+- Flux scaled by the star-planet distance in W m-2, time in yr, distance in AU
+```
+Flux = baraffe.BaraffeSolarConstant(time, distance)
+```
diff --git a/docs/How-to/activity_quantities.md b/docs/How-to/activity_quantities.md
new file mode 100644
index 0000000..3b70f88
--- /dev/null
+++ b/docs/How-to/activity_quantities.md
@@ -0,0 +1,129 @@
+## Rotation and activity quantities (How-to)
+
+### Goal
+Compute high-energy emission quantities (X-ray, EUV, Ly-α) from stellar **mass**, **age**, and **rotation**, optionally add variability/scatter, and (if needed) retrieve a larger set of model diagnostics via `ExtendedQuantities`.
+
+### Prerequisites
+Install MORS and download the required stellar evolution data:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Mstar` in **M☉**, `Age` in **Myr**, `Prot` in **days**, `Omega` in **Ω☉**. Luminosities are in **erg s⁻¹** and surface fluxes in **erg s⁻¹ cm⁻²**.
+
+---
+
+### Step 1: Get the full XUV dictionary (`Lxuv`)
+
+Use `Omega` (or `OmegaEnv`) **or** `Prot` to specify the surface rotation.
+
+**Using Ω (Ω☉):**
+```python
+import mors
+xuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+```
+
+**Using rotation period (days):**
+```python
+import mors
+xuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Prot=1.0)
+```
+
+The returned dictionary includes luminosities, surface fluxes, and bolometric-normalized values, e.g.:
+- Luminosities: `Lxuv`, `Lx`, `Leuv1`, `Leuv2`, `Leuv`, `Lly` (erg s⁻¹)
+- Fluxes: `Fxuv`, `Fx`, `Feuv1`, ... (erg s⁻¹ cm⁻²)
+- Normalized: `Rxuv`, `Rx`, `Reuv1`, ... (dimensionless)
+
+To see what keys are present:
+```python
+print(sorted(xuv.keys()))
+```
+
+---
+
+### Step 2: Retrieve a single quantity directly (`Lx`, `Leuv`, `Lly`)
+
+```python
+import mors
+
+Lx = mors.Lx(Mstar=1.0, Age=5000.0, Omega=10.0)
+Leu = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+Lly = mors.Lly(Mstar=1.0, Age=5000.0, Omega=10.0)
+```
+
+#### Choose an EUV sub-band with `band`
+`Leuv(..., band=...)` supports:
+- `band=0` → 10–92 nm (total EUV)
+- `band=1` → 10–36 nm (EUV1)
+- `band=2` → 36–92 nm (EUV2)
+
+```python
+import mors
+Leuv_total = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0, band=0)
+Leuv1 = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0, band=1)
+Leuv2 = mors.Leuv(Mstar=1.0, Age=5000.0, Omega=10.0, band=2)
+```
+
+---
+
+### Step 3: Add variability/scatter (optional)
+
+MORS provides random scatter terms intended to represent variability as a log-normal scatter.
+
+#### X-ray scatter (`XrayScatter`)
+You pass the mean value (e.g., `Lx`) and get a delta to add:
+```python
+import mors
+
+Lx = mors.Lx(Mstar=1.0, Age=5000.0, Omega=10.0)
+dLx = mors.XrayScatter(Lx)
+
+Lx_scattered = Lx + dLx
+```
+
+You can also pass `Fx` or `Rx` and receive `dFx` / `dRx`.
+
+#### Full XUV scatter (`XUVScatter`)
+`XUVScatter` takes the dictionary from `Lxuv` and returns a dictionary of deltas with the same keys:
+```python
+import mors
+
+xuv = mors.Lxuv(Mstar=1.0, Age=5000.0, Omega=10.0)
+dxuv = mors.XUVScatter(xuv)
+
+# Example: apply deltas to get one scattered realization
+xuv_scattered = {k: xuv[k] + dxuv[k] for k in xuv}
+```
+
+**Controlling the scatter width:** the X-ray scatter width is controlled by `sigmaXray` in the parameters file. You can override it by passing a custom `params` dictionary (see the parameter how-to):
+
+```python
+import mors
+
+params = mors.NewParams(sigmaXray=0.3) # example value
+Lx = mors.Lx(Mstar=1.0, Age=5000.0, Omega=10.0, params=params)
+dLx = mors.XrayScatter(Lx, params=params) if "params" in mors.XrayScatter.__code__.co_varnames else mors.XrayScatter(Lx)
+```
+*(If `XrayScatter` does not accept `params` directly, set `sigmaXray` via your model run parameters and use it consistently.)*
+
+---
+
+### Step 4: Get detailed diagnostics (`ExtendedQuantities`)
+
+If you need additional model internals (stellar properties, wind quantities, magnetic fields, torques), call `ExtendedQuantities` with envelope and core rotation rates:
+
+```python
+import mors
+
+q = mors.ExtendedQuantities(Mstar=1.0, Age=5000.0, OmegaEnv=10.0, OmegaCore=10.0)
+print(list(q))
+```
+
+---
+
+### Common pitfalls
+- Don’t mix `Omega` (Ω☉) and `Prot` (days) in the same call; pick one rotation representation.
+- Scatter functions are random; results differ each run unless you control the random seed in your workflow.
+- `ExtendedQuantities` requires **both** `OmegaEnv` and `OmegaCore`.
diff --git a/docs/How-to/activity_timelines.md b/docs/How-to/activity_timelines.md
new file mode 100644
index 0000000..d0f19ea
--- /dev/null
+++ b/docs/How-to/activity_timelines.md
@@ -0,0 +1,95 @@
+## Activity timelines (How-to)
+
+### Goal
+Compute how long a star (or a cluster of stars) stays above an activity threshold (e.g., X-ray luminosity) or how long it remains in the saturated regime.
+
+### Prerequisites
+Install MORS and download stellar evolution data:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** thresholds must match the quantity you use (e.g., `Lx` in **erg s⁻¹**). Ages returned are in **Myr**.
+
+---
+
+### Compute activity threshold
+
+#### Option 1: Use the standalone helper on an explicit track (basic)
+
+`mors.ActivityLifetime` takes an age array and a track array and returns the age when the track finally drops below the threshold.
+
+```python
+import mors
+
+AgeThreshold = mors.ActivityLifetime(
+ Age=star.AgeTrack,
+ Track=star.LxTrack,
+ Threshold=1.0e28
+)
+print(AgeThreshold)
+```
+
+**Behavior:**
+- If the track crosses the threshold multiple times → returns the **final** crossing time.
+- If it never goes below the threshold → returns the **final age** in the track.
+- If it is never above the threshold → returns `0.0`.
+
+Optional: limit the search to ages below `AgeMax`.
+
+---
+
+#### Option 2: Prefer the `Star.ActivityLifetime` method (recommended)
+
+This version selects the track internally; you only provide the quantity name and threshold.
+
+```python
+AgeThreshold = star.ActivityLifetime(Quantity="Lx", Threshold=1.0e28)
+print(AgeThreshold)
+```
+
+You can also set a maximum age to consider:
+
+```python
+AgeThreshold = star.ActivityLifetime(Quantity="Lx", Threshold=1.0e28, AgeMax=1000.0)
+```
+
+**Note**: The `Quantity` string must match one of the supported activity tracks:
+
+- `Lx`, `Fx`, `Rx`, `FxHZ`
+- `Leuv1`, `Feuv1`, `Reuv1`, `Feuv1HZ`
+- `Leuv2`, `Feuv2`, `Reuv2`, `Feuv2HZ`
+- `Leuv`, `Feuv`, `Reuv`, `FeuvHZ`
+- `Lly`, `Fly`, `Rly`, `FlyHZ`
+
+---
+
+### Compute cluster lifetimes (same interface, returns arrays)
+
+For a `mors.Cluster`, the call is the same but returns one value per star:
+
+```python
+AgeThreshold = cluster.ActivityLifetime(Quantity="Lx", Threshold=1.0e28)
+print(AgeThreshold) # numpy array
+```
+
+---
+
+### Compute saturation lifetime
+
+To get the time the star remains in the saturated regime, pass `Threshold="sat"`:
+
+```python
+AgeSaturated = star.ActivityLifetime(Quantity="Lx", Threshold="sat")
+print(AgeSaturated)
+```
+
+Note: the choice of `Quantity` should not change the saturation lifetime (XUV quantities saturate together), but it must be a valid option.
+
+---
+
+### Common pitfalls
+- Make sure your `Threshold` uses the same units as the chosen `Quantity` (e.g., `Lx` is erg s⁻¹).
+- The function returns the **final** threshold crossing time if there are multiple crossings.
diff --git a/docs/How-to/clusters.md b/docs/How-to/clusters.md
new file mode 100644
index 0000000..657cf0b
--- /dev/null
+++ b/docs/How-to/clusters.md
@@ -0,0 +1,145 @@
+## Cluster evolution calculation (How-to)
+
+### Goal
+Compute rotation/activity evolution tracks for a **cluster of stars**, then (a) inspect per-star tracks, (b) evaluate cluster distributions at a given age, and (c) save/reload the cluster to avoid recomputation.
+
+### Prerequisites
+Install MORS and download stellar evolution data:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Mstar` in **M☉**, `Age` in **Myr**, `Prot` in **days**, `Omega` in **Ω☉**.
+
+---
+
+### Step 1 — Create a cluster from arrays
+
+Create arrays/lists of masses and rotation rates (same length). If `Age` is omitted, MORS interprets `Omega` as the **initial (~1 Myr)** rotation rate.
+
+```python
+import numpy as np
+import mors
+
+Mstar = np.array([1.0, 0.5, 0.75])
+Omega = np.array([10.0, 10.0, 10.0])
+
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
+```
+
+To show progress while tracks are computed, enable `verbose`:
+
+```python
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega, verbose=True)
+```
+
+---
+
+### Step 2 — Fit tracks through a specified age (optional)
+
+If you provide `Age`, MORS fits tracks so that each star passes through the given rotation rate at that age.
+
+**One age applied to all stars**
+```python
+import numpy as np
+import mors
+
+Mstar = np.array([1.0, 0.5, 0.75])
+Omega = np.array([50.0, 30.0, 40.0])
+
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega, Age=100.0) # Myr
+```
+
+**Different ages for each star**
+```python
+import numpy as np
+import mors
+
+Mstar = np.array([1.0, 0.5, 0.75])
+Omega = np.array([50.0, 30.0, 40.0])
+Age = np.array([80.0, 120.0, 100.0]) # Myr; same length as Mstar/Omega
+
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega, Age=Age)
+```
+
+---
+
+### Step 3 — Access per-star tracks
+
+Each star is a `mors.Star` instance stored in `cluster.stars`.
+
+```python
+import matplotlib.pyplot as plt
+
+plt.plot(cluster.stars[0].AgeTrack, cluster.stars[0].LxTrack)
+plt.plot(cluster.stars[1].AgeTrack, cluster.stars[1].LxTrack)
+plt.plot(cluster.stars[2].AgeTrack, cluster.stars[2].LxTrack)
+plt.show()
+```
+
+Some versions also expose `stars0`, `stars1`, ... as attributes:
+
+```python
+plt.plot(cluster.stars0.AgeTrack, cluster.stars0.LxTrack)
+```
+
+---
+
+### Step 4 — Get cluster values at a fixed age
+
+Use `Values(Age=..., Quantity=...)` to retrieve an array across the cluster:
+
+```python
+Lx = cluster.Values(Age=100.0, Quantity="Lx")
+```
+
+Or use the convenience method named after the quantity:
+
+```python
+Lx = cluster.Lx(100.0)
+```
+
+Plot a distribution (e.g., Lx vs mass) at a given age:
+
+```python
+import matplotlib.pyplot as plt
+plt.scatter(cluster.Mstar, cluster.Lx(100.0))
+plt.xlabel("Mstar [Msun]")
+plt.ylabel("Lx [erg/s]")
+plt.show()
+```
+
+---
+
+### Step 5 — Save and reload (recommended for large clusters)
+
+Cluster calculations can be expensive for many stars, so saving is recommended.
+
+```python
+cluster.Save(filename="cluster.pickle")
+cluster2 = mors.Load("cluster.pickle")
+```
+
+---
+
+### Step 6 — Evolve the built-in “model cluster” (optional)
+
+MORS includes a composite “model cluster” distribution (derived from observed clusters at ~150 Myr evolved back to 1 Myr). You can load it and evolve it like any other cluster:
+
+```python
+import mors
+
+Mstar, Omega = mors.ModelCluster()
+cluster = mors.Cluster(Mstar=Mstar, Omega=Omega)
+```
+
+**Citation note:** If you use this model cluster distribution in research, cite the rotation-measurement sources listed in **Johnstone et al. (2020), Table 1** (150 Myr bin), in addition to the MORS model paper(s).
+
+---
+
+### Common pitfalls
+- Arrays for `Mstar`, `Omega` (and `Age` if provided) must have matching lengths.
+- `Age` is in **Myr**; `Omega` is in **Ω☉**.
+- Fitting tracks using `Age` + `Omega` can be unreliable at late ages where tracks converge.
diff --git a/docs/How-to/distribution_percentile.md b/docs/How-to/distribution_percentile.md
new file mode 100644
index 0000000..7b806af
--- /dev/null
+++ b/docs/How-to/distribution_percentile.md
@@ -0,0 +1,133 @@
+# Model distribution and percentiles (How-to)
+
+### Goal
+Pick an initial rotation rate from the built-in **model rotation distribution** (Johnstone et al. 2020), inspect a star’s inferred percentile, and convert between **rotation rate and percentile** for a given stellar mass.
+
+### Prerequisites
+Make sure MORS and the stellar evolution data are installed:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Omega` is in units of the current solar rotation rate (**Ω☉**). `Prot` is in **days**. `Mstar` is in **M☉**.
+
+---
+
+### Step 1: Get the built-in “model cluster” distribution
+
+This returns arrays of stellar masses and their **1 Myr** rotation rates (derived by evolving observed cluster distributions back to 1 Myr).
+
+```python
+import mors
+
+Mstar_dist, Omega_dist = mors.ModelCluster()
+print(Mstar_dist.shape, Omega_dist.shape)
+```
+
+**Citation note:** If you use the model cluster distribution directly in research, cite the rotation-measurement sources listed in **Johnstone et al. (2020), Table 1** (150 Myr bin), in addition to the MORS model paper(s).
+
+---
+
+### Step 2: Create a star by percentile
+
+You can set initial rotation using a numeric percentile (0–100) or the strings `'slow'`, `'medium'`, `'fast'` which map to 5th/50th/95th percentiles.
+
+```python
+import mors
+
+# numeric percentile
+star_p5 = mors.Star(Mstar=1.0, percentile=5.0)
+
+# equivalent string shortcut
+star_slow = mors.Star(Mstar=1.0, percentile="slow")
+```
+
+---
+
+### Step 3: Read back the star’s inferred percentile
+
+Regardless of how you create the star (Omega/Prot/percentile), MORS stores the inferred percentile after computing tracks:
+
+```python
+print(star_slow.percentile)
+```
+
+---
+
+### Step 4: Convert a rotation rate (Ω or Prot) to a percentile
+
+**From Ω (Ω☉) to percentile**
+```python
+import mors
+p = mors.Percentile(Mstar=1.0, Omega=10.0)
+print(p)
+```
+
+**From Prot (days) to percentile**
+```python
+import mors
+p = mors.Percentile(Mstar=1.0, Prot=1.0)
+print(p)
+```
+
+---
+
+### Step 5: Convert a percentile to a rotation rate
+
+```python
+import mors
+Omega_p10 = mors.Percentile(Mstar=1.0, percentile=10.0) # returns Ω in Ω☉
+print(Omega_p10)
+```
+
+---
+
+### Step 6: Control the mass window used for percentiles (`dMstarPer`)
+
+Percentiles are computed using stars within a mass window around `Mstar`. The width is controlled by `dMstarPer` (default: 0.1 M☉). To change it, modify parameters and pass them into `Star`:
+
+```python
+import mors
+
+params = mors.NewParams(dMstarPer=0.05)
+star = mors.Star(Mstar=1.0, percentile="medium", params=params)
+```
+
+---
+
+### Step 7: Use a custom rotation distribution (optional)
+
+If you have your own mass+rotation distribution (instead of the built-in 1 Myr model distribution), pass it into `Percentile`.
+
+**Using Ω distribution**
+```python
+import numpy as np
+import mors
+
+MstarDist = np.array([1.0, 1.0, 1.0, 0.9, 1.1])
+OmegaDist = np.array([2.0, 5.0, 10.0, 3.0, 7.0])
+
+p = mors.Percentile(Mstar=1.0, Omega=6.0, MstarDist=MstarDist, OmegaDist=OmegaDist)
+print(p)
+```
+
+**Using Prot distribution**
+```python
+import numpy as np
+import mors
+
+MstarDist = np.array([1.0, 1.0, 0.9, 1.1])
+ProtDist = np.array([12.0, 6.0, 9.0, 7.0]) # days
+
+p = mors.Percentile(Mstar=1.0, Prot=8.0, MstarDist=MstarDist, ProtDist=ProtDist)
+print(p)
+```
+
+---
+
+### Common pitfalls
+- Percentiles here refer to the **initial (~1 Myr)** distribution used by the model cluster.
+- `Omega` is **Ω☉**, not rad/s. `Prot` is **days**.
+- If your custom distribution is sparse or covers a narrow mass range, percentile estimates can become unstable.
diff --git a/docs/How-to/evolution.md b/docs/How-to/evolution.md
new file mode 100644
index 0000000..b1980de
--- /dev/null
+++ b/docs/How-to/evolution.md
@@ -0,0 +1,96 @@
+## Evolutionary calculations (How-to)
+
+### Goal
+Compute a star’s rotation and activity evolution tracks, then (a) plot a quantity, (b) query values at specific ages, and (c) save the result for reuse.
+
+### Prerequisites
+Make sure the package and stellar evolution data are installed:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Age` is in **Myr**, `Prot` is in **days**, and `Omega` is in units of the **current solar rotation rate (Ω☉)**.
+
+---
+
+### Step 1: Create a `Star`
+
+Pick **one** of these common ways:
+
+**A) Set initial rotation using Ω (at ~1 Myr)**
+
+```python
+import mors
+star = mors.Star(Mstar=1.0, Omega=10.0)
+```
+
+**B) Set initial rotation using rotation period (days)**
+
+```python
+import mors
+star = mors.Star(Mstar=1.0, Prot=2.7)
+```
+
+**C) Fit a track through a point (Age, Ω)**
+
+```python
+import mors
+star = mors.Star(Mstar=1.0, Age=100.0, Omega=50.0)
+```
+
+Use (C) only at ages where rotation tracks have not strongly converged for that mass; otherwise the fit may be ill-posed or fail if the requested rotation rate is unrealistic.
+
+---
+
+### Step 2: Inspect available tracks and units
+
+```python
+star.PrintAvailableTracks()
+print(star.Units['Lx'])
+```
+
+---
+
+### Step 3: Plot an evolutionary track
+
+Tracks are stored in `star.Tracks` and also exposed as `Track` arrays.
+
+```python
+import matplotlib.pyplot as plt
+
+plt.plot(star.Tracks['Age'], star.Tracks['Lx'])
+# or: plt.plot(star.AgeTrack, star.LxTrack)
+plt.xlabel(f"Age [{star.Units['Age']}]")
+plt.ylabel(f"Lx [{star.Units['Lx']}]")
+plt.show()
+```
+
+---
+
+### Step 4: Find a value at a given age
+
+```python
+print(star.Value(150.0, 'Lx'))
+# or:
+print(star.Value(Age=150.0, Quantity='Lx'))
+# or (preferred when available):
+print(star.Lx(150.0))
+```
+
+---
+
+### Step 5: Save and reload
+
+```python
+star.Save(filename="star.pickle")
+star2 = mors.Load("star.pickle")
+```
+
+---
+
+### Common pitfalls
+- `Age` is in **Myr** (not years).
+- `Omega` is in **Ω☉** (not rad/s).
+- Fitting with `Age=...` + `Omega=...` may fail at late ages or for extreme rotation rates.
diff --git a/docs/How-to/habitablezone.md b/docs/How-to/habitablezone.md
new file mode 100644
index 0000000..8fb915b
--- /dev/null
+++ b/docs/How-to/habitablezone.md
@@ -0,0 +1,79 @@
+## 8. Find habitable zone boundaries (How-to)
+
+### Goal
+Compute habitable zone (HZ) boundary orbital distances (AU) as a function of stellar **mass** and **age** using `mors.aOrbHZ`.
+
+### Prerequisites
+Install MORS and download stellar evolution data:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Mstar` in **M☉**, `Age` in **Myr**, returned distances in **AU**.
+
+**Citation note:** If you use these HZ boundaries in research, cite **Kopparapu et al. (2013)** (HZ prescription) and **Spada et al. (2013)** (stellar parameters used).
+
+---
+
+### Step 1: Compute HZ boundaries for a single star mass
+
+By default, `aOrbHZ` uses stellar parameters at **5000 Myr** if you don’t specify `Age`.
+
+```python
+import mors
+
+hz = mors.aOrbHZ(Mstar=1.0)
+```
+
+The returned dictionary contains these keys:
+- `RecentVenus`
+- `RunawayGreenhouse`
+- `MoistGreenhouse`
+- `MaximumGreenhouse`
+- `EarlyMars`
+- `HZ`
+
+Print them:
+
+```python
+for k in ["RecentVenus", "RunawayGreenhouse", "MoistGreenhouse",
+ "MaximumGreenhouse", "EarlyMars", "HZ"]:
+ print(k, hz[k])
+```
+
+`HZ` is defined (Johnstone et al. 2020) as the average of the `RunawayGreenhouse` and `MoistGreenhouse` orbital distances.
+
+---
+
+### Step 2: Specify a stellar age
+
+```python
+import mors
+hz_1gyr = mors.aOrbHZ(Mstar=1.0, Age=1000.0) # Myr
+print(hz_1gyr["HZ"])
+```
+
+---
+
+### Step 3: Compute HZ boundaries for multiple masses
+
+If you pass an array of masses, you get arrays back in the dictionary:
+
+```python
+import numpy as np
+import mors
+
+masses = np.array([0.3, 0.5, 0.8, 1.0])
+hz = mors.aOrbHZ(Mstar=masses, Age=5000.0)
+
+print(hz["HZ"]) # array of AU
+print(hz["RunawayGreenhouse"])
+```
+
+---
+
+### Common pitfalls
+- `Age` is in **Myr** (not years).
+- When passing arrays for `Mstar`, each returned dict entry becomes an array of the same length.
diff --git a/docs/How-to/index.md b/docs/How-to/index.md
new file mode 100644
index 0000000..5295db0
--- /dev/null
+++ b/docs/How-to/index.md
@@ -0,0 +1,14 @@
+# How-to guides
+
+Welcome to the how-to guides for MORS. These are task-focused recipes for common MORS workflows.
+
+1. [(Developer) installation](installation.md)
+2. [Calculate a star's evolution](evolution.md)
+3. [Calculate a cluster's evolution](clusters.md)
+4. [Using model percentiles for initial rotation](distribution_percentile.md)
+5. [Setting custom simulation parameters](set_parameters.md)
+6. [Find stellar rotation and activity quantities](activity_quantities.md)
+7. [Find stellar quantities with stellar evolution tracks](track_quantities.md)
+8. [Find habitable zone boundaries](habitablezone.md)
+9. [Find stellar activity timelines](activity_timelines.md)
+
diff --git a/docs/How-to/installation.md b/docs/How-to/installation.md
new file mode 100644
index 0000000..df45d2e
--- /dev/null
+++ b/docs/How-to/installation.md
@@ -0,0 +1,75 @@
+# Installation
+
+### Prerequisites
+- **Python:** >3.11 installed
+- **pip:** available (`python -m pip --version`)
+- **Git:** only needed for the developer install (`git --version`)
+- **Internet access:** required once to download the stellar evolution tracks
+- *(Optional)* **Conda/Anaconda/Miniconda:** only if you want to use a conda environment
+
+### 0. Optional: Conda/virtual environment
+
+Create and activate a Conda environment (requires `conda` installed):
+```bash
+conda create -n mors python=3.11 -y
+conda activate mors
+```
+
+No `conda`? create and activate a virtual environment (venv):
+```bash
+python -m venv .venv
+source .venv/bin/activate
+```
+
+### 1. Basic install
+
+The Forming Worlds Mors package is available on PyPI. Run the following command to install
+
+```sh
+pip install fwl-mors
+```
+### 2. Developer install
+
+You can alternatively download the source code from GitHub somewhere on your computer using
+
+```
+git clone git@github.com:FormingWorlds/MORS.git
+```
+
+Then run the following command inside the main directory to install the code (check the pyproject.toml file for dependencies)
+
+```
+pip install -e .
+```
+
+
+### 3. Stellar evolution tracks
+
+The code requires also a set of stellar evolution data, stored in the [OSF repository](https://osf.io/9u3fb/).
+
+You can use `mors download all` to download the data. This will download and extract package stellar evolution tracks data.
+
+By default, MORS stores the data in based on the [XDG specification](https://specifications.freedesktop.org/basedir-spec/latest/).
+You can check the location by typing `mors env` in your terminal.
+You can override the path using the `FWL_DATA` environment variable, e.g.
+
+```console
+export FWL_DATA=...
+```
+
+Where ... should be replaced with the path to your main data directory. To make this permanent on Ubuntu, use
+
+```console
+gedit ~/.profile
+```
+
+and add the export command to the bottom of the file.
+
+Alternatively, when creating a star object in your Python script, you can specify the path to a directory where evolution tracks are stored using the starEvoDir keyword
+
+```python
+import mors
+myStar = mors.StarEvo(starEvoDir=...)
+```
+
+where ... can be given as the path relative to the current directory. When this is done, no environmental variable needs to be set.
\ No newline at end of file
diff --git a/docs/How-to/set_parameters.md b/docs/How-to/set_parameters.md
new file mode 100644
index 0000000..0d33107
--- /dev/null
+++ b/docs/How-to/set_parameters.md
@@ -0,0 +1,101 @@
+## Setting custom simulation parameters (How-to)
+
+### Goal
+Create simulation parameters for a `Star` (or `Cluster`) run, and optionally restrict the ages saved in output tracks using `AgesOut`. This overrides MORS default simulation parameters.
+
+### Prerequisites
+Make sure MORS and stellar evolution data are installed:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+---
+
+### Step 1: Start from the default parameter set
+
+Create a parameter dictionary identical to the internal defaults (`paramsDefault` in `parameters.py`):
+
+```python
+import mors
+params = mors.NewParams()
+```
+
+This gives you a normal Python dictionary you can edit.
+
+---
+
+### Step 2: Change one or more parameters
+
+**Option A: edit the dictionary after creation**
+```python
+import mors
+params = mors.NewParams()
+
+params["param1"] = 1.5
+params["param2"] = 2.5
+```
+
+**Option B: set values when calling `NewParams` (recommended)**
+```python
+import mors
+params = mors.NewParams(param1=1.5, param2=2.5)
+```
+
+---
+
+### Step 3: Pass the parameter dictionary into `Star` (or `Cluster`)
+
+```python
+import mors
+star = mors.Star(Mstar=1.0, Omega=10.0, params=params)
+```
+
+(Use the same `params=` keyword when creating a `mors.Cluster`.)
+
+---
+
+### Step 4: Discover available parameters
+
+To print a complete list of parameters (and their meanings/values as exposed by MORS):
+
+```python
+import mors
+mors.PrintParams()
+```
+
+You can also inspect the source file `parameters.py` in your MORS installation.
+
+---
+
+### Step 5: Restrict output ages with `AgesOut` (optional)
+
+By default, MORS computes tracks on its internal age grid. If you only need values at specific ages, set `AgesOut` so the saved tracks contain only those ages **plus the starting age (1 Myr)**.
+
+**Single age**
+```python
+import mors
+star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=100.0) # Myr
+```
+
+**Multiple ages**
+```python
+import numpy as np
+import mors
+
+ages = np.array([100.0, 200.0, 300.0, 400.0]) # Myr (ascending)
+star = mors.Star(Mstar=1.0, Omega=10.0, AgesOut=ages)
+```
+
+Notes:
+- The simulation ends at the **largest** age in `AgesOut`.
+- Provide `AgesOut` in **ascending** order.
+- If you later call interpolating helpers like `star.Lx(Age)` at arbitrary ages, results can be inaccurate if `AgesOut` is too sparse.
+- Very large `AgesOut` grids can slow down calculations significantly.
+
+---
+
+### Common pitfalls
+- `AgesOut` is in **Myr**.
+- If you intend to query many arbitrary ages later, prefer the default age grid (don’t over-thin `AgesOut`).
diff --git a/docs/How-to/track_quantities.md b/docs/How-to/track_quantities.md
new file mode 100644
index 0000000..6fdf79e
--- /dev/null
+++ b/docs/How-to/track_quantities.md
@@ -0,0 +1,161 @@
+## Find stellar quantities using evolution tracks (How-to)
+
+### Goal
+Find basic stellar evolution properties (radius, luminosity, convective turnover time, moments of inertia, etc.) as a function of **stellar mass** and **age**, using the Spada et al. (2013) tracks bundled with MORS. Optionally, stellar evolution quantities according to the Baraffe model (Baraffe et al., 2002) can be found.
+
+### Prerequisites
+Install MORS and download the required data:
+
+```bash
+pip install fwl-mors
+mors download all
+```
+
+> **Units:** `Mstar` in **M☉**, `Age` in **Myr**. Output units depend on the quantity (listed below).
+
+## Spada tracks
+
+### Available quantities
+
+- `Rstar` — radius (**R☉**)
+- `Lbol` — bolometric luminosity (**L☉**)
+- `Teff` — effective temperature (**K**)
+- `Itotal` — total moment of inertia (**g cm²**)
+- `Icore` — core moment of inertia (**g cm²**)
+- `Ienv` — envelope moment of inertia (**g cm²**)
+- `Mcore` — core mass (**M☉**)
+- `Menv` — envelope mass (**M☉**)
+- `Rcore` — core radius (**R☉**)
+- `tau` — convective turnover time (**days**)
+- `dItotaldt` — d(Itotal)/dt (**g cm² Myr⁻¹**)
+- `dIcoredt` — d(Icore)/dt (**g cm² Myr⁻¹**)
+- `dIenvdt` — d(Ienv)/dt (**g cm² Myr⁻¹**)
+- `dMcoredt` — d(Mcore)/dt (**M☉ Myr⁻¹**)
+- `dRcoredt` — d(Rcore)/dt (**R☉ Myr⁻¹**)
+
+**Important definition note:** “core” and “envelope” here follow the rotation model definitions:
+- core = everything interior to the outer convective zone
+- envelope = the outer convective zone
+---
+
+### Option 1: Call a property function directly
+
+Each quantity has a dedicated function. For example, stellar radius:
+
+```python
+import mors
+Rstar = mors.Rstar(1.0, 1000.0) # Rsun
+print(Rstar)
+```
+
+---
+
+### Option 2: Use the generic accessor (`Value`)
+
+If you want to choose the quantity by name at runtime:
+
+```python
+import mors
+Rstar = mors.Value(1.0, 1000.0, "Rstar")
+print(Rstar)
+```
+
+---
+
+### Option 3: Vectorized inputs (arrays/lists)
+
+These functions accept lists or NumPy arrays for mass and/or age.
+
+**Multiple masses, one age → 1D array**
+```python
+import numpy as np
+import mors
+
+masses = np.array([0.8, 1.0, 1.2])
+R = mors.Rstar(masses, 1000.0)
+print(R.shape) # (3,)
+```
+
+**One mass, multiple ages → 1D array**
+```python
+import numpy as np
+import mors
+
+ages = np.array([10.0, 100.0, 1000.0])
+L = mors.Lbol(1.0, ages)
+print(L.shape) # (3,)
+```
+
+**Multiple masses and multiple ages → 2D array**
+```python
+import numpy as np
+import mors
+
+masses = np.array([0.8, 1.0, 1.2])
+ages = np.array([10.0, 100.0, 1000.0])
+T = mors.Teff(masses, ages)
+print(T.shape) # (len(masses), len(ages))
+```
+
+**Multiple quantities via `Value`**
+```python
+import numpy as np
+import mors
+
+masses = np.array([0.9, 1.0])
+ages = np.array([100.0, 1000.0])
+vals = mors.Value(masses, ages, ["Rstar", "Lbol", "tau"])
+# adds an extra dimension for the quantity list
+print(vals.shape)
+```
+
+---
+
+### Performance tip: pre-load a mass track (`LoadTrack`)
+
+The first time you call one of these functions, MORS loads and caches Spada tracks and writes a cache file (e.g., `SEmodels.pickle`) to speed up future runs. This cache can be deleted safely; it will be regenerated as needed.
+
+If you will query many ages for a particular mass, pre-load that mass track:
+
+```python
+import mors
+mors.LoadTrack(1.02)
+```
+
+If it’s already loaded, this call does nothing.
+
+---
+
+## Baraffe tracks
+
+MORS also provides access to Baraffe et al. (2002) tracks, which use **different units** than the Spada helpers above.
+
+> **Baraffe units:** `Mstar` in **M☉** (valid range: ~0.01–1.4), `time` in **years (yr)**.
+
+### Step 1. Load a Baraffe track (with interpolation)
+```python
+import mors
+
+Mstar = 0.5
+baraffe = mors.BaraffeTrack(Mstar)
+```
+
+`BaraffeTrack` performs mass interpolation (if needed) and interpolates time onto a fine grid.
+
+### Step 2. Query radius, luminosity, and “solar constant”
+```python
+time_yr = 1e7 # years
+
+Rstar = baraffe.BaraffeStellarRadius(time_yr) # Rsun
+Lbol = baraffe.BaraffeLuminosity(time_yr) # Lsun
+Flux = baraffe.BaraffeSolarConstant(time_yr, 1.0) # W m^-2 at 1 AU
+```
+
+`BaraffeSolarConstant(time, distance)` expects `distance` in **AU**.
+
+---
+
+### Common pitfalls
+- Spada helpers use `Age` in **Myr**; Baraffe helpers use `time` in **yr**.
+- “core” and “envelope” are defined by the rotation model (not the hydrogen-burning core definition).
+- The first Spada call can be slower due to model loading + cache generation.
diff --git a/docs/Tutorials/first_run.md b/docs/Tutorials/first_run.md
new file mode 100644
index 0000000..d315299
--- /dev/null
+++ b/docs/Tutorials/first_run.md
@@ -0,0 +1,151 @@
+# Tutorial: First run
+
+## What you’ll do
+By the end of this tutorial you will:
+1. Install MORS and download the required data
+2. Create a `Star` and compute evolutionary tracks
+3. Inspect available tracks + units
+4. Plot one track
+5. Find stellar values at specific ages
+6. Save and reload the result
+
+> **Time/units cheat sheet:** `Mstar` in **M☉**, `Age` in **Myr**, `Prot` in **days**, `Omega` in **Ω☉**.
+
+---
+
+## 0) Prerequisites
+You need:
+
+- Python 3 environment with `pip`
+- A working internet connection (for downloading data once)
+
+**Optional**: Create and activate a Conda environment (requires `conda` installed):
+```bash
+conda create -n mors python=3.11 -y
+conda activate mors
+```
+
+No `conda`? create and activate a virtual environment (venv):
+```bash
+python -m venv .venv
+source .venv/bin/activate
+```
+
+---
+
+## 1) Install MORS
+```bash
+pip install fwl-mors
+```
+
+Quick sanity check:
+```bash
+python -c "import mors; print('mors imported:', mors.__version__ if hasattr(mors,'__version__') else 'ok')"
+```
+
+---
+
+## 2) Download stellar evolution data
+MORS requires stellar evolution tracks (downloaded once):
+```bash
+mors download all
+```
+
+See where data are stored:
+```bash
+mors env
+```
+
+If you want to place data somewhere else, set `FWL_DATA` (optional):
+```bash
+export FWL_DATA=/path/to/your/data
+```
+
+---
+
+## 3) Create your first star
+Create a 1 M☉ star with an initial rotation period of 2.7 days (at ~1 Myr):
+```python
+import mors
+star = mors.Star(Mstar=1.0, Prot=2.7)
+```
+
+Alternative: specify initial rotation as Ω/Ω☉:
+```python
+star = mors.Star(Mstar=1.0, Omega=10.0)
+```
+
+---
+
+## 4) Inspect what tracks exist
+Print track names and units:
+```python
+star.PrintAvailableTracks()
+print("Lx units:", star.Units.get("Lx"))
+```
+
+You can access arrays either via the `Tracks` dict:
+```python
+age = star.Tracks["Age"]
+lx = star.Tracks["Lx"]
+```
+or via convenience attributes:
+```python
+age = star.AgeTrack
+lx = star.LxTrack
+```
+
+---
+
+## 5) Plot a track
+```python
+import matplotlib.pyplot as plt
+
+plt.plot(star.AgeTrack, star.LxTrack)
+plt.xlabel(f"Age [{star.Units['Age']}]")
+plt.ylabel(f"Lx [{star.Units['Lx']}]")
+plt.show()
+```
+
+If you see a plot and no errors, your installation + data are working.
+
+---
+
+## 6) Find stellar values at specific ages
+Use the generic accessor:
+```python
+print(star.Value(Age=150.0, Quantity="Lx"))
+```
+
+Or a convenience method (when available):
+```python
+print(star.Lx(150.0))
+```
+
+---
+
+## 7) (Optional) Try percentiles: slow/medium/fast rotators
+This uses the built-in model distribution:
+```python
+slow = mors.Star(Mstar=1.0, percentile="slow") # 5th percentile
+medium = mors.Star(Mstar=1.0, percentile="medium") # 50th percentile
+fast = mors.Star(Mstar=1.0, percentile="fast") # 95th percentile
+
+print("slow percentile:", slow.percentile)
+print("fast percentile:", fast.percentile)
+```
+
+---
+
+## 8) Save and reload (recommended)
+Save:
+```python
+star.Save(filename="star.pickle")
+```
+
+Reload later:
+```python
+import mors
+star2 = mors.Load("star.pickle")
+print(star2.Lx(150.0))
+```
diff --git a/docs/Tutorials/index.md b/docs/Tutorials/index.md
new file mode 100644
index 0000000..6b8f6bd
--- /dev/null
+++ b/docs/Tutorials/index.md
@@ -0,0 +1,5 @@
+# Tutorials
+
+Welcome to the tutorials for MORS. So far, there is only one available – but more to come!
+
+1. [First run](first_run.md)
\ No newline at end of file
diff --git a/docs/index.md b/docs/index.md
index 7d211af..a85034d 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -1 +1,23 @@
-{!../README.md!}
+
+
+# MODEL FOR ROTATION OF STARS (MORS)
+
+**This code is distributed as a python package for the purpose of the [PROTEUS framework](https://proteus-framework.org/PROTEUS/), a coupled simulation tool for the long-term evolution of atmospheres and interiors of rocky planets.
+The MORS package solves specifically the stellar rotation and evolution. It is based on the [original code](https://www.aanda.org/articles/aa/pdf/2021/05/aa38407-20.pdf) and model developed by Colin P. Johnstone.**
+
+**Original Author:** Colin P. Johnstone
+
+This code solves the stellar rotation and XUV evolution model presented in Johnstone et al. (2021). The package can be used to calculate evolutionary tracks for stellar rotaton and X-ray, EUV, and Ly-alpha emission for stars with masses between 0.1 and 1.25 Msun and has additional functionality such as allowing the user to get basic stellar parameters such as stellar radius and luminosity as functions of mass and age using the stellar evolution models of Spada et al. (2013). When publishing results that were calculated using this code, both the Johnstone et al. (2020) paper and Spada et al. (2013) should be cited.
+
+**NOTE:** This version contains the fix for the error in the equation converting EUV1 to EUV2.
+
+### Contributors
+
+| Name | Email address |
+| - | - |
+| Colin P. Johnstone | colinjohnstone@gmail.com |
+| Laurent Soucasse | l.soucasse@esciencecenter.nl |
+| Harrison Nicholls | harrison.nicholls@physics.ox.ac.uk |
+| Stef Smeets | s.smeets@esciencecenter.nl |
+| Tim Lichtenberg | tim.lichtenberg@rug.nl |
+| Karen Stuitje | e.k.e.stuitje@student.rug.nl |
\ No newline at end of file
diff --git a/mkdocs.yml b/mkdocs.yml
index 9490b0e..cb8c375 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -1,26 +1,45 @@
site_name: MORS
-site_url: https://fwl-mors.readthedocs.io
+site_url: https://proteus-framework.org/MORS/
repo_url: https://github.com/FormingWorlds/MORS
repo_name: GitHub
nav:
- Home: index.md
+
+ - How-to guides:
+ - Overview: How-to/index.md
+ - Installation: How-to/installation.md
+ - Calculate stellar evolution: How-to/evolution.md
+ - Calculate cluster evolution: How-to/clusters.md
+ - Using model percentiles for initial rotation: How-to/distribution_percentile.md
+ - Setting custom simulation parameters: How-to/set_parameters.md
+ - Find stellar rotation and activity quantities: How-to/activity_quantities.md
+ - Find stellar quantities with stellar evolution tracks: How-to/track_quantities.md
+ - Find habitable zone boundaries: How-to/habitablezone.md
+ - Find stellar activity timelines: How-to/activity_timelines.md
+
+ - Tutorials:
+ - Overview: Tutorials/index.md
+ - First run: Tutorials/first_run.md
+
- Contributing: CONTRIBUTING.md
- Code of Conduct: CODE_OF_CONDUCT.md
- Source code: https://github.com/FormingWorlds/MORS
- Issues page: https://github.com/FormingWorlds/MORS/issues
- - 🔗 PROTEUS: https://fwl-proteus.readthedocs.io/
- - 🔗 JANUS: https://fwl-proteus.readthedocs.io/projects/janus/en/latest/
- - 🔗 ZEPHYRUS: https://github.com/FormingWorlds/ZEPHYRUS/
- - 🔗 CALLIOPE: https://fwl-proteus.readthedocs.io/projects/calliope/en/latest/
+ - 🔗 PROTEUS: https://proteus-framework.org/PROTEUS/
+ - 🔗 JANUS: https://proteus-framework.org/JANUS/
+ - 🔗 ZEPHYRUS: https://proteus-framework.org/ZEPHYRUS/
+ - 🔗 CALLIOPE: https://proteus-framework.org/CALLIOPE/
- 🔗 AGNI: https://h-nicholls.space/AGNI/dev/
- - 🔗 LovePy: https://github.com/nichollsh/lovepy
- - 🔗 VULCAN: https://github.com/FormingWorlds/VULCAN
- - 🔗 Zalmoxis: https://zalmoxis.readthedocs.io/en/latest/
- - 🔗 Aragog: https://github.com/FormingWorlds/aragog
+ - 🔗 Obliqua: https://proteus-framework.org/Obliqua/
+ - 🔗 VULCAN: https://proteus-framework.org/VULCAN/
+ - 🔗 Zalmoxis: https://proteus-framework.org/ZALMOXIS/
+ - 🔗 Aragog: https://proteus-framework.org/aragog/
theme:
name: material
+ favicon: Images/PROTEUS_black_on_white_logo_only.png
+ logo: Images/PROTEUS_white_on_black.png
palette:
primary: black
accent: deep orange