Add intensity calibration

- Bump version
- Add intensity calibration
- Clean wavelength calibration (WIP)
- Add --intensity option
- Update --calibration option to --wavelength

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' )
+                    raise Exception( 'spectrum.py: wavelength have to take a value' )
-                argv[ i + 1 ] = '--calibration=' + argv[ i + 1 ]
+                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=':
+            elif len( arg ) > 13 and arg[ : 13 ] == '--wavelength=':
-                calibration = arg[ 14 : ]
+                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(
+    print( mean_data )
-        data[
-            spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ],
-            specturm[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ],
-        ]       ,
-        axis = 0,
-    ) )
 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(
+        output[ 'data' ] = mean_data
-            data[
-                spectrum[ 'y' ][ 'min' ] : spectrum[ 'y' ][ 'max' ],
-                spectrum[ 'x' ][ 'min' ] : spectrum[ 'x' ][ 'max' ],
-            ]       ,
-            axis = 0,
-        )
 if verbose:
     print( 'output writing finished' )
     print( '===========================================\