fix bug sls;add -2,2,-1 orders

ObservationSim/MockObject/Galaxy.py

@@ -3,7 +3,7 @@ import galsim
 import os, sys
 import astropy.constants as cons
 from astropy.table import Table
-from ._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
+from ._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders
 from .SpecDisperser import SpecDisperser
 from .MockObject import MockObject
 from scipy import interpolate
@@ -279,7 +279,7 @@ class Galaxy(MockObject):
             folding_threshold = 5.e-3
         gsp = galsim.GSParams(folding_threshold=folding_threshold)
         # nphotons_sum = 0
-
+        xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
 
@@ -314,8 +314,9 @@ class Galaxy(MockObject):
             spec_orders = sdp.compute_spec_orders()
             for k, v in spec_orders.items():
                 img_s = v[0]
-                origin_order_x = v[1]
-                origin_order_y = v[2]
+                img_s,orig_off = convolveGaussXorders(img_s,xOrderSigPlus[k])
+                origin_order_x = v[1]-orig_off
+                origin_order_y = v[2]-orig_off
                 specImg = galsim.ImageF(img_s)
                 photons = galsim.PhotonArray.makeFromImage(specImg)
                 photons.x += origin_order_x

ObservationSim/MockObject/MockObject.py

@@ -2,7 +2,7 @@ import galsim
 import numpy as np
 import astropy.constants as cons
 from astropy.table import Table
-from ._util import magToFlux, vc_A
+from ._util import magToFlux, vc_A, convolveGaussXorders
 from ._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG
 from .SpecDisperser import SpecDisperser
 
@@ -190,6 +190,8 @@ class MockObject(object):
         normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
                           names=('WAVELENGTH', 'FLUX'))
 
+
+        xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
         for i in range(len(bandpass_list)):
             bandpass = bandpass_list[i]
             psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)
@@ -211,8 +213,9 @@ class MockObject(object):
             spec_orders = sdp.compute_spec_orders()
             for k, v in spec_orders.items():
                 img_s = v[0]
-                origin_order_x = v[1]
-                origin_order_y = v[2]
+                img_s,orig_off = convolveGaussXorders(img_s,xOrderSigPlus[k])
+                origin_order_x = v[1]-orig_off
+                origin_order_y = v[2]-orig_off
                 specImg = galsim.ImageF(img_s)
                 photons = galsim.PhotonArray.makeFromImage(specImg)
                 photons.x += origin_order_x

ObservationSim/MockObject/SpecDisperser/SpecDisperser.py

@@ -101,6 +101,7 @@ class SpecDisperser(object):
         beam_names = self.grating_conf.beams
 
         for beam in beam_names:
+        #for beam in ['A','B']:
             all_orders[beam] = self.compute_spec(beam)
 
         return all_orders
@@ -123,7 +124,7 @@ class SpecDisperser(object):
         lam_index = argsort(lam_beam)
         conf_sens = self.grating_conf.sens[beam]
 
-        lam_intep = np.linspace(self.band_start, self.band_end, int((self.band_end - self.band_start) / 0.5))
+        lam_intep = np.linspace(self.band_start, self.band_end, int((self.band_end - self.band_start) / 0.1))
         
         thri = interpolate.interp1d(conf_sens['WAVELENGTH'], conf_sens['SENSITIVITY'])
         spci = interpolate.interp1d(self.spec['WAVELENGTH'], self.spec['FLUX'])
@@ -196,7 +197,7 @@ class SpecDisperser(object):
         return model, originOut_x, originOut_y
 
     def writerSensitivityFile(self, conffile = '', beam = '', w = None, sens = None):
-        orders={'A':'1st','B':'0st','C':'2st'}
+        orders={'A':'1st','B':'0st','C':'2st','D':'-1st','E':'-2st'}        
         sens_file_name = conffile[0:-5]+'_sensitivity_'+ orders[beam] + '.fits'
         if not os.path.exists(sens_file_name) == True:
             senstivity_out = Table(array([w,sens]).T, names=('WAVELENGTH', 'SENSITIVITY'))

ObservationSim/MockObject/SpecDisperser/disperse_c/interp.pyx

@@ -189,7 +189,7 @@ def interp_conserve_c(np.ndarray[DTYPE_t, ndim=1] x, np.ndarray[DTYPE_t, ndim=1]
             if (templmid[k+1] == tlam[i]) & (i < ntlam):
                 h = tlam[i]-tlam[i-1];
                 numsum+=h*(tf[i]+tf[i-1]);
-            else:
+            elif (i < ntlam):
                 i-=1;
                 h = templmid[k+1]-tlam[i];
                 numsum+=h*(tempfmid[k+1]+tf[i]);

ObservationSim/MockObject/_util.py

@@ -546,3 +546,21 @@ def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0):
     wave = np.array(h5file["wave"][model_tag_str][()]).ravel()
     flux = np.array(h5file["sed"][path][()]).ravel()
     return path, wave, flux
+
+from astropy.modeling.models import Gaussian2D
+from scipy import signal
+def convolveGaussXorders(img=None, sigma = 1):
+
+    offset = int(np.ceil(sigma*10))
+    g_size = 2*offset+1
+
+    m_cen = int(g_size/2)
+
+    g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma)
+    yp, xp = np.mgrid[0:g_size, 0:g_size]
+    g_PSF = g_PSF_(xp, yp)
+    psf = g_PSF/g_PSF.sum()
+    convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
+    return convImg, offset
+
+