C+
C NAME:
C	smei_skyd_combine
C PURPOSE:
C	Combines in into a single value.
C CATEGORY:
C	smei/camera/for/sky
C CALLING SEQUENCE:
	integer function smei_skyd_combine(iType,bKeepGlitch,iSilent,level,node_id,count,val,flag,pos,rgood,rbad,badpix)
	!integer function smei_skyd_combine(iType,bKeepGlitch,iSilent,count,val,flag,rgood,rbad,badpix)
C INPUTS:
C	iType			integer		defines type of final response
C						0: Mean
C						1: Median
C						2: Mean/median (ratio)
C						3: Mean-median (difference)
C						4: Standard deviation
C						5: Median/std (ratio)
C	bKeepGlitch		logical 	suppress glitch removal
C	iSilent			integer		higher values suppress more informational messages
C	count			integer		number of votes
C	level			integer
C	node_id			integer
C	val(*)			real		response for all votes
C	flag(*)			integer*1	bitwise flags (e.g. the 'near star' flag)
C	pos(*)			integer		linear array index for the originating SMEI frame
C						for each vote
C	badpix(*)		integer		'bad pixel' map		
C OUTPUTS:
C	smei_skyd_combine	integer		number of good votes from which final response
C						is calculated
C	rgood			real		final response for selected node,
C						calculated from all good votes
C	rbad			real		fraction of total adus in rejected votes: 1:nlow-1 where
C						too low; nhigh+1,count were too high
C	badpix(*)		integer		updated 'bad pixel' map
C INCLUDE:
	include		'filparts.h'
	include		'smei_skyd_dim.h'
