;+------------------------------------------------------------------------
; NAME: semiMinMax
; PURPOSE: For array y, return min and max values within 3 standard
;          deviations from the mean
; CALLING SEQUENCE:
;       semiMinMax, y, ymin, ymax
; INPUTS:
;       y = array of 4 or more elements
; KEYWORD PARAMETERS:
;       MODIFIED = use modified sigma which uses kurt. & skew. 
;                  results in a more liberal acceptance of data
; OUTPUTS:
;       ymin = minimum value within 3 stdev's of mean
;       ymax = maximum value within 3 stdev's of mean
; AUTHOR: ?
; MODIFICATION HISTORY:
;   10/96 R.Baldwin; replaced stdev w/ moment
;   03/99 R.Baldwin; added MODIFIED keyword
;
;Copyright 1996-2013 United States Government as represented by the 
;Administrator of the National Aeronautics and Space Administration. 
;All Rights Reserved.
;
;------------------------------------------------------------------
;
pro semiMinMax, y , ymin, ymax, MODIFIED=MODIFIED

;TJK modified on 12/18/2000 to put in a check so that ymin and ymax will never
;be adjusted any high/lower than what they start out as.  This modification 
;was requested by Bob, because of the new IMAGE data.

;TJK modified again on 1/10/2001 to check if ymin, ymax are actually defined
;since it looks like plottypes, except the image ones, do not have these set
;before entering this routine.

if (n_elements(ymin) eq 0) then minmax_flag = 0 else minmax_flag = 1

if (minmax_flag) then begin
  orig_min = ymin
  orig_max = ymax
endif

if keyword_set(MODIFIED) then mod_flag=1L else mod_flag=0L

if (n_elements(y) lt 3) then begin
  message, 'Not enough points in array', /info
  if (n_elements(y) lt 2) then begin
    Ymin = -1. & Ymax = -1.
  endif else Ymin = min(y, max=Ymax); 2 points only
  return
endif

Ymin = min(y, max=Ymax)
;TJK added check on 12/18/2000

if (minmax_flag) then begin
  if (ymin lt orig_min) then ymin = orig_min
  if (ymax gt orig_max) then ymax = orig_max
endif

if(Ymin eq Ymax) then return ;Avoid call to MOMENT where Std. Dev.=0

;as = stdev(y, mean) ;RTB replaced w/ moment
res=moment(y,sdev=as)
if(mod_flag) then begin
; Calculate Modified sigma taking into account the skewness and kurtosis 
; of the distribution
 as=res[1]*(1.+res[2]-0.01*(res[3]-3))  
 as=sqrt(abs(as))
endif

mean=res[0]
w = where(abs(y-mean) lt 3.*as, wc) ; use only points within 3 standard
;                                     deviations of the mean
if (wc gt 0) then begin
;print, 'time before for loop = ',systime()
;  Ymin = y(w(0)) & Ymax = y(w(0)) ;initialize
;  for i = 1L, wc-1 do begin
;	if (y(w(i)) lt Ymin) then Ymin =  y(w(i))
;	if (y(w(i)) gt Ymax) then Ymax =  y(w(i))
;  endfor
;print, 'time after for loop = ',systime()

  Ymin = min(y[w], max=Ymax) 

endif ; else begin
; Ymin = min(y, max=Ymax)
;endelse
;print, Ymin, Ymax

;TJK added check on 12/18/2000
if (minmax_flag) then begin
  if (ymin lt orig_min) then ymin = orig_min
  if (ymax gt orig_max) then ymax = orig_max
endif

return
end ; semiMinMax