;+
; NAME:
;	GAIN_PLOT4
; PURPOSE:
;	Plots the detector gain and energy resolution for the two 
;	BCS flight detectors from the data contained in gain log.
;	Plot is a panel of 4 plots with BCS-A resolution and gain on 
;	the (top and bottom) lefthand side, and BCS-B resolution and 
;	gain on the (top and bottom) righthand side.
; CALLING SEQUENCE:  
;	.run gain_plot4
; INPUT:
;       Reads from gain log table from the file GAIN_LOG.DAT. 
;       Looks FIRST for this on the current directory, OTHERWISE defaults 
;       to directory defined in logical/environment variable $DIR_BCS_LOGS.
; OUTPUT:
;	plots can be directed to PS if desired
; OPTIONAL OUTPUT:
;
; HISTORY
;	RDB   May 92	From GAIN_PLOT.PRO
;	RDB   7-Jun-92  Added polynomial fit to gain of BCS-A since launch
;	RDB  21-Jun-92	Capitalzed $DIR_BCS_LOGS for Unix
;			Used CONCAT_DIR
;	RDB  04-May-93	Changed limits for line fit plot
;	RDB  04-Nov-93	Searches current dir for file, else uses default.
;			Needed after centralizing to isass0::
;	RDB  13-Jan-94	Also do fit to BCS_B
;	RDB  23-Jan-94	Changed range of gain of BCS-A and deriv. of ytickv 
;	RDB  02-Nov-95	Changed range of gain, etc for BCS-A and BCS-B
;       RDB  07-Dec-95  Upped no. of points in arrays...
;
;-

yesnox,'* Plot fit to the slope of gain changes',ans_fit
print,' '

npnts = 200
centroid = fltarr(4,npnts)
fwhm = fltarr(4,npnts)
days = fltarr(npnts)
eggse = intarr(npnts)

qqfil = concat_dir([curdir(),'$DIR_BCS_LOGS'],'gain_log.dat')
exist = file_exist(qqfil)

if exist(0) then openr,dlun,qqfil(0),/get_lun  $
  else openr,dlun,qqfil(1),/get_lun
jjfil = fstat(dlun)
print,' '
print,'>>> Log file used: ',jjfil.name
print,' '

;;get_lun,dlun
;;openr,dlun,concat_dir('$DIR_BCS_LOGS','gain_log.dat')
;	skip title stuff
ss=''
for j=0,4 do begin readf,dlun,ss & print,ss & endfor

date=intarr(3)
vars=fltarr(8)
jj=0
ieg=''

while not eof(dlun) do begin
	readf,dlun,format='(i2,2(x,i2),a2,4(2x,2f7.3))',date,ieg,vars
;;	day = fix(days79([date(2),date(1),date(0)])) - fix(days79([90,1,1]))
	ex2int,[0,0,0,0,1,1,90],msody,days_90
	ex2int,[0,0,0,0,date],msody,days_curr
	day = fix(days_curr - days_90)
	print,format='(3i3,i5,8f8.3,a2)',date,day,vars,ieg
	qq = vars([0,2,4,6]) & centroid(4*jj) = qq
	qq1 = vars([1,3,5,7]) & fwhm(4*jj) = qq1
	days(jj) = day
	if strupcase(strmid(ieg,1,1)) eq 'E'  then eggse(jj)=1
	jj = jj+1
	endwhile

close,dlun
free_lun,dlun

centroid = centroid(*,0:jj-1)
fwhm = fwhm(*,0:jj-1)
days = days(0:jj-1)		; - fix(days79([90,1,1]))
eggse = eggse(0:jj-1)
qq = where(eggse gt 0)
print,'No of points read:',jj,n_elements(qq)
print,' '

;	define which point used as zero for gain plot
q0 = min(where(eggse eq 0))	;use first point when on s/c
q0 = min(where(days eq 632))	;use first point after launch
navg = 4
qx = where(days eq 669)		;for error bar

;	define if all points joined
ijoin = 1

ans=''
if !d.name ne 'PS' then  read,'Plot',ans
!p.multi = [0,2,2,0,1]
!x.margin = [9,2]
if !d.name eq 'PS' then begin
   device,/landscape
endif

