diff --git a/README.md b/README.md
index 11a233f..24b6c57 100644
--- a/README.md
+++ b/README.md
@@ -2,6 +2,7 @@
[Installation Guide](https://pyautolens.readthedocs.io/en/latest/installation/overview.html) |
[PyAutoLens readthedocs](https://pyautolens.readthedocs.io/en/latest/index.html) |
+[Browse Chapter 1 With Images](markdown/README.md) |
[autolens_workspace](https://github.com/PyAutoLabs/autolens_workspace)
diff --git a/config/build/markdown_examples.yaml b/config/build/markdown_examples.yaml
new file mode 100644
index 0000000..c34a517
--- /dev/null
+++ b/config/build/markdown_examples.yaml
@@ -0,0 +1,21 @@
+# Curated examples rendered to executed markdown pages (markdown/) with real
+# output images, for GitHub browsing. Built by PyAutoBuild's generate_markdown.py.
+# Batch 2b: chapter_1_introduction (no non-linear searches — fast).
+- script: scripts/chapter_1_introduction/tutorial_0_visualization.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_1_grids_and_galaxies.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_2_ray_tracing.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_3_more_ray_tracing.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_4_point_sources.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_5_lensing_formalism.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_6_data.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_7_fitting.py
+ max_minutes: 30
+- script: scripts/chapter_1_introduction/tutorial_8_summary.py
+ max_minutes: 30
diff --git a/markdown/README.md b/markdown/README.md
new file mode 100644
index 0000000..2580590
--- /dev/null
+++ b/markdown/README.md
@@ -0,0 +1,15 @@
+# HowToLens examples, executed — browse with output images
+
+Every page below is the corresponding example script **fully executed**, rendered to markdown with its real output images, so you can read the examples on GitHub exactly as they run. Each page links back to the `.py` script and Jupyter notebook it was generated from.
+
+- [Tutorial 0: Visualization](chapter_1_introduction/tutorial_0_visualization.md) — from `scripts/chapter_1_introduction/tutorial_0_visualization.py`
+- [HowToLens: Introduction](chapter_1_introduction/tutorial_1_grids_and_galaxies.md) — from `scripts/chapter_1_introduction/tutorial_1_grids_and_galaxies.py`
+- [Tutorial 2: Ray Tracing](chapter_1_introduction/tutorial_2_ray_tracing.md) — from `scripts/chapter_1_introduction/tutorial_2_ray_tracing.py`
+- [Tutorial 5: More Ray Tracing](chapter_1_introduction/tutorial_3_more_ray_tracing.md) — from `scripts/chapter_1_introduction/tutorial_3_more_ray_tracing.py`
+- [Tutorial 4: Point Sources](chapter_1_introduction/tutorial_4_point_sources.md) — from `scripts/chapter_1_introduction/tutorial_4_point_sources.py`
+- [Tutorial 5: Lensing Formalism](chapter_1_introduction/tutorial_5_lensing_formalism.md) — from `scripts/chapter_1_introduction/tutorial_5_lensing_formalism.py`
+- [Tutorial 6: Data](chapter_1_introduction/tutorial_6_data.md) — from `scripts/chapter_1_introduction/tutorial_6_data.py`
+- [Tutorial 7: Fitting](chapter_1_introduction/tutorial_7_fitting.md) — from `scripts/chapter_1_introduction/tutorial_7_fitting.py`
+- [Tutorial 9: Summary](chapter_1_introduction/tutorial_8_summary.md) — from `scripts/chapter_1_introduction/tutorial_8_summary.py`
+
+These pages are regenerated manually by PyAutoBuild's `generate_markdown.py` when a curated script changes.
diff --git a/markdown/chapter_1_introduction/tutorial_0_visualization.md b/markdown/chapter_1_introduction/tutorial_0_visualization.md
new file mode 100644
index 0000000..3b46d13
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_0_visualization.md
@@ -0,0 +1,193 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_0_visualization.py`](../../scripts/chapter_1_introduction/tutorial_0_visualization.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_0_visualization.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_0_visualization.ipynb).
+
+Tutorial 0: Visualization
+=========================
+
+In this tutorial, we quickly cover visualization in **PyAutoLens** and make sure images display
+clearly in your Jupyter notebook and on your computer screen.
+
+__Contents__
+
+- **Directories:** **PyAutoLens assumes** the working directory is `autolens_workspace` on your hard-disk.
+- **Dataset:** Load and plot the strong lens dataset.
+- **Subplots:** In addition to plotting individual figures, **PyAutoLens** can plot `subplots` which show multiple.
+- **Plot Customization:** Does the figure display correctly on your computer screen?
+- **Overlays:** Overlays such as critical curves and image positions are added using the `lines=` and `positions=`.
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+If the printed working directory does not match the workspace path on your computer, you can manually set it
+as follows (the example below shows the path I would use on my laptop. The code is commented out so you do not
+use this path in this tutorial!
+
+
+```python
+# workspace_path = "/Users/Jammy/Code/PyAuto/autolens_workspace"
+# #%cd $workspace_path
+# print(f"Working Directory has been set to `{workspace_path}`")
+```
+
+__Dataset__
+
+The `dataset_path` specifies where the dataset is located, which is the
+directory `autolens_workspace/dataset/imaging/simple__no_lens_light`.
+
+There are many example simulated images of strong lenses in this directory that will be used throughout the
+**HowToLens** lectures.
+
+
+```python
+from pathlib import Path
+
+import autolens as al
+import autolens.plot as aplt
+
+dataset_path = Path("dataset") / "imaging" / "simple__no_lens_light"
+```
+
+We now load this dataset from .fits files and create an instance of an `Imaging` object.
+
+
+```python
+dataset = al.Imaging.from_fits(
+ data_path=dataset_path / "data.fits",
+ noise_map_path=dataset_path / "noise_map.fits",
+ psf_path=dataset_path / "psf.fits",
+ pixel_scales=0.1,
+)
+```
+
+We can plot an image with `aplt.plot_array()`, passing the data array and a title.
+
+
+```python
+aplt.plot_array(array=dataset.data, title="Dataset Image")
+```
+
+
+
+
+
+
+
+__Subplots__
+
+In addition to plotting individual figures, **PyAutoLens** can plot `subplots` which show multiple
+views of the dataset at once.
+
+The `aplt.subplot_imaging_dataset()` function plots the data, noise-map and PSF together.
+
+
+```python
+aplt.subplot_imaging_dataset(dataset=dataset)
+```
+
+
+
+
+
+
+
+__Plot Customization__
+
+Does the figure display correctly on your computer screen?
+
+If not, the default matplotlib settings can be customized via the config files in:
+
+ autolens_workspace/config/visualize/
+
+Key config entries:
+
+ - `mat_wrap.yaml` -> Figure -> figure: -> figsize
+ - `mat_wrap.yaml` -> YLabel -> figure: -> fontsize
+ - `mat_wrap.yaml` -> XLabel -> figure: -> fontsize
+ - `mat_wrap.yaml` -> TickParams -> figure: -> labelsize
+ - `mat_wrap.yaml` -> YTicks -> figure: -> labelsize
+ - `mat_wrap.yaml` -> XTicks -> figure: -> labelsize
+
+For quick one-off adjustments you can pass `title=`, `colormap=`, and `use_log10=` directly:
+
+
+```python
+aplt.plot_array(array=dataset.data, title="Dataset Image (Log10)", use_log10=True)
+```
+
+
+
+
+
+
+
+__Overlays__
+
+Overlays such as critical curves and image positions are added using the `lines=` and `positions=`
+keyword arguments.
+
+For example, we can compute the critical curves of a tracer and overlay them on the image.
+
+
+```python
+grid = al.Grid2D.uniform(shape_native=(100, 100), pixel_scales=0.05)
+
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.0, 0.0)),
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.SersicCoreSph(
+ centre=(0.0, 0.0), intensity=1.0, effective_radius=0.5, sersic_index=2.0
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+
+tangential_critical_curve_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).tangential_critical_curve_list_from(grid=grid)
+
+aplt.plot_array(
+ array=tracer.image_2d_from(grid=grid),
+ title="Tracer Image with Critical Curves",
+ lines=tangential_critical_curve_list,
+)
+```
+
+
+
+
+
+
+
+__Wrap Up__
+
+Throughout the lectures you'll see lots more visuals plotted on figures and subplots.
+
+The key plotting functions you'll use are:
+
+ - `aplt.plot_array(array, title, ...)` — plot any 2D array.
+ - `aplt.plot_grid(grid, title, ...)` — plot a 2D grid of coordinates.
+ - `aplt.subplot_imaging_dataset(dataset)` — multi-panel dataset overview.
+ - `aplt.subplot_tracer(tracer, grid)` — multi-panel tracer overview.
+ - `aplt.subplot_fit_imaging(fit)` — multi-panel fit overview.
+
+Great! Hopefully, visualization in **PyAutoLens** is displaying nicely for us to get on with the
+**HowToLens** lecture series.
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_11_0.png b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_11_0.png
new file mode 100644
index 0000000..158af60
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_11_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_13_0.png b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_13_0.png
new file mode 100644
index 0000000..7ab7029
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_13_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_15_0.png b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_15_0.png
new file mode 100644
index 0000000..6d77181
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_15_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_9_0.png b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_9_0.png
new file mode 100644
index 0000000..f6e4fd0
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_0_visualization_files/tutorial_0_visualization_9_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies.md b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies.md
new file mode 100644
index 0000000..9b3a8ba
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies.md
@@ -0,0 +1,575 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_1_grids_and_galaxies.py`](../../scripts/chapter_1_introduction/tutorial_1_grids_and_galaxies.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_1_grids_and_galaxies.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_1_grids_and_galaxies.ipynb).
+
+HowToLens: Introduction
+=======================
+
+A strong gravitational lens is a system where two (or more) galaxies align perfectly down our line of sight from Earth
+such that the foreground galaxy's mass curves space-time in on itself, such that the light of a background source galaxy
+is deflected and magnified. This means we can see the background source galaxy multiple times, as multiple arcs or rings,
+because multiple paths through the foreground galaxy's mass are taken by the source's light.
+
+Here is a schematic of a strong gravitational lens:
+
+
+**Credit: F. Courbin, S. G. Djorgovski, G. Meylan, et al., Caltech / EPFL / WMKO**
+https://www.astro.caltech.edu/~george/qsolens/
+
+Today, Astronomers use computer software, statistical algorithms and image processing techniques analyse the light of
+these strong gravitational lenses. They are used for a wide range of scientific studies, including studying dark matter
+in the foreground lens galaxies, the morphology and structure of the background source galaxies, and even the expansion
+of the Universe itself.
+
+The **HowToLens** series of tutorials will teach you how to perform this analysis yourself, using the open-source
+software package **PyAutoLens**. By the end of the **HowToLens** series, you'll be able to take an image of
+strong lens and study it using the same techniques that professional astronomers use today.
+
+Tutorial 1: Grids And Galaxies
+==============================
+
+In this tutorial, we will introduce the first fundamental concepts and quantities used to study strong lenses.
+These concepts will enable us to create images of galaxies and analyze how their light is distributed across space.
+Additionally, we will explore how adjusting various properties of galaxies can alter their appearance. For instance,
+we can change the size of a galaxy, rotate it, or modify its brightness.
+
+To create these images, we first need to define 2D grids of \((y, x)\) coordinates. We will shift and rotate these
+grids to manipulate the appearance of the galaxy in the generated images. The grid will serve as the input for light
+profiles, which are analytic functions that describe the distribution of a galaxy's light. By evaluating these light
+profiles on the grid, we can effectively generate images that represent the structure and characteristics of galaxies.
+
+This tutorial won't yet perform any lensing calculations, which are introduced in the next tutorial.
+
+Here is an overview of what we'll cover in this tutorial:
+
+- **Grids**: We'll create a uniform grid of $(y,x)$ coordinates and show how it can be used to measure the light of a galaxy.
+- **Geometry**: We'll show how to shift and rotate a grid, and convert it to elliptical coordinates.
+- **Light Profiles**: We'll introduce light profiles, analytic functions that describe how a galaxy's light is distributed.
+- **Galaxies**: We'll create galaxies containing light profiles and show how to compute their image.
+- **Units**: We'll show how to convert the units of a galaxy's image to physical units like kiloparsecs.
+
+
+The imports below are required to run the HowToLens tutorials in a Jupiter notebook. They also import the
+`autolens` package and the `autolens.plot` module which are used throughout the tutorials.
+
+__Contents__
+
+- **Grids:** A `Grid2D` is a set of two-dimensional $(y,x)$ coordinates that represent points in space where we.
+- **Geometry:** The above grid is centered on the origin (0.0", 0.0").
+- **Light Profiles:** Galaxies are collections of stars, gas, dust, and other astronomical objects that emit light.
+- **One Dimension Projection:** We often want to calculative 1D quantities of a light profile, for example to plot how its light.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import autoarray as aa
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Grids__
+
+A `Grid2D` is a set of two-dimensional $(y,x)$ coordinates that represent points in space where we evaluate the
+light emitted by a galaxy.
+
+Each coordinate on the grid is referred to as a 'pixel'. This is because we use the grid to measure the brightness of a
+galaxy at each of these coordinates, allowing us to create an image of the galaxy.
+
+Grids are defined in units of 'arc-seconds' ("). An arc-second is a unit of angular measurement used by astronomers to
+describe the apparent size of objects in the sky.
+
+The `pixel_scales` parameter sets how many arc-seconds each pixel represents. For example, if `pixel_scales=0.1`,
+then each pixel covers 0.1" of the sky.
+
+We create a uniform 2D grid of 101 x 101 pixels with a pixel scale of 0.1", corresponding to an area
+of 10.1" x 10.1", spanning from -5.05" to 5.05" in both the y and x directions.
+
+
+```python
+grid = al.Grid2D.uniform(
+ shape_native=(
+ 101,
+ 101,
+ ), # The dimensions of the grid, which here is 100 x 100 pixels.
+ pixel_scales=0.1, # The conversion factor between pixel units and arc-seconds.
+)
+```
+
+We can visualize this grid as a uniform grid of dots, each representing a coordinate where the light is measured.
+
+
+```python
+aplt.plot_grid(grid=grid, title="Uniform Grid of Coordinates")
+```
+
+
+
+
+
+
+
+Each coordinate in the grid corresponds to an arc-second position. Below, we print a few of these coordinates to see
+the values.
+
+
+```python
+print("(y,x) pixel 0:")
+print(grid.native[0, 0]) # The coordinate of the first pixel.
+print("(y,x) pixel 1:")
+print(grid.native[0, 1]) # The coordinate of the second pixel.
+print("(y,x) pixel 2:")
+print(grid.native[0, 2]) # The coordinate of the third pixel.
+print("(y,x) pixel 100:")
+print(grid.native[1, 0]) # The coordinate of the 100th pixel.
+print("...")
+```
+
+ (y,x) pixel 0:
+ [ 5. -5.]
+ (y,x) pixel 1:
+ [ 5. -4.9]
+ (y,x) pixel 2:
+ [ 5. -4.8]
+ (y,x) pixel 100:
+ [ 4.9 -5. ]
+ ...
+
+
+Grids have two internal representations, `native` and `slim`:
+
+- `native`: A 2D array with shape [total_y_pixels, total_x_pixels, 2], where the 2 corresponds to the (y,x) coordinates.
+- `slim`: A 1D array with shape [total_y_pixels * total_x_pixels, 2], where the coordinates are 'flattened' into a single list.
+
+These formats are useful for different calculations and plotting. Here, we show the same coordinate using both formats.
+
+
+```python
+print("(y,x) pixel 0 (accessed via native):")
+print(grid.native[0, 0])
+print("(y,x) pixel 0 (accessed via slim 1D):")
+print(grid.slim[0])
+```
+
+ (y,x) pixel 0 (accessed via native):
+ [ 5. -5.]
+ (y,x) pixel 0 (accessed via slim 1D):
+ [ 5. -5.]
+
+
+We can also check the shapes of the `Grid2D` object in both `native` and `slim` formats. For this grid,
+the `native` shape is (101, 101, 2) and the `slim` shape is (10201, 2).
+
+
+```python
+print(grid.native.shape)
+print(grid.slim.shape)
+```
+
+ (101, 101, 2)
+ (10201, 2)
+
+
+For the HowToLens tutorials, you don't need to fully understand why grids have both native and slim representations.
+Just note that both are used for calculations and plotting.
+
+*Exercise*: Try creating grids with different shapes and pixel scales using the `al.Grid2D.uniform()` function above.
+Observe how the grid coordinates change when you adjust `shape_native` and `pixel_scales`.
+
+__Geometry__
+
+The above grid is centered on the origin (0.0", 0.0"). Sometimes, we need to shift the grid to be centered on a
+specific point, like the center of a galaxy.
+
+We can shift the grid to a new center, (y_c, x_c), by subtracting this center from each coordinate.
+
+
+```python
+centre = (0.3, 0.5) # Shifting the grid to be centered at y=1.0", x=2.0".
+
+grid_shifted = grid
+grid_shifted[:, 0] = grid_shifted[:, 0] - centre[0] # Shift in y-direction.
+grid_shifted[:, 1] = grid_shifted[:, 1] - centre[1] # Shift in x-direction.
+
+print("(y,x) pixel 0 After Shift:")
+print(grid_shifted.native[0, 0]) # The coordinate of the first pixel after shifting.
+```
+
+ (y,x) pixel 0 After Shift:
+ [ 4.7 -5.5]
+
+
+The grid is now centered around (0.3", 0.5"). We can plot the shifted grid to see this change.
+
+*Exercise*: Try shifting the grid to a different center, for example (0.0", 0.0") or (2.0", 3.0"). Observe how the
+center of the grid changes when you adjust the `centre` variable.
+
+
+```python
+aplt.plot_grid(grid=grid_shifted, title="Grid Centered Around (0.3, 0.5)")
+```
+
+
+
+
+
+
+
+Next, we can rotate the grid by an angle `phi` (in degrees). The rotation is counter-clockwise from the positive x-axis.
+
+To rotate the grid:
+
+1. Calculate the distance `radius` of each coordinate from the origin using $r = \sqrt{y^2 + x^2}$.
+2. Determine the angle `theta` counter clockwise from the positive x-axis using $\theta = \arctan(y / x)$.
+3. Adjust `theta` by the rotation angle and convert back to Cartesian coordinates via $y = r \sin(\theta)$ and $x = r \cos(\theta)$.
+
+
+```python
+angle_degrees = 60.0
+
+y = grid_shifted[:, 0]
+x = grid_shifted[:, 1]
+
+radius = np.sqrt(y**2 + x**2)
+theta = np.arctan2(y, x) - np.radians(angle_degrees)
+
+grid_rotated = grid_shifted
+grid_rotated[:, 0] = radius * np.sin(theta)
+grid_rotated[:, 1] = radius * np.cos(theta)
+
+print("(y,x) pixel 0 After Rotation:")
+print(grid_rotated.native[0, 0]) # The coordinate of the first pixel after rotation.
+```
+
+ (y,x) pixel 0 After Rotation:
+ [7.11313972 1.3203194 ]
+
+
+The grid has now been rotated 60 degrees counter-clockwise. We can plot it to see the change.
+
+*Exercise*: Try rotating the grid by a different angle, for example 30 degrees or 90 degrees. Observe how the grid
+changes when you adjust the `angle_degrees` variable.
+
+
+```python
+aplt.plot_grid(grid=grid_rotated, title="Grid Rotated 60 Degrees")
+```
+
+
+
+
+
+
+
+Next, we convert the rotated grid to elliptical coordinates using:
+
+$\eta = \sqrt{(x_r)^2 + (y_r)^2/q^2}$
+
+Where `q` is the axis-ratio of the ellipse and `(x_r, y_r)` are the rotated coordinates.
+
+Elliptical coordinates are a system used to describe positions in relation to an ellipse rather than a circle. They
+are particularly useful in astronomy when dealing with objects like galaxies, which often have elliptical shapes
+due to their inclination or intrinsic shape.
+
+*Exercise*: Try converting the grid to elliptical coordinates using a different axis-ratio, for example 0.3 or 0.8.
+What happens to the grid when you adjust the `axis_ratio` variable?
+
+
+```python
+axis_ratio = 0.5
+eta = np.sqrt((grid_rotated[:, 0]) ** 2 + (grid_rotated[:, 1]) ** 2 / axis_ratio**2)
+```
+
+Above, the angle $\phi$ (in degrees) was used to rotate the grid, and the axis-ratio $q$ was used to convert the grid
+to elliptical coordinates.
+
+From now on, we'll describe ellipticity using "elliptical components" $\epsilon_{1}$ and $\epsilon_{2}$, calculated
+from $\phi$ and $q$:
+
+$\epsilon_{1} = \frac{1 - q}{1 + q} \sin(2\phi)$
+$\epsilon_{2} = \frac{1 - q}{1 + q} \cos(2\phi)$
+
+We'll refer to these as `ell_comps` in the code for brevity.
+
+Future tutorials will explain why $\epsilon_{1}$ and $\epsilon_{2}$ are preferred over $q$ and $\phi$.
+
+*Exercise*: Try computing the elliptical components from the axis-ratio and angle above. What happens to the elliptical
+components when you adjust the `axis_ratio` and `angle_degrees` variables?
+
+
+```python
+fac = (1 - axis_ratio) / (1 + axis_ratio)
+epsilon_y = fac * np.sin(2 * np.radians(angle_degrees))
+epsilon_x = fac * np.cos(2 * np.radians(angle_degrees))
+
+ell_comps = (epsilon_y, epsilon_x)
+
+print("Elliptical Components:")
+print(ell_comps)
+
+```
+
+ Elliptical Components:
+ (np.float64(0.28867513459481287), np.float64(-0.16666666666666657))
+
+
+__Light Profiles__
+
+Galaxies are collections of stars, gas, dust, and other astronomical objects that emit light. Astronomers study this
+light to understand various properties of galaxies.
+
+To model the light of a galaxy, we use light profiles, which are mathematical functions that describe how a galaxy's
+light is distributed across space. By applying these light profiles to 2D grids of $(y, x)$ coordinates, we can
+create images that represent a galaxy's luminous emission.
+
+A commonly used light profile is the `Sersic` profile, which is widely adopted in astronomy for representing galaxy
+light. The `Sersic` profile is defined by the equation:
+
+$I_{\rm Ser} (\eta_{\rm l}) = I \exp \left\{ -k \left[ \left( \frac{\eta}{R} \right)^{\frac{1}{n}} - 1 \right] \right\}$
+
+In this equation:
+
+ - $\eta$ represents the elliptical coordinates of the profile in arc-seconds (refer to earlier sections for elliptical coordinates).
+ - $I$ is the intensity normalization of the profile, given in arbitrary units, which controls the overall brightness of the Sersic profile.
+ - $R$ is the effective radius in arc-seconds, which determines the size of the profile.
+ - $n$ is the Sersic index, which defines how 'steep' the profile is, influencing the concentration of light.
+ - $k$ is a constant that ensures half the light of the profile lies within the radius $R$, where $k = 2n - \frac{1}{3}$.
+
+We can evaluate this function using values for $(\eta, I, R, n)$ to calculate the intensity of the profile at
+a particular elliptical coordinate.
+
+
+```python
+elliptical_coordinate = (
+ 0.5 # The elliptical coordinate where we compute the intensity, in arc-seconds.
+)
+intensity = 1.0 # Intensity normalization of the profile in arbitrary units.
+effective_radius = 2.0 # Effective radius of the profile in arc-seconds.
+sersic_index = 1.0 # Sersic index of the profile.
+k = 2 * sersic_index - (
+ 1.0 / 3.0
+) # Calculating the constant k, note that this is an approximation.
+
+# Calculate the intensity of the Sersic profile at a specific elliptical coordinate.
+sersic_value = np.exp(
+ -k * ((elliptical_coordinate / effective_radius) ** (1.0 / sersic_index) - 1.0)
+)
+
+print("Intensity of Sersic Light Profile at Elliptical Coordinate 0.5:")
+print(sersic_value)
+```
+
+ Intensity of Sersic Light Profile at Elliptical Coordinate 0.5:
+ 3.4903429574618414
+
+
+The calculation above gives the intensity of the Sersic profile at an elliptical coordinate of 0.5.
+
+To create a complete image of the Sersic profile, we can evaluate the intensity at every point in our grid of
+elliptical coordinates.
+
+
+```python
+sersic_image = np.exp(-k * ((eta / effective_radius) ** (1.0 / sersic_index) - 1.0))
+```
+
+When we plot the resulting image, we can see how the properties of the grid affect its appearance:
+
+ - The peak intensity is at the position (0.3", 0.5"), where we shifted the grid.
+ - The image is elongated along a 60° counter-clockwise angle, corresponding to the rotation of the grid.
+ - The image has an elliptical shape, consistent with the axis ratio of 0.5.
+
+This demonstrates how the geometry of the grid directly influences the appearance of the light profile.
+
+*Exercise*: Try changing the values of `centre`, `ell_comps`, `effective_radius`, and `sersic_index` above.
+Observe how these adjustments change the Sersic profile image.
+
+
+```python
+aplt.plot_array(
+ array=aa.Array2D(values=sersic_image, mask=grid.mask), title="Sersic Image"
+)
+```
+
+
+
+
+
+
+
+Instead of manually handling these transformations, we can use `LightProfile` objects from the `light_profile`
+module (`lp`) for faster and more efficient calculations.
+
+Below, we define a `Sersic` light profile using the `Sersic` object. We can print the profile to display its parameters.
+
+
+```python
+sersic_light_profile = al.lp.Sersic(
+ centre=(0.0, 0.0),
+ ell_comps=(0.0, 0.1),
+ intensity=1.0,
+ effective_radius=2.0,
+ sersic_index=1.0,
+)
+
+print(sersic_light_profile)
+```
+
+ Sersic
+ centre: (0.0, 0.0)
+ ell_comps: (0.0, 0.1)
+ intensity: 1.0
+ effective_radius: 2.0
+ sersic_index: 1.0
+
+
+With this `Sersic` light profile, we can create an image by passing a grid to its `image_2d_from` method.
+
+The calculation will internally handle all the coordinate transformations and intensity evaluations we performed
+manually earlier, making it much simpler.
+
+The `Sersic` profile we created just above is different from the one we used to manually compute the image,
+so the image will look different. However, the process is the same.
+
+
+```python
+image = sersic_light_profile.image_2d_from(grid=grid)
+
+aplt.plot_array(array=image, title="Sersic Image via Light Profile")
+```
+
+
+
+
+
+
+
+The `image` is returned as an `Array2D` object. Similar to a `Grid2D`, it has two forms:
+
+ - `native`: A 2D array with shape [total_y_image_pixels, total_x_image_pixels].
+ - `slim`: A 1D array that flattens this data into shape [total_y_image_pixels * total_x_image_pixels].
+
+The `native` form is often used for visualizations, while the `slim` form can be useful for certain calculations.
+
+
+```python
+print("Intensity of pixel 0:")
+print(image.native[0, 0])
+print("Intensity of pixel 1:")
+print(image.slim[1])
+```
+
+ Intensity of pixel 0:
+ 0.01336649114209766
+ Intensity of pixel 1:
+ 0.014020623587576737
+
+
+To visualize the light profile's image, we use `aplt.plot_array`.
+
+We provide it with the light profile and the grid, which are used to create and plot the image.
+
+
+```python
+aplt.plot_array(
+ array=sersic_light_profile.image_2d_from(grid=grid), title="Image via Light Profile"
+)
+```
+
+
+
+
+
+
+
+__One Dimension Projection__
+
+We often want to calculative 1D quantities of a light profile, for example to plot how its light changes as
+a function of radius.
+
+To do this, we must still input a 2D grid into the `image_2d_from` method, therefore we create a project 2D
+radial grid as follows which has shape [Number_of_1d_coordinates, 2] and where all [:,0] entries are the same.
+
+A simple example of such a grid is as follows with 4 1D coordinates is:
+
+
+```python
+grid_2d_projected = al.Grid2DIrregular(
+ [
+ [1.000000e-06, 1.000000e-06],
+ [1.000000e-06, 1.000001e00],
+ [1.000000e-06, 2.000001e00],
+ [1.000000e-06, 3.000001e00],
+ ]
+)
+```
+
+As in this example, we often already have a 2D grid we are using to calculate images of a ligth profile
+and it would be convenient to simply create `grid_2d_projected` from that.
+
+For example, we may want the project grid which traces it major axis in uniform radial steps.
+
+This is easily computed using the `grid_2d_radial_project_from` function and passing the `centre` and `angle`
+of a light profile we can make it align with the light profile itself.
+
+Note how in this example the two galaxy bulges are not rotationally aligned but we aligned the projected
+grid with the first galaxy. The centres are aligned, but if they were not that would cause similar
+issues.
+
+
+```python
+grid_2d_projected = grid.grid_2d_radial_projected_from(
+ centre=sersic_light_profile.centre, angle=sersic_light_profile.angle()
+)
+
+image_1d = sersic_light_profile.image_2d_from(grid=grid_2d_projected)
+```
+
+We can now plot the 1D radial profile of the light profile. This profile shows how the intensity of the light
+changes as a function of distance from the profile's center. This is a more informative way to visualize the light p
+rofile's distribution.
+
+When we plot 1D quantities, we do not use built-in plotting functions as in 2D, but instead use standard
+matplotlib functionality.
+
+The reason is partly that 1D plotting is simple, but also because 1D plots have many different decisions
+about what is plotted and how they are computed, meaning its better to give the user full control.
+
+**Exercise**: Try plotting the 1D radial profile of Sersic profiles with different effective radii and Sersic indices.
+Does the 1D representation show more clearly how the light distribution changes with these parameters?
+
+
+```python
+plt.plot(grid_2d_projected[:, 1], image_1d)
+plt.xlabel("Radius (arcseconds)")
+plt.ylabel("Luminosity")
+plt.show()
+plt.close()
+
+```
+
+
+
+
+
+
+
+Since galaxy light distributions often cover a wide range of values, they are typically better visualized on a log10
+scale. This approach helps highlight details in the faint outskirts of a light profile.
+
+The `plot_array`/`subplot_\*` object has a `use_log10` option that applies this transformation automatically. Below, you can see
+that the image plotted in log10 space reveals more details.
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_15_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_15_0.png
new file mode 100644
index 0000000..a51d3a4
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_15_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_19_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_19_0.png
new file mode 100644
index 0000000..6524496
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_19_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_29_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_29_0.png
new file mode 100644
index 0000000..d8d2c2c
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_29_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_33_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_33_0.png
new file mode 100644
index 0000000..1047fe1
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_33_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_37_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_37_0.png
new file mode 100644
index 0000000..7393ae2
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_37_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_43_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_43_0.png
new file mode 100644
index 0000000..5945e1e
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_43_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_5_0.png b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_5_0.png
new file mode 100644
index 0000000..8487778
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_1_grids_and_galaxies_files/tutorial_1_grids_and_galaxies_5_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing.md b/markdown/chapter_1_introduction/tutorial_2_ray_tracing.md
new file mode 100644
index 0000000..bbe624c
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_2_ray_tracing.md
@@ -0,0 +1,298 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_2_ray_tracing.py`](../../scripts/chapter_1_introduction/tutorial_2_ray_tracing.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_2_ray_tracing.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_2_ray_tracing.ipynb).
+
+Tutorial 2: Ray Tracing
+=======================
+
+Strong gravitational lensing occurs when the mass of a foreground galaxy (or galaxies) curves space-time around it,
+causing light rays from a background source to appear deflected.
+
+The process of ray tracing calculates how much the path of each light ray is bent by the mass of the foreground galaxy.
+It then traces the paths of these light rays back to the observer, allowing us to determine how the source appears
+distorted. As a result, the images of the source may appear as arcs, multiple images, or other complex patterns.
+
+In the previous tutorial, we introduced **light profiles**, which are analytic functions that describe the distribution
+of light from a galaxy. In this tutorial, we will focus on **mass profiles**, which are analytic functions that
+describe the mass distribution within a galaxy. A key quantity derived from these profiles is the **deflection angle**,
+which quantifies how light is deflected by the mass of the galaxy at any point in space.
+
+Grids were crucial in the previous tutorial, as they enabled us to compute the light profile of a galaxy at every
+coordinate point. In ray tracing, grids are equally important because they are used to calculate the deflection
+angles of light rays caused by a mass profile and to map the source's light rays back to the observer.
+
+Let’s revisit the schematic of strong lensing:
+
+
+
+As observers, we do not see the true appearance of the source (e.g., a round blob of light). Instead, we only
+perceive its light after it has been deflected and lensed by the foreground galaxies. This resulting image is known
+as the **observed image** or **image-plane image**, which we will produce through the ray-tracing process.
+
+The schematic above uses the terms "image-plane" and "source-plane." In the context of gravitational lensing, a "plane"
+refers to a collection of galaxies located at the same redshift, meaning they are physically aligned parallel to one
+another.
+
+In this tutorial, we will create a strong lensing system consisting of planes similar to the one depicted above. While
+a plane can contain multiple galaxies, we will focus on a simple scenario with just one lens galaxy and one source
+galaxy.
+
+Here is an overview of what we'll cover in this tutorial:
+
+- **Grid**: How the 2D grids that were important for evaluating light profiles are equally important
+ for ray-tracing calculations.
+
+- **Mass Profiles**: Introduce mass profiles, which describe the mass distribution of galaxies and are used to
+ calculate deflection angles.
+
+- **Ray Tracing Grids**: How to map the light rays from the image-plane to the source-plane using the
+ deflection angles and **lens equation**.
+
+- **Ray Tracing Images**: How to evaluate the lensed image of a source galaxy after it has been
+ gravitationally lensed.
+
+- **Galaxies**: How to include both light and mass profiles in a single `Galaxy` object, and therefore
+construct realistic lens and source galaxies.
+
+- **Tracer**: Introduce the `Tracer` object, which automates the ray-tracing process and allows us to compute
+ images of the entire lens system.
+
+- **Mappings**: Visualize how image pixels map to the source plane and vice versa using the `lines=`/`positions=` overlays object.
+
+__Contents__
+
+- **Grid:** In the previous tutorial, we created 2D grids of (y,x) coordinates and showed how shifting and.
+- **Mass Profiles:** To perform lensing calculations, we use mass profiles available in the `mass_profile` module.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import matplotlib.pyplot as plt
+import autolens as al
+import autoarray as aa
+import autolens.plot as aplt
+
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Grid__
+
+In the previous tutorial, we created 2D grids of (y,x) coordinates and showed how shifting and rotating these grids
+is crucial for evaluating the light profiles of galaxies.
+
+Grids are also essential for performing ray-tracing calculations. The coordinates in the grid are deflected by the
+lens galaxy, and these deflected coordinates are used in ray-tracing calculations.
+
+Now, let’s create the grid for this tutorial, which we’ll call the `image_plane_grid`. It represents the grid of
+coordinates in the image plane, before the light is deflected by the lens galaxy. This grid is uniform, meaning
+every coordinate is evenly spaced. However, this uniformity will change after ray-tracing, as the light rays are
+mapped to the source plane.
+
+
+```python
+image_plane_grid = al.Grid2D.uniform(shape_native=(101, 101), pixel_scales=0.1)
+```
+
+__Mass Profiles__
+
+To perform lensing calculations, we use mass profiles available in the `mass_profile` module (accessible via `al.mp`).
+
+A mass profile is an analytic function that describes the mass distribution within a galaxy. It is used to calculate
+deflection angles and other quantities like surface density and gravitational potential.
+
+In gravitational lensing, deflection angles describe how a mass bends light by curving space-time.
+
+We will start with a simple mass profile, `IsothermalSph`, which represents a spherically symmetric isothermal
+mass distribution. This profile has two main properties: its center (the origin of the coordinate system) and its
+Einstein radius (which indicates the galaxy's mass and how much it bends light rays).
+
+
+```python
+sis_mass_profile = al.mp.IsothermalSph(
+ centre=(0.0, 0.0), # The (y,x) arc-second coordinates of the profile's center.
+ einstein_radius=1.6, # The Einstein radius of the profile in arc-seconds.
+)
+print(sis_mass_profile)
+```
+
+ IsothermalSph
+ centre: (0.0, 0.0)
+ ell_comps: (0.0, 0.0)
+ einstein_radius: 1.6
+ slope: 2.0
+ core_radius: 0.0
+
+
+In the previous tutorial, we used the `image_2d_from` method to compute the image of a light profile by evaluating
+its intensity at each (y,x) coordinate on the grid.
+
+Mass profiles have a similar method called `deflections_yx_2d_from`, which calculates the deflection angles at
+every (y,x) coordinate on the grid in units of arc-seconds.
+
+
+```python
+deflections = sis_mass_profile.deflections_yx_2d_from(grid=image_plane_grid)
+```
+
+Like grids and arrays, the deflection angles can be accessed using the `native` and `slim` attributes. These are
+structured similarly to a `Grid2D` object:
+
+- **native**: A 2D array with shape \([total_y_pixels, total_x_pixels, 2]\), where the last dimension represents the
+ (y,x) deflection components.
+
+- **slim**: A 1D array with shape \([total_y_pixels * total_x_pixels, 2]\), where the coordinates are flattened
+ into a single list.
+
+
+```python
+print("Deflection angles of pixel 0:")
+print(deflections.native[0, 0])
+print("Deflection angles of pixel 1:")
+print(deflections.slim[1])
+```
+
+ Deflection angles of pixel 0:
+ [ 1.13137085 -1.13137085]
+ Deflection angles of pixel 1:
+ [ 1.14274054 -1.11988573]
+
+
+There is an important difference between a grid and deflection angles. A `Grid2D` is a set of coordinates, while
+deflection angles are 2D vectors. This means that each deflection angle is defined at a specific (y,x) coordinate
+but has two components: a y and an x value, which vary across the grid.
+
+This is why the method is called `deflections_yx_2d_from`—the `yx` signifies that these are 2D vectors with both
+y and x components.
+
+The deflection angles are stored in a `VectorYX2D` data structure:
+
+
+```python
+print(type(deflections))
+```
+
+
+
+
+This structure includes a `grid`, which represents the `Grid2D` of coordinates where the deflection angles are
+calculated (in this case, the `image_plane_grid` we defined earlier). It also has vector-specific methods,
+such as `magnitude`, which calculates the magnitude of each deflection vector using \((x^2 + y^2)^{0.5}\).
+
+
+```python
+print("Deflection angle's `Grid2D` at pixel 0:")
+print(deflections.grid.native[0, 0])
+print("Deflection angle magnitude at pixel 0:")
+print(deflections.magnitudes.native[0, 0])
+```
+
+ Deflection angle's `Grid2D` at pixel 0:
+ [ 5. -5.]
+ Deflection angle magnitude at pixel 0:
+ 1.6
+
+
+We can use `aplt.plot_array` to visualize the deflection angles, which displays the y and x components separately.
+
+On this plot, you’ll see yellow and white lines called **critical curves**. These curves are important in lensing
+and will be explained in detail in the next tutorial.
+
+
+```python
+deflections = sis_mass_profile.deflections_yx_2d_from(grid=image_plane_grid)
+deflections_y = aa.Array2D(values=deflections.slim[:, 0], mask=image_plane_grid.mask)
+aplt.plot_array(array=deflections_y, title="Deflections Y")
+deflections = sis_mass_profile.deflections_yx_2d_from(grid=image_plane_grid)
+deflections_x = aa.Array2D(values=deflections.slim[:, 1], mask=image_plane_grid.mask)
+aplt.plot_array(array=deflections_x, title="Deflections X")
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+Mass profiles also have additional properties used in lensing calculations:
+
+- **convergence**: Represents the surface mass density of the profile in dimensionless units.
+- **potential**: Represents the "lensing potential" of the mass profile in dimensionless units.
+- **magnification**: Indicates how much brighter light rays appear due to the focusing effect of lensing.
+
+These quantities can be calculated using `*_from` methods and are returned as `Array2D` objects.
+
+
+```python
+convergence = sis_mass_profile.convergence_2d_from(grid=image_plane_grid)
+potential_2d = sis_mass_profile.potential_2d_from(grid=image_plane_grid)
+magnification_2d = al.LensCalc.from_mass_obj(
+ mass_obj=sis_mass_profile
+).magnification_2d_from(grid=image_plane_grid)
+```
+
+The same plotter API can be used to visualize these properties:
+
+
+```python
+aplt.plot_array(
+ array=sis_mass_profile.convergence_2d_from(grid=image_plane_grid),
+ title="Convergence",
+)
+aplt.plot_array(
+ array=sis_mass_profile.potential_2d_from(grid=image_plane_grid), title="Potential"
+)
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+One-dimensional plots can also be made using the same projection technique as in the previous tutorial:
+
+
+```python
+grid_2d_projected = image_plane_grid.grid_2d_radial_projected_from(
+ centre=sis_mass_profile.centre, angle=sis_mass_profile.angle()
+)
+
+convergence_1d = sis_mass_profile.convergence_2d_from(grid=grid_2d_projected)
+
+plt.plot(grid_2d_projected[:, 1], convergence_1d)
+plt.xlabel("Radius (arcseconds)")
+plt.ylabel("Luminosity")
+plt.show()
+plt.close()
+```
+
+
+
+
+
+
+
+The **convergence** and **potential** can be better understood when plotted in logarithmic space:
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_0.png b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_0.png
new file mode 100644
index 0000000..2acd4b7
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_1.png b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_1.png
new file mode 100644
index 0000000..f980698
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_15_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_0.png b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_0.png
new file mode 100644
index 0000000..fd284d6
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_1.png b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_1.png
new file mode 100644
index 0000000..c8858a0
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_19_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_21_0.png b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_21_0.png
new file mode 100644
index 0000000..c859b9f
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_2_ray_tracing_files/tutorial_2_ray_tracing_21_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing.md b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing.md
new file mode 100644
index 0000000..14aa99e
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing.md
@@ -0,0 +1,581 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_3_more_ray_tracing.py`](../../scripts/chapter_1_introduction/tutorial_3_more_ray_tracing.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_3_more_ray_tracing.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_3_more_ray_tracing.ipynb).
+
+Tutorial 5: More Ray Tracing
+============================
+
+We'll now reinforce the ideas that we learnt about ray-tracing in the previous tutorial and introduce the following
+new concepts:
+
+- What critical curves and caustics are.
+ - That by specifying redshifts and a cosmology, the results are converted from arc-second coordinates to physical
+ units of kiloparsecs (kpc). Again, if you're not an Astronomer, you may not be familiar with the unit of parsec, it
+ may be worth a quick Google!
+ - That a `Tracer` can be given any number of galaxies.
+
+Up to now, the planes have also had just one lens galaxy or source galaxy at a time. In this example, the tracer will
+have multiple galaxies at each redshift, meaning that each plane has more than one galaxy. In terms of lensing
+calculations:
+
+- If two or more lens galaxies are at the same redshift in the image-plane, the convergences, potentials and
+deflection angles of their mass profiles are summed when performing lensing calculations.
+
+- If two or more source galaxies are at the same redshift in the source-plane, their light can simply be summed before
+ray tracing.
+
+The `Tracer` fully accounts for this.
+
+__Contents__
+
+- **Initial Setup:** To begin, lets setup the grid we'll ray-trace using.
+- **Concise Code:** Lets set up the tracer used in the previous tutorial.
+- **Critical Curves:** To end, we can finally explain what the black lines that have appeared on many of the plots.
+- **Caustics:** In the previous tutorial, we plotted the critical curves of the mass profile on the image-plane.
+- **Units:** Lets plot the lensing quantities again.
+- **More Complexity:** We now make a lens with some attributes we didn`t in the last tutorial.
+- **Multi Galaxy Ray Tracing:** Now lets pass our 4 galaxies to a `Tracer`, which means the following will occur.
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import numpy as np
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Initial Setup__
+
+To begin, lets setup the grid we'll ray-trace using. But, lets do something crazy and use a higher resolution than
+the previous tutorials!
+
+Lets also stop calling it the `image_plane_grid`, and just remember from now on our `grid` is in the image-plane.
+
+
+```python
+grid = al.Grid2D.uniform(shape_native=(250, 250), pixel_scales=0.02)
+```
+
+The grid is now shape 250 x 250, which has more image-pixels than the 100 x 100 grid used previously.
+
+
+```python
+print(grid.shape_native)
+print(grid.shape_slim)
+```
+
+ (250, 250)
+ 62500
+
+
+__Concise Code__
+
+Lets set up the tracer used in the previous tutorial.
+
+Up to now, we have set up each profile one line at a time, making the code long and cumbersome to read.
+
+From here on, we'll set up galaxies in a single block of code, making it more concise and readable.
+
+
+```python
+lens = al.Galaxy(
+ redshift=0.5, mass=al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=1.6)
+)
+
+source = al.Galaxy(
+ redshift=1.0,
+ light=al.lp.SersicCoreSph(
+ centre=(0.0, 0.0),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=1.0,
+ radius_break=0.025,
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens, source])
+```
+
+__Critical Curves__
+
+To end, we can finally explain what the black lines that have appeared on many of the plots throughout this chapter
+actually are.
+
+These lines are called the 'critical curves', and they define line of infinite magnification due to a mass profile.
+They therefore mark where in the image-plane a mass profile perfectly `focuses` light rays such that if a source is
+located there, it will appear very bright: potentially 10-100x as brighter than its intrinsic luminosity.
+
+The outer white line is a `tangential_critical_curve`, because it describes how the image of the source galaxy is stretched
+tangentially. There is also an inner `radial_critical_curve` which appears in yellow on figures, which describes how the
+image of the source galaxy is stretched radially.
+
+However, a radial critical curve only appears when the lens galaxy's mass profile is shallower than isothermal (e.g.
+when its inner mass slope is less steep than a steep power-law). To make it appear below, we therefore change
+the mass profile of our lens galaxy to a `PowerLawSph` with a slope of 1.8.
+
+In the next tutorial, we'll introduce 'caustics', which are where the critical curves map too in the source-plane.
+
+
+```python
+mass_profile = al.mp.PowerLawSph(centre=(0.0, 0.0), einstein_radius=1.6, slope=1.8)
+
+lens = al.Galaxy(redshift=0.5, mass=mass_profile)
+
+tangential_critical_curve_list = al.LensCalc.from_mass_obj(
+ mass_obj=mass_profile
+).tangential_critical_curve_list_from(grid=grid)
+radial_critical_curves_list = al.LensCalc.from_mass_obj(
+ mass_obj=mass_profile
+).radial_critical_curve_list_from(grid=grid)
+
+```
+
+__Caustics__
+
+In the previous tutorial, we plotted the critical curves of the mass profile on the image-plane. We will now plot the
+'caustics', which correspond to each critical curve ray-traced to the source-plane. This is computed by using the
+lens galaxy mass profile's to calculate the deflection angles at the critical curves and ray-trace them to the
+source-plane.
+
+As discussed in the previous tutorial, critical curves mark regions of infinite magnification. Thus, if a source
+appears near a caustic in the source plane it will appear significantly brighter than its true luminosity.
+
+We again have to use a mass profile with a slope below 2.0 to ensure a radial critical curve and therefore radial
+caustic is formed. As above, the tangential critical curve is white and maps to the tangential caustic in the source-plane,
+which is also white. The radial critical curve is yellow and maps to the radial caustic, which is also yellow.
+
+We can plot both the tangential and radial critical curves and caustics using the `lines=`/`positions=` overlays object. The
+critical curves are plotted only for the image-plane grid, whereas the caustic in the source plane.
+
+
+```python
+sis_mass_profile = al.mp.IsothermalSph(centre=(0.0, 0.0), einstein_radius=1.6)
+
+tracer = al.Tracer(galaxies=[lens, source])
+
+tangential_critical_curve_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).tangential_critical_curve_list_from(grid=grid)
+radial_critical_curves_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).radial_critical_curve_list_from(grid=grid)
+
+
+traced_grid_list = tracer.traced_grid_2d_list_from(grid=grid)
+source_plane_grid = np.asarray(traced_grid_list[1])
+source_plane_grid = source_plane_grid[np.isfinite(source_plane_grid).all(axis=1)]
+
+aplt.plot_grid(grid=traced_grid_list[0], title="Plane 0 Grid")
+aplt.plot_grid(
+ grid=al.Grid2DIrregular(values=source_plane_grid),
+ title="Plane 1 Grid",
+)
+
+tangential_caustic_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).tangential_caustic_list_from(grid=grid)
+radial_caustics_list = al.LensCalc.from_tracer(tracer=tracer).radial_caustic_list_from(
+ grid=grid
+)
+
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+We can also plot the caustic on the source-plane image.
+
+
+```python
+aplt.plot_array(array=tracer.image_2d_list_from(grid=grid)[1], title="Plane 1 Image")
+```
+
+
+
+
+
+
+
+Caustics also mark the regions in the source-plane where the multiplicity of the strong lens changes. That is,
+if a source crosses a caustic, it goes from 2 images to 1 image. Try and show this yourself by changing the (y,x)
+centre of the source-plane galaxy's light profile!
+
+
+```python
+source = al.Galaxy(
+ redshift=1.0,
+ light=al.lp.SersicCoreSph(
+ centre=(0.0, 0.0),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=1.0,
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens, source])
+
+aplt.plot_array(array=tracer.image_2d_list_from(grid=grid)[1], title="Plane 1 Image")
+```
+
+
+
+
+
+
+
+__Units__
+
+Lets plot the lensing quantities again. However, we'll now use the `Units` object of the **PyAutoLens** plotter module
+to set `in_kpc=True` and therefore plot the y and x axes in kiloparsecs.
+
+This conversion is performed automatically, using the galaxy redshifts and cosmology.
+
+
+```python
+
+aplt.subplot_tracer(tracer=tracer, grid=grid)
+aplt.subplot_galaxies_images(tracer=tracer, grid=grid)
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+If you're too familiar with Cosmology, it will be unclear how exactly we converted the distance units from
+arcseconds to kiloparsecs. You'll need to read up on your Cosmology lecture to understand this properly.
+
+You can create a `Cosmology` object, which provides many methods for calculation different cosmological quantities,
+which are shown below (if you're not too familiar with cosmology don't worry that you don't know what these mean,
+it isn't massively important for using **PyAutoLens**).
+
+We will use a flat lambda CDM cosmology, which is the standard cosmological model often assumed in scientific studies.
+
+
+```python
+cosmology = al.cosmo.FlatLambdaCDM(H0=70, Om0=0.3)
+
+print("Image-plane arcsec-per-kpc:")
+print(cosmology.arcsec_per_kpc_from(redshift=0.5))
+print("Image-plane kpc-per-arcsec:")
+print(cosmology.kpc_per_arcsec_from(redshift=0.5))
+print("Angular Diameter Distance to Image-plane (kpc):")
+print(cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=0.5))
+
+print("Source-plane arcsec-per-kpc:")
+print(cosmology.arcsec_per_kpc_from(redshift=1.0))
+print("Source-plane kpc-per-arcsec:")
+print(cosmology.kpc_per_arcsec_from(redshift=1.0))
+print("Angular Diameter Distance to Source-plane:")
+print(cosmology.angular_diameter_distance_to_earth_in_kpc_from(redshift=1.0))
+
+print("Angular Diameter Distance From Image To Source Plane:")
+print(
+ cosmology.angular_diameter_distance_between_redshifts_in_kpc_from(
+ redshift_0=0.5, redshift_1=1.0
+ )
+)
+print("Lensing Critical convergence:")
+print(
+ cosmology.critical_surface_density_between_redshifts_solar_mass_per_kpc2_from(
+ redshift_0=0.5, redshift_1=1.0
+ )
+)
+```
+
+ Image-plane arcsec-per-kpc:
+ 0.1638289752817242
+ Image-plane kpc-per-arcsec:
+ 6.103926355398221
+ Angular Diameter Distance to Image-plane (kpc):
+ 1259025.187042171
+ Source-plane arcsec-per-kpc:
+ 0.12487552857882515
+ Source-plane kpc-per-arcsec:
+ 8.007974111346966
+ Angular Diameter Distance to Source-plane:
+ 1651763.228507974
+ Angular Diameter Distance From Image To Source Plane:
+ 707494.3382263458
+ Lensing Critical convergence:
+ 3083526795.6972556
+
+
+__More Complexity__
+
+We now make a lens with some attributes we didn`t in the last tutorial:
+
+ - A light profile representing a `bulge` of stars, meaning the lens galaxy's light will appear in the image for the
+ first time.
+ - An external shear, which accounts for the deflection of light due to line-of-sight structures.
+
+
+```python
+lens = al.Galaxy(
+ redshift=0.5,
+ bulge=al.lp.SersicSph(
+ centre=(0.0, 0.0), intensity=2.0, effective_radius=0.5, sersic_index=2.5
+ ),
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), ell_comps=(0.0, -0.111111), einstein_radius=1.6
+ ),
+ shear=al.mp.ExternalShear(gamma_1=0.05, gamma_2=0.0),
+)
+
+print(lens)
+```
+
+ Redshift: 0.5
+ Light Profiles:
+ SersicSph
+ centre: (0.0, 0.0)
+ ell_comps: (0.0, 0.0)
+ intensity: 2.0
+ effective_radius: 0.5
+ sersic_index: 2.5
+ Mass Profiles:
+ Isothermal
+ centre: (0.0, 0.0)
+ ell_comps: (0.0, -0.111111)
+ einstein_radius: 1.6
+ slope: 2.0
+ core_radius: 0.0
+ ExternalShear
+ centre: (0.0, 0.0)
+ ell_comps: (0.0, 0.0)
+ gamma_1: 0.05
+ gamma_2: 0.0
+
+
+Lets also create a small satellite galaxy nearby the lens galaxy and at the same redshift.
+
+
+```python
+lens_satellite = al.Galaxy(
+ redshift=0.5,
+ bulge=al.lp.DevVaucouleursSph(
+ centre=(1.0, 0.0), intensity=2.0, effective_radius=0.2
+ ),
+ mass=al.mp.IsothermalSph(centre=(1.0, 0.0), einstein_radius=0.4),
+)
+
+print(lens_satellite)
+```
+
+ Redshift: 0.5
+ Light Profiles:
+ DevVaucouleursSph
+ centre: (1.0, 0.0)
+ ell_comps: (0.0, 0.0)
+ intensity: 2.0
+ effective_radius: 0.2
+ sersic_index: 4.0
+ Mass Profiles:
+ IsothermalSph
+ centre: (1.0, 0.0)
+ ell_comps: (0.0, 0.0)
+ einstein_radius: 0.4
+ slope: 2.0
+ core_radius: 0.0
+
+
+Lets have a quick look at the appearance of our lens galaxy and its satellite.
+
+
+```python
+
+```
+
+And their deflection angles, noting that the satellite does not contribute as much to the deflections.
+
+
+```python
+
+
+# NOTE: In the new API, pass title directly as a string:
+# title="Lens Galaxy Deflections (x)"
+
+```
+
+Now, lets make two source galaxies at redshift 1.0. Instead of using the name `light` for the light profiles, lets
+instead use more descriptive names that indicate what morphological component of the galaxy the light profile
+represents. In this case, we'll use the terms `bulge` and `disk`, the two main structures that a galaxy can be made of
+
+
+```python
+source_0 = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.DevVaucouleursSph(
+ centre=(0.1, 0.2), intensity=0.3, effective_radius=0.3
+ ),
+ disk=al.lp.ExponentialCore(
+ centre=(0.1, 0.2),
+ ell_comps=(0.111111, 0.0),
+ intensity=3.0,
+ effective_radius=2.0,
+ ),
+)
+
+source_1 = al.Galaxy(
+ redshift=1.0,
+ disk=al.lp.ExponentialCore(
+ centre=(-0.3, -0.5),
+ ell_comps=(0.1, 0.0),
+ intensity=8.0,
+ effective_radius=1.0,
+ ),
+)
+
+print(source_0)
+print(source_1)
+```
+
+ Redshift: 1.0
+ Light Profiles:
+ DevVaucouleursSph
+ centre: (0.1, 0.2)
+ ell_comps: (0.0, 0.0)
+ intensity: 0.3
+ effective_radius: 0.3
+ sersic_index: 4.0
+ ExponentialCore
+ centre: (0.1, 0.2)
+ ell_comps: (0.111111, 0.0)
+ intensity: 3.0
+ effective_radius: 2.0
+ sersic_index: 1.0
+ radius_break: 0.01
+ alpha: 3.0
+ gamma: 0.25
+ Redshift: 1.0
+ Light Profiles:
+ ExponentialCore
+ centre: (-0.3, -0.5)
+ ell_comps: (0.1, 0.0)
+ intensity: 8.0
+ effective_radius: 1.0
+ sersic_index: 1.0
+ radius_break: 0.01
+ alpha: 3.0
+ gamma: 0.25
+
+
+Lets look at our source galaxies (before lensing)
+
+
+```python
+
+```
+
+__Multi Galaxy Ray Tracing__
+
+Now lets pass our 4 galaxies to a `Tracer`, which means the following will occur:
+
+ - Using the galaxy redshift`s, and image-plane and source-plane will be created each with two galaxies galaxies.
+
+We've also pass the tracer below a Planck15 cosmology, where the cosomology of the Universe describes exactly how
+ray-tracing is performed.
+
+
+```python
+tracer = al.Tracer(
+ galaxies=[lens, lens_satellite, source_0, source_1],
+ cosmology=al.cosmo.Planck15(),
+)
+```
+
+We can now plot the tracer`s image, which now there are two galaxies in each plane is computed as follows:
+
+ 1) First, using the image-plane grid, the images of the lens galaxy and its satellite are computed.
+
+ 2) Using the mass profiles of the lens and its satellite, their deflection angles are computed.
+
+ 3) These deflection angles are summed, such that the deflection of light due to the mass profiles of both galaxies in
+ the image-plane is accounted for.
+
+ 4) These deflection angles are used to trace every image-grid coordinate to the source-plane.
+
+ 5) The image of the source galaxies is computed by summing both of their images and ray-tracing their light back to
+ the image-plane.
+
+This process is pretty much the same as we have single in previous tutorials when there is one galaxy per plane. We
+are simply summing the images and deflection angles of the galaxies before using them to perform ray-tracing.
+
+
+```python
+aplt.plot_array(array=tracer.image_2d_from(grid=grid), title="Image")
+```
+
+
+
+
+
+
+
+As we did previously, we can plot the source plane grid to see how each coordinate was traced.
+
+
+```python
+aplt.plot_grid(grid=tracer.traced_grid_2d_list_from(grid=grid)[1], title="Plane 1 Grid")
+```
+
+
+
+
+
+
+
+We can zoom in on the source-plane to reveal the inner structure of the caustic.
+
+
+```python
+
+aplt.plot_grid(grid=tracer.traced_grid_2d_list_from(grid=grid)[1], title="Plane 1 Grid")
+```
+
+
+
+
+
+
+
+__Wrap Up__
+
+Tutorial 6 completed! Try the following:
+
+ 1) If you change the lens and source galaxy redshifts, does the tracer's image change?
+
+ 2) What happens to the cosmological quantities as you change these redshifts? Do you remember enough of your
+ cosmology lectures to predict how quantities like the angular diameter distance change as a function of redshift?
+
+ 3) The tracer has a small delay in being computed, whereas other tracers were almost instant. What do you think
+ is the cause of this slow-down?
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_0.png
new file mode 100644
index 0000000..96ec34a
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_1.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_1.png
new file mode 100644
index 0000000..7357030
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_11_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_13_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_13_0.png
new file mode 100644
index 0000000..c423256
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_13_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_15_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_15_0.png
new file mode 100644
index 0000000..3f1fadb
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_15_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_0.png
new file mode 100644
index 0000000..e474086
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_1.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_1.png
new file mode 100644
index 0000000..194e375
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_17_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_35_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_35_0.png
new file mode 100644
index 0000000..076478c
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_35_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_37_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_37_0.png
new file mode 100644
index 0000000..7e02c59
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_37_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_39_0.png b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_39_0.png
new file mode 100644
index 0000000..7e02c59
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_3_more_ray_tracing_files/tutorial_3_more_ray_tracing_39_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_4_point_sources.md b/markdown/chapter_1_introduction/tutorial_4_point_sources.md
new file mode 100644
index 0000000..d0c3761
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_4_point_sources.md
@@ -0,0 +1,38 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_4_point_sources.py`](../../scripts/chapter_1_introduction/tutorial_4_point_sources.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_4_point_sources.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_4_point_sources.ipynb).
+
+Tutorial 4: Point Sources
+=========================
+
+This tutorial is not wrriten yet, but will explain how point source lensing works.
+
+This tutorial is not necesary for using PyAutoLens or doing strong lens analysis, so don't worry that it is not
+written yet!
+
+Tutorial 8 summary is written and you should check that out instead!
+
+__Contents__
+
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Wrap Up__
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_5_lensing_formalism.md b/markdown/chapter_1_introduction/tutorial_5_lensing_formalism.md
new file mode 100644
index 0000000..93dcd00
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_5_lensing_formalism.md
@@ -0,0 +1,39 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_5_lensing_formalism.py`](../../scripts/chapter_1_introduction/tutorial_5_lensing_formalism.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_5_lensing_formalism.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_5_lensing_formalism.ipynb).
+
+Tutorial 5: Lensing Formalism
+=============================
+
+This tutorial is not wrriten yet, but will explain what all the different lens quantities are and give a more
+formal description of them.
+
+This tutorial is not necesary for using PyAutoLens or doing strong lens analysis, so don't worry that it is not
+written yet!
+
+Tutorial 8 summary is written and you should check that out instead!
+
+__Contents__
+
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Wrap Up__
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_6_data.md b/markdown/chapter_1_introduction/tutorial_6_data.md
new file mode 100644
index 0000000..62cac30
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_6_data.md
@@ -0,0 +1,540 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_6_data.py`](../../scripts/chapter_1_introduction/tutorial_6_data.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_6_data.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_6_data.ipynb).
+
+Tutorial 6: Data
+================
+
+In the last tutorials, we use tracers to create images of strong lenses. However, those images don't accurately
+represent what we would observe through a telescope.
+
+Real telescope images, like those taken with the Charge Coupled Device (CCD) imaging detectors on the Hubble Space
+Telescope, include several factors that affect what we see:
+
+**Telescope Optics:** The optical components of the telescope can blur the light, influencing the image's sharpness.
+
+**Exposure Time:** The time the detector collects light, affecting the clarity of the image. Longer exposure times
+gather more light, improving the signal-to-noise ratio and creating a clearer image.
+
+**Background Sky:** Light from the sky itself, such as distant stars or zodiacal light, adds noise to the image.
+adds additional noise to the image.
+
+In this tutorial, we'll simulate a strong lens image by applying these real-world effects to the light and mass
+profiles and images we created earlier.
+
+Here is an overview of what we'll cover in this tutorial:
+
+- **Optics Blurring:** We'll simulate how the telescope optics blur the galaxy's light, making the images appear blurred.
+- **Poisson Noise:** We'll add Poisson noise to the image, simulating the randomness in the photon-to-electron conversion process on the CCD.
+- **Background Sky:** We'll add a background sky to the image, simulating the light from the sky that adds noise to the image.
+- **Simulator:** We'll use the `SimulatorImaging` object to simulate imaging data that includes all these effects.
+
+__Contents__
+
+- **Initial Setup:** To create our simulated strong lens image, we first need a 2D grid.
+- **Optics Blurring:** All images captured using CCDs (like those on the Hubble Space Telescope or Euclid) experience some.
+- **Poisson Noise:** In addition to the blurring caused by telescope optics, we also need to consider Poisson noise when.
+- **Background Sky:** The final effect we will consider when simulating imaging data is the background sky.
+- **Simulator:** The `SimulatorImaging` object lets us create simulated imaging data while including the effects of.
+- **Output:** We will now save these simulated data to `.fits` files, the standard format used by astronomers for.
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+import numpy as np
+from pathlib import Path
+import autoarray as aa
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Initial Setup__
+
+To create our simulated strong lens image, we first need a 2D grid. This grid will represent the coordinate space over
+which we will simulate the strong lens's light distribution.
+
+
+```python
+grid = al.Grid2D.uniform(
+ shape_native=(
+ 101,
+ 101,
+ ), # The dimensions of the grid, which here is 100 x 100 pixels.
+ pixel_scales=0.1, # The conversion factor between pixel units and arc-seconds.
+)
+```
+
+Next, we define the properties of our strong lens. In this tutorial, we’ll represent the lens with no luminous
+emmission and an`Isothermal` mass profile. The source galaxy will be represented by a Sersic light profile.
+
+In the previous tutorial, the units of `intensity` were arbitrary. However, for this tutorial, where we simulate
+realistic imaging data, the intensity must have specific units. We’ll use units of electrons per second per pixel
+($e- pix^-1 s^-1$), which is standard for CCD imaging data.
+
+
+```python
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)
+ ),
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.Sersic(
+ centre=(0.1, 0.1),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0, # in units of e- pix^-1 s^-1
+ effective_radius=1.0,
+ sersic_index=2.5,
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+```
+
+Lets look at the tracer's image, which is the image we'll be simulating.
+
+
+```python
+aplt.plot_array(
+ array=tracer.image_2d_from(grid=grid), title="Tracer Image Before Simulating"
+)
+```
+
+
+
+
+
+
+
+__Optics Blurring__
+
+All images captured using CCDs (like those on the Hubble Space Telescope or Euclid) experience some level of blurring
+due to the optics of the telescope. This blurring occurs because the optical system spreads out the light from each
+point source (e.g., a star or a part of a galaxy).
+
+The Point Spread Function (PSF) describes how the telescope blurs the image. It can be thought of as a 2D representation
+of how a single point of light would appear in the image, spread out by the optics. In practice, the PSF is a 2D
+convolution kernel that we apply to the image to simulate this blurring effect.
+
+
+```python
+psf = al.Convolver.from_gaussian(
+ shape_native=(11, 11), # The size of the PSF kernel, represented as an 11x11 grid.
+ sigma=0.1, # Controls the width of the Gaussian PSF, which determines the level of blurring.
+ pixel_scales=grid.pixel_scales, # Maintains consistency with the scale of the image grid.
+ normalize=True, # Normalizes the PSF kernel so that its values sum to 1.
+)
+```
+
+We can visualize the PSF to better understand how it will blur the galaxy's image. The PSF is essentially a small
+image that represents the spreading out of light from a single point source. This kernel will be used to blur the
+entire tracer image when we perform the convolution.
+
+
+```python
+aplt.plot_array(array=psf.kernel, title="PSF 2D Kernel")
+```
+
+
+
+
+
+
+
+The PSF is often more informative when plotted on a log10 scale. This approach allows us to clearly observe values
+in its tail, which are much smaller than the central peak yet critical for many scientific analyses. The tail
+values may significantly affect the spread and detail captured in the data.
+
+
+```python
+aplt.plot_array(array=psf.kernel, title="PSF 2D Kernel")
+```
+
+
+
+
+
+
+
+Next, we'll manually perform a 2D convolution of the PSF with the image of the galaxy. This convolution simulates the
+blurring that occurs when the telescope optics spread out the galaxy's light.
+
+1. **Padding the Image**: Before convolution, we add padding (extra space with zero values) around the edges of the
+ image. This prevents unwanted edge effects when we perform the convolution, ensuring that the image's edges don't
+ become artificially altered by the process.
+
+2. **Convolution**: Using the `Convolver` object's `convolve` method, we apply the 2D PSF convolution to the padded
+ image. This step combines the PSF with the galaxy's light, simulating how the telescope spreads out the light.
+
+3. **Trimming the Image**: After convolution, we trim the padded areas back to their original size, obtaining a
+ convolved (blurred) image that matches the dimensions of the initial tracer image.
+
+
+```python
+image = tracer.image_2d_from(grid=grid) # The original unblurred image of the galaxy.
+padded_image = tracer.padded_image_2d_from(
+ grid=grid,
+ psf_shape_2d=psf.kernel.shape_native, # Adding padding based on the PSF size.
+)
+convolved_image = psf.convolved_image_from(
+ image=padded_image, blurring_image=None
+) # Applying the PSF convolution.
+blurred_image = convolved_image.trimmed_after_convolution_from(
+ kernel_shape=psf.kernel.shape_native
+) # Trimming back to the original size.
+```
+
+ .../PyAutoArray/autoarray/operators/convolver.py:1415: UserWarning: No blurring_image provided. Only the direct image will be convolved. This may change the correctness of the PSF convolution.
+ warnings.warn(
+
+
+We can now plot the original and the blurred images side by side. This allows us to clearly see how the PSF
+convolution affects the appearance of the galaxy, making the image appear softer and less sharp.
+
+
+```python
+aplt.plot_array(array=image, title="Tracer Image Before PSF")
+
+aplt.plot_array(array=blurred_image, title="")
+
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+__Poisson Noise__
+
+In addition to the blurring caused by telescope optics, we also need to consider Poisson noise when simulating imaging
+data.
+
+When a telescope captures an image of a galaxy, photons from the galaxy are collected by the telescope's mirror and
+directed onto a CCD (Charge-Coupled Device). The CCD is made up of a silicon lattice (or another material) that
+converts incoming photons into electrons. These electrons are then gathered into discrete squares, which form the
+pixels of the final image.
+
+The process of converting photons into electrons is inherently random, following a Poisson distribution. This randomness
+means that the number of electrons in each pixel can vary, even if the same number of photons hits the CCD. Therefore,
+the electron count per pixel becomes a Poisson random variable. For our simulation, this means that the recorded
+number of photons in each pixel will differ slightly from the true number due to this randomness.
+
+To replicate this effect in our simulation, we can add Poisson noise to the tracer image using NumPy’s random module,
+which generates values from a Poisson distribution.
+
+It's important to note that the blurring caused by the telescope optics occurs before the photons reach the CCD.
+Therefore, we need to add the Poisson noise after blurring the tracer image.
+
+We also need to consider the units of our image data. Let’s assume that the tracer image is measured in units of
+electrons per second ($e^- s^{-1}$), which is standard for CCD imaging data. To simulate the number of electrons
+actually detected in each pixel, we multiply the image by the observation’s exposure time. This conversion changes t
+he units to the total number of electrons collected per pixel over the entire exposure time.
+
+Once the image is converted, we add Poisson noise, simulating the randomness in the photon-to-electron conversion
+process. After adding the noise, we convert the image back to units of electrons per second for analysis, as
+this is the preferred unit for astronomers when studying their data.
+
+
+```python
+exposure_time = 300.0 # Units of seconds
+blurred_image_counts = (
+ blurred_image * exposure_time
+) # Convert to total electrons detected over the exposure time.
+blurred_image_with_poisson_noise = (
+ np.random.poisson(blurred_image_counts, blurred_image_counts.shape) / exposure_time
+) # Add Poisson noise and convert back to electrons per second.
+```
+
+Here is what the blurred image with Poisson noise looks like.
+
+
+```python
+aplt.plot_array(
+ array=aa.Array2D(values=blurred_image_with_poisson_noise, mask=blurred_image.mask),
+ title="Image With Poisson Noise",
+)
+```
+
+
+
+
+
+
+
+It is challenging to see the Poisson noise directly in the image above, as it is often subtle. To make the noise more
+visible, we can subtract the blurred image without Poisson noise from the one with noise.
+
+This subtraction yields the "Poisson noise realization" which highlights the variation in each pixel due to the Poisson
+distribution of photons hitting the CCD. It represents the noise values that were added to each pixel. We call
+it the realization because it is one possible outcome of the Poisson process, and the noise will be different each time
+we simulate the image.
+
+
+```python
+poisson_noise_realization = blurred_image_with_poisson_noise - blurred_image
+
+aplt.plot_array(
+ array=aa.Array2D(values=poisson_noise_realization, mask=blurred_image.mask),
+ title="Poisson Noise Realization",
+)
+```
+
+
+
+
+
+
+
+__Background Sky__
+
+The final effect we will consider when simulating imaging data is the background sky.
+
+In addition to light from the strong lens, the telescope also picks up light from the sky. This background sky light is
+primarily due to two sources: zodiacal light, which is light scattered by interplanetary dust in the solar system,
+and the unresolved emission from distant stars and tracer.
+
+For our simulation, we'll assume that the background sky has a uniform brightness across the image, measured at
+0.1 electrons per second per pixel. The background sky is added to the image before applying the PSF convolution
+and adding Poisson noise. This is important because it means that the background contributes additional noise to the
+image.
+
+The background sky introduces noise throughout the entire image, including areas where the strong lens is not present.
+This is why CCD images often appear noisy, especially in regions far from where the strong lens signal is detected.
+The sky noise can make it more challenging to observe faint details of the lens and source galaxies.
+
+To simulate this, we add a constant background sky to the tracer image and then apply Poisson noise to create the
+final simulated image as it would appear through a telescope.
+
+
+```python
+background_sky_level = 0.1
+
+# Add background sky to the blurred tracer image.
+blurred_image_with_sky = blurred_image + background_sky_level
+blurred_image_with_sky_counts = blurred_image_with_sky * exposure_time
+
+# Apply Poisson noise to the image with the background sky.
+blurred_image_with_sky_poisson_noise = (
+ np.random.poisson(
+ blurred_image_with_sky_counts, blurred_image_with_sky_counts.shape
+ )
+ / exposure_time
+)
+
+# Visualize the image with background sky and Poisson noise.
+aplt.plot_array(
+ array=aa.Array2D(
+ values=blurred_image_with_sky_poisson_noise, mask=blurred_image.mask
+ ),
+ title="Image With Background Sky",
+)
+
+# Create a noise map showing the differences between the blurred image with and without noise.
+poisson_noise_realization = (
+ blurred_image_with_sky_poisson_noise - blurred_image_with_sky
+)
+
+aplt.plot_array(
+ array=aa.Array2D(values=poisson_noise_realization, mask=blurred_image.mask),
+ title="Poisson Noise Realization",
+)
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+__Simulator__
+
+The `SimulatorImaging` object lets us create simulated imaging data while including the effects of PSF blurring,
+Poisson noise, and background sky all at once:
+
+
+```python
+simulator = al.SimulatorImaging(
+ exposure_time=300.0,
+ psf=psf,
+ background_sky_level=0.1,
+ add_poisson_noise_to_data=True,
+)
+
+dataset = simulator.via_tracer_from(tracer=tracer, grid=grid)
+```
+
+ .../PyAutoArray/autoarray/operators/convolver.py:1415: UserWarning: No blurring_image provided. Only the direct image will be convolved. This may change the correctness of the PSF convolution.
+ warnings.warn(
+
+
+By plotting the `data` from the dataset, we can see that it matches the image we simulated earlier. It includes
+the effects of PSF blurring, Poisson noise, and noise from the background sky. This image is a realistic
+approximation of what a telescope like the Hubble Space Telescope would capture.
+
+
+```python
+aplt.plot_array(array=dataset.data, title="Simulated Imaging Data")
+```
+
+
+
+
+
+
+
+The dataset also includes the `psf` (Point Spread Function) used to blur the strong lens image.
+
+For actual telescope data, the PSF is determined during data processing and is provided along with the observations.
+It's crucial for accurately deconvolving the PSF from the strong lens image, allowing us to recover the true properties
+of the strong lens. We'll explore this further in the next tutorial.
+
+
+```python
+aplt.plot_array(array=dataset.psf.kernel, title="Simulated PSF")
+```
+
+
+
+
+
+
+
+The dataset includes a `noise_map`, which represents the Root Mean Square (RMS) standard deviation of the noise
+estimated for each pixel in the image. Higher noise values mean that the measurements in those pixels are
+less certain, so those pixels are given less weight when analyzing the data.
+
+This `noise_map` is different from the Poisson noise arrays we plotted earlier. The Poisson noise arrays show the
+actual noise added to the image due to the random nature of photon-to-electron conversion on the CCD, as calculated
+using the numpy random module. These noise values are theoretical and cannot be directly measured in real telescope data.
+
+In contrast, the `noise_map` is our best estimate of the noise present in the image, derived from the data itself
+and used in the fitting process.
+
+
+```python
+aplt.plot_array(array=dataset.noise_map, title="Simulated Noise Map")
+```
+
+
+
+
+
+
+
+The `signal-to-noise_map` shows the ratio of the signal in each pixel to the noise level in that pixel. It is
+calculated by dividing the `data` by the `noise_map`.
+
+This ratio helps us understand how much of the observed signal is reliable compared to the noise, allowing us to
+see where we can trust the detected signal from the strong lens and where the noise is more significant.
+
+In general, a signal-to-noise ratio greater than 3 indicates that the signal is likely real and not overwhelmed by
+noise. For our datasets, the signal-to-noise ratio peaks at ~70, meaning we can trust the signal detected in the
+image.
+
+
+```python
+aplt.plot_array(array=dataset.signal_to_noise_map, title="Signal-To-Noise Map")
+```
+
+
+
+
+
+
+
+The `aplt.subplot_imaging_dataset` object can display all of these components together, making it a powerful tool for visualizing
+simulated imaging data.
+
+It also shows the Data and PSF on a logarithmic (log10) scale, which helps highlight the faint details in these
+components.
+
+The "Over Sampling" plots on the bottom of the figures display advanced features that can be ignored for now.
+
+
+```python
+aplt.subplot_imaging_dataset(dataset=dataset)
+```
+
+
+
+
+
+
+
+__Output__
+
+We will now save these simulated data to `.fits` files, the standard format used by astronomers for storing images.
+Most imaging data from telescopes like the Hubble Space Telescope (HST) are stored in this format.
+
+The `dataset_path` specifies where the data will be saved, in this case, in the directory
+`autolens_workspace/dataset/imaging/howtolens/`, which contains many example images distributed with
+the `autolens_workspace`.
+
+The files are named `data.fits`, `noise_map.fits`, and `psf.fits`, and will be used in the next tutorial.
+
+
+```python
+dataset_path = Path("dataset") / "imaging" / "howtolens"
+print("Dataset Path: ", dataset_path)
+
+aplt.fits_imaging(
+ dataset=dataset,
+ data_path=dataset_path / "data.fits",
+ noise_map_path=dataset_path / "noise_map.fits",
+ psf_path=dataset_path / "psf.fits",
+ overwrite=True,
+)
+```
+
+ Dataset Path: dataset/imaging/howtolens
+
+
+__Wrap Up__
+
+In this tutorial, you learned how CCD imaging data of a lens is collected using real telescopes like the
+Hubble Space Telescope, and how to simulate this data using the `SimulatorImaging` object.
+
+Let's summarise what we've covered:
+
+- **Optics Blurring**: The optics of a telescope blur the light from tracer, reducing the clarity and sharpness of
+the images.
+
+- **Poisson Noise**: The process of converting photons to electrons on a CCD introduces Poisson noise, which is random
+variability in the number of electrons collected in each pixel.
+
+- **Background Sky**: Light from the sky is captured along with light from the lens, adding a layer of noise across
+the entire image.
+
+- **Simulator**: The `SimulatorImaging` object enables us to simulate realistic imaging data by including all of
+these effects together and contains the `data`, `psf`, and `noise_map` components.
+
+- **Output**: We saved the simulated data to `.fits` files, the standard format used by astronomers for storing images.
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_11_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_11_0.png
new file mode 100644
index 0000000..eb6a698
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_11_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_13_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_13_0.png
new file mode 100644
index 0000000..eb6a698
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_13_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_0.png
new file mode 100644
index 0000000..2e8a101
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_1.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_1.png
new file mode 100644
index 0000000..cacc07e
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_17_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_21_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_21_0.png
new file mode 100644
index 0000000..0d7e951
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_21_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_23_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_23_0.png
new file mode 100644
index 0000000..41ba8fb
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_23_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_0.png
new file mode 100644
index 0000000..13a499f
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_1.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_1.png
new file mode 100644
index 0000000..bf23b9e
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_25_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_29_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_29_0.png
new file mode 100644
index 0000000..3b6b60b
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_29_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_31_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_31_0.png
new file mode 100644
index 0000000..0f6fa08
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_31_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_33_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_33_0.png
new file mode 100644
index 0000000..0d78811
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_33_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_35_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_35_0.png
new file mode 100644
index 0000000..38b6203
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_35_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_37_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_37_0.png
new file mode 100644
index 0000000..072cc53
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_37_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_7_0.png b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_7_0.png
new file mode 100644
index 0000000..5d76d26
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_6_data_files/tutorial_6_data_7_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting.md b/markdown/chapter_1_introduction/tutorial_7_fitting.md
new file mode 100644
index 0000000..23519e7
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_7_fitting.md
@@ -0,0 +1,1065 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_7_fitting.py`](../../scripts/chapter_1_introduction/tutorial_7_fitting.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_7_fitting.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_7_fitting.ipynb).
+
+Tutorial 7: Fitting
+===================
+
+In previous tutorials, we used light profiles to create simulated images of tracer and visualized how these images
+would appear when captured by a CCD detector on a telescope like the Hubble Space Telescope.
+
+However, this simulation process is the reverse of what astronomers typically do when analyzing real data. Usually,
+astronomers start with an observation—an actual image of a strong lens - and aim to infer detailed information about the
+lens’s properties, such as its mass and unlensed source properties.
+
+To achieve this, we must fit the observed image data with a model, identifying the combination of light and mass
+profiles that best matches the lens's appearance in the image. In this tutorial, we'll illustrate this process using
+the imaging data simulated in the previous tutorial. Our goal is to demonstrate how we can recover the parameters of
+the light profiles that we used to create the original simulation, as a proof of concept for the fitting procedure.
+
+The process of fitting data introduces essential statistical concepts like the `model`, `residual_map`, `chi-squared`,
+`likelihood`, and `noise_map`. These terms are crucial for understanding how fitting works, not only in astronomy but
+also in any scientific field that involves data modeling. This tutorial will provide a detailed introduction to these
+concepts and show how they are applied in practice to analyze astronomical data.
+
+Here is an overview of what we'll cover in this tutorial:
+
+- **Dataset**: Load the imaging dataset that we previously simulated, consisting of the image, noise map, and PSF.
+- **Mask**: Apply a mask to the data, excluding regions with low signal-to-noise ratios from the analysis.
+- **Masked Grid**: Create a masked grid, which contains only the coordinates of unmasked pixels, to evaluate the
+ galaxy's light profile in only unmasked regions.
+- **Fitting**: Fit the data with a galaxy model, computing key quantities like the model image, residuals,
+ chi-squared, and log likelihood to assess the quality of the fit.
+- **Bad Fits**: Demonstrate how even small deviations from the true parameters can significantly impact the fit.
+- **Model Fitting**: Perform a basic model fit on a simple dataset, adjusting the model parameters to improve the
+ fit quality.
+
+__Contents__
+
+- **Dataset & Mask:** Standard set up of the dataset and mask that is fitted.
+- **Masked Grid:** In tutorials 1 and 2, we emphasized that the `Grid2D` object is crucial for evaluating a lens's.
+- **Fitting:** Fit the lens model to the dataset and inspect the results.
+- **Incorrect Fit:** In the previous section, we successfully created and fitted a lens model to the image data.
+- **Model Fitting:** In the previous sections, we used the true model to fit the data, which resulted in a high log.
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import setup_notebook; setup_notebook()
+
+import numpy as np
+from pathlib import Path
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Dataset__
+
+We begin by loading the imaging dataset that we will use for fitting in this tutorial. This dataset is identical to the
+one we simulated in the previous tutorial, representing how a lens would appear if captured by a CCD camera.
+
+In the previous tutorial, we saved this dataset as .fits files in the `autolens_workspace/dataset/imaging/howtolens`
+folder. The `.fits` format is commonly used in astronomy for storing image data along with metadata, making it a
+standard for CCD imaging.
+
+The `dataset_path` below specifies where these files are located: `autolens_workspace/dataset/imaging/howtolens/`.
+
+
+```python
+dataset_path = Path("dataset") / "imaging" / "howtolens"
+```
+
+__Dataset Auto-Simulation__
+
+If the dataset does not already exist on your system, it will be created by running the corresponding
+simulator script. This ensures that all example scripts can be run without manually simulating data first.
+
+
+```python
+if not dataset_path.exists():
+ import subprocess
+ import sys
+
+ subprocess.run(
+ [sys.executable, "scripts/simulator/no_lens_light__mass_sis.py"],
+ check=True,
+ )
+
+dataset = al.Imaging.from_fits(
+ data_path=dataset_path / "data.fits",
+ noise_map_path=dataset_path / "noise_map.fits",
+ psf_path=dataset_path / "psf.fits",
+ pixel_scales=0.1,
+)
+```
+
+The `Imaging` object contains three key components: `data`, `noise_map`, and `psf`:
+
+- `data`: The actual image of the lens, which we will analyze.
+
+- `noise_map`: A map indicating the uncertainty or noise level in each pixel of the image, reflecting how much the
+ observed signal in each pixel might fluctuate due to instrumental or background noise.
+
+- `psf`: The Point Spread Function, which describes how a point source of light is spread out in the image by the
+ telescope's optics. It characterizes the blurring effect introduced by the instrument.
+
+Let's print some values from these components and plot a summary of the dataset to refresh our understanding of the
+imaging data.
+
+
+```python
+print("Value of first pixel in imaging data:")
+print(dataset.data.native[0, 0])
+print("Value of first pixel in noise map:")
+print(dataset.noise_map.native[0, 0])
+print("Value of first pixel in PSF:")
+print(dataset.psf.kernel.native[0, 0])
+
+aplt.subplot_imaging_dataset(dataset=dataset)
+```
+
+ Value of first pixel in imaging data:
+ 0.01999999999999999
+ Value of first pixel in noise map:
+ 0.02
+ Value of first pixel in PSF:
+ 0.0
+
+
+
+
+
+
+
+
+__Mask__
+
+The signal-to-noise map of the image highlights areas where the signal (light from the lens and source tracer)
+is detected above the background noise. Values above 3.0 indicate regions where the light is detected with a
+signal-to-noise ratio of at least 3, while values below 3.0 are dominated by noise, where the light is not
+clearly distinguishable.
+
+To ensure the fitting process focuses only on meaningful data, we typically mask out regions with low signal-to-noise
+ratios, removing areas dominated by noise from the analysis. This allows the fitting process to concentrate on the
+regions where the lens is clearly detected.
+
+Here, we create a `Mask2D` to exclude certain regions of the image from the analysis. The mask defines which parts of
+the image will be used during the fitting process.
+
+For our simulated image, a circular 3" mask centered at the center of the image is appropriate, since the simulated
+lens was positioned at the center.
+
+
+```python
+mask = al.Mask2D.circular(
+ shape_native=dataset.shape_native,
+ pixel_scales=dataset.pixel_scales,
+ radius=3.0, # The circular mask's radius in arc-seconds
+ centre=(0.0, 0.0), # center of the image which is also the center of the lens
+)
+
+print(mask) # 1 = True, meaning the pixel is masked. Edge pixels are indeed masked.
+print(mask[48:53, 48:53]) # Central pixels are `False` and therefore unmasked.
+```
+
+ Mask2D([[False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False]])
+ []
+
+
+We can visualize the mask over the strong lens image using an `aplt.subplot_imaging_dataset`, which helps us adjust the mask as needed.
+This is useful to ensure that the mask appropriately covers the lens and source light and does not exclude important
+regions.
+
+To overlay objects like a mask onto a figure, we use the `lines=`/`positions=` overlays object. This tool allows us to add custom
+visuals to any plot, providing flexibility in creating tailored visual representations.
+
+
+```python
+
+aplt.plot_array(array=dataset.data, title="Imaging Data With Mask")
+```
+
+
+
+
+
+
+
+Once we are satisfied with the mask, we apply it to the imaging data using the `apply_mask()` method. This ensures
+that only the unmasked regions are considered during the analysis.
+
+
+```python
+dataset = dataset.apply_mask(mask=mask)
+```
+
+ 2026-07-11 16:29:25,502 - autoarray.dataset.imaging.dataset - INFO - IMAGING - Data masked, contains a total of 225 image-pixels
+
+
+When we plot the masked imaging data again, the mask is now automatically included in the plot, even though we did
+not explicitly pass it using the `lines=`/`positions=` overlays object. The plot also zooms into the unmasked area, showing only the
+region where we will focus our analysis. This is particularly helpful when working with large images, as it centers
+the view on the regions where the strong lens's signal is detected.
+
+
+```python
+aplt.plot_array(array=dataset.data, title="Masked Imaging Data")
+```
+
+
+
+
+
+
+
+The mask is now stored as an additional attribute of the `Imaging` object, meaning it remains attached to the
+dataset. This makes it readily available when we pass the dataset to a `FitImaging` object for the fitting process.
+
+
+```python
+print("Mask2D:")
+print(dataset.mask)
+```
+
+ Mask2D:
+ Mask2D([[False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False],
+ [False, False, False, False, False, False, False, False, False,
+ False, False, False, False, False, False]])
+
+
+In earlier tutorials, we discussed how grids and arrays have `native` and `slim` representations:
+
+- `native`: Represents the original 2D shape of the data, maintaining the full pixel array of the image.
+- `slim`: Represents a 1D array containing only the values from unmasked pixels, allowing for more efficient
+ processing when working with large images.
+
+After applying the mask, the `native` and `slim` representations change as follows:
+
+- `native`: The 2D array keeps its original shape, [total_y_pixels, total_x_pixels], but masked pixels (those where
+ the mask is True) are set to 0.0.
+- `slim`: This now only contains the unmasked pixel values, reducing the array size
+ from [total_y_pixels * total_x_pixels] to just the number of unmasked pixels.
+
+Let's verify this by checking the shape of the data in its `slim` representation.
+
+
+```python
+print("Number of unmasked pixels:")
+print(dataset.data.native.shape)
+print(
+ dataset.data.slim.shape
+) # This should be lower than the total number of pixels, e.g., 100 x 100 = 10,000
+```
+
+ Number of unmasked pixels:
+ (15, 15)
+ (225,)
+
+
+The `mask` object also has a `pixels_in_mask` attribute, which gives the number of unmasked pixels. This should
+match the size of the `slim` data structure.
+
+
+```python
+print(dataset.data.mask.pixels_in_mask)
+```
+
+ 225
+
+
+We can use the `slim` attribute to print the first unmasked values from the image and noise map:
+
+
+```python
+print("First unmasked image value:")
+print(dataset.data.slim[0])
+print("First unmasked noise map value:")
+print(dataset.noise_map.slim[0])
+```
+
+ First unmasked image value:
+ 0.01999999999999999
+ First unmasked noise map value:
+ 0.02
+
+
+Additionally, we can verify that the `native` data structure has zeros at the edges where the mask is applied and
+retains non-zero values in the central unmasked regions.
+
+
+```python
+print("Example masked pixel in the image's native representation at its edge:")
+print(dataset.data.native[0, 0])
+print("Example unmasked pixel in the image's native representation at its center:")
+centre = tuple(s // 2 for s in dataset.data.shape_native)
+print(dataset.data.native[centre])
+```
+
+ Example masked pixel in the image's native representation at its edge:
+ 0.01999999999999999
+ Example unmasked pixel in the image's native representation at its center:
+ 0.6233333333333334
+
+
+__Masked Grid__
+
+In tutorials 1 and 2, we emphasized that the `Grid2D` object is crucial for evaluating a lens's light profile. This grid
+contains (y, x) coordinates for each pixel in the image and is used to ray-trace to the source plane and map out the
+positions where the source galaxy's light is calculated.
+
+From a `Mask2D`, we derive a `masked_grid`, which consists only of the coordinates of unmasked pixels. This ensures
+that light profile calculations focus exclusively on regions where the strong lens's light is detected, saving
+computational time and improving efficiency.
+
+Below, we plot the masked grid:
+
+
+```python
+masked_grid = mask.derive_grid.unmasked
+
+aplt.plot_grid(grid=masked_grid, title="Masked Grid2D")
+```
+
+
+
+
+
+
+
+By plotting this masked grid over the lens image, we can see that the grid aligns with the unmasked pixels of the
+image.
+
+This alignment **is crucial** for accurate fitting because it ensures that when we evaluate a strong lens's light
+profile, the calculations occur only at positions where we have real data from.
+
+
+```python
+aplt.plot_array(array=dataset.data, title="Image Data With 2D Grid Overlaid")
+```
+
+
+
+
+
+
+
+__Fitting__
+
+Now that our data is masked, we are ready to proceed with the fitting process.
+
+Fitting the data is done using the `Galaxy` and `Tracer objects that we introduced in previous tutorials. We will start by
+setting up a `Tracer`` object, using the same galaxy configuration that we previously used to simulate the
+imaging data. This setup will give us what is known as a 'perfect' fit, as the simulated and fitted models are identical.
+
+
+```python
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)
+ ),
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.Sersic(
+ centre=(0.1, 0.1),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=2.5,
+ ),
+)
+
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+
+aplt.plot_array(array=tracer.image_2d_from(grid=dataset.grid), title="Image")
+```
+
+
+
+
+
+
+
+Next, let's plot the image of the tracer. This should look familiar, as it is the same image we saw in
+previous tutorials. The difference now is that we use the dataset's `grid`, which corresponds to the `masked_grid`
+we defined earlier. This means that the tracer image is only evaluated in the unmasked region, skipping calculations
+in masked regions.
+
+
+```python
+aplt.plot_array(
+ array=tracer.image_2d_from(grid=dataset.grid), title="Tracer Image To Be Fitted"
+)
+```
+
+
+
+
+
+
+
+Now, we proceed to fit the image by passing both the `Imaging` and `Tracer` objects to a `FitImaging` object.
+This object will compute key quantities that describe the fit’s quality:
+
+`image`: Creates an image of the tracer using their image_2d_from() method.
+`model_data`: Convolves the tracer image with the data's PSF to account for the effects of telescope optics.
+`residual_map`: The difference between the model data and observed data.
+`normalized_residual_map`: Residuals divided by noise values, giving units of noise.
+`chi_squared_map`: Squares the normalized residuals.
+`chi_squared` and `log_likelihood`: Sums the chi-squared values to compute chi_squared, and converts this into
+a log_likelihood, which measures how well the model fits the data (higher values indicate a better fit).
+
+Let's create the fit and inspect each of these attributes:
+
+
+```python
+fit = al.FitImaging(dataset=dataset, tracer=tracer)
+```
+
+The `model_data` represents the tracer's image after accounting for effects like PSF convolution.
+
+An important technical note is that when we mask data, we discussed above how the image of the tracer is not evaluated
+outside the mask and is set to zero. This is a problem for PSF convolution, as the PSF blurs light from these regions
+outside the mask but at its edge into the mask. They must be correctly evaluated to ensure the model image accurately
+represents the image data.
+
+The `FitImaging` object handles this internally, but evaluating the model image in the additional regions outside the mask
+that are close enough to the mask edge to be blurred into the mask.
+
+
+```python
+print("First model image pixel:")
+print(fit.model_data.slim[0])
+aplt.plot_array(array=fit.model_data, title="Model Image")
+```
+
+ First model image pixel:
+
+
+ 1.8009715550322956
+
+
+
+
+
+
+
+
+Even before computing other fit quantities, we can normally assess if the fit is going to be good by visually comparing
+the `data` and `model_data` and assessing if they look similar.
+
+In this example, the tracer used to fit the data are the same as the tracer used to simulate it, so the two
+look very similar (the only difference is the noise in the image).
+
+
+```python
+aplt.plot_array(array=fit.data, title="Data")
+aplt.plot_array(array=fit.model_data, title="Model Image")
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+The `residual_map` is the different between the observed image and model image, showing where in the image the fit is
+good (e.g. low residuals) and where it is bad (e.g. high residuals).
+
+The expression for the residual map is simply:
+
+\[ \text{residual} = \text{data} - \text{model\_data} \]
+
+The residual-map is plotted below, noting that all values are very close to zero because the fit is near perfect.
+The only non-zero residuals are due to noise in the image.
+
+
+```python
+residual_map = dataset.data - fit.model_data
+print("First residual-map pixel:")
+print(residual_map.slim[0])
+
+print("First residual-map pixel via fit:")
+print(fit.residual_map.slim[0])
+
+aplt.plot_array(array=fit.residual_map, title="Residual Map")
+```
+
+ First residual-map pixel:
+ -1.7809715550322955
+ First residual-map pixel via fit:
+ -1.7809715550322955
+
+
+
+
+
+
+
+
+Are these residuals indicative of a good fit to the data? Without considering the noise in the data, it's difficult
+to ascertain. That is, its hard to ascenrtain if a residual value is large or small because this depends on the
+amount of noise in that pixel.
+
+The `normalized_residual_map` divides the residual-map by the noise-map, giving the residual in units of the noise.
+Its expression is:
+
+\[ \text{normalized\_residual} = \frac{\text{residual\_map}}{\text{noise\_map}} = \frac{\text{data} - \text{model\_data}}{\text{noise\_map}} \]
+
+If you're familiar with the concept of standard deviations (sigma) in statistics, the normalized residual map represents
+how many standard deviations the residual is from zero. For instance, a normalized residual of 2.0 (corresponding
+to a 95% confidence interval) means that the probability of the model underestimating the data by that amount is only 5%.
+
+
+```python
+normalized_residual_map = residual_map / dataset.noise_map
+
+print("First normalized residual-map pixel:")
+print(normalized_residual_map.slim[0])
+
+print("First normalized residual-map pixel via fit:")
+print(fit.normalized_residual_map.slim[0])
+
+aplt.plot_array(array=fit.normalized_residual_map, title="Normalized Residual Map")
+```
+
+ First normalized residual-map pixel:
+ -89.04857775161477
+ First normalized residual-map pixel via fit:
+ -89.04857775161477
+
+
+
+
+
+
+
+
+Next, we define the `chi_squared_map`, which is obtained by squaring the `normalized_residual_map` and serves as a
+measure of goodness of fit.
+
+The chi-squared map is calculated as:
+
+\[ \chi^2 = \left(\frac{\text{data} - \text{model\_data}}{\text{noise\_map}}\right)^2 \]
+
+Squaring the normalized residual map ensures all values are positive. For instance, both a normalized residual of -0.2
+and 0.2 would square to 0.04, indicating the same quality of fit in terms of `chi_squared`.
+
+As seen from the normalized residual map, it's evident that the model provides a good fit to the data, in this
+case because the chi-squared values are close to zero.
+
+
+```python
+chi_squared_map = (normalized_residual_map) ** 2
+print("First chi-squared pixel:")
+print(chi_squared_map.slim[0])
+
+print("First chi-squared pixel via fit:")
+print(fit.chi_squared_map.slim[0])
+
+aplt.plot_array(array=fit.chi_squared_map, title="Chi Squared Map")
+```
+
+ First chi-squared pixel:
+ 7929.649199585381
+ First chi-squared pixel via fit:
+ 7929.649199585381
+
+
+
+
+
+
+
+
+Now, we consolidate all the information in our `chi_squared_map` into a single measure of goodness-of-fit
+called `chi_squared`.
+
+It is defined as the sum of all values in the `chi_squared_map` and is computed as:
+
+\[ \chi^2 = \sum \left(\frac{\text{data} - \text{model\_data}}{\text{noise\_map}}\right)^2 \]
+
+This summing process highlights why ensuring all values in the chi-squared map are positive is crucial. If we
+didn't square the values (making them positive), positive and negative residuals would cancel each other out,
+leading to an inaccurate assessment of the model's fit to the data.
+
+The lower the `chi_squared`, the fewer residuals exist between the model's fit and the data, indicating a better
+overall fit!
+
+
+```python
+chi_squared = np.sum(chi_squared_map)
+print("Chi-squared = ", chi_squared)
+print("Chi-squared via fit = ", fit.chi_squared)
+```
+
+ Chi-squared = 645295.8813889123
+ Chi-squared via fit = 645295.8813889123
+
+
+The reduced chi-squared is the `chi_squared` value divided by the number of data points (e.g., the number of pixels
+in the mask).
+
+This quantity offers an intuitive measure of the goodness-of-fit, as it normalizes the `chi_squared` value by the
+number of data points. That is, a reduced chi-squared of 1.0 indicates that the model provides a good fit to the data,
+because every data point is fitted with a chi-squared value of 1.0.
+
+A reduced chi-squared value significantly greater than 1.0 indicates that the model is not a good fit to the data,
+whereas a value significantly less than 1.0 suggests that the model is overfitting the data.
+
+
+```python
+reduced_chi_squared = chi_squared / dataset.mask.pixels_in_mask
+print("Reduced Chi-squared = ", reduced_chi_squared)
+```
+
+ Reduced Chi-squared = 2867.9816950618324
+
+
+Another quantity that contributes to our final assessment of the goodness-of-fit is the `noise_normalization`.
+
+The `noise_normalization` is computed as the logarithm of the sum of squared noise values in our data:
+
+\[
+\text{{noise\_normalization}} = \sum \log(2 \pi \text{{noise\_map}}^2)
+\]
+
+This quantity is fixed because the noise-map remains constant throughout the fitting process. Despite this,
+including the `noise_normalization` is considered good practice due to its statistical significance.
+
+Understanding the exact meaning of `noise_normalization` isn't critical for our primary goal of successfully
+fitting a model to a dataset. Essentially, it provides a measure of how well the noise properties of our data align
+with a Gaussian distribution.
+
+
+```python
+noise_normalization = np.sum(np.log(2 * np.pi * dataset.noise_map**2))
+print("Noise Normalization = ", noise_normalization)
+print("Noise Normalization via fit = ", fit.noise_normalization)
+```
+
+ Noise Normalization = -1023.7670414997456
+ Noise Normalization via fit = -1023.7670414997456
+
+
+From the `chi_squared` and `noise_normalization`, we can define a final goodness-of-fit measure known as
+the `log_likelihood`.
+
+This measure is calculated by taking the sum of the `chi_squared` and `noise_normalization`, and then multiplying the
+result by -0.5:
+
+\[ \text{log\_likelihood} = -0.5 \times \left( \chi^2 + \text{noise\_normalization} \right) \]
+
+Don't worry about why we multiply by -0.5; it's a standard practice in statistics to ensure the log likelihood is
+defined correctly.
+
+
+```python
+log_likelihood = -0.5 * (chi_squared + noise_normalization)
+print("Log Likelihood = ", log_likelihood)
+print("Log Likelihood via fit = ", fit.log_likelihood)
+```
+
+ Log Likelihood = -322136.0571737063
+ Log Likelihood via fit = -322136.0571737063
+
+
+In the previous discussion, we noted that a lower \(\chi^2\) value indicates a better fit of the model to the
+observed data.
+
+When we calculate the log likelihood, we take the \(\chi^2\) value and multiply it by -0.5. This means that a
+higher log likelihood corresponds to a better model fit. Our goal when fitting models to data is to maximize the
+log likelihood.
+
+The **reduced \(\chi^2\)** value provides an intuitive measure of goodness-of-fit. Values close to 1.0 suggest a
+good fit, while values below or above 1.0 indicate potential underfitting or overfitting of the data, respectively.
+In contrast, the log likelihood values can be less intuitive. For instance, a log likelihood value printed above
+might be around 5300.
+
+However, log likelihoods become more meaningful when we compare them. For example, if we have two models, one with
+a log likelihood of 5300 and the other with 5310 we can conclude that the first model fits the data better
+because it has a higher log likelihood by 10.0.
+
+In fact, the difference in log likelihood between models can often be associated with a probability indicating how
+much better one model fits the data compared to another. This can be expressed in terms of standard deviations (sigma).
+
+As a rule of thumb:
+
+- A difference in log likelihood of **2.5** suggests that one model is preferred at the **2.0 sigma** level.
+- A difference in log likelihood of **5.0** indicates a preference at the **3.0 sigma** level.
+- A difference in log likelihood of **10.0** suggests a preference at the **5.0 sigma** level.
+
+All these metrics can be visualized together using the `aplt.subplot_fit_imaging` object, which offers a comprehensive
+overview of the fit quality. It also shows separate model images for the lens and source galaxies, and the appearance
+of the source galaxy in the image and source planes.
+
+
+```python
+fit = al.FitImaging(dataset=dataset, tracer=tracer)
+
+aplt.subplot_fit_imaging(fit=fit)
+```
+
+
+
+
+
+
+
+If you're familiar with model-fitting, you've likely encountered terms like 'residuals', 'chi-squared',
+and 'log_likelihood' before.
+
+These metrics are standard ways to quantify the quality of a model fit. They are applicable not only to 1D data but
+also to more complex data structures like 2D images, 3D data cubes, or any other multidimensional datasets.
+
+__Incorrect Fit___
+
+In the previous section, we successfully created and fitted a lens model to the image data, resulting in an
+excellent fit. The residual map and chi-squared map showed no significant discrepancies, indicating that the
+strong lens's light was accurately captured by our model. This optimal solution translates to one of the highest log
+likelihood values possible, reflecting a good match between the model and the observed data.
+
+Now, let's modify our lens model to create a fit that is close to the correct solution but slightly off.
+Specifically, we will slightly offset the center of the source galaxy by half a pixel (0.05") in both the x and y
+directions. This change will allow us to observe how even small deviations from the true parameters can impact the
+quality of the fit.
+
+
+```python
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)
+ ),
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.Sersic(
+ centre=(0.15, 0.15),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=2.5,
+ ),
+)
+
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+
+aplt.plot_array(array=tracer.image_2d_from(grid=dataset.grid), title="Image")
+```
+
+
+
+
+
+
+
+After implementing this slight adjustment, we can now plot the fit. In doing so, we observe that residuals have
+emerged at the multiple images of the lensed source, which indicates a mismatch between our model and the data.
+Consequently, this discrepancy results in increased chi-squared values, which in turn affects our log likelihood.
+
+
+```python
+fit_bad = al.FitImaging(dataset=dataset, tracer=tracer)
+
+aplt.subplot_fit_imaging(fit=fit_bad)
+```
+
+
+
+
+
+
+
+Next, we can compare the log likelihood of our current model to the log likelihood value we computed previously.
+
+
+```python
+print("Previous Likelihood:")
+print(fit.log_likelihood)
+print("New Likelihood:")
+print(fit_bad.log_likelihood)
+```
+
+ Previous Likelihood:
+ -322136.0571737063
+ New Likelihood:
+ -401645.6055797826
+
+
+As expected, we observe that the log likelihood has decreased! This decline confirms that our new model is indeed a
+worse fit to the data compared to the original model.
+
+Now, let’s change our lens model once more, this time setting it to a position that is far from the true parameters.
+We will offset the source's center significantly to see how this extreme deviation affects the fit quality.
+
+
+```python
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)
+ ),
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.Sersic(
+ centre=(0.5, 0.5),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=2.5,
+ ),
+)
+
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+
+fit_very_bad = al.FitImaging(dataset=dataset, tracer=tracer)
+
+aplt.subplot_fit_imaging(fit=fit_very_bad)
+```
+
+
+
+
+
+
+
+It is now evident that this model provides a terrible fit to the data. The tracer do not resemble a plausible
+representation of our simulated strong lens dataset, which we already anticipated given that we generated the data ourselves!
+
+As expected, the log likelihood has dropped dramatically with this poorly fitting model.
+
+
+```python
+print("Previous Likelihoods:")
+print(fit.log_likelihood)
+print(fit_bad.log_likelihood)
+print("New Likelihood:")
+print(fit_very_bad.log_likelihood)
+```
+
+ Previous Likelihoods:
+ -322136.0571737063
+ -401645.6055797826
+ New Likelihood:
+ -2703441.378104259
+
+
+__Model Fitting__
+
+In the previous sections, we used the true model to fit the data, which resulted in a high log likelihood and minimal
+residuals. We also demonstrated how even small deviations from the true parameters can significantly degrade the fit
+quality, reducing the log likelihood.
+
+In practice, however, we don't know the "true" model. For example, we might have an image of a strong lens observed with
+the Hubble Space Telescope, but the values for parameters like its `einstein_radius` and others are
+unknown. The process of determining the best-fit model is called model fitting, and it is the main topic of
+Chapter 2 of *HowToGalaxy*.
+
+To conclude this section, let's perform a basic, hands-on model fit to develop some intuition about how we can find
+the best-fit model. We'll start by loading a simple dataset that was simulated without any lens light, using
+an `IsothermalSph` lens mass profile and `ExponentialCoreSph` source light profile, but the true parameters of these
+profiles are unknown.
+
+
+```python
+dataset_name = "simple__no_lens_light__mass_sis"
+dataset_path = Path("dataset") / "imaging" / dataset_name
+
+dataset = al.Imaging.from_fits(
+ data_path=dataset_path / "data.fits",
+ psf_path=dataset_path / "psf.fits",
+ noise_map_path=dataset_path / "noise_map.fits",
+ pixel_scales=0.1,
+)
+
+mask = al.Mask2D.circular(
+ shape_native=dataset.shape_native,
+ pixel_scales=dataset.pixel_scales,
+ radius=3.0,
+)
+
+dataset = dataset.apply_mask(mask=mask)
+
+aplt.subplot_imaging_dataset(dataset=dataset)
+```
+
+ 2026-07-11 16:29:33,067 - autoarray.dataset.imaging.dataset - INFO - IMAGING - Data masked, contains a total of 225 image-pixels
+
+
+
+
+
+
+
+
+Now, you'll try to determine the best-fit model for this image, corresponding to the parameters used to simulate the
+dataset.
+
+We'll use the simplest possible approach: try different combinations of light and mass profile parameters and adjust
+them based on how well each model fits the data. You’ll quickly find that certain parameters produce a much better fit
+than others. For example, determining the correct values of the `centre` should not take too long.
+
+Pay attention to the `log_likelihood` and the `residual_map` as you adjust the parameters. These will guide you in
+determining if your model is providing a good fit to the data. Aim to increase the log likelihood and reduce the
+residuals.
+
+Keep experimenting with different values for a while, seeing how small you can make the residuals and how high you
+can push the log likelihood. Eventually, you’ll likely reach a point where further improvements become difficult,
+even after trying many different parameter values. This is a good point to stop and reflect on your first experience
+with model fitting.
+
+
+```python
+
+lens_galaxy = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.IsothermalSph(
+ centre=(1.0, 1.0), einstein_radius=1.0
+ ), # These are the lens parameters you need to adjust
+)
+
+source_galaxy = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.ExponentialCoreSph(
+ centre=(1.0, 1.0),
+ intensity=1.0,
+ effective_radius=1.0,
+ radius_break=0.025, # These are the source parameters you need to adjust
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens_galaxy, source_galaxy])
+
+fit = al.FitImaging(dataset=dataset, tracer=tracer)
+
+aplt.subplot_fit_imaging(fit=fit)
+
+print("Log Likelihood:")
+print(fit.log_likelihood)
+
+```
+
+
+
+
+
+
+
+ Log Likelihood:
+ -91864.03686663607
+
+
+Manually guessing model parameters repeatedly is a very inefficient and slow way to find the best fit. If the model
+were more complex—say, if the source galaxy had additional light profile components beyond just its `bulge` (like a
+second `Sersic` profile representing a `disk`)—the model would become so intricate that this manual approach
+would be practically impossible. This is definitely not how model fitting is done in practice.
+
+However, this exercise has given you a basic intuition for how model fitting works. The statistical inference tools
+that are actually used for model fitting will be introduced in Chapter 2. Interestingly, these tools are not entirely
+different from the approach you just tried. Essentially, they also involve iteratively testing models until those
+with high log likelihoods are found. The key difference is that a computer can perform this process thousands of
+times, and it does so in a much more efficient and strategic way.
+
+__Wrap Up__
+
+In this tutorial, you have learned how to fit a lens model to imaging data, a fundamental process in astronomy
+and statistical inference.
+
+Let's summarise what we have covered:
+
+- **Dataset**: We loaded the imaging dataset that we previously simulated, consisting of the tracer image, noise map,
+ and PSF.
+
+- **Mask**: We applied a circular mask to the data, excluding regions with low signal-to-noise ratios from the analysis.
+
+- **Masked Grid**: We created a masked grid, which contains only the coordinates of unmasked pixels, to evaluate the
+ tracer's light profile.
+
+- **Fitting**: We fitted the data with a lens model, computing key quantities like the model image, residuals,
+ chi-squared, and log likelihood to assess the quality of the fit.
+
+- **Bad Fits**: We demonstrated how even small deviations from the true parameters can significantly impact the fit
+ quality, leading to decreased log likelihood values.
+
+- **Model Fitting**: We performed a basic model fit on a simple dataset, adjusting the model parameters to improve the
+ fit quality.
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_11_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_11_0.png
new file mode 100644
index 0000000..ce921f7
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_11_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_15_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_15_0.png
new file mode 100644
index 0000000..2425188
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_15_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_27_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_27_0.png
new file mode 100644
index 0000000..e84c182
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_27_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_29_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_29_0.png
new file mode 100644
index 0000000..afa9d16
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_29_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_31_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_31_0.png
new file mode 100644
index 0000000..42ec776
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_31_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_33_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_33_0.png
new file mode 100644
index 0000000..c12c345
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_33_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_37_2.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_37_2.png
new file mode 100644
index 0000000..b51b224
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_37_2.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_0.png
new file mode 100644
index 0000000..e7649d6
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_1.png
new file mode 100644
index 0000000..b51b224
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_39_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_41_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_41_1.png
new file mode 100644
index 0000000..93e0681
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_41_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_43_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_43_1.png
new file mode 100644
index 0000000..80ab621
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_43_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_45_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_45_1.png
new file mode 100644
index 0000000..52e8560
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_45_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_55_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_55_0.png
new file mode 100644
index 0000000..e163d68
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_55_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_57_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_57_0.png
new file mode 100644
index 0000000..c607d80
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_57_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_59_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_59_0.png
new file mode 100644
index 0000000..1b5f28c
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_59_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_63_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_63_0.png
new file mode 100644
index 0000000..75c9523
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_63_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_67_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_67_1.png
new file mode 100644
index 0000000..553d39a
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_67_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_69_0.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_69_0.png
new file mode 100644
index 0000000..274acc1
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_69_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_7_1.png b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_7_1.png
new file mode 100644
index 0000000..707c207
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_7_fitting_files/tutorial_7_fitting_7_1.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_8_summary.md b/markdown/chapter_1_introduction/tutorial_8_summary.md
new file mode 100644
index 0000000..71d55da
--- /dev/null
+++ b/markdown/chapter_1_introduction/tutorial_8_summary.md
@@ -0,0 +1,284 @@
+> ✏️ **This page is auto-generated from [`scripts/chapter_1_introduction/tutorial_8_summary.py`](../../scripts/chapter_1_introduction/tutorial_8_summary.py) — do not edit it directly.**
+> It shows the example fully executed, with its real output images.
+> Run it yourself via the [Python script](../../scripts/chapter_1_introduction/tutorial_8_summary.py) or the [Jupyter notebook](../../notebooks/chapter_1_introduction/tutorial_8_summary.ipynb).
+
+Tutorial 9: Summary
+===================
+
+In this chapter, we have learnt that:
+
+ 1) **PyAutoLens** uses Cartesian `Grid2D`'s of $(y,x)$ coordinates to perform ray-tracing.
+ 2) These grids are combined with light and mass profiles to compute images, deflection angles and other quantities.
+ 3) Profiles are grouped together to make galaxies.
+ 4) Collections of galaxies (at the same redshift) form a plane.
+ 5) A `Tracer` can make an image-plane + source-plane strong lens system.
+ 6) The Universe's cosmology can be input into this `Tracer` to convert its units to kiloparsecs.
+ 7) The tracer's image can be used to simulate strong lens `Imaging` like it was observed with a real telescope.
+ 8) This data can be fitted, so to as quantify how well a model strong lens system represents the observed image.
+
+In this summary, we'll go over all the different Python objects introduced throughout this chapter and consider how
+they come together as one.
+
+__Contents__
+
+- **Start:** Below, we do all the steps we have learned this chapter, making profiles, galaxies, a tracer, etc.
+- **Object Composition:** Lets now consider how all of the objects we've covered throughout this chapter (`LightProfile`'s.
+- **Visualization:** Furthermore, using the `MatPLot2D` and `lines=`/`positions=` overlays objects we can visualize any.
+- **Code Design:** To end, I want to quickly talk about the **PyAutoLens** code-design and structure, which was really.
+- **Source Code:** If you do enjoy code, variables, functions, and parameters, you may want to dig deeper into the.
+- **Wrap Up:** Summary of the script and next steps.
+
+
+```python
+
+from autoconf import jax_wrapper # Sets JAX environment before other imports
+
+from autoconf import setup_notebook; setup_notebook()
+
+from pathlib import Path
+import autolens as al
+import autolens.plot as aplt
+```
+
+ Working Directory has been set to `HowToLens`
+
+
+__Start__
+
+Below, we do all the steps we have learned this chapter, making profiles, galaxies, a tracer, etc.
+
+Note that in this tutorial, we omit the lens galaxy's light and include two light profiles in the source representing a
+`bulge` and `disk`.
+
+
+```python
+grid = al.Grid2D.uniform(shape_native=(100, 100), pixel_scales=0.05)
+
+lens = al.Galaxy(
+ redshift=0.5,
+ mass=al.mp.Isothermal(
+ centre=(0.0, 0.0), einstein_radius=1.6, ell_comps=(0.17647, 0.0)
+ ),
+)
+
+source = al.Galaxy(
+ redshift=1.0,
+ bulge=al.lp.SersicCore(
+ centre=(0.1, 0.1),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=4.0,
+ ),
+ disk=al.lp.SersicCore(
+ centre=(0.1, 0.1),
+ ell_comps=(0.0, 0.111111),
+ intensity=1.0,
+ effective_radius=1.0,
+ sersic_index=1.0,
+ ),
+)
+
+tracer = al.Tracer(galaxies=[lens, source])
+```
+
+__Object Composition__
+
+Lets now consider how all of the objects we've covered throughout this chapter (`LightProfile`'s, `MassProfile`'s,
+`Galaxy`'s, `Plane`'s, `Tracer`'s) come together.
+
+The `Tracer`, which contains planes of `Galaxies`, which contains the `Galaxy`'s which contains the `Profile`'s:
+
+
+```python
+print(tracer)
+print()
+print(tracer.planes[0])
+print()
+print(tracer.planes[1])
+print()
+print(tracer.planes[0])
+print()
+print(tracer.planes[1])
+print()
+print(tracer.planes[0][0].mass)
+print()
+print(tracer.planes[1][0].bulge)
+print()
+print(tracer.planes[1][0].disk)
+print()
+```
+
+
+
+ [Redshift: 0.5
+ Mass Profiles:
+ Isothermal
+ centre: (0.0, 0.0)
+ ell_comps: (0.17647, 0.0)
+ einstein_radius: 1.6
+ slope: 2.0
+ core_radius: 0.0]
+ ... [60 lines of output truncated] ...
+ centre: (0.1, 0.1)
+ ell_comps: (0.0, 0.111111)
+ intensity: 1.0
+ effective_radius: 1.0
+ sersic_index: 4.0
+ radius_break: 0.025
+ alpha: 3.0
+ gamma: 0.25
+
+ SersicCore
+ centre: (0.1, 0.1)
+ ell_comps: (0.0, 0.111111)
+ intensity: 1.0
+ effective_radius: 1.0
+ sersic_index: 1.0
+ radius_break: 0.025
+ alpha: 3.0
+ gamma: 0.25
+
+
+
+Once we have a tracer we can therefore use any of the plotting function objects throughout this chapter to plot
+any specific aspect, whether it be a profile, galaxy, galaxies or tracer.
+
+For example, if we want to plot the image of the source galaxy's bulge and disk, we can do this in a variety of
+different ways.
+
+
+```python
+aplt.plot_array(array=tracer.image_2d_from(grid=grid), title="Image")
+
+source_plane_grid = tracer.traced_grid_2d_list_from(grid=grid)[1]
+
+```
+
+
+
+
+
+
+
+Understanding how these objects decompose into the different components of a strong lens is important for general
+**PyAutoLens** use.
+
+As the strong lens systems that we analyse become more complex, it is useful to know how to decompose their light
+profiles, mass profiles, galaxies and planes to extract different pieces of information about the strong lens. For
+example, we made our source-galaxy above with two light profiles, a `bulge` and `disk`. We can plot the lensed image of
+each component individually, now that we know how to break-up the different components of the tracer.
+
+
+```python
+aplt.plot_array(
+ array=tracer.planes[1][0].bulge.image_2d_from(grid=source_plane_grid),
+ title="Bulge Image",
+)
+
+aplt.plot_array(
+ array=tracer.planes[1][0].disk.image_2d_from(grid=source_plane_grid),
+ title="Disk Image",
+)
+```
+
+
+
+
+
+
+
+
+
+
+
+
+
+__Visualization__
+
+Furthermore, using the `MatPLot2D` and `lines=`/`positions=` overlays objects we can visualize any aspect we're interested
+in and fully customize the figure.
+
+Before beginning chapter 2 of **HowToLens**, you should checkout the package `autolens_workspace/plot`. This provides a
+full API reference of every plotting option in **PyAutoLens**, allowing you to create your own fully customized
+figures of strong lenses with minimal effort!
+
+
+```python
+
+tangential_critical_curve_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).tangential_critical_curve_list_from(grid=grid)
+radial_critical_curve_list = al.LensCalc.from_tracer(
+ tracer=tracer
+).radial_critical_curve_list_from(grid=grid)
+
+
+aplt.plot_array(
+ array=tracer.image_2d_from(grid=grid),
+ title="Tracer Image with Critical Curves",
+ lines=tangential_critical_curve_list,
+)
+```
+
+
+
+
+
+
+
+And, we're done, not just with the tutorial, but the chapter!
+
+__Code Design__
+
+To end, I want to quickly talk about the **PyAutoLens** code-design and structure, which was really the main topic of
+this tutorial.
+
+Throughout this chapter, we never talk about anything like it was code. We didn`t refer to 'variables', 'parameters`'
+'functions' or 'dictionaries', did we? Instead, we talked about 'galaxies', 'planes' a 'Tracer', etc. We discussed
+the objects that we, as scientists, think about when we consider a strong lens system.
+
+Software that abstracts the underlying code in this way follows an `object-oriented design`, and it is our hope
+with **PyAutoLens** that we've made its interface (often called the API for short) very intuitive, whether you were
+previous familiar with gravitational lensing or a complete newcomer!
+
+__Source Code__
+
+If you do enjoy code, variables, functions, and parameters, you may want to dig deeper into the **PyAutoLens** source
+code at some point in the future. Firstly, you should note that all of the code we discuss throughout the **HowToLens**
+lectures is not contained in just one project (e.g. the **PyAutoLens** GitHub repository) but in fact four repositories:
+
+**PyAutoFit** - Everything required for lens modeling (the topic of chapter 2): https://github.com/PyAutoLabs/PyAutoFit
+
+**PyAutoArray** - Handles all data structures and Astronomy dataset objects: https://github.com/PyAutoLabs/PyAutoArray
+
+**PyAutoGalaxy** - Contains the light profiles, mass profiles and galaxies: https://github.com/PyAutoLabs/PyAutoGalaxy
+
+**PyAutoLens** - Everything strong lensing: https://github.com/PyAutoLabs/PyAutoLens
+
+Instructions on how to build these projects from source are provided here:
+
+https://pyautolens.readthedocs.io/en/latest/installation/source.html
+
+We take a lot of pride in our source code, so I can promise you its well written, well documented and thoroughly
+tested (check out the `test` directory if you're curious how to test code well!).
+
+__Wrap Up__
+
+You`ve learn a lot in this chapter, but what you have not learnt is how to 'model' a real strong gravitational lens.
+
+In the real world, we have no idea what the 'correct' combination of light profiles, mass profiles and galaxies are
+that will give a good fit to a lens. Lens modeling is the process of finding the lens model which provides a good fit
+and it is the topic of chapter 2 of **HowToLens**.
+
+Finally, if you enjoyed doing the **HowToLens** tutorials please git us a star on the **PyAutoLens** GitHub
+repository:
+
+ https://github.com/PyAutoLabs/PyAutoLens
+
+Even the smallest bit of exposure via a GitHub star can help our project grow!
+
+
+```python
+
+```
diff --git a/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_11_0.png b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_11_0.png
new file mode 100644
index 0000000..97ce0e7
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_11_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_7_0.png b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_7_0.png
new file mode 100644
index 0000000..82f87ab
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_7_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_0.png b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_0.png
new file mode 100644
index 0000000..79b5f73
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_0.png differ
diff --git a/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_1.png b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_1.png
new file mode 100644
index 0000000..d80e430
Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_8_summary_files/tutorial_8_summary_9_1.png differ
diff --git a/notebooks/chapter_1_introduction/tutorial_7_fitting.ipynb b/notebooks/chapter_1_introduction/tutorial_7_fitting.ipynb
index d9b07ce..7e38d99 100644
--- a/notebooks/chapter_1_introduction/tutorial_7_fitting.ipynb
+++ b/notebooks/chapter_1_introduction/tutorial_7_fitting.ipynb
@@ -87,6 +87,8 @@
"cell_type": "code",
"metadata": {},
"source": [
+ "\n",
+ "from autoconf import setup_notebook; setup_notebook()\n",
"\n",
"import numpy as np\n",
"from pathlib import Path\n",
diff --git a/scripts/chapter_1_introduction/tutorial_7_fitting.py b/scripts/chapter_1_introduction/tutorial_7_fitting.py
index eb01303..23f6596 100644
--- a/scripts/chapter_1_introduction/tutorial_7_fitting.py
+++ b/scripts/chapter_1_introduction/tutorial_7_fitting.py
@@ -42,6 +42,8 @@
"""
+# from autoconf import setup_notebook; setup_notebook()
+
import numpy as np
from pathlib import Path
import autolens as al