;+
; NAME:
;       CUMULATIVE_CORRELATION
;
; PURPOSE:
;
;       The purpose of this procedure is to provide
;	a graphical method of determining correlation
;	or covariance between two time-series of data.
;	Unlike a 'traditional' FFT, this method adds
;	each correlation or covariance so that the
;	areas where the steepest slope occur is where
;	there is maximum correlation/covariance.
;
; AUTHOR:
;
;       Nathan Johnson
;       University of Arizona
;       Department of Atmospheric Sciences
;       johnson AT atmo DOT arizona DOT edu
;
; CATEGORY:
;
;       Graphics
;
; CALLING SEQUENCE (WITH REQUIRED INPUT PARAMETERS):
;
;       CUMULATIVE_CORRELATION, var1_in, var2_in: 
;                   var1_in is a time series of data
;                   var2_in is a time series of data
;
;	var1_in and var2_in should be the same length.
;	If they aren't, the longer one will be truncated.
;
; OPTIONAL INPUT PARAMETERS:
;
;       mean_velocity: if this parameter is specified, instead of displaying the result
;		       as a function of time, it is displayed as a function of distance
;
; INPUT KEYWORD PARAMETERS:
;
;	print: Boolean keyword. If true (1), then the result will plot on the Postscript device.
;	range: The subscript range over which to plot. Defaults to the whole range.
;	oplot: Boolean keyword. If true (1), then the result is plotted over the current plot using OPLOT.
;	color: The color to use for the plot. The background is currently always white, and the axes are black.
;		This defaults to red if using the oplot keyword, otherwise it defaults to green. Note: this
;		is the integer index of the color, and not a string.
;	covariance: Boolean keyword. If true (1), then the result shows the cumulative covariance, not correlation.
;	filename: If using /print, the file which to save. Defaults to 'a.eps'
;	title, xtitle, ytitle : Same as for the PLOT procedure.
;       xscale: The amount to scale the x-axis by. Useful for unit conversions (ie. m -> km)
;       badvalue: Used for interpolating the bad data points. Default value is -9999.
;
; OUTPUT KEYWORD PARAMETERS:
;       None.
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None known.
;
; RESTRICTIONS:
;
;   Required programs from the Coyote Library:
;   (Written by David Fanning)
;
;      http://www.dfanning.com/programs/error_message.pro
;      http://www.dfanning.com/programs/pickcolorname.pro
;      http://www.dfanning.com/programs/fsc_color.pro
;
;   I suggest downloading the entire library from:
;      http://www.dfanning.com/documents/programs.html
;
; EXAMPLE:
;
;      cumulative_correlation, x, y, 100, /covariance
;
;	This command will find the cumulative covariance between
;	the two time series x and y, figure out the corresponding
;	x-axis (since a mean_velocity was specified, a spatial
;	axis), and plot it to the screen.
;
;
; MODIFICATION HISTORY:
;
;       Written by: Nathan Johnson, 16 January 2008.
;-

PRO CUMULATIVE_CORRELATION, var1_in, var2_in, mean_velocity, print=print, range=range, oplot=oplot, $
     color=color, covariance=covariance, filename=filename, title=title, xtitle=xtitle, ytitle=ytitle, $
     xscale=xscale, badvalue=badvalue

  on_error, 2

  IF N_PARAMS() LT 2 THEN MESSAGE, 'This procedure requires at least 2 parameters.'
  IF N_PARAMS() LT 3 THEN no_v = 1 ELSE no_v = 0

; Check that the two input time series are the same length. If not, truncate the longest.
  var1 = var1_in
  var2 = var2_in

  N1 = N_ELEMENTS(var1)
  N2 = N_ELEMENTS(var2)

  IF N1 GT N2 THEN var1 = var1[0:N2-1]
  IF N2 GT N1 THEN var2 = var2[0:N1-1]

; Remember current device
  thisdevice = !d.name
  device, get_decomposed = get_decomposed

