Fluence calibration

Add routine to calculate fluence from laser power and beam size measurements

src/toolbox_scs/misc/laser_utils.py

@@ -7,7 +7,8 @@ __all__ = [
 import numpy as np
 import matplotlib.pyplot as plt
 
-def positionToDelay(pos, origin=0, invert = False, reflections=1):
+
+def positionToDelay(pos, origin=0, invert=True, reflections=1):
     ''' converts a motor position in mm into optical delay in picosecond
         Inputs:
             pos: array-like delay stage motor position
@@ -18,29 +19,32 @@ def positionToDelay(pos, origin=0, invert = False, reflections=1):
         Output:
             delay in picosecond
     '''
-    c_ = 299792458 *1e-9 # speed of light in mm/ps
+    c_ = 299792458 * 1e-9  # speed of light in mm/ps
     x = -1 if invert else 1
     return 2*reflections*(pos-origin)*x/c_
-    
+
+
 def degToRelPower(x, theta0=0):
     ''' converts a half-wave plate position in degrees into relative power
         between 0 and 1.
         Inputs:
             x: array-like positions of half-wave plate, in degrees
             theta0: position for which relative power is zero
-            
+
         Output:
             array-like relative power
     '''
     return np.sin(2*(x-theta0)*np.pi/180)**2
 
+
 def fluenceCalibration(hwp, power_mW, npulses, w0x, w0y=None,
                        train_rep_rate=10, fit_order=1,
                        plot=True, xlabel='HWP [%]'):
     """
-    Given a series of relative powers or half wave plate angles and
-    averaged powers in mW, this routine calculates the corresponding
+    Given a measurement of relative powers or half wave plate angles
+    and averaged powers in mW, this routine calculates the corresponding
     fluence and fits a polynomial to the data.
+
     Parameters
     ----------
     hwp: array-like (N)
@@ -50,9 +54,9 @@ def fluenceCalibration(hwp, power_mW, npulses, w0x, w0y=None,
     npulses: int
         number of pulses per train during power measurement
     w0x: float
-        radius at 1/e^2 in x-axis
+        radius at 1/e^2 in x-axis in meter
     w0y: float, optional
-        radius at 1/e^2 in y-axis. If None, w0y=w0x is assumed.
+        radius at 1/e^2 in y-axis in meter. If None, w0y=w0x is assumed.
     train_rep_rate: float
         repetition rate of the FEL, by default equals to 10 Hz.
     fit_order: int
@@ -66,49 +70,42 @@ def fluenceCalibration(hwp, power_mW, npulses, w0x, w0y=None,
     ------
     F: ndarray (N)
         fluence in mJ/cm^2
-    fit: ndarray
-        coefficients of the polynomial fit
+    fit_F: ndarray
+        coefficients of the fluence polynomial fit
+    E: ndarray (N)
+        pulse energy in microJ
+    fit_E: ndarray
+        coefficients of the fluence polynomial fit
     """
     power = np.array(power_mW)
     hwp = np.array(hwp)
-    E = power/(train_rep_rate*npulses)*1e-3 # pulse energy in J
+    E = power/(train_rep_rate*npulses)*1e-3  # pulse energy in J
     if w0y is None:
         w0y = w0x
-    F = 2*E/(np.pi*w0x*w0y) # fluence in J/m^2
-    
+    F = 2*E/(np.pi*w0x*w0y)  # fluence in J/m^2
+
     fit_E = np.polyfit(hwp, E*1e6, fit_order)
     fit_F = np.polyfit(hwp, F*1e-1, fit_order)
     x = np.linspace(hwp.min(), hwp.max(), 100)
     if plot:
-        plt.figure(figsize=(7,4))
-        ax = plt.subplot(121)
-        plt.title(f'w0x = {w0x*1e6:.0f} $\mu$m, w0y = {w0y*1e6:.0f} $\mu$m')
-        plt.plot(hwp, E*1e6, 'o', label='data')
-        fit_label=''
-        for i in range(len(fit_E)-1, 1, -1):
-            fit_label += f'{fit_E[i]:.2g}x$^{i}$ + '
-            if i%2 == 0:
-                fit_label +='\n'
-        fit_label += f'{fit_E[-2]:.2g}x + {fit_E[-1]:.2g}'
-        plt.plot(x, np.poly1d(fit_E)(x), label=fit_label)
-        plt.ylabel('Pulse energy [$\mu$J]')
-        plt.xlabel('HWP [%]')
-        plt.legend()
-        plt.grid()
-        plt.subplot(122, sharex=ax)
-        plt.title(f'w0x = {w0x*1e6:.0f} $\mu$m, w0y = {w0y*1e6:.0f} $\mu$m')
-        plt.plot(hwp, F*1e-1, 'o', label='data')
-        fit_label=''
+        fig, ax = plt.subplots(figsize=(6, 4))
+        ax.set_title(f'w0x = {w0x*1e6:.0f} $\mu$m, w0y = {w0y*1e6:.0f} $\mu$m')
+        ax.plot(hwp, F*1e-1, 'o', label='data')
+        fit_label = 'F = '
         for i in range(len(fit_F)-1, 1, -1):
             fit_label += f'{fit_F[i]:.2g}x$^{i}$ + '
-            if i%2 == 0:
-                fit_label +='\n'
+            if i % 2 == 0:
+                fit_label += '\n'
         fit_label += f'{fit_F[-2]:.2g}x + {fit_F[-1]:.2g}'
-        #fit_label = '+'.join([f'{f:.2g}x$^{i}$' for i, f in enumerate(fit_F[:-1])])
-        plt.plot(x, np.poly1d(fit_F)(x), label=fit_label)
-        plt.ylabel('Fluence [mJ/cm$^2$]')
-        plt.xlabel('HWP [%]')
-        plt.legend()
-        plt.grid()
-        plt.tight_layout()
-    return F*1e-1, fit_F, E*1e6, fit_E
+        ax.plot(x, np.poly1d(fit_F)(x), label=fit_label)
+        ax.set_ylabel('Fluence [mJ/cm$^2$]')
+        ax.set_xlabel(xlabel)
+        ax.legend()
+        ax.grid()
+        ax2 = ax.secondary_yaxis(
+            'right',
+            functions=(lambda x: 1e7*x*np.pi*w0x*w0y/2,
\ No newline at end of file
+                       lambda x: 1e-7*2*x/(np.pi*w0x*w0y)))
+        ax2.set_ylabel(r'Pulse energy [$\mu$J]')
+        fig.tight_layout()
+    return F*1e-1, fit_F, E*1e6, fit_E