Implemented error code SAS_SOL_TOO_FAR

i4Ds · Oct 6, 2023 · 00614f6 · 00614f6
1 parent 042e485
commit 00614f6
Show file tree

Hide file tree

Showing 3 changed files with 49 additions and 30 deletions.
diff --git a/stix/idl/processing/aspect/SAS_test.pro b/stix/idl/processing/aspect/SAS_test.pro
@@ -24,6 +24,9 @@ print,"Reading L1 data file..."
 
 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
+
 in_file = "solo_L1_stix-hk-maxi_"+one_day+"_V01.fits"
 data = prepare_aspect_data(data_dir + in_file, quiet=0)
 show_info, data
@@ -54,8 +57,10 @@ data.CHB_DIODE1 *= cal_corr_factor
 print,"Computing aspect solution..."
 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.ERROR eq '')
+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
+;;;;;;;;;;;;;;;;
 
 end
diff --git a/stix/idl/processing/aspect/stx_calib_sas_data.pro b/stix/idl/processing/aspect/stx_calib_sas_data.pro
@@ -50,7 +50,7 @@ pro stx_calib_sas_data, data, calibfile, factor=factor
   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.5,nb_bad)
-    for i=0,nb_bad-1 do data[ind_bad[i]].ERROR = "SAS_DATA_WRONG_DURATION"
+    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

diff --git a/stix/idl/processing/aspect/stx_derive_aspect_solution.pro b/stix/idl/processing/aspect/stx_derive_aspect_solution.pro
@@ -37,6 +37,9 @@ function solve_aspect_one_plane, inputA_B, inputC_D, plane_AB, plane_CD, all_X,
   default, max_iter, 10      ; stop after iterations...
   default, delta_conv, 10.   ; or if successive solutions don't differ by more than 10 mic
 
+  ; Store number of elements along main direction in simulated data
+  nb_X = n_elements(all_X)
+
   ; Test: if inputA_B or inputC_D is not within the range covered by simulated data,
   ; then we cannot derive any solution: set results to NaN
   if inputA_B lt min(plane_AB) or inputC_D lt min(plane_CD) or $
@@ -82,30 +85,33 @@ function solve_aspect_one_plane, inputA_B, inputC_D, plane_AB, plane_CD, all_X,
       if n_iter ge max_iter or sol_diff lt delta_conv then do_more = 0
     endwhile
 
-    if keyword_set(interpol_XY) then begin
-      ; refine solution by interpolating between two nearest points on the grid
-      delta_AB = inputA_B - plane_AB[tmpAB,ind_Y_AB]
-      ; by construction, plane_AB[*,ind_Y] is monotonically decreasing, therefore:
-      if delta_AB lt 0 then begin
-        ind_AB_pos = max([0,tmpAB-1])  &  ind_AB_neg = tmpAB
-      endif else begin
-        ind_AB_pos = tmpAB  &  ind_AB_neg = min([tmpAB+1,1000])
-      endelse
-      delta_pos = inputA_B - plane_AB[ind_AB_pos,ind_Y_AB]
-      delta_neg = plane_AB[ind_AB_neg,ind_Y_AB] - inputA_B
-      x_AB = delta_pos / (delta_pos+delta_neg) * all_X[ind_AB_neg] + delta_neg / (delta_pos+delta_neg) * all_X[ind_AB_pos]
+    if keyword_set(interpol_XY) then $
+      ; refine solution by interpolating between two nearest points on the grid,
+      ; but only if not on the edge of the parameter space:
+      if tmpAB gt 0 and tmpAB lt nb_X-1 and tmpCD gt 0 and tmpCD lt nb_X-1 then begin
+        delta_AB = inputA_B - plane_AB[tmpAB,ind_Y_AB]
+        ; At least near Sun centre, plane_AB[*,ind_Y] is monotonically increasing, 
+        ; therefore delta_AB is decreasing. Thus:
+        if delta_AB lt 0 then begin
+          ind_AB_pos = tmpAB-1  &  ind_AB_neg = tmpAB
+        endif else begin
+          ind_AB_pos = tmpAB  &  ind_AB_neg = tmpAB+1
+        endelse
+        delta_pos = inputA_B - plane_AB[ind_AB_pos,ind_Y_AB]
+        delta_neg = plane_AB[ind_AB_neg,ind_Y_AB] - inputA_B
+        x_AB = delta_pos / (delta_pos+delta_neg) * all_X[ind_AB_neg] + delta_neg / (delta_pos+delta_neg) * all_X[ind_AB_pos]
 
-      ; Same game for x_CD
-      delta_CD = inputC_D - plane_CD[tmpCD,ind_Y_CD]
-      if delta_CD lt 0 then begin
-        ind_CD_pos = tmpCD-1  &  ind_CD_neg = tmpCD
-      endif else begin
-        ind_CD_pos = tmpCD  &  ind_CD_neg = tmpCD+1
-      endelse
-      delta_pos = inputC_D - plane_CD[ind_CD_pos,ind_Y_CD]
-      delta_neg = plane_CD[ind_CD_neg,ind_Y_CD] - inputC_D
-      x_CD = delta_pos / (delta_pos+delta_neg) * all_X[ind_CD_neg] + delta_neg / (delta_pos+delta_neg) * all_X[ind_CD_pos]
-    endif
+        ; Same game for x_CD
+        delta_CD = inputC_D - plane_CD[tmpCD,ind_Y_CD]
+        if delta_CD lt 0 then begin
+          ind_CD_pos = tmpCD-1  &  ind_CD_neg = tmpCD
+        endif else begin
+          ind_CD_pos = tmpCD  &  ind_CD_neg = tmpCD+1
+        endelse
+        delta_pos = inputC_D - plane_CD[ind_CD_pos,ind_Y_CD]
+        delta_neg = plane_CD[ind_CD_neg,ind_Y_CD] - inputC_D
+        x_CD = delta_pos / (delta_pos+delta_neg) * all_X[ind_CD_neg] + delta_neg / (delta_pos+delta_neg) * all_X[ind_CD_pos]
+      endif
   endelse
 
   result = {x_AB:x_AB, x_CD:x_CD}
@@ -179,12 +185,20 @@ pro stx_derive_aspect_solution, data, simu_data_file, interpol_r=interpol_r, int
       ; convert to SAS frame
       x_sas[i] = -1.*(x_AB - x_CD) / sqrt(2.) * 1.e-6
       y_sas[i] = -1.*(x_AB + x_CD) / sqrt(2.) * 1.e-6
-      if (~finite(x_AB)  or ~finite(x_CD)) then data[i].ERROR = 'NO_ASPECT_SOL'
+
+      ; Test whether the solution can be used:
+      if (~finite(x_AB) or ~finite(x_CD)) then data[i].ERROR = 'NO_ASPECT_SOL' $
+        else begin
+          ; also not good if too far off the solar limb:
+          dist_center = norm([x_AB,x_CD]) * 1.e-6
+          if dist_center gt 1.1 * rsol[i] then data[i].ERROR = 'SAS_SOL_TOO_FAR'
+        endelse
     endif
   endfor
 
   ; Store results as arcsec in SRF in the data structure
-  data.y_srf = y_sas * 0.375e6
-  data.z_srf = x_sas * 0.375e6
+  linear_to_asec = (1./foclen) * 180./!pi * 3600.  ; conversion factor from m to arcsec
+  data.y_srf = y_sas * linear_to_asec
+  data.z_srf = x_sas * linear_to_asec
 
 end