;==========================================================================
PRO coread,profile,Data3,Wavel,WaveT,norm=norm
;PRO coread,profile,Data3,Sky,Wavel,WaveT,norm=norm
;
;Procedure for reading in the data from either .img files or from
;  a vector list of FITS format images.
;If the switch NORM is set, images are normalized to their image mean
;The parameters for the the .dat file are:
;	imagefile is the image filename or the keyword vector for a string
;		vector list of FITS images
;	Dim1,Dim2 are the image dimensions
;	WB is the number of wavelength bands and lgrid is the wavelength grid
;		spacing (normally set to 1 for all wavelengths and 3 to speed
;		things up, taking every third wavelength image).
;	x,y,t,b margins are the left,right, top, and bottom margins to be
;		excluded; xygrid is the x,y sampling grid of the image (normally
;		set to 1 for all or 3 to speed thing up, taking every third pixel
;		in both x and y)
;   badmin and badmax are the intensity limits to valid pixels
;	normal is set to 1 for the hyperspectra to be scaled (to band 1)
;		if so, band 1 is set to 1000.
;	InX,InY are the dimensions of the final, reduced image
;;	Sky are the values of the Sky pixels (0 for *.img images)
;	swapbyte is a flag for byte swapping the data; set to 1 if .img file
;		was created under UNIX but pcatool is running under WIN (and vice
;		versa); else set t 0
;	fltdata is a flag indicating the data format; 1 = floating point, else
;		the data are assumed to be integer
;	WaveT is a title for the wavelength axis in plots, e.g. Wavelength (nm)
;	Wavel is a vector of wavelength values for each image
COMMON GLOBE, WB, InX, InY, XCenter, YCenter, XSize, YSize, XHigh, XLow, $
          YHigh, YLow, High, Low, scl, good_pix, exclude, include
COMMON VECT, vector
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
openr,1,profile
imagefile = '              '
WaveT='           '
readf,1,imagefile
readf,1,Dim1,Dim2
readf,1,WB,lgrid
Wavel=intarr(WB)
readf,1,xmargin,ymargin,tmargin,bmargin,xygrid
readf,1,badmin,badmax,normal
readf,1,swapbytes, fltdata
readf,1,WaveT
for r=0,WB-1 do begin
	readf,1,wave
	Wavel(r)=wave
endfor ;r
close,1
dWB=WB/lgrid
doDim1=Dim1-xmargin-ymargin
doDim2=Dim2-tmargin-bmargin
dDim1=doDim1/xygrid
dDim2=doDim2/xygrid
good_pix=findgen(dDim1*dDim2)
if (fltdata eq 1) then begin
	Data=fltarr(dWB,doDim1,doDim2)
	Data2=fltarr(Dim1,Dim2)
endif else begin
	Data=lonarr(dWB,doDim1,doDim2)
	Data2=lonarr(Dim1,Dim2)
endelse
if imagefile eq 'vector' and n_elements(vector) gt 1 then begin
	l=0
	for v=0,WB-1,lgrid do begin
		Data2=readfits(vector(v))
		if (swapbytes eq 1) then byteorder,Data2,/LSWAP
		Data(l,*,*)=Data2(xmargin:Dim1-1-ymargin,bmargin:Dim2-1-tmargin)
		if keyword_set(norm) then Data(l,*,*)=Data(l,*,*)/mean(Data(l,*,*))
		l=l+1
		for vn=1,lgrid-1 do data2=readfits(vector(vn))
	endfor ;v
endif else if imagefile eq 'vector' and n_elements(vector) eq 1 then begin
	filenamevect=strarr(WB)
	openr,99,vector
	readf,99,filenamevect
	close,99
	l=0
	for v=0,WB-1,lgrid do begin
		Data2=readfits(filenamevect(v))
		if (swapbytes eq 1) then byteorder,Data2,/LSWAP
		Data(l,*,*)=Data2(xmargin:Dim1-1-ymargin,bmargin:Dim2-1-tmargin)
		if keyword_set(norm) then Data(l,*,*)=Data(l,*,*)/mean(Data(l,*,*))
		l=l+1
		for vn=1,lgrid-1 do data2=readfits(vector(vn))
	endfor ;v
