import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt from .io import load_spectrum from .utils import OptimizeResult, setup_matplotlib from .analysis import finds_peak def _thicknesses_scheludko_at_order(wavelengths, intensities, interference_order, refractive_index, intensities_void=None): """ Compute thicknesses vs wavelength for a given interference order. Parameters ---------- wavelengths : array Wavelength values in nm. intensities : array Intensity values. interference_order : int Interference order. refractive_index : array_like (or float) Refractive index. intensities_void : array, optional Intensities of void. Returns ------- thicknesses : array """ if intensities_void is None: Imin = np.min(intensities) else: Imin = intensities_void n = refractive_index m = interference_order I_norm = (np.asarray(intensities) - Imin) / (np.max(intensities) - Imin) prefactor = wavelengths / (2 * np.pi * n) argument = np.sqrt(I_norm / (1 + (1 - I_norm) * (n**2 - 1)**2 / (4 * n**2))) if m % 2 == 0: term1 = (m / 2) * np.pi else: term1 = ((m+1) / 2) * np.pi term2 = (-1)**m * np.arcsin(argument) return prefactor * (term1 + term2) def _Delta(wavelengths, thickness, interference_order, refractive_index): """ Compute the Delta values. Parameters ---------- wavelengths : array Wavelength values in nm. thickness : array_like (or float) Film thickness. interference_order : int Interference order. refractive_index : array_like (or float) Refractive index. Returns ------- ndarray Delta values. """ # ensure that the entries are numpy arrays wavelengths = np.asarray(wavelengths) h = np.asarray(thickness) n = np.asarray(refractive_index) m = interference_order # Calculation of p as a function of the parity of m if m % 2 == 0: p = m / 2 else: p = (m + 1) / 2 # Calculation of alpha alpha = ((n**2 - 1)**2) / (4 * n**2) # Argument of sinus angle = (2 * np.pi * n * h / wavelengths) - p * np.pi # A = sin²(argument) A = np.sin(angle)**2 # Final calcuation of Delta return (A * (1 + alpha)) / (1 + A * alpha) def _Delta_fit(xdata, thickness, interference_order): """ Wrapper on Delta() for curve_fit. Parameters ---------- xdata : tuple (wavelengths, refractive_index) thickness : array_like (or float) Film thickness. interference_order : int Interference order. Returns ------- ndarray Delta values. """ lambdas, r_index = xdata return _Delta(lambdas, thickness, interference_order, r_index) def get_default_start_stop_wavelengths(wavelengths, intensities, refractive_index, min_peak_prominence, plot=None): """ Returns the start and stop wavelength values of the last monotonic branch. Parameters ---------- wavelengths : array Wavelength values in nm. intensities : array Intensity values. refractive_index : scalar, optional Value of the refractive index of the medium. min_peak_prominence : scalar Required prominence of peaks. plot : bool, optional Display a curve, useful for checking or debuging. The default is None. Raises ------ RuntimeError if at least one maximum and one minimum are not detected. Returns ------- wavelength_start : scalar wavelength_stop : scalar """ # idx_min idx max idx_peaks_min, idx_peaks_max = finds_peak(wavelengths, intensities, min_peak_prominence=min_peak_prominence, plot=plot) failure, message = False, '' if len(idx_peaks_min) == 0: message += 'Failed to detect at least one minimum. ' failure = True if len(idx_peaks_max) == 0: message += 'Failed to detect at least one maximum. ' failure = True if failure: raise RuntimeError(message) # Get the last oscillation peaks lambda_min = wavelengths[idx_peaks_min[-1]] lambda_max = wavelengths[idx_peaks_max[-1]] # Order them wavelength_start = min(lambda_min, lambda_max) wavelength_stop = max(lambda_min, lambda_max) return wavelength_start, wavelength_stop def thickness_from_scheludko(wavelengths, intensities, refractive_index, wavelength_start=None, wavelength_stop=None, interference_order=None, intensities_void=None, plot=None): """ Compute the film thickness based on Scheludko method. Parameters ---------- wavelengths : array Wavelength values in nm. intensities : array Intensity values. refractive_index : scalar, optional Value of the refractive index of the medium. wavelength_start : scalar, optional Starting value of a monotonic branch. Mandatory if interference_order != 0. wavelength_start : scalar, optional Stoping value of a monotonic branch. Mandatory if interference_order != 0. interference_order : scalar, optional Interference order, zero or positive integer. If set to None, the value is guessed. intensities_void : array, optional Intensity in absence of a film. Mandatory if interference_order == 0. plot : bool, optional Display a curve, useful for checking or debuging. The default is None. Returns ------- results : Instance of `OptimizeResult` class. The attribute `thickness` gives the thickness value in nm. """ if plot: setup_matplotlib() if interference_order is None or interference_order > 0: if wavelength_stop is None or wavelength_start is None: raise ValueError('wavelength_start and wavelength_stop must be passed for interference_order != 0.') else: if wavelength_start > wavelength_stop: raise ValueError('wavelength_start and wavelength_stop are swapped.') r_index = refractive_index # Handle the interference order if interference_order is None: # A bit extreme... max_tested_order = 12 # mask input data mask = (wavelengths >= wavelength_start) & (wavelengths <= wavelength_stop) wavelengths_masked = wavelengths[mask] r_index_masked = r_index[mask] intensities_masked = intensities[mask] min_difference = np.inf thickness_values = None if plot: plt.figure() plt.ylabel(r'$h$ ($\mathrm{{nm}}$)') plt.xlabel(r'$\lambda$ ($ \mathrm{nm} $)') for _order in range(0, max_tested_order+1): h_values = _thicknesses_scheludko_at_order(wavelengths_masked, intensities_masked, _order, r_index_masked) difference = np.max(h_values) - np.min(h_values) print(f"h-difference for m={_order}: {difference:.1f} nm") if difference < min_difference: min_difference = difference interference_order = _order thickness_values = h_values if plot: plt.plot(wavelengths_masked, h_values, '.', markersize=3, label=f"Order={_order}") elif interference_order == 0: min_peak_prominence = 0.02 peaks_min, peaks_max = finds_peak(wavelengths, intensities, min_peak_prominence=min_peak_prominence, plot=plot) if len(peaks_max) != 1: raise RuntimeError('Failed to detect a single maximum peak.') lambda_unique = wavelengths[peaks_max[0]] # keep rhs from the maximum mask = wavelengths >= lambda_unique wavelengths_masked = wavelengths[mask] r_index_masked = r_index[mask] intensities_masked = intensities[mask] intensities_I_min_masked = intensities_void[mask] interference_order = 0 thickness_values = _thicknesses_scheludko_at_order(wavelengths_masked, intensities_masked, interference_order, r_index_masked, Imin=intensities_I_min_masked) elif interference_order > 0: h_values = _thicknesses_scheludko_at_order(wavelengths_masked, intensities_masked, interference_order, r_index_masked) thickness_values = h_values else: raise ValueError('interference_order must be >= 0.') # Compute the thickness for the selected order # Delta if interference_order == 0: num = intensities_masked - np.min(intensities_I_min_masked) denom = np.max(intensities_masked) - np.min(intensities_I_min_masked) else: num = intensities_masked - np.min(intensities_masked) denom = np.max(intensities_masked) - np.min(intensities_masked) Delta_from_data = num / denom # Delta_from_data = (intensities_masked -np.min(intensities_masked))/(np.max(intensities_masked) -np.min(intensities_masked)) # Delta_from_data = (intensities_raw_masked -np.min(intensities_raw_masked))/(np.max(intensities_raw_masked) -np.min(intensities_raw_masked)) DeltaScheludko = _Delta(wavelengths_masked, np.mean(thickness_values), interference_order, r_index_masked) xdata = (wavelengths_masked, r_index_masked) popt, pcov = curve_fit(lambda x, h: _Delta_fit(x, h, interference_order), xdata, Delta_from_data, p0=[np.mean(thickness_values)]) fitted_h = popt[0] std_err = np.sqrt(pcov[0][0]) if plot: Delta_values = _Delta(wavelengths_masked, fitted_h, interference_order, r_index_masked) plt.figure() plt.plot(wavelengths_masked, Delta_from_data, 'bo-', markersize=2, label=r'$\mathrm{{Smoothed}}\ \mathrm{{Data}}$') # Scheludko label = rf'$\mathrm{{Scheludko}}\ (h = {np.mean(thickness_values):.1f} \pm {np.std(thickness_values):.1f}\ \mathrm{{nm}})$' plt.plot(wavelengths_masked, DeltaScheludko, 'go-', markersize=2, label=label) # Fit label = rf'$\mathrm{{Fit}}\ (h = {fitted_h:.1f}\pm {std_err:.1f} \ \mathrm{{nm}})$' plt.plot(wavelengths_masked, Delta_values, 'ro-', markersize=2, label=label) plt.legend() plt.ylabel(r'$\Delta$') plt.xlabel(r'$\lambda$ ($ \mathrm{nm} $)') import inspect plt.title(inspect.currentframe().f_code.co_name) return OptimizeResult(thickness=fitted_h, stderr=std_err)