;==========================================================================
PRO in_data,datafile=df,norm=norm
;
;This procedure asks the user for a filename containing the list of images,
;preferably in order of wavelength.  Images should bin in one .img file, a
;headerless, binary format with each image written sequentially into the file.
;The images are then read into the array Data1.
;
;     Data1    =  Array containing image.
;     NImages  =  Number of images (Wavebands).
;
;If keyword NORM is set, then each images is normalized by its mean
;
COMMON PCA1, wavelength, eigvecs, wave_title
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 WIND, clean, w0x, w1x, w2x, w0y, w1y, w2y
;COMMON SPEC, Sky, Sky_pca, Sky_pca2
COMMON VECT, vector
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
;Reading in name list of files containing data at various wavebands.
if n_elements(df) lt 1 then begin
	profile = '                '
	read,'Profile name? ',profile
endif else begin
	profile=df
endelse
;coread, profile, Data1,Sky,wavelength,wave_title,$
coread, profile, Data1,wavelength,wave_title,$
	norm=keyword_set(norm) ;in pcalib3.pro
NImages=WB
return
end

;==========================================================================
PRO resize_im
;
;Interactive method of selecting the size of the images to be run through
;PCAINIT.
;
;     width  =  array containing X and Y dimensions of box.
;     X      =  X position for the center of the box.
;     Y      =  Y position for the center of the box.
;     XSize  =  the final resized X dimension of the image.
;     YSize  =  the final resized Y dimension of the image.
;     tvbox  =  display routine that draws box.
;
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 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
;
;Entering center point for image box.
tst = 0
while (tst ne 1) do begin
	print
	print, '====================================================='
	print, 'Click mouse, first to choose UL, then LR of image.'
	print, '====================================================='
	device, /cursor_crosshair
	lct_dmw ;in pcalib1.pro
	color = 1
	rubber_box,color,wid,X,Y,InX*Scl,InY*Scl ;in pcalib3.pro
	Choice = 0
	xmenu, ['REDEFINE IMAGE BOX',$
	        'CONTINUE'], BASE=base, BUTTONS=B, title='Resize Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = Choice, event.id
	WIDGET_CONTROL, event.top, /DESTROY
	case Choice of
	;Redefining image box
	0: begin
		tst = 0
		tvscl, Image < High(Im_Num) > Low(Im_Num)
	end ;0
	;Continue with initialization.
	1:  begin
		tst = 1
		XSize = wid(0) / scl
		YSize = wid(1) / scl
	end ;1
	endcase
endwhile
;
XSize = long(XSize)
YSize = long(YSize)
XCenter = X/scl                   ;modified 3/30/94
YCenter = Y/scl                   ;modified 3/30/94
X2 = XSize
Y2 = YSize
if (XSize ge InX) then begin
	XHigh = InX - 1
	XLow = 0
	XSize = InX
endif else begin
	X2 = XSize / 2
	XHigh = XCenter + X2
	XLow = XHigh - XSize + 1
endelse
if (YSize ge InY) then begin
	YHigh = InY - 1
	YLow = 0
	YSize = InY
endif else begin
	Y2 = YSize / 2
	YHigh = YCenter + Y2
	YLow = YHigh - YSize + 1
endelse
res_pix=make_array(Xsize*Ysize,/LONG,/INDEX)
good_pix=make_array(Xsize*Ysize,/LONG,/INDEX)
;
return
end

;==========================================================================
PRO get_box, color, boxsz
;
;This procedure displays a box on an image or a PCA Plot and determines a
;box size and position so that specific pixels or eigenpoints may be
;extracted.
;
;     win    =  Window box is to be displayed in.
;     boxsz  =  X and Y size of the box displayed [X,Y].
;     width  =  array containing X and Y dimensions of box.
;     X_im   =  X position for the center of the box (image).
;     Y_im   =  Y position for the center of the box (image).
;     X_pca  =  X position for the center of the box (PCA Plot).
;     Y_pca  =  Y position for the center of the box (PCA Plot).
;     tvbox  =  display routine that draws box.
;
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 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
;
tst1 = 0
;Entering center point for box.
while (tst1 eq 0) do begin
;	if (win lt 3) then begin
		if (win eq 0) then begin
			Element1 = Elem1
			Element2 = Elem2
			!x=w0x
			!y=w0y
		endif
		if (win eq 1) then begin
			Element1 = Elem3
			Element2 = Elem4
			!x=w1x
			!y=w1y
		endif
		if (win eq 2) then begin
			Element1 = Elem5
			Element2 = Elem6
			!x=w2x
			!y=w2y
		endif
;	endif
	if (clean(3) eq 1) then begin
		wset,3
		tvscl, Image < High(Im_Num) > Low(Im_Num)
	endif
	wset,win
	print
	print, '=================================================='
	print, 'Click mouse, first to choose UL, then LR of image.'
	print, '=================================================='
	device, /cursor_original
	if (win eq 3) then device, /cursor_crosshair
	lct_dmw ;in pcalib1.pro
	cursor,x0,y0,/dev,/DOWN
	cursor,x1,y1,/dev,/DOWN
	X = (x0+x1)/2.
	Y = (y0+y1)/2.
	width = [x1-x0,y0-y1]
	tvbox,width,X,Y,color ;in pcalib3.pro
;
	Chose = 0
	xmenu, ['REDEFINE IMAGE BOX','DISPLAY Different PCA Plot', $
		'DISPLAY Different Image','CONTINUE'], BASE=base, BUTTONS=B, $
		title='Pixel Display Menu'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = Chose, event.id
	WIDGET_CONTROL, event.top, /DESTROY
