Working beta version.

README.md

@@ -1,32 +1,108 @@
 # tresmerge
 Pseudo-flux calibrate and merge TRES echelle spectra
 
-Requires Python 3.
+![Sample](sample.png)
 
-Required packages:
+
+This package requires Python 3 and the following packages:
 
 * numpy
 * scipy
 * matplotlib
 * astropy
 * astropy/specutils: https://github.com/astropy/specutils
 
+You can try installing these individually yourself beforehand, or the installation process should be able to install them all for you.
+
 ## Installation
 
 Download this package, `cd` to the directory, and run
 
-  $python setup.py install
+    $python setup.py install
 
 or
 
-  $python3 setup.py install
+    $python3 setup.py install
 
 Make sure that `python` or `python3` points to the Python 3 version you want to use. You can check that this is the correct version by first entering the Python interpreter, and you should see something like
 
-  $ python
-  Python 3.6.3 |Anaconda custom (64-bit)| (default, Nov  3 2017, 19:19:16)
-  [GCC 7.2.0] on linux
-  Type "help", "copyright", "credits" or "license" for more information.
-  >>>
+    $ python
+    Python 3.6.3 |Anaconda custom (64-bit)| (default, Nov  3 2017, 19:19:16)
+    [GCC 7.2.0] on linux
+    Type "help", "copyright", "credits" or "license" for more information.
+    >>>
 
 These scripts have not been tested with Python 2.x and will likely not work. It's recommended you upgrade to Python 3.
+
+The main script that you should use to run everything is called `tresmerge-process`, e.g.,
+
+    $ tresmerge-process --help
+    usage: tresmerge-process [-h] [--outfile OUTFILE] [--clobber] [--plot]
+                         [-t TRIM] [--shift SHIFT] [--poly-order POLY_ORDER]
+                         rfits bfits template
+
+    Merge TRES echelle orders.
+
+    positional arguments:
+      rfits                 Name of the FITS file containing the RAW spectrum to
+                            merge.
+      bfits                 Name of the FITS file containing the BLAZE-corrected
+                            spectrum to merge.
+      template              Name of the FITS file containing the Kurucz template
+                            spectrum.
+
+    optional arguments:
+      -h, --help            show this help message and exit
+      --outfile OUTFILE     Name of the output file to write the merged echelle
+                            spectrum to.
+      --clobber             Overwrite any existing output file with new data.
+      --plot                Make a set of plots of the merged spectra.
+      -t TRIM, --trim TRIM  How many pixels to trim from the front of the file.
+                            Default is 6
+      --shift SHIFT         Doppler shift the synthetic spectrum by this amount
+                            (in km/s) before doing the merge process. This may
+                            help if the target star has an exceptionally high
+                            radial velocity. Positive velocities correspond to
+                            redshifting the template.
+      --poly-order POLY_ORDER
+                            The order Chebyshev polynomial used to flatten each
+                            echelle order. 0 = constant, 1 = line, 2 = parabola, 3
+                            = cubic, ... etc.
+
+## Updating your package
+
+If you've previously installed this package, and there's been a new update, here is how you can proceed. If you used `git` to download this package, first `cd` to the repository and run
+
+    $ git pull
+
+This command will use git to bring down all of the new changes. Then, run
+
+    $python setup.py install
+
+again to install these changes to your system. If you didn't use `git` to download the package in the first place, you can delete the repository, download a new `.zip` package, and install it as before.
+
+## Citations
+
+If you make use of this code, please cite this repository as
+
+
+
+## How this package works
+
+This purpose of this code is to merge overlapping echelle orders of TRES spectrograph data into a single spectrum. Where there is no data in the redder echelle gaps, the code leaves no-data (despite what may appear as plotting artifacts).
+
+The data in each echelle order is compared to a synthetic template spectrum that the user provides. The mismatch is accounted for by fitting an `n`-th order (user specifies `n`) polynomial to the data and flattening it until it matches the template.
+
+The output is written to an enhanced-CSV file, which can be read by `astropy.io.ascii`, or any CSV reader. The header values in the original FITS files are copied as comments.
+
+### Usage
+
+For example, if you had two directories each containing the raw and blaze-corrected spectra, and a template, then you would run the code like this
+
+    $ tresmerge-process raw/hii468_2017-11-23_04h35m56s_cb.spec.fits blaze/hii468_2017-11-23_04h35m56s_cb.spec.fits t12000g40p00v180s.fits --plot
+
+If you enabled plotting (via `--plot`), then you should see an `output.pdf` with the spectra. You should also see an outfile containing the merged spectrum as a text file.
+
+## Notes about quality of merging
+
+For orders with deep absorption lines, the quality of fit of the synthetic spectra to the data makes a difference. Because the polynomial-fitting process seeks to minimize the distance between the synthetic spectra and data for each order, if the template has a number of inaccurate lines, then the polynomial fit will be biased. Future plans for this package may include adding in a masking option to do the polynomial fitting while excluding user-specified regions of the template, so please contact me or raise an issue on this repo if such a feature extension would be helpful.

tresmerge/common.py

@@ -1,26 +1,26 @@
 import numpy as np
 from scipy.interpolate import interp1d
 
-
+from scipy.linalg import cho_factor, cho_solve
 from numpy.polynomial import Chebyshev as Ch
 
 
-import pkg_resources
-
-# load the sensitivity curves into an easy-to-use interpolation object.
-sensfuncs_fname = pkg_resources.resource_filename("tresmerge", "data/sensfuncs.npy")
-wls, sens = np.load(sensfuncs_fname)
+# import pkg_resources
+#
+# # load the sensitivity curves into an easy-to-use interpolation object.
+# sensfuncs_fname = pkg_resources.resource_filename("tresmerge", "data/sensfuncs.npy")
+# wls, sens = np.load(sensfuncs_fname)
 # the sensfuncs is a (2,51,2998) shaped object, whereby the first is the wavelengths, and the second is the sensitivity.
 
 # the number of echelle orders in TRES
