adding a more complex example to the docs

The [basic usage](basic_usage.md) page shows how to run a simulation and get outputs, here we want to show how a script with some more advanced reading in and writing out of data would look like with some customized plotting.

If you wish to make your own plots we recommend the user to learn how to use [matplotlib](https://matplotlib.org/) and/or [xarrays](https://docs.xarray.dev/en/stable/).

The [basic usage](basic_usage.md) page demonstrates how to run a simulation and retrieve results. This page provides an example of a more advanced workflow using **pyThermoGIS**, including some plotting examples to illustrate the outputs.

**Note: pythermogis has been made to enable users to run doublet simulations; it is not our goal to make lots of customised plotting functions. 

The reason is that when it comes to plotting users have a specific idea of how they want their plots to look and it is hard to cater to everyone. It is better for us to give you the tools to run the simulations and let you design your plots 

yourselves.**

!!! info "Plotting"

    pyThermoGIS is designed to enable users to run geothermal doublet simulations. Customized plotting is intentionally limited, as users often have specific preferences for visual presentation. Rather than attempt to accommodate all 

    potential styles, our focus is on providing reliable simulation tools. For creating your own visualizations, we recommend using libraries such as [matplotlib](https://matplotlib.org/) and/or [xarray](https://docs.xarray.dev/en/stable/).

Here is an example which:

1. reads in grids of: mean thickness, thickness standard deviation, net-to-gross, porosity, mean permeability, ln(permeability) standard deviation

2. runs a doublet simulations and calculates economics across all the non-nan cells of the grids for p-values of 10%-90% (with increments of 10%)

3. writes the results out to file

4. plots maps of power, capex and npv for the pvalues: 10%, 50%, 90%

5. plots the npv curve for a single location on the grids

### This example demonstrates the following steps:

1. **Reading input grids** for:

     - Mean thickness  

     - Thickness standard deviation  

     - Net-to-gross ratio  

     - Porosity  

     - Mean permeability  

     - Log(permeability) standard deviation  

2. **Running doublet simulations** and calculating economic indicators (e.g., power production, CAPEX, NPV) across all non-NaN cells in the input grids for a range of p-values (10% to 90% in 10% increments).

3. **Exporting results** to file.

4. **Plotting result maps** for selected p-values (10%, 50%, and 90%) for:

     - Power output  

     - Capital expenditure (CAPEX)  

     - Net Present Value (NPV)

5. **Plotting an NPV curve** at a single specified location on the grid.

Example input data for this workflow is available in the `/resources/example_data` directory of the repository.

The example input data can be found in the `/resoureces/example_data` directory in the repository.

```python

from pythermogis import calculate_doublet_performance

from pygridsio import read_grid, plot_grid

from matplotlib import pyplot as plt

import xarray as xr

from pygridsio import read_grid, plot_grid # for reading and plotting grids

from matplotlib import pyplot as plt # for more customized plotting

import xarray as xr # to handle the in and output data

import numpy as np

from pathlib import Path

from os import path

from pathlib import Path # to handle pathing

# the location of the input data: the data can be found in the resources/example_data directory of the repo

input_data_path = Path(path.dirname(__file__), "resources") / "test_input" / "example_data"

input_data_path = Path("path/to/example_data")

# create a directory to write the output files to

output_data_path = Path(path.dirname(__file__), "resources") / "test_output" / "example_data"

output_data_path = Path("path/to/output_data")

output_data_path.mkdir(parents=True, exist_ok=True)

# if set to True then simulation is always run, otherwise pre-calculated results are read (if available)

    results_p_value = results.sel(p_value=p_value)

    for i, variable in enumerate(variables_to_plot):

        plot_grid(results_p_value[variable], axes=axes[i,j], add_netherlands_shapefile=True)  # See documentation on plot_grid in pygridsio, you can also provide your own shapefile

plt.tight_layout()  # ensure there is enough spacing

plt.tight_layout()  # ensure there is enough spacing between subplots

plt.savefig(output_data_path / "maps.png")

# plot net-present value at a single location as a function of p-value and find the probability of success

# plot map

plot_grid(results.sel(p_value=50).npv,axes=axes[1], add_netherlands_shapefile=True)

axes[1].scatter(x,y,marker="x",color="tab:red")

plt.tight_layout()  # ensure there is enough spacing

plt.tight_layout()  # ensure there is enough spacing between subplots

plt.savefig(output_data_path / "npv.png")

```

And the output plots are:

And the output plots produced by the code are:

![maps.png](../images/maps.png)

![npv.png](../images/npv.png)