endif else if imagefile eq 'cube' then begin
	l=0
	dl=0
	Data2=readfits(vector)
	if (swapbytes eq 1) then byteorder,Data2,/LSWAP
	for v=0,WB-1,lgrid do begin
		Data(l,*,*)=Data2(xmargin:Dim1-1-ymargin,bmargin:Dim2-1-tmargin,dl)
		if keyword_set(norm) then Data(l,*,*)=Data(l,*,*)/mean(Data(l,*,*))
		l=l+1
		dl=dl+lgrid
	endfor ;v
endif else begin
	openr,1,imagefile
	l=0
	for ii=0,WB-1,lgrid do begin
		readu,1,Data2
		if (swapbytes eq 1) then byteorder,Data2,/LSWAP
		Data(l,*,*)=Data2(xmargin:Dim1-1-ymargin,bmargin:Dim2-1-tmargin)
		if keyword_set(norm) then Data(l,*,*)=Data(l,*,*)/mean(Data(l,*,*))
		l=l+1
		for ij=1,lgrid-1 do readu,1,Data2
	endfor ;jj
	close,1
endelse
WB=dWB
Data=rebin(Data,l,dDim1,dDim2)
;Sky = fltarr(WB)
if (normal eq 1) then begin
	for j=1,WB-1 do Data(j,*,*)=Data(j,*,*)/Data(0,*,*)
	Data=Data*1000
endif
good=make_array(Dim1,Dim2,value=1.0)
if n_elements(exclude) gt 1 and exclude(0) ne -1 then good(exclude)=0.0
if include(0) ne -1 then good(include)=1.0
good=good(xmargin:Dim1-1-ymargin,bmargin:Dim2-1-tmargin)
good=rebin(good,dDim1,Ddim2)
for i = 0,WB-1 do begin
	bad=where(Data(i,*,*) le badmin or Data(i,*,*) ge badmax)
;	if (n_elements(bad) gt 1) then good(bad)=-1.
	if (n_elements(bad) gt 1) then good(bad)=0.0
	Data(i,*,*)=Data(i,*,*)*good
endfor ;i
good_pix=where(good eq 1.0)
InX = dDim1
InY = dDim2
print,InX,InY
Data3=Data
;
return
end

;==========================================================================
Pro Rubber_box,color,wid,x,y,InX,InY
cursor,x0,y0,/dev,/DOWN
XC = (float(x0) / !D.X_SIZE)*InX + 1
YC = (float(y0) / !D.Y_SIZE)*InY + 1
x0 =  long(XC)
y0 = long(YC)
cursor,x1,y1,/dev,/DOWN
XC = (float(x1) / !D.X_SIZE)*InX + 1
YC = (float(y1) / !D.Y_SIZE)*InY + 1
x1 =  long(XC)
y1 =  long(YC)
x = (x0 + x1)/2.
y = (y0 + y1)/2.
wid = [x1-x0, y0-y1]
print,!D.X_SIZE,!D.Y_SIZE,InX,InY
print,x0,y0,x1,y1
print,wid,x,y
tvbox,wid,x,y,color ;in pcalib3.pro
return
end

;==========================================================================
pro tvbox,width,x,y,color,THICK=thick
;+
; NAME:
;   TVBOX
; PURPOSE:
;   Draw a box(es) or rectangle(s) of specified width centered on the
;   cursor, or at specified position(s).
; CALLING SEQUENCE:
;   TVBOX,WIDTH	        ;Draw box of width WIDTH @ current cursor position
;   TVBOX,WIDTH,X,Y     ;Draw box at position X,Y
;   TVBOX,WIDTH,X,Y,COLOR,THICK = thick      ;Color, thickness options
; INPUTS:
;   WIDTH -  either a scalar giving the width of a box, or a 2 element
;            vector giving the length and width of a rectangle.
; OPTIONAL INPUTS:
;   X  -  x position for box center, scalar or vector
;   Y  -  y position for box center, scalar or vector.   If vector, then Y
;         must have the same number of elements as X
;         If X and Y are not specified then program will draw
;         box at current cursor position
;   COLOR - color intensity (0-255).  Default=!P.COLOR
; OUTPUTS:
;   None
; OPTIONAL KEYWORD INPUT:
;   THICK - scalar specifying thickness of the drawn box, default = 1
; SIDE EFFECTS:
;   A box or rectangle will be drawn on screen
;   For best results WIDTH should be odd.  (If WIDTH is even, the actual
;   size of the box will be WIDTH + 1.)
; RELATED PROCEDURES:
;   The color of a box in the graphics overlay can be set with GRCOL.
; RESTRICTIONS:
;   TVBOX does not check whether box is off the edge of the screen
; REVISON HISTORY:
;   Written, W. Landsman   STX Co.           10-6-87
;   Modified to take vector arguments. Greg Hennessy Mar 1991
;   Fixed centering of odd width    W. Landsman    Sep. 1991
;-
On_error,2
npar = N_params()                         ;Get number of parameters
if (npar LT 1) then begin
     print,'CALLING SEQUENCE - tvbox,width,[x,y,color,THICK = ]'
     return
