; plotspatial.pro
; IDL program written by Tim Mitchell on 17.05.01
; last modified on 31.07.03
; permits b/w line coutours, or colour contours, or grid cells

; ******************************************************************************
; unified entry

pro PlotSpatial

print, ""
print, "  > ***** PlotSpatial.pro: plots .glo(s) *****"
print, ""

Black = 0 & Grey = 254 & White = 255
MissVal  = -999.0 & MissCol = White
ScaleSegN = 10 & ModelChosen = -999 & BatchFileN = 0
InputText = "" & MainTitle=""

MonthNames = strarr (12)
MonthNames [*] = ['January','February','March','April','May','June','July','August', $
		  'September','October','November','December']
SeasonNames = strarr (4)
SeasonNames = ['winter (DJF)','spring (MAM)','summer (JJA)','autumn (SON)']
ColorBounds  = [1,253]
; ColorBounds = [128,253]			; ######################################

print, "  > Use standard (=0) or custom (=1) specifications ? "
read, QStanCust

; ******************************************************************************
; get basic stuff

print, "  > Enter the main title of the plot: "
read, MainTitle

print, "  > Enter the number of .glo to plot: "
read, PlotN

if            (PlotN EQ 1) then begin
  ColN = 1 & RowN = 1
endif else if (PlotN EQ 4) then begin
  ColN = 2 & RowN = 2
endif else if (PlotN EQ 6) then begin
  ColN = 3 & RowN = 2
endif else if (PlotN EQ 9) then begin
  ColN = 3 & RowN = 3
endif else if (PlotN EQ 12) then begin
  ColN = 3 & RowN = 4
endif else begin
  print, "  > Specify the number of columns: "
  read, ColN  
  print, "  > Specify the number of rows: "
  read, RowN
endelse

if (PlotN EQ 4) then begin
  print, "  > Do the 4 plots correspond to the four seasons ? (0=no,1=yes)"
  read, QSeasons
endif else if (PlotN EQ 12) then begin
  print, "  > Do the 12 plots correspond to the 12 months ? (0=no,2=yes)"
  read, QSeasons
endif else begin
  QSeasons = 0
endelse

; ******************************************************************************
; get model stuff

Models = intarr (PlotN) & ModelsLongN = intarr (PlotN) & ModelsLatN = intarr (PlotN)

if (PlotN EQ 1) then begin
  SelectModel, ModelChosen, ModelFilePath, LongN, LatN, /Extend
  Models(0)=ModelChosen & ModelsLongN(0)=LongN & ModelsLatN(0)=LatN
  QOneGrid = 1
endif else begin
  print, "  > Are all the .glo files on the same model grid ? (1=yes,2=no)"
  read, QOneGrid
  
  if (QOneGrid EQ 1) then begin
    SelectModel, ModelChosen, ModelFilePath, LongN, LatN, /Extend
    Models(*)=ModelChosen & ModelsLongN(*)=LongN & ModelsLatN(*)=LatN
  endif
endelse

GetViewBounds, ViewN, ViewName, ViewKeyWords, ViewBounds, ViewFullName=ViewFullName
GetScales, ScaleN, ScaleName, ScaleColTab, ScaleSeg, ScaleLimits 

FilePaths = strarr (PlotN)
FileInfos = strarr (PlotN)

if (QSeasons NE 1 AND PlotN GT 1 AND QOneGrid EQ 1) then begin
  print, "  > Load the .glo files as a batch (=1) or singly (=2) ?"
  read, QBatchSingle
endif else begin
  QBatchSingle = 2
endelse

if (QBatchSingle EQ 1) then begin
 while (BatchFileN NE PlotN) do begin
   GetFileBatch, Files, BatchFileN
   if (BatchFileN NE PlotN) then print, "   > Try again. This filter gives a file total of: ", BatchFileN
 endwhile
 
 for XPlot = 0, (PlotN-1) do begin
  CheckGlo, ModelChosen, FilePath, FileInfo, /NoPrompt, CallFile=Files(XPlot)
  FilePaths(XPlot) = FilePath
  FileInfos(XPlot) = FileInfo
 endfor
