; plotglofour.pro
; written by Tim Mitchell on 17.12.99 using plotglo.pro as the basis
; last modified on 01.08.00
; program that loads multiple .glo files and plots the spatial slices

; *********************************************************************
; initialising of variables

initx

MissVal  = -999.0
DatumFlt  = 0.0D
DatumInt  = 0L & DatumInt1 = 0L & DatumInt2 = 0L & DatumInt3 = 0L
DatumCha  = "bunnywhoops"
ImageLatMin  = 0.0 & ImageLatMax  = 0.0 & ImageLongMin = 0.0 & ImageLongMax = 0.0

; *********************************************************************
; select and load multiple .glo files from a single model

SelectModel, ModelChosen, ModelFilePath, ModelLongN, ModelLatN
  
GloFileN = 0
print, "  > Enter the number of .glo files to load (up to four): "
read, GloFileN

Slice         = fltarr (4,ModelLongN,ModelLatN)
Slice [*,*,*] = MissVal

for XGlo = 0, (GloFileN-1) do begin
  print, "  > Load .glo file number: ", (XGlo+1)
    
  LoadGlo, ModelFilePath, GloSlice
    
  for XLong = 0, (ModelLongN-1) do begin
      for XLat = 0, (ModelLatN-1) do begin
        Slice (XGlo,XLong,XLat) = GloSlice (XLong,XLat)
      endfor
  endfor    
endfor

GetImageBounds, ModelChosen, ImageLatMin, ImageLatMax, ImageLongMin, ImageLongMax
GetViewBounds, ViewN, ViewName, ViewKeyWords, ViewBounds
GetScales, ScaleN, ScaleName, ScaleColTab, ScaleSeg, ScaleLimits

; *********************************************************************
; initialising of variables

Black = 0
Grey  = 254
White = 255

LevelIncrement = 0.0D
ColorIncrement = 0.0D

ColourSlice       = intarr (ModelLongN,ModelLatN)
ColourSlice (*,*) = Grey

ImageSlice4        = intarr (4,ModelLongN,ModelLatN)
ImageSlice4(*,*,*) = MissVal

RangeFract = 0.0D
ColSelect  = 0
	
SetGlo        = intarr (4)
SetView       = intarr (4)
SetScale      = intarr (4)
SetGlo    [*] = MissVal
SetView   [*] = MissVal
SetScale  [*] = MissVal

CurrentScale = fix(MissVal)

Compression = 1

;******************************************************************************
; main loop

Choice = 78

while (Choice LT 99) do begin
	
case 1 of
	
;******************************************************************************
; select defaults for SetGlo SetView SetScale

(Choice EQ 1): begin	

QGlo   = 0
QView  = 0
QScale = 0

QGlo = -2
print, ""
print, "  > Select the default assignment of .glo files (-1=list): "
while (QGlo LT 1 OR QGlo GT 3) do begin
 read, QGlo
 if (QGlo EQ -1) then begin
  print, 1, " : ", "each file to a plot", format='(i4,a3,a25)'
  print, 2, " : ", "first 2 files to 4 plots", format='(i4,a3,a25)'
  print, 3, " : ", "first file for all plots", format='(i4,a3,a25)'
 endif
endwhile

if 		(QGlo EQ 1) then begin
  for XGloFile = 0, (GloFileN-1) do begin
    SetGlo [XGloFile]    = XGloFile
  endfor
endif else if 	(QGlo EQ 2) then begin
  for XPlot = 0, 1 do begin
    SetGlo [XPlot]    = XPlot
    SetGlo [XPlot+2]  = XPlot
  endfor
endif else if 	(QGlo EQ 3) then begin
  for XPlot = 0, 3 do begin
    SetGlo [XPlot]    = 0
  endfor
endif

QView = -2
print, "  > Select the default view (-1=list): "
while (QView LT 0 OR QView GT (ViewN-1)) do begin
 read, QView
 if (QView EQ -1) then begin
  for XView = 0, (ViewN-1) do begin
   print, XView, " : ", ViewName[XView], format='(i4,a3,a10)'
  endfor
 endif
endwhile

QScale = -2
print, "  > Select the default scale (-1=list): "
while (QScale LT 0 OR QScale GT (ScaleN-1)) do begin
 read, QScale
 if (QScale EQ -1) then begin
  for XScale = 0, (ScaleN-1) do begin
    print, XScale, " : ", ScaleName[XScale], ScaleLimits[XScale,0], ScaleLimits[XScale,1], $
  		 format='(i4,a3,a10,2f10.2)'
  endfor
 endif
endwhile

for XPlot = 0, 3 do begin
   SetView   [XPlot] = QView
   SetScale  [XPlot] = QScale
endfor

end

;******************************************************************************
; select .glo files to plot

(Choice EQ 2): begin	

print, ""
print, "  > Select a plot, .glo file (99,99=end): "

QPlot   = -1
QGlo    = -1