;
    case Chose of
	;Changing center point for box.
	0: begin
		tvbox,width,X,Y,0 ;in pcalib3.pro
		tst1 = 0
	end ;0
	1: begin
		print
		d_pcaplot ;in pcalib1.pro
		tst1 = 0
	end ;1
	2: begin
		XDim = [XLow,XHigh]
		YDim = [YLow,YHigh]
		d_image,Xdim,Ydim,scl ;in pcalib1.pro
		tst1 = 0
	end ;2
	;Continue with pixel coloring.
	3: begin
		tst1 = 1
		if (win eq 3) then begin
			X_im = X
			Y_im = Y
			print,' Box Center coordinates = ',X_im,Y_im
		endif else begin
			X_pca = X
			Y_pca = Y
			print,' Box Center coordinates = ',X_pca,Y_pca
		endelse
	end ;3
	endcase
endwhile
;
;tvbox,width,X,Y,0 ;in pcalib3.pro
return
end

;==========================================================================
PRO d_c_pca,new_color
;
;This procedure takes box information from get_box and image_pix to display
;colored pixels on the PCA Plots from the image frame.
;
;     X_pcapl# =  X PCA values for selected pixels in the image frame.
;     Y_pcapl# =  Y PCA values for selected pixels in the image frame.

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 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
;
refr_pca,/color ;in pcalib1.pro
X_pca0 = new_vecs(c_pix,Elem1)
Y_pca0 = new_vecs(c_pix,Elem2)
X_pca1 = new_vecs(c_pix,Elem3)
Y_pca1 = new_vecs(c_pix,Elem4)
X_pca2 = new_vecs(c_pix,Elem5)
Y_pca2 = new_vecs(c_pix,Elem6)
;
numb1 = n_elements(c_pix)
wset,0
!x = w0x
!y = w0y
for i = 0l,numb1-1l do begin
	oplot,[X_pca0(i)],[Y_pca0(i)],psym=1,color=new_color
endfor
wset,1
!x = w1x
!y = w1y
for i = 0l,numb1-1l do begin
	oplot,[X_pca1(i)],[Y_pca1(i)],psym=1,color=new_color
endfor
wset,2
!x = w2x
!y = w2y
for i = 0l,numb1-1l do begin
	oplot,[X_pca2(i)],[Y_pca2(i)],psym=1,color=new_color
endfor
;
clean(*)= 0
;
return
end

;==========================================================================
PRO s_endm
;
;Procedure for interactively selecting end-member spectra from either the
;original image, or PCA Plots.  End-members are then used to determine
;fractional abundance information.
;
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 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, v_names, StoN
;
refr_pca,/color ;in pcalib1.pro
xmenu, ['2 Endmembers','           3 Endmembers           ','4 Endmembers', $
	'5 Endmembers'], BASE=base,BUTTONS=B,title='Number of Endmemebers?'
WIDGET_CONTROL,base,/REALIZE
event = WIDGET_EVENT(base)
WIDGET_CONTROL, GET_UVALUE = Choice, event.id
WIDGET_CONTROL, event.top, /DESTROY
;
case Choice of
	0: N_EndM=2
	1: N_EndM=3
	2: N_EndM=4
	3: N_EndM=5
endcase
;print
;print, 'What is the number of end-members? '
;read, N_EndM
v_names=strarr(5)
print, 'The number of end-members is', N_EndM
End_M = fltarr(N_EndM-1,N_EndM)
;
xmenu, ['SELECT Endmembers from PCA plot',' SELECT Endmembers from IMAGE ', $
	'  Enter Endmember Pixles at KEYBOARD  ',$
	'Enter Endmember PCA Coords at KEYBOARD', 'EXIT End-member selection'],$
	BASE=base, BUTTONS=B, title='Endmember Selction Menu'
WIDGET_CONTROL,base,/REALIZE
event = WIDGET_EVENT(base)
WIDGET_CONTROL, GET_UVALUE = Choice, event.id
WIDGET_CONTROL, event.top, /DESTROY
;
case Choice of
0: begin
	Chose_pca = 1
	s_pca_em2 ;in pcalib2.pro
end ;0
1: begin
	Chose_im = 1
	s_im_em2 ;in pcalib2.pro
end ;1
2: begin
	Chose_kb = 1
	s_kb_em2 ;in pcalib2.pro
end ;2
3: begin
	Chose_pcakb = 1
	s_pcakb_em2 ;in pcalib2.pro
end ;3
else:
endcase
;
return
end

