pro varadi,file,rmax,M,K,filters,seed,S_lag,eigen_values
;+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; NAME:
;	varadi
; CALLING SEQUENCE:
;	varadi,file,rmax,M,K,filters,seed,S_lag,eigen_values
; PURPOSE:
;	Compute the RLSSA from Varadi et al, ApJ, 1999, 526, 1052
; INPUTS:
;	file		filename
;	rmax		number of largest eigenvalues for reconstruction
;	M		number of random lags
;	K		maximum lag
;	filters		reconstructed filters of size fltarr(K)
;	seed		dummy argument
;	S_lag		return set of lags fltarr(M)
;	eigen_values	return eigenvalues fltarr(M)
; OPTIONAL KEYWORDS:
; 	none
; OUTPUTS:
;	none
; OPTIONAL OUTPUT:
;	none
; EXAMPLE:
;	thomson,file
; LIMITATIONS:
;	None yet
; COMMONS:
;	None
; PROCEDURES USED:
;       fft, readfits
; MODIFICATION HISTORY:
;	Beta 1: WRITTEN, Thierry Appourchaux, November 14, 2000
;	Beta 2: Return also eigenvalues, November 17, 2000
;------------------------------------------------------------------------------
;
; Read file
	print,'r_max (10 is a good number)=',rmax
	print,'M (number of different lags, 250 is OK)',M
	print,'K (maximum lag, 25000 is OK)',K
	a=readfits(file,/silent)
;	a=randomn(seed,131072)+cos(2.*!pi*findgen(131072)/3.)
	m_2=moment(a,/double)
	a=(a-m_2(0))/sqrt(m_2(1))
;	a=randomn(seed,131072)*1.d0
;
;
	N_a=N_elements(a)
;
	k_lag=K
;
; Here we assume that the time series is sampled at 60 s
;
	xfreq=findgen(N_a)*1.e06/60./N_a
	a_rev=reverse(a)
;
;	Pad the time series with zero
;
	a_new=dblarr(2*N_a)
	a_new_rev=dblarr(2*N_a)
	a_new(0:N_a-1)=a
	a_new_rev(0:N_a-1)=a_rev
;
;	Then compute the correlation using fft (55)
;	Please note the use of the different normalization
;	It was calibrated using the real formula (53)
;

	q=2.d0*N_a*fft(fft(a_new,-1,/double)*fft(a_new_rev,-1,/double),1,/double)
;
; Shift to put the max at 0
;
	q=shift(double(q),-(N_a-1))
	q=q(0:k_lag-1)
	print,q(0:10)
;
; Calibration
;
;
;	q_clas=fltarr(k_lag)
;	for k=0,k_lag-1 do begin
;		q_clas(k)=total(a(0:N_a-1-k)*a(k:N_a-1))
;	endfor
;
; Equation (57)
;
	q=q/double(N_a-findgen(k_lag))
;
; normalize, equation (58)
;
	c=q/q(0)
;
; Set randomly the lag and makes sure that they are all different
;
	C_diff=dblarr(M,M)
	S=floor(k_lag*randomu(seed,M))+1
	for i=0,M-1 do begin
		for j=i,M-1 do begin
			;print,i,j
			C_diff(i,j)=abs(S(i)-S(j))
			if (i ne j) then begin
				repeat begin
					nodiff=1
					if (C_diff(i,j) eq 0.) then begin
						S(j)=floor(k_lag*randomu(seed))+1
						C_diff(i,j)=abs(S(i)-S(j))
						;print,C_diff(i,j),i,j
						nodiff=0
					endif
				endrep until nodiff
			endif
		endfor
	endfor
;	print,S
	S_lag=S
;
; Set the random matrix
;
	C_rand=dblarr(M,M)
	for i=0,M-1 do begin
		for j=i,M-1 do begin
			C_rand(i,j)=c(C_diff(i,j))
			C_rand(j,i)=C_rand(i,j)
		endfor
	endfor
	C_buffer=C_rand
	Print,'Eigenvalues...'
;
; Compute eigenvalues and eigenvectors
;
	evecs=1.d0
	eigenvalues=eigenql(C_buffer,eigenvectors=evecs,/double)
;
; Set the number of eigenvalues
;
	r=rmax
;
; Set filter vector to 0
;
	f1=dblarr(k_lag)
;
; Compute filter
;
	Print,'Filter...'
	for mm=0,r-1 do begin
		for i=0,M-1 do begin
			nn=S(i)-S
			index=where(nn ge 0.)
			f1(nn(index))=f1(nn(index))+evecs(i,mm)*evecs(index,mm)/M
		endfor
	endfor
	f1(0)=0.
	filters=f1
	eigen_values=eigenvalues
end