Add & Fix final format
This commit is contained in:
parent
5e8e21b66b
commit
32403fdce4
1 changed files with 62 additions and 23 deletions
85
spectrum.py
85
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_stair , index_stair - 1 ],
|
||||
[ stair_before , stair_after ]
|
||||
)
|
||||
|
||||
stair_polyval = np.polyfit(
|
||||
index_polyval = np.polyfit( # fraction stair index from intensity value
|
||||
[ stair_before , stair_after ],
|
||||
[ 0 , step ]
|
||||
[ index_stair , index_stair - 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 = 0
|
||||
|
||||
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( '===========================================\
|
||||
|
|
|
|||
Loading…
Reference in a new issue