endif else begin
 for XPlot = 0, (PlotN-1) do begin
  if (QSeasons EQ 1) then begin
    print, "  > Enter .glo for season: ", SeasonNames [XPlot]
  endif else if (QSeasons EQ 2) then begin
    print, "  > Enter .glo for month: ", MonthNames [XPlot]
  endif else if (PlotN EQ 1) then begin
    print, "  > Enter the .glo file."
  endif else if (QOneGrid NE 1) then begin
    print, "  > Plot: ", XPlot
    ModelChosen = MissVal
    SelectModel, ModelChosen, ModelFilePath, LongN, LatN, /Extend
    Models(XPlot)=ModelChosen & ModelsLongN(XPlot)=LongN & ModelsLatN(XPlot)=LatN
    print, "  > Enter the .glo file: ", XPlot
  endif else if (XPlot EQ 0) then begin
    print, "  > Enter each .glo file in turn."
  endif

  CheckGlo, ModelChosen, FilePath, FileInfo, /NoPrompt
  FilePaths(XPlot) = FilePath
  FileInfos(XPlot) = FileInfo
 endfor
endelse

; ******************************************************************************
; plot choices

print, "  > Do we label the plots ? (0=no,1=yes)" 
read, QPlotName
 
if (QPlotName EQ 1 AND QSeasons EQ 0) then begin
   PlotName = strarr (PlotN)
   print, "  > Enter the title for each plot in turn: "
   for XPlot = 0, (PlotN-1) do begin
      read, InputText
      PlotName[XPlot] = InputText
      GiveSymbol, PlotName[XPlot]
   endfor
endif

QView = MissVal & InputInt = 0
while (QView LT 0) do begin
     print, "  > Identify the region to plot (-1=list): "
     read, InputInt
    
     if (InputInt EQ -1) then begin
      for XView = 0, (ViewN-1) do begin
        print, "  > ", XView, " : ", ViewName[XView], format='(a4,i4,a3,a10)'
      endfor
     endif else begin
      if (InputInt GE 0 AND InputInt LT ViewN) then QView = InputInt
     endelse
endwhile
Limits = [ViewBounds[QView,2],ViewBounds[QView,0],ViewBounds[QView,3],ViewBounds[QView,1]]
       
print, "  > Plot onto map (=0) or grid (=1) ? "
read, QMapGrid

print, "  > Plot contours (=0), cells (=1), or b/w lines (=2) ?"
read, QContourCell

print, "  > Select the contour density (1=min,5=max): "
read, ContourFactor

if (QMapGrid EQ 0) then begin
  print, "  > Plot grid ? (0=no,1=yes)"
  read, QGrid

  print, "  > Plot outlines ? (-2=IDLcoasts,-1=select,0=no,>0=predefined)"
  read, QOutline
  if (QOutLine EQ -1 OR QOutLine GT 0) then SelectPoly, QOutline
endif else begin
  QGrid = -999 & QOutline = -999
endelse

if (QcontourCell NE 2) then begin
  print, "  > Side bar: 0=none,1=all,2=attrib,3=info (scale IS plotted)"
endif else begin
  print, "  > Side bar: 0=none,1=all,2=attrib,3=info (scale NOT plotted)"
endelse
read, QBarInfo
QBottomRight = QBarInfo
if (QBarInfo eq 1) then QBottomRight = 3
if (QBarInfo eq 2) then QBottomRight = 0

if (QcontourCell EQ 2) then begin
  QSameScale = 0
endif else if (PlotN NE 1) then begin
  print, "  > Do the plots have the same scale ? (1=yes,2=no)"
  read, QSameScale
endif else begin
  QSameScale = 1
endelse

if (QSameScale EQ 1) then begin
  Info = strarr (10) & Info [*] = "missing"
  GrabInfo, "@@@@@@@@", "@@@@@@@@", Info, -999, 0
  
  if (QSeasons EQ 0) then begin  
    SeasName = "missing"
    SearchSeas, FilePaths(0), FileInfos(0), SeasName, SeasIndex=SeasIndex
  endif else begin
    SeasName = "annual" & SeasIndex=17
  endelse
