pro roll_plot,duration=duration,endtime,outfil=outfil,average=av,orbit=orbit,calc=calc,p=p,ecliptic=ecliptic,outplot=a,plotsize=plotsize,debug=debug,yrange=yrange,despike=despike,nsigma=nsigma,smooth=sm
;+
; NAME:
;	ROLL_PLOT
; 
; PURPOSE:
;	Generate a plot of the SOHO roll angle versus time
;
; CALLING SEQUENCE:
;	ROLL_PLOT, [duration], [endtime], [keywords=keywords]
;	
; INPUTS:
; 	ENDTIME:  Time the plot should *END* (default is Right NOW!)
;
; KEYWORD INPUTS:
;	DURATION: Duration of the plot in hours (default is 24)
;	OUTFILE: If specified, the output is saved here.
;	AV:	If set, averaging is done.  If set to 
;		a value greater than 5, averaging is done over that
;		many minutes (values less than five cause five-minute
;		averaging).
;	SMOOTH:	Synonym to AV.
;	ORBIT:	If set, use only minutes for which orbit data were available
;	CALC:	If set, use only minutes for which calculated orbit was used.
;
;	P:	If set, report solar "P" angle (this is also the default)
; 	ECLIPTIC: If set, report angle from ecliptic north
;
;	YRANGE: If set, specify the output Y range in degrees  (otherwise,
;		ROLL_PLOT guesses, using an algorith that loses some outlying
;		points in order to see small variations in the roll).
;
;	DESPIKE: If set, throw out all the points that look "spiky" -- ie
;		those that exceed a threshhold number of standard deviations
;		from the average from both their neighbors.
;	
;	NSIGMA: Number of standard deviations to use in autoscaling and in
;		despike threshholding (default is 3).
;	
; KEYWORD OUTPUTS:
;	OUTPLOT : If set, gets the plot itself as a byte array.
;
; HISTORY:
;	Written 20-May-99 CED
;	Added time-of-plot 20-May-99 CED
;	Fixed missing-day contingency, 20-May-99 CED
;	Fixed plotsize bug, 21-May-99 CED
; 	Added dual-trace plotting for times when both orbital and calculated
;		rolls are listed in the file (to avoid unsightly jumps!) 
;		8-Jun-99 CED
;	Added std-deviation-based scaling and despiking, 2-Aug-99, CED
;	Added rejection of outlying dates; added SMOOTH synonymo for AV.
;	Added simple singularity avoidance near the -180/180 branch jump. CED.
;       Modified handling of zbuffer to speed up display/output. 3-Dec-99 RIB.
;-
if keyword_set(orbit) then teststr='(ORBIT)'
if keyword_set(calc) then teststr='(CALC)'

; Set the index column to plot...
tfile_index = 2
if keyword_set(ecliptic) then tfile_index = 3
if keyword_set(p) then tfile_index = 2
columns = ' '+['Day','Time','Apparent P angle','Ecliptic angle','Method']


if not isvalid(endtime) then endtime = anytim(ut_time())
if not isvalid(duration) then duration = 24
if not isvalid(nsigma) then nsigma = 3

startime = endtime-duration*3600.0

if not isvalid(plotsize) then begin
	if !D.name eq 'X' then begin
		a = tvrd()
		plotsize = [n_elements(a(*,0)),n_elements(a(0,*))]
	end else plotsize = [600,400]
end

if isvalid(outfil) then begin
	old_device = !D.name
	set_plot,'Z'
	tv2, plotsize(0), plotsize(1), /init
	loadct, 0
	erase
end

; Assemble list of days for which we need data...
for i=startime,endtime,24*3600.0 do begin
	a = anytim(i,/yymmdd)
	a = strmid(a,0,2)+strmid(a,3,2)+strmid(a,6,2)
	if not isvalid(daylist) then daylist = a else daylist = [daylist,a]
end

; Get list of relevant files...
; This is SOOO tedious!  Give me perl any day!
; At the end, 
f=findfile(getenv('MDI_MON5K_DIR')) ; Read all filenames
r = strpos(f,'.roll')		    ; Find filenames with ".roll"
l = strlen(f) - strlen('.roll')	    ; Lengths of the filenames
w = where((l ne -1) and (r eq l))
if (w(0) eq -1) then message,'No roll files found for the selected dates!'
f = f(w)

days = strmid(f,0,6) ; This should contain the "yymmdd" part of each file.