endif
zparcheck,'TVBOX',width,1,[1,2,3,4,5],[0,1],'Box Width' ;in pcalib3.pro
if ( N_elements(width) EQ 2 ) then w = width/2. else w = [width,width]/2.
if ( npar LT 3 ) then begin
    cursor,x,y,/DEVICE,/NOWAIT    ;Get unroamed,unzoomed position
    if (x LT 0) or (y LT 0) then begin
       message,'Position cursor in window ' + strtrim(!D.WINDOW,2) + $
              ' -- then hit mouse button',/INF
       cursor,x,y,/DEVICE,/WAIT
    endif
endif
if npar LT 4 then color = !P.COLOR              ;Set color index.
if not keyword_set(THICK) then thick = !P.THICK
;color = color mod 256                      ;Make sure the color index is legal.
nbox = N_elements(x)                      ;Number of boxes to draw
if ( nbox NE N_elements(Y) ) then $
       message,'ERROR - X and Y positions must have same number of elements'
xs = fix(x+0.5)  &  ys = fix(y+0.5)
for i=0,nbox-1 do begin
  xt = xs(i) + fix( [1,1,-1,-1,1]*w(0) )      ;X edges of rectangle
  yt = ys(i) + fix( [-1,1,1,-1,-1]*w(1) )     ;Y edges of rectangle
  plots, xt, yt, /DEVICE,COLOR=color,THICK=thick  ;Plot the box on the window.
endfor
return
end

;==========================================================================
pro zparcheck,progname,parameter,parnum,types,dimens,message
;+
; NAME:
;     ZPARCHECK
; PURPOSE:
;	Routine to check user parameters to a procedure
; CALLING SEQUENCE:
;	zparcheck,progname,parameter,parnum,types,dimens,message
; INPUTS:
;	progname  - string name of calling procedure
;	parameter - parameter passed to the routine
;	parnum    - integer parameter number
;	types     - integer scalar or vector of valid types
;		 1 - byte        2 - integer  3 - int*4
;		 4 - real*4      5 - real*8   6 - complex
;		 7 - string      8 - structure
;	ndimens   - integer scalar or vector giving number
;		      of allowed dimensions.
;	message - (optional string message to be printed if an
;		      error is found).
; OUTPUTS:
;	none
; SIDE EFFECTS:
;	If an error in the parameter is a message is printed
;	a RETALL issued
; HISTORY
;	version 1  D. Lindler  Dec. 86
;	documentation updated.  M. Greason, May 1990.
;-
;----------------------------------------------------------
;
; convert types and ndimens to vectors if scalars supplied
;
  vtypes = intarr( N_elements(types) ) + types
  vndims = intarr( N_elements(dimens) ) + dimens
;
; get type and size of parameter
;
  s = size(parameter)
  ndim = s(0)
  type = s(ndim+1)
;
; check if parameter defined.
;
  if type eq 0 then begin
	err=' is undefined.'
	goto,abort
  endif
;
; check for valid dimensions
;
  valid = where(ndim eq vndims,nvalid)
  if nvalid lt 1 then begin
	err = 'has wrong number of dimensions'
	goto,abort
  endif
;
; check for valid type
;
  valid = where(type eq vtypes,ngood)
  if ngood lt 1 then begin
	err = 'is an invalid data type'
	goto,abort
  endif
;
  return
;
; bad parameter
;
abort:
  mess = ' '
  if N_params() lt 6 then message=''
  if message ne '' then mess=' ('+message+') '
  print,string(7b) + 'Parameter '+strtrim(parnum,2)+mess,$
	' of routine ',strupcase(progname)+' ',err
  sdim = ' '
  for i = 0,N_elements(vndims)-1 do begin
	if vndims(i) eq 0 then sdim=sdim+'scalar' $
			  else sdim=sdim+string(vndims(i),'(i3)')
  end
  print,'Valid dimensions are:'+sdim
