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( '===========================================\