pro calhrs_dio,fname,udls,ih,data,eps,err,log
;+
;*NAME:
;			calhrs_dio
;
;*PURPOSE:
; Perform correction for diode non-linearities
;
;*PARAMETERS:
; CALLING SEQUENCE:
;	calhrs_dio,fname,udls,ih,data,eps,err,log
;
; INPUTS:
;	fname - reference file name
;	udls - array of udls for the observation
;	ih - headers for the data (128 x N)
;
; INPUT/OUTPUTS:
;	data - data arrays (512 x N)
;	eps - epsilon values
;
; OPTIONAL INPUT/OUTPUTS:
;	err - propagated statistical errors
;		If not supplied or set to a scalar, no error propagation
;		is done.
;	log - processing log
;
; HISTORY:
;	version 1  D. Lindler April 1989
;       Apr 17 1991      JKF/ACC    - handle degenerating array (N)
;       Apr 19 1991      JKF/ACC    - avoid bad UDL data (incorrect IH)
;       Apr 22 1991      DJL/ACC    - new version of CALHRS
;	Apr 24 1991	 DJL/ACC    - replaced NBADM with BADM
;-
;------------------------------------------------------------------------

	VERSION = 1.0
	DEAD = 400			;dead diode data quality flag
	MINDIO = 0.02			;threshold for dead diodes

	if n_elements(err) gt 1 then error_process=1 else error_process=0
;
; Read reference file
;
	sxopen,8,strtrim(fname,2)
	resp=sxread(8)
	close,8
;
; Extract combaddition pattern from first UDL
;
	ixdef = udls(37:41)
	rcodes = [1,udls(63:68)]
	max_rcode=udls(34)>1				;AVOID BAD IH VALUES
	last=where(ixdef eq 0)
	if !err lt 1 then ninit=5 else ninit=last(0)	;number of def. pairs	
	ixdef=ixdef(0:ninit-1)
	offsets=fix((ixdef-ixdef(0))/8.0+0.5)		;combaddition offsets
							; in diode units
	ddlink = (ih(0) mod 10) eq 0
	if ddlink then ninit=1
	s=size(ih) & nreads=n_elements(ih)/128
;
; History processing
;
	hist=strarr(3)
	hist(0)='CALHRS_DIO version '+string(version,'(f5.2)')+ $
		': Removal of diode non-uniformities'
	hist(1)='    Reference file = '+strtrim(fname,2)
	hist(2)='    Response computed for combaddition over '+ $
		string(ninit,'(i1)')+' diodes'
	if !dump gt 0 then printf,!textunit,hist
	if n_elements(log) gt 0 then sxaddhist,hist,log

;
; Easy case (no combaddition) -----------------------------------------------
;
	if ddlink or (ninit eq 1) then begin
		bad=where(resp lt mindio)
		nbad=!err
		if nbad gt 0 then resp(bad)=1
		for i=0,nreads-1 do begin
			d=data(*,i)/resp
			e=eps(*,i)
			if nbad gt 0 then begin
				e(bad)=e(bad)>dead
				d(bad)=0
			endif
			data(0,i)=d
			eps(0,i)=e
			if error_process then begin
				e=err(*,i)
				if nbad gt 0 then e(bad)=0
				err(0,i)=e/resp
			endif
			ih(66,i)=ih(66,i) or 8
		endfor
;
; hard case (combaddition) ---------------------------------------------------
;
	   end else begin
		main=resp(6:505)			;main array diodes
;
; compute multiplexed array for complete substep pattern
;
		cresp=main
		for i=1,ninit-1 do $
			cresp=cresp+[fltarr(offsets(i)),main(0:499-offsets(i))]
		cresp=[resp(0:5),cresp/ninit,resp(506:511)]
		badc=where(cresp lt mindio) & nbadc=!err
		if nbadc gt 0 then cresp(badc)=1
;
; loop on readouts
;
		for i=0,nreads-1 do begin
		    d=data(*,i)
		    e=eps(*,i)
		    if error_process then er=err(*,i)
		    ncoadds=ih(46,i)
		    bin=ih(0,i) mod 10	;bin number
		    rcode=rcodes(bin-1)	;repeat code for the bin
		    coadds_per_deflection=max_rcode/rcode
		    coadds_per_pattern=coadds_per_deflection*ninit
;
; completed pattern
;
		    if (ncoadds mod coadds_per_pattern eq 0) then begin
			d=d/cresp
			if nbadc gt 0 then begin
				d(badc)=0
				e(badc)=dead
			endif
			if error_process then begin
				er=er/cresp
				if nbadc gt 0 then er(badc)=0
			endif
		    end else begin
;
; Worse case (combadditon and incomplete pattern
;
		    ncoadds_per_def=lonarr(ninit)	;coadds for each def.
							;	pair
;
;	complete patterns
;
		    complete=ncoadds/coadds_per_pattern ;completed patterns
		    ncoadds_per_def(*)=complete*coadds_per_deflection
		    left=ncoadds-complete*coadds_per_pattern
;
;	completed deflection pairs
;
		    ndef=left/coadds_per_deflection
		    if ndef gt 0 then ncoadds_per_def(0)= $
					ncoadds_per_def(0:ndef-1)+ $
					coadds_per_deflection
		    left = left-ndef*coadds_per_deflection
;
;	coadds for last incomplete def. pair
;
		    ncoadds_per_def(ndef)=ncoadds_per_def(ndef)+left
;
; compute fraction time on each deflection pair
;
		    frac=ncoadds_per_def/total(ncoadds_per_def)
;
; compute multiplexed response
;
		    mresp=main*frac(0)
		    for j=1,ninit-1 do $
			mresp=mresp + frac(j)* $
				[fltarr(offsets(j)),main(0:499-offsets(j))]
		    mresp=[resp(0:5),mresp,resp(506:511)]
		    badm=where(mresp lt mindio)
		    nbadm=!err
		    if nbadm gt 0 then mresp(badm)=1.0
		    d=d/mresp
		    if nbadm gt 0 then begin
				d(badm)=0
				e(badm)=dead
		    endif
		    if error_process then begin
				er=er/mresp
				if nbadm gt 0 then er(badm)=0
		    endif
		    endelse
;
; place into output
;
		ih(66,i)=ih(66,i) or 8
		data(0,i)=d
		eps(0,i)=e
		if error_process then err(0,i)=er
	   endfor
	endelse
return
end