;The comusiclag function obtains the Cospectrum and its corresponding Ogive
;for two time series where one may lag the other by some number of points. 
;No special windows are used. Use is as follows:
;       comusiclag,xin,yin,lag,fs,avg,dtrnd,cxy,og,f,psisq
;  WHERE
;          |xin & yin = variable names of data vectors
;          |      lag = amount y lags x in points.
;    input |       fs = sampling frequency of data (hz)
;          |      avg = 1 to subt. mean from data before calc., 0 otherwise
;          |    dtrnd = 1 to detrend input data, 0 otherwise
;
;          |      cxy = cospectrum
;   output |       og = ogive function
;          |        f = frequency vector corresponding to cxy and og
;          |    psisq = mean sq. if mean removed; =covariance or flux value
;                       Note: if set~=1 psisq is avg of psisq from each set
;Written by BWBerger 2/1/00 based heavily on comusiclag.m for matlab.
;Reference: Bendat & Piersol as well as D.E. Newland.
;******************************************************************************

pro comusiclag,xin,yin,lag,fs,avg,dtrnd,cxy,og,f,psisq

x=xin
y=yin

if dtrnd eq 1 then begin
   print,'Detrending selected'
   detren,x,x,xdmn,xtren,-999.
   detren,y,y,ydmn,ytren,-999.
endif

x=x(0:len(x)-1-lag)
y=y(lag:len(y)-1)

; Remove mean from time series if requested.
if avg eq 1 then begin
   print,'Remove mean selected'
   x=x-mean(x)
   y=y-mean(y)
endif

if len(x) ne len(y) then begin
   print,'Input vectors not the same size.'
   goto,fini
endif

; Put data vectors into column format.
rw=row(x)
cl=col(x)

if cl gt rw then begin
   x=transpose(x)
end

rw=row(y)
cl=col(y)

if cl gt rw then begin
   y=transpose(y)
   rw=row(y)
   cl=col(y)
end

q=rw

ev=2*floor(q/2);
if q gt ev then begin
   q=q-1
endif
n=q

;Do fft as requested.
a=fft(x(0:q-1),/double)
b=fft(y(0:q-1),/double)

s=double(conj(a)*b)

t=(n-1)/double(fs)
s=2*t*s
cxy=float(s)

;Want only half of cospectrum (remove frequency values larger than Nyquist).
half=n/2+1
cxy=cxy(0:half-2)
f=findgen(1,half-1)*fs/n

; Find mean square value (assuming avg=1), and calculate ogive function.
psisq=(cxy(0)/(t*2))+(1/t)*total(cxy(1:half-2))
if avg ne 1 then begin
  print,'You did not remove mean from time series. psisq may ~= mean square.'
endif

og=fltarr(1,len(cxy))
og(0)=psisq
og(1)=og(0)-cxy(0)*(1/t)*(1/2)
for i=2,half-2 do begin
   og(i)=og(i-1)-cxy(i-1)*(1/t)
endfor

fini:
return
end