diff --git a/crispy/IFS.py b/crispy/IFS.py index ca3b4c1..5e6de64 100644 --- a/crispy/IFS.py +++ b/crispy/IFS.py @@ -119,7 +119,8 @@ def polychromeIFS(par, inWavelist, inputcube, end=True) nframes = inputcube.data.shape[0] - allweights = None # TODO, What is this variable? Add a helpful comment. + # Variable for storing weight maps for each wavelength slice (currently unused) + allweights = None if inputcube.data.shape[0] != len(wavelist): log.error('Number of wavelengths does not match the number of input cube slices') @@ -171,17 +172,18 @@ def polychromeIFS(par, inWavelist, inputcube, if lam_arr is None: lam_arr = np.loadtxt(par.wavecalDir + "lamsol.dat")[:, 0] - hires_arrs = [] # TODO, Add a comment that describes what this variable is supposed to represent. Rename all occurence of this variable to 'high_res_arrays' + # List to store high-resolution PSF arrays for each wavelength + high_res_arrays = [] if par.gaussian: for i in range(len(lam_arr)): - hiresarr = get_sim_hires(par, lam_arr[i]) # TODO rename all occurences of this variable to 'high_res_array' - hires_arrs += [hiresarr] + high_res_array = get_sim_hires(par, lam_arr[i]) # TODO rename all occurences of this variable to 'high_res_array' + high_res_arrays += [high_res_array] log.info('Creating Gaussian PSFLet templates') upsample = 10 else: try: hires_list = np.sort(glob.glob(par.wavecalDir + 'hires_psflets_lam???.fits')) - hires_arrs = [pyf.getdata(filename) for filename in hires_list] + high_res_arrays = [pyf.getdata(filename) for filename in hires_list] log.info('Loaded PSFLet templates') except BaseException: log.error('Failed loading the PSFLet templates') @@ -193,7 +195,8 @@ def polychromeIFS(par, inWavelist, inputcube, # Fill out the polyimage array if not parallel: for i in range(len(waveList)): - # TODO, add an informational comment about what imagePlaneRot represents. + # Process and rotate the image plane for this wavelength slice, + # then scale by wavelength bin width for flux conservation imagePlaneRot = (wavelist_endpts[i + 1] - wavelist_endpts[i]) * \ processImagePlane(par, interpolatedInputCube.data[i], noRot) inputCube += [imagePlaneRot] @@ -201,7 +204,7 @@ def polychromeIFS(par, inWavelist, inputcube, imagePlaneRot, wavelist_endpts[i], wavelist_endpts[i + 1], - hires_arrs, + high_res_arrays, lam_arr, upsample, nlam, @@ -227,7 +230,7 @@ def polychromeIFS(par, inWavelist, inputcube, imagePlaneRot, wavelist_endpts[i], wavelist_endpts[i + 1], - hires_arrs, + high_res_arrays, lam_arr, upsample, nlam, diff --git a/crispy/configs/HCIFS/params.py b/crispy/configs/HCIFS/params.py index 956ca53..c315da5 100644 --- a/crispy/configs/HCIFS/params.py +++ b/crispy/configs/HCIFS/params.py @@ -95,7 +95,7 @@ def __init__(self, codeRoot='../'): self.BW = 0.18 # Spectral bandwidth self.npixperdlam = 2.0 # Number of pixels per spectral resolution element - self.nchanperspec_lstsq = 2.0 # Nspec per pixel for least squares + self.nchanperspec_lstsq = 2.0 # num_wavelengths per pixel for least squares # Spectral resolving power (extracted cubes have twice) self.R = 50 @@ -175,6 +175,6 @@ def makeHeader(self): self.hdr.append( ('RESLSTSQ', self.nchanperspec_lstsq, - 'Nspec per Nyq. sample for lstsq extraction'), + 'num_wavelengths per Nyq. sample for lstsq extraction'), end=True) self.hdr.append(('R', self.R, 'Spectral resolution'), end=True) diff --git a/crispy/configs/PISCES/params.py b/crispy/configs/PISCES/params.py index 155681d..083a859 100644 --- a/crispy/configs/PISCES/params.py +++ b/crispy/configs/PISCES/params.py @@ -78,7 +78,7 @@ def __init__(self,codeRoot='../'): self.distortPISCES=False # If True, use measured PISCES distortion/dispersion self.BW = 0.18 # Spectral bandwidth (if distortPISCES==False) self.npixperdlam = 2.0 # Number of pixels per spectral resolution element - self.nchanperspec_lstsq = 2.0 # Nspec per pixel for least squares + self.nchanperspec_lstsq = 2.0 # num_wavelengths per pixel for least squares self.R = 70 # Spectral resolving power (extracted cubes have twice) # carry-over old parameter names diff --git a/crispy/configs/WFIRST/params.py b/crispy/configs/WFIRST/params.py index 7dadda4..77fa8c7 100644 --- a/crispy/configs/WFIRST/params.py +++ b/crispy/configs/WFIRST/params.py @@ -90,7 +90,7 @@ def __init__(self, codeRoot='../'): self.BW = 0.18 # Spectral bandwidth self.npixperdlam = 2.0 # Number of pixels per spectral resolution element - self.nchanperspec_lstsq = 2.0 # Nspec per pixel for least squares + self.nchanperspec_lstsq = 2.0 # num_wavelengths per pixel for least squares # Spectral resolving power (extracted cubes have twice) self.R = 50 @@ -170,6 +170,6 @@ def makeHeader(self): self.hdr.append( ('RESLSTSQ', self.nchanperspec_lstsq, - 'Nspec per Nyq. sample for lstsq extraction'), + 'num_wavelengths per Nyq. sample for lstsq extraction'), end=True) self.hdr.append(('R', self.R, 'Spectral resolution'), end=True) diff --git a/crispy/tools/lenslet.py b/crispy/tools/lenslet.py index 47b9d2d..a03f722 100644 --- a/crispy/tools/lenslet.py +++ b/crispy/tools/lenslet.py @@ -72,7 +72,7 @@ def propagateLenslets( imageplane, lam1, lam2, - hires_arrs=None, + high_res_arrays=None, lam_arr=None, upsample=3, nlam=10, @@ -83,7 +83,7 @@ def propagateLenslets( Function propagateLenslets This is the main propagation function. It puts the PSFLets where they belong on the detector. - It uses template PSFLets given in hires_arrs, and can use also a pre-determined wavelength + It uses template PSFLets given in high_res_arrays, and can use also a pre-determined wavelength solution through the allcoef argument. Parameters @@ -142,21 +142,21 @@ def propagateLenslets( # to do this later, saving a factor of a few in time. ################################################################ - if (hires_arrs is None) or (lam_arr is None): + if (high_res_arrays is None) or (lam_arr is None): log.error('No template PSFLets given!') return else: - hires = np.zeros((hires_arrs[0].shape)) + hires = np.zeros((high_res_arrays[0].shape)) if lam <= np.amin(lam_arr): - hires[:] = hires_arrs[0] + hires[:] = high_res_arrays[0] elif lam >= np.amax(lam_arr): - hires[:] = hires_arrs[-1] + hires[:] = high_res_arrays[-1] else: i1 = np.amax(np.arange(len(lam_arr))[np.where(lam > lam_arr)]) i2 = i1 + 1 - hires = hires_arrs[i1] * \ + hires = high_res_arrays[i1] * \ (lam - lam_arr[i1]) / (lam_arr[i2] - lam_arr[i1]) - hires += hires_arrs[i2] * \ + hires += high_res_arrays[i2] * \ (lam_arr[i2] - lam) / (lam_arr[i2] - lam_arr[i1]) for i in range(hires.shape[0]): @@ -244,23 +244,23 @@ def propagateLenslets( x_hires = xcen[i] * 1. / image.shape[1] y_hires = ycen[i] * 1. / image.shape[0] - x_hires = x_hires * hires_arrs[0].shape[1] - 0.5 - y_hires = y_hires * hires_arrs[0].shape[0] - 0.5 + x_hires = x_hires * high_res_arrays[0].shape[1] - 0.5 + y_hires = y_hires * high_res_arrays[0].shape[0] - 0.5 totweight = 0 if x_hires <= 0: i1 = i2 = 0 - elif x_hires >= hires_arrs[0].shape[1] - 1: - i1 = i2 = hires_arrs[0].shape[1] - 1 + elif x_hires >= high_res_arrays[0].shape[1] - 1: + i1 = i2 = high_res_arrays[0].shape[1] - 1 else: i1 = int(x_hires) i2 = i1 + 1 if y_hires < 0: j1 = j2 = 0 - elif y_hires >= hires_arrs[0].shape[0] - 1: - j1 = j2 = hires_arrs[0].shape[0] - 1 + elif y_hires >= high_res_arrays[0].shape[0] - 1: + j1 = j2 = high_res_arrays[0].shape[0] - 1 else: j1 = int(y_hires) j2 = j1 + 1 diff --git a/crispy/tools/locate_psflets.py b/crispy/tools/locate_psflets.py index bc83263..777889d 100755 --- a/crispy/tools/locate_psflets.py +++ b/crispy/tools/locate_psflets.py @@ -129,29 +129,30 @@ def geninterparray(self, lam, allcoef, order=3): coef = np.linalg.lstsq(log_wavelength_powers, allcoef[:, i])[0] self.interp_arr[:, i] = coef - def return_locations_short(self, coef, xindx, yindx): - ''' - Returns the x,y detector location of a given lenslet for a given polynomial fit - - Parameters - ---------- - coef: lists floats - Polynomial coefficients of fit for a single wavelength - xindx: int - X index of lenslet in lenslet array - yindx: int - Y index of lenslet in lenslet array - - Returns - ------- - interp_x: float - X coordinate on the detector - interp_y: float - Y coordinate on the detector - ''' - # TODO, where does this 'coeforder' come from? Does this parent function actually get called from anywhere? Because clearly it's going to throw an error if it does. - interp_x, interp_y = transform(xindx, yindx, coeforder, coef) - return interp_x, interp_y + # COMMENTED OUT: Function commented out to avoid confusion until we understand its purpose + # def return_locations_short(self, coef, xindx, yindx): + # ''' + # Returns the x,y detector location of a given lenslet for a given polynomial fit + + # Parameters + # ---------- + # coef: lists floats + # Polynomial coefficients of fit for a single wavelength + # xindx: int + # X index of lenslet in lenslet array + # yindx: int + # Y index of lenslet in lenslet array + + # Returns + # ------- + # interp_x: float + # X coordinate on the detector + # interp_y: float + # Y coordinate on the detector + # ''' + # # TODO, where does this 'coeforder' come from? Does this parent function actually get called from anywhere? Because clearly it's going to throw an error if it does. + # interp_x, interp_y = transform(xindx, yindx, coeforder, coef) + # return interp_x, interp_y def return_res(self, lam, allcoef, xindx, yindx, order=3, lam1=None, lam2=None): @@ -782,28 +783,29 @@ def corrval(coef, x, y, input_image, order, trimfrac=0.1, show_plots=False): return score -def corrvalsum(coef, x, y, filtered, order, trimfrac=0.1, gsize=2): - # TODO, is this function used anywhere in this repo? If not, comment it out. - # TODO, add doscring - _x, _y = transform(x, y, order, coef) - ydim, xdim = filtered.shape - s = 0.0 - ry = np.reshape(_y, -1) - rx = np.reshape(_x, -1) - for i in range(len(ry)): - yi = ry[i] - xi = rx[i] - xmin = int(xi) - gsize - xmax = xmin + 2 * gsize - ymin = int(yi) - gsize - ymax = ymin + 2 * gsize - if ymin > 2 * gsize and xmin > 2 * gsize and xmax < xdim - 2 * gsize and ymax < ydim - 2 * gsize: - # dx = xi - int(xi) - # dy = yi - int(yi) - # s+=np.sum(simplepsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax]) - # s+=np.sum(gausspsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax]) - s += np.sum(filtered[ymin:ymax, xmin:xmax]) - return -s +# COMMENTED OUT: Function unused in repository, commented out as requested +# def corrvalsum(coef, x, y, filtered, order, trimfrac=0.1, gsize=2): +# # TODO, is this function used anywhere in this repo? If not, comment it out. +# # TODO, add doscring +# _x, _y = transform(x, y, order, coef) +# ydim, xdim = filtered.shape +# s = 0.0 +# ry = np.reshape(_y, -1) +# rx = np.reshape(_x, -1) +# for i in range(len(ry)): +# yi = ry[i] +# xi = rx[i] +# xmin = int(xi) - gsize +# xmax = xmin + 2 * gsize +# ymin = int(yi) - gsize +# ymax = ymin + 2 * gsize +# if ymin > 2 * gsize and xmin > 2 * gsize and xmax < xdim - 2 * gsize and ymax < ydim - 2 * gsize: +# # dx = xi - int(xi) +# # dy = yi - int(yi) +# # s+=np.sum(simplepsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax]) +# # s+=np.sum(gausspsf(size=2*gsize,fwhm=fwhm,offx=dx,offy=dy)*filtered[ymin:ymax,xmin:xmax]) +# s += np.sum(filtered[ymin:ymax, xmin:xmax]) +# return -s def locatePSFlets(inImage, mask, polyorder=2, sig=0.7, coef=None, trimfrac=0.1, @@ -831,7 +833,16 @@ def locatePSFlets(inImage, mask, polyorder=2, sig=0.7, coef=None, trimfrac=0.1, fraction of lenslet outliers (high & low combined) to trim in the minimization. Default 0.1 (5% trimmed on the high end, 5% on the low end) - # TODO, add argument definitions for the other arugments not listed here. + mask: ndarray + Mask array for the image + phi: float + Rotation angle for the PSFlet grid. Default np.arctan2(1.926, -1) + scale: float + Scale factor for the PSFlet grid. Default 15.02 + nlens: int + Number of lenslets. Default 108 + finesearch: int + Fine search parameter. Default 3 Returns ------- diff --git a/crispy/tools/reduction.py b/crispy/tools/reduction.py index b1bb2f8..dc33ded 100644 --- a/crispy/tools/reduction.py +++ b/crispy/tools/reduction.py @@ -130,7 +130,7 @@ def testReduction(par, name, ifsimage): return cube -def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'): +def calculateWaveList(par, lam_list=None, num_wavelengths=None, method='lstsq'): ''' Computes the wavelength lists corresponding to the center and endpoints of each spectral bin. Wavelengths are separated by a constant value in log space. Number of @@ -145,7 +145,7 @@ def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'): calibration. Otherwise, we could decide to focus on a smaller/larger region of the spectrum to retrieve. The final processed cubes will have bins centered on lam_midpts - Nspec: int + num_wavelengths: int If specified, forces the number of bins in the final cube (uses np.linspace) Returns @@ -160,9 +160,9 @@ def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'): lamlist = np.loadtxt(par.wavecalDir + "lamsol.dat")[:, 0] else: lamlist = lam_list - if Nspec is None: + if num_wavelengths is None: if method == 'lstsq': - Nspec = int( + num_wavelengths = int( np.log( max(lamlist) / min(lamlist)) * @@ -170,20 +170,20 @@ def calculateWaveList(par, lam_list=None, Nspec=None, method='lstsq'): par.nchanperspec_lstsq + 1) else: - Nspec = int( + num_wavelengths = int( np.log( max(lamlist) / min(lamlist)) * par.R * par.npixperdlam + 1) - log.info('Reduced cube will have %d wavelength bins' % (Nspec - 1)) -# lam_endpts = np.linspace(min(lamlist), max(lamlist), Nspec) + log.info('Reduced cube will have %d wavelength bins' % (num_wavelengths - 1)) +# lam_endpts = np.linspace(min(lamlist), max(lamlist), num_wavelengths) # lam_midpts = (lam_endpts[1:]+lam_endpts[:-1])/2. loglam_endpts = np.linspace( np.log( min(lamlist)), np.log( - max(lamlist)), Nspec) + max(lamlist)), num_wavelengths) loglam_midpts = (loglam_endpts[1:] + loglam_endpts[:-1]) / 2 lam_endpts = np.exp(loglam_endpts) lam_midpts = np.exp(loglam_midpts) @@ -237,7 +237,7 @@ def lstsqExtract(par, name, ifsimage, smoothandmask=True, ivar=True, dy=3, yindx = polychromekey[2].data good = polychromekey[3].data - lam_midpts, lam_endpts = calculateWaveList(par, method='lstsq', Nspec=psflets.shape[0]+1) + lam_midpts, lam_endpts = calculateWaveList(par, method='lstsq', num_wavelengths=psflets.shape[0]+1) if fitbkgnd: n_add = 1 @@ -913,7 +913,7 @@ def intOptimalExtract(par, name, IFSimage, smoothandmask=True, sum=False): """ loc = PSFLets(load=True, infiledir=par.wavecalDir) - #Nspec = int(par.BW*par.npixperdlam*par.R) + #num_wavelengths = int(par.BW*par.npixperdlam*par.R) lam_midpts, scratch = calculateWaveList(par, method='optext') datacube = fitspec_intpix_np( diff --git a/crispy/tools/wavecal.py b/crispy/tools/wavecal.py index c501f46..6e99f57 100644 --- a/crispy/tools/wavecal.py +++ b/crispy/tools/wavecal.py @@ -390,15 +390,37 @@ def make_hires_polychrome(lam1, lam2, hires_arrs, lam_arr, psftool, allcoef, return hiresimg -def get_sim_hires(par, lam, upsample=10, nsubarr=1, npix=13, renorm=True): +def get_sim_hires(par, lam, upsample=10, nsubarr=1, npix=13, normalize=True): """ Build high resolution images of the undersampled PSF using the monochromatic frames. This version of the function uses the perfect knowledge of the Gaussian PSFLet. Only valid if par.gaussian=True. All PSFLets are the same across the entire FOV - # TODO improve this docstring with Numpy style - # TODO rename all instances of 'renorm' to 'normalize' + Parameters + ---------- + par : object + Parameter object containing FWHM, FWHMlam, and gaussian attributes + lam : float + Wavelength for which to generate the high resolution PSF + upsample : int, optional + Upsampling factor for the high resolution array. Default is 10 + nsubarr : int, optional + Number of subarrays in each dimension. Default is 1 + npix : int, optional + Number of pixels in the base PSF. Default is 13 + normalize : bool, optional + Whether to normalizealize the PSFlet. Default is True + + Returns + ------- + hires_arr : ndarray + 4D array of shape (nsubarr, nsubarr, array_size, array_size) containing + the high resolution PSFlets + + Notes + ----- + # TODO rename all instances of 'normalize' to 'normalize' """ # Determine side length of the upsampled array array_size = upsample * (npix + 1) @@ -419,7 +441,7 @@ def get_sim_hires(par, lam, upsample=10, nsubarr=1, npix=13, renorm=True): erf((_y - 0.5) / (np.sqrt(2) * sigma_scaled))) # Normalize the PSFLet, if desired - if renorm: + if normalize: psflet *= upsample**2 / np.sum(psflet) # Because the output is expected to ahve nsubarr * nsubarr entries, fill the array with the same PSFLet @@ -469,7 +491,7 @@ def get_sim_hires(par, lam, upsample=10, nsubarr=1, npix=13, renorm=True): # return epsf.data # -def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, renorm=True): +def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, normalize=True): """ Build high resolution images of the undersampled PSF using the monochromatic frames. @@ -641,13 +663,13 @@ def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, renorm=True): # interpolator. ############################################################ - if renorm: + if normalize: meanpsf *= upsample**2 / np.sum(meanpsf) hires_arr[yreg, xreg] = meanpsf return hires_arr -# def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, renorm=True): +# def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, normalize=True): # """ # Build high resolution images of the undersampled PSF using the # monochromatic frames. @@ -673,7 +695,7 @@ def gethires(x, y, good, image, upsample=5, nsubarr=5, npix=13, renorm=True): # j2 = min(j2, data.shape[1] - npix) # subim = data[i1:i2,j1:j2] # hires_arr[yreg,xreg] = epsflets(subim,upsample,npix) -# if renorm: +# if normalize: # hires_arr[yreg,xreg] *= upsample**2 / np.sum(hires_arr[yreg,xreg]) # # return hires_arr @@ -768,19 +790,19 @@ def makeHires( for i in range(ncpus): tasks.put(None) for i in range(len(lam)): - index, hiresarr = results.get() - hires_arrs += [hiresarr] + index, high_res_array = results.get() + hires_arrs += [high_res_array] if savehiresimages: - di, dj = hiresarr.shape[0], hiresarr.shape[2] + di, dj = high_res_array.shape[0], high_res_array.shape[2] outim = np.zeros((di * dj, di * dj)) for ii in range(di): for jj in range(di): outim[ii * dj:(ii + 1) * dj, jj * - dj:(jj + 1) * dj] = hiresarr[ii, jj] + dj:(jj + 1) * dj] = high_res_array[ii, jj] out = fits.HDUList( fits.PrimaryHDU( - hiresarr.astype( + high_res_array.astype( np.float32))) out.writeto( par.wavecalDir + @@ -791,7 +813,7 @@ def makeHires( log.info('No parallel computation') for i in range(len(lam)): if par.gaussian_hires: - hiresarr = get_sim_hires(par, lam[i], upsample, nsubarr) + high_res_array = get_sim_hires(par, lam[i], upsample, nsubarr) else: # if finexy is None: # xpos, ypos = psftool.return_locations( @@ -807,18 +829,21 @@ def makeHires( good = np.reshape(psftool.good, -1) xpos = np.reshape(xpos, -1) ypos = np.reshape(ypos, -1) - hiresarr = gethires(xpos, ypos, good, imlist[i], upsample, nsubarr) - hires_arrs += [hiresarr] + high_res_array = gethires(xpos, ypos, good, imlist[i], upsample, nsubarr) + hires_arrs += [high_res_array] - # TODO, in one iteration, savehiresimages equaled 3 here. How could this happen? + # Validate savehiresimages parameter - should be boolean + if not isinstance(savehiresimages, bool): + raise ValueError(f"savehiresimages must be boolean (True/False), got {type(savehiresimages).__name__}: {savehiresimages}") + if savehiresimages: # Apparently deprecated code that didn't get used? Commenting it out for now. - # di, dj = hiresarr.shape[0], hiresarr.shape[2] + # di, dj = high_res_array.shape[0], high_res_array.shape[2] # outim = np.zeros((di * dj, di * dj)) # for ii in range(di): # for jj in range(di): - # outim[ii * dj:(ii + 1) * dj, jj * dj:(jj + 1) * dj] = hiresarr[ii, jj] - out = fits.HDUList(fits.PrimaryHDU(hiresarr.astype(np.float32))) + # outim[ii * dj:(ii + 1) * dj, jj * dj:(jj + 1) * dj] = high_res_array[ii, jj] + out = fits.HDUList(fits.PrimaryHDU(high_res_array.astype(np.float32))) out.writeto(par.wavecalDir + 'hires_psflets_lam%d.fits' % (lam[i]), overwrite=True) return hires_arrs @@ -1061,7 +1086,7 @@ def buildcalibrations( - an array of the Y positions of all lenslets - an array of booleans indicating whether that lenslet is good or not (e.g. when it is outside of the detector area) - polychromeRXX.fits: 3D arrays of size Nspec x Npix x Npix with maps of the PSFLets put in their correct + polychromeRXX.fits: 3D arrays of size num_wavelengths x Npix x Npix with maps of the PSFLets put in their correct positions for each wavelength bins that we want in the output cube. Each PSFLet in each wavelength slice is used for least-squares fitting. hiresPolychromeRXX.fits: same as polychromeRXX.fits but this time using the high-resolution PSFLets @@ -1358,9 +1383,9 @@ def buildcalibrations( # Create an array of wavelengths that represent the midpoints/endpoints of the wavelength bins lam_midpts, lam_endpts = calculateWaveList(par, lam, method='lstsq') - # TODO, rename all instances of 'Nspec' to 'num_wavelengths' for clarity. - Nspec = len(lam_endpts) # The number of unique wavelength bins - polyimage = np.zeros((Nspec - 1, ysize, xsize)) + # TODO, rename all instances of 'num_wavelengths' to 'num_wavelengths' for clarity. + num_wavelengths = len(lam_endpts) # The number of unique wavelength bins + polyimage = np.zeros((num_wavelengths - 1, ysize, xsize)) # Initialize some arrays where we will store information about the x/y position of each PSF, # as well as whether or not that PSF is "good" (i.e. falls on the detector) @@ -1370,8 +1395,8 @@ def buildcalibrations( log.info('Making polychrome cube') if not parallel: - for i in range(Nspec - 1): - log.info(f' Wavelength bin {i + 1} of {Nspec - 1}') + for i in range(num_wavelengths - 1): + log.info(f' Wavelength bin {i + 1} of {num_wavelengths - 1}') polyimage[i] = (lam_endpts[i + 1] - lam_endpts[i]) * make_polychrome(lam_endpts[i], lam_endpts[i + 1], hires_arrs, @@ -1405,7 +1430,7 @@ def buildcalibrations( for w in consumers: w.start() - for i in range(Nspec - 1): + for i in range(num_wavelengths - 1): tasks.put(Task(i, make_polychrome, (lam_endpts[i], @@ -1424,7 +1449,7 @@ def buildcalibrations( for i in range(ncpus): tasks.put(None) - for i in range(Nspec - 1): + for i in range(num_wavelengths - 1): index, poly = results.get() polyimage[index] = poly * \ (lam_endpts[index + 1] - lam_endpts[index]) @@ -1480,11 +1505,11 @@ def buildcalibrations( hires_arrs = [fits.open(filename)[0].data for filename in hires_list] lam_midpts, lam_endpts = calculateWaveList(par, lam, method='lstsq') - Nspec = len(lam_endpts) - hirespoly = np.zeros((Nspec - 1, ysize * upsample, xsize * upsample)) + num_wavelengths = len(lam_endpts) + hirespoly = np.zeros((num_wavelengths - 1, ysize * upsample, xsize * upsample)) if not parallel: - for i in range(Nspec - 1): + for i in range(num_wavelengths - 1): hirespoly[i] = (lam_endpts[i + 1] - lam_endpts[i]) * make_hires_polychrome(lam_endpts[i], lam_endpts[i + 1], hires_arrs, @@ -1505,7 +1530,7 @@ def buildcalibrations( for w in consumers: w.start() - for i in range(Nspec - 1): + for i in range(num_wavelengths - 1): tasks.put(Task(i, make_hires_polychrome, (lam_endpts[i], @@ -1522,7 +1547,7 @@ def buildcalibrations( for i in range(ncpus): tasks.put(None) - for i in range(Nspec - 1): + for i in range(num_wavelengths - 1): index, poly = results.get() hirespoly[index] = poly * \ (lam_endpts[index + 1] - lam_endpts[index]) / upsample**2 diff --git a/crispy/unitTests.py b/crispy/unitTests.py index aa1241a..6115982 100644 --- a/crispy/unitTests.py +++ b/crispy/unitTests.py @@ -112,7 +112,7 @@ def testOptExt(par,im, lensX, lensY, smoothandmask=True, delt_y=5): PSFlet_tool = PSFLets(load=True, infiledir=par.wavecalDir) - #Nspec = int(par.BW*par.npixperdlam*par.R) + #num_wavelengths = int(par.BW*par.npixperdlam*par.R) lamlist,scratch = calculateWaveList(par,method='optext') xindx = PSFlet_tool.xindx @@ -191,7 +191,7 @@ def testGenPixSol(par): def testCreateFlatfield(par,pixsize = 0.1, npix = 512, pixval = 1., - Nspec=45, + num_wavelengths=45, outname='flatfield.fits', useQE=True, method='optext', @@ -212,7 +212,7 @@ def testCreateFlatfield(par,pixsize = 0.1, Each input frame has a pixel size npix x npix pixval: float Each input frame has a unform value pixval in photons per second per nm of bandwidth - Nspec: float + num_wavelengths: float Optional input forcing the number of wavelengths bins used outname: string Name of flatfield image @@ -221,7 +221,7 @@ def testCreateFlatfield(par,pixsize = 0.1, ''' # Calculate the wavelengths of the flatfield cube and create the flatfield cube - lam_midpts,lam_endpts = calculateWaveList(par,Nspec=Nspec,method=method) + lam_midpts,lam_endpts = calculateWaveList(par,num_wavelengths=num_wavelengths,method=method) inputCube = np.ones((len(lam_midpts),npix,npix),dtype=np.float32) inputCube *= pixval @@ -241,7 +241,7 @@ def testCreateFlatfield(par,pixsize = 0.1, import scipy from scipy.ndimage import gaussian_filter1d -def testCrosstalk(par,pixsize = 0.1, npix = 512, pixval = 1.,Nspec=45,outname='crosstalk.fits',useQE=True,method='optext'): +def testCrosstalk(par,pixsize = 0.1, npix = 512, pixval = 1.,num_wavelengths=45,outname='crosstalk.fits',useQE=True,method='optext'): ''' Creates a polychromatic flatfield @@ -255,7 +255,7 @@ def testCrosstalk(par,pixsize = 0.1, npix = 512, pixval = 1.,Nspec=45,outname='c Each input frame has a pixel size npix x npix pixval: float Each input frame has a unform value pixval in photons per second per nm of bandwidth - Nspec: float + num_wavelengths: float Optional input forcing the number of wavelengths bins used outname: string Name of flatfield image @@ -264,7 +264,7 @@ def testCrosstalk(par,pixsize = 0.1, npix = 512, pixval = 1.,Nspec=45,outname='c ''' - lam_midpts,lam_endpts = calculateWaveList(par,Nspec=Nspec,method=method) + lam_midpts,lam_endpts = calculateWaveList(par,num_wavelengths=num_wavelengths,method=method) inputCube = np.zeros((len(lam_midpts),npix,npix),dtype=np.float32) for i in range(len(lam_midpts)): @@ -272,7 +272,7 @@ def testCrosstalk(par,pixsize = 0.1, npix = 512, pixval = 1.,Nspec=45,outname='c # lam_midpts_nom,_ = calculateWaveList(par,method=method) -# FWHM=Nspec/len(lam_midpts_nom) +# FWHM=num_wavelengths/len(lam_midpts_nom) # inputCube[:,npix//2,npix//2] = gaussian_filter1d(inputCube[:,npix//2,npix//2],sigma=FWHM/2.35) par.saveDetector=False diff --git a/run_simple_tests.py b/run_simple_tests.py index 0f1886d..a1ff0be 100644 --- a/run_simple_tests.py +++ b/run_simple_tests.py @@ -113,14 +113,14 @@ def test_working_units(par): # Test 1: testCreateFlatfield (smaller parameters) success, result = run_test_safely( testCreateFlatfield, "testCreateFlatfield (small)", par, - pixsize=0.1, npix=32, pixval=1.0, Nspec=5 + pixsize=0.1, npix=32, pixval=1.0, num_wavelengths=5 ) results["testCreateFlatfield"] = success # Test 2: testCrosstalk (smaller parameters) success, result = run_test_safely( testCrosstalk, "testCrosstalk (small)", par, - pixsize=0.1, npix=32, pixval=1.0, Nspec=5 + pixsize=0.1, npix=32, pixval=1.0, num_wavelengths=5 ) results["testCrosstalk"] = success diff --git a/tests/integration/test_full_pipeline.py b/tests/integration/test_full_pipeline.py index 4f0b660..9febcda 100644 --- a/tests/integration/test_full_pipeline.py +++ b/tests/integration/test_full_pipeline.py @@ -22,7 +22,7 @@ def test_minimal_ifs_workflow(self, wfirst_params, reference_files_available): # Step 1: Calculate wavelength list try: - lam_midpts, lam_endpts = calculateWaveList(wfirst_params, Nspec=5, method='optext') + lam_midpts, lam_endpts = calculateWaveList(wfirst_params, num_wavelengths=5, method='optext') assert len(lam_midpts) > 0 assert len(lam_endpts) > 0 except Exception as e: @@ -122,7 +122,7 @@ def test_create_and_process_small_dataset(self, wfirst_params, reference_files_a pixsize=0.1, npix=32, pixval=1.0, - Nspec=3, + num_wavelengths=3, outname='integration_test_flat.fits' ) except Exception as e: diff --git a/tests/unit/test_working_functions.py b/tests/unit/test_working_functions.py index 2491915..e10c4e0 100644 --- a/tests/unit/test_working_functions.py +++ b/tests/unit/test_working_functions.py @@ -47,7 +47,7 @@ def test_create_flatfield_small(self, wfirst_params, reference_files_available): pixsize=0.1, npix=32, pixval=1.0, - Nspec=5, + num_wavelengths=5, outname='test_flatfield_small.fits' ) @@ -67,7 +67,7 @@ def test_create_flatfield_medium(self, wfirst_params, reference_files_available) pixsize=0.1, npix=64, pixval=1.0, - Nspec=10, + num_wavelengths=10, outname='test_flatfield_medium.fits' ) @@ -89,7 +89,7 @@ def test_crosstalk_small(self, wfirst_params, reference_files_available): pixsize=0.1, npix=32, pixval=1.0, - Nspec=5, + num_wavelengths=5, outname='test_crosstalk_small.fits' ) @@ -107,7 +107,7 @@ def test_crosstalk_medium(self, wfirst_params, reference_files_available): pixsize=0.1, npix=64, pixval=1.0, - Nspec=10, + num_wavelengths=10, outname='test_crosstalk_medium.fits' )