;
  stype = ' '
  for i = 0,N_elements(vtypes)-1 do begin
	case vtypes(i) of
		1: stype = stype + ' byte'
		2: stype = stype + ' integer'
		3: stype = stype + ' longword'
		4: stype = stype + ' real*4'
		5: stype = stype + ' real*8'
		6: stype = stype + ' complex'
		7: stype = stype + ' string'
                8: stype = stype + ' structure'
	endcase
  endfor
  print,'Valid types are:'+stype
  ;if !debug then stop
  retall  ; zparcheck
  end

;==========================================================================
pro pcainit,aa,wavl,tek=tek,xtitle=xt,subtitle=subt,norm=norm,object=object,$
	scale=scale,dim_image=dpix
;
;Written by: T.N.Titus, Univ of Wyo, 1992
;
;Propose:
;       This is the PCA Initialization Routine.
;       This routine must be run before other PCA Routine will operate.
;
;Input:
;       aa:     The image that is to be operated on.  The format must
;               be such that columns are objects or samples, and
;               rows are the coordinates or wavelengths.
;Optional Input:
;       wavl:   wavelength
;       tek:    if tek is set, then all pca routines will plot in tek.
;       xtitle: the title for the x axis in all future plots.
;       subtitle: defines subtitle for all future pca plots.
;       norm:   this will normalize each coord to its st dev before
;               performing PCA on the data set.
;       object: Default is all objects else only objects listed will
;               be used to calculate the variance, ie the evecs.
;       scale: normalizes the intensity of each column (object). A
;               mean col value less than one is set to one.
;       dim_image: a 2 componet vector containing x,y image size for
;                       use with 2-D spacial images.
;
COMMON PCA0, sub_title, tek_term
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA3, no_obj
COMMON PCA4, image, rot_image
COMMON PCA5, objects
COMMON PCA6, dim_image
COMMON PCA7, cm_vector
;
image=aa
fnddim,aa,c,r
if n_elements(dpix) ne 0 then dim_image=dpix else dim_image=[c,1]
if dim_image(0)*dim_image(1) ne c then message,'*** Error in image size ***'
if n_elements(object) eq  0 then no_obj=(c<r) else no_obj=n_elements(object)
;if n_elements(xt) ne 0 then wave_title=xt else wave_title='n''th Coordinate'
if n_elements(object) ne 0 then objects=object else objects=$
   indgen(c)
if n_elements(subt) ne 0 then  sub_title=subt else sub_title='PCA Plot'
if n_elements(tek) eq 0 then tek=(!d.name eq 'TEK')
tek_term=tek
if keyword_set(scale) then aaa=aa/(transpose(rmnarr(transpose(aa)))>1.) $
	else aaa=aa
var=varcalc(aaa,cc,cm,norm=norm,object=object) ;in pcalib3.pro
eigsol,var,a,lam ;in pcalib3.pro
print,'Eigenvalues'
printf,5,'Eigenvalues'
maxlam=n_elements(wavelength)-1<10
print, lam[0:maxlam]
printf,5,lam[0:maxlam]
sumeign = total(lam)
cumvar = lam/sumeign
print,'Relative Variance per Axis'
printf,5,'Relative Variance per Axis'
print, cumvar[0:maxlam]
printf,5,cumvar[0:maxlam]
;if n_elements(wavl) eq 0 then wavl=indgen(n_elements(lam))
xx=cc#a
rot_image=cc
new_vecs=xx
eigvals=lam
;wavelength=wavl
eigvecs=a
cm_vector=cm
return
end
;==========================================================================
pro fnddim,a,c,r
;finds the dimensions of a 2-d array
r=n_elements(a(0,*))
c=n_elements(a(*,0))
return
end
;==========================================================================
function mean,a
ave=total(a)/n_elements(a)
return,ave
end
;==========================================================================
function varcalc,aa,cc,cm_vect,norm=norm,object=object
;
;Purpose:  To calculate the variance matrix of the rows of a matrix.
;
;	Input:
;		aa:  the data matrix.
;
;	Keywords:
;		norm:  if set, the variance is normalized to one.
;		object:  if set,  the variance is calculated only using
;			the columns contained in the vector.
;		cm_vect: the vector that contains the center of mass.
;
;	Output:  The variance matrix.
;
if n_elements(object) eq 0 then bb=rmnarr(aa) $
   else bb=rmnarr(aa(object,*),n_col=n_elements(aa(*,0)))
