Update & fix intensity calibration

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
+        indexes_stair_lower = np.where( # stairs lower than intensity value
-        index_stair = np.where(
             intensity_stairs[
                 :              ,
                 intensity_index,
-                0              ,
+                0
             ] < intensity_value
-        )
+        )[0]
-        if len( index_stair[0] ) == 0:
-            index_stair = intensity_stairs.shape[0] - 1
-        else:
-            index_stair = index_stair[0][0]
 
-        stair_before = intensity_stairs[ index_stair     , intensity_index , 0 ]
+        if len( indexes_stair_lower ) == 0: # intensity value outside ETA (below)
-        stair_after  = intensity_stairs[ index_stair - 1 , intensity_index , 0 ]
+            final_intensity[ index ] = 0
-        step         = 0.2
+            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      ],
+            [ stair_intensity[ 'higher' ] , stair_intensity[ 'lower' ] ],
-            [ index_stair  , index_stair - 1  ],
+            [     index_stair[ 'higher' ] ,     index_stair[ 'lower' ] ],
-            1                                  ,
+            1                                                           ,
         )
-        true_intensity_value = ( intensity_stairs.shape[0] - np.polyval( index_polyval , intensity_value ) - 1 ) * step
+
-        if true_intensity_value < 0:
+        true_intensity_value = ( intensity_stairs.shape[0] - np.polyval( index_polyval , intensity_value ) ) * step
-            true_intensity_value = 0
 
         final_intensity[index] = true_intensity_value