endif

print, "  > Plot shape: standard (=0) or fit to .glo (=1) ?"
read, QStanFit

if (QStanFit EQ 1 AND QOneGrid EQ 2) then begin
  print, "  > Plot shape: enter LongN,LatN: "
  read, LongN,LatN
endif

PlotSizeFactor = 1.0
print, "  > Plot  size factor (standard=1.0)  ?"
read, PlotSizeFactor

LabelSizeFactor = 1.0
print, "  > Label size factor (standard=1.0)  ?"
read, LabelSizeFactor

TextSize = 0.8 & TextThick = 2.0

print, "  > Plot to screen (=21) or .eps (=22) ?"
read, Choice

; ******************************************************************************
; determine plot size and identify plot segments
; note that when multiple grids are used in a single execution, this will go wonky

TitlePixels = 2.0

if (QBarInfo EQ 0) then begin 				; this allows width for just the scale
  BarPixels = 3.0
endif else begin					; this allows width for the logo and info as well
  BarPixels = 4.0
endelse

if (QPlotName EQ 0) then begin 				; this allows height for the label, if required
  NamePixels = 0.0
endif else begin					
  NamePixels = 1.0
endelse

if (QStanFit EQ 0) then begin				; this makes any PS plot have a size 16cm * 12cm
							;   so fits into Word doc without any resizing
  XPixels = (16.0*PlotSizeFactor) & YPixels = (12.0*PlotSizeFactor)  	
endif else begin					; this sizes the plot according to the .glo sizes
					  		;   so minimises white space in the plot
  YPixels = (float(RowN)*PlotSizeFactor * ((float(LatN)/20.0) + (LabelSizeFactor*NamePixels))) $
  		+ (TitlePixels*LabelSizeFactor)
  
  if (QSameScale EQ 2) then begin
    XPixels =  float(ColN) * ( (float(LongN)*PlotSizeFactor/20.0) + BarPixels )
  endif else begin
    XPixels = (float(ColN) * float(LongN) * PlotSizeFactor / 20.0)  + BarPixels
  endelse
endelse

BarWidth   = BarPixels / XPixels			; this converts the key width into normalised units
NameHeight = LabelSizeFactor * NamePixels / YPixels	; ... and the name height likewise
TopBarHeight= TitlePixels / XPixels			; ... and the title bar

							; this identifies the normalised key boundaries
if (QSameScale EQ 2) then begin				; here, one scale per .glo
  GetCorners, [0.0,0.0,1.0,(1.0-TopBarHeight)], PlotCorners, ColN, RowN, 0
  ScaleCorners = PlotCorners
  PlotCorners  [*,2] = PlotCorners  [*,2] - BarWidth
  ScaleCorners [*,0] = ScaleCorners [*,2] - BarWidth
endif else begin					; here one scale for the entire plot
  GetCorners, [0.0,0.0,(1.0-BarWidth),(1.0-TopBarHeight)], PlotCorners,  ColN, RowN, 0
  GetCorners, [(1.0-BarWidth),0.0,1.0,(1.0-TopBarHeight)], ScaleCorners, 1,	  1, 0
endelse

if (QPlotName EQ 1) then begin				; this provides space for a label for each .glo
  LabelCorners = PlotCorners
  
  PlotCorners  [*,3] = PlotCorners  [*,3] - NameHeight
  LabelCorners [*,1] = LabelCorners [*,3] - NameHeight

  if (QSameScale EQ 2) then ScaleCorners [*,3] = ScaleCorners [*,3] - NameHeight
endif else begin
  LabelCorners = PlotCorners * 0.0
endelse

