;+
; NAME:
;       PLOT2AXIS
;
; PURPOSE:
;
;       The purpose of this procedure is to provide
;	a simple method of plotting 2 time series
; 	of data on a single graph.
;
; 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):
;
;       PLOT2AXIS, x1data_in, y1data_in: 
;		Depending upon how many parameters are specified,
;		x2data_in and y1data_in could represent several
;		different things. See below for details.
;
; OPTIONAL INPUT PARAMETERS:
;
;       x2data_in: If y2data_in is not specified, this will be the second dependent data set
;		corresponding to the independent data found in x1data_in. If y2data_in is
;		specified, this will be the dependent data for y2data_in.
;	y2data_in: If specified, will always be indepenent data for the dependent data set in x2data_in.
;
; 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.
;	color1: The color to be used for the first data set. This is a string.
;	color2: Same as color1, except for the second data set.
;	filename: If using /print, the file which to save. Defaults to 'a.eps'
;	title, xtitle: Same as for the PLOT procedure.
;	y1title, y2title: Similar to ytitle, except for use on the two different y axes.
;       badvalue: Used for interpolating the bad data points. Default value is -9999.
;	noclose, noopen: If using print, and you don't want to open or close the file - for whatever reason.
;	trend: Boolean keyword. If true (1), trend lines will be plotted on the graph.
;	psym, psym1, psym2: Same as for the PLOT procedure. psym1 corresponds to the first time series, and
;		psym2 to the second. Use psym to specify both at the same time.
;	smooth: If set to a value greater than one, applies the IDL smooth() function to the data sets.
;	sameaxis: Boolean keyword. If true (1), both of the y axes will have the same maximum and minimum.
;	noiter: The default is to interpolate bad data points. Set this keyword if you don't want this.
;	bad: The bad value to use for iterating bad data points.
;	window: Specify which window to plot to. Default is 0.
;
; OUTPUT KEYWORD PARAMETERS:
;       None.
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None known.
;
; RESTRICTIONS:
;
;   Required programs from my library:
;     interpolate_bad_points()
;
;   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
;
; IMPORTANT NOTE:
;	This procedure will not work as currently written unless you have
;	a function called numtostring() which takes a number (integer, float,
;	or double), and converts it into a string variable. This is a function
;	I wrote that I have not released yet. If you wish to use this procedure,
;	just replace all occurrences of "numtostring(" with string("; without
;	the quotes.
;
; EXAMPLE:
;
;      plot2axis, x, y, /print, filename='helloworld.ps'
;
;	This command will take the variables x and y, and
;	plot them on the two different y axes, using the 
;	same x axis. It will print it to a file called
;	helloworld.ps.
;
;
; MODIFICATION HISTORY:
;
;       Written by: Nathan Johnson, 17 January 2008.
;-

pro plot2axis, x1data_in, y1data_in, x2data_in, y2data_in, xtitle=xtitle, y1title=y1title, y2title=y2title, $
               filename=filename, title=title, range=range, print=print, noclose=noclose, $
               y2range=y2range, sameaxis=sameaxis, smooth=smooth, noopen=noopen, xstyle=xstyle, $
               trend=trend, color1=color1, color2=color2, noiter=noiter,  bad=bad, psym=psym, $
               p1sym=psym1, p2sym=psym2, window=window

  on_error, 0

; Remember this. To be used later.
  thisdevice = !d.name
  device, get_decomposed = get_decomposed

; Preserve the data that is passed to this procedure.
  IF KEYWORD_SET(x1data_in) THEN x1data = x1data_in
  IF KEYWORD_SET(y1data_in) THEN y1data = y1data_in
  IF KEYWORD_SET(x2data_in) THEN x2data = x2data_in
  IF KEYWORD_SET(y2data_in) THEN y2data = y2data_in

; This procedure can accept from 2 to 4 parameters.
; If 2 parameters are specified, the first is taken to be
; the first dependent variable, and the second is taken to
; be the second dependent variable. They are both given
; the same number of independent values as there are
; dependent values.

  if N_PARAMS() eq 2 then begin
    if not keyword_set(range) then range = [0, N_ELEMENTS(x1data)]
    p_x1 = indgen(N_ELEMENTS(x1data))
    p_x2 = indgen(N_ELEMENTS(y1data))
    p_y1 = x1data
    p_y2 = y1data
  endif

; If 4 parameters are specified, the first is taken to be the
; first independent variable, the second is the first dependent
; variable, the third is the second independent variable, and 
; the fourth is the second dependent variable.

  if N_PARAMS() eq 4 then begin
      if not keyword_set(range) then range = [min(x1data), max(x1data)]
      p_x1 = x1data
      p_x2 = x2data
      p_y1 = y1data
      p_y2 = y2data
  endif

; If 3 parameters are specified, use the first for the independent
; variable for each of the next two parameters, which are used as
; the dependent variables. First make sure that the data are the
; same length. Chop off the ends of the dependent variables if they
; are longer than the dependent.

  if N_PARAMS() eq 3 then begin
    if not keyword_set(range) then range = [min(x1data), max(x1data)]
    IF N_ELEMENTS(y1data) gt N_ELEMENTS(x1data) THEN $
      y1data = y1data[0:N_ELEMENTS(x1data)-1] ELSE $
      x1data = x1data[0:N_ELEMENTS(y1data)-1]

    IF N_ELEMENTS(x2data) gt N_ELEMENTS(x1data) THEN $
      x2data = x2data[0:N_ELEMENTS(x1data)-1] ELSE $
      x1data = x1data[0:N_ELEMENTS(y1data)-1]

    p_x1 = x1data
    p_x2 = x1data
    p_y1 = y1data
    p_y2 = x2data
  endif

; Set some variables if they are not specified by the user.
  if not keyword_set(noopen) then noopen=0 else noopen = 1
  if not keyword_set(smooth) then smooth=0
  if not keyword_set(sameaxis) then sameaxis=0 else sameaxis=1
  if not keyword_set(doclose) then doclose=0 else doclose=1
  if not keyword_set(trend) then trend=0 else trend=1
  if not keyword_set(color1) then color1 = "navy"
  if not keyword_set(color2) then color2 = "orange"
  if not keyword_set(bad) then noiter=1
  if not keyword_set(filename) then fname='a.eps'
  if not keyword_set(window) then window = 0

; Check to make sure color is a string
  varInfo = Size(color1)
    IF varInfo[varInfo[0] + 1] NE 7 THEN color1="navy"
  varInfo = Size(color2)
    IF varInfo[varInfo[0] + 1] NE 7 THEN color2="orange"

; Default is to do a linear interpolation to find
; values for bad points. Use the keyword /noiter to
; disable interpolation. 
  if not keyword_set(noiter) then begin
    p_x1 = interpolate_bad_points(p_x1, bad)
    p_x2 = interpolate_bad_points(p_x2, bad)
    p_y1 = interpolate_bad_points(p_y1, bad)
    p_y2 = interpolate_bad_points(p_y2, bad)
  endif

; Set the appropriate symbol for the plots if necessary.
  if not keyword_set(psym1) then $
    if keyword_set(psym) then psym1 = psym

  if not keyword_set(psym2) then $
    if keyword_set(psym) then psym2 = psym

; Find the appropriate range, if not set by the user. 
  if not keyword_set(y2range) then begin
     y2min = MIN(smooth(p_y2,smooth), MAX=y2max, /NAN) 
  endif else begin
     y2min = y2range[0]
     y2max = y2range[1]
  endelse

  y1min=MIN(smooth(p_y1,smooth), MAX=y1max,/NAN)

; If the user wants the 2 y-axes to be the same, figure out
; the values to use.
  IF sameaxis eq 1 then begin
    ymax = max([y1max, y2max], /NAN)
    ymin = min([y1min, y2min], /NAN)
    y1max=ymax
    y2max=ymax
    y1min=ymin
    y2min=ymin
  ENDIF

; If the arrays are of equal length, find the correlation between them.
; To be later used as the subtitle of the plot.
  IF Array_Equal(p_x1, p_x2) eq 1 THEN BEGIN
    corr = real_part(correlate(smooth(p_y1,smooth), smooth(p_y2,smooth)))
    if abs(corr) lt 0.001 then corr=0.d
    if corr ge 0 then subtitle='R: '+strmid(numtostring(corr), 0, 5) $
    else subtitle='R: '+strmid(numtostring(corr), 0, 6)       ; 3 sig figs for both + and - values
  ENDIF

; Set the device based upon the /PRINT keyword
  if keyword_set(print) then begin
    set_plot, 'ps'
    !p.font=0
    if not keyword_set(noopen) then 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, window, xsize=xsize, ysize=ysize
  endelse

; Use fsc_color() to get the numerical color based on the current device.
  color1 = fsc_color(color1)
  color2 = fsc_color(color2)

; Plot the data, omit left and right axes: .
  PLOT, p_x1, smooth(p_y1, smooth), TITLE = title, SUBTITLE = subtitle, /YNOZERO, XTITLE = xtitle, YTITLE = y1title, $
    YSTYLE=5, XRANGE=[range[0],range[1]], YRANGE=[y1min, y1max],charsize = 1.5, XSTYLE=xstyle, XMARGIN=[8, 8], YMARGIN=[4, 4], $
    psym=psym1, background=fsc_color('white'), color=fsc_color('black')

; Draw the left axis and data in the specified color.
  AXIS, YAXIS=0, YRANGE =[y1min,y1max], YSTYLE = 1, YTITLE = y1title, /SAVE, color=color1, charsize = 1.5
  OPLOT, p_x1, smooth(p_y1, smooth), color=color1, psym=psym1

; Add a trendline, if desired.
  if trend eq 1 then begin
    trend1 = linfit(p_x1, p_y1)
    oplot, [min(p_x1),max(p_x1)],[trend1[0]+trend1[1]*min(p_x1),trend1[0]+trend1[1]*max(p_x1)], linestyle=2, color=color1
  endif
  
; Draw the right axis and data in the specified color.  
  AXIS, YAXIS=1, YRANGE =[y2min,y2max], YSTYLE = 1, YTITLE = y2title, /SAVE, color=color2, charsize = 1.5
  OPLOT, p_x2, smooth(p_y2, smooth), color=color2, psym=psym2

; Add a trendline, if desired.
  if trend eq 1 then begin
    trend2 = linfit(p_x2, p_y2)
    oplot, [min(p_x2),max(p_x2)],[trend2[0]+trend2[1]*min(p_x2),trend2[0]+trend2[1]*max(p_x2)], linestyle=2, color=color2
  endif

; 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) AND not keyword_set(noclose) then device, /close
  set_plot, thisdevice
  device, decomposed=get_decomposed

END