; Get a list of the relevant files...
for i=0,n_elements(daylist)-1 do begin
	w = where(days eq daylist(i))
	if(w(0) ne -1) then if isvalid(flist) then $
		flist=[flist,f(w)] else flist=f(w)
end

; Read in the relevant files and assemble a master time/roll array
fug = 0
for i=0,n_elements(flist) - 1 do begin
	new = rd_tfile(getenv('MDI_MON5K_DIR')+'/'+flist(i),5)

	if n_elements(new) gt 1 then begin
		if isvalid(teststr) then begin	
			fug = 0
			w = where(new(4,*) eq teststr)
			if w(0) eq -1 then $
				fug = 1 else $
				new = new(*,w)
		end
	
		if not isvalid(teststr) or (fug eq 0) then begin
			newtimes = (new(0,*)+' '+new(1,*))(*)
			newroll = (float(new(tfile_index,*)))(*)
			
			if isvalid(types) then $
				types = [types,(new(4,*))(*)] else $
				types = (new(4,*))(*)

			if isvalid(times) then $
				times = [times,newtimes] else $
				times = newtimes
	
			if isvalid(roll) then $
				roll = [roll,newroll] else $
				roll = newroll
		end
	end
end

;
; WHen we get bad telemetry the date field gets munged.  Filter out lines
; with munged dates.  This is a quick-n-dirty way to do it:  just filter out
; times outside the selected range.
;
t2 = anytim(times)
wok = where((t2 ge startime) and (t2 le endtime))
types = types(wok)
times = times(wok)
roll = roll(wok)

; Plot calculated and orbit-file determinations separately if necessary...
; wcalc gets the locations of calculated-orbit rolls; wfile gets the locations
; of the orbit-file rolls.
wcalc = where(types eq '(CALC)')
wfile = where(types eq '(ORBIT)')

plotmode = 2*(wcalc(0) ne -1) + (wfile(0) ne -1)



; Check for singularities.  If we're right on the singularity, pick the 
; side where the latest data or the bulk of the data are.  If neither 
; heuristic works, use positive numbers.
;
; This is a pretty stoopid heuristic -- we shove the singularity over to 
; 0 and pick (-360 - 0) or (0 - 360) as our range; we don't pick the 
; singularity properly by examining the data in detail.
if(max(roll) gt 170) and (min(roll) lt -170) and $
   ( (where((roll gt -10) and (roll lt 10)))(0) eq -1) then begin
	
	branch='POSITIVE' ; default branch to take

	if(nlm(roll) gt 30) then begin
	  if(total(roll(nlm(roll)-20:nlm(roll)-1) gt 0) gt 15) then begin
	    	branch='POSITIVE' 
	  end else begin
	  	if(total(roll(nlm(roll)-20:nlm(roll)-1) lt 0) gt 15) then begin
			branch='NEGATIVE' 
	    	end else begin
			if(total(roll gt 0) gt (nlm(roll)/2)) then begin
				branch='POSITIVE' 
			end else branch='NEGATIVE'
		end
	  end
	end

;Explicitly check each case -- leaves room for expansion later
 	if(branch eq 'POSITIVE') then begin
		w = where(roll lt 0);guaranteed not to fail so no check needed
		roll(w) = roll(w)+360
	end else if(branch eq 'NEGATIVE') then begin
		w = where(roll gt 0);guaranteed not to fail so no check needed
		roll(w) = roll(w)-360
	end
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; Pick ranges

if(plotmode and 1) ne 0 then begin
	froll = roll(wfile)
	ftimes = times(wfile)

	av1 = total(froll) / nlm(froll)
	stdev1 = total(abs(froll - av1))/nlm(froll)

	if keyword_set(despike) and (nlm(froll) gt 3) then begin
		w = where((abs(froll(1:nlm(froll)-2) - (froll(0:nlm(froll)-3) + froll(2:nlm(froll)-1))/2) gt stdev1 * nsigma) eq 0)
		w = [0,w+1,nlm(froll)-1]
;		print,'despiking: kept ',nlm(w),' of ',nlm(froll),' points.'
		froll = froll(w)
		ftimes = ftimes(w)
	
		av1 = total(froll) / nlm(froll)
		stdev1 = total(abs(froll - av1))/nlm(froll)
	end

	pyrange1 = [min((smooth([froll,froll],10))(nlm(froll)/2:nlm(froll)*3/2.0))<(av1-stdev1*nsigma/2.0) , $
		    max((smooth([froll,froll],10))(nlm(froll)/2:nlm(froll)*3/2.0))>(av1+stdev1*nsigma/2.0) ]
end

