pro calfos_ppc,ppcoef,h,data,err,eps
;+
;			CALFOS_PPC
;
; FOS paired pulse correction
;
; CALLING SEQUENCE:
;	calfos_ppc,ppcoef,h,data,err,eps
;
; INPUTS:
;	ppcoef - paired pulse coef. (from routine CALFOS_CCG2)
;
; INPUT/OUTPUTS:
;
;	h - FITS header for the data
;	data - data array
;	err - propagated error array
;	eps - data quality array
;
; HISTORY:
;	version 1  D. Lindler  Jan 1990
;-
;-------------------------------------------------------------------------
;
; extract coefficients from PPCOEF
;
	tau1 = ppcoef(0)
	threshold = ppcoef(1)>0.01	;0.01 needed to avoid divide by zero
	q0 = ppcoef(2)
	q1 = ppcoef(3)
	F = ppcoef(4)
	rsat = ppcoef(5)
	r20 = ppcoef(6)
	r5 = ppcoef(7)
	if n_elements(err) gt 1 then propagate = 1 else propagate = 0
;
; set up history
;
	hist = strarr(3)
	hist(0) = 'Deadtime correction performed for count rates >'+ $
			strtrim (float(threshold))
	if q0 eq 0 then begin			;formula 1
		hist(1) = '  TAU1 = '+string(tau1)
	      end else begin				;formula 2
		hist(1) = '    Q0 ='+string(q0)+'    Q1 ='+string(Q1) + $
			'    F ='+ string(F)
	endelse
;
; flag bad and sort of bad data
;
	hist(2) = '  RSAT ='+string(rsat)+'   R20 ='+string(r20)+'   R5 ='+ $
			string(R5)
	bad = where(data gt rsat) & nbad = !err		;saturated
	if nbad gt 0 then begin
		data(bad) = 0.0
		eps(bad) = 300
		if propagate then err(bad) = 0.0
	endif

	bad = where(data gt r20) & nbad = !err		;20% error
	if nbad gt 0 then eps(bad) = eps(bad)>190

	bad = where(data gt r5) & nbad = !err		;5% error
	if nbad gt 0 then eps(bad) = eps(bad)>130
;
; perform correction on data greater than THRESHOLD
;
	perform = where(data gt THRESHOLD) & n=!err
	if n gt 0 then begin
	    values = data(perform)
	    if q0 eq 0 then begin			;formula 1
		x = (1.0-values*tau1) > 1.0e-12
		new_values = alog(x)/(-tau1)
		data(perform) = new_values
		if propagate then err(perform) = err(perform)*new_values/values
	      end else begin				;formula 2
		t = q0 + (values gt F)*(values - F)*q1
		x = (1.0 - values*t)>1.0e-12
	 	data(perform) = values/x
		if propagate then err(perform) = err(perform)/x
	    endelse
	endif

	if !dump gt 0 then print,hist
	sxaddhist,hist,h
return
end