PRO pfold, time,rate,profile,period=period,nbins=nbins,time0=time0, $ proferr=proferr,gap=gap,nogap=nogap,tolerance=tolerance, $ phbin=phbin,ttot=ttot, raterr=raterr,chatty=chatty,npts=npts, $ quantiles=quantiles, statquant=quant,dt=dt,gti=gti, $ pdot=pdot,pddot=pddot,freq=nu,nudot=nudot,nuddot=nuddot, $ mjd=mjd ;+ ; NAME: ; pfold ; ; ; PURPOSE: ; folds a lightcurve with a given period or frequency and returns ; the resulting profile ; ; ; CATEGORY: ; lightcurve ; ; ; CALLING SEQUENCE: ; pfold, time,rate,profile,period=period,nbins=nbins,time0=time0, $ ; proferr=proferr,gap=gap,nogap=nogap,tolerance=tolerance, $ ; phbin=phbin,ttot=ttot, raterr=raterr,chatty=chatty,npts=npts, $ ; quantiles=quantiles, statquant=quant,dt=dt,gti=gti, $ ; pdot=pdot,pddot=pddot ; ; INPUTS: ; time : the starting time of each rate bin (seconds, unless ; mjd keyword is given) ; rate : the count rate of each bin (counts/second, ; independent of the mjd keyword!) ; if the rate argument is omitted, the code assumes ; EVENT data instead, and time is the time at which ; each event was detected) ; ; ; OPTIONAL INPUTS: ; period : period on which to fold (seconds) ; freq : frequency on which to fold (Hz) ; -- only one of period and freq can be given ; pdot, pddot: Pdot and Pdotdot (pulse period derivatives) ; if pddot is set, pdot needs to be defined as well ; nudot, nuddot: nu_dot and nu_dotdot (frequency derivatives) ; if nuddot is set, nudot needs to be defined as well ; nbins : the number of phase bins to use ; time0 : use time0 instead of the first entry in the ; timecolumn as zerotime of phase bin, same units as ; time ; gti : an array of good time intervals ; gti(0,*) contains the starttimes, ; gti(1,*) contains the stoptimes. ; if given, the exposure time per phase bin (ttot) is ; computed from this information. alternatively, for ; lightcurves, the gap detection can be used. ; raterr : uncertainty of the counting rates in rate(*) ; dt : width of each time bin ; gap : An array containing the indices of gaps in the ; lightcurve (i.e. time(gap(2)) is the starting time ; of the 3rd gap in the lightcurve). These bins are ; ignored in the epoch folding. ; If nogap is not set and no gti is given, the code determines ; gap using the subroutine timegap and returns gap ; If nogap is not set and a gti is known, the gti is ; used to determines the bad bins and returns these ; in gap ; A result of these options is that for a given ; lightcurve the gaps are only determined once. ; ; tolerance: parameter defining the lower limit for the gap ; length; the reference is the time difference ; between the first and second entry in the time ; array; tolerance defines the maximum allowed relative ; deviation from this reference bin length; ; default: 1e-8; this parameter is passed to timegap ; (see timegap.pro for further explanation) ; Be aware that tolerance needs to be set to a higher ; value if a lightcurve is barycentered after it is ; binned. Only used if gti is not given. ; ; ; ; ; KEYWORD PARAMETERS: ; nogap : If set, data are assumed to NOT contain any ; gaps. Has only influence for lightcurves. ; quantiles: If set, return the 1st,2nd, and 3rd quantile in ; keyword statquant (only for lightcurves) ; mjd : if set, then time argument is in days, not seconds. ; ; OUTPUTS: ; profile : the pulse profile (counts/sec) ; ; OPTIONAL OUTPUTS: ; proferr : Array containing the statistical uncertainty of each ; bin of the profile. computed from raterr for ; lightcurves, for event data, Poissonian errors are ; assumed. ; phbin : array containing the start-phase of each bin of the ; profile ; ttot : array containing the total observing time ; accumulated in each bin of the profile ; note that ttot is not 100% exact if pdot/pddot are ; nonzero ; npts : total number of rate points in each bin of the ; profile, or total number of events in each bin. ; statquant: a 2d array, quant(0,*) contains the 1st, quant(1,*) ; the 2nd, and quant(2,*) the 3rd quantile for each ; pulse point; only for lightcurve data ; ; RESTRICTIONS: ; * it is assumed that the period/frequency are only changing ; slowly, such that it can be assumed that during one gti ; interval the period is close to constant. For all relevant ; astronomical objects this seems to be the case. In other ; words, when the GTI is known, then ttot (the exposure time of ; each phase bin) is computed using the period as determined ; for the middle of the GTI interval as the correct period. ; Note that the phase of each phase bin is determined correctly ; (using getphase) ; ; ; PROCEDURE: ; ; ; EXAMPLE: ; nbins=128 ; period=1.7 ; pfold,time,counts,profile,283.,nbins=nbins,period=period ; phase=findgen(nbins)/(nbins-1) ; plot,phase,profile ; ; ; MODIFICATION HISTORY: ; $Log: pfold.pro,v $ ; Revision 1.16 2006/02/14 17:54:47 wilms ; now also works in the case that the period is an integer variable ; ; Revision 1.15 2005/05/14 20:55:39 wilms ; added mjd keyword ; ; Revision 1.14 2004/11/22 07:59:45 wilms ; further improvement on handling of phase bins with no exposure (now ; also returns NaN for binned data). ; ; Revision 1.13 2004/11/21 19:10:10 kuster ; - fixed bug for phase bins with zero integration time (event list mode) ; - change print to message in chatty mode ; ; Revision 1.12 2004/11/01 21:05:57 wilms ; timezero warning only if chatty is set ; ; Revision 1.11 2004/09/13 16:02:56 wilms ; cosmetic changes only ; ; Revision 1.10 2004/09/07 17:26:50 wilms ; * added options to fold on frequency instead of period ; * added ability to have time dependent period (pdot and pddot and ; nudot, nuddot options) ; * switched phase handling over to getphase ; ; Revision 1.9 2004/09/06 20:54:04 wilms ; added treatment of GTIs and also greatly improved general treatment of ; gaps for the case of lightcurves. ; ; Revision 1.8 2004/05/24 18:03:30 wilms ; simpler computation of npts for events ; ; Revision 1.7 2003/06/30 10:26:25 wilms ; complete rewrite (Ljuba/JW), now also allows event lists ; ; ; Version 1.0: Derived from foldlc, Joern Wilms, 1998/05/14 ; Version 1.1: J.W., 1998/08/09, added quantile keyword ; Version 2.0: J.W., 1999/08/04: ; * now use the histogram function for determining the ; countrates contributing to each phase bin ; (speed improvement roughly a factor of 20%) ; * compute width of time bins only if desired ; (slight speed improvement) ; CVS Version 1.2: JW/PR, 2001.07.20 ; code now produces exact results in case time-time0 is negative ; CVS Version 1.3: MK, 2001.10.18 ; added keyword tolerance ; Version 1.4: JW, 2002.02.12 ; repaired tolerance keyword. now calls with tolerance eq ; 0. are possible ; Version 1.5: JW/SF, 2002.08.15 ; correction of code in the example ; Version 1.6: JW: CHANGE TO CVS LOGGING ;- ;; two non-keyword arguments: we're working on event data ;; three non-keyword arguments: we're working on a light curve eventdata=0 IF (n_params() EQ 2) THEN eventdata=1 IF (n_params() LT 2 OR n_params() GT 3) THEN BEGIN message,'syntax: pfold,time,profile for event data',/informational message,' pfold,time,rate,profile for a light curve' ENDIF IF (keyword_set(chatty) AND time[0] EQ 0. AND n_elements(time0) NE 0) THEN BEGIN message,'time[0]=0 and you have set a time0; this will',/informational message,'most probably not give the result you expect.',/informational ENDIF IF (n_elements(time) LE 1) THEN BEGIN message,'Need more than one lightcurve bin or event' ENDIF IF (n_elements(time0) EQ 0) THEN time0=time[0] IF (n_elements(nbins) EQ 0) THEN nbins=20 IF ((n_elements(pddot) NE 0) AND (n_elements(pdot) EQ 0)) THEN BEGIN message,'Pddot set: need pdot as well' ENDIF ;; for the case that pdot is set, we always also assume that ;; pddot is given as well. if not: set it to zero, to make ;; the equations below somewhat easier IF ((n_elements(pdot) NE 0) AND (n_elements(pddot) EQ 0)) THEN BEGIN pddot=0d0 ENDIF ;; event data: no gti file known IF (eventdata EQ 1 AND n_elements(gti) EQ 0) THEN BEGIN message,'eventdata: GTI not known, using time of first',/informational message,' and last event. This is most probably NOT',/informational message,' what you want!',/informational gti=fltarr(2,1) gti[0,0]=time[0] gti[1,0]=time[n_elements(time)-1] ENDIF IF (eventdata EQ 1 AND keyword_set(quantiles)) THEN BEGIN message,'quantiles can only be computed for lightcurves' ENDIF ;; Set gap tolerance to default value (only for lightcurves) ttol=1d-7 IF (n_elements(tolerance) NE 0) THEN ttol=tolerance ;; disable chatty if not set IF (n_elements(chatty) EQ 0) THEN chatty=0 ;; Width of each time-bin, take care of wrap-around IF (n_elements(dt) EQ 0) THEN BEGIN dt=temporary(shift(time,-1))-time dt(n_elements(time)-1)=dt(0) ENDIF ;; ;; start of phase bins, including phase 1 (is removed at end of ;; code and makes code more elegant) ;; phbin=dindgen(nbins+1)/nbins ;; ;; If the gti file on which the lightcurve or event data are based ;; on is known, we can compute the exact exposure for each bin ;; (if not: compute ttot from the bin size) ;; ttot=dblarr(nbins) IF (n_elements(gti) NE 0) THEN BEGIN ;; ;; Loop over all good time intervals ;; FOR i=0, n_elements(gti[0,*])-1 DO BEGIN ;; ;; compute information for start and stop times ;; starttime=gti[0,i] stoptime=gti[1,i] ;; period number of start period nstart = long(getphase(starttime,time0=time0,/numpulse, $ period=period,pdot=pdot,pddot=pddot,$ freq=nu,nudot=nudot,nuddot=nuddot, $ mjd=mjd)) IF (starttime LT time0) THEN nstart = nstart - 1L ;; period number of stop period nstop = long(getphase(stoptime,time0=time0,/numpulse, $ period=period,pdot=pdot,pddot=pddot,$ freq=nu,nudot=nudot,nuddot=nuddot, $ mjd=mjd)) IF (stoptime LT time0) THEN nstop = nstop - 1L phasestart = getphase(starttime,time0=time0, $ period=period,pdot=pdot,pddot=pddot,$ freq=nu,nudot=nudot,nuddot=nuddot, $ mjd=mjd) ;; phase of stop period phasestop = getphase(stoptime,time0=time0, $ period=period,pdot=pdot,pddot=pddot,$ freq=nu,nudot=nudot,nuddot=nuddot, $ mjd=mjd) ;; nbstart -- number of start bin ;; nbstop -- nubber of stop bin nbstart = long(phasestart*nbins) nbstop = long(phasestop*nbins) ;; period for this phase interval ;; ;; we will assume that the GTI is always short with ;; respect to changes in the period, if we would not do ;; this then the code would be a real pain... ;; IF (n_elements(period) NE 0) THEN BEGIN IF (n_elements(pdot) EQ 0) THEN BEGIN perint=period ENDIF ELSE BEGIN tmid=(starttime+stoptime)/2. ddd=tmid-time0 IF (keyword_set(mjd)) THEN ddd=ddd*86400d0 perint=(0.5*pddot*ddd+pdot)*ddd+period ENDELSE ENDIF ELSE BEGIN ;; compute period from frequency IF (n_elements(nudot) EQ 0) THEN BEGIN perint=1d0/nu ENDIF ELSE BEGIN tmid=(starttime+stoptime)/2. ddd=tmid-time0 IF (keyword_set(mjd)) THEN ddd=ddd*86400d0 nuint=(0.5*nuddot*ddd+nudot)*ddd+nu perint=1d0/nuint ENDELSE ENDELSE ;; width of each time bin (also works if perint is ;; accidentally an integer [yes, this HAS happened on ;; 2006-02-14]) size_bin=double(perint)/nbins IF (nstart NE nstop) THEN BEGIN ;; ;; Case that the start and stop time of the GTI do ;; NOT fall into the same period bin ;; ;; ;; Take care of the periods between ;; the start and the stop interval ;; ttot[*] = ttot[*] + (nstop - (nstart+1L))*size_bin ;; ;; Take care of the first and the last period interval ;; ;; ... first period interval ttot[nbstart]=ttot[nbstart] + $ (phbin[nbstart+1] - phasestart)*perint IF (nbstart LT nbins-1) THEN BEGIN ttot[nbstart+1:nbins-1]=ttot[nbstart+1:nbins-1]+size_bin ENDIF ;; ... last period interval IF (nbstop NE 0) THEN BEGIN ttot[0:nbstop-1] = ttot[0:nbstop-1]+size_bin ENDIF ttot[nbstop]=ttot[nbstop] + $ (phasestop - phbin[nbstop])*perint ENDIF ELSE BEGIN ;; ;; GTI is shorter than one period ;; and falls into the same period interval ;; IF (nbstart EQ nbstop) THEN BEGIN ;; ;; GTI only falls into ONE phase bin ;; ttot[nbstart]=ttot[nbstart]+(phasestop-phasestart)*perint ENDIF ELSE BEGIN ;; ;; GTI falls into more than one phase bin ;; IF (nbstart+1 NE nbstop) THEN BEGIN ;; ... phase bins between start and stop bin ttot[nbstart+1:nbstop-1] = $ ttot[nbstart+1:nbstop-1] + size_bin ENDIF ;; ... start and stop bin ttot[nbstart]=ttot[nbstart] + $ (phbin[nbstart+1] - phasestart)*perint ttot[nbstop]=ttot[nbstop] + $ (phasestop - phbin[nbstop])*perint ENDELSE ENDELSE ENDFOR ENDIF ;; ... computation of the pulse profile ;; ;; Compute phase of each timebin ;; phase=getphase(time,time0=time0, $ period=period,pdot=pdot,pddot=pddot,$ freq=nu,nudot=nudot,nuddot=nuddot,mjd=mjd) ;; If lc contains gaps or gti is known: mark them invalid ;; IF ((eventdata EQ 0) AND NOT keyword_set(nogap)) THEN BEGIN IF (n_elements(gti) EQ 0 OR n_elements(gap) NE 0) THEN BEGIN ;; ;; no GTI known -> search for gaps and mark ;; time bins on gaps as invalid ;; IF (n_elements(gap) EQ 0) THEN BEGIN timegap,time,gap,dblock,tolerance=ttol,chatty=chatty IF (keyword_set(chatty) AND (n_elements(gap) GE 0.1*n_elements(time))) THEN BEGIN message,'WARNING: More than 10% gaps in the lightcurve',/informational message,'Could be due to too stringent tolerance setting',/informational message,'and barycentering the lightcurve',/informational ENDIF ENDIF IF (gap[0] NE -1) THEN phase[gap]=2. ENDIF ELSE BEGIN ;; ;; if we know the gti, then mark data outside of the GTI ;; or time bins crossing the GTI as invalid ;; (should in principle also take care of bins COVERING ;; a whole gti, but we will forget this for the moment) ;; note that the following is very inefficient and not ;; necessary where GTIs fall on time bin boundaries, but in ;; general we do not have this (unfortunately), so it is ;; better to keep this in. ;; ;; note that this returns gap, i.e., if pfold is called ;; multiple times, then the following code is executed ;; only once (optimization for pfold) ;; ;; ;; step 1: data clearly outside of the bins ;; ngti=n_elements(gti[0,*]) FOR i=1,ngti-1 DO BEGIN ndx=where(time GT gti[1,i-1] AND time LT gti[0,i]) IF (ndx[0] NE -1) THEN phase[ndx]=2. ENDFOR IF (time[0] LT gti[0,0]) THEN BEGIN ndx=where(time LT gti[0,0]) phase[ndx]=2. ;; NB we are guaranteed that ndx is not -1 ENDIF IF (time[n_elements(time)-1] GE gti[1,ngti-1]) THEN BEGIN ndx=where(time GE gti[1,ngti-1]) phase[ndx]=2. ENDIF ;; ;; step 2: cut out time bins crossing a gti ;; tstop=shift(time,-1) tstop[n_elements(time)-1]=time[n_elements(time)-1]+dt[0] FOR i=0,ngti-1 DO BEGIN ndx=where(time LT gti[0,i] AND tstop GT gti[0,i]) IF (ndx[0] NE -1) THEN phase[ndx]=2. ndx=where(time LT gti[1,i] AND tstop GT gti[1,i]) IF (ndx[0] NE -1) THEN phase[ndx]=2. ENDFOR ;; remember bad bins IF (NOT keyword_set(nogap)) THEN BEGIN gap=where(phase EQ 2.) ENDIF ENDELSE END profile=fltarr(nbins) proferr=fltarr(nbins) npts=lonarr(nbins) ;; ;; Get indices to all times falling into each phase bin ;; dummy=histogram(phase,binsize=1./nbins,min=0.,max=1., $ reverse_indices=r) ;; ;; Compute the profiles ;; IF (eventdata EQ 1) THEN BEGIN ;; ;; Event data ;; FOR i=0,nbins-1 DO BEGIN IF (r[i] EQ r[i+1]) THEN BEGIN IF (keyword_set(chatty)) THEN BEGIN message, 'Warning: phase '+strtrim(string(i),2)+ $ ' contains no bins',/informational ENDIF npts[i]=0L ;; not really needed... ENDIF ELSE BEGIN ;; jw, old code: ;; npts[i]=n_elements(r[r[i]:r[i+1]-1]) ;; .. which is equivalent to npts[i]=r[i+1]-r[i] ENDELSE ENDFOR ;; avoid division by zero ;; get index to phase bins with zero integration time ;; and apply a NAN to these bins ndx=where(ttot NE 0., complement=zerondx) IF (zerondx[0] NE -1) THEN BEGIN profile[zerondx]=!values.f_nan proferr[zerondx]=!values.f_nan ENDIF IF (ndx[0] NE -1) THEN BEGIN profile[ndx]=npts[ndx]/ttot[ndx] proferr[ndx]=sqrt(npts[ndx])/ttot[ndx] ENDIF ENDIF ELSE BEGIN ;; ;; Lightcurve data ;; IF (keyword_set(quantiles)) THEN BEGIN quant=fltarr(3,nbins) ENDIF FOR i=0,nbins-1 DO BEGIN IF (r[i] EQ r[i+1]) THEN BEGIN IF (keyword_set(chatty)) THEN BEGIN message, 'Warning: phase '+strtrim(string(i),2)+ $ ' contains no bins',/informational ENDIF profile[i]=!values.f_nan npts[i]=0 IF (n_elements(gti) EQ 0) THEN ttot[i]=0. END ELSE BEGIN ndx=r[r[i]:r[i+1]-1] dte=dt[ndx] IF (keyword_set(mjd)) THEN dte=dte*86400.d0 IF (n_elements(gti) EQ 0) THEN ttot[i]=total(dte,/double) npts[i]=n_elements(ndx) ;; need dt because of boundaries and GTI's profile[i]=total(rate[ndx]*dte,/double)/ttot[i] IF (n_elements(raterr) NE 0) THEN BEGIN proferr[i]=sqrt(total((raterr[ndx]*dte)^2.,/double))/ttot[i] ENDIF ;; compute statistical quantiles if desired IF (keyword_set(quantiles)) THEN BEGIN dist=rate[ndx] nnn=sort(dist) n25=fix(npts[i]*0.25) n50=fix(npts[i]*0.50) n75=fix(npts[i]*0.75) quant(0,i)=dist[nnn[n25]] quant(1,i)=dist[nnn[n50]] quant(2,i)=dist[nnn[n75]] ENDIF END ENDFOR ENDELSE ;; truncate last phase bin phbin=phbin[0:nbins-1] ;; event data: return profile in the 2nd argument IF (eventdata EQ 1) THEN rate=profile END