epferror

             estimate epoche folding error with monte-carlo
             simulation approach.

             timing tools

             error=epferror(time,rate,period=period, raterr=raterr,
                    pstart=pstart,pstop=pstop, ntrials=ntrials,
                    nbins=nbins,sampling=sampling,gti=gti,chatty=chatty,tolerance=tolerance)

             time : a vector containing the time in arbitary units
             rate : a vector containing the countrate; if not given,
                    the times are assumed to come from individual
                    events.
            pstart: lowest period to be considered.
                    Propagated to epfold when analysing simulated
                    lightcurve.
             pstop: highest period to be considered.
                    Propagated to epfold when analysing simulated
                    lightcurve.

          ntrial: number of trials to apply epoche folding on
                  a simulated lightcurve. Default: 1000.
             gti: at the moment for event data only: gti[2,*] is an
                  array with the good times from the event selection,
                  gti[0,*] are the start times, gti[1,*] the stop
                  times. Obviously, use the same time system for time
                  and gti (i.e., also barycenter gtis!). Needed for
                  the computation of the exposure time of the phase
                  bins, if not given, the time of the first and the
                  last event is used instead.
                  Propagated to epfold when analysing simulated
                  lightcurve.
          raterr: a given error vector for the countrate. In
                  this case monte carlo simulation will be performed
                  with gaussian distribution using this error.
                  THEREFORE THIS ERROR IS INTERPRETED AS GAUSSIAN
                  STANDARD DEVIATION!!!
                  If not given and we are not working on events,
                  raterror is computed using poissonian statistics.
                  This should be the standard case.
                  REMARK: For high countrates/photon numbers (> 50)
                  it is recommended to use this keyword with
                  raterr=sqrt(rate) to avoid numerical problems with
                  poisson statistics.
             nbins:    number of phase-bins to be used in creating
                  the epoche folded trial profile. Default: 20.
                  Propagated to epfold when analysing simulated
                  lightcurve.
          sampling: how many periods per peak to use (default=10)
                  Propagated to epfold when analysing simulated
                  lightcurve.
       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 timgap.pro for further explanation)
                  Propagated to epfold when analysing simulated
                  lightcurve.
       seed     : initial seed number for monte carlo
                  simulation. Must be negative number when starting
                  new. Refer ran3 lib function.
                  Per default seed is obtained from /dev/random, if
                  this does not exist, it is set to -systime (seconds since
                  1970, jan 1.)
       chatty   : if 1 (or /chatty): report progress. If 2, report
                  single epfold progress also.
       numinst  : input lightcurve is normalized to one instrument
                  and was observed with numinst independent instruments
                  (currently only for Poisson case; default:1)

             fitchi   : fit a Gauss distribution to the chi^2 and use
                        the center of the gauss to determine the
                        period, instead of using the maximum of the
                        chi^2 distribution.
             debug    : Use rate as the estimated data set, i.e. rate
                        contains the actual period without noise. Can
                        be used to debug the profile multiplying
                        method.
             noabort  : keyboard has no influence on aborting

             returns: estimated error value for given period applied
                      to the input rate data set.

          maxperlist: (unsorted) list of periods found to be a
                      maximum using a simulated input lightcurve. Contains
                      ntrial entries.

             none !!!

             none

             The processing may last a long time (be prepared to wait
             hours to weeks). To abort the main monte carlo loop one
             may enter "x" on the keyboard unless the noabort keyword
             is set. Already computed values are returned in maxperlist.

             the input lightcurve has to be given in count rates (not
             in photon numbers).

             This routine tries to estimate the error of a previous
             received period using the epoche folding approach. For
             this we do:
             1.) calculate a mean profile with given period.
             2.) compute the intensity for all times applying the
                 period multiplied profile.
             3.) simulate an error for all times (standard: using
                 Poisson statistics, or use the error for normal
                 statistics).
             4.) Perform epoch folding for that simulated lightcurve.
             5.) Go to step 2.) Ntrial times, sum up the maximum of epoche
                 folding found.
             6.) Compute the standard deviation of the Ntrial maxima
                 obtained and take these as the error.

                 Davies, S.R., 1990, MNRAS 244, 93
                 Larsson, S., 1996, A&AS 117, 197
                 Leahy, D.A., Darbro, W., Elsner, R.F., et al., 1983,
                    ApJ 266, 160-170
                 Schwarzenberg-Czerny, A., 1989, MNRAS 241, 153

            epfold,time,rate,chierg=chierg,maxchierg=maxchierg,pstart=9,pstop=11
            print,"Period:", maxchierg[0],$
                  "Error:", epferror(time,rate,pstart=9,pstop=11,period=maxchierg[0])

 $Log: epferror.pro,v $
 Revision 1.10  2006/02/24 13:40:08  kreyken
 Added tolerance to the function call and promoted it to epfold.
 The tolerance was alreday used to determine the profile, but as it could not be set, it had no effect result in crappy simulted profiles!
 Revision 1.9  2004/10/26 14:03:51  goehler
 fix of raterr description.
 Revision 1.8  2004/09/06 22:56:22  wilms
 made compatibile with SGE batch jobs (which do not have a TTY)
 Revision 1.7  2004/09/06 22:22:26  wilms
 * corrected major bug in Poisson simulation to now correctly draw
  number of photons in a bin and only after this compute the
  simulated count rate.
 * Switch from ran3 to ran2 because of issues with ran3.
 * If present, read from /dev/random to initialize the RNGs (better
   granularity than using systime(1).
 * Added numinst keyword, used in the Poisson simulation only, to take
   care of case of lightcurves normalized to numinst instruments.
 * Further cleanups and optimizations
 Revision 1.6  2004/03/31 12:24:41  goehler
 implementation of beta state error computation for real event lists
 using the appropriate statistics. *this is very cpu consuming*.
 Revision 1.5  2004/03/01 16:30:29  goehler
 fix: allow event lists without rate info also. In this case the event time
      will directly passed to the pfold routine, and simulations are done
      for *all* event times given.
      This approach is not fully statistically funded  and should be
      replaced with a more sensible algorithm which is able to simulate
      events for given rates within the good time intervals.
 Revision 1.4  2003/11/04 13:01:01  rodina
 fix: in case of /fitchi option maxchierg[0] also contains the best
      period found (and *not* chierg[2,0]).
 Revision 1.3  2003/10/22 12:16:53  rodina
 Increased chatty period precission (e.g.)
 Revision 1.2  2003/10/11 19:46:31  goehler
 replaced random generator with ran3_poisson (misc library)
 Revision 1.1  2003/10/10 12:14:19  goehler
 initial, still buggy version

[Home Page] [Software, Documentation] [IDL Documentation] [Quick Reference] [Feedback]