if((plotmode and 2) ne 0)then begin
	croll = roll(wcalc)
	ctimes = times(wcalc)

	av2 = total(croll) / nlm(croll)
	stdev2 = total(abs(croll - av2))/nlm(croll)

	if keyword_set(despike) and (nlm(croll) gt 3) then begin
		w = where((abs(croll(1:nlm(croll)-2) - (croll(0:nlm(croll)-3) + croll(0:nlm(croll)-1))/2) gt stdev2 * nsigma) eq 0)
		w = [0,w+1,nlm(croll)-1]
;		print,'despiking: kept ',nlm(w),' of ',nlm(croll),' calculated points.'
		croll = croll(w)
		ctimes = ctimes(w)
	
		av2 = total(croll) / nlm(croll)
		stdev2 = total(abs(croll - av2))/nlm(croll)
	end

	pyrange2 = [min((smooth([croll,croll],10))(nlm(croll)/2:nlm(croll)*3/2.0))<(av2-stdev2*nsigma/2.0) , $
		    max((smooth([croll,croll],10))(nlm(croll)/2:nlm(croll)*3/2.0))>(av2+stdev2*nsigma/2.0)]
end

; Do smoothing if necessary.
; This is quick and dirty, assuming that most of the data are 
; separated by only one minute.  We don't deal with gaps at all!
if not isvalid(av) and isvalid(sm) then av = sm
if keyword_set(av) then begin
	if av gt 5 then window=av else window=5

	if (plotmode and 1) ne 0 then begin
		froll = smooth(froll,window)
		if(av lt n_elements(ftimes)-1) then begin
			ftimes = ftimes(av/2:n_elements(ftimes)-av/2-1)
			froll  = froll(av/2:n_elements(froll)-av/2-1)
		end
	end

	if (plotmode and 2) ne 0 then begin
		croll = smooth(croll,window)

		if(av lt n_elements(ctimes)-1) then begin
			ctimes = ctimes(av/2:n_elements(ctimes)-av/2-1)
			croll = croll(av/2:n_elements(croll)-av/2-1)
		end
	end

	avstr= "("+strtrim(window,2)+" min av)"
end else avstr=""

;
; Initialize plotting
;device,set_resolution=plotsize ; Initialize the output size
ptitle="SOHO Roll Angle from MDI/LOI-V data "+avstr
pytit="Roll Angle ("+columns(tfile_index)+")"

delvarx,pyrange
if isvalid(yrange) then if nlm(yrange) eq 2 then pyrange=yrange

if not isvalid(pyrange) then begin
	if isvalid(froll) then begin
	end 
	if isvalid(croll) then begin
	end
	if isvalid(pyrange1) and not isvalid(pyrange2) then pyrange = pyrange1
	if isvalid(pyrange2) then begin
		if isvalid(pyrange1) then $
			pyrange = [min([pyrange1(0),pyrange2(0)]), $
			  	   max([pyrange1(1),pyrange2(1)])] $ 
		else $
			pyrange = pyrange2
	end 
end

x0 = !D.X_SIZE*0.15
y0 = !D.Y_SIZE*0.2


; Do the plot!
if((plotmode and 1) ne 0) then begin
	utplot,yrange=pyrange,timerange=anytim([startime,endtime],/atime),ftimes,froll,ytit=pytit,title=ptitle,linestyle=0
	
	plots,x0+[0,20],y0+[5,5],linestyle=0,color=255,/device
	xyouts2,/device,x0+25,y0,color=255,charsize=0.75,"From orbital file"

	if(plotmode and 2) ne 0 then begin
		outplot,ctimes,croll,linestyle=1
		plots,x0+120+[0,20],y0+[5,5],linestyle=1,color=255,/device
		xyouts2,x0+120+25,y0,/device,charsize=0.75,"From calc. orbit"
	end
end else if(plotmode and 2) ne 0 then begin
	utplot,yrange=pyrange,timerange=anytim([startime,endtime],/atime),ctimes,croll,ytit=pytit,title=ptitle,linestyle=0
	plots,x0+[0,20],y0+[5,5],linestyle=0,color=255,/device
	xyouts2,x0+25,y0,color=255,/device,charsize=0.75,"From calc. orbit"
end
	
; Note the time it was made...
xyouts2,3,3,'Plot made '+anytim(ut_time(),/atime)+' UT',charsize=0.67,/device,color=255

if isvalid(outfil) then begin
	zbuff2file, outfil, /gif
	set_plot,old_device
end

if keyword_set(debug) then stop
end