Fixed factors of 2pi, dPhidT calculation

Author: Thomas Satterthwaite

sotodlib/toast/ops/mumux_crosstalk_util.py

@@ -181,6 +181,11 @@ def pos_to_chi(focalplane, dets, alpha=9.64e-30, tol=1.0, collision=1.0):
                     is_north_subset.append(None)
                     mux_band_subset.append(None)
                     bond_pad_subset.append(None)
+        # Convert to numpy objects
+        freq_subset     = np.array(freq_subset)
+        is_north_subset = np.array(is_north_subset)
+        mux_band_subset = np.array(mux_band_subset)
+        bond_pad_subset = np.array(bond_pad_subset)
         # Compute chi for all detector pairs in this subset
         ndet = len(detector_subset)
         for idet1, det1 in enumerate(detector_subset):
@@ -211,6 +216,14 @@ def pos_to_chi(focalplane, dets, alpha=9.64e-30, tol=1.0, collision=1.0):
                     continue
                 if np.abs(bond_pad1 - bond_pad2) != 4:
                     continue
+                # Check that this truly is the nearly frequency neighbor
+                neighbors = np.argwhere((is_north_subset == is_north1) & \
+                                        (mux_band_subset == mux_band1)).flatten()
+                freq_neighbors = np.argsort(np.abs(freq_subset[neighbors] - freq1)).flatten()
+                if not idet2 in neighbors[freq_neighbors[1:3]]:
+                    # freq_neighbors[0] is idet1
+                    # frequency neighbors on either side
+                    continue
                 x2, y2 = position_subset[idet2]
                 # Translate frequencies to chi
                 df = freq1 - freq2

sotodlib/toast/ops/sim_mumux_crosstalk.py

@@ -81,6 +81,8 @@
 R_BOLO = 0.008
 # Readout noise fraction
 R_FRAC = 0.5
+# Shunt Resistance [Ohm]
+R_SHUNT = 400e-6
 
 
 @trait_docs
@@ -152,7 +154,7 @@ def _squid_phase_to_temperature(self, input_signal, dPhi0dT):
         output_signal = input_signal / dPhi0dT
         return output_signal
 
-    def _draw_Phi0(self, obs, focalplane, detectors, vmin=0.3, vmax=1.3):
+    def _draw_Phi0(self, obs, focalplane, detectors, vmin=0.0, vmax=1.0):
         """ Draw initial SQUID phases from a flat distribution
         """
         Phi0 = {}
@@ -173,6 +175,8 @@ def _draw_Phi0(self, obs, focalplane, detectors, vmin=0.3, vmax=1.3):
                 counter=(counter1, counter2),
             )[0]
             v = vmin + x * (vmax - vmin)
+            # Convert from Phi0 to Phi
+            v *= 2 * np.pi
             Phi0[det] = v
 
         return Phi0
@@ -209,9 +213,12 @@ def _evaluate_dPhi0dT(self, obs, signal, detectors, rows, Phi0):
             P_atm = efficiency * bandpass.optical_loading(det, median_signal)  # W
             P_opt += P_atm - P_atm_ref
             dPdT = bandpass.kcmb2w(det)  # K_CMB -> W
-            dIdP = 1 / np.sqrt((P_sat - P_opt) * R_FRAC * R_BOLO)  # W -> A
+            R_TES = R_FRAC * R_BOLO
+            I_TES = np.sqrt((P_sat - P_opt) / R_TES)
+            dIdP = -1 / (I_TES * (R_TES - R_SHUNT))  # W -> A
             dPhi0dI = 1 / 9e-6  # A -> [rad]
-            dPhi0dT[det] = dPdT * dIdP * dPhi0dI  # K_CMB -> [rad]
+            dPhidPhi0 = 2 * np.pi
+            dPhi0dT[det] = dPdT * dIdP * dPhi0dI * dPhidPhi0  # K_CMB -> [rad]
 
         return dPhi0dT
 
@@ -290,10 +297,11 @@ def _exec(self, data, detectors=None, **kwargs):
                         crosstalk += chi * np.sin(
                             source_squid_phase - target_squid_phase
                         )
-                    else:
-                        # Otherwise flag
-                        temp_obs.detdata[self.det_flags][det_source] |= self.det_flag_mask
-                        temp_obs.detdata[self.det_flags][det_target] |= self.det_flag_mask
+                    #else:
+                    #    # Otherwise flag
+                    #    print('Flagging!')
+                    #    temp_obs.detdata[self.det_flags][det_source] |= self.det_flag_mask
+                    #    temp_obs.detdata[self.det_flags][det_target] |= self.det_flag_mask
 
                 # Translate crosstalk into temperature units and scale to
                 # match input data