;tek
set_xy
!xtitle = 'Days Since 1/1/90' + $
  string("   (Last point: ",date,fix(days(jj-1)),")",'(a,i2,2(1h/,i2),i5,a)')
!p.charsize=0.80
!x.style=1 & !y.style=1

!grid=0
!mtitle = 'BCS-A Gain History  (Detector F2)'
!y.range = [13,23]	;!ymin=13 & !ymax=23
!ytitle = 'Resolution (percent)'

;;!x.range = [300.,(fix(days(jj-1))+110)/50*50]
xtoday = anytim2ints(today())
xtoday = xtoday.day - days_90
!x.range = [300.,(xtoday+110)/50*50]

cfwhm0 = fwhm(0,*) & vv = where(cfwhm0 gt 0) & !c=0
wid0 = cfwhm0/centroid(0,*)*100.
cfwhm1 = fwhm(1,*) & vv = where(cfwhm1 gt 0) & !c=0
wid1 = cfwhm1/centroid(1,*)*100.

vv = where(cfwhm0 gt 0 and cfwhm1 gt 0) & !c=0
vve = where(cfwhm0 gt 0 and cfwhm1 gt 0 and eggse eq 0) & !c=0
wids = (wid0+wid1)/2.
!c=0
!p.psym=1
plot,days(vv),wids(vv)
!p.psym=-1
oplot,days(vve),wids(vve)
!p.thick=2
oploterr,[days(qx),days(qx)],[wids(qx),wids(qx)],[0.75,0.75],1
;;oplot,[days(qx),days(qx)],[wids(qx)+0.75,wids(qx)-0.75]
!p.thick=0

if !d.name ne 'PS' then  read,ans

!grid=0		;1
!ytitle = 'Gain'
gsta = 0.84 & gra = .18 & ngra = fix(100*gra)/2
!yticks=ngra
!y.range = [gsta,gsta+gra]	;0.95,1.05]
y_tickv = str2arr(string(findgen(ngra+1)*0.02+gsta,format='(f4.2,20(1h,,f4.2))'))

cgain0 = centroid(0,*) & gain00 = total(cgain0(q0:q0+navg))/n_elements(where(cgain0(q0:q0+navg) gt 0))
ddgain0 = (cgain0 - gain00)/gain00 + 1.
vv = where(cgain0 gt 0) & !c=0
vve = where((cgain0 gt 0) and (eggse eq 0))

cgain1 = centroid(1,*) & gain00 = total(cgain1(q0:q0+navg))/n_elements(where(cgain1(q0:q0+navg) gt 0))
ddgain1 = (cgain1 - gain00)/gain00 + 1.
vv = where(cgain1 gt 0) & !c=0
vve = where((cgain1 gt 0) and (eggse eq 0))

!mtitle=''
vvx  = where((cgain0 gt 0) and (cgain1 gt 0))
vvxe = where((cgain0 gt 0) and (cgain1 gt 0) and (eggse eq 0))
ddgain = (ddgain0+ddgain1)/2.
!c=0
!psym = 1
plot,days(vvx),ddgain(vvx),ytickname=y_tickv
!psym = -1
oplot,days(vvxe),ddgain(vvxe)
!p.thick=2
oploterr,[days(qx),days(qx)],[ddgain(qx),ddgain(qx)],[0.01,0.01],1
;;oplot,[days(qx),days(qx)],[ddgain(qx)-.01,ddgain(qx)+.01]
!p.thick=0

;	do a fit to the slope on BCS-A
since_day = 680
qq = where(days gt since_day and ddgain gt 0)
cc = poly_fit(days(qq),ddgain(qq),1,yfit,yband,sigma)
print,'Fit to BCS-A gain since day',since_day
print,'Fit coefficents =',cc,format='(a,f10.6,e16.6)
print,'Sigma of fit    =',sigma

if ans_fit then begin
   print,' '
   percent_gain = string(cc(1)*100*365,format='(f6.2)')+'%'
   print,'Decline per year for BCS-B = ',percent_gain
   fdays = [610,!x.range(1)]		;launch day = 606
   !p.psym = 0
   oplot,fdays,poly(fdays,cc),linestyle=1
   set_xy,0,1,0,1
   xyouts,.70,.87,percent_gain + ' p.a.'
   set_xy
