From 7579a3dde15f4f4f27b8ebf8635696b46b84227e Mon Sep 17 00:00:00 2001 From: linarphy Date: Tue, 18 Jul 2023 14:45:00 +0200 Subject: [PATCH] Update intensity calibration --- spectrum.py | 304 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 182 insertions(+), 122 deletions(-) diff --git a/spectrum.py b/spectrum.py index 7c2d3b2..e2b8589 100755 --- a/spectrum.py +++ b/spectrum.py @@ -3,6 +3,7 @@ import numpy as np import matplotlib.pyplot as plt from scipy.signal import medfilt , find_peaks +from scipy.signal import convolve as sp_convolve from scipy.optimize import curve_fit import utils import sys @@ -11,7 +12,8 @@ import pathlib from scipy.ndimage import rotate from astropy.io import fits -cache , filename , output , calibration , intensity_calibration , verbose , no_cache = None , None , None , None , None , False, False +cache , filename , output , calibration , intensity_calibration +, verbose , no_cache = None , None , None , None , None , False, False if len( sys.argv ) < 2: raise Exception( 'spectrum.py: type \'spectrum.py -h\' for more information' ) @@ -143,7 +145,9 @@ if 'cache' in files and files[ 'cache' ].is_file() and not no_cache: with shelve.open( str( files[ 'cache' ] ) ) as cache: for key in [ 'data' , 'border' , 'calibrations' ]: if key not in cache: - raise Exception( 'spectrum.py: missing data in cache file' ) + raise Exception( + 'spectrum.py: missing data in cache file' + ) data = cache[ 'data' ] border = cache[ 'border'] spectrum = cache[ 'spectrum' ] @@ -282,7 +286,11 @@ else: ] number = 1000 position_peaks = np.zeros( number ) - indexes = np.linspace( border[ 'x' ][ 'min' ] , border[ 'x' ][ 'max' ] , number , dtype = int ) + indexes = np.linspace( + border[ 'x' ][ 'min' ], + border[ 'x' ][ 'max' ], + number , dtype = int , + ) for i in range( number ): x = np.arange( border[ 'y' ][ 'min' ], @@ -338,7 +346,7 @@ else: axis = 1, ) indexes = utils.near_value( - list_ , + list_ , np.mean( list_ ) + ( np.max( list_ ) - np.mean( list_ ) ) / 10, ) + border[ 'y' ][ 'min' ] @@ -381,13 +389,29 @@ else: while len( indexes ) == 2: factor += 1 indexes = utils.near_value( - list_ , - np.min( list_ ) + ( np.mean( list_ ) - np.min( list_ ) ) / factor, + list_ , + np.min( + list_, + ) + ( + np.mean( + list_, + ) - np.min( + list_, + ) + ) / factor, ) factor -= 1 indexes = utils.near_value( - list_ , - np.min( list_ ) + ( np.mean( list_ ) - np.min( list_ ) ) / factor, + list_ , + np.min( + list_, + ) + ( + np.mean( + list_, + ) - np.min( + list_, + ) + ) / factor, ) + border[ 'x' ][ 'min' ] + 100 # valid convolution only spectrum[ 'x' ] = { @@ -418,7 +442,8 @@ else: list_ , height = np.mean( list_ ) + ( 1 - np.mean( list_ ) ) / 2, )[0] - number = 0.01 + abs( len( peaks ) - 90 ) # the lower the better + number = 0.01 + abs( len( peaks ) - 90 ) # the lower the + # better intensity = np.sum( list_[ peaks ] ) return intensity / ( amplitude ** 2 * number ) @@ -439,22 +464,38 @@ else: 'valid' , ) indicators /= np.max( indicators ) - calibration_areas = utils.consecutive( np.where( indicators > 1 / 1000 )[0] ) + calibration_areas = utils.consecutive( + np.where( + indicators > 1 / 1000, + )[0], + ) calibration_areas = [ - [ calibration_area for calibration_area in calibration_areas if ( - calibration_area[0] < ( - border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] - ) / 2 - ) ], - [ calibration_area for calibration_area in calibration_areas if ( - calibration_area[0] > ( - border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] - ) / 2 - ) ], + [ + calibration_area for calibration_area in calibration_areas if ( + calibration_area[0] < ( + border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] + ) / 2 + ) + ], + [ + calibration_area for calibration_area in calibration_areas if ( + calibration_area[0] > ( + border[ 'y' ][ 'max' ] - border[ 'y' ][ 'min' ] + ) / 2 + ) + ], ] calibration_sizes = [ - [ len( calibration_area ) for calibration_area in calibration_areas[0] ], - [ len( calibration_area ) for calibration_area in calibration_areas[1] ], + [ + len( + calibration_area + ) for calibration_area in calibration_areas[0] + ], + [ + len( + calibration_area + ) for calibration_area in calibration_areas[1] + ], ] calibrations_y = [ calibration_areas[0][ @@ -506,7 +547,9 @@ else: Calibration """ -wavelengths = np.arange( spectrum[ 'x' ][ 'max' ] - spectrum[ 'x' ][ 'min' ] ) +wavelengths = np.arange( + spectrum[ 'x' ][ 'max' ] - spectrum[ 'x' ][ 'min' ] +) if calibration != None: if verbose: @@ -530,19 +573,48 @@ if calibration != None: ] ) * 1e-10 start , end = 5000 , 18440 - polyval_before = np.polyfit( abciss[ : start ] , mean_data[ : start ] , 2 ) - polyval_middle = np.polyfit( abciss[ start : end ] , mean_data[ start : end ] , 2 ) - polyval_end = np.polyfit( abciss[ end : ] , mean_data[ end : ] , 2 ) + polyval_before = np.polyfit( + abciss[ : start ] , + mean_data[ : start ], + 2 , + ) + polyval_middle = np.polyfit( + abciss[ start : end ] , + mean_data[ start : end ], + 2 , + ) + polyval_end = np.polyfit( + abciss[ end : ] , + mean_data[ end : ], + 2 , + ) - mean_data[ : start ] = mean_data[ : start ] - np.polyval( polyval_before , abciss[ : start ] ) - mean_data[ start : end ] = mean_data[ start : end ] - np.polyval( polyval_middle , abciss[ start : end ] ) - mean_data[ end : ] = mean_data[ end : ] - np.polyval( polyval_end , abciss[ end : ] ) + mean_data[ : start ] = mean_data[ : start ] - np.polyval( + polyval_before , + abciss[ : start ], + ) + mean_data[ start : end ] = mean_data[ start : end ] - np.polyval( + polyval_middle , + abciss[ start : end ], + ) + mean_data[ end : ] = mean_data[ end : ] - np.polyval( + polyval_end , + abciss[ end : ], + ) mean_data_normalized = mean_data.copy() # normalization mean_data_normalized -= np.min( mean_data_normalized ) mean_data_normalized /= np.max( mean_data_normalized ) - lines = [ np.mean( cons ) for cons in utils.consecutive( np.where( mean_data_normalized < 0.3 )[0] ) ] + lines = [ + np.mean( + cons, + ) for cons in utils.consecutive( + np.where( + mean_data_normalized < 0.3, + )[0], + ) + ] lines = np.array( lines[ : len( ref ) - 1 ] ) # Balmer discontinuity ref = ref[ 1 : ] # start with H-beta @@ -578,17 +650,50 @@ bias = { ), } -print( bias[ 'top' ] , bias[ 'down' ] ) -mean_bias = np.mean( [ bias[ 'top' ] , bias[ 'down' ] ] , axis = 0 ) / 2 -print( 'mean_bias spectrum: ' + str( mean_bias ) ) +mean_bias = sp_convolve( + np.mean( + [ + bias[ 'top' ], + bias[ 'down' ], + ] , + axis = 0, + ) , + np.ones( + ( + 50, + ), + ) , + 'same', +) / 50 + +plt.plot( bias[ 'top' ] , label = 'top' ) +plt.plot( bias[ 'down' ] , label = 'down' ) +plt.plot( mean_bias , label = 'mean' ) +plt.legend() +plt.show() + +data[ + : , + spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ] +] = data[ + :, + spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ] +] - mean_bias if verbose: print( 'bias substraction finished' ) +ETA_bias = np.load( 'ETA_bias.npy' ) + +plt.plot( mean_bias , label = 'spectrum' ) +plt.plot( ETA_bias , label = 'ETA' ) +plt.legend() +plt.show() + mean_data = np.mean( data[ spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ], spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ] -] , axis = 0 ) - mean_bias +] , axis = 0 ) if intensity != None: if verbose: @@ -597,121 +702,76 @@ if intensity != None: intensity_file = pathlib.Path( intensity ) with shelve.open( str( intensity_file ) ) as storage: - intensity_stairs = storage[ 'data' ] - intensity_wavelengths = storage[ 'wavelength' ] * 1e-10 + ETA_stairs = storage[ 'data' ] + ETA_wavelengths = storage[ 'wavelength' ] * 1e-10 wavelengths = wavelengths[ # remove wavelengths outside range np.where( np.logical_and( - wavelengths > np.min( intensity_wavelengths ), - wavelengths < np.max( intensity_wavelengths ), + wavelengths > np.min( ETA_wavelengths ), + wavelengths < np.max( ETA_wavelengths ), ), ) ] - intensity_wavelengths = intensity_wavelengths[ # remove intensity_wavelengths outside range + ETA_wavelengths = ETA_wavelengths[ # remove intensity_wavelengths + # outside range np.where( np.logical_and( - intensity_wavelengths > np.min( wavelengths ), - intensity_wavelengths < np.max( wavelengths ), + ETA_wavelengths > np.min( wavelengths ), + ETA_wavelengths < np.max( wavelengths ), ), ) ] if len( wavelengths ) == 0: - raise Exception( 'spectrum.py: spectrum and ETA does not share any common wavelengths' ) + raise Exception( + 'spectrum.py: spectrum and ETA does not share any common' + + 'wavelengths' + ) step = 0.2 # depends of ETA source #TODO final_intensity = np.zeros( len( wavelengths ) ) for index in range( len( wavelengths ) ): - intensity_value = mean_data[ index ] - intensity_wavelength = wavelengths[ index ] - intensity_index = utils.near_value( # list of index corresponding to index for intensity_stairs - intensity_wavelengths, - intensity_wavelength , + wavelength = wavelengths[ index ] + current_ETA_index = np.argmin( + ETA_wavelengths - wavelength ) + ETA_wavelength = ETA_wavelengths[ current_ETA_index ] - if len( intensity_index ) != 1: # too much or no intensity found near value - final_intensity[ index ] = - 1 - continue - intensity_index = intensity_index[0] + current_stair = ETA_stairs[ + : , + current_ETA_index, + 0 + ] - indexes_stair_lower = np.where( # stairs lower than intensity value - intensity_stairs[ - : , - intensity_index, - 0 - ] < intensity_value + abciss_intensity_curve = np.arange( + ( + ETA_stairs.shape[0] - 1 # we want 11 value, no 12 + ) * step , + - step / 2, # a little more than - step to remove negative + - step , + ) + function_intensity_curve = lambda x , a , b : a * np.log( # log + # because the curve is inverted (exp on the + # other side ) + x + ) + b + + results = curve_fit( + function_intensity_curve, + current_stair , + abciss_intensity_curve , )[0] - if len( indexes_stair_lower ) == 0: # intensity value outside ETA (below) - final_intensity[ index ] = 0 - continue - if len( indexes_stair_lower ) != intensity_stairs.shape[0] - indexes_stair_lower[0]: # stairs intensity does not decrease with index as it should - # could indicate an artefact in ETA - indexes_stair_lower = [ - int( np.mean( - [ - intensity_stairs.shape[0] - - indexes_stair_lower[0] , - len( indexes_stair_lower ) - ] - ) ) - ] - - indexes_stair_higher = np.where( # stairs higher than intensity value - intensity_stairs[ - : , - intensity_index, - 0 - ] > intensity_value - )[0] - - if len( indexes_stair_higher ) == 0: # intensity value outside ETA (upper) - final_intensity[ index ] = np.exp( - step * intensity_stairs.shape[0] - ) - continue - if len( indexes_stair_higher ) - 1 != indexes_stair_higher[-1]: # stairs intensity does not decrease with index as it should - indexes_stair_higher = [ - int( np.mean( - [ - indexes_stair_higher[ - 1 ] , - len( indexes_stair_higher ) - 1, - ], - ) ), - ] - indexes_stair_higher = [ - indexes_stair_lower[ 0 ] - 1, - ] - - index_stair = { - 'higher': indexes_stair_higher[-1], - 'lower' : indexes_stair_lower[0] , - } - - if index_stair[ 'lower' ] - index_stair[ 'higher' ] != 1: # ETA curve should decrease - raise Exception( 'spectrum.py: given intensity stairs (from ETA) are missformed' ) - - stair_intensity = { - 'higher': intensity_stairs[ index_stair[ 'higher' ] , intensity_index , 0 ], - 'lower' : intensity_stairs[ index_stair[ 'lower' ] , intensity_index , 0 ], - } - - index_polyval = np.polyfit( # fraction stair index from intensity value - [ stair_intensity[ 'higher' ] , stair_intensity[ 'lower' ] ], - [ index_stair[ 'higher' ] , index_stair[ 'lower' ] ], - 1 , + final_intensity[ index ] = np.exp( + function_intensity_curve( + mean_data[ index ], + *results , + ), ) - true_intensity_value = ( intensity_stairs.shape[0] - np.polyval( index_polyval , intensity_value ) ) * step - - final_intensity[index] = np.exp( true_intensity_value ) - - if final_intensity[ index ] > np.exp( 2.2 ): - final_intensity[ index ] = np.exp( 2.2 ) - if verbose: print( 'intensity calibration finished' )