pro ipsplot, geomagnetic=geomagnetic,position=position

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	IPSPLOT
; PURPOSE:
;	Plot time series for IMP and/or G-IPS data for specified Carrington
;	rotation
; CALLING SEQUENCE:
;	IPS_PLOT, /geomagnetic
; INPUTS:
;	User prompts. EGIPSY output files located in X$DATA directory
; OUTPUTS:
;	Plot to screen and/or printer
; CALLS:
; SIDE EFFECTS:
;	The time axis runs from right to left
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;	September 1992, Kari Winfield & Susan Rappoport (UCSD)
;	JUN-1994, Paul Hick (UCSD), rewrite
;-

if keyword_set(geomagnetic) then	$
    ips_geo, title, titles, tags, data,	$
	logpos, logshift, invalid, ranges, gparm, nparm, yparm, y	$
else					$
    ips_imp, title, titles, tags, data,	$
	logpos, logshift, invalid, ranges, gparm, nparm, yparm, y

elongstr = ['30-50','30-80','50-70']

rot = 0  &  echo,'Carrington Rotation?', rot

print, '          Elongation (Inner-Outer)?: '+	$
	'(0='+elongstr(0)+		$
	',1='+elongstr(1)+		$
	',2='+elongstr(2)+')'
l = 1  &  echo,' ',l,0,2
elstr = elongstr(l)

l = n_elements(tags)
print, '          ', data
echo, 'First time series', yparm, 0,l-1

print, '          ', data+','+strcompress(l,/rem)+'=None'
echo, 'Second time series', y, 0,l

if y ne yparm and y ne l then yparm = [yparm,y]

ylog = logpos(yparm)				; = 1 if log has to be taken
rparm = yparm-ylog*logshift(yparm)		; Shift to argument of log

dayparm = -8  &  echo,'Day: Aft,Fore,Current(-8 to 8)',dayparm,-8,8

MJDtoJD = 2400000.5D0

carrtime, file='X$AUX:GIPSY.CAR', rot, [0.,360.],jd=jdrot

i5 = 5
julian, 1, iyr, doy5, jdrot(1)
doy5 = int(doy5)/i5*i5
julian, 0, iyr, doy5, jd5

jd5 = jd5+i5*indgen(10)
jd5 = jd5( where(jdrot(1) lt jd5 and jd5 lt jdrot(0)) )
julian, 1, iyr, doy5, jd5		; Label time axis every 5 days
doy5 = round(doy5)
jd5 = jd5-MJDtoJD

frst = long( 1+long(jdrot(1)) -MJDtoJD)	; MJD containing first noon in rotation
last = long( long(jdrot(0)) -MJDtoJD)	; MJD containing last noon in rotation
xrange = jdrot-MJDtoJD			; MJD time of 0 and 360 degree passage

span = last-frst+1			; # MJDs to be read
tt = frst+indgen(span)			; List of (long int) MJDs to be read
yy = fltarr(span, n_elements(yparm) )	; Time series

for i = frst,last do begin
					; Read file offset by -`dayparm' days
	if flt_read( filepath(root=getenv('DAT'),strcompress(i-dayparm,/rem)+'.dat'),f, /silent) then begin

		; fc(*,i) = [day,doy,gval,delg,dens,deldens,vel,delvel]
		; day = offset -8,-6,..0,..,6,8
		; doy = day of year
		; Get data offset by `dayparm' days
		f = f(*, where( f(0,*) eq dayparm ) )
		yy( where( tt eq i) ,*) = f(2+rparm)
	endif else	$
		yy( where( tt eq i) ,*) = invalid(rparm)

endfor

n = where( yy(*,0) ne invalid(rparm(0)))	; Identify valid data
if n(0) eq -1 then message, 'First time series must contain valid data

subtitle = 'Starting at MJD:'+strcompress(frst,/rem)+	$
	'  Elongation:'+elstr+' Degrees'+	$
	'  Day:'+strcompress(dayparm,/rem)+		$
	'  Carrington Rotation:'+strcompress(rot,/rem)

if not keyword_set(position) then position = [0.1,0.15, 0.91, 0.93]

set_plot, !TheTerminal
outpr = 0

AGAIN:

y = yy(*,0)
n = where( y ne invalid(rparm(0)))

t = tt(n)+.5				; All data refer to noon
y = y(n)
if ylog(0) then y = alog10(y)

if n_elements(yparm) eq 1 then	$
    ystyle = 0			$	; One time series: both y-axes
else				$
    ystyle = 8				; Two time series: only left axis

yrange = [min(y),max(y)]

plot,	xrange,yrange,		$	; Plot bottom x-axis and left y-axis
	/nodata,		$
	/ynozero,		$
	xrange = xrange,	$
;	yrange = yrange,	$
	position = position,	$
	xstyle = 1,		$	; 1=Exact scaling, 8=bottom axis only
	ystyle = ystyle,	$
	xtitle = 'DOY',		$
	ytitle = titles(yparm(0)),	$
	subtitle = subtitle,		$
	xticks = n_elements(jd5)-1,	$
	xtickv = jd5,			$
	xtickname = strcompress(doy5,/rem) ; Day of year labels for x-axis

;axis,	xaxis  = 1,	$		; Plot top axis (same as bottom axis)
;	xstyle = 1,	$
;	xtickname = strcompress(doy5,/rem) ; Separate axis call to get labels

; Plot time series 1 (connect only data points which are one day apart)

n = n_elements(n)
m = 0
while m lt n do begin
    d = where ( t(m:n-1) eq t(m)+indgen(n-m), i )
    oplot, [t(m+d)],[y(m+d)],psym = -2
    m = m+i
endwhile
xyouts, t(n-1)+1, y(n-1), tags(yparm(0))	; Tag last data points

if n_elements(yparm) eq 1 then goto, SKIP

; Plot right y-axis for the second time series.
; Save the user scale for the subsequent `oplot' calls

y = yy(*,1)
n = where( y ne invalid(rparm(1)))
if n(0) eq -1 then goto, SKIP

t = tt(n)+.5					; All data refer to noon
y = y(n)
if ylog(1) then y = alog10(y)

axis,	yaxis  =1,			$
;	ystyle =1,			$
	yrange = [min(y), max(y)],	$	; Needed to change the y-range
	ytitle = titles(yparm(1)),	$
	/save					; Save the user scaling

; Plot time series 2 (connect only data points which are one day apart)

n = n_elements(n)
m = 0
while m lt n do begin
    d = where ( t(m:n-1) eq t(m)+indgen(n-m) , i )
    oplot, [t(m+d)],[y(m+d)],psym = -6
    m = m+i
endwhile
xyouts, t(n-1)+1, y(n-1), tags(yparm(1))	; Tag last data point

SKIP:

if outpr then begin				; Hardcopy
    spitplot
    set_plot, !TheTerminal
endif else begin
    echo, 'Hardcopy', outpr, /log
    if outpr then begin
	set_plot, !ThePrinter
	goto, again
    endif
endelse

twin,/hide
return &  end