Merge branch 'i4ds' into time-shift-keyword

* i4ds: Imaging pipeline improvements (i4Ds#176) Check online version of ELUT and science energy channel configurations. (i4Ds#175) Fix elut correction factor (i4Ds#178)
grazwegian · Oct 23, 2023 · 23b8738 · 23b8738
2 parents a8cdf39 + d87eb57
commit 23b8738
Show file tree

Hide file tree

Showing 12 changed files with 299 additions and 37 deletions.
diff --git a/.gitignore b/.gitignore
@@ -39,3 +39,4 @@ stx_scenario_*
 /*.fits
 
 .vscode
+
diff --git a/stix/idl/io/stx_check_fits_compatibility.pro b/stix/idl/io/stx_check_fits_compatibility.pro
@@ -30,7 +30,7 @@
 ;-
 function stx_check_fits_compatibility, fits_path
 
-  fits_data = mrdfits( fits_path, 0, primary_header, silent = silent, /unsigned )
+  fits_data = mrdfits( fits_path, 0, primary_header, /silent, /unsigned )
 
   processing_level = sxpar(primary_header,'level')
 

diff --git a/stix/idl/io/stx_read_fits.pro b/stix/idl/io/stx_read_fits.pro
@@ -23,7 +23,7 @@
 ;
 ; :keywords:
 ;
-;    silent : in, type="int", default="0"
+;    silent : in, type="int", default="1"
 ;             If set prevents informational messages being displayed. Passed through to mrdfits.
 ;
 ; :returns:
@@ -45,10 +45,10 @@
 ;-
 function stx_read_fits, fits_path, extension, header, silent = silent, mversion_full = mversion_full
 
-  default, silent, 0
+  default, silent, 1
 
   if  ~keyword_set(mversion_full) then begin
-    ver = mrdfits(/version)
+    ;; ver = mrdfits(/version, /silent)  ; commented out (FSc, 2023-09-26) since not used anywhere
     mversion_full = stx_get_mrd_version()
     mversion = mversion_full.split('\.')
     if ~(fix(mversion[0]) ge 2 and fix(mversion[1]) ge 27) and ~silent then begin

diff --git a/stix/idl/processing/aux_data/stx_create_auxiliary_data.pro b/stix/idl/processing/aux_data/stx_create_auxiliary_data.pro
@@ -83,6 +83,8 @@ if ~nb_within then begin
            (time_middle-time_data[t_before]) * aux_data_str[t_after].Z_SRF) / $
            (time_data[t_after]-time_data[t_before])
   sigma_X = 0.  &  sigma_Y = 0.
+  ; convert to single precision
+  X_SAS = float(X_SAS)  &  Y_SAS = float(Y_SAS)
 
   ; Apparent solar radius (arcsec)
   RSUN = ((time_data[t_after]-time_middle) * aux_data_str[t_before].spice_disc_size + $
@@ -167,6 +169,7 @@ if ~X_SAS.isnan() and ~Y_SAS.isnan() then begin
     print, X_SAS, Y_SAS, format='(" --- found (Y_SRF, -Z_SRF) = ", F7.1,",",F7.1)'
     print, sigma_X, sigma_Y, format='("                 std. dev. = ", F7.1,",",F7.1)'
   endif
+
   ; Correct SAS solution for systematic error measured in 2021
   X_SAS += offset_X  &  Y_SAS += offset_Y
   sas_pointing = [X_SAS, Y_SAS]

diff --git a/stix/idl/processing/imaging/fit_map_2dgauss.pro b/stix/idl/processing/imaging/fit_map_2dgauss.pro
@@ -3,10 +3,12 @@
 ;    fit_map_2dgauss
 ;
 ; PURPOSE:
-;    Fits a 2D-Gaussian on a previously computed map
+;    Fits a 2D-Gaussian on a previously computed map. The centroid positions obtained from Gaussian fitting
+;    is printed on screen, both in helioprojective Cartesian coordinates and in heliographic (Stonyhurst)
+;    longitude and latitude.
 ;
 ; CALLING SEQUENCE:
-;    fit_map_2dgauss, in_map [, /do_plot]
+;    fit_map_2dgauss, in_map, hlon=hlon, hlat=hlat [, /do_plot]
 ;
 ; INPUTS:
 ;    in_map : a map structure, as returned by imaging functions
@@ -15,15 +17,16 @@
 ;   do_plot : if set, overplots contours at 25%, 50% and 75% of the Gaussian peak on the map
 ;             (note that the map has to be already plotted in the active window) 
 ;
-; OUTPUTS:
-;   None. The centroid positions obtained from Gaussian fitting is printed on screen, in both
-;         helioprojective cartesian coordinates and heliographic (Stonyhurst) longitude and latitude
+; OPTIONAL OUTPUTS:
+;   hlon, hlat       : Heliographic Stonyhurst coordinates of source center
+;   gauss_X, gauss_Y : Helioprojective Cartesian coordinates of source center
 ;
 ; MODIFICATION HISTORY:
 ;    2022-01-07: F. Schuller (AIP, Germany): created
+;    2023-02-21, FSc (AIP): added hlon, hlat as optional outputs
 ;
 ;-
-pro fit_map_2dgauss, in_map, do_plot=do_plot
+pro fit_map_2dgauss, in_map, do_plot=do_plot, hlon=hlon, hlat=hlat, gauss_X=gauss_X, gauss_Y=gauss_Y
   ; Dimensions of the map array
   if not is_struct(in_map) then begin
     print,"ERROR: input variable is not a structure."
@@ -42,7 +45,8 @@ pro fit_map_2dgauss, in_map, do_plot=do_plot
   roll_xy, gauss_X/60.,gauss_Y/60., -1.*in_map.roll_angle, p1_x, p1_y
   coord = arcmin2hel(p1_x,p1_y, date=in_map.time, rsun=in_map.rsun, b0=in_map.b0, l0=in_map.l0)
   print,coord[1],coord[0],format='("   --> Heliographic LON = ",F7.2," ; LAT = ",F6.2)'
-
+  hlon = coord[1]  &  hlat = coord[0]
+
   ; Overplot contours at 25%, 50% and 75% of peak
   if keyword_set(do_plot) then begin
     ; arrays of X and Y coordinates

diff --git a/stix/idl/processing/imaging/stx_imaging_pipeline.pro b/stix/idl/processing/imaging/stx_imaging_pipeline.pro
@@ -18,34 +18,63 @@
 ;   imsize       : size (in pixels) of the image to be generated (default: [128, 128])
 ;   pixel        : size (in arcsec) of one pixel in the map (default: [2.,2.])
 ;   x_ptg, y_ptg : if provided, use these values instead of those found in the auxiliary file to correct for pointing
+;   force_sas    : if set, uses SAS solution even if it's far off SolO's pointing
+;   no_sas       : if set, don't use SAS solution but rely on SolO's pointing
+;   subc_labels  : list of sub-collimators to be used in imaging algorithm
 ;
 ; OPTIONAL KEYWORDS:
 ;   force_sas    : if set, use SAS pointing solution even if very different from s/c pointing
 ;   no_sas       : if set, bypass SAS solution and use spacecraft pointing (corrected for systematics) instead
 ;   no_small     : if set, don't use small pixels data to generate the map
-;
+;   method       : select imaging algorithm; should be one of: "MEM" [default], "EM", or "clean"
+;   w_clean      : for clean method, choose between uniform weighting (w_clean=1, default) and natural weighting (w_clean=0)
+;   clean_beam_width : for clean method, size of the beam to convolve the clean components with
+;   set_clean_boxes  : should the user define the clean boxes interactively? (default: NO)
+;   
 ; OUTPUTS:
 ;   Returns a map object that can be displayed with plot_map
+;   
+; OPTIONAL OUTPUT:
+;   path_sci_file : contains the full path to the L1 SCI data file used as input
 ;
 ; EXAMPLES:
-;   mem_ge_map = stx_imaging_pipeline('2109230031', ['23-Sep-2021 15:20:30', '23-Sep-2021 15:22:30'], [18,28], [650.,-650.])
-;   map_1 = stx_imaging_pipeline('2110090002', ['2021-10-09T06:29:50','2021-10-09T06:32:30'], [18,28], [20., 420.])
-;   map_2 = stx_imaging_pipeline('2110090002', ['2021-10-09T06:29:50','2021-10-09T06:32:30'], [4,8], [20., 420.])
+;   mem_ge_map = stx_imaging_pipeline('2109230031', ['23-Sep-2021 15:20:30', '23-Sep-2021 15:22:30'], [18,28])
+;   map_1 = stx_imaging_pipeline('2110090002', ['2021-10-09T06:29:50','2021-10-09T06:32:30'], [18,28])
+;   map_2 = stx_imaging_pipeline('2110090002', ['2021-10-09T06:29:50','2021-10-09T06:32:30'], [4,8])
+;; Example with user-provided pointing correction:
+;   map1 = stx_imaging_pipeline('2204020888', ['2022-04-02T13:18:10','2022-04-02T13:20:40'], [25,50], $
+;                               x_ptg=-1901.0, y_ptg=871.3)
+;; User-defined map center and dimensions:
+;   map2 = stx_imaging_pipeline('2204020888', ['2022-04-02T13:22:30','2022-04-02T13:26:50'], [25,50], $
+;                               x_ptg=-1900.0, y_ptg=871.2, xy_flare=[-1900.,550.], imsize=[261,221], pixel=[2.5,2.5])
+;; Using only a sub-set of collimators:
+;   low_res = ['10a','10b','10c','9a','9b','9c','8a','8b','8c','7b','7c','6c']
+;   map_th = stx_imaging_pipeline('2204020888', ['2022-04-02T13:37','2022-04-02T13:44'], [4,8], $
+;                                 x_ptg=-1897., y_ptg=870.3, xy_flare=[-2000.,600.], imsize=[201,201], pixel=[3.,3.], subc_labels = low_res)
 ;
 ; MODIFICATION HISTORY:
 ;    2022-05-19: F. Schuller (AIP, Germany): created
 ;    2022-08-30, FSc: use stx_estimate_location to find source position if not given
 ;    2022-09-09, FSc: added optional argument bkg_uid
 ;    2022-10-06, FSc: adapted to recent changes in other procedures
+;    2022-11-16, FSc: added optional argument subc_labels
 ;    2023-02-24, FSc: added optional keyword no_small
+;    2023-09-06, FSc: added optional keyword method
+;    2023-10-16, FSc: added optional keywords clean_beam_width and set_clean_boxes
 ;
 ;-
-function stx_imaging_pipeline, stix_uid, time_range, energy_range, bkg_uid=bkg_uid, $
-                               xy_flare=xy_flare, imsize=imsize, pixel=pixel, x_ptg=x_ptg, y_ptg=y_ptg, $
-                               force_sas=force_sas, no_sas=no_sas, no_small=no_small
+function stx_imaging_pipeline, stix_uid, time_range, energy_range, bkg_uid=bkg_uid, xy_flare=xy_flare, $
+                               imsize=imsize, pixel=pixel, x_ptg=x_ptg, y_ptg=y_ptg, force_sas=force_sas, no_sas=no_sas, $
+                               subc_labels=subc_labels, no_small=no_small, method=method, $
+                               w_clean=w_clean, clean_beam_width=clean_beam_width, set_clean_boxes=set_clean_boxes, $
+                               path_sci_file=path_sci_file
+
   if n_params() lt 3 then begin
     print, "STX_IMAGING_PIPELINE"
-    print, "Syntax: result = stx_imaging_pipeline(stix_uid, time_range, energy_range [, xy_flare=xy_flare, imsize=imsize, pixel=pixel, x_ptg=x_ptg, y_ptg=y_ptg])"
+    print, "Syntax: result = stx_imaging_pipeline(stix_uid, time_range, energy_range [, bkg_uid=bkg_uid, xy_flare=xy_flare, $"
+    print, "                 imsize=imsize, pixel=pixel, x_ptg=x_ptg, y_ptg=y_ptg, force_sas=force_sas, no_sas=no_sas, $"
+    print, "                 subc_labels=subc_labels, no_small=no_small, method=method, w_clean=w_clean, $
+    print, "                 clean_beam_width=clean_beam_width, set_clean_boxes=set_clean_boxes, path_sci_file=path_sci_file])"
     return, 0
   endif
 
@@ -54,8 +83,24 @@ function stx_imaging_pipeline, stix_uid, time_range, energy_range, bkg_uid=bkg_u
 ;   l1a_data_folder = '/store/data/STIX/L1A_FITS/L1/'
   l1a_data_folder = '/store/data/STIX/L1_FITS_SCI/'
 
+
+  ; sub-collimator labels
+  default, subc_labels, ['10a','10b','10c','9a','9b','9c','8a','8b','8c','7a','7b','7c','6a','6b','6c','5a','5b','5c','4a','4b','4c','3a','3b','3c']
+  subc_index = stx_label2ind(subc_labels)
+
   default, imsize, [128, 128]
   default, pixel,  [2.,2.]
+
+  ; Imaging algorithm to be used: make sure that the method is implemented
+  default, method, "MEM"
+  method = strupcase(method)
+  known_methods = ["MEM", "EM", "CLEAN"]
+  tst = where(known_methods eq method, i_tst)
+  if ~i_tst then begin
+    print, method, format='("Method ",A," not known. Please use one of the following:")'
+    print, known_methods
+    return, 0
+  endif
 
   ;;;;
 
@@ -93,7 +138,8 @@ function stx_imaging_pipeline, stix_uid, time_range, energy_range, bkg_uid=bkg_u
   ; If not given, try to estimate the location of the source from the data
   !p.background=0
   if not keyword_set(xy_flare) then begin
-    stx_estimate_flare_location, path_sci_file, time_range, aux_data, flare_loc=xy_flare, energy_range=energy_range
+    stx_estimate_flare_location, path_sci_file, time_range, aux_data, $
+                                 flare_loc=xy_flare, energy_range=energy_range, subc_index=subc_index
     print, xy_flare, format='(" *** INFO: Estimated flare location = (",F7.1,", ",F7.1,") arcsec")'
     print, xy_flare / aux_data.rsun, format='(" ... in units of solar radius = (",F6.3,", ",F6.3,")")'
 
@@ -104,14 +150,40 @@ function stx_imaging_pipeline, stix_uid, time_range, energy_range, bkg_uid=bkg_u
   mapcenter = stx_hpc2stx_coord(xy_flare, aux_data)
   flare_loc  = mapcenter
 
-  ; Compute calibrated visibilities
-  if keyword_set(no_small) then $
-     vis = stx_construct_calibrated_visibility(path_sci_file, time_range, energy_range, mapcenter, $
-                                               path_bkg_file=path_bkg_file, xy_flare=flare_loc, /no_small, sumcase='TOP+BOT') $
-     else vis = stx_construct_calibrated_visibility(path_sci_file, time_range, energy_range, mapcenter, $
-                                               path_bkg_file=path_bkg_file, xy_flare=flare_loc)
+  ; Next, we call the functions that take the data as input and generate the emission map. The functions
+  ; to be called depend on the imaging algorithm.
+
+  if method eq "EM" then begin
+    pixel_data_summed = stx_construct_pixel_data_summed(path_sci_file, time_range, energy_range, $
+                                                        path_bkg_file=path_bkg_file, xy_flare=xy_flare)
+
+    out_map = stx_em(pixel_data_summed, aux_data, imsize=imsize, pixel=pixel,mapcenter=mapcenter)
+
+  endif else begin
+    ; Compute calibrated visibilities
+    if keyword_set(no_small) then $
+      vis = stx_construct_calibrated_visibility(path_sci_file, time_range, energy_range, mapcenter, subc_index=subc_index, $
+      path_bkg_file=path_bkg_file, xy_flare=flare_loc, /no_small, sumcase='TOP+BOT') $
+    else vis = stx_construct_calibrated_visibility(path_sci_file, time_range, energy_range, mapcenter, subc_index=subc_index, $
+      path_bkg_file=path_bkg_file, xy_flare=flare_loc)
+
+    case method of
+      "MEM": out_map = stx_mem_ge(vis,imsize,pixel,aux_data,total_flux=max(abs(vis.obsvis)), /silent)
+      "CLEAN": begin
+        default, w_clean, 1  ; 1 = uniform weighting, 0 = natural weighting
+        niter  = 100   ; Number of iterations
+        gain   = 0.1   ; Gain used in each clean iteration
+        nmap   = 10    ; Plot clean components and cleaned map every 10 iterations
+        if not keyword_set(clean_beam_width) then clean_beam_width = 14. ; clean components are convolved with this beam
+        if not keyword_set(set_clean_boxes) then set_clean_boxes = 0
+        clean_map=stx_vis_clean(vis, aux_data, niter=niter, image_dim=imsize[0], PIXEL=pixel[0], $
+                                uniform_weighting = w_clean, gain=gain, nmap=nmap, $
+                                /plot, set_clean_boxes = set_clean_boxes, beam_width=clean_beam_width)
+
+        out_map = clean_map[0]   ; contains the CLEAN map
+      end
+    endcase
+  endelse
 
-  ; Finally, generate the map using MEM_GE
-  out_map = stx_mem_ge(vis,imsize,pixel,aux_data,total_flux=max(abs(vis.obsvis)), /silent)
   return, out_map
 end
diff --git a/stix/idl/processing/imaging/stx_vis_clean.pro b/stix/idl/processing/imaging/stx_vis_clean.pro
@@ -52,6 +52,7 @@
 ;     input visibility bag
 ;   - spatial_frequency_weight - an array with the same number of elements as the input vis bag computed
 ;     by the user's preference.
+;   - set_clean_boxes  If set, let the user define a clean box or polygon
 ;   - box_map       map that is displayed for selecting the clean box
 ;
 ;
@@ -363,7 +364,7 @@ function stx_vis_clean, vis, aux_data, niter = niter, image_dim = image_dim_in,
 
   if keyword_set(plot) then begin
     window,6,xsize=5*this_disp,ysize=2*this_disp
-    cleanplot
+    cleanplot, /silent
     !p.multi=[0,3,1]
     chs2=2.
     plot_map,out0,charsize=chs2

diff --git a/stix/idl/processing/pixel_data/stx_plot_selected_time_range.pro b/stix/idl/processing/pixel_data/stx_plot_selected_time_range.pro
@@ -52,6 +52,10 @@
 pro stx_plot_selected_time_range, tim_axis, energy_ind, time_ind, counts, live_time_bins, subc_index, sumcase, energy_range, $
                                   time_range, counts_bkg=counts_bkg, live_time_bkg=live_time_bkg
 
+  ; re-initialise UTBASE to handle UTC times properly when calling utplot
+  common utcommon
+  utbase = 0
+
 ;;********** Exclude first and last energy bin
 
 lightcurve = counts[1:30,*,*,*]
@@ -169,7 +173,7 @@ wset,2
 chsize=1.2
 chsize_leg=1.8
 clearplot
-utplot,tim_axis,lightcurve, psym=10,ytitle='STIX count rate [s!U-1!N cm!U-2!N keV!U-1!N]', charsize=chsize,$
+utplot,tim_axis,lightcurve, psym=10, /xs, ytitle='STIX count rate [s!U-1!N cm!U-2!N keV!U-1!N]', charsize=chsize,$
        title = this_date + ' ' + this_start_time + '-' + this_end_time + ' UT, ' + $
        trim(energy_range[0],'(f12.1)') + '-' + trim(energy_range[1],'(f12.1)') + ' keV'
 outplot,tim_axis[time_ind],lightcurve[time_ind],psym=10,color=122