C CALLS:
C	ArrR4Total, BadR4, Say, cInt2Str, Str2Str, Int2Str
C	IndexR4, ArrR4Mean, ArrR4Median, ArrR4Stdev, smei_skyd_node2skyloc
C PROCEDURE:
C MODIFICATION HISTORY:
C	JAN-2006, Paul Hick (UCSD/CASS)
C	JUN-2011, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
C	    Added arguments level and node_id (used in constructing diagnostic messages)
C	    Made diagnostic messages more intelligible. All votes are now
C	    printed if iSilent <=1.
C-
	    integer	iType
	    logical	bKeepGlitch
	    integer	iSilent
	    integer	count
	    real	val (*)
	    integer*1	flag(*)
	    integer	pos (*)
	    real	rgood
	    real	rbad
	    integer	badpix(*)

	character	cSay*12			/'skyd_combine'/
	character	cStr*(FIL__LENGTH)

	integer		Str2Str
	integer		Flt2Str

	integer		nvote_min		/1/

	integer		indx	 (SKYD__VOTES_PER_NODE_MAX)

	real		val_sort (SKYD__VOTES_PER_NODE_MAX)
	integer		flag_sort(SKYD__VOTES_PER_NODE_MAX)
	integer		pos_sort (SKYD__VOTES_PER_NODE_MAX)

	character	cInt2Str*14

	badval = BadR4()

	if (count .gt. SKYD__VOTES_PER_NODE_MAX)
     &	    call Say('skyd_combine','E','SKYD__VOTES_PER_NODE_MAX='//cInt2Str(SKYD__VOTES_PER_NODE_MAX),
     &		 'is too small#must be at least '//cInt2Str(count))

	call smei_skyd_node2skyloc(level,node_id,kp,ix,jy)

	istr = 0
	istr = istr+Str2Str('('	, cStr(istr+1:))
	istr = istr+Int2Str(kp	, cStr(istr+1:))
	istr = istr+Str2Str(','	, cStr(istr+1:))
	istr = istr+Int2Str(ix	, cStr(istr+1:))
	istr = istr+Str2Str(','	, cStr(istr+1:))
	istr = istr+Int2Str(jy	, cStr(istr+1:))
	istr = istr+Str2Str(')'	, cStr(istr+1:))+1
	
	! Sort arrays

	call IndexR4(1,count,1,count,val,indx)
	do i=1,count
	    j = indx(i)
	    val_sort (i) = val (j)
	    flag_sort(i) = flag(j)
	    pos_sort (i) = pos (j)
	end do

	! Count number of votes with the near_star flag set.

	near_star = 0
	do i=1,count
	    if (btest(flag_sort(i),0)) near_star = near_star+1
	end do

	if (near_star .ne. 0) then
	    i = istr
	    i = i+Int2Str(near_star	, cStr(i+1:))
	    i = i+Str2Str('/'		, cStr(i+1:))
	    i = i+Int2Str(count		, cStr(i+1:))+1
	    i = i+Str2Str('"near star"'	, cStr(i+1:))
	    if (near_star .lt. 0.5*count) then
		!if (iSilent .le. 0) call Say(cSay,'W','marked', cStr(:i)//', IGNORED !')
		near_star = 0
	    else
		!if (iSilent .le. 0) call Say(cSay,'W','marked', cStr)
		near_star = 1
	    end if
	end if

	if (count .lt. nvote_min) then		! Not enough votes: reject node

	    rmedian = badval

	    nlow    = 1
	    nhigh   = 0				! Effectively discards all votes
	    nvote   = 0

	    nlow_old  = 1
	    nhigh_old = count
	    nvote_old = count

	else if (bKeepGlitch) then		! Cmd line switch -keepglitch set: accept node

	    imedian = min(max(count/2,2),count)
	    rmedian = val_sort(imedian)

	    nlow    = 1
	    nhigh   = count
	    nvote   = count

	    nlow_old  = 1
	    nhigh_old = count
	    nvote_old = count

	else if (near_star .eq. 1) then		! Close to bright star

	    imedian = min(max(count/2,2),count)
	    rmedian = val_sort(imedian)
	    !rstddev = ArrR4Stdev(0,imedian,val_sort,rmedian,i)

	    nlow    = 1
	    nhigh   = count
	    nvote   = count

	    nlow_old  = nlow
	    nhigh_old = nhigh
	    nvote_old = nvote

	    ! Note that this test gest ugly if the median is negative

	    dtest = 3.3*rmedian
	    !dtest = rmedian+3.0*rstddev
	    !print *, '2', rmedian, dtest, dtest1, count
	    do while (nhigh .gt. imedian .and. val_sort(nhigh) .gt. dtest)
		nhigh = nhigh-1
	    end do

	    if (nhigh .lt. nlow) then
		print *, (val_sort(i),i=1,count)
		print *, nlow_old, nhigh_old, nvote_old, imedian, rmedian
		print *, nlow, nhigh
		stop 'skyd_combine, oops 1'
	    end if

	    dtest = min(0.1*rmedian,rmedian-50.0)
	    !dtest = rmedian-3.0*rstddev
	    !print *, '2', rmedian, dtest, dtest1, count
	    do while (nlow .lt. imedian .and. val_sort(nlow) .lt. dtest .and. btest(flag_sort(nlow),1))
		nlow = nlow+1
	    end do

	    if (nhigh .lt. nlow) then
		print *, (val_sort(i),i=1,count)
		print *, nlow_old, nhigh_old, nvote_old, imedian, rmedian
		print *, nlow, nhigh
		stop 'skyd_combine, oops 2'
	    end if

	    nvote = nhigh-nlow+1

	    if (nvote .ne. nvote_old) rmedian = ArrR4Median(nvote,val_sort(nlow),0.0,i,val_sort(nlow))

	else

	    imedian = min(max(count/2,2),count)
	    rmedian = val_sort(imedian)
	    !rstddev = ArrR4Stdev(0,imedian,val_sort,rmedian,i)

	    nlow    = 1
	    nhigh   = count
	    nvote   = count

	    nlow_old  = nlow
	    nhigh_old = nhigh
	    nvote_old = nvote

	    dtest = max(1.25*rmedian,rmedian+50.0)
	    !dtest = rmedian+3.0*rstddev
	    !print *, '2', rmedian, dtest, dtest1, count
	    do while (nhigh .gt. imedian .and. val_sort(nhigh) .gt. dtest)
		nhigh = nhigh-1
	    end do

	    if (nhigh .lt. nlow) then
		print *, (val_sort(i),i=1,count)
		print *, nlow_old, nhigh_old, nvote_old, imedian, rmedian
		print *, nlow, nhigh
		stop 'skyd_combine, oops 3'
	    end if

	    dtest = min(0.8*rmedian,rmedian-50.0)
	    !dtest = rmedian-3.0*rstddev
	    !print *, '2', rmedian, dtest, dtest1, count
	    do while (nlow .lt. imedian .and. val_sort(nlow) .lt. dtest .and. btest(flag_sort(nlow),1))
		nlow = nlow+1
	    end do

	    if (nhigh .lt. nlow) then
		print *, (val_sort(i),i=1,count)
		print *, nlow_old, nhigh_old, nvote_old, imedian, rmedian
		print *, nlow, nhigh
		stop 'skyd_combine, oops 4'
	    end if

	    nvote = nhigh-nlow+1

	    if (nvote .ne. nvote_old) rmedian = ArrR4Median(nvote,val_sort(nlow),0.0,i,val(nlow))

	end if

	if (iType .eq. 0) then					! Mean
	    rgood = ArrR4Mean(nvote,val_sort(nlow),i)
	else if (iType .eq. 1) then				! Median
	    rgood = rmedian
	else if (iType .eq. 2) then				! Mean/median (ratio)
	    rgood = rmedian
	    if (rgood .ne. badval) rgood = ArrR4Mean(nvote,val_sort(nlow),i)/rgood
	else if (iType .eq. 3) then				! Mean-median (difference)
	    rgood = rmedian
	    if (rgood .ne. badval) rgood = ArrR4Mean(nvote,val_sort(nlow),i)-rgood
	else if (iType .eq. 4) then
	    rgood = ArrR4Stdev(0,nvote,val_sort(nlow),rgood,i)	! "Std"
	else if (iType .eq. 5) then				! Median/std (ratio)
	    rgood = rmedian
	    if (rgood .ne. badval) then
		rstd = ArrR4Stdev(0,nvote,val_sort(nlow),rgood,i)
		if (rstd .ne. badval) then
		    rgood = rgood/rstd
		else
		    rgood = badval
		end if
	    end if
	end if

	if (iSilent .le. 1) then

	    if (nlow .ne. nlow_old) then
		i = iStr
		i = i+Int2Str(nlow-nlow_old	, cStr(i+1:))
		i = i+Str2Str('/'		, cStr(i+1:))
		i = i+Int2Str(count		, cStr(i+1:))+1
		i = Str2Str('low votes rejected', cStr(i+1:))
		call Say(cSay,'I','bin',cStr)
		write (*,*) (val_sort(i),i=nlow_old,nlow-1)
	    end if

	    if (nvote .ne. 0) then
		i = iStr
		i = i+Int2Str(nvote	    , cStr(i+1:))
		i = i+Str2Str('/'	    , cStr(i+1:))
		i = i+Int2Str(count	    , cStr(i+1:))+1
		i = Str2Str('votes accepted', cStr(i+1:))
		call Say(cSay,'I','bin',cStr)
		write (*,*) (val_sort(i),i=nlow,nhigh)
	    end if

	    if (nhigh .ne. nhigh_old) then
		i = iStr
		i = i+Int2Str(nhigh_old-nhigh    , cStr(i+1:))
		i = i+Str2Str('/'		 , cStr(i+1:))
		i = i+Int2Str(count		 , cStr(i+1:))+1
		i = Str2Str('high votes rejected', cStr(i+1:))
		call Say(cSay,'I','bin',cStr)
		write (*,*) (val_sort(i),i=nhigh+1,nhigh_old)
	    end if

	    i = iStr
	    i = i+Str2Str('Final value:', cStr(i+1:))+1
	    i = i+Flt2Str(rgood, 3	, cStr(i+1:))+1
	    if (near_star .eq. 1) i = i+Str2Str('(near star)',cStr(i+1:))+1
	    i = i+Str2Str('<=='		, cStr(i+1:))

	    call Say(cSay,'I','==>',cStr)

	end if

	smei_skyd_combine = nvote

	! Fraction of total adus in rejected votes: 1:nlow-1 where too low; nhigh+1,count were too high

	rbad = (ArrR4Total(nlow-1,val_sort,i)+ArrR4Total(count-nhigh,val_sort(nhigh+1),i))/ArrR4Total(count,val_sort,i)

	! Build a 'bad pixel' map: the counter for the originating pixel in the SMEI frame is incremented
	! Rejected votes are 1:nlow-1 (too low) and nhigh+1,count (too high).

	do i=1,nlow-1
	    j = pos(indx(i))	    		
	    badpix(j) = badpix(j)+1
	end do

	do i=nhigh+1,count
	    j = pos(indx(i))	    		
	    badpix(j) = badpix(j)+1
	end do

	return
	end