; Set some variables if the user didn't specify them as keywords.
  if not keyword_set(xscale) then xscale = 1.0
  xscale = double(xscale)
  if not keyword_set(filename) then filename = 'a.eps'
  if not keyword_set(range) then range=[1, N_ELEMENTS(var1)/2]
  if range[0] lt 1 then range[0] = 1
  if range[1] gt N_ELEMENTS(var1)/ 2 then range[1]=N_ELEMENTS(var1)/2
  if keyword_set(mean_velocity) then mean_velocity = double(mean_velocity)
  if not keyword_set(badvalue) then badvalue = -9999

; Look for bad points and interpolate them. Otherwise you'll get 'bad' answers.
; Linear interpolation may not be the best solution, but it'll do for now.

  var1 = interpolate_bad_points(var1, badvalue)
  var2 = interpolate_bad_points(var2, badvalue)


; Figure out the x-axis.
  data_x = dindgen(n_elements(var1)+1)
  n_samples = n_elements(var1)
  n_nu_points = n_samples/2 + 1
  x_range = max(data_x) - min(data_x)

; Temporal x-axis. If user specified a mean_velocity, then the x-axis
; will be changed later to use a spacial domain.
  nu = dindgen(n_nu_points)/(n_nu_points-1) / (x_range / n_samples) / 2	

; Remove the means before doing the correlation/covariance.
  x=fft(var1 - mean(var1))
  y=conj(fft(var2 - mean(var2)))
  z=x*y
  z=z(0:N_ELEMENTS(Z)/2+1)

; Which do you want, correlation or covariance?
  IF KEYWORD_SET(covariance) THEN BEGIN
    a=total(z, /cumulative) 
    IF NOT KEYWORD_SET(ytitle) THEN ytitle='Cumulative Covariance'
  ENDIF ELSE BEGIN
    a=2.*total(z, /cumulative)/(stdev(var1)*stdev(var2))
    IF NOT KEYWORD_SET(ytitle) THEN ytitle='Cumulative Correlation'
  ENDELSE

; Set the device based upon the /PRINT keyword
  if keyword_set(print) then begin
    set_plot, 'ps'
    !p.font=0
    device, filename=filename, /color
  endif else begin
    if thisdevice eq 'PS' then begin
      ; this will only work on linux (including new macs) and windows machines!
      if !version.os eq 'linux' then set_plot, 'X' else set_plot, 'win'
    endif else begin
      set_plot, thisdevice
    endelse
    !p.font=-1
    device, decomposed=0
    window, xsize=xsize, ysize=ysize
  endelse

; If the user didn't specify a color, choose one.
  IF NOT KEYWORD_SET(color) THEN $
    IF KEYWORD_SET(oplot) THEN color=fsc_color('red') ELSE color=fsc_color('green') 

; If a mean_velocity is specified, make the x-axis into a spatial domain,
; otherwise, it will be a temporal domain.
  IF no_v eq 1 THEN BEGIN
    x_axis = nu[1:*] * xscale
    IF NOT KEYWORD_SET(xtitle) then xtitle = 'Temporal Scale'
  ENDIF ELSE BEGIN
    x_axis = mean_velocity/nu[1:*] * xscale
    IF NOT KEYWORD_SET(xtitle) then xtitle = 'Spatial Scale'
  ENDELSE

; Plot the result to the device specified.
  IF NOT KEYWORD_SET(oplot) THEN $
    PLOT, x_axis, a, /YNOZERO, charsize = 1.5, /xlog, background=fsc_color('white'), $
    color=fsc_color('black'), title=title, xtitle=xtitle, ytitle=ytitle

  OPLOT, x_axis, a[range[0]-1:range[1]], color=color

; Clean up. Close the postscript device if necessary, and set the plot back to 
; whatever it was set to before this procedure was called.
  if keyword_set(print) then device, /close
  set_plot, thisdevice
  device, decomposed=get_decomposed

END