;==========================================================================
pro eigpltsetup
; Setup for eigplot.pro
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA4, image, rot_image
ntot=n_elements(eigvals)
if ntot gt 16 then ntot=16
print,'----------------------'
print,'EIGENVECTOR PLOT SETUP'
print,'----------------------'
print,'This routine plots the first n eigenvectors'
restart: read,nplts,prompt='How many eigenvectors do you want to plot? '
if nplts le 0 or nplts gt ntot then begin
	print,'Number out of range: 1 to ',ntot
	print,'Please choose again'
	goto, restart
endif
n=make_array(nplts,/INT,/INDEX)
ans=' '
read,ans,prompt='Also direct spectra to ASCII file (Y/n)? '
!p.charsize=2
!p.color=255
!p.linestyle=0
if ans ne 'n' then begin
	pltname='     '
	read,pltname,prompt='Enter filename [evs.txt]:'
	if pltname eq '' then pltname='evs.txt'
	window,15,title='Eigenvector Spectral Plots',xs=1000,ys=700
	eigplot,n,eigwin=15,zero=1,filename=pltname
endif else begin
	window,15,title='Eigenvector Spectral Plots',xs=1000,ys=700
	eigplot,n,eigwin=15,zero=1
endelse
for imgn=0,nplts-1 do begin
	print,'  Eignevector',imgn
	print,'    Eigenvalue =',eigvals(imgn)
	percentval=eigvals(imgn)*100./total(eigvals)
	print,'               =',percentval,'% of total variance'
endfor ;imgn
return
end

;==========================================================================
pro eigplot,n,vector,ps=ps,noplot=noplot,eigwin=eigwin,mouse=mou,$
	histogram=histo,power=power,oplot=oplot,fill=fill,subtitle=sub,$
	top=top,zero=zero,filename=filename
;
;Written by T.N. Titus, Univ of WYO.
;Last Modification: 23 Sept 1993.
;
;Purpose:  This routine plots the n'th eignvector.
;
;Input:
;	n:  The eigenfunction to be plotted.  NOTE n can be a vector.
;		However, vector will be a matrix: vector(*,[n])
;	noplot: suppress plot of eigenvector.
;	win: window to plot function. Default is zero.
;	mouse: enable mouse routines: location of min & max values.
;		The value of mouse determines window for extra plots.
;	histogram: Plot in histogram mode instead of psym=0.
;	power: if set, the plot will be of power, not weight.
;	oplot: overplot eigenfunction in oplot color.
;	fill: a 2 element vector specifying a region to be filled.
;
;Optional Output:
;	vector: The n'th eigenvector.
;
;
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA4, image, rot_image
;
on_error,2
if keyword_set(power) then begin
	b=eigvecs*eigvecs
	ytt=' Component Power'
	endif else begin
	b=eigvecs
	ytt=' Component Weight'
	endelse
if keyword_set(oplot) then begin
	oplot,wavelength,b(*,n),color=oplot
	return
	endif
if keyword_set(histo) then !p.psym=10
if n_elements(sub) eq 0 then sub=sub_title
if n_elements(n) eq 0 then n=0
no_n=n_elements(n)
if n_elements(top) eq 0 then top=0
set_mp,no_n,/top
if keyword_set(ps) then begin
	ps
	if no_n ge 2 then device,/land
	for ii=0,no_n-1 do plot,wavelength,b(*,n(ii)),$
		title=strith(n(ii))+' Eigenvector',xstyle=3,$
		xtitle=wave_title,ytitle=ytt,subtitle=sub
	endif
if keyword_set(tek_term) then tek else begin
	x
	if n_elements(eigwin) ne 0 then wset,eigwin else wset,0
	endelse
if not keyword_set(noplot) then $
	for ii=0,no_n-1 do begin
	plot,wavelength,b(*,n(ii)),title=strith(n(ii))+' Eigenvector',xstyle=3,$
		xtitle=wave_title,ytitle=ytt,subtitle=sub,color=255,$
		xrange=[min(wavelength),max(wavelength)],$
		yrange=[min(b[*,n[ii]]),max(b[*,n[ii]])]
	if keyword_set(zero) then oplot,wavelength,wavelength*0,color=zero
	if n_elements(fill) ne 0 then for jj=0,(n_elements(fill)-1)/2 do begin
		lft=(!d.name eq 'TEK')
		wpos=where((wavelength ge fill(0+jj*2)) and (wavelength le fill(1+jj*2)))
		polyfill,[fill(0+jj*2),wavelength(wpos),fill(1+jj*2),fill(0+jj*2)],$
			[0,b(wpos,n(ii)),0,0],$
			line_fill=lft,color=6-jj,orientation=45,spac=lft+0.01
		endfor
	endfor
if keyword_set(mou) then $
	mouse,pict=image,absc=wavelength,pixel=indgen(n_elements(wavelength)),$
		ordin=b(*,n),win=mou
if no_n ge 2 then !p.multi=0
if keyword_set(histo) then !p.psym=0
if n_elements(filename) ne 0 then begin
	namearr=strarr(no_n+1)
	namearr(0)=wave_title
	namearr(1:*)='Eigenvector'+strno(n(*),remove=1)
	fwave=transpose(wavelength)
	fb=transpose(b(*,n))
	farray=[fwave,fb]
	openw,10,filename,width=160
	printf,10,namearr
	printf,10,farray
	close,10
endif
set_mp
return
end

