Changeset 45

Timestamp:

Sep 19, 2014 4:52:22 PM (10 years ago)

Author:

JohnLightfoot

Message:

almost working

File:

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

trunk/skygenerator.py

@@ -51 +51 @@
 
 
-class TestLine(object):
+class linespectrum(object):
     """Class to return spectrum with line at wn.
     """
-    def __init__(self, line_wn, frequency_axis):
-        self.line_wn = line_wn
+    def __init__(self, temperature, linefreq, frequency_axis):
+        self.temperature = temperature
+        self.linefreq = linefreq
         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
+        # spectrum is a sync function centred on the specified line freq
+        x = (self.frequency_axis - self.linefreq) / \
+          (self.frequency_axis[1] - self.frequency_axis[0])
+        spectrum = np.sinc(x) * self.temperature
+
         return spectrum
 
@@ -78 +82 @@
         fts = self.previous_results['fts']
         fts_wn_truncated = fts['fts_wn_truncated']
+        cutoffmin = fts['wnmin']
+        cutoffmax = fts['wnmax']
 
         beamsgenerator = self.previous_results['beamsgenerator']
@@ -88 +94 @@
 
         sky = self.parameters['substages']['Sky']
-        columns = sky['SourceNum'].keys()
+        columns = sky['sourcenum'].keys()
 
         self.result['sources'] = collections.OrderedDict()
 
         for column in columns:
-            temp = sky['SourceNum'][column]
+            temp = sky['sourcenum'][column]
             sourcenum = int(round(temp))
             if sourcenum not in self.result['sources'].keys():
                 self.result['sources'][sourcenum] = {}
 
-            type = sky['Type'][column]
+            type = sky['type'][column]
 
             temp = sky['x pos [asec]'][column]
@@ -106 +112 @@
             ypos = float(temp)
 
-            temp = sky['Temp'][column]
+            temp = sky['xwidth'][column]
+            try:
+                xwidth = float(temp)
+            except:
+                xwidth = None
+
+            temp = sky['ywidth'][column]
+            try:
+                ywidth = float(temp)
+            except:
+                ywidth = None
+
+            spectrum = sky['spectrum'][column]
+
+            temp = sky['temperature'][column]
             temperature = float(temp)
 
-            temp = sky['cutoffmin'][column]
-            cutoffmin = float(temp)
-
-            temp = sky['cutoffmax'][column]
-            cutoffmax = float(temp)
+            temp = sky['linefreq'][column]
+            try:
+                linefreq = float(temp)
+            except:
+                linefreq = None
 
             temp = sky['emissivity'][column]
             emissivity = float(temp)
 
-            print 'generating source:%s type:%s xpos:%s ypos:%s temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
-              sourcenum, type, xpos, ypos, temperature, cutoffmin, cutoffmax,
-              emissivity)
-
-            spectrum_func = BB_spectrum(temperature, fts_wn_truncated,
-              cutoffmin, cutoffmax, emissivity)
-#            spectrum_func = TestLine(50.0, fts_wn_truncated)
+            print 'generating source:%s type:%s xpos:%s ypos:%s' % (sourcenum, 
+              type, xpos, ypos)
+            if type.upper().strip() == 'GAUSSIAN':
+                print '  fwhmx:%s fwhmy:%s' % (xwidth, ywidth)
+
+            if spectrum.upper().strip() == 'BLACKBODY':
+                print '  blackbody spectrum temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
+                  temperature, cutoffmin, cutoffmax, emissivity)
+                spectrum_func = BB_spectrum(temperature, fts_wn_truncated,
+                  cutoffmin, cutoffmax, emissivity)
+
+            elif spectrum.upper().strip() == 'LINE':
+                print '  line spectrum brightness temperature:%s cutoffmin:%s cutoffmax:%s e:%s' % (
+                  temperature, cutoffmin, cutoffmax, emissivity)
+                spectrum_func = linespectrum(temperature, linefreq, fts_wn_truncated)
 
             if type.upper().strip() == 'POINT':
                 source_spectrum = self._create_point_source(xpos, ypos,
                   skymodel, spatial_axis, fts_wn_truncated, spectrum_func)
+            elif type.upper().strip() == 'GAUSSIAN':
+                source_spectrum = self._create_gaussian_source(xpos, ypos,
+                  xwidth, ywidth, skymodel, spatial_axis, fts_wn_truncated,
+                  spectrum_func)
             else:
                 source_spectrum = None
@@ -140 +172 @@
 
         return self.result
+
+    def _create_gaussian_source(self, xpos, ypos, xwidth, ywidth,
+      skymodel, spatial_axis, frequency_axis, spectrum_function):
+        """Create a point source.
+        """
+        # calculate spectrum
+        spectrum = spectrum_function.calculate()
+
+        # calculate Gaussian profile - a cunning use of array slicing found
+        # on the web
+        x = spatial_axis
+        y = x[:,np.newaxis]
+        profile = np.exp(-4*np.log(2) * ((x-xpos)**2/xwidth**2 + (y-ypos)**2/ywidth**2))
+
+#        # apodise
+#        profile[((x-xpos)**2 + (y-ypos)**2) < 4] = 1.0
+
+        # go through freq planes and add source to each
+        for iwn,wn in enumerate(frequency_axis):
+            # add to sky model
+            skymodel[:,:,iwn] += profile * spectrum[iwn]
+
+        # return spectrum
+        axis = co.Axis(data=frequency_axis, title='wavenumber', units='cm-1')
+        spectrum = co.Spectrum(data=spectrum, axis=axis,
+          title='Source spectrum', units='W sr-1 m-2 Hz-1')
+
+        return spectrum
 
     def _create_point_source(self, xpos, ypos, skymodel, spatial_axis,