diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index e5a3674..1c3d283 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -5,7 +5,9 @@ on:
   push:
     branches:
     - main
-
+  pull_request:
+    branches:
+    - main
 env:
   BUILD_NAME: 'github-actions[bot]'
   BUILD_EMAIL: '<41898282+github-actions[bot]@users.noreply.github.com>'
diff --git a/build-requirements.txt b/build-requirements.txt
index c93aa03..863d134 100644
--- a/build-requirements.txt
+++ b/build-requirements.txt
@@ -1,3 +1,4 @@
+fitsio
 numpy
 healpy
 camb
diff --git a/examples/advanced/DefineVariables.py b/examples/advanced/DefineVariables.py
new file mode 100644
index 0000000..61dfca1
--- /dev/null
+++ b/examples/advanced/DefineVariables.py
@@ -0,0 +1,90 @@
+'''
+Define Variables
+=======================
+
+To avoid large example files filled with variable definitions we define all relevant variables
+for the advanced exampels here
+
+'''
+
+
+import numpy as np
+import camb
+import glass.ext.camb
+from cosmology import Cosmology
+from glass.fields import generate_lognormal, lognormal_gls
+from glass.lensing import MultiPlaneConvergence
+from glass.shells import distance_grid, tophat_windows
+from glass.observations import gaussian_nz, smail_nz, equal_dens_zbins
+
+
+def get_common_data():
+    global lmax, nside, ncorr, rng, n_arcmin2
+    lmax = nside = 256
+    ncorr = 3
+    # photometric redshift distribution
+    gal_mean_z = [0.5, 1.0]
+    gal_sigma_z = 0.125
+    # total number of galaxies per arcmin2 in each bin
+    gal_dens = 2.0
+    # galaxy number density in units of galaxies/arcmin2/dz
+    z = np.linspace(0., 2., 100)
+    dndz = gaussian_nz(z, gal_mean_z, gal_sigma_z, norm=gal_dens)
+    # ellipticity distribution standard deviation
+    sigma_e = 0.26
+    # photometric redshift error
+    sigma_z0 = 0.03
+    # galaxy bias < 1 so that we can compute expectation from theory
+    beff = 0.9
+    nbin = len(dndz)
+    # labels for tomographic bins
+    binlabels = np.arange(nbin)
+    # name of the catalog we will create
+    catalogName = "myCatalog.FITS"
+    # random number generator, fix seed to make reproducible
+    rng = np.random.default_rng()
+    # galaxy density (using 1/100 of the expected galaxy number density for Stage-IV)
+    n_arcmin2 = 0.3
+    return lmax, nside, sigma_e, sigma_z0, beff, binlabels, catalogName, rng, n_arcmin2
+
+
+def get_glass_data(*, loadCls: bool = False):
+    # cosmology for the simulation
+    h = 0.7
+    Oc = 0.25
+    Ob = 0.05
+    # set up CAMB parameters for matter angular power spectrum
+    pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,
+                           NonLinear=camb.model.NonLinear_both)
+    # get the cosmology from CAMB
+    cosmo = Cosmology.from_camb(pars)
+    # shells of 200 Mpc in comoving distance spacing
+    zgrid = distance_grid(cosmo, 0., 3., dx=200.)
+    # shells from windows
+    shells = tophat_windows(zgrid)
+
+    ws = glass.shells.tophat_windows(zgrid, weight=glass.ext.camb.camb_tophat_weight)
+    if loadCls is True:
+        print("Loading Cls")
+        cls = np.load('../basic/cls.npy')
+    else:
+        print("Calculating Cls (Consider saving and loading them)")
+        cls = glass.ext.camb.matter_cls(pars, lmax, ws)
+
+    gls = lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=ncorr)
+    # assert len(gls) == len(shells) * (len(shells) + 1) // 2
+    # footprint
+
+    print("simulating fields using NSIDE = ", str(nside))
+    # generator for lognormal fields
+    matter = generate_lognormal(gls, nside, ncorr=ncorr, rng=rng)
+    # iterative computation of convergence
+    convergence = MultiPlaneConvergence(cosmo)
+    # true redshift distribution following a Smail distribution
+    z = np.arange(0., 3., 0.01)
+    dndz = smail_nz(z, z_mode=0.9, alpha=2., beta=1.5)
+    dndz *= n_arcmin2
+    # compute tomographic redshift bin edges with equal density
+    nbins = 10
+    zbins = equal_dens_zbins(z, dndz, nbins=nbins)
+    return pars, matter, shells, convergence, ws, z, dndz, zbins
diff --git a/examples/advanced/FITS_IO.py b/examples/advanced/FITS_IO.py
new file mode 100644
index 0000000..cbd4567
--- /dev/null
+++ b/examples/advanced/FITS_IO.py
@@ -0,0 +1,97 @@
+'''
+FITS I/O example
+================
+
+This example creates a galaxy catalogue from a Stage IV Space Satellite Galaxy
+Survey such as *Euclid* and *Roman* combining the :doc:`/basic/plot_density` and
+:doc:`/basic/plot_lensing` examples with galaxy ellipticities and galaxy shears,
+as well as using some auxiliary functions.
+
+The focus in this example is mock catalogue generation using auxiliary functions
+built for simulating Stage IV galaxy surveys.
+
+The catalog is then read and the redshifts plotted
+
+'''
+from itertools import count
+import numpy as np
+from glass.galaxies import galaxy_shear, redshifts_from_nz, gaussian_phz
+from glass.points import positions_from_delta
+from glass.shapes import ellipticity_intnorm
+from glass.lensing import from_convergence
+from glass.user import write_catalog
+from glass.shells import restrict
+from glass.observations import vmap_galactic_ecliptic, tomo_nz_gausserr
+from DefineVariables import get_common_data, get_glass_data
+import healpy as hp
+import matplotlib.pyplot as plt
+import fitsio
+
+
+# Import common variables from create_FITS_helper module
+lmax, nside, sigma_e, sigma_z0, beff, binlabels, catalogName, rng, n_arcmin2 = get_common_data()
+pars, matter, shells, convergence, ws, z, dndz, zbins = get_glass_data(loadCls=False)
+
+# Make a visibility map typical of a space telescope survey, seeing both
+# hemispheres, and low visibility in the galactic and ecliptic bands.
+vis = vmap_galactic_ecliptic(nside)
+# checking the mask:
+hp.mollview(vis, title='Stage IV Space Survey-like Mask', unit='Visibility')
+plt.show()
+
+print("Creating the catalog " + catalogName)
+# open a glass defined write context with a HDU extension name 'CATALOG'
+with write_catalog(catalogName, ext="CATALOG") as out:
+    # go through the shells and simulate
+    for i, delta, shell in zip(count(), matter, shells):
+        z_i, dndz_i = restrict(z, dndz, ws[i])
+        ngal = np.trapz(dndz_i, z_i)
+        print("shell ", i)
+        print("galaxy density: ", ngal)
+        convergence.add_window(delta, shell)
+        kappa = convergence.kappa
+        gamma, = from_convergence(kappa, lmax, shear=True)
+
+        for gal_lon, gal_lat, gal_count in positions_from_delta(ngal, delta, beff, vis):
+            print("galaxies sampled: ", gal_count)
+            gal_eps = ellipticity_intnorm(gal_count, sigma_e, rng=rng)
+            gal_she = galaxy_shear(gal_lon, gal_lat, gal_eps,
+                                   kappa, gamma.real, gamma.imag,
+                                   reduced_shear=False)
+            # some variance in shear weights
+            gal_w = 10**rng.uniform(-2, 2, size=gal_she.shape)
+            gal_z = redshifts_from_nz(gal_count, z_i, dndz_i)
+            # generator photometric redshifts using a Gaussian model
+            gal_phz = gaussian_phz(gal_z, sigma_z0)
+            # attach tomographic bin IDs to galaxies, based on photometric redshifts
+            gal_zbin = np.digitize(gal_phz, np.unique(zbins)) - 1
+
+            out.write(
+                RA=gal_lon,
+                DEC=gal_lat,
+                E1=gal_she.real,
+                E2=gal_she.imag,
+                W=gal_w,
+                Z_TRUE=gal_z,
+                PHZ=gal_phz,
+                ZBIN=gal_zbin,
+            )
+
+# split dndz using the same Gaussian error model assumed in the sampling
+tomo_nz = tomo_nz_gausserr(z, dndz, sigma_z0, zbins)
+nbins = 10
+print("Reading the catalog " + catalogName)
+with fitsio.FITS(catalogName, "r") as fits:
+    data = fits[1].read()
+    plt.figure()
+    plt.title('redshifts in catalogue')
+    plt.ylabel('dN/dz - normalised')
+    plt.xlabel('z')
+    for i in range(nbins):
+        in_bin = (data['ZBIN'] == i)
+        plt.hist(data['Z_TRUE'][in_bin], histtype='stepfilled', edgecolor='none', alpha=0.5, bins=50, density=1, label=f'cat. bin {i}')
+    for i in range(nbins):
+        plt.plot(z, (tomo_nz[i]/n_arcmin2)*nbins, alpha=0.5, label=f'inp. bin {i}')
+    plt.plot(z, dndz/n_arcmin2*nbins, ls='--', c='k')
+    plt.legend(ncol=2)
+    plt.show()
diff --git a/examples/advanced/plot_s4_galaxies.py b/examples/advanced/plot_s4_galaxies.py
deleted file mode 100644
index 2227cdd..0000000
--- a/examples/advanced/plot_s4_galaxies.py
+++ /dev/null
@@ -1,202 +0,0 @@
-'''
-Stage IV Galaxy Survey
-======================
-
-This example simulates a galaxy catalogue from a Stage IV Space Satellite Galaxy
-Survey such as *Euclid* and *Roman* combining the :doc:`/basic/plot_density` and
-:doc:`/basic/plot_lensing` examples with galaxy ellipticities and galaxy shears,
-as well as using some auxiliary functions.
-
-The focus in this example is mock catalogue generation using auxiliary functions
-built for simulating Stage IV galaxy surveys.
-
-'''
-
-# %%
-# Setup
-# -----
-# The setup is essentially the same as in the :doc:`/advanced/plot_shears`
-# example.
-#
-# In addition to a generator for intrinsic galaxy ellipticities,
-# following a normal distribution, we also show how to use auxiliary functions
-# to generate tomographic redshift distributions and visibility masks.
-#
-# Finally, there is a generator that applies the reduced shear from the lensing
-# maps to the intrinsic ellipticities, producing the galaxy shears.
-
-import numpy as np
-import healpy as hp
-import matplotlib.pyplot as plt
-
-# use the CAMB cosmology that generated the matter power spectra
-import camb
-from cosmology import Cosmology
-
-# GLASS modules: cosmology and everything in the glass namespace
-import glass.shells
-import glass.fields
-import glass.points
-import glass.shapes
-import glass.lensing
-import glass.galaxies
-import glass.observations
-import glass.ext.camb
-
-
-# cosmology for the simulation
-h = 0.7
-Oc = 0.25
-Ob = 0.05
-
-# basic parameters of the simulation
-nside = lmax = 256
-
-# set up CAMB parameters for matter angular power spectrum
-pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,
-                       NonLinear=camb.model.NonLinear_both)
-
-# get the cosmology from CAMB
-cosmo = Cosmology.from_camb(pars)
-
-# %%
-# Set up the matter sector.
-
-# shells of 200 Mpc in comoving distance spacing
-zb = glass.shells.distance_grid(cosmo, 0., 3., dx=200.)
-
-# tophat window functions for shells
-ws = glass.shells.tophat_windows(zb)
-
-# compute the angular matter power spectra of the shells with CAMB
-cls = glass.ext.camb.matter_cls(pars, lmax, ws)
-
-# compute Gaussian cls for lognormal fields for 3 correlated shells
-# putting nside here means that the HEALPix pixel window function is applied
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
-
-# generator for lognormal matter fields
-matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
-
-# %%
-# Set up the lensing sector.
-
-# this will compute the convergence field iteratively
-convergence = glass.lensing.MultiPlaneConvergence(cosmo)
-
-# %%
-# Set up the galaxies sector.
-
-# galaxy density (using 1/100 of the expected galaxy number density for Stage-IV)
-n_arcmin2 = 0.3
-
-# true redshift distribution following a Smail distribution
-z = np.arange(0., 3., 0.01)
-dndz = glass.observations.smail_nz(z, z_mode=0.9, alpha=2., beta=1.5)
-dndz *= n_arcmin2
-
-# compute tomographic redshift bin edges with equal density
-nbins = 10
-zbins = glass.observations.equal_dens_zbins(z, dndz, nbins=nbins)
-
-# photometric redshift error
-sigma_z0 = 0.03
-
-# constant bias parameter for all shells
-bias = 1.2
-
-# ellipticity standard deviation as expected for a Stage-IV survey
-sigma_e = 0.27
-
-
-# %%
-# Make a visibility map typical of a space telescope survey, seeing both
-# hemispheres, and low visibility in the galactic and ecliptic bands.
-vis = glass.observations.vmap_galactic_ecliptic(nside)
-
-# checking the mask:
-hp.mollview(vis, title='Stage IV Space Survey-like Mask', unit='Visibility')
-plt.show()
-
-
-# %%
-# Simulation
-# ----------
-# Simulate the galaxies with shears.  In each iteration, get the quantities of
-# interest to build our mock catalogue.
-
-# we will store the catalogue as a structured numpy array, initially empty
-catalogue = np.empty(0, dtype=[('RA', float), ('DEC', float),
-                               ('Z_TRUE', float), ('PHZ', float), ('ZBIN', int),
-                               ('G1', float), ('G2', float)])
-
-# simulate the matter fields in the main loop, and build up the catalogue
-for i, delta_i in enumerate(matter):
-
-    # compute the lensing maps for this shell
-    convergence.add_window(delta_i, ws[i])
-    kappa_i = convergence.kappa
-    gamm1_i, gamm2_i = glass.lensing.shear_from_convergence(kappa_i)
-
-    # the true galaxy distribution in this shell
-    z_i, dndz_i = glass.shells.restrict(z, dndz, ws[i])
-
-    # integrate dndz to get the total galaxy density in this shell
-    ngal = np.trapz(dndz_i, z_i)
-
-    # generate galaxy positions from the matter density contrast
-    for gal_lon, gal_lat, gal_count in glass.points.positions_from_delta(ngal, delta_i, bias, vis):
-
-        # generate random redshifts from the provided nz
-        gal_z = glass.galaxies.redshifts_from_nz(gal_count, z_i, dndz_i)
-
-        # generator photometric redshifts using a Gaussian model
-        gal_phz = glass.galaxies.gaussian_phz(gal_z, sigma_z0)
-
-        # attach tomographic bin IDs to galaxies, based on photometric redshifts
-        gal_zbin = np.digitize(gal_phz, np.unique(zbins)) - 1
-
-        # generate galaxy ellipticities from the chosen distribution
-        gal_eps = glass.shapes.ellipticity_intnorm(gal_count, sigma_e)
-
-        # apply the shear fields to the ellipticities
-        gal_she = glass.galaxies.galaxy_shear(gal_lon, gal_lat, gal_eps,
-                                              kappa_i, gamm1_i, gamm2_i)
-
-        # make a mini-catalogue for the new rows
-        rows = np.empty(gal_count, dtype=catalogue.dtype)
-        rows['RA'] = gal_lon
-        rows['DEC'] = gal_lat
-        rows['Z_TRUE'] = gal_z
-        rows['PHZ'] = gal_phz
-        rows['ZBIN'] = gal_zbin
-        rows['G1'] = gal_she.real
-        rows['G2'] = gal_she.imag
-
-        # add the new rows to the catalogue
-        catalogue = np.append(catalogue, rows)
-
-print(f'Total number of galaxies sampled: {len(catalogue):,}')
-
-# %%
-# Catalogue checks
-# ----------------
-# Here we can perform some simple checks at the catalogue level to
-# see how our simulation performed.
-
-# split dndz using the same Gaussian error model assumed in the sampling
-tomo_nz = glass.observations.tomo_nz_gausserr(z, dndz, sigma_z0, zbins)
-
-# redshift distribution of tomographic bins & input distributions
-plt.figure()
-plt.title('redshifts in catalogue')
-plt.ylabel('dN/dz - normalised')
-plt.xlabel('z')
-for i in range(nbins):
-    in_bin = (catalogue['ZBIN'] == i)
-    plt.hist(catalogue['Z_TRUE'][in_bin], histtype='stepfilled', edgecolor='none', alpha=0.5, bins=50, density=1, label=f'cat. bin {i}')
-for i in range(nbins):
-    plt.plot(z, (tomo_nz[i]/n_arcmin2)*nbins, alpha=0.5, label=f'inp. bin {i}')
-plt.plot(z, dndz/n_arcmin2*nbins, ls='--', c='k')
-plt.legend(ncol=2)
-plt.show()
diff --git a/examples/advanced/plot_shears.py b/examples/advanced/plot_shears.py
index 8ed3781..1b05be2 100644
--- a/examples/advanced/plot_shears.py
+++ b/examples/advanced/plot_shears.py
@@ -22,7 +22,6 @@
 
 # use the CAMB cosmology that generated the matter power spectra
 import camb