;==========================================================================
pro lampltsetup
; Setup for lamplot.pro
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA4, image, rot_image
ntot=n_elements(eigvals)
print,'----------------------'
print,'EIGENVALUE PLOT SETUP'
print,'----------------------'
print,'This routine plots the first n eigenvaluess'
restart: 
read,nevs,prompt='How many eigenvalues do you want in the plot? '
if nevs gt ntot then goto, restart
ans=' '
read,ans,prompt='Include 0th eigenvalue (Y/n)? '
if ans eq 'n' then nz=1 else nz=0
if nz and nevs+1 gt ntot then goto, restart
ans=' '
read,ans,prompt='Liner plot (Y=linear/n=log)? '
if ans eq 'n' then lgplot=1 else lgplot=0
ans=' '
read,ans,prompt='Also direct values to ASCII file (Y/n)? '
lct_dmw
!p.charsize=2
!p.color=255
!p.linestyle=0
if ans ne 'n' then begin
	pltname='     '
	read,pltname,prompt='Enter filename [evals.txt]:'
	if pltname eq '' then pltname='evals.txt'
	window,11,title='Eigenvalue Plot',xs=400,ys=400
	lamplot,dim=nevs,symbol=4,log=lgplot,nozero=nz,filename=pltname
endif else begin
	window,11,title='Eigenvalue Plot',xs=400,ys=400
	lamplot,dim=nevs,symbol=4,log=lgplot,nozero=nz
endelse
for imgn=nz,nevs+nz-1 do begin
	print,'  Eignevector',imgn
	print,'    Eigenvalue =',eigvals(imgn)
	percentval=eigvals(imgn)*100./total(eigvals)
	print,'               =',percentval,'% of total variance'
endfor ;imgn
return
end

;==========================================================================
pro lamplot,evals,ps=ps,noplot=noplot,arm=arm,dim=dim,symbol=sym,$
        broken_stick=broke,pie=pie,label=label,log=log,$
		nozero=nozero,filename=filename
;
;Written by: T.N. Titus, Univ of Wyo, 29 Jun 1992.
;Last Modification: 2 Nov 95.
;
;This procedure plots the eigenvalues calculated by PCAINIT.
;
;Input: ps - plots to a post script file called idl.ps.
;       noplot - do not plot eigenvalues on screen.
;       arm - overplot the aveage root * arm.
;       symbol: use symbol with lines connecting.
;       broken_stick: another overplot. This is the Horn Method (Jolliffe,86)
;   label: print out the eigenvalues on the plot
;
;Keywords: -pie - dispaly as a pie graph
;
;Output: Evals contains the eigenvalues (Optional)
;
common pca0,subt,tek
common pca2,ev,vecs
;
if keyword_set(pie) then begin
        pie,ev,/lab,log=log
        return
        endif
if keyword_set(nozero) then evals=ev[1:*] else evals=ev
if n_elements(dim) ne 0 then evals=evals(0:dim-1)
idx=indgen(n_elements(evals))
if keyword_set(nozero) then idx=idx+1
if n_elements(sym) eq 0 then sym=0
fidx=findgen(n_elements(ev))
if keyword_set(broke) then horn=(1./(fidx+1.))#trimat(n_elements(fidx))
;
if keyword_set(ps) then begin
        ps
;        plot,idx,evals,title='Eigenvalues',xtitle='Principle Component',$
;			ytitle='Eigenvalue',xstyle=2,ystyle=2,subt=subt,$
;			psym=sym,xrange=[0,max(idx)+2]
		plot,idx,evals,xrange=[0,max(idx)+1],ylog=log, color=255, $
			title='Eigenvalues',xtitle='Principle Component',$
			ytitle='Eigenvalue',xstyle=2,ystyle=2,subt=subt,psym=sym, $
			yrange=[min(evals),max(evals)]
		oplot,idx,evals
        if keyword_set(arm) then oplot,idx,idx*0+arm*mean(evals)
        if keyword_set(broke) then oplot,idx,horn
        if keyword_set(tek) then tek else x
 endif
;
if (not keyword_set(noplot)) then begin
	plot,idx,evals,xrange=[0,max(idx)+1],ylog=log, color=255, $
		title='Eigenvalues',xtitle='Principle Component',$
		ytitle='Eigenvalue',xstyle=2,ystyle=2,subt=subt,psym=sym, $
		yrange=[min(evals),max(evals)]
	oplot,idx,evals
	if keyword_set(arm) then oplot,idx,idx*0+arm*mean(evals)
	if keyword_set(broke) then oplot,idx,horn
	if keyword_set(label) then xyouts,idx+0.85,evals,evals
endif
if n_elements(filename) ne 0 then begin
	namearr=['      Num','      Eval','        %Var']
	farray=strarr(3,n_elements(idx)+1)
	farray[*,0]=namearr
	farray[0,1:n_elements(idx)]=string(transpose(idx))
	farray[1,1:n_elements(idx)]=string(transpose(evals[idx[0]:$
		idx[n_elements(idx)-1]]))
	farray[2,1:n_elements(idx)]=string(transpose(evals[idx[0]:$
		idx[n_elements(idx)-1]]/total(evals)*100))
	openw,10,filename
	printf,10,farray
	close,10
endif
;
return
end

;==========================================================================
pro set_mp,n,top=top
;
;Written by T.N. Titus, 29 Sept 93.
;Last Modified: 25 Oct 93.
;
;Purpose: To set up the screen for multiple plots.
;
;Input:
;	n: the number of plots. [default is n=0 which resets screen]
;Keyword:
;	top: plot top down. [default is side ways.]
;
on_error,2
if n_elements(n) eq 0 then n=0
if n gt 16 then message,'N to large: Max value is 16.'
v0=intarr(1,17)
v1=[0,1,1,1,2,2,2,3,3,3,3,3,3,4,4,4,4]
v2=[0,1,2,3,2,3,3,3,3,3,4,4,4,4,4,4,4]
v1=transpose(v1)
v2=transpose(v2)
if keyword_set(top) then arr=[v0,v1,v2] else arr=[v0,v2,v1]
!p.multi=arr(*,n)
return
end

