;+
; PROJECT:
;   HESSI
;
; NAME:
;   hsi_vis_normalize
;
; PURPOSE:
;   Corrects visibility amplitudes to compensate for presumed relative detector calibration errors
;   Visibilities are in the form of an array of 'standard' RHESSI visibility structures.
;
; CATEGORY:
;   \hessi\idl\imaging\visibilities\gh
;
; CALLING SEQUENCE:
;   visout = hsi_vis_normalize(visin, MAX_CORR=max_corr, CORR_FACTORS=corr_factors)
;
; CALLS:
;   hsi_conditional_average
;
; INPUTS:
;   visin = a bag of visibilities to be normalized
;
; OPTIONAL INPUTS:
;   none
;
; OUTPUTS:
;   visout = normalized version of input visibilities.

; OPTIONAL OUTPUTS:
;   See keyword, CORR_FACTORS.
;
; KEYWORDS:
;   MAX_CORR        = maximum allowable fractional difference of 1 detector's total flux from the average.  (Default = 0.25)
;	 VALID_DET_INDEX_MASK = 9 element mask (0's and 1's) showing which detectors where included in the average total flux
;   CORR_FACTORS    = 9-element vector (one element per detector) of factors that were used to multiply visibility amplitudes, total flux and sigamp
;	CORR_FACTORS of 0 means that this detector did not have an average positive TOTFLUX
;	QUIET - if set, don't print info message
;
; COMMON BLOCKS:
;   none
;
; METHOD:
;   Method is based on the expectation that the average values of TOTFLUX will be the same for all detectors.
;   Calculates the ratio of the detector-associated values of TOTFLUX to the average and then multiplies these values,
;       the visibility amplitudes and their statistical errors by CORR_FACTORS to make <TOTFLUX> the same for all detectors.
;   Detectors whose initial <TOTFLUX> differs from the average by more than a factor of (1+-MAX_CORR) are
;       excluded from participation in calculating the average, but they are still modified as above for use downstream
;
; SIDE EFFECTS:
;   None -  visin is not modified.
;
; RESTRICTIONS:
;   It is assumed that visin values have been run through hsi_vis_edit() to remove bad values.
;   It is assumed that all visibilities correspond to the same energy- and time-range.
;   Returns with visout=visin and an error message if this is not the case.
;
; MODIFICATION HISTORY:
;   24-Nov-08       Initial version (ghurford@ssl.berkeley.edu);
;   29-Dec-08 gh    Added keywords, MAX_CORR, CORR_FACTORS.
;                   Excluded outlier detectors from calculation of <TOTFLUX>
;                   Added ssw header.
;    8-Jan-09 gh    Add printed warning if any RMC's are excluded from the averaging calculation
;   21-Jan-09 gh    Correct indexing bug.
;   30-Apr-09 richard,kim  Return valid_det_index_mask. Also in visout, only include vis for detectors
;                   that passed max_corr test, and added quiet keyword
;  14-sep-2009 ras, gh wants no detectors removed, just outliers excluded from the average
;  02-Jun-2017 Kim  Added to print statement about inconsistent RMCs
;  14-Dec-2018 Kim  If visin contains only one visibility, just return visin
;  28-Aug-2019, Kim. Changed 'message' to 'print' (because message outputs multiple lines in log file)
;-
;
FUNCTION hsi_vis_normalize, visin, MAX_CORR=max_corr, $
  CORR_FACTORS=corr_factors,$
  VALID_DET_INDEX_MASK = valid_det_index_mask, quiet=quiet
  ;
  DEFAULT, max_corr, 0.25
  corr_factors     = FLTARR(9) + 0.
  avtotflux        = FLTARR(9);
  default, quiet, 0
  ;
  ; Check that all time and energy intervals are the same.
  IF  N_ELEMENTS(UNIQ(visin.trange[0]))       NE 1 OR $
    N_ELEMENTS(UNIQ(visin.trange[1]))       NE 1 OR $
    N_ELEMENTS(UNIQ(visin.erange[0])) NE 1 OR $
    N_ELEMENTS(UNIQ(visin.erange[1])) NE 1 THEN BEGIN
    PRINT, 'WARNING: hsi_vis_normalize aborting due to inconsistent energy- and/or time-ranges.'
    RETURN, visin
  ENDIF
  
  if n_elements(visin) eq 1 then begin
    message, 'Only one visibility - not normalizing.', /info
    return, visin
  endif
  
  ;
  ; Calculate the average totflux value for each detector; then the grand average, excluding outliers
  ;FOR isc = 0,8 DO BEGIN
  ;    iok = WHERE(visin.isc EQ isc, niok)
  ;    IF niok GT 0 THEN avtotflux[isc] = AVERAGE(visin[iok].totflux)                      ; average for each detector
  ;ENDFOR
  h = histogram( visin.isc, min=0, max=8, rev=rr)
  selisc = where( h, nselisc)
  avtotflux = fltarr(9)
  for i = 0, nselisc-1 do begin
    isc = selisc[i]
    avtotflux[isc] = average( visin[ rr[rr[isc]: rr[isc+1]-1]].totflux)
  endfor
  valid_det_index_mask = intarr(9)
  j               = WHERE(avtotflux GT 0, n_ok)                                                 ; Indices of detectors for which there are visibiilities.
  valid_det_index_mask[j] = 1
  grandavtotflux  = hsi_conditional_average(avtotflux[j], max_corr, OK_FLAG=ok_flag)      ; calculates the average of avtotflux, excluding outliers

  not_ok          = WHERE(ok_flag EQ 0, nnotok)                                           ; indices of outlier detectors referenced to j
  if nnotok gt 0 then begin
    valid_det_index_mask[j[not_ok]] = 0 							     							;indices of outliers referenced to whole set
    if not quiet then print, 'HSI_VIS_NORMALIZE: WARNING - RMC ' + arr2str(trim(1+j[not_ok])) + $
      ' total flux was not consistent with the others. Not used in calculating average.'
  endif
  ;
  ; Normalize visibilities and return.
  corr_factors[j] = grandavtotflux / avtotflux[j]
  ;Next line removed by ras, 14-sep-2009, gh wants no detectors removed, just outliers excluded from the average
  ;visout          = visin[ where_arr( visin.isc, where( valid_det_index_mask))]  ; include only the good detectors
  visout = visin
  isc = visout.isc
  visout.obsvis   = visout.obsvis  * corr_factors[isc]
  visout.totflux  = visout.totflux * corr_factors[isc]
  visout.sigamp   = visout.sigamp  * corr_factors[isc]

  RETURN, visout
END