From fb0e5482437b84d1976880d2e8b1474e9cd91a6e Mon Sep 17 00:00:00 2001 From: linarphy Date: Fri, 16 Jun 2023 17:44:17 +0200 Subject: [PATCH] Add intensity calibration - Bump version - Add intensity calibration - Clean wavelength calibration (WIP) - Add --intensity option - Update --calibration option to --wavelength --- spectrum.py | 148 +++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 123 insertions(+), 25 deletions(-) diff --git a/spectrum.py b/spectrum.py index e91cb1e..612d2d1 100644 --- a/spectrum.py +++ b/spectrum.py @@ -8,7 +8,7 @@ import shelve import pathlib from scipy.ndimage import rotate -cache , filename , output , calibration , verbose , no_cache = '' , None , None , None , False, False +cache , filename , output , calibration , intensity_calibration , verbose , no_cache = '' , None , None , None , None , False, False if len( sys.argv ) < 2: raise Exception( 'spectrum.py: type \'spectrum.py -h\' for more information' ) @@ -39,10 +39,16 @@ while i < len( argv ): argv[ i + 1 ] = '--output=' + argv[ i + 1 ] i += 1 continue - elif arg == '-a': + elif arg == '-w': if i == len( sys.argv ) - 1: - raise Exception( 'spectrum.py: calibration have to take a value' ) - argv[ i + 1 ] = '--calibration=' + argv[ i + 1 ] + raise Exception( 'spectrum.py: wavelength have to take a value' ) + argv[ i + 1 ] = '--wavelength=' + argv[ i + 1 ] + i += 1 + continue + elif arg == '-i': + if i == len( sys.argv ) - 1: + raise Exception( 'spectrum.py: intensity have to take a value' ) + argv[ i + 1 ] = '--intensity=' + argv[ i + 1 ] i += 1 continue else: @@ -50,8 +56,10 @@ while i < len( argv ): if arg[1] == '-': # not elif because arg can change after last if if arg == '--help': print( 'spectrum.py [options...] filename\ - \n -a --calibration calibration file, default to no calibration.\ + \n -w --wavelength wavelength calibration file, default to no calibration.\ \n No calibration means no wavelength interpolation\ + \n -i --intensity intensity calibration file, default to no calibration.\ + \n No calibration means no intensity interpolation\ \n -c --cache use given cache\ \n -h --help show this help and quit\ \n -n --no-cache do not use cache and rewrite it\ @@ -62,7 +70,7 @@ while i < len( argv ): \nParse a naroo ETA fits' ) exit() elif arg == '--version': - print( '0.2' ) + print( '0.3' ) exit() elif arg == '--verbose': verbose = True @@ -72,8 +80,10 @@ while i < len( argv ): cache = arg[ 8 : ] elif len( arg ) > 9 and arg[ : 9 ] == '--output=': output = arg[ 9 : ] - elif len( arg ) > 14 and arg[ : 14 ] == '--calibration=': - calibration = arg[ 14 : ] + elif len( arg ) > 13 and arg[ : 13 ] == '--wavelength=': + calibration = arg[ 13 : ] + elif len( arg ) > 12 and arg[ : 12 ] == '--intensity=': + intensity = arg[ 12 : ] else: raise Exception( 'spectrum.py: unknown argument "' + arg + '", type \'ETA.py -h\' for more information' ) else: @@ -89,7 +99,8 @@ if verbose: print( f'spectrum.py: launching now with parameters:\ \n --filename: {filename}\ \n --cache: {cache} ( default: \'\' )\ - \n --calibration: {calibration} ( default to None )\ + \n --wavelength: {calibration} ( default to None )\ + \n --intensity: {intensity} ( default to None )\ \n --output: {output} ( default to None )\ \n --verbose: True ( default to False)\ \n\ @@ -471,10 +482,52 @@ else: Calibration """ +wavelengths = np.arange( spectrum[ 'x' ][ 'max' ] - spectrum[ 'x' ][ 'min' ] ) + if calibration != None: if verbose: - print( 'starting calibration' ) + print( 'starting wavelength calibration' ) + mean_data = np.mean( data[ + spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ], + spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ] + ] , axis = 0 ) + abciss = np.arange( len( mean_data ) ) + + ref = np.array( [ + 6562.79, + 4861.35, + 4340.47, + 4101.73, + 3970.08, + 3889.06, + 3835.40, +# 3646 , + ] ) * 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 ) + + 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.array( lines[ : len( ref ) - 1 ] ) # Balmer discontinuity + + ref = ref[ 1 : ] # start with H-beta + + wavelength_polyval = np.polyfit( lines , ref , 1 ) + + wavelengths = np.polyval( wavelength_polyval , abciss ) + + """ mean_up = np.mean( data[ calibrations[ 'top' ][ 'y' ][ 'min' ] : calibrations[ 'top' ][ 'y' ][ 'max' ], calibrations[ 'top' ][ 'x' ][ 'min' ] : claibrations[ 'top' ][ 'x' ][ 'max' ] @@ -499,33 +552,78 @@ if calibration != None: ) , dtype = int, ) + """ if verbose: - print( 'calibration finished' ) + print( 'wavelength calibration finished' ) + +mean_data = np.mean( data[ + spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ], + spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ] +] , axis = 0 ) + +if intensity != None: + if verbose: + print( 'starting intensity calibration' ) + + intensity_file = pathlib.Path( intensity ) + + with shelve.open( str( intensity_file ) ) as storage: + intensity_stairs = storage[ 'data' ] + intensity_wavelengths = storage[ 'wavelength' ] + + 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( + intensity_wavelengths, + intensity_wavelength , + ) + index_stair = np.where( + intensity_stairs[ + : , + intensity_index, + 0 , + ] < intensity_value + )[0][-1] + stair_before = intensity_stairs[ index_stair , intensity_index ] + stair_after = intensity_stairs[ index_stair - 1 , intensity_index ] + step = 0.2 + start_step = step * ( stairs.shape[0] - index_stair - 1 ) + + index_polyval = np.polyfit( + [ index_stair , index_stair - 1 ], + [ stair_before , stair_after ] + ) + + stair_polyval = np.polyfit( + [ stair_before , stair_after ], + [ 0 , step ] + ) + + final_intensity[index] = start_step + np.polyval( # the step fraction + stair_polyval, + np.polyval( index_polyval , intensity_index ), # the stair "fraction" + ) + + plt.plot( final_intensity ) + plt.show() + + if verbose: + print( 'intensity calibration finished' ) if verbose: print( 'starting output writing' ) if output == None: - print( np.mean( - data[ - spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ], - specturm[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ], - ] , - axis = 0, - ) ) + print( mean_data ) else: if verbose: print( 'storing result in ' + output ) output_file = pathlib.Path( output ) with shelve.open( str( output_file ) ) as output: - output[ 'data' ] = np.mean( - data[ - spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ], - spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ], - ] , - axis = 0, - ) + output[ 'data' ] = mean_data if verbose: print( 'output writing finished' ) print( '===========================================\