for XPlot = 0, (PlotN-1) do begin			; this ensures each .glo has correct width/length
  XGloPixels = (PlotCorners[XPlot,2] - PlotCorners[XPlot,0]) * XPixels
  YGloPixels = (PlotCorners[XPlot,3] - PlotCorners[XPlot,1]) * YPixels
  
  if (YGloPixels/XGloPixels LT float(LatN)/float(LongN)) then begin
    XNewPixels = YGloPixels * (float(LongN)/float(LatN))
    PlotCorners[XPlot,0] = PlotCorners[XPlot,0] + ((XGloPixels-XNewPixels)/(2.0*XPixels))
    PlotCorners[XPlot,2] = PlotCorners[XPlot,2] - ((XGloPixels-XNewPixels)/(2.0*XPixels))
  endif else begin
    YNewPixels = XGloPixels * (float(LatN)/float(LongN))
    PlotCorners[XPlot,1] = PlotCorners[XPlot,1] + ((YGloPixels-YNewPixels)/(2.0*YPixels))
    PlotCorners[XPlot,3] = PlotCorners[XPlot,3] - ((YGloPixels-YNewPixels)/(2.0*YPixels))
  endelse
endfor

TopBarSpec = fltarr(3,4) & TopBarSpec(*,*)=-1.0			; gets title bar positions
TopBarIndent = TopBarHeight/5.0

if (QBarInfo eq 1 or QBarInfo eq 2) then begin				; logo x:y = 834:184
	TopBarSpec(0,1) = (1.0-(TopBarHeight-TopBarIndent))
	TopBarSpec(0,3) = (TopBarHeight-(TopBarIndent*2.0))

	TopBarSpec(0,0) = TopBarIndent
	TopBarSpec(0,2) = TopBarSpec(0,3)*834.0/184.0
endif

TopBarSpec(1,*) = [0.5,(1.0-(TopBarHeight/1.4)),-1,-1]

if (QBarInfo eq 1 or QBarInfo eq 2) then begin
	TopBarSpec(2,0) = 1.0-TopBarIndent
	TopBarSpec(2,1) = 1.0-(TopBarIndent*2.5)
	TopBarSpec(2,2) = 1.0-TopBarIndent
	TopBarSpec(2,3) = 1.0-(TopBarHeight-TopBarIndent)
endif

; ******************************************************************************
; plot each .glo

SeasName = "" & SeasClash = 0

if (Choice EQ 21) then begin
  set_plot, 'X'	
  if (!D.Window EQ -1) then window, /free, xsize=(XPixels*32.0), ysize=(YPixels*32.0), $
  		xpos=(1000-(XPixels*32.0)), title="box plot"
endif else begin
  PSPath = ""
  print, "  > Enter the filepath of the .eps file: "
  read, PSPath

  Set_Plot, 'ps', /copy
  device, filename=PSPath, bits_per_pixel=8, xsize=XPixels, ysize=YPixels,/Color,/encapsulated
endelse

print, "  > Plotting..."

for XPlot = 0, (PlotN-1) do begin
;  if (QSeasons EQ 0) then begin				; these are not the four seasons
;    print, "  > Plot: ", (XPlot+1)
;  endif else if (QSeasons EQ 1) then begin		; these are the four seasons
;    print, "  > Plot: ", SeasonNames (XPlot)
;  endif else if (QSeasons EQ 2) then begin		; these are the 12 months
;    print, "  > Plot: ", MonthNames (XPlot)
;  endif
  
; ##### do standard or custom formatting #####

  if (QStanCust EQ 0) then begin
    if (QSeasons EQ 0 OR XPlot EQ 0) then begin
      if (QSameScale NE 1) then begin
        Info = strarr (10) & Info [*] = "missing"
        if (QSeasons GT 0) then Info[5] = "various"

        GrabInfo, FilePaths(XPlot), FileInfos(XPlot), Info, QView, QBarInfo
        
        if (Info[5] EQ "missing") then begin
          SeasName = "missing"
          SearchSeas, FilePaths(XPlot), FileInfos(XPlot), SeasName, SeasIndex=SeasIndex
        endif else begin
          SeasName = "missing"
          SearchSeas, FilePaths(XPlot), Info[5], SeasName, SeasIndex=SeasIndex
        endelse
      endif
      
      LoadVariCT, Info, ContourFactor, DataBounds, ContourN, TickVals, TickFormat, SeasIndex=SeasIndex  
    endif
  endif else begin
    Info = strarr (10) & Info [*] = "missing"
  
    LoadVariCT, Info, ContourFactor, DataBounds, ContourN, TickVals, TickFormat, SeasIndex=SeasIndex 
  endelse

