Viewing contents of file '../idllib/contrib/groupk/avg_lc.pro'
;+
; NAME:
; AVG_LC
;
; PURPOSE:
; Plots the average light curves for a source, customized to
; the HEAO A-1 data.
;
; CATEGORY:
; HEAO A-1.
;
; INPUTS:
;
; Time: Array of times in SECONDS since 1/1/77, [dblarr( npts )].
;
; Flux: Array of corresponding intensities for the source,
; [fltarr( npts )] in units specified by the FLUX_UNITS keyword.
;
; Err: Array of corresponding uncertainties for each intensity
; [fltarr( npts )] in units specified by the FLUX_UNITS keyword.
;
; SrcName: The name of the source, [string].
;
; OPTIONAL INPUT PARAMETERS:
;
; FLUX_UNITS: The units of the given intensities, [string], (''=Default).
;
; BINSIZE: Time binsize of the average light curve in DAYS.
;
; OUTPUTS:
;
; This routine makes light curves of the HEAO A-1 data and sends
; them to the graphic device if the USER selects it from the widget
; menu.
;
; MODIFICATION HISTORY:
; September, 1994 H.C. Wen
; 12-JAN-1995 Added the BINSIZE keyword
; 01-FEB-1995 Implemented XPLOT, XOPLOT
; 08-FEB-1995 Eliminated the HALF keyword for simplicity.
;-
pro AVG_LC, Time, Flux, Err, SrcName, FLUX_UNITS=Flux_units, $
BINSIZE=Binsize
JD_1_1_1977 = 2443145L
xtitle = 'Modified Julian Date ( JD - 2443145 )'
tday = Time/86400.D0
Pmulti = !P.MULTI
!P.MULTI = [0,1,2]
units = ''
if keyword_set( Flux_units ) then units=Flux_units
flux1 = flux
sigF = Err
avgf = AVG( flux1 )
siga = sqrt( TOTAL( sigF^2. ) )/N_ELEMENTS( flux1 )
print,'Avg. Flux: ',arr2str(avgf),plusminus(),arr2str(siga)
; Plot the fitted intensities or "flux" as a function of time
MINy = MIN( flux1-sigF )
MAXy = MAX( flux1+sigF )
MINy = FLOOR( (0.9*MINy)*100 ) / 100
MAXy = CEIL( (1.1*MAXy)*100 ) / 100
yrng1= [ MINy, MAXy ]
mintd= MIN( tday, MAX=maxtd )
xrng1= [mintd-0.2,maxtd+0.2]
xplot,tday,flux1,sigF,psym=3,/xstyle,$
xrange=xrng1, $
yrange=yrng1, $
XTITLE=xtitle,$
YTITLE='Intensity '+units,$
TITLE ='Light Curve of '+srcname,$
/PORTRAIT, WINDOW=i
; Now do a weighted histogram of these intensities
wgt = 1./(sigF^2.)
hist_wgt = wgt * flux1
if N_ELEMENTS(Binsize) eq 0 then $
binsize = ((maxtd - mintd)/10.)*1.01
sumf = HIST1D( tday, hist_wgt, BINSIZE=Binsize, $
MIN=mintd, MAX=maxtd,$
OBIN=tdbins, DENSITY=dens )
wgtf = HIST1D( tday, wgt, BINSIZE=Binsize, $
MIN=mintd, MAX=maxtd)
igt0 = where( wgtf gt 0, ngt0 )
avgf = sumf
avgf(igt0) = avgf(igt0)/wgtf(igt0)
sigf = 0.*avgf
sigf(igt0) = 1./sqrt( wgtf(igt0) )
sigtd = replicate( (binsize/2.)*.95, n_elements( tdbins ))
; plot results
xplot,tdbins,avgf,sigf,psym=3,/xstyle,$
xrange=xrng1,$
yrange=yrng1,$
XTITLE=xtitle,$
YTITLE='Weighted Average Intensity '+units,$
TITLE ='Average Light Curve of '+srcname,$
/PORTRAIT
xoplot,tdbins,avgf,psym=10,/PORTRAIT
!P.MULTI = Pmulti
end