idlastro / Math and Statistics: PENT

[Source code]

NAME
PENT
PURPOSE
Return the information entropy of a time series
EXPLANATION
This function will return S, the information entropy of a time series
for a set of trial periods 
CATEGORY
Time series analysis, period finding, astronomical utilities.
CALLING SEQUENCE
Result = PENT(P, T, X, [N, M ] )
INPUTS
P - array of trial period values.
T - array of observation times (same units as P).
X - array of observations.
OPTIONAL INPUTS
N   - If  four parameters are given then the 4th parameter is assumed
        to be N. Then NxN boxes are used to calculate S.
M,N - If five parameters are given then parameter 4 is M and parameter
        5 is N. S is then calculated using MxN boxes - M partitions for the
        phase and N partitions for the data.
OUTPUTS
This function returns S, the information entropy of the time series for
the periods given in P as defined by Cincotta, Me'ndez & Nu'n~ez
(Astrophysical Journal 449, 231-235, 1995). The minima of S occur at
values of P where X shows periodicity.
PROCEDURE
The procedure involves dividing the phase space into N^2 partitions 
(NxN boxes) and then calculating:
        __ N^2
  S = - \        mu_i . ln(mu_i)  for all mu_i <> 0
        /_  
          i = 1 
where  mu_i is the number of data points in partition i normalised by 
the number of partitions.
The option of using MxN boxes is an additional feature of this routine.
EXAMPLE
To generate a similar synthetic data set to Cincotta et al. we
 do the following:
IDL> P0 = 173.015                        ; Fundamental period
IDL> T = randomu(seed,400)*15000         ; 400 random observation times
IDL> A0 = 14.0                           ; Mean magnitude
IDL> M0 = -0.5  * sin(2*!pi*T/P0)        ; Fundamental mode
IDL> M1 = -0.15 * sin(4*!pi*T/P0)        ; 1st harmonic
IDL> M2 = -0.05 * sin(6*!pi*T/P0)        ; 2nd harmonic
IDL> sig = randomu(seed,400)*0.03        ; noise
IDL> U = A0 + M0 + M1 + M2 + sig         ; Synthetic data
IDL> Ptest = 100. + findgen(2000)/2.     ; Trial periods 
IDL> S = pent(Ptest,T,U)                 ; Calculate S
        ... this takes a few seconds ...
IDL> plot,Ptest,S,xtitle="P",ytitle="S"  ; plot S v. P
IDL> print,Ptest(where(S eq min(S)))     ; Print best period (+/- 0.5)
The plot produced should be similar to Fig. 2 of Cincotta et al.
RESTRICTIONS
My own (limited) experience with this routine suggests that it is not
as good as other techniques for finding  weak,  multi-periodic signals in 
poorly sampled  data, but is good for establishing periods of eclipsing
binary stars when M is quite large (try MxN = 64x16, 128x16 or even 
256x16).  This suggests it may be good for other periodic light curves 
(Cepheids, RR Lyrae etc.).
I would be glad to receive reports of other peoples experience with
this technique (e-mail pflm@bro730.astro.ku.dk).
MODIFICATION HISTORY
Written by:   Pierre Maxted, 14Sep95
Modifications:
Normalisation of S corrected, T-min(T) taken out of loop.
        -  Pierre Maxted, 15Sep95
Converted to IDL V5.0   W. Landsman   September 1997