;==========================================================================
PRO s_pca_em2
;
;Procedure for interactively selecting end-member spectra from PCA Plot.
;This procedure uses a box to extract the data and then takes the median in
;order to find the resulting end-member.
;
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 PCA7, cm_vector
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, v_names, StoN
;
EM_pix = intarr(N_EndM)
EM_pix = long(EM_pix)
EM_pix(0) = -1
EM_pix_set=make_array(5,XSize*YSize,/long,value=(-1))
;for j = 0,N_EndM-2 do begin
;  End_M(j,0) = Sky_pca(j)
;endfor
Left = 1
i = 0
XBox = 24
YBox = 24
color=1
while (i lt N_EndM) do begin
	s_pcaplot ;in pcalib1.pro
	print
	print, '==================================================='
	print, 'Click LEFT button on PCA Plot to choose end-member.'
	print, '==================================================='
	device, /cursor_original
	cursor, X, Y, /device
	if(!err eq Left) then begin
		lct_dmw ;in pcalib1.pro
		width=[XBox,YBox]
		tvbox,width,X,Y,color ;in pcalib3.pro
	endif
	width(0) = XBox
	width(1) = YBox
	X_pca = X
	Y_pca = Y
	pca_pix, ZC ;in pcalib1.pro
	Numb_pix = n_elements(c_pix1)
	if c_pix1[0] ne -1 then begin ; That is, if the pixel is in the image set
		med1 = MEDIAN(new_vecs(c_pix1,Element1))
		med2 = MEDIAN(new_vecs(c_pix1,Element2))
		res1 = fltarr(Numb_pix)
		for kk = 0l,Numb_pix-1l do $
			res1(kk) = (new_vecs(c_pix1(kk),Element1)-med1)^2 + $
			(new_vecs(c_pix1(kk),Element2)-med2)^2
		tmp1 = MIN(res1,elm)
		EM_pix(i) = c_pix1(elm)
	endif else begin ; That is, if pixel is not in the image set
		k=n_oem_pix+i
		zvec=fltarr(WB)
		zvec[Element1]=ZC[0]
		zvec[Element2]=ZC[1]
		print,'PCA vector: ',zvec
		ans=' '
		read,ans,prompt='Would you like to add values to any other dimensions (y/n)? '
		while ans eq 'y' do begin
				print,'PCA vector: ',zvec
				read,l,zl,prompt='Enter dimension number and value: '
				zvec[l]=zl
				read,ans,prompt='Again (y/n)? '
		endwhile
;		new_vecz = zvec ## eigvecs
;		new_vecs[exclude[k],*]=new_vecz
		new_vecs[exclude[k],*]=zvec
