;------------------------------------------------------------------
;------------------------------------------------------------------
;------------------------------------------------------------------
;------------------------------------------------------------------
;------------------------------------------------------------------
;------------------------------------------------------------------
;------------------------------------------------------------------
; Main Program:

;--------------------------
; Variables to Set

mr = 40.0  ; this is the Maximum RANGE (i.e. 90 - Minimum Latitude)

psfile_base = "test" ; this is the base name of the postscript files

iFirst = 60  ; the first time to plot (0-1439)
iLast  = 63  ; the last time to plot  (iFirst - 1439)
iStep  = 1   ; The cadence of the plots

; Set this to the maximum potential for the color bar
;  - put 0 if you want an automatic selection
Max_Potential = 30.0  

; Set this to 1 if you want a zero line
Contour_Zero_Line = 1

;--------------------------

filelist = findfile("b*.save")

filename = filelist(0)

print, "restoring file : ", filename

restore, filename

potential = reform(data(*,0,*,*))

nlats = n_elements(lats)
nmlts = n_elements(mlts)

lat2d = fltarr(nmlts,nlats)
lon2d = fltarr(nmlts,nlats)

for i=0,nlats-1 do lon2d(*,i) = mlts*!pi/12.0 - !pi/2.0
for i=0,nmlts-1 do lat2d(i,*) = lats

x = (90.0-lat2d)*cos(lon2d)
y = (90.0-lat2d)*sin(lon2d)

nFields = n_elements(fields)

;for i=0,nFields-1 do print, string(i,format="(i2)")+". ",fields(i)

;------------------------------------------------
; Figure out levels:

joule = 0
nLevels = 31

if (Max_Potential eq 0) then begin

    mini_1  = min(potential)
    maxi_1  = max(potential)
    if (mini_1 ge 0.0) then begin
        mini_1  = 0.0
        range_1 = maxi_1
    endif else begin
        if (not joule) then range_1 = max([abs(mini_1),maxi_1])*2.0 $
        else range_1 = maxi_1 - mini_1
    endelse
    dc      = 10.0^fix(alog10(range_1/100.0))
    factor  = 1.0
    while (range_1 gt dc*(nLevels-1)*factor) do factor=factor+0.1
    dc = factor*dc

endif else begin

  mini_1 = -Max_Potential
  dc = 2*Max_Potential / (nLevels-1)

endelse

if (mini_1) eq 0.0 then 						$
  levels = findgen(nLevels)*dc						$
else levels = (findgen(nLevels) - (nLevels-1)/2)*dc

makect,'mid'
ncolors = 255
clevels = (ncolors-1)*findgen(nLevels)/(nLevels-1) + 1

nLevels = (nLevels-1)/4+1
if (mini_1) eq 0.0 then 					$
  levels_lines = findgen(nLevels)*dc*4			$
else begin
    levels_lines = (findgen(nLevels/2)+1)*dc*4
    if (Contour_Zero_Line) then begin
        levels_lines = [-reverse(levels_lines),0.0,levels_lines]
    endif else begin
        levels_lines = [-reverse(levels_lines),levels_lines]
    endelse
endelse

loc = where(lat2d(0,*) gt 90.0-mr)

for n=iFirst, iLast, iStep do begin

    c_r_to_a, itime, time(n)
    c_a_to_s, itime, stime

    sn = chopr('0000'+tostr(n),4)
    psfile = psfile_base+'_'+sn+'.ps'

    setdevice, psfile, 'p', 5, 0.95

    !p.charsize = 1.2
    !p.thick = 2

    ppp = 1
    space = 0.05
    pos_space, ppp, space, sizes
    get_position, ppp, space, sizes, 0, pos

    contour, potential(n,*,loc),x(*,loc),y(*,loc),	$
      /follow, xstyle = 5, ystyle = 5,				$
      xrange = [-mr,mr],yrange=[-mr,mr], levels = levels, 		$
      pos = pos, /noerase, /cell_fill, c_color = clevels

    contour, potential(n,*,loc),x(*,loc),y(*,loc),	$
      /follow, xstyle = 5, ystyle = 5,				$
      xrange = [-mr,mr],yrange=[-mr,mr], levels = levels_lines,		$
      pos = pos, /noerase, c_linestyle = 3.0*(levels_lines lt 0.0)

;--------------------------------------------------------------
; Draw the MLT grid
;--------------------------------------------------------------

    plotmlt, mr, /no06

;--------------------------------------------------------------
; Tell us what time it is
;--------------------------------------------------------------

    dy  = float(!d.y_ch_size)/float(!d.y_size)

    xyouts, pos(0), pos(3)-0.75*dy, strmid(stime,10,5)+' UT', 		$
      /norm

    ct_left = pos(2)+0.01
    ct_right = ct_left + 0.02
    ctpos = [ct_left, pos(1), ct_right, pos(3)]
    plotct, ncolors, ctpos, mm(levels), 'kV', /right

    mini = min(potential(n,*,loc))
    maxi = max(potential(n,*,loc))

    mini_string = 'Min:'+tostr(round(mini))
    maxi_string = 'Max:'+tostr(round(maxi))

    up_shift    = 0.5*dy
    xyouts, pos(0), pos(1)+up_shift, mini_string, /norm
    xyouts, pos(2), pos(1)+up_shift, maxi_string, /norm, 	$
      alignment=1.0

    closedevice

endfor

end