Merge branch 'master' of github.com:mikeireland/pyxao

Almost have chara.py working again.

Author: mikeireland

chara.py

@@ -4,7 +4,12 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
-import pdb
+try:
+    import ipdb
+except:
+    #The following ling is dodgy, but still enables set_trace()
+    import pdb as ipdb
+    
 plt.ion()
 #from PyQt4 import QtGui
 #from PyQt4 import QtCore
@@ -32,7 +37,7 @@
 #dm  = pyxao.DeformableMirror(wavefronts=[wf_sense,wf_image],actuator_pitch=0.135,geometry='hexagonal', plotit=True,central_actuator=True)
 #wfs = pyxao.ShackHartmann(wavefronts=[wf_sense],lenslet_pitch = 0.17,plotit=True, central_lenslet=True)
 
-#pdb.set_trace()
+#ipdb.set_trace()
 #print("Click to continue")
 #dummy=plt.ginput(1)
 
@@ -48,16 +53,17 @@
 #See if the reconstructor works!
 #sensors,ims  = aos.correct_twice()
 
-im_mn, im_perfect = aos.run_loop(dt=0.001,plotit=True,niter=260,nframesbetweenplots=50)
+im_mn= aos.run_loop(dt=0.001,plotit=True,niter=260,nframesbetweenplots=50)
 
+#im_perfect = ??? #How do we compute Strehl now?
 
 #Uncomment the line below instead to see uncorrected seeing.
 #im_mn = aos.run_loop(plotit=True,gain=0.0)
 
 #Strehl: Regrid in order to sample finely the peak.
-strehl = np.max(ot.utils.regrid_fft(im_mn,(1024,1024)))/np.max(ot.utils.regrid_fft(im_perfect,(1024,1024)))
+#strehl = np.max(ot.utils.regrid_fft(im_mn,(1024,1024)))/np.max(ot.utils.regrid_fft(im_perfect,(1024,1024)))
 
