gdl-astrolib 2018.02.16+dfsg-1 / usr / share / gnudatalanguage / astrolib / pent.pro

This file is indexed.

/usr/share/gnudatalanguage/astrolib/pent.pro is in gdl-astrolib 2018.02.16+dfsg-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below. 

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
 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 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
;-

 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

APT Browse - Built by Thomas Orozco.