diff --git a/tools/pymultispec.py b/tools/pymultispec.py
new file mode 100644
index 00000000..b86f7e66
--- /dev/null
+++ b/tools/pymultispec.py
@@ -0,0 +1,270 @@
+# -*- coding: utf-8 -*-
+"""readmultispec.py
+
+Read IRAF (echelle) spectrum in multispec format from a FITS file.
+Can read most multispec formats including linear, log, cubic spline,
+Chebyshev or Legendre dispersion spectra.
+
+Usage: retdict = readmultispec(fitsfile, reform=True)
+
+Inputs:
+
+fitfile     Name of the FITS file
+reform      If true (the default), a single spectrum dimensioned
+            [4,1,NWAVE] is returned as flux[4,NWAVE].  If false,
+            it is returned as a 3-D array flux[4,1,NWAVE].
+
+Returns a dictionary with these entries:
+flux        Array dimensioned [NCOMPONENTS,NORDERS,NWAVE] with the spectra.
+            If NORDERS=1, array is [NCOMPONENTS,NWAVE]; if NCOMPONENTS is also
+            unity, array is [NWAVE].  (This can be changed
+            using the reform keyword.)  Commonly the first dimension
+            is 4 and indexes the spectrum, an alternate version of
+            the spectrum, the sky, and the error array.  I have also
+            seen examples where NCOMPONENTS=2 (probably spectrum and
+            error).  Generally I think you can rely on the first element
+            flux[0] to be the extracted spectrum.  I don't know of
+            any foolproof way to figure out from the IRAF header what the
+            various components are.
+wavelen     Array dimensioned [NORDERS,NWAVE] with the wavelengths for
+            each order.
+header      The full FITS header from pyfits.
+wavefields  [NORDERS] List with the analytical wavelength
+            description (polynomial coefficients, etc.) extracted from
+            the header.  This is probably not very useful but is
+            included just in case.
+
+History:
+Created by Rick White based on my IDL readechelle.pro, 2012 August 15
+Apologies for any IDL-isms that remain!
+"""
+
+import numpy as np
+from astropy.io import fits as pyfits
+
+
+def nonlinearwave(nwave, specstr, verbose=False):
+    """Compute non-linear wavelengths from multispec string
+
+    Returns wavelength array and dispersion fields.
+    Raises a ValueError if it can't understand the dispersion string.
+    """
+
+    fields = specstr.split()
+    if int(fields[2]) != 2:
+        raise ValueError("Not nonlinear dispersion: dtype=" + fields[2])
+    if len(fields) < 12:
+        raise ValueError("Bad spectrum format (only %d fields)" % len(fields))
+    wt = float(fields[9])
+    w0 = float(fields[10])
+    ftype = int(fields[11])
+    if ftype == 3:
+
+        # cubic spline
+
+        if len(fields) < 15:
+            raise ValueError("Bad spline format (only %d fields)" % len(fields))
+        npieces = int(fields[12])
+        pmin = float(fields[13])
+        pmax = float(fields[14])
+        if verbose:
+            print("Dispersion is order-%d cubic spline" % npieces)
+        if len(fields) != 15 + npieces + 3:
+            raise ValueError(
+                "Bad order-%d spline format (%d fields)" % (npieces, len(fields))
+            )
+        coeff = np.asarray(fields[15:], dtype=float)
+        # normalized x coordinates
+        s = (np.arange(nwave, dtype=float) + 1 - pmin) / (pmax - pmin) * npieces
+        j = s.astype(int).clip(0, npieces - 1)
+        a = (j + 1) - s
+        b = s - j
+        x0 = a ** 3
+        x1 = 1 + 3 * a * (1 + a * b)
+        x2 = 1 + 3 * b * (1 + a * b)
+        x3 = b ** 3
+        wave = coeff[j] * x0 + coeff[j + 1] * x1 + coeff[j + 2] * x2 + coeff[j + 3] * x3
+
+    elif ftype == 1 or ftype == 2:
+
+        # chebyshev or legendre polynomial
+        # legendre not tested yet
+
+        if len(fields) < 15:
+            raise ValueError("Bad polynomial format (only %d fields)" % len(fields))
+        order = int(fields[12])
+        pmin = float(fields[13])
+        pmax = float(fields[14])
+        if verbose:
+            if ftype == 1:
+                print("Dispersion is order-%d Chebyshev polynomial" % order)
+            else:
+                print("Dispersion is order-%d Legendre polynomial (NEEDS TEST)" % order)
+        if len(fields) != 15 + order:
+            # raise ValueError('Bad order-%d polynomial format (%d fields)' % (order, len(fields)))
+            if verbose:
+                print(
+                    "Bad order-%d polynomial format (%d fields)" % (order, len(fields))
+                )
+                print("Changing order from %i to %i" % (order, len(fields) - 15))
+            order = len(fields) - 15
+        coeff = np.asarray(fields[15:], dtype=float)
+        # normalized x coordinates
+        pmiddle = (pmax + pmin) / 2
+        prange = pmax - pmin
+        x = (np.arange(nwave, dtype=float) + 1 - pmiddle) / (prange / 2)
+        p0 = np.ones(nwave, dtype=float)
+        p1 = x
+        wave = p0 * coeff[0] + p1 * coeff[1]
+        for i in range(2, order):
+            if ftype == 1:
+                # chebyshev
+                p2 = 2 * x * p1 - p0
+            else:
+                # legendre
+                p2 = ((2 * i - 1) * x * p1 - (i - 1) * p0) / i
+            wave = wave + p2 * coeff[i]
+            p0 = p1
+            p1 = p2
+
+    else:
+        raise ValueError("Cannot handle dispersion function of type %d" % ftype)
+
+    return wave, fields
+
+
+def readmultispec(fitsfile, reform=True, quiet=False):
+    """Read IRAF echelle spectrum in multispec format from a FITS file
+
+    Can read most multispec formats including linear, log, cubic spline,
+    Chebyshev or Legendre dispersion spectra
+
+    If reform is true, a single spectrum dimensioned 4,1,NWAVE is returned
+    as 4,NWAVE (this is the default.)  If reform is false, it is returned as
+    a 3-D array.
+    """
+
+    fh = pyfits.open(fitsfile)
+    try:
+        header = fh[0].header
+        flux = fh[0].data
+    finally:
+        fh.close()
+    temp = flux.shape
+    nwave = temp[-1]
+    if len(temp) == 1:
+        nspec = 1
+    else:
+        nspec = temp[-2]
+
+    # first try linear dispersion
+    try:
+        crval1 = header["crval1"]
+        crpix1 = header["crpix1"]
+        cd1_1 = header["cd1_1"]
+        ctype1 = header["ctype1"]
+        if ctype1.strip() == "LINEAR":
+            wavelen = np.zeros((nspec, nwave), dtype=float)
+            ww = (np.arange(nwave, dtype=float) + 1 - crpix1) * cd1_1 + crval1
+            for i in range(nspec):
+                wavelen[i, :] = ww
+            # handle log spacing too
+            dcflag = header.get("dc-flag", 0)
+            if dcflag == 1:
+                wavelen = 10.0 ** wavelen
+                if not quiet:
+                    print("Dispersion is linear in log wavelength")
+            elif dcflag == 0:
+                if not quiet:
+                    print("Dispersion is linear")
+            else:
+                raise ValueError("Dispersion not linear or log (DC-FLAG=%s)" % dcflag)
+
+            if nspec == 1 and reform:
+                # get rid of unity dimensions
+                flux = np.squeeze(flux)
+                wavelen.shape = (nwave,)
+            return {
+                "flux": flux,
+                "wavelen": wavelen,
+                "header": header,
+                "wavefields": None,
+            }
+    except KeyError:
+        pass
+
+    # get wavelength parameters from multispec keywords
+    try:
+        wat2 = header["wat2_*"]
+        count = len(wat2)
+    except KeyError:
+        raise ValueError("Cannot decipher header, need either WAT2_ or CRVAL keywords")
+
+    # concatenate them all together into one big string
+    watstr = []
+    for i in range(len(wat2)):
+        # hack to fix the fact that older pyfits versions (< 3.1)
+        # strip trailing blanks from string values in an apparently
+        # irrecoverable way
+        # v = wat2[i].value
+        v = wat2[i]
+        v = v + (" " * (68 - len(v)))  # restore trailing blanks
+        watstr.append(v)
+    watstr = "".join(watstr)
+
+    # find all the spec#="..." strings
+    specstr = [""] * nspec
+    for i in range(nspec):
+        sname = "spec" + str(i + 1)
+        p1 = watstr.find(sname)
+        p2 = watstr.find('"', p1)
+        p3 = watstr.find('"', p2 + 1)
+        if p1 < 0 or p1 < 0 or p3 < 0:
+            raise ValueError("Cannot find " + sname + " in WAT2_* keyword")
+        specstr[i] = watstr[p2 + 1 : p3]
+
+    wparms = np.zeros((nspec, 9), dtype=float)
+    w1 = np.zeros(9, dtype=float)
+    for i in range(nspec):
+        w1 = np.asarray(specstr[i].split(), dtype=float)
+        wparms[i, :] = w1[:9]
+        if w1[2] == -1:
+            raise ValueError(
+                "Spectrum %d has no wavelength calibration (type=%d)" % (i + 1, w1[2])
+            )
+            # elif w1[6] != 0:
+            #    raise ValueError('Spectrum %d has non-zero redshift (z=%f)' % (i+1,w1[6]))
+
+    wavelen = np.zeros((nspec, nwave), dtype=float)
+    wavefields = [None] * nspec
+    for i in range(nspec):
+        # if i in skipped_orders:
+        #    continue
+        verbose = (not quiet) and (i == 0)
+        if wparms[i, 2] == 0 or wparms[i, 2] == 1:
+            # simple linear or log spacing
+            wavelen[i, :] = np.arange(nwave, dtype=float) * wparms[i, 4] + wparms[i, 3]
+            if wparms[i, 2] == 1:
+                wavelen[i, :] = 10.0 ** wavelen[i, :]
+                if verbose:
+                    print("Dispersion is linear in log wavelength")
+            elif verbose:
+                print("Dispersion is linear")
+        else:
+            # non-linear wavelengths
+            wavelen[i, :], wavefields[i] = nonlinearwave(
+                nwave, specstr[i], verbose=verbose
+            )
+        wavelen *= 1.0 + wparms[i, 6]
+        if verbose:
+            print("Correcting for redshift: z=%f" % wparms[i, 6])
+    if nspec == 1 and reform:
+        # get rid of unity dimensions
+        flux = np.squeeze(flux)
+        wavelen.shape = (nwave,)
+    return {
+        "flux": flux,
+        "wavelen": wavelen,
+        "header": header,
+        "wavefields": wavefields,
+    }
diff --git a/tools/wavecal_creator_from_existing.py b/tools/wavecal_creator_from_existing.py
new file mode 100644
index 00000000..aa9cc16e
--- /dev/null
+++ b/tools/wavecal_creator_from_existing.py
@@ -0,0 +1,106 @@
+# -*- coding: utf-8 -*-
+"""
+This script creates a new linelist file based on an existing wavelength solution and an atlas for a specific element of the gas lamp
+used in the wavelength calibration.
+"""
+
+from os.path import dirname, join
+
+import matplotlib.pyplot as plt
+import numpy as np
+from pymultispec import readmultispec
+
+from pyreduce.configuration import get_configuration_for_instrument
+from pyreduce.instruments import instrument_info
+from pyreduce.wavelength_calibration import LineAtlas, LineList, WavelengthCalibration
+
+# Load Existing wavelength solution and ThAr file
+# Here from IRAF
+fname = join(dirname(__file__), "wef10006.ec.fits")
+spec = readmultispec(fname)
+wave = spec["wavelen"]
+flux = np.nansum(spec["flux"], axis=0)
+nord, ncol = flux.shape
+
+# Load the lineatlas for the used element
+# Here we use ThAr in vaccuum
+atlas = LineAtlas("thar", "vac")
+
+# Fill the linelist with the expected values
+linelist = {
+    "wlc": [],
+    "wll": [],
+    "posc": [],
+    "posm": [],
+    "xfirst": [],
+    "xlast": [],
+    "approx": [],
+    "width": [],
+    "height": [],
+    "order": [],
+    "flag": [],
+}
+# Adjust this threshold to what makes sense in your spectra
+# it is used to differentiate between background and line
+threshold = 1000
+
+for i in range(nord):
+    wmin, wmax = wave[i].min(), wave[i].max()
+    # Find only useful lines
+    idx_list = (atlas.linelist.wave > wmin) & (atlas.linelist.wave < wmax)
+    height = np.interp(atlas.linelist[idx_list].wave, wave[i], flux[i])
+    idx_list[idx_list] &= height > threshold
+    nlines = np.sum(idx_list)
+
+    linelist["wlc"] += [atlas.linelist[idx_list].wave]
+    linelist["wll"] += [atlas.linelist[idx_list].wave]
+    linelist["posc"] += [
+        np.interp(atlas.linelist[idx_list].wave, wave[i], np.linspace(0, ncol, ncol))
+    ]
+    linelist["posm"] += [
+        np.interp(atlas.linelist[idx_list].wave, wave[i], np.linspace(0, ncol, ncol))
+    ]
+    linelist["xfirst"] += [np.clip(linelist["posc"][-1] - 5, 0, ncol).astype(int)]
+    linelist["xlast"] += [np.clip(linelist["posc"][-1] + 5, 0, ncol).astype(int)]
+    linelist["approx"] += [np.full(nlines, "G")]
+    linelist["width"] += [np.full(nlines, 1, float)]
+    linelist["height"] += [
+        np.interp(atlas.linelist[idx_list].wave, wave[i], flux[i] / np.nanmax(flux[i]))
+    ]
+    linelist["order"] += [np.full(nlines, i)]
+    linelist["flag"] += [np.full(nlines, True)]
+
+# Combine the data from the different orders
+linelist = {k: np.concatenate(v) for k, v in linelist.items()}
+linelist = np.rec.fromarrays(
+    list(linelist.values()), names=list(linelist.keys()), dtype=LineList.dtype
+)
+
+# Run the lines through the wavecal
+# This updates the linelist inplace and flags bad ones
+# You need to check if this is a good solution
+# And update the settings in the config file accordingly
+# here: pyreduce/settings/settings_MCDONALD.json
+# or after the
+# config = get_configuration_for_instrument(instrument, plot=1)
+# line in your script
+
+# Setup the wavelength calibration module of PyReduce
+instrument = instrument_info.load_instrument("MCDONALD")
+module = WavelengthCalibration(
+    plot=1,
+    manual=True,
+    degree=8,
+    threshold=500,
+    iterations=3,
+    dimensionality="1D",
+    nstep=0,
+    shift_window=0.1,
+    element="thar",
+    medium="vac",
+)
+result = module.execute(flux, linelist)
+
+# Save the linelist
+linelist = LineList(linelist)
+linelist.save("mcdonald.npz")