-norders = 51
-
-# create order-by-order interpolation objects
-interps = [interp1d(wls[i], sens[i], kind="linear", fill_value="extrapolate", assume_sorted=True) for i in range(norders)]
+# norders = 51
+#
+# # create order-by-order interpolation objects
+# interps = [interp1d(wls[i], sens[i], kind="linear", fill_value="extrapolate", assume_sorted=True) for i in range(norders)]
 
 
 
-def optimize_calibration_static(wl_cal, fl_cal, sigma_cal, wl_fixed, fl_fixed, sigma_fixed, amp, l_f, order=1, mu_GP=1.0):
+def solve_flux(wl, fl, sigma, fl_synth, order=3):
     '''
     Determine the calibration parameters for this epoch of observations. Assumes all wl, fl arrays are 1D, and that the relative velocities between all epochs are zero.
 
@@ -33,8 +33,6 @@ def optimize_calibration_static(wl_cal, fl_cal, sigma_cal, wl_fixed, fl_fixed, s
         wl_fixed (np.array) : the 1D (flattened) array of the reference wavelengths
         fl_fixed (np.array) : the 1D (flattened) array of the reference fluxes
         sigma_fixed (np.array) : the 1D (flattened) array of the reference sigmas
-        amp (float): the GP amplitude
-        l_f (float): the GP length
         order (int): the Chebyshev order to use
         mu_GP (optional): the mean of the GP to assume.
 
@@ -45,90 +43,32 @@ def optimize_calibration_static(wl_cal, fl_cal, sigma_cal, wl_fixed, fl_fixed, s
     wl0 = np.min(wl)
     wl1 = np.max(wl)
 
-    N_A = len(wl)
-    A = np.empty((N_A, N_A), dtype=np.float64)
-
-    N_B = len(wl_fixed)
-    B = np.empty((N_B, N_B), dtype=np.float64)
-
-    C = np.empty((N_A, N_B), dtype=np.float64)
-
-    matrix_functions.fill_V11_f(A, wl_cal, amp, l_f)
-    matrix_functions.fill_V11_f(B, wl_fixed, amp, l_f)
-    matrix_functions.fill_V12_f(C, wl_cal, wl_fixed, amp, l_f)
-
-    # Add in sigmas
-    A[np.diag_indices_from(A)] += sigma_cal**2
-    B[np.diag_indices_from(B)] += sigma_fixed**2
-
-
     # Get a clean set of the Chebyshev polynomials evaluated on the input wavelengths
     T = []
     for i in range(0, order + 1):
         coeff = [0 for j in range(i)] + [1]
         Chtemp = Ch(coeff, domain=[wl0, wl1])
-        Ttemp = Chtemp(wl_cal)
+        Ttemp = Chtemp(wl)
         T += [Ttemp]
 
     T = np.array(T)
 
-    D = fl_cal[:,np.newaxis] * T.T
-
-
-    # Solve for the calibration coefficients c0, c1
+    # AKA, dcal
+    D = fl[:,np.newaxis] * T.T
 
-    # Find B^{-1}, fl_prime, and C_prime
-    try:
-        B_cho = cho_factor(B)
-    except np.linalg.linalg.LinAlgError:
-        print("Failed to solve matrix inverse. Calibration not valid.")
-        raise
-
-    fl_prime = mu_GP + np.dot(C, cho_solve(B_cho, (fl_fixed.flatten() - mu_GP)))
-    C_prime = A - np.dot(C, cho_solve(B_cho, C.T))
-
-    # Find {C^\prime}^{-1}
-    CP_cho = cho_factor(C_prime)
+    # Inverse covariance matrix for data
+    C_inv = np.diag(1/sigma**2)
 
     # Invert the least squares problem
-    left = np.dot(D.T, cho_solve(CP_cho, D))
-    right = np.dot(D.T, cho_solve(CP_cho, fl_prime))
+    left = np.dot(D.T, np.dot(C_inv, D))
+    right = np.dot(D.T, np.dot(C_inv, fl_synth))
 
     left_cho = cho_factor(left)
 
     # the coefficents, X = [c0, c1]
     X = cho_solve(left_cho, right)
 
-    # Apply the correction
+    # Apply the correction to the data
     fl_cor = np.dot(D, X)
 
     return fl_cor, X
-
-
-            # optimize the calibration of "tweak" with respect to all other orders
-            fl_cor, X = covariance.optimize_calibration_ST1(lwl0, lwl1, lwl_tweak, fl_tweak, fl_remain[indsort], gp, A, C, order=config["order_cal"])
-
-            # since fl_cor may have actually have fewer pixels than originally, we can't just
-            # stuff the corrected fluxes directly back into the array.
-            # Instead, we need to re-evaluate the line on all the wavelengths
-            # (including those that may have been masked)
-            # using the chebyshev coefficients, and apply this.
-
-            # Here is where we need to make sure that wl0 and wl1 are the same.
-            T = []
-            for k in range(0, config["order_cal"] + 1):
-                pass
-                coeff = [0 for j in range(k)] + [1]
-                Chtemp = Ch(coeff, domain=[lwl0, lwl1])
-                Ttemp = Chtemp(lwl_tweak_unmasked)
-                T += [Ttemp]
-
-            T = np.array(T)
-
-            Q = fl_tweak_unmasked[:,np.newaxis] * T.T
-
-            # Apply the correction
-            fl_cor = np.dot(Q, X)
-
-            # replace this epoch with the corrected fluxes
-            fl_out[i] = fl_cor[:,0]