; ##### prepare plot #####
  
  GetContours, DataBounds, ColorBounds, Contours, ContourColors, ContourN
  
  ExeRef = ((LabelCorners[XPlot,2] - LabelCorners[XPlot,0]) / 2.0) + LabelCorners[XPlot,0]
  WyeRef = ((LabelCorners[XPlot,3] - LabelCorners[XPlot,1]) / 3.0) + LabelCorners[XPlot,1]
  
  LoadGlo, Models(XPlot), GloSlice, FilePaths(XPlot), FileInfo, /Silent
  if (QPlotName EQ 1) then begin
    if (QSeasons EQ 0) then xyouts, ExeRef,WyeRef,PlotName[XPlot],font=5,charsize=(1.2*LabelSizeFactor), $
    					charthick=2.0,/normal,alignment=0.5
    if (QSeasons EQ 1) then xyouts, ExeRef,WyeRef,SeasonNames(XPlot),font=5,charsize=(1.2*LabelSizeFactor), $
    					charthick=2.0,/normal,alignment=0.5
    if (QSeasons EQ 2) then xyouts, ExeRef,WyeRef,MonthNames(XPlot),font=5,charsize=(1.2*LabelSizeFactor), $
    					charthick=2.0,/normal,alignment=0.5
  endif
  
  if            (QContourCell EQ 0) then begin
    TrimToScale, GloSlice, TrimSlice, DataBounds
  endif else if (QContourCell EQ 1) then begin
    TrimSlice = DataToColours (GloSlice, [DataBounds[0],DataBounds[1],ColorBounds[0],ColorBounds[1]], $
    		ModelsLongN(XPlot), ModelsLatN(XPlot), White)    		
  endif else if (QContourCell EQ 2) then begin
    TrimSlice = GloSlice
  endif
  
