pro selplot,nom,cube,whichplot,lvlflag=lvlflag,nplots=nplots
;------------------------------------------
;Selects greyscale and plots contour maps for a data cube.
;The inputs to the routine include:
;	CUBE - the data cube to be analyzed
;	NOM - The name of the output postscript file
;	BOXP - A two element vector: [resolution_size,pixel_size]
;	WHICHPLOT - A vector containing the indexes of the plots 
;		to be plotted
;	NPLOTS - The number of plots on the page in each direction.
;------------------------------------------
if n_elements(lvlflag) eq 0 then lvlflag=0
if n_elements(nplots) eq 0 then nplots=1
boxp=[0.,0.]
set_plot,'ps'
device,file=nom,/inches,xoffset=0.5,yoffset=0.5,xsize=7.5,ysize=10,$
/color,bits_per_pixel=4,/portrait
erase

numberplot=n_elements(whichplot)
a=replicate(255,256)-findgen(256)
b=findgen(256)

sz=size(cube)
sx=float(sz(1))
sy=float(sz(2))
spam=16
;float(sz(3))

pasp=11/8.5
asp=(sy/sx)
fillfac=0.8
xpl=float(nplots)
ypl=float(nplots)
scale=1./(nplots)
dpx=1./xpl
dpy=1./ypl
x0=0.0
y0=1-dpy
arat=asp/pasp

plotcounter=0

if (asp le pasp) then begin
s=1
if asp gt 1. then s=1/asp 
pos=[(1-s)/2,1-s*arat,1-(1-s)/2,1]
posbar=[0.1,0.1*(1-arat),0.9,0.5*(1-arat)]
endif else begin
pos=[0,0,1/arat,1]
xend=1-1/arat
posbar=[1-0.7*xend,0.1,1-0.3*xend,0.9]
endelse
pos0=fillfac*pos+replicate((1-fillfac)/2,4)
posbar0=fillfac*posbar+replicate((1-fillfac)/2,4)
ttles=strarr(16)
ttles=['!8S!Iij!3','!8l!8!Iij!3','!8s!Iij!3','!8S!Iij!El!3',$
'!8l!8!Iij!El!3','!8S!Iij!Es!N!3','!8s!Iij!Es!N!3','!8S!Iij!E0!N!3',$
'!8T!IA!E*!N!3','!4D!8v!3','!8v!ILSR!3','!8T!Iint!3','!4r!8!Irms!3',$
'!8T!E*!IA!N/!4r!8!Irms!3','!10"!8l!8!Iij!N!10"!3',$
'!10"!8l!8!El!Iij!N!10"!3','!10"!8l!8!Iij!N!10"!3',$
'!10"!8l!8!El!Iij!N!10"!3']

newpage:

for j=0,ypl-1 do begin
;for i=0,xpl-1 do begin

;px0=x0+i*dpx
;px1=px0+dpx

py0=y0-j*dpy
py1=py0+dpy

for i=0,xpl-1 do begin
q=whichplot(plotcounter)
cubei=cube(*,*,q)
pointer=where(cubei ne 99999.0)
if pointer(0) eq -1 then begin
	xyouts,pos(0),pos(1),'No Valid Data',charsize=2
	goto,skipout
endif
mival=min((cubei ne 99999.0)*cube)
maval=max((cubei ne 99999.0)*cube)
mini=mival*(0.9*(mival gt 0.)+1.1*(mival lt 0.))
maxi=maval*(1.1*(maval gt 0.)+0.9*(maval lt 0.))
px0=x0+i*dpx
px1=px0+dpx
;py0=y0-j*dpy
;py1=py0+dpy
offst=[px0,py0,px0,py0]
pos=scale*pos0+offst
posbar=scale*posbar0+offst

;ylen=py1-py0-2*dply0
;pos=[px0+dplx0,py0+dply0+offy,px1-dplx0,py1-dply0+offy]
;posbar=pos-[0.,(dply0+offy)*0.7,0.,(dply0+offy)*0.5+ylen]
;print,pos
cub=cube(*,*,q)
lvl=levs(cub,15)
lvl=lvl(uniq(lvl))
pointer=where(cub ne 99999.0)
dit=cub(pointer)
mini=min(cub)
maxi=max(cub)
if lvlflag ne 0 then begin
	lvl=findgen(15)*(maxi-mini)/15.+mini
	print,'Difference between contours:',(maxi-mini)/15.
