diff --git a/wavelength_calibration.py b/wavelength_calibration.py index d09dcfc..0965ce7 100644 --- a/wavelength_calibration.py +++ b/wavelength_calibration.py @@ -1,30 +1,10 @@ import numpy as np import scipy.signal as sig -import matplotlib.backends.backend_qtagg as qt -import matplotlib.figure as fig -import utils - -calib_top = np.load( 'asset/calib_top.npy' ) -calib_down = np.load( 'asset/calib_down.npy' ) - -peaks_reference_sorted = np.loadtxt( 'calibration/OHP_Hg.calib' ) -peaks_reference = np.sort( peaks_reference_sorted ) - -mean_calib_top = np.mean( calib_top , axis = 0 ) -mean_calib_down = np.mean( calib_down , axis = 0 ) - -peaks_calib_top = np.array( - utils.retrieve_peaks( mean_calib_top ) -).astype( int ) -peaks_calib_down = np.array( - utils.retrieve_peaks( mean_calib_down ) -).astype( int ) - -""" -Signal without peaks -""" def remove_peaks( signal , peaks ): + """ + remove peaks from a signal + """ peakless_signal = signal.copy() for peak in peaks: first = peak @@ -41,27 +21,30 @@ def remove_peaks( signal , peaks ): peakless_signal[ first : peak ] = peakless_signal[ peak ] * np.ones( peak - first ) return sig.medfilt( peakless_signal , 111 ) -peakless_calib_top = remove_peaks( mean_calib_top , peaks_calib_top ) -def get_extremities( signal ): +def get_extremities( signal , peaks , temp ): """ It's possible to have an idea of the part the spectrum begins with the small large peak at the start of it, the peak at 3810 A should be inside it. The same goes for the area at the end. ONLY TRUE FOR THE CURRENT CALIBRATION (Hg) """ - argmin = np.argmin( - signal[ : len( signal ) // 2 ] - ) - argmax = np.argmax( peakless_calib_top[ : argmin ] ) - peaks_inside = np.where( - np.logical_and( - argmax < peaks_calib_top, - argmin > peaks_calib_top, + peaks_inside , i = [] , 0 + argmin = len( signal ) // 2 + while len( peaks_inside ) != 1: + if i > 50 or len( signal[ : argmin ] ) == 0: + raise Exception( 'unknown plage, cannot autocalibrate' ) + argmin = np.argmin( + signal[ : argmin ] ) - )[0] - if len( peaks_inside ) < 1: - raise Exception( 'unknown plage, cannot autocalibrate' ) + argmax = np.argmax( signal[ : argmin ] ) + peaks_inside = np.where( + np.logical_and( + argmax < peaks, + argmin > peaks, + ) + )[0] + i += 1 first_peak = peaks_inside[0] """ @@ -70,85 +53,48 @@ def get_extremities( signal ): ONLY TRUE FOR THE CURRENT CALIBRATION (Hg) """ argmin_1 = np.argmin( - peakless_calib_top[ len( peakless_calib_top ) // 2 : ] - ) + len( peakless_calib_top ) // 2 - + signal[ + len( signal ) // 2 : + - int( + 0.1 * len( signal ) // 2 + ) + ] # not at the end + ) + len( signal ) // 2 peaks_inside = np.where( - argmin_1 < peaks_calib_top + argmin_1 < peaks, )[0] if len( peaks_inside ) < 1: raise Exception( 'unknown plage, cannot autocalibrate' ) return ( first_peak , peaks_inside[0] ) -first_peak_ref = 3 # hard-coded -last_peak_ref = 20 # hard-coded -first_peak_cal , last_peak_cal = get_extremities( peakless_calib_top ) +def only_keep_calib( peaks_data , peaks_calib ): + """ + only keep data peaks corresponding to calibration + """ + diff_calib = ( peaks_calib[ 1 : ] - peaks_calib[ : -1 ] ).astype( float ) + diff_calib -= np.min( diff_calib ) + diff_calib /= np.max( diff_calib ) -""" -delete if too much peak in calib -""" + diff_data = ( peaks_data[ 1 : ] - peaks_data[ : -1 ] ).astype( float ) + diff_data -= np.min( diff_data ) + diff_data /= np.max( diff_data ) -peaks_calib_top = peaks_calib_top[ first_peak_cal : last_peak_cal ] -peaks_reference = peaks_reference[ first_peak_ref : last_peak_ref ] + peaks = [ -1 ] + for i in range( len( diff_calib ) ): + good , sum_ = -1 , 0 + for j in range( peaks[ - 1 ] + 1 , len( diff_data ) ): + sum_ += diff_data[ j ] + if sum_ - diff_calib[ i ] > 0.002: + print( sum_ - diff_calib[ i ] ) + raise Exception( 'reference peak not found' ) + if sum_ - diff_calib[ i ] > - 0.002: + good = j + break + if good == -1: + raise Exception( 'reference peak not found and not exceeded' ) + peaks.append( good ) + peaks.append( peaks[-1] + 1 ) # append the last peak -diff_ref = ( peaks_reference[ 1 : ] - peaks_reference[ : -1 ] ).astype( float ) -diff_ref -= np.min( diff_ref ) -diff_ref /= np.max( diff_ref ) -diff_cal = ( peaks_calib_top[ 1 : ] - peaks_calib_top[ : -1 ] ).astype( float ) -diff_cal -= np.min( diff_cal ) -diff_cal /= np.max( diff_cal ) - -peaks = [ -1 ] -for i in range( len( diff_ref ) ): - good , sum_ = -1 , 0 - for j in range( peaks[ - 1 ] + 1 , len( diff_cal ) ): - sum_ += diff_cal[ j ] - if sum_ - diff_ref[ i ] > 0.002: - print( sum_ - diff_ref[ i ] ) - raise Exception( 'reference peak not found' ) - if sum_ - diff_ref[ i ] > - 0.002: - good = j - break - if good == -1: - raise Exception( 'reference peak not found and not exceeded' ) - peaks.append( good ) -peaks.append( peaks[-1] + 1 ) - -peaks_calib_top = np.array( - [ peaks_calib_top[ i ] for i in peaks[ 1 : ] ] -).astype( int ) - -""" -Calibration -""" - -polyval_wavelength = np.polyfit( - peaks_calib_top, - peaks_reference, - 1 , -) - -wavelength = np.polyval( - polyval_wavelength, - np.arange( - 0 , - len( mean_calib_top ), - 1 , - ) -) - -manager = qt.FigureManager( - qt.FigureCanvas( - fig.Figure( - figsize = ( 10 , 5 ) - ), - ), - 0, -) -manager.canvas.figure.gca().plot( - peaks_reference, - peaks_reference - np.polyval( polyval_wavelength , peaks_calib_top ), - 'o-' , -) -manager.show() -manager.start_main_loop() + return np.array( + [ peaks_data[ i ] for i in peaks[ 1 : ] ] # remove the first -1 value + ).astype( int )