; ##### do plot #####
  
  if (QMapGrid EQ 0) then begin
    if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits 
    if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits
    if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits
  endif
  
  if (QOutline NE 0 AND QMapGrid EQ 0) then begin
   if (QOutline EQ -2) then begin
    if (QContourCell EQ 2) then begin
      if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Grey, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits 
      if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Grey, /noerase, $
  	position=PlotCorners[XPlot,*], /noborder, limit=Limits
      if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Grey, /noerase, $
   	position=PlotCorners[XPlot,*], /noborder, limit=Limits
      Map_Continents, /coasts, /fill, color=Grey
    endif else begin
      if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /continents, /noborder, limit=Limits 
      if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /continents, /noborder, limit=Limits
      if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
   	position=PlotCorners[XPlot,*], /continents, /noborder, limit=Limits
    endelse
   endif else begin
    PlotPoly, QOutline, linestyle=0, color=Black
   endelse
  endif
  
  if (QGrid EQ 1 AND QMapGrid EQ 0) then begin
    if (ViewKeyWords[QView,0] EQ 'Lamb') then Map_set, 0, 0, /lambert, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /grid, glinestyle=0, /noborder, limit=Limits 
    if (ViewKeyWords[QView,0] EQ 'Cyli') then Map_set, 0, 0, /cylindrical, color=Black, /noerase, $
  	position=PlotCorners[XPlot,*], /grid, glinestyle=0, /noborder, limit=Limits
    if (ViewKeyWords[QView,0] EQ 'Hamm') then Map_set, 0, 0, /hammer, color=Black, /noerase, $
   	position=PlotCorners[XPlot,*], /grid, glinestyle=0, /noborder, limit=Limits
  endif
  
  GetContourVec, Models(XPlot), Lons, Lats
  
  if 		(QContourCell EQ 0) then begin		; colour contours
   if (QMapGrid EQ 0) then begin
    contour, TrimSlice, Lons, Lats, levels=Contours, c_colors=ContourColors, /cell_fill, $
  	/overplot,/noerase
   endif else begin
    contour, TrimSlice, xrange=[Limits(1),Limits(3)], xstyle=5, ystyle=5, xticks=1, yticks=1, $
    	yrange=[Limits(0),Limits(2)], color=Black, $
    	/noerase, position=PlotCorners[XPlot,*], /nodata
    contour, TrimSlice, Lons, Lats, levels=Contours, c_colors=ContourColors, /cell_fill, $
  	/overplot,/noerase
   endelse
  endif else if (QContourCell EQ 1) then begin		; boxes
   if (QMapGrid EQ 0) then begin
     Map = map_image (TrimSlice, startx, starty, xs, ys, $
	lonmin=Lons[0], lonmax=Lons[ModelsLongN(XPlot)-1], latmin=Lats[0], latmax=Lats[ModelsLatN(XPlot)-1], $
	compress=1, min_value=(MissVal+1), missing=white)
	
     TV, Map, startx, starty, xsize=xs, ysize=ys
   endif else begin
     TV,TrimSlice,PlotCorners[XPlot,0],PlotCorners[XPlot,1],/normal, $
      	xsize=(PlotCorners[XPlot,2]-PlotCorners[XPlot,0]),ysize=(PlotCorners[XPlot,3]-PlotCorners[XPlot,1])
   endelse
  endif else if (QContourCell EQ 2) then begin  	; b/w contours
   if (QMapGrid EQ 0) then begin
    contour, TrimSlice, Lons, Lats,levels=Contours,c_labels=Contours, $
    		min_value=(MissVal+1),/overplot,/noerase
   endif else begin
    contour, TrimSlice, xrange=[Limits(1),Limits(3)], xstyle=5, ystyle=5, xticks=1, yticks=1, $
    	yrange=[Limits(0),Limits(2)], color=Black, $
    	/noerase, position=PlotCorners[XPlot,*], /nodata
    contour, TrimSlice, Lons, Lats,levels=Contours,c_labels=Contours, $
    		min_value=(MissVal+1),/overplot,/noerase
   endelse
  endif
  
  if (QSameScale EQ 2) then DrawConScale,ScaleCorners[XPlot,*],Contours,ContourColors,$
    	TickVals,TickFormat,TextSize,TextThick,Info,XPixels,YPixels,QBottomRight
endfor

if (QSameScale EQ 1) then DrawConScale,ScaleCorners[0,*],Contours,ContourColors,$
    	TickVals,TickFormat,TextSize,TextThick,Info,XPixels,YPixels,QBottomRight

if (QBarInfo eq 1 or QBarInfo eq 2) then begin
 if (Choice EQ 22) then begin
  read_jpeg, '/cru/u2/f709762/goglo/ref/tyn-logo-torok.jpg', Logo, LogoCT, colors=256, dither=1, /two_pass_quantize
  tvlct, Red,Green,Blue, /get						; get old colors
  tvlct, LogoCT(*,0), LogoCT(*,1), LogoCT(*,2)
  tv, Logo, TopBarSpec[0,0],TopBarSpec[0,1],xsize=TopBarSpec[0,2],ysize=TopBarSpec[0,3],/normal
  tvlct, Red,Green,Blue							; restore old colors
 endif else begin
  xyouts,TopBarSpec[0,0],TopBarSpec[0,1],'logo',/normal,alignment=0,charsize=1,charthick=1,font=5
 endelse
 
 Author = 'Tim Mitchell'
 DateTime  = systime ()
 DateStamp = strmid(DateTime,4,7) + strmid(DateTime,20,4)
 xyouts,TopBarSpec[2,0],TopBarSpec[2,1],Author,/normal,alignment=1,charsize=1,charthick=1,font=5
 xyouts,TopBarSpec[2,2],TopBarSpec[2,3],DateStamp,/normal,alignment=1,charsize=1,charthick=1,font=5
endif

xyouts,TopBarSpec[1,0],TopBarSpec[1,1],MainTitle,/normal,alignment=0.5,charsize=2.5,charthick=3,font=5


; ******************************************************************************
; end plot

if (Choice EQ 22) then device, /close

set_plot, 'X'	

; ******************************************************************************
; end program

print, ""

end