Changeset 38

Timestamp:

Jun 10, 2014 10:40:25 AM (10 years ago)

Author:

JohnLightfoot

Message:

TestLine added

File:

• trunk/skygenerator.py (modified) (9 diffs)

trunk/skygenerator.py

@@ -14 +14 @@
     """
     freq = wavenumber * sc.c * 100.0
-    jnu = 2.0 * sc.h * pow(freq,3) / (pow(sc.c,2) *
-      (math.exp((sc.h * freq) / (sc.k * temperature)) - 1.0))
+    if freq > 0:
+        jnu = 2.0 * sc.h * pow(freq,3) / (pow(sc.c,2) *
+          (math.exp((sc.h * freq) / (sc.k * temperature)) - 1.0))
+    else:
+        jnu = 0.0
 
     return jnu
+
+class BB_spectrum(object):
+    """Class to generate BB spectrum.
+    """
+    def __init__(self, temperature, frequency_axis, cutoffmin, cutoffmax,
+      emissivity):
+        self.temperature = temperature
+        self.frequency_axis = frequency_axis
+        self.cutoffmin = cutoffmin
+        self.cutoffmax = cutoffmax
+        self.emissivity = emissivity
+
+    def calculate(self):
+        spectrum = np.zeros(np.shape(self.frequency_axis))
+        # ignore floating point errors
+        old_settings = np.seterr(all='ignore')
+
+        for iwn,wn in enumerate(self.frequency_axis):
+            # simulate real-life 'rounded' cutoffs numerically  
+            f1 = 1.0 / (1.0 + pow(self.cutoffmin/wn, 18) + 
+              pow(wn/self.cutoffmax, 24))
+            spectrum[iwn] = black_body(self.temperature, wn) * f1 * \
+              self.emissivity
+
+        # restore fp behaviour
+        ignore = np.seterr(**old_settings)
+
+        return spectrum
+
+
+class TestLine(object):
+    """Class to return spectrum with line at wn.
+    """
+    def __init__(self, line_wn, frequency_axis):
+        self.line_wn = line_wn
+        self.frequency_axis = frequency_axis
+
+    def calculate(self):
+        spectrum = np.zeros(np.shape(self.frequency_axis))
+        spectrum[self.frequency_axis==self.line_wn] = 1.0e-15
+        return spectrum
 
 
@@ -33 +77 @@
 
         fts = self.previous_results['fts']
-        frequency_axis = fts['fts_wn']
-        nspec = len(frequency_axis)
+        fts_wn_truncated = fts['fts_wn_truncated']
 
         beamsgenerator = self.previous_results['beamsgenerator']
@@ -42 +85 @@
         # skymodel is complex so that its fft can hold truncated version
         # of infinitesimally sampled map - does that make sense?
-        skymodel = np.zeros([npix, npix, nspec], np.complex)
+        skymodel = np.zeros([npix, npix, len(fts_wn_truncated)], np.complex)
 
         sky = self.parameters['substages']['Sky']
@@ -79 +122 @@
               emissivity)
 
+            spectrum_func = BB_spectrum(temperature, fts_wn_truncated,
+              cutoffmin, cutoffmax, emissivity)
+#            spectrum_func = TestLine(50.0, fts_wn_truncated)
+
             if type.upper().strip() == 'POINT':
                 source_spectrum = self._create_point_source(xpos, ypos,
-                  temperature, 
-                  cutoffmin, cutoffmax, emissivity, skymodel, 
-                  spatial_axis, frequency_axis)
+                  skymodel, spatial_axis, fts_wn_truncated, spectrum_func)
             else:
                 source_spectrum = None
@@ -92 +137 @@
         self.result['sky model'] = skymodel
         self.result['spatial axis'] = spatial_axis
-        self.result['frequency axis'] = frequency_axis 
+        self.result['frequency axis'] = fts_wn_truncated 
 
         return self.result
 
-    def _create_point_source(self, xpos, ypos, temperature, cutoffmin,
-      cutoffmax, emissivity, skymodel, spatial_axis, frequency_axis):
+    def _create_point_source(self, xpos, ypos, skymodel, spatial_axis,
+      frequency_axis, spectrum_function):
         """Create a point source.
         """
@@ -128 +173 @@
             shift[:,i] *= shifty
 
-        # calculate black-body spectrum, modified for cutoffs in 
-        # sensitivity and source emissivity
-        bbmod = np.zeros(np.shape(frequency_axis))
-        # ignore floating point errors
-        old_settings = np.seterr(all='ignore')
-        for iwn,wn in enumerate(frequency_axis):
-            # simulate real-life 'rounded' cutoffs numerically  
-            f1 = 1.0 / (1.0 + pow(cutoffmin/wn, 18) + pow(wn/cutoffmax, 24))
-            bbmod[iwn] = black_body(temperature, wn) * f1 * emissivity
-        # restore fp behaviour
-        ignore = np.seterr(**old_settings)
+        # calculate spectrum
+        spectrum = spectrum_function.calculate()
 
         # go through freq planes and add point source to each
@@ -144 +180 @@
             # create point in frequency space
             temp = np.zeros([nx,nx])
-            temp[0,0] = bbmod[iwn]
+            temp[0,0] = spectrum[iwn]
             # 2d fft
             temp = np.fft.fft2(temp)
@@ -151 +187 @@
             # transform back
             temp = np.fft.ifft2(temp)
-            temp = np.real(temp)
 
             # add to sky model
@@ -158 +193 @@
         # return spectrum
         axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
-        spectrum = co.Spectrum(data=bbmod, axis=axis, title='Source spectrum',
-          units='W sr-1 m-2 Hz-1')
+        spectrum = co.Spectrum(data=spectrum, axis=axis,
+          title='Source spectrum', units='W sr-1 m-2 Hz-1')
 
         return spectrum