diff --git a/spectrum.py b/spectrum.py index ac46b52..d1973e2 100644 --- a/spectrum.py +++ b/spectrum.py @@ -566,7 +566,7 @@ if verbose: 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: @@ -578,43 +578,110 @@ if intensity != None: intensity_stairs = storage[ 'data' ] intensity_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 ), + ), + ) + ] + intensity_wavelengths = intensity_wavelengths[ # remove intensity_wavelengths outside range + np.where( + np.logical_and( + intensity_wavelengths > np.min( wavelengths ), + intensity_wavelengths < np.max( wavelengths ), + ), + ) + ] + + if len( wavelengths ) == 0: + 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( + intensity_index = utils.near_value( # list of index corresponding to index for intensity_stairs intensity_wavelengths, intensity_wavelength , ) - if len( intensity_index ) != 1: - break + + if len( intensity_index ) != 1: # too much or no intensity found near value + final_intensity[ index ] = - 1 + continue intensity_index = intensity_index[0] - if abs( intensity_wavelength - intensity_wavelengths[ intensity_index ] ) > 5e-10: - break - index_stair = np.where( + + indexes_stair_lower = np.where( # stairs lower than intensity value intensity_stairs[ : , intensity_index, - 0 , + 0 ] < intensity_value - ) - if len( index_stair[0] ) == 0: - index_stair = intensity_stairs.shape[0] - 1 - else: - index_stair = index_stair[0][0] + )[0] - stair_before = intensity_stairs[ index_stair , intensity_index , 0 ] - stair_after = intensity_stairs[ index_stair - 1 , intensity_index , 0 ] - step = 0.2 + 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 ] = intensity_stairs.shape[0] * step + 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_before , stair_after ], - [ index_stair , index_stair - 1 ], - 1 , + [ stair_intensity[ 'higher' ] , stair_intensity[ 'lower' ] ], + [ index_stair[ 'higher' ] , index_stair[ 'lower' ] ], + 1 , ) - true_intensity_value = ( intensity_stairs.shape[0] - np.polyval( index_polyval , intensity_value ) - 1 ) * step - if true_intensity_value < 0: - true_intensity_value = 0 + + true_intensity_value = ( intensity_stairs.shape[0] - np.polyval( index_polyval , intensity_value ) ) * step final_intensity[index] = true_intensity_value