if n_elements(object) eq 0 and keyword_set(norm) then sigma=rsdarr(aa)
if n_elements(object) ne 0 and keyword_set(norm) $
   then sigma=rsdarr(aa(object,*),n_col=n_elements(aa(*,0)))
if not keyword_set(norm) then sigma=1.0
pos=where(sigma eq 0)
if pos(0) ne -1 then sigma(pos)=1.0
cm_vect=bb(0,*)
cc=(aa-bb)/sigma
; the following line is the beginnings of a test to do PCA across
; the spatial dimensionality instead of the wavelength dimensionality
;cc=transpose(cc)
if n_elements(object) eq 0 then var=transpose(cc)#cc else var=$
   transpose(cc(object,*))#cc(object,*)
return,var
end
;==========================================================================
pro eigsol,var,aa,lam
a=eigvec(var,d) ;in pcalib3.pro
pos=sort(-d)
lam=d(pos)
aa=a(*,pos)
return
end
;==========================================================================
function eigvec,a,dd
aa=a
trired,aa,d,e
;tred2,aa,d,e
triql,d,e,aa
;tqli,d,e,aa
if n_params() eq 2 then dd=d
return,aa
end
;==========================================================================
function rmnarr,a,n_col=n_col
;
;Purpose: To calculate a matrix where each element is the ave of the row.
;
;Input:
;	a: the matrix to be operated on.
;
;Optional Input:
;	n_col:  to specify the number of columns of output matrix.
;		The default is the same size as a.
;
fnddim,a,c,r ;in pcalib3.pro
if n_elements(n_col) ne 0 then c=n_col
mu=fltarr(c,r)
for i=0,r-1 do mu(*,i)=mean(a(*,i)) ;in pcalib3.pro
return,mu
end
;==========================================================================
function rsdarr,a,n_col=n_col
;
;Purpose: To calculate a matrix whose elements have the st dev
;	of the row of matrix a.
;
fnddim,a,c,r
if n_elements(n_col) ne 0 then c=n_col
mu=fltarr(c,r)
for i=0,r-1 do mu(*,i)=stdev(a(*,i))
return,mu
end

;==========================================================================
function newcoo,coords
;
;Purpose: To take coordinates in "normal" space and
;	transform them into pca space.
;
;Input:
;	coords:  N objs x M dims matrix.
;
;Output:  N obj x M dim matrix.
;
COMMON PCA1, wavelength, eigvecs, wave_title
COMMON PCA7, cm_vector
;
fnddim,coords,col,row ;in pcalib3.pro
cm_arr=rebin(cm_vector,col,row)
xx=(coords-cm_arr)#eigvecs
;
return,xx
end
;==========================================================================
pro pcaplot,i,j,k,sym=sym,title=titl,ps=ps,iter=iter,outfile=outfile,$
	label=label,xrange=xrange,yrange=yrange,all=all,zoomplt=zoomplt,$
	thrd=thrd
;
;PCA j'th Eigenvector vs i'th Eigenvector Scatter Plot.
;
COMMON PCA0, sub_title, tek_term
COMMON PCA2, eigvals, new_vecs, PCArr
COMMON PCA5, objects
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
;
!p.multi=0
!p.charsize=1
xr=xrange
yr=yrange
n_data=n_elements(new_vecs(*,0))
if keyword_set(all) then list=lindgen(n_data) else list=objects
n_list=n_elements(list)
if n_elements(new_vecs) eq 0 then print,'PCA not Init!'
if n_elements(new_vecs) eq 0 then return
xt=strith(i)+' Eigenfunction' ;in pcalib3.pro
yt=strith(j)+' Eigenfunction' ;in pcalib3.pro
x=new_vecs(list,i)
y=new_vecs(list,j)
if n_elements(k) eq 1 then begin
	zt=strith(k)+' Eigenfunction' ;in pcalib3.pro
	z=new_vecs(list,k)
	zrange=minmax(z)