;==========================================================================
pro endspcsetup
; Setup for endspec.pro
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON RES2, v_names, StoN
COMMON RES4, verts, verts_set
print,'--------------------'
print,'ENDMEMBER PLOT SETUP'
print,'--------------------'
print,'This routine plots the spectra of your chosen endmembers'
verts=EM_pix
verts_set=EM_pix_set
;v_names=strarr(N_EndM)
n=make_array(N_EndM,/INT,/INDEX)
ans=' '
read,ans,prompt='Also direct spectra to ASCII file (Y/n)? '
!p.charsize=2
!p.color=255
!p.linestyle=0
if ans ne 'n' then begin
	pltname='     '
	read,pltname,prompt='Enter filename [ems.txt]:'
	if pltname eq '' then pltname='ems.txt'
	window,16,title='Endmember Spectral Plots',xs=1000,ys=700
	endspec,n,reswin=16,filename=pltname
endif else begin
	window,16,title='Endmember Spectral Plots',xs=1000,ys=700
	endspec,n,reswin=16
endelse
return
end

;==========================================================================
pro endspec, n, filename=filename, ps=ps, reswin=reswin
;+
; NAME:
;	endspec
; PURPOSE:
;	To plot (and write to an ASCII file) calculated endmemeber spectra
;	from T. Titus' PCA code.  Note, this program will only work after
;	variables have been set using pcainit.pro and resinit.pro
; INPUTS:
;	n = Endmember number of endmember to plot.  This is the nth
;	    element in the vector verts
;	filename = file to write extracted spectra
;	ps = 1 (or /ps) to make a postscript plot.  Note that this file must
;	     be closed before printing/ftp'ing.  Do this by typing:
;	     ps & device, /close & x
; OUTPUTS:
;	None.  ASCII file is written if filename is given
; PROCEDURES USED:
;	common blocks from pcainit and resinit, strith
; MODIFICATION HISTORY:
;	David R. Klassen, April 1996 - written
;-
	On_error, 1
if n_params() lt 1 then message, /NONAME, $
	'Syntax:  endspec, n[, filename=filename'
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON RES2, v_names, StoN
COMMON RES4, verts, verts_set
;
no_n=n_elements(n)
set_mp,no_n,/top
imgerr=fltarr(no_n,n_elements(wavelength))
if keyword_set(ps) then begin
	ps
	if no_n ge 2 then device, /land
	for i=0,no_n-1 do begin
			print,verts_set(n(i),*)
		errset=verts_set(n(i),where(verts_set(n(i),*) ge 0))
			print,errset
		if n_elements(errset) gt 1 then $
			for l=0,n_elements(wavelength)-1 do imgerr(i,l)=stdev(PCArr(errset,l))
		plot,wavelength,PCArr(verts(n(i)),*),title=$
			strith(n(i))+'Endmember: '+v_names(n(i)),xstyle=3,$
			xtitle=wave_title,ytitle='Reflectance',subtitle=sub_title
		errplot,wavelength,PCArr(verts(n(i)),*)+imgerr(i,*),PCArr(verts(n(i)),*)-imgerr(i,*)
	endfor ;i
endif
x
if n_elements(reswin) ne 0 then wset,reswin else wset,0
for i=0,no_n-1 do begin
	errset=verts_set(n(i),where(verts_set(n(i),*) ge 0))
	if n_elements(errset) gt 1 then $
		for l=0,n_elements(wavelength)-1 do imgerr(i,l)=stdev(PCArr(errset,l))
	plot,wavelength,PCArr(verts(n(i)),*),title=strith(n(i))+' Endmemeber:  '+$
		v_names(n(i)),xstyle=3,xtitle=wave_title,ytitle='Reflectance',$
		subtitle=sub_title
	errplot,wavelength,PCArr(verts(n(i)),*)+imgerr(i,*),PCArr(verts(n(i)),*)-imgerr(i,*)
endfor ;i
if n_elements(filename) ne 0 then begin
	namearr=make_array(2*no_n+1,/string,value='test')
	namearr(0)='Wavelength'
	namearr(1:no_n)=v_names(n(*))
	errnames='sigma_'+v_names(n(*))
	namearr(no_n+1:no_n+no_n)=errnames
	fwave=transpose(wavelength)
	fspec=PCArr(verts(n),*)
	ferror=imgerr
	farray=[fwave,fspec,ferror]
	openw, 10, filename,width=160
	printf,10,namearr
	printf,10,farray
	close,10
endif
set_mp
return
end

;==========================================================================
pro pamtvsetup
; Setup for pamtv.pro
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON PCA0, sub_title, tek_term
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA5, objects
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
ntot=n_elements(eigvals)
print,'-----------------------'
print,'EIGENVECTOR IMAGE SETUP'
print,'-----------------------'
print,'This routine will display an eigen-image'
ans=''
read,ans,$
	prompt='   FAST: Display 0--4 and output to default files [Y/n]? '
if ans ne 'n' then begin
	for imgn=0,4 do begin
		deffname='ev'+strno(imgn,/remove)
		title=strith(imgn)+' Eigenvector Image'
		pamtv,imgn,pamwin=17,zoom=scl,filename=deffname,wtitle=title
		print,'   Eigenvalue =',eigvals(imgn)
		print,'              =',eigvals(imgn)*100./total(eigvals),$
			'% of total variance'
	endfor
endif else begin
	restart: read,imgn,prompt='Which eigen-image do you want to display? '
	if imgn lt 0 or imgn ge ntot then begin
		print,'Number out of range: 0 to ',ntot-1
		print,'Please choose again'
		goto, restart
	endif
	imgn=fix(imgn)
	deffname='ev'+strno(imgn,/remove)
	ans=' '
	read,ans,prompt='Also also write to GIF/FITS files (Y/n)? '
	;window,17,xs=InX*scl,ys=InY*scl,title=strith(imgn)+' Eigenvector Image'
	title=strith(imgn)+' Eigenvector Image'
	if ans ne 'n' then begin
		imgname='     '
		read,imgname,prompt='Enter filename w/o extensions ['+deffname+']:'
		if imgname eq '' then imgname=deffname
		pamtv,imgn,pamwin=17,zoom=scl,filename=imgname,wtitle=title
	endif else $
		pamtv,imgn,pamwin=17,zoom=scl,wtitle=title
	print,'   Eigenvalue =',eigvals(imgn)
	print,'              =',eigvals(imgn)*100./total(eigvals),$
		'% of total variance'
