From 32403fdce47b4308b0e93011b8e26c324373bcc5 Mon Sep 17 00:00:00 2001
From: linarphy <linarphy@linarphy.net>
Date: Mon, 19 Jun 2023 17:31:14 +0200
Subject: [PATCH] Add & Fix final format

---
 spectrum.py | 85 ++++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 62 insertions(+), 23 deletions(-)

diff --git a/spectrum.py b/spectrum.py
index 3b57daa..5fc2869 100644
--- a/spectrum.py
+++ b/spectrum.py
@@ -7,6 +7,7 @@ import sys
 import shelve
 import pathlib
 from scipy.ndimage import rotate
+from astropy.io import fits
 
 cache , filename , output , calibration , intensity_calibration , verbose , no_cache = '' , None , None , None , None , False, False
 if len( sys.argv ) < 2:
@@ -107,7 +108,10 @@ if verbose:
            \n===========================================' )
 # TODO: check in advance file to check if exists or writeable
 
-data = utils.load( filename )
+hdul = fits.open( filename )
+data = hdul[0].data
+head = hdul[0].header
+hdul.close()
 if verbose:
     print( 'data loaded' )
 
@@ -570,7 +574,7 @@ if intensity != None:
 
     with shelve.open( str( intensity_file ) ) as storage:
         intensity_stairs      = storage[ 'data'       ]
-        intensity_wavelengths = storage[ 'wavelength' ]
+        intensity_wavelengths = storage[ 'wavelength' ] * 1e-10
 
     final_intensity = np.zeros( len( wavelengths ) )
     for index in range( len( wavelengths ) ):
@@ -580,35 +584,37 @@ if intensity != None:
             intensity_wavelengths,
             intensity_wavelength ,
         )
+        if len( intensity_index ) != 1:
+            break
+        intensity_index = intensity_index[0]
+        if abs( intensity_wavelength - intensity_wavelengths[ intensity_index ] ) > 5e-10:
+            break
         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 ]
+        )
+        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 ]
+        stair_after  = intensity_stairs[ index_stair - 1 , intensity_index , 0 ]
         step         = 0.2
-        start_step   = step * ( stairs.shape[0] - index_stair - 1 )
 
-        index_polyval = np.polyfit(
+        index_polyval = np.polyfit(            # fraction stair index from intensity value
+            [ stair_before , stair_after      ],
             [ index_stair  , index_stair - 1  ],
-            [ stair_before , stair_after      ]
+            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
 
-        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()
+        final_intensity[index] = true_intensity_value
 
     if verbose:
         print( 'intensity calibration finished' )
@@ -621,9 +627,42 @@ if output == None:
 else:
     if verbose:
         print( 'storing result in ' + output )
-    output_file = pathlib.Path( output )
-    with shelve.open( str( output_file ) ) as output:
-        output[ 'data' ] = mean_data
+    main_hdu  = fits.PrimaryHDU( final_intensity )
+    main_hdu.header[ 'CRVAL1'   ] = wavelengths[0]
+    main_hdu.header[ 'CDELT1'   ] = wavelengths[1] - wavelengths[0]
+    main_hdu.header[ 'CTYPE1'   ] = 'Wavelength'
+    main_hdu.header[ 'CUNIT1'   ] = 'Angstrom'
+    main_hdu.header[ 'OBJNAME'  ] = head[ 'OBJECT' ]
+    # missing from Naroo
+    #main_hdu.header[ 'OBSERVER' ] = head[ '' ]
+    #main_hdu.header[ 'DATE-OBS' ] = head| '' ]
+    #main_hdu.header[ 'EXPTIME'  ] = head[ '' ]
+    #main_hdu.header[ 'TELESCOP' ] = head[ '' ]
+    #main_hdu.header[ 'INSTRUME' ] = head[ '' ]
+    #main_hdu.header[ 'RA'       ] = head[ '' ]
+    #main_hdu.header[ 'DEC'      ] = head[ '' ]
+    #main_hdu.header[ 'EQUINOX'  ] = head[ '' ]
+    #main_hdu.header[ 'MID-HJD'  ] = head[ '' ]
+    #main_hdu.header[ 'DETNAME'  ] = head[ '' ]
+    #main_hdu.header[ 'INSTRUME' ] = head[ '' ]
+    main_hdu.header[ 'RADECSYS' ] = 'FKS'
+    main_hdu.header[ 'OBS-ID'   ] = head[ 'OBS_ID' ]
+
+    # BSS keywords
+    main_hdu.header[ 'BSS_VHEL' ] = 0
+    main_hdu.header[ 'BSS_TELL' ] = 'None'
+    main_hdu.header[ 'BSS_COM'  ] = 'None'
+    # Not understood
+    # main_hdu.header[ 'BSS_ORD'  ] = ?
+    # missing from Naroo
+    #main_hdu.header[ 'BSS_INST' ] = head[ '' ]
+    #main_hdu.header[ 'BSS_SITE' ] = head[ '' ]
+    #main_hdu.header[ 'BSS_LAT'  ] = head[ '' ]
+    #main_hdu.header[ 'BSS_LONG' ] = head[ '' ]
+    #main_hdu.header[ 'BSS_ELEV' ] = head[ '' ]
+
+    hdul = fits.HDUList( [ main_hdu ] )
+    hdul.writeto( output , overwrite = True )
 if verbose:
     print( 'output writing finished' )
     print( '===========================================\