endif
;if n_elements(xrange) eq 0 then xrange=[min(new_vecs(list,i)), $
;   max(new_vecs(list,i))]
;if n_elements(yrange) eq 0 then yrange=[min(new_vecs(list,j)), $
;   max(new_vecs(list,j))]
if keyword_set(thrd) then window,19,xs=400,ys=400,title='PCA Plot 3',retain=2
;if keyword_set(thrd) then window,19,xs=600,ys=400,title='PCA Plot 3',retain=2
case n_elements(k) of
	0: begin
		k_col=bytscl(new_vecs(*,0),top=6)+1
		k_col(list)=255
	end
	1: begin
		k_title=strith(k) +' Eigenfunction'
		k_col=bytscl(new_vecs(*,k),top=6)+1
		k_col(list)=bytscl(new_vecs(list,k),top=6)+1
	end
	else: begin ; for coloring individual points!
		if n_elements(k_title) eq 0 then k_title='Parameter'
		k_col=k
;		k_col=bytscl(k,top=6)+1
;		k_col(list)=bytscl(k(list),top=6)+1
	end
endcase
if n_elements(titl) eq 0 then titl=yt+' vs '+xt
if n_data ge 1000. then def_sym=3 else def_sym=1
;if n_elements(sym) eq 0 then sym=fltarr(n_data)+def_sym
if n_elements(sym) eq 0 then sym=def_sym
if n_elements(label) eq 1 then label=strno(indgen(n_data)) ;in pcalib3.pro
if keyword_set(ps) then begin
	symang=findgen(9)*(!pi*2/8.)
	usersym,0.25*cos(symang),0.25*sin(symang),/fill
	if sym eq 3 then pssym=8 else pssym=sym
	iter = strcompress(string(iter), /remove_all)
	outfile = strcompress(string(outfile), /remove_all)
	filenm='pca'+outfile+'_'+iter+'.eps'
	if keyword_set(zoomplt) then filenm='pca'+outfile+'_'+iter+'_zoom.eps'
	zero = long(0)
	set_plot, 'ps'
	device, /encapsulated
	device, FILENAME = filenm
	device, /color
	device, xsize=20, ysize=20
	colortest=where(k_col eq 255 or k_col eq 7,whitecount)
	if whitecount ne 0 then k_col(colortest)=0
	lct_dmw ;in pcalib1.pro
	if not keyword_set(thrd) then begin
		plot,new_vecs(list,i),new_vecs(list,j),/nodata,xstyle=3,ystyle=3,$
			title=titl, subtitle=sub_title, xtitle=xt, ytitle=yt,$
			xrange=xrange,yrange=yrange,charthick=3,charsize=1.2,$
			xthick=2,ythick=2
		for ii=0l,n_list-1 do oplot,$
			[new_vecs(list(ii),i)],[new_vecs(list(ii),j)],$
				psym=pssym,color=k_col(list(ii))
;				psym=sym(list(ii)),color=k_col(list(ii))
		if n_elements(label) ne 0 then for ii=0,n_list-1 do xyouts,$
			new_vecs(list(ii),i),new_vecs(list(ii),j),label(list(ii))
	endif
	if keyword_set(thrd) then begin
		SURFACE, DIST(5), /NODATA, /SAVE, XRANGE=xr,YRANGE=yr,$
			ZRANGE=zrange,XSTYLE=4,YSTYLE=4, ZSTYLE=4, CHARSIZE=1.5, $
       			POSITION=[0.1, 0.1, 0.95, 0.95, 0.1, 0.95], $
       			xtitle=xt,ytitle=yt,ztitle=zt,charthick=3, $
					xthick=5, ythick=5, zthick=5
		Triangulate, x, y, triangles, boundaryPts