endelse
return
end

;==========================================================================
pro pamtv,n,ps=ps,pamwin=pamwin,surface=sur,contour=con,zoom=zoom,neg=neg,$
	eight=eight,log=log,filename=filename,wtitle=wtitle
;
;Lasted updated 8 November 1995
;Updated 22 APR 96 by David R. Klassen - added filename option to output
;               a FITS image with informational header
;
;This procedure will plot n'th eigenfunction values vs param.
;
;Input: n - the component to be viewed. (the eigenvector)
;
;Options:
;	ps: Output a postscript plot to idl.ps
;	window: display image in window no.
;	surface: do a surface plot.
;	contour: do a contour plot [default for tek]
;	zoom: mainly for use with tvscl to enlarge the image.
;	neg: use the negative of the image.
;	eight: color code image as in pca3plot - black are masked out.
;	log: use a logrithmic display
;	filename:  name of file to save eigenfunction image in FITS format
;	windowtitle: title of the display window for gifbar
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON PCA0, sub_title, tek_term
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA5, objects
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
image_size=[XSize,YSize]
type=(keyword_set(sur)*2+keyword_set(con))
if n_elements(n) eq 0 then n=0
n_obj=n_elements(new_vecs(*,n))
picture=reform(new_vecs(*,n),image_size(0),image_size(1))
if keyword_set(neg) then picture=-picture
if keyword_set(eight) then begin
	colours=bytscl(picture,top=6,max=max(picture(objects)),$
		min=min(picture(objects)))+1
	blank=picture*0
	blank(objects)=1
	picture=blank*colours
endif
if n_elements(wtitle) eq 0 then wtitle='Eigen-map #'+strno(n)
if n_elements(zoom) ne 0 then picture=rebin(picture,image_size(0)*zoom,$
	image_size(1)*zoom,/sample)
if keyword_set(log) then picture=alog(picture-min(picture)+1)
if keyword_set(ps) then begin
	ps
	case type of
		0: imgprt,picture,wtitle,/border
		1: contour,picture
		2: surface,picture
		3: message,'Incompatable keywords'
		endcase
	endif
if keyword_set(tek_term) then tek $
else begin
	x
	if n_elements(pamwin) eq 0 then pamwin=0
endelse
case type of
;	0: tvscl,picture
	0: gifbar,picture,dispwin=pamwin,wintitle=wtitle
	1: contour,picture
	2: surface,picture
	3: message,'Incompatable keywords'
	endcase
