fixed bug in W0 optimization

fast/fast.py

@@ -293,8 +293,13 @@ def init_pupil_mask(self):
             else:
                 ptype = 'gauss'
 
-            self.pupil = funcs.compute_pupil(self.Npxls, self.dx, self.D_ground, 
+            puptmp = funcs.compute_pupil(self.Npxls, self.dx, self.D_ground, 
                 self.W0, self.obsc_ground, ptype=ptype)
+
+            if self.W0 == "opt":
+                self.pupil, self.W0 = puptmp
+            else:
+                self.pupil = puptmp
 
             # Circ aperture at satellite
             self.pupil_sat = funcs.compute_pupil(32, self.dx_sat, self.D_sat,
@@ -307,8 +312,13 @@ def init_pupil_mask(self):
                 Tx_obsc=self.obsc_ground, ptype=ptype)
 
             # Gaussian aperture at satellite (NOTE hard coded 32 pxls)
-            self.pupil_sat = funcs.compute_pupil(32, self.dx_sat, self.D_sat, 
+            pupsattmp = funcs.compute_pupil(32, self.dx_sat, self.D_sat, 
                 W0=self.W0, Tx_obsc=self.obsc_sat, ptype='gauss')
+
+            if self.W0 == "opt":
+                self.pupil_sat, self.W0 = pupsattmp
+            else:
+                self.pupil_sat = pupsattmp
 
         self.pupil_filter = funcs.pupil_filter(self.freq.main, self.pupil, spline=False)
 

fast/funcs.py

@@ -299,7 +299,7 @@ def compute_pupil(N, dx, Tx, W0=None, Tx_obsc=0, Raxicon=None, ptype='gauss'):
         if W0 == "opt":
             g, opt = optimize_fibre(circ_ap, dx, return_size=True)
             logger.debug(f"Optimised gaussian size: {opt}")
-            return g * circ_ap
+            return g * circ_ap, opt
         else:
             I0 = 2 / (numpy.pi * W0**2)
             return gaussian2d(N, W0/dx/numpy.sqrt(2)) * circ_ap * numpy.sqrt(I0)
@@ -343,11 +343,20 @@ def _opt_func(W):
         return coupling_loss(W, N, pupil, dx)
 
     opt = minimize_scalar(_opt_func, bracket=[size_min, size_max]).x
-
+
+    # the optimization can fail for some (seemingly random) parameter combinations, 
+    # producing a very small value of opt. Try again but change size_max, if still 
+    # doesn't work, then we are stuffed and raise an Exception.
+    if abs(opt) < dx:
+        logger.info("Gaussian mode optimisation failed, trying with different parameters")
+        opt = minimize_scalar(_opt_func, bracket=[size_min, 2*size_max]).x
+        if abs(opt) < dx:
+            raise Exception("Cannot optimise gaussian mode, try changing DX?")
+
     g = gaussian2d(N, opt/dx/numpy.sqrt(2)) * numpy.sqrt(2./(numpy.pi * opt**2))
 
     if return_size:
-        return g, opt
+        return g, numpy.abs(opt)
     else:
         return g