pent.pro

basic median filter simulation

 function pent,p,t,x,m,n
;+
; 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 sythetic 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 [e-mail]).
;
; 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
;-

 on_error,2 ; return to caller

; Check suitable no. of parameters have been entered.

  case N_params() of 
   3 : begin
        n = 8.0 
        m = 8.0
       end
   4 : begin
        n = float(fix(m)) 
        m = n
       end
   5 : begin
        m = float(fix(m)) 
        n = float(fix(n))
       end
  else : message,/noname,' Syntax - Result = ( P, T, X [ [,M ] ,N ])'
  endcase

  nbox = m*n
  np = n_elements(p)
  npts = n_elements(x)

  if n_elements(t) ne  npts  then message , $
     'Input arrays T and X  must have same number of elements'

  if npts lt 3 then message,' Insufficient data in input arrays'

  npts = float(npts)

  S = fltarr(np)
 
   norm = (X - min(X))/(max(x) - min(x))   ; normalised data
   norm = norm - (norm eq 1.0)*(0.1/n) ; norm = 1 -> norm = 0.99..
   ni = 1 + n*(floor(norm*n))

   Tplus = T-min(T)  ; take this operation out of the loop
 
  for j = 0l,np - 1l do begin
 
   phi = ( Tplus / P[j] ) mod  1.0
 
   mu = histogram(floor(phi*m) + ni,max=nbox,min=0.0)/(npts)
 
   mu = mu[where(mu gt 0.0)]
   S[j] = -total(mu*alog(mu))
   
  endfor
  
  S = S/alog(nbox) ; normalise S

  return,S

end   ; That's all folks