Changeset 43

Timestamp:

Jun 10, 2014 10:43:56 AM (10 years ago)

Author:

JohnLightfoot

Message:

calculates beams only over truncated fts_wn

File:

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

trunk/beamsgenerator.py

@@ -7 +7 @@
 import common.commonobjects as co
 
-def calculate_dirty_beam(npix, rpix, mult, bmax, bxbylist, wn_max, wn):
+def calculate_dirty_beam(npix, rpix, mult, bmax, bxbylist, wnmax, wn):
     lamb = 1.0 / (wn * 100.0)
-    lambda_min = 1.0 / (wn_max * 100.0)
+    lambda_min = 1.0 / (wnmax * 100.0)
     maxbx = bmax / lambda_min
 
@@ -39 +39 @@
 
 def calculate_primary_beam(npix, rpix, mult, mx, my, bmax, m1_diameter,
-  wn_max, wn):
+  wnmax, wn):
     # assume uniform illumination of primary. uv goes out to
     # +- bmax / lambda_min, fill uv plane appropriate to primary
     # size
     mirror = numpy.ones([mult*npix,mult*npix])
-    lambda_min = 1.0 / (wn_max * 100.0)
+    lambda_min = 1.0 / (wnmax * 100.0)
     lamb = 1.0 / (wn * 100.0)
     mirror[(numpy.sqrt(mx*mx + my*my) / (mult*rpix)) * 
@@ -99 +99 @@
         nsample = fts['ftsnsample']
         fts_wn = fts['fts_wn']
-        fts_lambda = fts['fts_lambda']
+        fts_wn_truncated = fts['fts_wn_truncated']
       
         # max pixel size (radians) that will fully sample the image.
-        # Sampling freq is 2 * Nyquist freq. So sampling freq = 2 * b_max / lambda_min.
+        # Sampling freq is 2 * Nyquist freq.
+        # So sampling freq = 2 * b_max / lambda_min.
         bmax = uvmapgen['bmax']
         self.result['max baseline [m]'] = bmax
-        max_pixsize = np.min(fts_lambda) / (2.0 * bmax)
+
+        lambda_min = 1.0 / (np.max(fts_wn) * 100.0)
+        max_pixsize = lambda_min / (2.0 * bmax)
         self.result['max pixsize [rad]'] = max_pixsize
 
@@ -111 +114 @@
         # radius of first null of Airy disk is theta=1.22lambda/d. Number of
         # pixels is beam diameter/max_pixsize. Calculate this for the 
-        # longest wavelength, largest beam
-        max_beam_radius = 1.22 * np.max(fts_lambda) /\
+        # longest wavelength that has flux (at wnmin), largest beam
+        
+        max_beam_radius = 1.22 * (1.0 / (np.min(fts_wn_truncated) * 100.0)) /\
           telescope['Primary mirror diameter']
         self.result['beam diam [rad]'] = 2.0 * max_beam_radius
@@ -128 +132 @@
         axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec') 
         axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec') 
-        axis3 = co.Axis(data=fts_wn, title='Frequency', units='cm-1')
-
-        # calculate beams for each point on fts_wn
+        axis3 = co.Axis(data=fts_wn_truncated, title='Frequency', units='cm-1')
+
+        # calculate beams for each point on fts_wn_truncated
         # oversample uv grid by mult to minimise aliasing
         mult = 5
@@ -139 +143 @@
         self.result['dirty beam'] = collections.OrderedDict()
 
-        wn_max = np.max(fts_wn)
+        wnmax = np.max(fts_wn_truncated)
         m1_diameter = telescope['Primary mirror diameter']
 
         # first calculate primary beam for each wn
         jobs = {}
-        for wn in fts_wn:
+        for wn in fts_wn_truncated:
             # submit jobs
-            indata = (npix, rpix, mult, mx, my, bmax, m1_diameter, wn_max,
+            indata = (npix, rpix, mult, mx, my, bmax, m1_diameter, wnmax,
               wn,)
             jobs[wn] = self.job_server.submit(calculate_primary_beam,
               indata, (), ('numpy',))
 
-        primary_amplitude_beam = np.zeros([npix,npix,len(fts_wn)], np.complex)
-        for iwn,wn in enumerate(fts_wn):
+        primary_amplitude_beam = np.zeros([npix,npix,len(fts_wn_truncated)],
+          np.complex)
+        for iwn,wn in enumerate(fts_wn_truncated):
             # collect and store results
             primary_intensity_beam, primary_amplitude_beam[:,:,iwn] = jobs[wn]()
@@ -164 +169 @@
 
         # second calculate dirty beam, same use of pp
-        for wn in fts_wn:
+        for wn in fts_wn_truncated:
             #
             # 1. uv plane covered by discrete points so should use a
@@ -176 +181 @@
             # positions. Fastest but adequate?
 
-            indata = (npix, rpix, mult, bmax, uvmapgen['bxby'], wn_max, wn,)
+            indata = (npix, rpix, mult, bmax, uvmapgen['bxby'], wnmax, wn,)
             jobs[wn] = self.job_server.submit(calculate_dirty_beam,
               indata, (), ('numpy',))
@@ -183 +188 @@
         # axes for dirty beam
         axis = np.arange(-rpix, rpix, dtype=np.float)
-        axis *= (np.min(fts_lambda) / bmax)
+        axis *= (lambda_min / bmax)
         axis *= 206264.806247
         axis1 = co.Axis(data=-axis, title='RA offset', units='arcsec') 
         axis2 = co.Axis(data=axis, title='Dec offset', units='arcsec') 
 
-        for wn in fts_wn:
+        for wn in fts_wn_truncated:
             dirty_beam = jobs[wn]()
             image = co.Image(data=dirty_beam, axes=[axis1, axis2],