endif


if !d.name ne 'PS' then  read,ans

!mtitle = 'BCS-B Gain History  (Detector F4)'

!grid=0
!yticks=0
!y.range = [14,24]	;13,23]	;!ymin=13 & !ymax=23
!ytitle = 'Resolution (percent)'

;;!x.range = [300.,(fix(days(jj-1))+110)/50*50]
!x.range = [300.,(xtoday+110)/50*50]

cfwhm0 = fwhm(2,*) & vv = where(cfwhm0 gt 0) & !c=0
wid0 = cfwhm0/centroid(2,*)*100.
cfwhm1 = fwhm(3,*) & vv = where(cfwhm1 gt 0) & !c=0
wid1 = cfwhm1/centroid(3,*)*100.

vv = where(cfwhm0 gt 0 and cfwhm1 gt 0) & !c=0
vve = where(cfwhm0 gt 0 and cfwhm1 gt 0 and eggse eq 0) & !c=0
wids = (wid0+wid1)/2.
!c=0
!p.psym=1
plot,days(vv),wids(vv)
!p.psym=-1
oplot,days(vve),wids(vve)
!p.thick=2
oploterr,[days(qx),days(qx)],[wids(qx),wids(qx)],[0.75,0.75],1
!p.thick=0

if !d.name ne 'PS' then  read,ans

!grid=0		;1
!ytitle = 'Gain'
gsta = 0.84 & gra = .18 & ngra = fix(100*gra)/2
!yticks=ngra
!y.range = [gsta,gsta+gra]	;0.95,1.05]
y_tickv = str2arr(string(findgen(ngra+1)*0.02+gsta,format='(f4.2,20(1h,,f4.2))'))

cgain0 = centroid(2,*) & gain00 = total(cgain0(q0:q0+navg))/n_elements(where(cgain0(q0:q0+navg) gt 0))
ddgain0 = (cgain0 - gain00)/gain00 + 1.
vv = where(cgain0 gt 0) & !c=0
vve = where((cgain0 gt 0) and (eggse eq 0))

cgain1 = centroid(3,*) & gain00 = total(cgain1(q0:q0+navg))/n_elements(where(cgain1(q0:q0+navg) gt 0))
ddgain1 = (cgain1 - gain00)/gain00 + 1.
vv = where(cgain1 gt 0) & !c=0
vve = where((cgain1 gt 0) and (eggse eq 0))

!mtitle=''
vvx  = where((cgain0 gt 0) and (cgain1 gt 0))
vvxe = where((cgain0 gt 0) and (cgain1 gt 0) and (eggse eq 0))
ddgain = (ddgain0+ddgain1)/2.
!c=0
!psym = 1
plot,days(vvx),ddgain(vvx),ytickname=y_tickv
!psym = -1
oplot,days(vvxe),ddgain(vvxe)
qx = where(days eq 669)
!p.thick=2
oploterr,[days(qx),days(qx)],[ddgain(qx),ddgain(qx)],[0.01,0.01],1
!p.thick=0

;	do a fit to the slope on BCS-B
since_day = 850
qq = where(days gt since_day and ddgain gt 0)
cc = poly_fit(days(qq),ddgain(qq),1,yfit,yband,sigma)
print,'Fit to BCS-B gain since day',since_day
print,'Fit coefficents =',cc,format='(a,f10.6,e16.6)
print,'Sigma of fit    =',sigma

if ans_fit then begin
;;   print,'Decline per year for BCS-B = ',cc(1)*100*365,' %',format='(a,f6.2,a)'
   print,' '
   percent_gain = string(cc(1)*100*365,format='(f6.2)')+'%'
   print,'Decline per year for BCS-B = ',percent_gain
   fdays = [610,!x.range(1)]		;launch day = 606
   !p.psym = 0
   oplot,fdays,poly(fdays,cc),linestyle=1
   set_xy,0,1,0,1
   xyouts,.70,.87,percent_gain + ' p.a.'
endif


;;if !d.name ne 'PS' then  read,ans

!type=0 & !yticks=0
!grid=0 & !p.psym=0

print,' '
print,'>>> Log file used: ',jjfil.name
print,' '

print,'>>> Days since last measurement =',fix(xtoday-days(jj-1))
print,' '

end