;		gridSpace = [0.01, 0.05]
;		griddedData = TriGrid(x, y, z, triangles, gridSpace, $
		griddedData = TriGrid(x, y, z, triangles, $
			XGrid=xvector, YGrid=yvector)
		SHADE_SURF, griddeddata, Shades=BYTSCL(griddeddata, Top=6)+1,$
			xvector,yvector, xtitle=xt, ytitle=yt, ztitle=zt, /T3, $
			title='3D PCA Plot: 2 vs. 1 vs. 0',xrange=xr,yrange=yr
		AXIS, /XAxis,0,0,0, /T3D,CharSize=1.5,Color=2,xtitle=xt
		AXIS, /YAxis,0,0,0, /T3D,CharSize=1.5,Color=2,ytitle=yt
		AXIS, /ZAxis,0,0,0, /T3D,CharSize=1.5,Color=2,ztitle=zt
;		PLOTS, x, y, z, PSYM=3, /T3D
	endif
	if whitecount ne 0 then k_col(colortest)=255
	device, /close
endif
x ;in pcalib3.pro
if n_elements(win) ne 0 then wset,win
lct_dmw ;in pcalib1.pro
if not keyword_set(thrd) then begin
	plot,new_vecs(list,i),new_vecs(list,j),/nodata,xstyle=3,ystyle=3,$
		title=titl, subtitle=sub_title, xtitle=xt, ytitle=yt,$
		xrange=xrange,yrange=yrange,charthick=2,charsize=1.2
	for ii=0l,n_list-1 do oplot,$
		[new_vecs(list(ii),i)],[new_vecs(list(ii),j)],psym=sym,$
;		[new_vecs(list(ii),i)],[new_vecs(list(ii),j)],psym=sym(list(ii)),$
		color=k_col(list(ii))
	if n_elements(label) ne 0 then for ii=0,n_list-1 do xyouts,$
		new_vecs(list(ii),i),new_vecs(list(ii),j),label(list(ii))
endif
if keyword_set(thrd) then begin
;	window,19,xs=600,ys=400,title='PCA Plot 3',retain=2
	SURFACE, DIST(5), /NODATA, /SAVE, XRANGE=xr,YRANGE=yr,$
		ZRANGE=zrange,XSTYLE=4,YSTYLE=4, ZSTYLE=4, CHARSIZE=1.5, $
       		POSITION=[0.1, 0.1, 0.95, 0.95, 0.1, 0.95], $
       		xtitle=xt,ytitle=yt,ztitle=zt
	Triangulate, x, y, triangles, boundaryPts
;	gridSpace = [0.01, 0.05]
;	griddedData = TriGrid(x, y, z, triangles, gridSpace, $
	griddedData = TriGrid(x, y, z, triangles, $
		XGrid=xvector, YGrid=yvector)
	SHADE_SURF, griddeddata, Shades=BYTSCL(griddeddata, Top=6)+1,xvector, $
		yvector, xtitle=xt, ytitle=yt, ztitle=zt, /T3, $
		title='3D PCA Plot: 2 vs. 1 vs. 0',xrange=xr,yrange=yr
	AXIS, /XAxis, 0, 0, 0, /T3D, CharSize=1.5, Color=2,xtitle=xt
	AXIS, /YAxis, 0, 0, 0, /T3D, CharSize=1.5, Color=2,ytitle=yt
	AXIS, /ZAxis, 0, 0, 0, /T3D, CharSize=1.5, Color=2,ztitle=zt
;	PLOTS, x, y, z, PSYM=3, /T3D
endif
return
end
;==========================================================================
function strno,n,remove=remove
if keyword_set(remove) then s=strcompress(string(n),/rem) $
   else s=strcompress(string(n))
return,s
end
;==========================================================================
function strith,n
; produces labels for pcaplots
if n_elements(n) gt 1 then return,strno(n)
th=strarr(20)+'th'
th(1)='st'
th(2)='nd'
th(3)='rd'
if n lt 20 then m=n else m=n mod 10
nth=strcompress(string(n))+th(m)
return,nth
end
;==========================================================================
PRO r_endm,a,nendm,nbands
;used by PCABAT.PRO
;reads endmember spectra from open unit 1
;user must open and close unit 1 outside subroutine
; NO INPUT
; OUTPUT:
; nendm= number of endmembers read
; nbands= number of wavelengths
; a=array of results
readf,1,nendm,nbands
print,nendm,nbands
a=fltarr(nbands,nendm)
for i=0,nendm-1 do begin
	for ii = 0,nbands-1 do begin
		readf,1,w,x
		a(ii,i)=x
	endfor ;ii
endfor ;i
return
end
;==========================================================================
pro x
; returns plotting to xterminal
set_plot,'x'
return
end
;==========================================================================
pro ps,land_scape=land
;
;Pupose: to set plot name to post script.
;
;Keyword:
;	land_scape: Set output to landscape, [Default: Portrait]
;
set_plot,'ps'
if keyword_set(land) then device,/land else device,/port
return
end