endif
;print,q,mini,maxi
colmatrix=((max(lvl)-min(lvl))/100*findgen(100)+min(lvl))#[1.,1.,1.,1.]
if arat ge 1 then colmatrix=transpose(colmatrix)

tvlct,a,a,a
;loadct,6
contour,cube(*,*,q),$
position=pos,$
levels=lvl,/cell_fill,max_value=99998.0,xrange=[0,sx-1],yrange=[0,sy-1],$
/noerase,xstyle=5,ystyle=5

contour,colmatrix,/fill,$
position=posbar,levels=lvl,/noerase,xstyle=5,ystyle=5


tvlct,b,b,b

if q eq 1 or q eq 4 or q eq 10 then begin
ts=(1-total(lvl ge 1)/n_elements(lvl))*255
print,ts
if ts lt 135 and ts gt 120 then ts=255
abba=(max(lvl)-min(lvl))/25
contour,cube(*,*,q),position=pos,levels=[0.],color=ts,xstyle=5,$
ystyle=5,/noerase,max_value=99998.0,xrange=[0,sx-1],yrange=[0,sy-1]
contour,colmatrix,position=posbar,levels=[0.],xstyle=5,c_color=ts,$
ystyle=5,/noerase
endif
if q eq 2 or q eq 6 then begin
ts=(1-total(lvl ge 1)/n_elements(lvl))*255
if ts lt 140 and ts gt 115 then ts=255
print,ts
contour,cube(*,*,q),position=pos,levels=[1.],color=ts,xstyle=5,$
ystyle=5,/noerase,max_value=99998.0,xrange=[0,sx-1],yrange=[0,sy-1]
contour,colmatrix,position=posbar,levels=[1.],xstyle=5,c_color=ts,$
ystyle=5,/noerase
;stop
endif

contour,cube(*,*,q),position=pos,/nodata,$
title=ttles(q),xtitle='!3x (pixels)',ytitle='!3y (pixels)',$
xticks=1,yticks=1,xcharsize=0.8,ycharsize=0.8,$
xtickname=[' ',' '],ytickname=[' ',' '],/noerase,charsize=1.

xtn=strarr(3)
xtn(0)=sigfig(min(lvl),3)
xtn(1)=sigfig((max(lvl)-min(lvl))/2+min(lvl),3)
xtn(2)=sigfig(max(lvl),3)
if total(boxp) ne 0 then begin
res=boxp(0)
scl=boxp(1)

bdx0=scl*(pos(2)-pos(0))/sx
bdy0=scl*(pos(3)-pos(1))/sy
bx0=pos(0)+bdx0
by0=pos(1)+bdy0
by1=by0+(res-1)/2*bdy0
bx1=bx0+(res-1)/2*bdx0
plots,bx0,by0,/normal
plots,bx0+res*bdx0,by0,/continue,/normal
plots,bx0+res*bdx0,by0+res*bdy0,/continue,/normal
plots,bx0,by0+res*bdy0,/continue,/normal
plots,bx0,by0,/continue,/normal
plots,bx1,by1,/normal
plots,bx1+bdx0,by1,/normal,/continue
plots,bx1+bdx0,by1+bdy0,/normal,/continue
plots,bx1,by1+bdy0,/normal,/continue
plots,bx1,by1,/normal,/continue
endif


if asp le pasp then begin
;contour,colmatrix,/nodata,position=posbar,$
;xstyle=9,ystyle=1,xticks=2,yticks=1,xticklen=-0.25,$
;xtickname=[mini,(maxi-mini)/2+mini,maxi],$
;ytickname=[' ',' '],/noerase,xcharsize=0.7

contour,colmatrix,/nodata,xstyle=1,ystyle=1,/noerase,$
position=posbar,$
xticks=1,yticks=1,xtickname=[' ',' '],ytickname=[' ',' ']
axis,0,0,xax=0,/data,xticks=2,$
xtickname=xtn,xticklen=-0.25,xcharsize=0.75
endif else begin

contour,colmatrix,/nodata,xstyle=1,ystyle=1,/noerase,$
position=posbar,$
xticks=1,yticks=1,xtickname=[' ',' '],ytickname=[' ',' ']
axis,3,0,yax=3,/ynozero,/data,yticks=2,$
ytickname=xtn,yticklen=-0.15,ycharsize=0.75
endelse

skipout:
;q=q+1
;if q eq spam then goto,ed
plotcounter=plotcounter+1
if plotcounter eq numberplot then goto,ed
endfor
endfor
if plotcounter lt numberplot then begin
erase
goto,newpage
endif


ed:
device,/close
set_plot,'x'
end