while (QPlot NE 99 AND QGlo NE 99) do begin
  if (QPlot GE 0 AND QPlot LT 4) then begin
    if (QGlo GE 0 AND QGlo LT GloFileN) then begin
      SetGlo [QPlot] = QGlo
    endif
  endif
  
  read, QPlot, QGlo
endwhile

end

;******************************************************************************
; select view

(Choice EQ 3): begin	

print, ""
print, "  > Select a plot, view (-1,-1=list, 99,99=end): "

QPlot = -2
QView = -2

while (QPlot NE 99 AND QView NE 99) do begin
  if (QPlot EQ -1 AND QView EQ -1) then begin
    for XView = 0, (ViewN-1) do begin
      print, XView, " : ", ViewName[XView], format='(i4,a3,a10)'
    endfor  
  endif
  
  if (QPlot GE 0 AND QPlot LT 4) then begin
    if (QView GE 0 AND QView LT ViewN) then begin
      SetView [QPlot] = QView
    endif
  endif
  
  read, QPlot, QView
endwhile

end

;******************************************************************************
; select scale

(Choice EQ 4): begin	

print, ""
print, "  > Select a plot, scale (-1,-1=list, 99,99=end): "

QPlot  = -2
QScale = -2

while (QPlot NE 99 AND QScale NE 99) do begin
  if (QPlot EQ -1 AND QScale EQ -1) then begin
    for XScale = 0, (ScaleN-1) do begin
      print, XScale, " : ", ScaleName[XScale], ScaleLimits[XScale,0], ScaleLimits[XScale,1], $
  		 format='(i4,a3,a10,2f10.2)'
  endfor
  
  endif
  
  if (QPlot GE 0 AND QPlot LT 4) then begin
    if (QScale GE 0 AND QScale LT ScaleN) then begin
      SetScale [QPlot] = QScale
    endif
  endif
  
  read, QPlot, QScale
endwhile

end

;******************************************************************************
; reverse all scales

(Choice EQ 5): begin	

print, "  > Reverse all scales ? (1=no,2=yes)"
read, QReverseAll  

if (QReverseAll EQ 2) then begin

 for XScale = 0, (ScaleN-1) do begin
  DatumFlt = ScaleLimits[XScale,2]
  ScaleLimits[XScale,2] = ScaleLimits[XScale,3]
  ScaleLimits[XScale,3] = DatumFlt
 endfor

endif else begin

 for XScale = 0, (ScaleN-1) do begin
  print, XScale, " : ", ScaleName[XScale], ScaleLimits[XScale,0], ScaleLimits[XScale,1], $
  		 format='(i4,a3,a10,2f10.2)'
 endfor
 
 QReverseOne = -1
 print, "  > Select scales to reverse (99=end): "
 while (QReverseOne NE 99) do begin
   if (QReverseOne GE 0 AND QReverseOne LT ScaleN) then begin
    DatumFlt                   = ScaleLimits[QReverseOne,2]
    ScaleLimits[QReverseOne,2] = ScaleLimits[QReverseOne,3]
    ScaleLimits[QReverseOne,3] = DatumFlt
   endif
   
   read, QReverseOne 
 endwhile
 
endelse

end

;******************************************************************************
; construct image

(Choice EQ 10): begin	

for XPlot = 0, 3 do begin
 if (SetGlo[XPlot] NE MissVal) then begin
  if (SetScale[XPlot] NE MissVal) then begin
   
   ColourSlice = DataToColours (reform(Slice(SetGlo[XPlot],*,*)), $
   			reform(ScaleLimits(SetScale[XPlot],*)), ModelLongN, ModelLatN, Grey)
   
   for XLong = 0, (ModelLongN-1) do begin
     for XLat = 0, (ModelLatN-1) do begin
       ImageSlice4 (XPlot,XLong,XLat) = ColourSlice (XLong,XLat)  
     endfor   
   endfor
   
   print, "  > Image prepared, scale, .glo: ", XPlot, " ", ScaleName[SetScale[XPlot]], " ", $
    			SetGlo[XPlot], format="(A32,I3,A1,A10,A1,I3)"
  
  endif else begin
   print, "  > No scale assigned to plot:   ", XPlot, format="(A32,I3)"
  endelse
 
 endif else begin
   print, "  > No .glo assigned to plot:    ", XPlot, format="(A32,I3)"
 endelse
endfor

end

;******************************************************************************
; plot image to screen (11) or .eps (12)

(Choice EQ 11 OR Choice EQ 12): begin	

if (Choice EQ 11) then begin
  set_plot, 'X'	
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, /Color, /encapsulated
endelse

!P.Multi = [0, 2, 2, 0, 0]

