at least runs to completion, probably doesn't get the right answer yet

py/specter/psf/psf.py

@@ -77,15 +77,15 @@ 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):
         """
         Fit the cross-sectional Gaussian sigma of PSF spots vs. wavelength.
         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