;		endmz=cm_vector*(invert(eigvecs) ## transpose(new_vecz))
		endmz=cm_vector*(invert(eigvecs) ## transpose(zvec))
		c_pix1=exclude[k]
		EM_pix[i]=exclude[k]
		PCArr[exclude[k],*]=endmz
	endelse
	EM_pix_set(i,0:n_elements(c_pix1)-1)=c_pix1
	for j = 0,N_EndM-2 do End_M(j,i) = new_vecs(EM_pix(i),j)
	lct_dmw ;in pcalib1.pro
	pix = EM_pix(i)
	wset,3
	XX = (pix mod XSize) * scl
	YY = (pix / XSize) * scl
	print,'Endmember X,Y,pix = ',XX/scl,YY/scl,pix
	print,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
	xyouts,XX,YY,'x',/dev,size=2,color=color,align=0.5, charthick=2
	wset,0
	!x = w0x
	!y = w0y
	oplot,[new_vecs(pix,Elem1)],[new_vecs(pix,Elem2)],$
		psym=6,color=color,thick=2
	wset,1
	!x = w1x
	!y = w1y
	oplot,[new_vecs(pix,Elem3)],[new_vecs(pix,Elem4)],$
		psym=6,color=color,thick=2
	wset,2
	!x = w2x
	!y = w2y
	oplot,[new_vecs(pix,Elem5)],[new_vecs(pix,Elem6)],$
		psym=6,color=color,thick=2
;
	Choice = 0
	xmenu, ['    CHANGE Box Center Point    ', 'CHANGE Box Size',$
		'SELECT End-member'], BASE=base, BUTTONS=B, $
		title='PCA Endmember Selection'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = Choice, event.id
	WIDGET_CONTROL, event.top, /DESTROY
;
	case Choice of
	;Changing center point for box.
	0: begin
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		redraw,i ;in pcalib1.pro
	end ;0
	;Changing the size of the box.
	1: begin
		print
		print, '=================================================='
		print, '              CHANGING BOX SIZE'
		print, '=================================================='
		print, 'Length =',XBox
		print, 'Width =',YBox
		print, '=================================================='
		print, 'INPUT the new length, width of the box.'
		print, '(Length, Width)? '
		read, XBox, YBox
		width = [XBox,YBox]
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		redraw,i ;in pcalib1.pro
	end ;1
	2: begin
		v_names(i)='Endmember'+strno(i)
		junk=''
		read,junk,prompt='Enter endmember name ['+v_names(i)+']:  '
		if (junk ne '') then v_names(i)=junk
		printf,5,'Endmember X,Y,pix = ',XX/scl,YY/scl,pix
		printf,5,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
		i = i + 1
		clean(*) = 1
		color=color+1
	end ;2
	endcase
;	endif else begin
;		oplot,ZC[0],ZC[1],psym=6,color=color,thick=2
;
;		Choice=0
;		xmenu, ['    CHANGE Box Center Point    ', $
;				'Synthesize End-member'], BASE=base, BUTTONS=B, $
;				title='No Image Points in Box!'
;		WIDGET_CONTROL,base,/REALIZE
;		event = WIDGET_EVENT(base)
;		WIDGET_CONTROL, GET_UVALUE = Choice, event.id
;		WIDGET_CONTROL, event.top, /DESTROY
;;
;		case Choice of
;		;Changing center point for box.
;		0: begin
;			clean(*)=1
;			refr_pca,/color ;in pcalib1.pro
;			redraw,i ;in pcalib1.pro
;		end ;0
;		;Synthesizing an endmemeber spectrum
;		1: begin
;			k=n_oem_pix+i
;			zvec=fltarr(WB)
;			zvec[Element1]=ZC[0]
;			zvec[Element2]=ZC[1]
;			print,'PCA vector: ',zvec
;			ans=' '
;			read,ans,prompt='Would you like to add values to any other dimensions? '
;			while ans eq 'y' do begin
;				print,'PCA vector: ',zvec
;				read,l,zl,prompt='Enter dimension number and value: '
;				zvec[l]=zl
;				read,ans,prompt='Again (y/n)? '
;			endwhile
;			new_vecz = zvec ## eigvecs
;			new_vecs[exclude[k],*]=new_vecz
;			endmz=cm_vector*(invert(eigvecs) ## transpose(new_vecz))
;			PCArr[exclude[k],*]=endmz
;			v_names(i)='Endmember'+strno(i)
;			junk=''
;			read,junk,prompt='Enter endmember name ['+v_names(i)+']:  '
;			if (junk ne '') then v_names(i)=junk
;			printf,5,'Endmember pix = ',pix
;			i = i + 1
;			clean(*) = 1
;			color=color+1
;		end ;1
;		endcase
;	endelse
endwhile
clean(*)=1
;
return
end

;==========================================================================
PRO s_pcakb_em2
;
;Procedure for interactively selecting end-member spectra from
;the keyboard in PCA space
;
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 PCA7, cm_vector
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, v_names, StoN
;
EM_pix = intarr(N_EndM)
EM_pix = long(EM_pix)
EM_pix(0) = -1
EM_pix_set=make_array(5,XSize*YSize,/long,value=(-1))
color=1
i=0
while (i lt N_EndM) do begin
	print,'Enter endmember 0th through 4th PC coordinates: '
	read,ZC0,ZC1,ZC2,ZC3,ZC4
	k=n_oem_pix+i
	zvec=fltarr(WB)
	zvec[0]=ZC0
	zvec[1]=ZC1
	zvec[2]=ZC2
	zvec[3]=ZC3
	zvec[4]=ZC4
	print,'PCA vector: ',zvec
	ans=' '
	read,ans,prompt='Add values to any other dimensions (y/n)? '
	while ans eq 'y' do begin
		print,'PCA vector: ',zvec
		read,l,zl,prompt='Enter dimension number and value: '
		zvec[l]=zl
		read,ans,prompt='Again (y/n)? '
	endwhile
	new_vecs[exclude[k],*]=zvec
;	endmz=(invert(eigvecs) ## transpose(zvec))
	endmz=cm_vector*(invert(eigvecs) ## transpose(zvec))
	c_pix1=exclude[k]
	EM_pix[i]=exclude[k]
	EM_pix_set[i,0:n_elements(c_pix1)-1]=c_pix1
	PCArr[exclude[k],*]=endmz
	for j = 0,N_EndM-2 do End_M(j,i) = new_vecs(EM_pix(i),j)
	lct_dmw ;in pcalib1.pro
	pix = EM_pix(i)
	wset,3
	XX = fix((pix mod XSize) * scl)
	YY = fix((pix / XSize) * scl)
	print,'Endmember X,Y,pix = ',XX/scl,YY/scl,pix
	print,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
	printf,5,'Endmember X, Y, pix = ',XX/scl,YY/scl,pix
	printf,5,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
	xyouts,XX,YY,'x',/dev,size=2,color=color,align=0.5, charthick=2
	wset,0
	!x = w0x
	!y = w0y
	oplot,[new_vecs(pix,Elem1)],[new_vecs(pix,Elem2)],$
		psym=6,color=color,thick=2
	wset,1
	!x = w1x
	!y = w1y
	oplot,[new_vecs(pix,Elem3)],[new_vecs(pix,Elem4)],$
		psym=6,color=color,thick=2
	wset,2
	!x = w2x
	!y = w2y
	oplot,[new_vecs(pix,Elem5)],[new_vecs(pix,Elem6)],$
		psym=6,color=color,thick=2
;
	v_names[i]='Endmember'+strno(i)
	junk=''
	read,junk,prompt='Enter endmember name ['+v_names[i]+']:  '
	if junk ne '' then v_names[i]=junk
	i=i+1
	color=color+1
endwhile
clean(*)=1
;
return
end

;==========================================================================
PRO s_im_em2
;
;Procedure for interactively selecting end-member spectra from the Image.
;This procedure uses a box to extract the data and then takes the median in
;order to find the resulting end-member.
;
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 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, v_names, StoN
;
EM_pix = intarr(N_EndM)
EM_pix = long(EM_pix)
dispwidth=intarr(N_EndM,2)
XBox = scl * 3
YBox = scl * 3
EM_pix_set=make_array(5,1024,/long,value=(-1))
Left = 1
;for j = 0,N_EndM-2 do begin
;  End_M(j,0) = Sky_pca(j)
;endfor
win = 3
color=1
i = 0
while (i lt N_EndM) do begin
	wset,win
	print
	print, '================================================='
	print, 'Click LEFT button on Image to choose center point'
	print, 'for box to select end-member.'
	print, '================================================='
	device, /cursor_crosshair
	cursor, X, Y, /dev
;
	if(!err eq Left) then begin
		lct_dmw ;in pcalib1.pro
		width = [XBox,YBox]
		tvbox, width, X, Y, color ;in pcalib3.pro
		clean(3) = 1
	endif
	width(0) = Xbox
	width(1) = Ybox
	dispwidth(i,0)=Xbox
	dispwidth(i,1)=Ybox
	X_im = X
	Y_im = Y
	image_pix ;in pcalib1.pro
	Numb_pix = n_elements(c_pix)
	med1 = MEDIAN(new_vecs(c_pix,Elem1))
	med2 = MEDIAN(new_vecs(c_pix,Elem2))
	res1 = fltarr(Numb_pix)
	for kk = 0,Numb_pix-1 do $
		res1(kk) = (new_vecs(c_pix(kk),Elem1)-med1)^2 + $
		(new_vecs(c_pix(kk),Elem2)-med2)^2
	tmp1 = MIN(res1,elm)
	EM_pix(i) = c_pix(elm)
	EM_pix_set(i,0:n_elements(c_pix)-1)=c_pix
	for j = 0,N_EndM-2 do End_M(j,i) = new_vecs(EM_pix(i),j)
	lct_dmw ;in pcalib1.pro
	pix = EM_pix(i)
	print,'Endmember X,Y,pix = ',(pix mod XSize),(pix / XSize),pix
	print,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
	wset,0
	!x = w0x
	!y = w0y
	oplot,[new_vecs(pix,Elem1)],[new_vecs(pix,Elem2)],psym=6, $
		color=color,thick=2
 	wset,1
 	!x = w1x
 	!y = w1y
 	oplot,[new_vecs(pix,Elem3)],[new_vecs(pix,Elem4)],psym=6, $
		color=color,thick=2
 	wset,2
 	!x = w2x
 	!y = w2y
 	oplot,[new_vecs(pix,Elem5)],[new_vecs(pix,Elem6)],psym=6, $
		color=color,thick=2
;
	Choice = 0
	xmenu, ['CHANGE BOX CENTER POINT',$
		'CHANGE BOX SIZE', 'SELECT ENDMEMBER'], BASE=base, BUTTONS=B, $
		title='Image Endmember Selection'
	WIDGET_CONTROL,base,/REALIZE
	event = WIDGET_EVENT(base)
	WIDGET_CONTROL, GET_UVALUE = Choice, event.id
	WIDGET_CONTROL, event.top, /DESTROY
;
	case Choice of
	;Changing center point for box.
	0: begin
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		redraw,i ;in pcalib1.pro
	end ;0
	;Changing the size of the box.
	1: begin
		print
		print, '=================================================='
		print, '              CHANGING BOX SIZE'
		print, '=================================================='
		print, 'Length =',XBox
		print, 'Width =',YBox
		print, '=================================================='
		print, 'INPUT the new length, width of the box.'
		print, '(Length, Width)? '
		read, XBox, YBox
		width = [XBox,YBox]
		clean(*)=1
		refr_pca,/color ;in pcalib1.pro
		redraw,i ;in pcalib1.pro
	end ;1
	2: begin
		v_names(i)='Endmember'+strno(i)
		junk=''
		read,junk,prompt='Enter endmember name ['+v_names(i)+']:  '
		if (junk ne '') then v_names(i)=junk
		printf,5,'Endmember X,Y,pix = ',(pix mod XSize),(pix / XSize),pix
		printf,5,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
		i = i + 1
		clean(*) = 1
		color=color+1
	end ;2
	else:
	endcase
endwhile
;
clean(*)=1
return
end

;==========================================================================
PRO s_kb_em2
;
;Procedure for interactively selecting end-member spectra from the Keyboard.
;
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 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, v_names, StoN
;
EM_pix = intarr(N_EndM)
EM_pix = long(EM_pix)
XBox = scl * 0
YBox = scl * 0
EM_pix_set=make_array(5,1024,/long,value=(-1))
i = 0
EM_pix(0) = -1
win = 3
color=1
while (i lt N_EndM) do begin
	print,'Enter endmember X,Y coordinate:'
	read,X,Y
	X=X*scl
	Y=Y*scl
;
	lct_dmw ;in pcalib1.pro
	wset,3
	width = [XBox,YBox]
	tvbox, width, X, Y, color, thick=2 ;in pcalib3.pro
	clean(3) = 1
	width(0) = Xbox
	width(1) = Ybox
	X_im = X
	Y_im = Y
	image_pix ;in pcalib1.pro
	Numb_pix = n_elements(c_pix)
	med1 = MEDIAN(new_vecs(c_pix,Elem1))
	med2 = MEDIAN(new_vecs(c_pix,Elem2))
	res1 = fltarr(Numb_pix)
	for kk = 0,Numb_pix-1 do $
		res1(kk) = (new_vecs(c_pix(kk),Elem1)-med1)^2 + $
		(new_vecs(c_pix(kk),Elem2)-med2)^2
	tmp1 = MIN(res1,elm)
	EM_pix(i) = c_pix(elm)
	EM_pix_set(i,0:n_elements(c_pix)-1)=c_pix
	for j = 0,N_EndM-2 do End_M(j,i) = new_vecs(EM_pix(i),j)
	lct_dmw ;in pcalib1.pro
	pix = EM_pix(i)
	print,'Endmember X,Y,pix = ',X_im,Y_im,pix
	print,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
	wset,0
	!x = w0x
	!y = w0y
	oplot,[new_vecs(pix,Elem1)],[new_vecs(pix,Elem2)],psym=6, $
		color=color,thick=2
 	wset,1
 	!x = w1x
 	!y = w1y
 	oplot,[new_vecs(pix,Elem3)],[new_vecs(pix,Elem4)],psym=6, $
		color=color,thick=2
 	wset,2
 	!x = w2x
 	!y = w2y
 	oplot,[new_vecs(pix,Elem5)],[new_vecs(pix,Elem6)],psym=6, $
		color=color,thick=2
;
;	Choice = 0
;	xmenu, ['CHANGE BOX CENTER POINT',$
;		'SELECT ENDMEMBER'], BASE=base, BUTTONS=B, $
;		title='Keyboard Endmember Selection'
;	WIDGET_CONTROL,base,/REALIZE
;	event = WIDGET_EVENT(base)
;	WIDGET_CONTROL, GET_UVALUE = Choice, event.id
;	WIDGET_CONTROL, event.top, /DESTROY
;
;	case Choice of
;	;Changing center point for box.
;	0: begin
;		wset,win
;		tvscl, Image < High(Im_Num) > Low(Im_Num)
;		clean(3) = 0
;	end ;0
;	1: begin
		v_names(i)='Endmember'+strno(i)
		junk=''
		read,junk,prompt='Enter endmember name ['+v_names(i)+']:  '
		if (junk ne '') then v_names(i)=junk
		printf,5,'Endmember X,Y,pix = ',X_im,Y_im,pix
		printf,5,'PCA coordinates = ',transpose(new_vecs[EM_pix(i),*])
		i = i + 1
		clean(*)=1
		color=color+1
;	end ;1
;	endcase
endwhile
;
clean(*)=1
return
end

;==========================================================================
PRO frac_ab
;
;This procedure calculates the fractional abundances of endmembers contained
;in each pixel of the image.   The algorithm used to determine the abundances
;comes from Salmon, 'Treatise on the Analytic Geometry of Three Dimensions',
;p. 19.
;
;
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,wave,a,wavtitle
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, v_names, StoN
;
res1 = 1
N_pix = n_elements(res_pix)
if (N_pix eq XSize*YSize) then begin
	res1 = 0
endif
if (N_EndM le 0) then s_endm
N_dim = N_EndM - 1
Mix_EM = fltarr(N_dim,2*N_EndM+1)
Frac1 = fltarr(N_pix,N_EndM+1)
Frac = fltarr(N_pix,N_EndM)
Frac2 = fltarr(N_pix,N_EndM)
FracTot = fltarr(N_pix)
FracT = fltarr(N_pix)
FracTot(*) = 0.0
Work = fltarr(N_EndM,N_EndM)
;
;Calculating abundances over all pixels, BEGINNING OF BIG LOOP over pix
for i = 0l,N_pix-1 do begin
;
	if  (res1 eq 0) then begin
		;Form Mix_EM from mixtures and end_members.
		Mix_EM(*,1:N_EndM) = End_M(0:N_dim-1,0:N_EndM-1)
		Mix_EM(*,1+N_EndM:N_EndM+N_EndM) = End_M(0:N_dim-1,0:N_EndM-1)
		Mix_EM(*,0)=new_vecs(i,0:N_dim-1)
	endif else begin
		;Form Mix_EM from mixtures and end_members for residuals.
		Mix_EM(*,1:N_EndM) = End_M(0:N_dim-1,0:N_EndM-1)
		Mix_EM(*,1+N_EndM:N_EndM+N_EndM) = End_M(0:N_dim-1,0:N_EndM-1)
		Mix_EM(*,0)=new_vecs(res_pix(i),0:N_dim-1)
	endelse
;
;Create work array and find determinant.
	Work(N_dim,*) = 1.0
	Work(0:N_dim-1,0) = Mix_EM(*,0)
;
	for ii = 0,N_dim do begin
		Work(0:N_dim-1,1:N_dim) = Mix_EM(*,1+ii:N_dim+ii)
		Frac1(i,ii) = DETERM(Work)
	endfor ;ii
	FracTot(i) = FracTot(i) + Total(Frac1(i,0:N_dim))  ;modified
;
endfor ;i
;END OF BIG LOOP
;
for kk = 0,N_Dim do Frac1(*,kk) = Frac1(*,kk) / FracTot(*)
Frac1(*,N_EndM) = Frac1(*,0)    ;new
Frac = Frac1(*,1:N_EndM)   ;new
;
;New method of calculating fractions for all but sky.
;for kk = 0l,N_pix-1 do begin
 ; for mm = 0,N_EndM-1 do begin
  ;  FracT(kk) = FracT(kk) + Frac1(kk,mm)
 ; endfor
;endfor
;
;for n = 0l,N_pix-1 do begin
 ; for o = 0,N_EndM-1 do begin
  ;  Frac2(n,o) = Frac(n,o) / (FracT(n) - Frac(n,0))
 ; endfor
;endfor
;Frac2(*,0) = Frac(*,0)
;Frac2 = Frac
;
return
end

;==========================================================================
PRO d_frac
;
;Procedure for displaying fractional abundance information.
;
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, v_names, StoN
;
print
print, 'Displaying fractional abundance images.'
print
;
;titl = strarr(5)
;titl = ['Component 1','Component 2','Component 3','Component 4','Component 5']
;
Frac1 = fltarr(XSize*YSize,N_EndM)
N_pix = n_elements(res_pix)
for j = 0l,N_pix-1 do begin
	for i = 0,N_EndM-1 do begin
		Frac1(res_pix(j),i) = Frac(j,i)
	endfor ;i
endfor ;j
loadct,0
for i = 0,N_EndM-1 do begin
	window,i+4,xs=scl*XSize,ys=scl*YSize,title=v_names(i),retain=2
	b = Frac1(*,i)
	d = reform(b,XSize,YSize)
	c = rebin(d,XSize*scl,YSize*scl,/sa)
	tvscl, c < 1.0 > 0.0
endfor ;i
;
return
end

;==========================================================================
PRO d_resid
;
;Procedure for displaying the residual information.
;
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, v_names, StoN
COMMON RESIM, Res_Spec
COMMON RESCALE, xygrid, Dim1, Dim2, xmargin, ymargin, bmargin, tmargin
;
N_pix = n_elements(res_pix)
sigma = fltarr(N_pix)
sigma1 = fltarr(N_pix)
sigma2 = fltarr(N_pix)
sigma3 = fltarr(N_pix)
Th_Im = fltarr(N_pix,NImages)
Data2 = fltarr(N_pix,NImages)
Fin_Ims = fltarr(N_pix,NImages)
Fin_Im = fltarr(N_pix)
EM_flux = fltarr(N_EndM,NImages)
Noise1 = fltarr(NImages)
Noise2 = fltarr(NImages)
Res_Spec=fltarr(XSize*YSize,NImages)
;
print, 'Enter three values for sigma multiplicative factors'
read, mult1, mult2, mult3
printf,5,'sigma',mult1,mult2,mult3
;
mult1 = mult1*mult1
mult2 = mult2*mult2
mult3 = mult3*mult3
;
for i = 0l,N_pix-1 do Data2(i,*) = PCArr(res_pix(i),*)
;
Noisestimate,Data2,StoN,sigma ;in pcalib2.pro
sigma1 = sigma * mult1
sigma2 = sigma * mult2
sigma3 = sigma * mult3
;
;Calculating the flux from end_member pixels.
for i = 0,N_EndM-1 do EM_flux(i,*) = PCArr(EM_pix(i),*)
;
;Constructing theoretical image from fractional abundances and end-members.
   Th_Im =  Frac # EM_flux   ;new
;
;Subtracting theoretical images from original images.
   Res_Buffer = (Data2 - Th_Im)
   Fin_Ims = Res_Buffer * Res_Buffer
   for i=0l,N_pix-1 do Res_Spec(res_pix(i),*)=Res_Buffer(i,*)
;
;Find final residual image and noise associated with each filter.
for i = 0 ,NImages-1 do Noise1(i) = Total(Fin_Ims(*,i))
Noise2 = sqrt(Noise1)
for j = 0l  ,N_pix-1 do Fin_Im(j) = Total(Fin_Ims(j,*))
;
;Selecting pixels for display purposes by number of sigma from rms value.
col1 = where(Fin_Im  lt sigma1)
col2 = where(Fin_Im  gt sigma1 and Fin_Im  lt sigma2)
col3 = where(Fin_Im  gt sigma2 and Fin_Im  lt sigma3)
col4 = where(Fin_Im  gt sigma3)
col1 = col1 -2 *(col1 < 0)
col2 = col2 -2 *(col2 < 0)
col3 = col3 -2 *(col3 < 0)
col4 = col4 -2 *(col4 < 0)
print, 'Number of  BLACK pixels is ',n_elements(col1)
print, 'Number of    RED pixels is ',n_elements(col2)
print, 'Number of YELLOW pixels is ',n_elements(col3)
print, 'Number of  GREEN pixels is ',n_elements(col4)
printf,5, 'Number of  BLACK pixels is ',n_elements(col1)
printf,5, 'Number of    RED pixels is',n_elements(col2)
printf,5, 'Number of YELLOW pixels is',n_elements(col3)
printf,5, 'Number of  GREEN pixels is',n_elements(col4)
;
b = Data1(Im_Num,XLow:XHigh,YLow:YHigh)
c = reform(b,XSize*YSize)
lct_dmw ;in pcalib1.pro
;
;Displaying residual data in window 8
c_pix2 = bytscl(c) > 8b
c_pix2(res_pix(col1)) = 0b
c_pix2(res_pix(col2)) = 1b
c_pix2(res_pix(col3)) = 2b
c_pix2(res_pix(col4)) = 3b
c_pix3 = reform(c_pix2,XSize,YSize)
c_pix4 = rebin(c_pix3,XSize*scl,YSize*scl,/sa)
tm = where(c_pix4 gt 3)
if tm(0) ne -1 then c_pix4(tm) = 0
window,12,xs=scl*XSize,ys=scl*YSize,title='Residual Image',retain=2
tv, c_pix4
;
;Writing residual image to GIF/FITS files.
resname = strarr(10)
resname = ['resp0','resp1','resp2','resp3','resp4',$
           'resp5','resp6','resp7','resp8','resp9']
;write_gif,resname(N_iter)+'.gif',bytscl(c_pix4)
lct_dmw
tvlct,rvect,gvect,bvect,/get
write_png,resname(N_iter)+'.png',bytscl(c_pix4),rvect,gvect,bvect
;gifbar,c_pix4,0,4,fname=resname(N_iter)+'.disp.gif',ct=0,dispwin=20
;gifbar,c_pix4,0,4,fname=resname(N_iter)+'.disp.c.gif',ct=39,dispwin=20
gifbar,c_pix4,0,4,fname=resname(N_iter)+'.disp',ct=0,dispwin=20
gifbar,c_pix4,0,4,fname=resname(N_iter)+'.disp',ct=39,dispwin=20
BufImage=rebin(c_pix3,Xsize*xygrid,YSize*xygrid)
tm = where(BufImage gt 3)
if tm(0) ne -1 then BufImage(tm) = 0
OutImage=fltarr(XSize,YSize)
OutImage(xmargin:XSize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
writefits,resname(N_iter)+'.fit',OutImage
;
;Displaying residual data on the PCAPlot.
loadct,0
window,13,xs=400,ys=400,title='PCAPlot of Residuals',retain=2
;window,13,xs=600,ys=400,title='PCAPlot of Residuals',retain=2
!x = w0x
!y = w0y
lct_dmw ;in pcalib1.pro
pcaplot,Elem1,Elem2,ps=1,iter=N_iter,outfile=13, $
	xrange=xrange(*,0),yrange=yrange(*,0) ;in pcalib3.pro
!x = w0x
!y = w0y
X_pca1 = new_vecs(res_pix(col2),Elem1)
Y_pca1 = new_vecs(res_pix(col2),Elem2)
if n_elements(X_pca1) gt 1 then $
	oplot,X_pca1,Y_pca1,psym=3,color=1 $
else print, 'No pixels in sigma range: ',sqrt(mult1),' to ',sqrt(mult2)
X_pca2 = new_vecs(res_pix(col3),Elem1)
Y_pca2 = new_vecs(res_pix(col3),Elem2)
if n_elements(X_pca2) gt 1 then $
	oplot,X_pca2,Y_pca2,psym=3,color=2 $
else print, 'No pixels in sigma range: ',sqrt(mult2),' to ',sqrt(mult3)
X_pca3 = new_vecs(res_pix(col4),Elem1)
Y_pca3 = new_vecs(res_pix(col4),Elem2)
if n_elements(X_pca3) gt 1 then $
	oplot,X_pca3,Y_pca3,psym=3,color=3 $
else print, 'No pixels in sigma range: greater than ',sqrt(mult3)
print,'Total Error per Filter',Noise2
printf,5,'Total Error per Filter',Noise2
;
;Writing frac maps of modeled pixels to GIF/FITS files
endfn = strarr(5,11)
endfn(0,*) = ['end00','end10','end20','end30','end40','end50',$
         'end60','end70','end80','end90','end100']
endfn(1,*) = ['end01','end11','end21','end31','end41','end51',$
         'end61','end71','end81','end91','end101']
endfn(2,*) = ['end02','end12','end22','end32','end42','end52',$
         'end62','end72','end82','end92','end102']
endfn(3,*) = ['end03','end13','end23','end33','end43','end53',$
         'end63','end73','end83','end93','end103']
endfn(4,*) = ['end04','end14','end24','end34','end44','end54',$
         'end64','end74','end84','end94','end104']
;
Frac1 = fltarr(XSize*YSize,N_EndM)
for j = 0l,N_pix-1 do begin
	for i = 0,N_EndM-1 do begin
		Frac1(res_pix(j),i)=Frac(j,i)
	endfor ;i
endfor ;j
for ii=0,N_EndM-1 do begin
	FracD = fltarr(XSize,YSize)
	FracD = reform(Frac1(*,ii),XSize,YSize)
	FracD(res_pix(col2)) = 0
	FracD(res_pix(col3)) = 0
	FracD(res_pix(col4)) = 0
	printf,5,'>1.2',n_elements(where(FracD gt 1.2))
	printf,5,'<-0.1',n_elements(where(FracD lt -0.1))
	FracDgif = rebin(FracD,XSize*scl,YSize*scl,/sa)
;	gifbar,FracDgif,0,1,fname=endfn(ii,N_iter)+'.gif',ct=0,dispwin=20
;	gifbar,FracDgif,0,1,fname=endfn(ii,N_iter)+'.c.gif',dispwin=20
	gifbar,FracDgif,0,1,fname=endfn(ii,N_iter),ct=0,dispwin=20
	gifbar,FracDgif,0,1,fname=endfn(ii,N_iter),dispwin=20
	BufImage=rebin(FracD,Xsize*xygrid,YSize*xygrid)
	OutImage=fltarr(XSize,YSize)
	OutImage(xmargin:XSize-1-ymargin,bmargin:YSize-1-tmargin)=BufImage
	mkhdr,hdr,OutImage
	vertx=EM_pix(ii) mod Xsize
	vertx=vertx*xygrid+xmargin
	comment='x-coord of endmember pixel'
	sxaddpar,hdr,'x_coord',vertx,comment
	verty=EM_pix(ii)/Xsize
	verty=verty*xygrid+bmargin
	comment='y-coord of endmember pixel'
	sxaddpar,hdr,'y_coord',verty,comment
	comment='Endmember name'
	sxaddpar,hdr,'vname',v_names(ii),comment
	writefits,endfn(ii,N_iter)+'.fit',OutImage,hdr
endfor ;ii
wdelete,20
;
n_c2 = n_elements(col2)
n_c3 = n_elements(col3)
n_c4 = n_elements(col4)
N_pix2 = n_c2 + n_c3 + n_c4
res_pix2 = intarr(N_pix2)
res_pix2 = long(res_pix2)
for i = 0l,n_c2-1 do res_pix2(i) = res_pix(col2(i))
for i = 0l,n_c3-1 do res_pix2(i+n_c2) = res_pix(col3(i))
for i = 0l,n_c4-1 do res_pix2(i+n_c2+n_c3) = res_pix(col4(i))
res_pix = intarr(N_pix2)
res_pix = long(res_pix)
for i = 0l,N_pix2-1 do res_pix(i)=res_pix2(i)
;
for i = 0,N_EndM-1 do old_em_pix(i+n_oem_pix) = EM_pix(i)
n_oem_pix=n_oem_pix+N_EndM
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
N_EM(0) = N_EndM-1
;for j = 0,NImages-1 do printf,4,Sky(j)
;
printf,4,N_EndM,NImages
for i = 0,N_EndM-1 do begin
	for j = 0,NImages-1 do begin
		printf,4,j,PCArr(EM_pix(i),j)
	endfor ;j
endfor ;i
;
return
end

;==========================================================================
PRO Noisestimate,Image,StoN,Noisarray
; This array creates a noise array from the data.  The output from
; this routine is the sum (across wavelength) of the square errors
; so to get RMS error it must be divided by number of wavelengths
; and the square-root taken.  Both of these must be done outside of
; this routine.
temparray=float(Image)
temparray=temparray/StoN
temparray=temparray*temparray
Noisarray = total(temparray,2)
;
return
end