-sz = im_mn.shape[0]
-plt.clf()
-plt.imshow(im_mn[sz//2-20:sz//2+20,sz//2-20:sz//2+20],interpolation='nearest', cmap=cm.gist_heat)
-print("Strehl: {0:6.3f}".format(strehl))
+#sz = im_mn.shape[0]
+#plt.clf()
+#plt.imshow(im_mn[sz//2-20:sz//2+20,sz//2-20:sz//2+20],interpolation='nearest', cmap=cm.gist_heat)
+#print("Strehl: {0:6.3f}".format(strehl))

minerva_red.py

@@ -4,7 +4,10 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
-import pdb
+try:
+    import ipdb
+except:
+    import pdb
 plt.ion()
 np.random.seed(1)
 
@@ -13,43 +16,53 @@
 #Strehl:
 #np.max(ot.utils.regrid_fft(im_mn_on,(1024,1024)))/np.max(ot.utils.regrid_fft(im_mn,(1024,1024)))
 
+"""
+	Observations:
+		- it's different for different sampling values because of how ot.kmf() computes the phase screen
+"""
+
 #Initialize some wavefronts. 0.56 to 0.76 microns sense (average 0.66 microns)
 #0.77 to 0.93 microns science (average 0.85 microns). For such a narrow field of
 #view, we can assume monochromatic 
 pupil={'type':"annulus", 'dout':0.7,'din':0.32}
 # Wavefront sensor wavelength
-wf_sense = pyxao.Wavefront(wave=0.66e-6,m_per_px=0.01,sz=128,pupil=pupil)
+scale_factor = 1
+wf_sense = pyxao.Wavefront(wave=0.66e-6,m_per_px=0.01/scale_factor,sz=256*scale_factor,pupil=pupil)
 # Science wavelength
-wf_image = pyxao.Wavefront(wave=0.85e-6,m_per_px=0.01,sz=128,pupil=pupil)
+wf_image = pyxao.Wavefront(wave=0.85e-6,m_per_px=0.01/scale_factor,sz=256*scale_factor,pupil=pupil)
 
 dm  = pyxao.DeformableMirror(wavefronts=[wf_sense,wf_image],actuator_pitch=0.115,geometry='square', plotit=False,central_actuator=True)
-wfs = pyxao.ShackHartmann(wavefronts=[wf_sense],lenslet_pitch = 0.105,central_lenslet=False,sampling=1,plotit=False,geometry='hexagonal')
+wfs = pyxao.ShackHartmann(wavefronts=[wf_sense],lenslet_pitch = 0.105,central_lenslet=False,sampling=1,plotit=False,geometry='hexagonal',N_phot = 1600 * 1e4 * 1/2000 * 0.90 * 1000)
 
 #Add an atmosphere model to our wavefronts. 
-atm = pyxao.Atmosphere(sz=1024, m_per_px=wf_sense.m_per_px,r_0=[0.1,0.1,0.1]) #For 1.7" seeing, try: ,r_0=[0.1,0.1,0.1])
+atm = pyxao.Atmosphere(sz=512*scale_factor, m_per_px=wf_sense.m_per_px,r_0=[0.1,0.1,0.1],v_wind=[0,0,0],seed=5) #For 1.7" seeing, try: ,r_0=[0.1,0.1,0.1])
 
 aos = pyxao.SCFeedBackAO(dm,wfs,atm,image_ixs=1)
-aos.find_response_matrix()
-aos.compute_reconstructor(threshold=0.1)
+# aos.find_response_matrix()
+# aos.compute_reconstructor(threshold=0.1)
+# np.save('minerva_red_response_matrix', aos.response_matrix)
+# np.save('minerva_red_reconstructor', aos.reconstructor)
 
+# ipdb.set_trace()
+aos.reconstructor = np.load('minerva_red_reconstructor.npy')
+aos.response_matrix = np.load('minerva_red_response_matrix.npy')
 
-wf_sense.add_atmosphere(atm)
-wf_image.add_atmosphere(atm)
+# Want 1600 photons/cm2/s for the LGS. 
+# = photons/m/s = 1600 * 1e4
+# photons in im = 1600 * 1e4 * D * 1/2000
 
 #See if the reconstructor works!
 #sensors,ims  = aos.correct_twice()
 
-# pdb.set_trace()
-psfs = aos.run_loop(plotit=True,dt=0.002,niter=50,mode='integrator',plate_scale_as_px=0.1,psf_ix = 1,nframesbetweenplots=10)
+# ipdb.set_trace()
+psfs = aos.run_loop(plotit=True,dt=0.002,niter=1,mode='open loop',plate_scale_as_px=0.1,psf_ix = 1,nframesbetweenplots=1)
 psf_mn = np.mean(psfs,axis=0)
-#Uncomment the line below instead to see uncorrected seeing.
-# im_mn = aos.run_loop(plotit=True,dt=0.002,nphot=1e4,niter=50,mode='dodgy_damping',dodgy_damping=0.0,gain=0.0,plate_scale_as_px=0.125)[-2]
 
 #Strehl: Regrid in order to sample finely the peak.
 psf_dl = aos.psf_dl(plate_scale_as_px=0.1,psf_ix=1)
 strehl = np.max(ot.utils.regrid_fft(psf_mn,(1024,1024)))/np.max(ot.utils.regrid_fft(psf_dl,(1024,1024)))
 
-sz = psf_mn.shape[0]
-plt.figure()
-plt.imshow(im_mn[sz//2-20:sz//2+20,sz//2-20:sz//2+20],interpolation='nearest', cmap=cm.gist_heat)
-print("Strehl: {0:6.3f}".format(strehl))
+# sz = psf_mn.shape[0]
+# plt.figure()
+# plt.imshow(psf_mn[sz//2-20:sz//2+20,sz//2-20:sz//2+20],interpolation='nearest', cmap=cm.gist_heat)
+# print("Strehl: {0:6.3f}".format(strehl))

pyxao/ao_system.py

@@ -4,7 +4,11 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.ticker import FormatStrFormatter
-import pdb
+try:
+    import ipdb
+except:
+    #The following ling is dodgy, but still enables set_trace()
+    import pdb as ipdb
 import scipy.ndimage as nd
 import scipy.linalg as la
 import matplotlib.cm as cm
@@ -87,6 +91,8 @@ def __init__(self,
         self.wfs = wfs        
         self.atm = atm
         self.dm_poke_scale = dm_poke_scale
+        self.response_matrix = None
+        self.reconstructor = None
 
         #Add an atmosphere to all wavefronts.
         for wf in self.wavefronts:
@@ -95,7 +101,6 @@ def __init__(self,
     #############################################################################################################
     def psf_dl(self, plate_scale_as_px,
         psf_ix = None,
-        band = None,
         crop = True,    # Whether or not to crop the PSF. If set to False, the below two arguments are irrelevant.
         psf_sz_cropped = None,          # Size of the PSF. By default cropped at the 10th Airy ring
         psf_sigma_limit_N_os = TENTH_AIRY_RING,     # Corresponds to 10 Airy rings
@@ -163,7 +168,7 @@ def find_response_matrix(self,mode='onebyone',amplitude=1e-7):
                 #Check that poking in both directions is equivalent!
                 if (np.sum(wfs_plus*wfs_minus)/np.sum(wfs_plus*wfs_plus) > -0.9):
                     print("WARNING: Poking the DM is assymetric!!!")
-                    # pdb.set_trace()
+                    # ipdb.set_trace()
 
                 # Taking the mean response value.
                 self.response_matrix[i] = 0.5*(wfs_plus - wfs_minus).flatten()
@@ -264,8 +269,6 @@ def correct_twice(self,plotit=False):
         return [measurements0,measurements1,measurements2],[im0,im1,im2] 
 
 
-
-
     #############################################################################################################
     def run_loop(self, dt, 
                 mode = 'integrator',
@@ -339,19 +342,30 @@ def run_loop(self, dt,
             else:
                 psf_sz_cropped = psf_sz
 
+        #Make plate_scale_as_px for Nyquist sampling of the full size.
+        #See XXX below for a bug.
+        if not plate_scale_as_px:
+            plate_scale_rad_px = self.wavefronts[psf_ix].wave / (self.wavefronts[psf_ix].sz * self.wavefronts[psf_ix].m_per_px)
+            plate_scale_as_px = np.degrees(plate_scale_rad_px) * 3600
+            plate_scale_as_px_in = None
+        else:
+            plate_scale_as_px_in = plate_scale_as_px
+
         # Arrays to hold images
         psfs_cropped = np.zeros((niter, psf_sz_cropped, psf_sz_cropped))# at the psf_ix wavelength
 
         """ AO Control loop """
         print("Starting the AO control loop with control logic mode '%s'..." % mode)
         for k in range(niter):  
-            #------------------ AO CONTROL ------------------#
-            print("Iteration %d..." % (k+1))
-            # Evolve the atmosphere & update the wavefront fields to reflect the new atmosphere.                
+            #------------------ EVOLVING THE ATMOSPHERE ------------------#
+            if ((k / niter) * 100) % 10 == 0:
+                print("{:d}% done...".format(int((k / niter) * 100)))
+            # Evolve the atmosphere & update the wavefront fields to reflect the new atmosphere.
             self.atm.evolve(dt * k)
             for wf in self.wavefronts:   
                 wf.atm_field() 
 
+            #------------------ AO CONTROL ------------------#
             if mode == 'open loop':
                 # We still measure the wavefront for plotting.
                 self.wfs.sense()
@@ -374,13 +388,14 @@ def run_loop(self, dt,
                     raise UserWarning            
                 coefficients_current = coefficients_next;   # Update the DM coefficients.
 
+            #-------------------- SAVING PSF ----------------------# 
             # Create the PSF            
-            psf = self.wavefronts[psf_ix].image(plate_scale_as_px = plate_scale_as_px) 
-            psf /= sum(psf.flatten())             
-            # Saving the PSF.            
+            # NB if the following line has non-none, then we have a uint 8 error XXX
+            psf = self.wavefronts[psf_ix].image(plate_scale_as_px = plate_scale_as_px_in) 
+            psf /= np.sum(psf.flatten())             
             psfs_cropped[k] = centre_crop(psf, psf_sz_cropped) 
 
-            # Plotting
+            #------------------ PLOTTING ------------------#
             if plotit & ((k % nframesbetweenplots) == 0):
                 if k == 0:
                     axes = []
@@ -405,6 +420,7 @@ def run_loop(self, dt,
                     plt.gca().xaxis.set_major_formatter(FormatStrFormatter('%.2f\"'))
                     plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%.2f\"'))
                     plt.colorbar(plots[-1],fraction=COLORBAR_FRACTION, pad=COLORBAR_PAD)
+                    ipdb.set_trace()
                 else:
                     # Update the plots
                     plot_wfs.set_data(self.wfs.im)  
@@ -416,51 +432,5 @@ def run_loop(self, dt,
 
         return psfs_cropped
 
-    #############################################################################################################
-    def _plot_ao_loop_iteration(self, k, plot_sz_px):
-        """ Plot the WFS detector image, phase and psf in an iteration of the AO loop. """
-        if k == 0:
-            axes = []
-            plots = []
-            plt.rc('text', usetex=True)
-
-            # WFS plot                    
-            fig_wfs = mu.newfigure(2,2)
-            ax_wfs = fig_wfs.add_subplot(111)  # WFS detector image
-            ax_wfs.title.set_text(r'WFS detector')
-            ax_wfs.axis( [0,self.wavefronts[0].sz,0,self.wavefronts[0].sz] )
-            plot_wfs = ax_wfs.imshow(self.wfs.im,interpolation='nearest',cmap=cm.gray)
-
-            fig = mu.newfigure(width=2,height=1)
-        
-            # Phase
-            axes.append(fig.add_subplot(N_science_imgs,2,2*j+1)) 
-            axes[-1].title.set_text(r'Corrected phase ($\lambda$ = %d nm)' % (self.wavefronts[self.image_ixs[j]].wave*1e9))
-            plots.append(axes[-1].imshow(np.angle(self.wavefronts[self.image_ixs[j]].field)*self.wavefronts[self.image_ixs[j]].pupil,interpolation='nearest', cmap=cm.gist_rainbow, vmin=-np.pi, vmax=np.pi))
-            mu.colourbar(plots[-1])
-
-            # Science image
-            axes.append(fig.add_subplot(N_science_imgs,2,2*j+2))  
-            axes[-1].title.set_text(r'Science Image ($\lambda$ = %d nm)' % (self.wavefronts[self.image_ixs[j]].wave*1e9))
-            pad_x = max((imsize - plot_sz_px[0]) // 2, 0)
-            pad_y = max((imsize - plot_sz_px[1]) // 2, 0)
-            im_cropped = psf_science[j,pad_x:pad_x + min(imsize, plot_sz_px[0]),pad_y:pad_y + min(imsize, plot_sz_px[1])]
-            plots.append(axes[-1].imshow(im_cropped,interpolation='nearest', cmap=cm.gist_heat, norm=LogNorm()))
-            mu.colourbar(plots[-1])
-        else:
-            # Update the plots
-            plot_wfs.set_data(self.wfs.im)  
-            fig.suptitle(r'AO-corrected phase and science images, $k = %d, K_i = %.2f, K_{leak} = %.2f$' % (k, K_i, K_leak))
-            for j in range(N_science_imgs):
-                plots[2*j].set_data(np.angle(self.wavefronts[self.image_ixs[j]].field)*self.wavefronts[self.image_ixs[j]].pupil)
-                pad_x = max((imsize - plot_sz_px[0]) // 2, 0)
-                pad_y = max((imsize - plot_sz_px[1]) // 2, 0)
-                im_cropped = im_science[j,pad_x:pad_x + min(imsize, plot_sz_px[0]),pad_y:pad_y + min(imsize, plot_sz_px[1])]
-                plots[2*j+1].set_data(im_cropped)
-
-        plt.draw()
-        plt.pause(0.00001)  # Need this to plot on some machines.
-
-    
 
 

pyxao/atmosphere.py

@@ -4,7 +4,10 @@
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
 import scipy.ndimage as nd
-import pdb
+try:
+    import ipdb
+except:
+    import pdb
 
 class Atmosphere():
     """A model of the atmosphere, including a fixed phase screen. 
@@ -59,7 +62,7 @@ def __init__(self,sz = 1024,
         if ( (len(elevations) != len(v_wind)) |
             (len(v_wind) != len(r_0)) |
             (len(r_0) != len(angle_wind)) ):
-            pdb.set_trace()
+            ipdb.set_trace()
             print("ERROR: elevations, v_wind, r_0 and angle_wind must all be the same length")
             raise UserWarning 
 
@@ -73,7 +76,7 @@ def __init__(self,sz = 1024,
         self.delays0 = np.empty( (len(r_0),sz,sz) )
         for i in range(self.nlayers): 
             self.delays0[i] = ot.kmf(sz,seed) * np.sqrt(6.88*(m_per_px/r_0[i])**(5.0/3.0)) * wave_ref / 2 / np.pi
-        
+
         #Delays in meters at another time.
         self.delays  = self.delays0.copy()
 

pyxao/deformable_mirror.py

@@ -3,7 +3,10 @@
 import numpy as np
 import matplotlib.pyplot as plt
 import scipy.ndimage as nd
-import pdb
+try:
+    import ipdb
+except:
+    import pdb
 
 class DeformableMirror():
     """A deformable mirror.