for XPlot = 0, 3 do begin
 if (SetGlo(XPlot) NE MissVal) then begin
  if (SetScale(XPlot) NE MissVal) then begin
   if (SetView(XPlot) NE MissVal) then begin
    
    if (SetScale[XPlot] NE CurrentScale OR XPlot EQ 0) then begin
     
     MyColTab = ScaleColTab (SetScale[XPlot])
     LoadAnyCT, MyColTab

     GetScaleVectors,reform(ScaleLimits(SetScale[XPlot],*)),ScaleLevels,ScaleColors,ScaleSeg[SetScale[XPlot]]
     
     CurrentScale = SetScale[XPlot]
     
    endif
    
    !P.Multi = [(4-XPlot), 2, 2, 0, 0]

    if (Choice EQ 11) then begin
      Set_Plot, 'X'
    endif else if (Choice EQ 12) then begin
      Set_Plot, 'ps', /copy
    endif
    
    scale_vert, levels=ScaleLevels, c_colors=ScaleColors, noextremes=[" "," "], chwid=2, ytickformat='(f7.2)'

    if (ViewKeyWords[SetView[XPlot],0] EQ 'Cyli') then begin
      Map_set, 0, 20, /cylindrical, color=White, /noerase, /advance, $
	limit=[ViewBounds[SetView[XPlot],2],ViewBounds[SetView[XPlot],0],$
	       ViewBounds[SetView[XPlot],3],ViewBounds[SetView[XPlot],1]] 
    endif else if (ViewKeyWords[SetView[XPlot],0] EQ 'Lamb') then begin
      Map_set, 0, 20, /lambert, color=White, /noerase, /advance, $
	limit=[ViewBounds[SetView[XPlot],2],ViewBounds[SetView[XPlot],0],$
	       ViewBounds[SetView[XPlot],3],ViewBounds[SetView[XPlot],1]] 
    endif else if (ViewKeyWords[SetView[XPlot],0] EQ 'Hamm') then begin
      Map_set, 0, 20, /hammer, color=White, /noerase, /advance, $
	limit=[ViewBounds[SetView[XPlot],2],ViewBounds[SetView[XPlot],0],$
	       ViewBounds[SetView[XPlot],3],ViewBounds[SetView[XPlot],1]] 
    endif
    
    ImageSlice4 [*,48,57] = Grey
    PlotImageSlice = fltarr (ModelLongN,ModelLatN)    
    PlotImageSlice = reform (ImageSlice4[XPlot,*,*])
    
    Map = map_image (PlotImageSlice, startx, starty, xs, ys, $
	lonmin=ImageLongMin, lonmax=ImageLongMax, $
	latmin=ImageLatMin,  latmax=ImageLatMax,  $
	compress=Compression, min_value=(MissVal+1), missing=grey)
	
    TV, Map, startx, starty, xsize=xs, ysize=ys

    if 		  (ViewKeyWords[SetView[XPlot],1] EQ 'Cont' AND $
    			ViewKeyWords[SetView[XPlot],2] EQ 'Coun') then begin
      Map_continents, /continents, /countries, color=Black, /noerase
    endif else if (ViewKeyWords[SetView[XPlot],1] EQ 'Coas' AND $
    			ViewKeyWords[SetView[XPlot],2] EQ 'Coun') then begin
      Map_continents, /coasts, /countries, color=Black, /noerase
    endif else if (ViewKeyWords[SetView[XPlot],1] EQ 'Coas' AND $
    			ViewKeyWords[SetView[XPlot],2] EQ 'Bare') then begin
      Map_continents, /coasts, color=Black, /noerase
    endif
    
    scale_vert, /off

    print, "  > Plot, scale, region, .glo: ", XPlot, " ", ScaleName[SetScale[XPlot]], " ", $
    			ViewName[SetView[XPlot]], " ", SetGlo[XPlot], $
    			format="(A30,I3,A1,A10,A1,A10,A1,I3)"
    
    LastPlot = XPlot
    
   endif else begin
    print, "  > No view assigned to plot:  ", XPlot, format="(A30,I3)"
   endelse

  endif else begin
    print, "  > No scale assigned to plot: ", XPlot, format="(A30,I3)"
  endelse

 endif else begin
    print, "  > No .glo assigned to plot:  ", XPlot, format="(A30,I3)"
 endelse

endfor

if (Choice EQ 12) then 	device, /close

if (Choice EQ 11) then begin
  UseMouse, reform(ViewBounds[LastPlot,*]), ModelChosen, XLong, XLat, XButton
  print, XLong, XLat, XButton
endif

set_plot, 'X'	

end

; *********************************************************************
; end of main loop

else: begin
	print, ""
	print, "  >  1 : select defaults"
	print, "  >  2 : select glo files to plot"
	print, "  >  3 : select regions to view"
	print, "  >  4 : select scales"
	print, "  >  5 : reverse scales"
	print, "  > 10 : construct image"
	print, "  > 11 : plot image to screen"
	print, "  > 12 : plot image to .eps"
	print, "  > 99 : end"	
end

endcase
        
        print, ""
	print, "  > Main menu: Enter a number:"
	read, Choice
	
endwhile

print, ""
!P.Multi = [0, 1, 1, 0, 0]

end