From e7f733ccbd414ffbf2d2fba4ba6656a5477d91c2 Mon Sep 17 00:00:00 2001
From: lastephey <lastephey@github.com>
Date: Tue, 19 Jun 2018 09:51:54 -0700
Subject: [PATCH] at least runs to completion, probably doesn't get the right
 answer yet

---
 py/specter/psf/psf.py | 32 ++++++++++++++++++++++----------
 1 file changed, 22 insertions(+), 10 deletions(-)

diff --git a/py/specter/psf/psf.py b/py/specter/psf/psf.py
index 2ea3ce1..e04deca 100644
--- a/py/specter/psf/psf.py
+++ b/py/specter/psf/psf.py
@@ -77,7 +77,7 @@ def __init__(self, filename):
         #- Filled only if needed
         self._xsigma = None
         self._ysigma = None
-        
+       
     #- Utility function to fit spot sigma vs. wavelength
     def _fit_spot_sigma(self, ispec, axis=0, npoly=5):
         """
@@ -85,7 +85,7 @@ def _fit_spot_sigma(self, ispec, axis=0, npoly=5):
         Return callable Legendre object.
         
         Inputs:
-            ispec : spectrum number
+           ispec : spectrum number
             axis  : 0 or 'x' for cross dispersion sigma;
                     1 or 'y' or 'w' for wavelength dispersion
             npoly : order of Legendre poly to fit to sigma vs. wavelength
@@ -216,7 +216,7 @@ def _xypix(self, ispec, wavelength):
         """
         raise NotImplementedError
         
-    def xypix(self, ispec, wavelength, xmin=0, xmax=None, ymin=0, ymax=None):
+    def xypix(self, ispec, wavelength, xmin=0, xmax=None, ymin=0, ymax=None, iwave_cache=None):
         """
         Evaluate PSF for spectrum[ispec] at given wavelength
         
@@ -242,11 +242,11 @@ def xypix(self, ispec, wavelength, xmin=0, xmax=None, ymin=0, ymax=None):
             if key in self._cache:
                 xx, yy, ccdpix = self._cache[key]
             else:
-                xx, yy, ccdpix = self._xypix(ispec, wavelength)
+                xx, yy, ccdpix = self._xypix(ispec, wavelength, iwave_cache)
                 self._cache[key] = (xx, yy, ccdpix)
         except AttributeError:
             self._cache = CacheDict(2500)
-            xx, yy, ccdpix = self._xypix(ispec, wavelength)
+            xx, yy, ccdpix = self._xypix(ispec, wavelength, iwave_cache)
             
         xlo, xhi = xx.start, xx.stop
         ylo, yhi = yy.start, yy.stop
@@ -612,14 +612,17 @@ def wmax_all(self):
         """Maximum wavelength seen by all spectra"""
         return self._wmax_all
     
-    def projection_matrix(self, spec_range, wavelengths, xyrange):
+    def projection_matrix(self, spec_range, wavelengths, xyrange, iwave_cache=None):
         """
         Returns sparse projection matrix from flux to pixels
     
-        Inputs:
+        Args:
             spec_range = (ispecmin, ispecmax) or scalar ispec
             wavelengths = array_like wavelengths
             xyrange  = (xmin, xmax, ymin, ymax)
+        Options:
+            iwave_cache: index of wavelengths[0] in the possibly larger
+                wavelengths array previously passed to self.cache_xypix()
             
         Usage:
             xyrange = xmin, xmax, ymin, ymax
@@ -645,14 +648,23 @@ def projection_matrix(self, spec_range, wavelengths, xyrange):
         A = np.zeros( (ny*nx, nspec*nflux) )
         tmp = np.zeros((ny, nx))
         for ispec in range(specmin, specmax):
-            for iflux, w in enumerate(wavelengths):
+            for iw, w in enumerate(wavelengths):
+                #- Are use using a pre-cached wavelength?
+                if iwave_cache is not None:
+                    iwave = iwave_cache + iw
+                else:
+                    iwave = None
+
                 #- Get subimage and index slices
-                xslice, yslice, pix = self.xypix(ispec, w, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
+                xslice, yslice, pix = self.xypix(ispec, w,
+                    xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
+                    iwave_cache = iwave
+                    )
                 
                 #- If there is overlap with pix_range, put into sub-region of A
                 if pix.shape[0]>0 and pix.shape[1]>0:
                     tmp[yslice, xslice] = pix
-                    ij = (ispec-specmin)*nflux + iflux
+                    ij = (ispec-specmin)*nflux + iw
                     A[:, ij] = tmp.ravel()
                     tmp[yslice, xslice] = 0.0