;+
; NAME:
;	PPHA
; PURPOSE:
;	Extract BCS PHA data from the DP_SYNC data stream
;	plots PHA spectra of selected time interval
; CALLING SEQUENCE:  
;	.run PPHA
; INPUTS:
;
; OUTPUTS:
;	Plot of PHA spectra
; OPTIONAL OUTPUTS:
;	hardcopy is optional
; RESTRICTIONS:
;	Only works on and X-terminal
;	This programme should be run after PP_CARTEX or PLOT_ERATE
; PROCEDURE:
;	The user is allowed to select a region from the existing 
;	light curve. The PHA data corresponding to this interval
;	is extracted from the DP_SYNC data and is thrown up on the 
;	screen for the user to further select the required time interval 
;	The PHA spectra of the sum of this period is plotted and can
;	optionally be sent to the hardcopy.
; MODIFICATION HISTORY:
;	RDB	1991/92
;-

;	  reset to window zero

print,'** Light-curve programme must have been run with Window 0 **'
wset,0
!x = pp_x & !y = pp_y & !z = pp_z & !p = pp_p
!fancy=2

print,' '
print,'** Set Cursor to required region on Light Curve'
get_limits,xx,yy	;limits from time plot
lim_region = where((ttime ge xx(0)) and (ttime le xx(1)))
ista = min(lim_region) & ista = ista(0)
istp = max(lim_region) & istp = istp(0)
print,'ista, istp =',ista,istp

fioff = 0
fi = dp_sync.dp_time(3)
istax = min(where(fi(ista:*) eq fioff)) + ista & istax = istax(0)
nnpha = (istp - istax)/4
istp  = istax + nnpha*4
print,'ista, istp =',ista,istp
pha   = byte(dp_sync.pha_counts,istax*8,32,nnpha)	;8 PHA bytes/rec
print,'nnpha =',nnpha
help,pha


;	Plot contour map of the PHA data to aid channel selection
window,1,xsiz=300,ysiz=600,retain=2
next_one:
wset,1
set_xy
!p.multi = 0
contour,pha,levels=[2,5,10,20,50,100,200] $
	,xtitle='PHA bins',ytitle='PHA frames (Time)'	$
	,xticks=4,xrange=[0,32],xstyle=1

print,' '
print,'** Set Cursor to required region on contour plot'
get_limits,xx1,yy1			;required region on contour plot

jyy = [0,0]  &  jyy4 = [0,0]
jyy = fix(yy1+0.5)			;try to synch locally
jkkk = min(where(fi(jyy(0)*4+istax:jyy(1)*4+istax) eq 0)) 
jyy4(0) = jyy(0)*4 + jkkk + istax	;normal records
jyy4(1) = jyy4(0) + (jyy(1)-jyy(0)-1)*4
jyy = (jyy4 - istax)/4 + 1 		;?? should this be +1 ??
;;print,jyy,jyy4

print,' '
phach = (dp_sync(jyy4(0):jyy4(1)).pha_control and 'c0'x)/64
phach = all_vals(phach)
print,'PHA Channel   = ',phach
phatm = (dp_sync(jyy4(0):jyy4(1)).pha_control and '3f'x)
phatm = all_vals(phatm)*0.5 - 0.125
print,'PHA Int. time = ',string(phatm,'(f7.3)')

;         start of PHA accumulation
int2ex,dp_sync(jyy4(0)).time,dp_sync(jyy4(0)).day,tarr
sstm =  fmt_tim(tarr,sday,stime)

;	  SCA levels during this interval
scaidx = sel_qs(qs.instr,dp_sync(jyy4(0):jyy4(1)))
scaidx = all_vals(scaidx)
if n_elements(scaidx) gt 1 then print,'** ERROR - SCA values changed **',scaidx
sca = byte(qs.instr(scaidx(0)).scaval,0,2,4)	;extract and redimension
print,'SCAs= ',string(sca(*,phach(0)),'(2(x,Z2.2))')

;	  get PHA again locally and sum the data over the interval
ppha = dp_sync(jyy4(0):jyy4(1)).pha_counts
nppha = (jyy4(1)-jyy4(0))/4
ppha = byte(ppha,0,32,nppha)
phas = float(sum(ppha))/float(nppha)/phatm(0)

window,2,xsiz=300,ysiz=300,retain=2
ans=''

;	plot the extract and summed PHA data
rrplot:
!psym=10
set_xy
!xtitle = 'PHA Bins'
!ytitle = 'Counts/Bin/Sec'
cymax=max(phas(0:30))
!mtitle = string('BCS Chan.',phach(0),'  Int.',phatm(0),'s','(a,i2,a,f8.3,a)')
plot,phas,xrange=[0,32],xstyle=1,yrange=[0,1.1*cymax],ystyle=1
xyouts,0.9,0.05,filnam+ext,/device,siz=1.2
set_xy,0,1,0,1
xyouts,0.05,0.95,string(sstm,nppha,'(a,2x,i3)')
xyouts,0.05,0.90,'SCAs= '+string(sca(*,phach(0)),'(x,Z2.2,2h, ,Z2.2)')
;;;; 	NO! need to be in the correct bin space
;;;oplot,[sca(0,phach(0)),sca(0,phach(0))],[0,0.8*!y.crange(1)]
;;;oplot,[sca(1,phach(0)),sca(1,phach(0))],[0,0.8*!y.crange(1)]


;	Output can be copied to paper, or new region selected, or quit
if ans ne 'Y' then begin
  read,'* Plot to paper? [q to quit]: ',ans
  ans = strupcase(strmid(ans,0,1))
  if ans eq 'Y' then begin
	set_plot,'ps'
	goto, rrplot
	endif
  if ans eq 'Q' then goto, endit
  endif

if ans eq 'Y' and !d.name eq 'PS' then pltdv2

set_plot,'x'
goto, next_one

endit:
set_plot,'x'
end