-from cosmology import Cosmology
 
 # GLASS modules: cosmology and everything in the glass namespace
 import glass.shells
@@ -32,49 +31,12 @@
 import glass.lensing
 import glass.galaxies
 from glass.core.constants import ARCMIN2_SPHERE
+from DefineVariables import get_common_data, get_glass_data
 
 
-# cosmology for the simulation
-h = 0.7
-Oc = 0.25
-Ob = 0.05
-
-# basic parameters of the simulation
-nside = lmax = 256
-
-# set up CAMB parameters for matter angular power spectrum
-pars = camb.set_params(H0=100*h, omch2=Oc*h**2, ombh2=Ob*h**2,
-                       NonLinear=camb.model.NonLinear_both)
-
-# get the cosmology from CAMB
-cosmo = Cosmology.from_camb(pars)
-
-# shells of 200 Mpc in comoving distance spacing
-zb = glass.shells.distance_grid(cosmo, 0., 1., dx=200.)
-
-# tophat window function for shells
-ws = glass.shells.tophat_windows(zb)
-
-# load the angular matter power spectra previously computed with CAMB
-cls = np.load('../basic/cls.npy')
-
-# %%
-# Matter
-# ------
-
-# compute Gaussian cls for lognormal fields for 3 correlated shells
-# putting nside here means that the HEALPix pixel window function is applied
-gls = glass.fields.lognormal_gls(cls, nside=nside, lmax=lmax, ncorr=3)
-
-# generator for lognormal matter fields
-matter = glass.fields.generate_lognormal(gls, nside, ncorr=3)
-
-# %%
-# Lensing
-# -------
-
-# this will compute the convergence field iteratively
-convergence = glass.lensing.MultiPlaneConvergence(cosmo)
+# Import common variables from create_FITS_helper module
+lmax, nside, sigma_e, sigma_z0, beff, binlabels, catalog, rng, n_arcmin2 = get_common_data()
+pars, matter, shells, convergence, ws, z, dndz, zbins = get_glass_data(loadCls=True)
 
 # %%
 # Galaxy
@@ -136,7 +98,7 @@
 
         # apply the shear fields to the ellipticities
         gal_she = glass.galaxies.galaxy_shear(gal_lon, gal_lat, gal_eps,
-                                            kappa_i, gamm1_i, gamm2_i)
+                                              kappa_i, gamm1_i, gamm2_i)
 
         # map the galaxy shears to a HEALPix map; this is opaque but works
         gal_pix = hp.ang2pix(nside, gal_lon, gal_lat, lonlat=True)