Some Improvement in aspect processing

stix/idl/io/stx_read_aux_fits.pro

@@ -68,6 +68,8 @@ pro stx_read_aux_fits, fits_path, aux_data=aux_data, primary_header = primary_he
   aux_data.spice_disc_size            = data.spice_disc_size
   aux_data.y_srf                      = data.y_srf
   aux_data.z_srf                      = data.z_srf
+  if tag_exist(data, "sas_ok") then $
+    aux_data.sas_ok                   = data.sas_ok
   aux_data.solo_loc_carrington_lonlat = data.solo_loc_carrington_lonlat
   aux_data.solo_loc_carrington_dist   = data.solo_loc_carrington_dist
   aux_data.solo_loc_heeq_zxy          = data.solo_loc_heeq_zxy

stix/idl/processing/aspect/SAS_test.pro

@@ -3,67 +3,64 @@
 ; 
 ; F. Schuller (AIP Potsdam, Germany)
 ;    Created:     2021-08-06
-;    Last modif.: 2022-01-31
+;    Last modif.: 2023-09-19
 ; 
 
 ;;;;
 ; First, let's define paths to input directories and parameter files
 ;
 data_dir = '/store/data/STIX/L1_FITS_HK/'
-param_dir = '/home/fschuller/Documents/IDL/STIX-GSW/stix/idl/processing/aspect/SAS_param/'
-calib_file = param_dir + 'SAS_calib_20211005.sav'
+param_dir = getenv('SAS_PARAM_DIR')
+def_calibfile = getenv('SAS_CALIBFILE')
+calib_file = param_dir +def_calibfile
 aperfile = param_dir + 'apcoord_FM_circ.sav'
-simu_data_file = param_dir + 'SAS_simu.sav'   ; this can be a link to the actual file containing simulated data
-
-;;;
-; Read SPICE kernel: needed here to compute SPICE_DISC_SIZE
-; NB: This won't be necessary once the data reading and formatting is done within STIXcore
-spice_dir = getenv('SPICE_DIR')
-spice_kernel = getenv('SPICE_KERNEL')
-cspice_kclear
-add_sunspice_mission, 'solo'
-load_sunspice_solo
-cspice_furnsh, spice_dir + spice_kernel
+; simu_data_file = param_dir + 'SAS_simu.sav'   ; this can be a link to the actual file containing simulated data
+simu_data_file = param_dir + getenv('SAS_DATA_SIMU')
 
 ; Read data from L1 FITS file and store them in an array of stx_aspect_dto structures
 ; NB: This won't be necessary once the data reading and formatting is done within STIXcore
 ;
 print,"Reading L1 data file..."
-in_file = "solo_L1_stix-hk-maxi_20210401_V01.fits"
+
+; one_day = '20230330'    ; case where the HK file contains too many rows (duplicate entries, many with duration close to 0)
+; one_day = '20220509'    ; solar distance gets beyond 0.75 AU at the end of that day
+; one_day = '20230321'    ; pointing mostly at Sun centre, with a flat-field calib. from 19:00 to 20:00
+one_day = '20230329'    ; includes pointing at pole, and other off-centre pointings, and too many rows in HK file
+
+in_file = "solo_L1_stix-hk-maxi_"+one_day+"_V01.fits"
 data = prepare_aspect_data(data_dir + in_file, quiet=0)
-show_info, data
+stx_show_sas_info, data
 
 print,"Calibrating data..."
 ; First, substract dark currents and applies relative gains
-calib_sas_data, data, calib_file
+stx_calib_sas_data, data, calib_file
+; copy result in a new object
+data_calib = data
+; Added 2023-09-18: remove data points with some error detected during calibration
+stx_remove_bad_sas_data, data_calib
+
 ; Now automatically compute global calibration correction factor and applies it
 ; Note: this takes a bit of time
-auto_scale_sas_data, data, simu_data_file, aperfile
+stx_auto_scale_sas_data, data_calib, simu_data_file, aperfile
 
 print,"Plotting the signals..."
-show_info, data
-plot4sig, data
+stx_show_sas_info, data_calib
+plot4sig, data_calib
+
+; apply same calibration correction factor to all data (including data points that were removed due to errors)
+cal_corr_factor = data_calib[0].calib
+data.CHA_DIODE0 *= cal_corr_factor
+data.CHA_DIODE1 *= cal_corr_factor
+data.CHB_DIODE0 *= cal_corr_factor
+data.CHB_DIODE1 *= cal_corr_factor
 
 print,"Computing aspect solution..."
-derive_aspect_solution, data, simu_data_file, interpol_r=1, interpol_xy=1
+stx_derive_aspect_solution, data, simu_data_file, interpol_r=1, interpol_xy=1
 !p.multi = [0,1,2]
-utplot, data.TIME, data.y_srf, /xs, /ynoz, ytit='!6Y!dSRF !n [arcsec]',chars=1.4
-utplot, data.TIME, data.z_srf, /xs, /ynoz, ytit='!6Z!dSRF !n [arcsec]',chars=1.4
-;;;;;;;;;;;;;;;;;
-
+; Only plot solution with no error message
+good = where(data.sas_ok eq 1)
+utplot, data[good].TIME, data[good].y_srf, /xs, /ynoz, ytit='!6Y!dSRF !n [arcsec]',chars=1.4,/psym
+utplot, data[good].TIME, data[good].z_srf, /xs, /ynoz, ytit='!6Z!dSRF !n [arcsec]',chars=1.4,/psym
+;;;;;;;;;;;;;;;;
 