pcmap=reform(new_vecs(*,n),image_size(0),image_size(1))
if n_elements(filename) ne 0 then begin
	BufImage=rebin(pcmap,Xsize*xygrid,YSize*xygrid)
	OutImage=fltarr(XSize,YSize)
	OutImage(xmargin:XSize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
	mkhdr,hdr,OutImage
	comment='Eigenvector number'
	sxaddpar,hdr,'eig_num',n,comment
	comment='Eigenvalue'
	sxaddpar,hdr,'eig_val',eigvals(n),comment
	comment='% of Varience'
	sxaddpar,hdr,'var_%',(eigvals(n)/total(eigvals))*100.,comment
	writefits,filename+'.fit',OutImage,hdr
;	gifbar,OutImage,fname=filename+'.gif',ct=0,dispwin=20
;	gifbar,OutImage,fname=filename+'.c.gif',dispwin=20
	gifbar,OutImage,fname=filename,ct=0,dispwin=20
	gifbar,OutImage,fname=filename,dispwin=20
	wdelete,20
endif
return
end

;==========================================================================
pro dispresidsetup
;
; The procedure get some of the variables setup for restv to display a
; residual spectral image
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
             X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
print,'-----------------------------'
print,'RESIDUAL SPECTEAL IMAGE SETUP'
print,'-----------------------------'
print,'This routine will display a residual spectral image'
dispunit=0
xmenu,['Absolute Image Units','Relative Image Units'], base=base, $
	buttons=B,title='Image Unit Choices'
WIDGET_CONTROL,base,/realize
event=WIDGET_EVENT(base)
WIDGET_CONTROL, get_uvalue=dispunit, event.id
WIDGET_CONTROL, event.top,/destroy
case dispunit of
	0: begin
		print,'   Absolute image units chosen'
		dfn2='abs'
	end ;0
	1: begin
		print,'   Relative image units chosen'
		dfn2='rel'
	end ;1
endcase
read,imgn,prompt='Which residual spectral image to display (-1 for RMS)? '
imgn=fix(imgn)
ans=' '
if imgn ne -1 then begin
	title=strith(imgn)+' Residual Spectral Image'
	dfn1='r'+strno(imgn,/remove)
endif else begin
	title=' RMS Residual Spectral Image'
	dfn1='rms'
endelse
read,ans,prompt='Also also write to GIF/FITS files (Y/n)? '
if ans ne 'n' then begin
	deffname=dfn1+dfn2+'resid'
	imgname='     '
	read,imgname,prompt='Enter filename w/o extensions ['+deffname+']:'
	if imgname eq '' then imgname=deffname
	restv,imgn,reswin=17,zoom=scl,filename=imgname,/verbose,wtitle=title,$
		disprelunit=dispunit
endif else $
	restv,imgn,reswin=17,zoom=scl,/verbose,wtitle=title,disprelunit=dispunit
return
end

;==========================================================================
pro restv,n,resmap,ps=ps,reswin=reswin,surface=sur,contour=con,$
	zoom=zoom,neg=neg,filename=filename,abs=abs,verbose=verbose,$
	wtitle=wtitle,disprelunit=disprelunit
;
;Written by T.N. Titus, 19 December 1995.
;
;This procedure will image the n'th filter residue.
;	The default is residual magnitude.
;
;Input: n - the filter to be viewed.
;
;Options:
;	ps: Output a postscript plot to idl.ps
;	window: display image in window no.
;	surface: do a surface plot.
;	contour: do a contour plot [default for tek]
;	zoom: mainly for use with tvscl to enlarge the image.
;	neg: use the negative of the image.
;	windowtitle: title for display window for gifbar
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA5, objects
COMMON PCA6, image_size
COMMON RES1, fr, recon
COMMON RES2, vname, StoN
COMMON RES3, orig_mag
COMMON RESIM, resid
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
fnddim,resid,c,r
image_size=[XSize,YSize]
type=(keyword_set(sur)*2+keyword_set(con))
if n_elements(n) eq 0 then n=-1
if n ne -1 then picture=resid(*,n) else picture=sqrt(rebin(resid*resid,c,1))
if keyword_set(disprelunit) then begin
	if n ne -1 then begin
		print,'Displaying Residual/Original'
		print,'    Note: original image pixels that had a value of 0.0'
		print,'    were set to -30000.0 before division to avoid errors.'
		relimg=PCArr(*,n)
		ztest=where(relimg eq 0.0)
		if ztest(0) ne -1 then relimg(ztest)=-30000.0
		picture=picture/relimg
	endif else begin
		print,'Displaying  Wavelength_RMS(Residual)/Wavelength_Mean(Original)'
		print,'    Note: original image pixels that had a value of 0.0'
		print,'    were set to -30000.0 before division to avoid errors.'
		relimg=rebin(PCArr,c,1)
		ztest=where(relimg eq 0.0)
		if ztest(0) ne -1 then relimg(ztest)=-30000.0
		picture=picture/relimg
	endelse
endif else print,'Displaying Residuals in Image Intesity Units'
picture=reform(picture,image_size(0),image_size(1))
resmap=picture
if keyword_set(neg) then picture=-picture
if keyword_set(abs) then picture=abs(picture)
if n_elements(wtitle) eq 0 then begin
	strtitle=['',': filter#'+strno(indgen(15))]
	wtitle='Residual Map'+strtitle(n+1)
endif
if n_elements(zoom) ne 0 then picture=rebin(picture,image_size(0)*zoom,$
	image_size(1)*zoom,/sample)
if keyword_set(ps) then begin
	case type of
		0: imgprt,picture,wtitle,/border
		1: contour,picture
		2: surface,picture
		3: message,'Incompatable keywords'
		endcase
	endif
if keyword_set(tek_term) then tek $
else begin
	x
	if n_elements(reswin) eq 0 then reswin=0
endelse
case type of
;	0: tvscl,picture
	0: gifbar,picture,dispwin=reswin,wintitle=wtitle
	1: contour,picture
	2: surface,picture
	3: message,'Incompatable keywords'
	endcase
resstdev=stdev(resmap,resmean)
if n_elements(filename) ne 0 then begin
	BufImage=rebin(resmap,Xsize*xygrid,YSize*xygrid)
	OutImage=fltarr(XSize,YSize)
	OutImage(xmargin:XSize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
	mkhdr,hdr,OutImage
	if n ne -1 then begin
		comment=wave_title
		sxaddpar,hdr,'res_wave',wavelength(n),comment
	endif else begin
		comment='RMS Residual over all wavelengths'
		sxaddpar,hdr,'res_wave',0,comment
	endelse
	comment='Image Mean'
	sxaddpar,hdr,'mean',resmean,comment
	comment='Image Standard Deviation'
	sxaddpar,hdr,'std_dev',resstdev,comment
	writefits,filename+'.fit',OutImage,hdr
;	gifbar,OutImage,fname=filename+'.gif',ct=0,dispwin=20
;	gifbar,OutImage,fname=filename+'.c.gif',dispwin=20
	gifbar,OutImage,fname=filename,ct=0,dispwin=20
	gifbar,OutImage,fname=filename,dispwin=20
	wdelete,20
endif
if keyword_set(verbose) then begin
	if n ne -1 then print,'Wavelength =',wavelength(n) $
		else print,'RMS over all wavelengths'
	print,'Mean = ',resmean,'     St Dev = ',resstdev
endif
return
end

;==========================================================================
pro mrestv,reswin=reswin,zoom=zoom,neg=neg,abs=abs,verbose=verbose,wtitle=wtitle
;
;Written by T.N. Titus, 19 December 1995.  Modified to show all residual
;images in a 5x5 display by David Klassen, Nov. 1996.
;
;This procedure will display all residual images as a function of
;wavelength 25 at a time in a 5x5 window.
;
;Input: none
;
;Options
;	reswin: display image in window no.
;	zoom: mainly for use with tvscl to enlarge the image.
;	neg: use the negative of the image.
;	wtitle: title for display window
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA5, objects
COMMON PCA6, image_size
COMMON RES1, fr, recon
COMMON RES2, vname, StoN
COMMON RES3, orig_mag
COMMON RESIM, resid
;
fnddim,resid,c,r
image_size=[XSize,YSize]
if n_elements(wtitle) eq 0 then wtitle='Residual Maps'
if n_elements(reswin) eq 0 then reswin=0
if n_elements(zoom) ne 0 then begin
	winx=XSize*5*zoom
	winy=Ysize*5*zoom
endif else begin
	winx=XSize*6
	winy=Ysize*6
endelse
x
window,reswin,xs=winx,ys=winy,title=wtitle
sets=fix(NImages/36)+1
step=0
while step lt sets do begin
	erase
	n=0
	while step*36+n lt NImages and n lt 36 do begin
		picture=reform(resid(*,step*36+n),image_size(0),image_size(1))
		if keyword_set(neg) then picture=-picture
		if keyword_set(abs) then picture=abs(picture)
		if n_elements(zoom) ne 0 then picture=rebin(picture,image_size(0)*zoom,$
			image_size(1)*zoom,/sample)
		tvscl,picture,n
		resstdev=stdev(picture,resmean)
		if keyword_set(verbose) then begin
			print,'  Image ',n,'  Wavelength =',wavelength(n)
			print,'    Mean = ',resmean,'     St Dev = ',resstdev
		endif
		n=n+1
	endwhile; n
	continue=' '
	read,continue,prompt='Hit RETURN to continue'
	step=step+1
endwhile; step
return
end

;==========================================================================
PRO an_subset
;
;Procedure for analyzing a subset of non-contiguous pixels chosen from
;the fractional abundance maps
;
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
              YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON IMAG, Data1, NImages, Im_Num, Hi, Lo, Image, X_im, Y_im
COMMON PCAP, Elem1, Elem2, Elem3, Elem4, Elem5, Elem6, Element1, Element2,$
             X_pca, Y_pca, width, c_pix1, win, xrange, yrange
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON COLOR, c_pix, Num_pix
COMMON ENDM, N_EndM, End_M, EM_pix, EM_pix_set, EM_npca
COMMON FRACAB, Frac, Frac2
COMMON RES, res_pix, res_pix2
COMMON WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON SPEC2, old_vecs, n_oem_pix, N_iter, old_em_pix, em_pca_pix,N_EM
COMMON RES2, vname, StoN
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
print,'=================================================='
print,'   DOING PCA/LMM ON NON-CONTIGUOUS PIXEL SUBSET'
print,'=================================================='
res_pix2=res_pix
Continue=0
while (Continue ne 1) do begin
	menu=0
	xmenu, ['0 - '+vname(0),'1 - '+vname(1),'2 - '+vname(2), $
		'3 - '+vname(3),'4 - '+vname(4)], base=base, button=B, $
		title='Frac Map?'
	widget_control, base, /realize
	event = widget_event(base)
	widget_control, get_uvalue=menu, event.id
	widget_control, event.top, /destroy
	case menu of
	0: begin
		f=0
		continue=1
	end ;0
	1: begin
		f=1
		continue=1
	end ;1
	2: begin
		if N_Endm gt 2 then begin
			f=2
			continue=1
		endif else begin
			print,'Endmember out of range - rechoose'
			continue=0
		endelse
	end ;2
	3: begin
		if N_Endm gt 3 then begin
			f=3
			continue=1
		endif else begin
			print,'Endmember out of range - rechoose'
			continue=0
		endelse
	end ;3
	4: begin
		if N_Endm gt 4 then begin
			f=4
			continue=1
		endif else begin
			print,'Endmember out of range - rechoose'
			continue=0
		endelse
	end ;4
	endcase
endwhile
print,' '
print,'  Using Fractional Abundance Map',f
read,cutoff,prompt='Enter CUTOFF value: '
continue=0
while (continue ne 1) do begin
	menu=0
	xmenu,['<CUTOFF','       <=CUTOFF       ','=CUTOFF','>=CUTOFF','>CUTOFF'], $
		base=base,title='Use pixels that are...'
	widget_control, base, /realize
	event = widget_event(base)
	widget_control, get_uvalue=menu, event.id
	widget_control, event.top, /destroy
	case menu of
	0: begin
		res_pix=where(Frac(*,f) lt cutoff)
		continue=1
	end; 0
	1: begin
		res_pix=where(Frac(*,f) le cutoff)
		continue=1
	end; 1
	2: begin
		res_pix=where(Frac(*,f) eq cutoff)
		continue=1
	end; 2
	3: begin
		res_pix=where(Frac(*,f) ge cutoff)
		continue=1
	end; 3
	4: begin
		res_pix=where(Frac(*,f) gt cutoff)
		continue=1
	end; 4
	endcase
endwhile
print,'  Using Cutoff = ',cutoff
print,' '
n_r_pix = n_elements(res_pix)
print, 'Number of Residual pixels',n_r_pix
printf,5,'Number of Residual pixels',n_r_pix
;
;Saving old PCA coordinates.
;old_vecs = new_vecs
;
;Finding subset of good_pix and res_pix for PCAINIT.
good_pix2 = intarr(XSize*YSize)
good_pix2 = long(good_pix2)
good_pix3 = intarr(XSize*YSize)
good_pix3 = long(good_pix3)
good_pix4 = intarr(XSize*YSize)
good_pix4 = long(good_pix4)
good_pix2(res_pix) = 1
good_pix3(good_pix) = 1
good_pix4 = good_pix3 + good_pix2
good_pix = where(good_pix4 eq 2)
print,'# good elements = ', n_elements(good_pix)
;
;Running high residual pixels through PCA again.
pcainit,PCArr,object=good_pix ,tek=0 ;in pcalib3.pro
;
;Finding pca coordinate for the sky pixel.
;Sky_pca2 = newcoo(Sky)
;
;Finding pca coordinates for old end-members.
EM_npca = fltarr(n_oem_pix,NImages)
for i = 0, n_oem_pix-1 do begin
  EM_npca(i,*) = newcoo(PCArr(old_em_pix(i),*)) ;in pcalib3.pro
endfor ;i
;Opening up windows in order to display new PCAPLOTs.
window,2,xs=400,ys=400,title='New PCA Plot 3',retain=2
window,1,xs=400,ys=400,title='New PCA Plot 2',retain=2
window,0,xs=400,ys=400,title='New PCA Plot 1',retain=2
;window,2,xs=600,ys=400,title='New PCA Plot 3',retain=2
;window,1,xs=600,ys=400,title='New PCA Plot 2',retain=2
;window,0,xs=600,ys=400,title='New PCA Plot 1',retain=2
;
;Displaying PCAPLOTs.
Elem1 = 0
Elem2 = 1
Elem3 = 0
Elem4 = 2
Elem5 = 2
Elem6 = 1
wset,2
lct_dmw
pcaplot,Elem5,Elem6,0,ps=1,iter=N_iter+1,outfile=2 ;in pcalib3.pro
;oplot,[Sky_pca2(Elem5)],[Sky_pca2(Elem6)],psym=1,color=1
for jj = 0, n_oem_pix-1 do begin
  oplot,[EM_npca(jj,Elem5)],[EM_npca(jj,Elem6)],psym=1,color=4
endfor
w2x = !x
w2y = !y
wset,1
pcaplot,Elem3,Elem4,1,ps=1,iter=N_iter+1,outfile=1 ;in pcalib3.pro
;oplot,[Sky_pca2(Elem3)],[Sky_pca2(Elem4)],psym=1,color=1
for jj = 0, n_oem_pix-1 do begin
  oplot,[EM_npca(jj,Elem3)],[EM_npca(jj,Elem4)],psym=1,color=4
endfor
w1x = !x
w1y = !y
wset,0
pcaplot,Elem1,Elem2,2,ps=1,iter=N_iter+1,outfile=0 ;in pcalib3.pro
;oplot,[Sky_pca2(Elem1)],[Sky_pca2(Elem2)],psym=1,color=1
for jj = 0, n_oem_pix-1 do begin
	oplot,[EM_npca(jj,Elem1)],[EM_npca(jj,Elem2)],psym=1,color=4
endfor ;jj
w0x = !x
w0y = !y
clean = intarr(4)
clean(*) = 0
;
;Select new set of end_members for residual pixels
Choose = 0
while (Choose ne 1) do begin;
	xmenu, ['DISPLAY PCA pixels','New PCA Plot','REFRESH Image and PCA plots',$
		'SELECT Endmembers','MEASURE Fractional Abundances', $
		'ITERATE Residuals Again','Change COLOR PALETTE','STOP'], $
		BASE=base, BUTTONS=B,title='Residual Analysis Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = menu, event.id
	WIDGET_CONTROL,event.top,/DESTROY
;
	case menu of
	;Select pca pixels to be displayed
	0: begin
		refr_pca,/color ;in pcalib1.pro
		Ans = 'y'
		color=1
		b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
		c = reform(b,Long(XSize)*Long(YSize))
		c_pix2 = bytscl(c) > 8b
		spec = fltarr(10,WB)
		specerr = fltarr(10,WB)
		window,4,xs=XSize*scl,ys=YSize*scl,title='PCA pixels',retain=2,$
			yp=0,xp=300
		while (Ans ne 'n') do begin
			s_pcaplot ;in pcalib1.pro
			get_box,color ;in pcalib2.pro
			;Display image with pixels highlighted in color from color_im
			pca_pix ;in pcalib1.pro
			if (Num_pix le 0) then print, 'No points inside boxed area.'
			;Display selected points, in color, in window 0.
			if (Num_pix gt 0) then begin
				c_pix2(c_pix1) = color
				c_pix3 = reform(c_pix2,XSize,YSize)
				c_pix4 = rebin(c_pix3,XSize*scl,YSize*scl,/sa)
				specextract,spec,specerr,color ;in pcatool.pro
				wset,4
				tv, c_pix4
			endif
			c_pix=c_pix1
			d_c_pca,color
			color = color+1
			xmenu,['    No    ','    Yes    '],base=base, $
				buttons=b,title='Again?'
			WIDGET_CONTROL,base,/REALIZE
			event = WIDGET_EVENT(base)
			WIDGET_CONTROL, GET_UVALUE = dispmore, event.id
			WIDGET_CONTROL, event.top, /DESTROY
			case dispmore of
			0: Ans='n'
			1: Ans='y'
			endcase
		endwhile
		specwin,spec,specerr,color ;in pcatool.pro
		clean(*)=1
        print,'*****   Display Routine Complete    *****'
        print,'***** Ready for Endmember Selection *****'
	end ;0
	; New PCA Plot
	1: begin
		N_iter=N_iter+1
		d_pcaplot ;in pcalib1.pro
		N_iter=N_iter-1
        print,'*****    New PCA Plot Displayed     *****'
        print,'***** Ready for Endmember Selection *****'
	end ;1
	; Refreshing Image and PCA plots
	2: begin
		print,'Refreshing Image and PCA plots...'
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
        print,'***** Image and PCA Plot Refreshed  *****'
        print,'***** Ready for Endmember Selection *****'
	end ;2
	;Choose endmembers
	3: begin
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		s_endm ;in pcalib2.pro
        print,'*****     Endmembers Selection Complete      *****'
        print,'***** Ready to Calcualte Frational Abuncaces *****'
 	end ;3
	;Measure endmember abundances
	4: begin
		N_iter = N_iter+1
		frac_ab ;in pcalib2.pro
		d_frac ;in pcalib2.pro
		d_resid ;in pcalib2.pro
        print,'*****   Frational Abundances and
        print,'*****   Frational Abundances and Residuals Calculated    *****'
        print,'*****  For Endmember and Eigenvector Display Functions   *****'
        print,'***** Choose ITERATE RESIDUALS AGAIN and make selections *****'
        print,'*****                 From Analysis Menu                 *****'
    end ;4
    ;Repeat iteration
    5: begin
        choose=1
        print,'=============================================='
        print,'CHOOSE ANALYZE RESIDUALS FROM MENU'
        print,'OR CHOOSE ONE OF THE DISPLAY FUNCTIONS'
        print,'=============================================='
        print
    end ;5
    ;Change color palette
    6: begin
        read,colorvalue,prompt='Enter color palette number (0-41): '
        if colorvalue le 40 then loadct,colorvalue $
        else if colorvalue eq 41 then loadct,0
        print,'***** Color Palette Changed *****'
        print,'***** Continue with analysis *****'
    end ;6
    ;Save and Stop
    7: begin
        loadct,0
        choose =1
    end ;7
    endcase
endwhile
;
for i = 1,N_EndM-1 do old_em_pix(i-1+n_oem_pix) = EM_pix(i)
n_oem_pix = n_oem_pix + N_EndM - 1
N_EM(N_iter) = N_EndM-1
;
return
end
;==========================================================================
;The following sections deaing with plotting a single pixel's original,
;model, and residual spectra isn't quite working.  The option was removed
;from the Analysis Tools Menu.

;==========================================================================
pro spcpltsetup
; Setup for specplot.pro
common pca0,subt,tek
common pca1,wave,eigvect,wavetitle
COMMON PCA2,eigvals,new_vecs,PCArr
common pca6,dim_image
COMMON FRACAB, Frac,Frac2
COMMON ENDM, N_EndM,End_M,EM_pix, EM_pix_set, EM_npca
COMMON RES, res_pix,res_pix2
print,'-------------------------'
print,'PIXEL SPECTRAL PLOT SETUP'
print,'-------------------------'
print,'This routine plots spectra for a single pixel.  As specified it'
print,'will plot the original, model, and/or residual spectra.'
read,pixel,prompt='Enter pixel to be plotted:'
pxl=where(res_pix eq pixel)
read,pltorig,prompt='Plot original spectrum (0=n,1=y)? '
read,pltmod,prompt='Plot original spectrum (0=n,1=y)? '
read,pltres,prompt='Plot original spectrum (0=n,1=y)? '
ans=' '
read,ans,prompt='Also direct spectra to ASCII file (y/n)? '
window,17,title='Specific Pixel Plots'
if ans eq 'y' then begin
        pltname='     '
        read,pltname,prompt='Enter filename:'
        specplot,pxl,orig=pltorig,model=pltmod,res=pltres,filename=pltname,$
        win=17
endif else $
        specplot,pxl,orig=pltorig,model=pltmod,res=pltres,win=17
return
end

;==========================================================================
pro specplot,pixel,orig=orig,model=model,res=res,ps=ps,filename=filename,$
        win=win
;+
; NAME:
;   specplot
; PURPOSE:
;   To exract spectra from the pca common block images produced by the
;   various programs written by Tim Titus.  These spectra can come from
;   the original image, the model image, and/or the residual image.
;   Spectra are plotted to the screen in the current active window
;   and can also be plotted to the postscript file idl.ps (which, as per
;   Titus' programs, must be manually closed before printing).  The
;   spectral data can also be output to an ASCII file.
; INPUTS:
;   pixel = the pixel to be plotted; either a single pixel number or a
;           (c,r) ordered pair (which is then turned into a pixel number
;           which corresponds to the column number for Titus' image arrays)
;   filename = name of file for ASCII output.  (Optional input)
; KEYWORDS:
;   /orig = extract spectrum from original image
;   /model = extract spectrum from model image
;   /res = extract spectrum from residual image
;   /ps = create a postscript plot
; OUTPUTS:
;   Plots are output to screen and, optionally a postscript file and ASCII
;   file in the form of a wavelength column and columns of reflectance.
; MODIFICATION HISTORY:
;   Written  David R. Klassen, May 1996
;-
        on_error,1
;
; SET UP COMMON BLOCKS
;
common pca0,subt,tek
common pca1,wave,eigvect,wavetitle
COMMON PCA2,eigvals,new_vecs,PCArr
common pca6,dim_image
COMMON FRACAB, Frac,Frac2
COMMON ENDM, N_EndM,End_M,EM_pix, EM_pix_set, EM_npca
COMMON RES, res_pix,res_pix2
;
no_plt=keyword_set(res)+keyword_set(model)+keyword_set(orig)
if keyword_set(win) then wset,win
set_mp,no_plt,/top
dim_wave=n_elements(wave)
;
N_pix=n_elements(res_pix)
EM_flux=fltarr(N_EndM,dim_wave)
orig_spectra=fltarr(N_pix,dim_wave)
for i=0l,N_pix-1 do orig_spectra(i,*) = PCArr(res_pix(i),*)
recon_spectra=PCArr*0
for i=0,N_EndM-1 do EM_flux(i,*) = PCArr(EM_pix(i),*)
recon_spectra = Frac # EM_flux
res_spectra=orig_spectra-recon_spectra
if n_elements(pixel) eq 0 then pixel=-1
if n_elements(pixel) gt 1 then pix=dim_image(0)*pixel(1)+pixel(0) $
        else pix=pixel
pixx=pix mod dim_image(0)
pixy=pix/dim_image(0)
if pixel eq -1 then begin
        original=rebin(orig_spectra,1,dim_wave)
        modeled=rebin(recon_spectra,1,dim_wave)
        residual=rebin(res_spectra*res_spectra,1,dim_wave)
endif else begin
        original=orig_spectra(pix,*)
        modeled=recon_spectra(pix,*)
        residual=res_spectra(pix,*)
endelse
;
; SET COMMON PLOT KEYWORDS
;
!x.title=wavetitle
!y.title='Reflectance'
!p.subtitle=subt
!x.style=3
tstring=' of Pixel '+strno(pix)+' ['+strno(pixx)+','+strno(pixy)+']'
;
if pixel eq -1 then !p.title='Average Original Spectrum' $
        else !p.title='Original Spectrum'+tstring
if keyword_set(orig) then plot, wave, original
if pixel eq -1 then !p.title='Average Modeled Spectrum' $
        else !p.title='Modeled Spectrum'+tstring
if keyword_set(model) then plot, wave, modeled
if pixel eq -1 then !p.title='Average RMS Residual Spectrum' $
        else !p.title='Residual Spectrum'+tstring
if keyword_set(res) then plot, wave, residual
;
; Check for PostScript plot
;
if keyword_set(ps) then begin
        ps
        if no_plt ge 2 then device, /land
        if keyword_set(orig) then plot, wave, original, $
                title='Original Spectrum'+tstring
        if keyword_set(model) then plot, wave, modeled, $
                title='Modeled Spectrum'+tstring
        if keyword_set(res) then plot, wave, residual, $
                title='Residual Spectrum'+tstring
        x
endif
;
; Check for ASCII file output
if n_elements(filename) ne 0 then begin
        namearr=[wavetitle,'Original','Modeled','Residual']
        fwave=transpose(wave)
        farray=[fwave,original,modeled,residual]
        openw,10,filename
        printf,10,namearr
        printf,10,farray
        close,10
endif
set_mp
return
end