-;;;;
-; Alternatively, call "process_SAS_data" to do at once:
-;  - calib_sas_data
-;  - auto_scale_sas_data
-;  - derive_aspect_solution
-;
-print,"End-to-end processing test..."
-data_dir = '/store/data/STIX/L1_FITS_HK/'
-in_file = "solo_L1_stix-hk-maxi_20210401_V01.fits"
-data = prepare_aspect_data(data_dir + in_file, quiet=0)
-cal_factor = 1.10   ; starting value, roughly OK for 2021 data [NB: this parameter is fully optional]
-process_SAS_data, data, calib_file, simu_data_file, aperfile, cal_factor=cal_factor
-; Now the results are stored in data.y_srf and data.z_srf
-;
-;;;;
+end

stix/idl/processing/aspect/prepare_aspect_data.pro

@@ -76,7 +76,11 @@ function prepare_aspect_data, infile, quiet=quiet
     a =  {stx_aspect_dto, $
           cha_diode0: tbl[i].hk_asp_photoa0_v, cha_diode1: tbl[i].hk_asp_photoa1_v, $
           chb_diode0: tbl[i].hk_asp_photob0_v, chb_diode1: tbl[i].hk_asp_photob1_v, $
-          time: utc_str[i], duration : duration[i], spice_disc_size : r_sol[i], y_srf : 0.0, z_srf : 0.0, calib : 0.0,  error : ""}
+          time: utc_str[i], $
+           scet_time_c: 0LL, scet_time_f: 0LL, $
+          duration : duration[i], spice_disc_size : r_sol[i], $
+          y_srf : 0.0, z_srf : 0.0, calib : 0.0,  sas_ok : 1, error : "", $
+          control_index: 0LL, parentfits: 0LL}
     if i eq 0 then result = [a] else result = [result, a]
   endfor
 

stix/idl/processing/aspect/process_sas_data.pro

@@ -24,10 +24,26 @@
 ;
 ;-
 
-pro process_SAS_data, data, calibfile, simu_data_file, aperfile, cal_factor=cal_factor
+pro process_SAS_data, data, calib_file, simu_data_file, aperfile, cal_factor=cal_factor
   default, cal_factor, 1.
 
-  calib_sas_data, data, calibfile, factor=cal_factor
-  auto_scale_sas_data, data, simu_data_file, aperfile
-  derive_aspect_solution, data, simu_data_file, /interpol_r, /interpol_xy
+  ; First, substract dark currents and applies relative gains
+  stx_calib_sas_data, data, calib_file, factor=cal_factor
+  ; copy result in a new object
+  data_calib = data
+  ; remove data points with some error detected during calibration
+  stx_remove_bad_sas_data, data_calib
+
+  ; Now automatically compute global calibration correction factor and applies it
+  stx_auto_scale_sas_data, data_calib, simu_data_file, aperfile
+
+  ; apply same calibration correction factor to all data (including data points that were removed due to errors)
+  cal_corr_factor = data_calib[0].calib
+  data.CHA_DIODE0 *= cal_corr_factor
+  data.CHA_DIODE1 *= cal_corr_factor
+  data.CHB_DIODE0 *= cal_corr_factor
+  data.CHB_DIODE1 *= cal_corr_factor
+
+  ; Compute aspect solution
+  stx_derive_aspect_solution, data, simu_data_file, /interpol_r, /interpol_xy
 end

stix/idl/processing/aspect/stx_calib_sas_data.pro

@@ -4,10 +4,10 @@
 ;
 ; Category    : analysis
 ;
-; Syntax      : calib_sas_data, data , calibfile  [, factor=factor ]
+; Syntax      : stx_calib_sas_data, data , calibfile  [, factor=factor ]
 ;
 ; Inputs      :
-;   data      = a structure as returned by read_hk_data
+;   data      = a structure as returned by prepare_aspect_data
 ;   calibfile = name of the file with calibration parameters (gains and bias values), including full absolute path
 ;
 ; Output      : The input structure is modified
@@ -27,7 +27,7 @@ pro stx_calib_sas_data, data, calibfile, factor=factor
 
   default, factor, 1.
 
-  if n_params() lt 2 then message,' SYNTAX: calib_sas_data, data, calibfile [, factor=factor]'
+  if n_params() lt 2 then message,' SYNTAX: stx_calib_sas_data, data, calibfile [, factor=factor]'
 
   if not is_struct(data) then message," ERROR: input variable is not a structure."
 
@@ -49,8 +49,8 @@ pro stx_calib_sas_data, data, calibfile, factor=factor
   nb_rows = n_elements(data)
   if nb_ok lt nb_rows then begin
     print,nb_rows - nb_ok,format='(" CALIB_SAS_DATA Warning: Found",I5," entries with wrong duration.")'
-    ind_bad = where(abs(data.duration-64.) ge 0.1,nb_bad)
-    for i=0,nb_bad-1 do data[ind_bad[i]].ERROR = "SAS_DATA_WRONG_DURATION"
+    ind_bad = where(abs(data.duration-64.) ge 0.5,nb_bad)
+    for i=0,nb_bad-1 do data[ind_bad[i]].ERROR = "SAS_WRONG_DUR"
   endif
   data[ind_ok].CHA_DIODE0 = (data[ind_ok].CHA_DIODE0 / 16. - V_base) / R_m
   data[ind_ok].CHA_DIODE1 = (data[ind_ok].CHA_DIODE1 / 16. - V_base) / R_m