; MODIFICATION HISTORY:
;	01-Oct-2008 Removed mdsconnect [BD]
;--------------------------------------------------------------
PRO oplot_bloops_subset,s,size,loop_indices
isub	= loop_indices-1

x0	= [-1.00,+0.50,+0.50,+1.00,+0.50,+0.50,-1.00,-1.00]
y0	= [+0.05,+0.05,-0.25,+0.00,+0.25,-0.05,-0.05,+0.05] * 2.0
r0	= SQRT(x0^2+y0^2)*size
ang0	= ATAN(y0,x0)

nels	= N_ELEMENTS(s.bloop.data.r[isub])
FOR i=0,nels-1 DO BEGIN
  angle	= s.bloop.data.pang(isub[i])*!PI/180.0
  xc	= s.bloop.data.r(isub[i])
  yc	= s.bloop.data.z(isub[i])
  x	= xc+r0*COS(ang0+angle)
  y	= yc+r0*SIN(ang0+angle)
;  POLYFILL,x,y,color=s.bloop.data.color
  POLYFILL,x,y,color=35

ENDFOR

RETURN
END

;---------------------------------------------------------------------------------------
PRO oplot_bloops,s,size
oplot_bloops_subset,s,size,LINDGEN(N_ELEMENTS(s.bloop.data.r))+1
RETURN
END

;---------------------------------------------------------------------------------------
PRO plot_nstx_poloidal_mirnovs

file = '/p/nstxusr1/nstx-users/jmenard/analysis/plot_efit/poloidal_array_data.sav'
RESTORE,file
s	= {bloop:{data:position}}

siz	= AVERAGE(s.bloop.data.r)/12.0

oplot_bloops,s,siz

RETURN
END


;--------------------------------------------------------------------------------------
pro plot_nstx_hhfw_antenna
box_z = [-.00635,.00635,0.00635,-.00635,-.00635]
box_r = [0.,0.,.0508,.0508,0.]
r_curv = 23.15*.0254
r_cent = 40.15*.0254
t_boron = -.03
d_strap = 0.625*.0254
t_strap = 0.375*.0254
theta_max = !pi*11.25/180
r_wall = 1.7
;
; Boron
;
theta = theta_max*(findgen(21)/10-1.)
offset_parallel = (0.5+findgen(13))*.0254
offset_r = -1.0*sin(theta_max)*offset_parallel
offset_z = cos(theta_max)*offset_parallel
r_start = (r_curv + t_boron)*cos(theta_max)+r_cent
z_start = (r_curv + t_boron)*sin(theta_max)
r_boron = [r_wall,r_start+reverse(offset_r), $
  (r_curv+t_boron)*cos(theta)+r_cent,r_start+offset_r,r_wall]
z_boron = [-1.*z_start-max(offset_z),-1.*z_start-reverse(offset_z), $
  (r_curv+d_strap)*sin(theta),z_start+offset_z,z_start+max(offset_z)]
polyfill,r_boron,z_boron,color=13
oplot,r_boron,z_boron
;
;Antenna strap
;
theta = theta_max*(findgen(21)/10-1.)
offset_parallel = (0.5+findgen(12))*.0254
offset_r = -1.0*sin(theta_max)*offset_parallel
offset_z = cos(theta_max)*offset_parallel
r_start = (r_curv + d_strap)*cos(theta_max)+r_cent
z_start = (r_curv + d_strap)*sin(theta_max)
r_antenna = [r_start+reverse(offset_r), $
  (r_curv+d_strap)*cos(theta)+r_cent,r_start+offset_r, $
   r_start+reverse(offset_r)+t_strap, $
   (r_curv+d_strap+t_strap)*cos(theta)+r_cent, $
   r_start+offset_r+t_strap]
z_antenna = [-1.*z_start-reverse(offset_z), $
  (r_curv+d_strap)*sin(theta),z_start+offset_z, $
   z_start+reverse(offset_z), $
   (r_curv+d_strap+t_strap)*reverse(sin(theta)), $
    -1.*z_start-offset_z]
polyfill,r_antenna,z_antenna,color=11
oplot,r_antenna,z_antenna
;
; Horizontal part of straps
;
polyfill,r_start+min(offset_r)+[0.,0.,.14,.14,0.], $
 z_start+max(offset_z)+[0.,.01,.01,0.,0.],color=11
oplot,r_start+min(offset_r)+[0.,0.,.14,.14,0.], $
 z_start+max(offset_z)+[0.,.01,.01,0.,0.]
polyfill,r_start+min(offset_r)+[0.,0.,.13,.13,0.], $
 -1.0*(z_start+max(offset_z)+[0.,.01,.01,0.,0.]),color=11
oplot,r_start+min(offset_r)+[0.,0.,.13,.13,0.], $
 -1.0*(z_start+max(offset_z)+[0.,.01,.01,0.,0.])
;
;Faraday shields on curved part of antenna
;
theta = !pi*2.5*(findgen(9) - 4.0)/180
box2_r = cos(theta)#(r_curv+box_r)-sin(theta)#box_z
box2_z = sin(theta)#(r_curv+box_r)+cos(theta)#box_z
for i=0,8 do polyfill,box2_r(i,*)+r_cent,box2_z(i,*),$
 color=14
;
;Faraday shields on straight parts of antenna
;
box3_r = cos(theta_max)*(r_curv+box_r)-sin(theta_max)*box_z
box3_z = sin(theta_max)*(r_curv+box_r)+cos(theta_max)*box_z
offset_parallel = (0.5+findgen(11))*.0254
offset_r = -1.0*sin(theta_max)*offset_parallel
offset_z = cos(theta_max)*offset_parallel
for i = 0,10 do polyfill,box3_r+r_cent+offset_r[i],box3_z+offset_z[i], $
  color=14
for i = 0,10 do polyfill,box3_r+r_cent+offset_r[i], $
 -1.0*(box3_z+offset_z[i]),color=14
end

;--------------------------------------------------------------------------------------
pro plot_nstx_tiles,reverse=reverse
if not keyword_set(reverse) then s=1. else s=-1.
;
; Centerstack tiles
;
polyfill,s*[.1711,.1836,.1836,.1711,.1711], $
 [-1.019,-1.019,1.019,1.019,-1.019],$
 color=0
;
; Center Stack Angle flanges
;
polyfill,s*[.1811,.1811,.282,.282,.1811], $
 [1.004,1.034,1.209,1.179,1.004],color=0
polyfill,s*[.1811,.1811,.282,.282,.1811], $
 [-1.004,-1.034,-1.209,-1.179,-1.004], $
 color=0
;
; PF1A casings
;
polyfill,s*[.26,.285,.285,.26,.26]-.004, $
 [1.164,1.169,1.6,1.6,1.164],color=0
polyfill,s*[.26,.285,.285,.26,.26]-.004, $
 -1.*[1.164,1.169,1.6,1.6,1.164],color=0
;
; Beveled corner tiles between CS and inner divertor
;
polyfill,s*[.266,.266,.287,.318,.266], $
 [1.56,1.6,1.63,1.63,1.56],color=0
polyfill,s*[.266,.266,.287,.318,.266], $
 -1.0*[1.56,1.6,1.63,1.63,1.56],color=0
;
; Inner divertor
;
polyfill,s*[.287,.572,.572,.288,.288], $
 [1.605,1.605,1.655,1.655,1.605],color=0
polyfill,s*[.287,.572,.572,.288,.288], $
 -1.0*[1.605,1.605,1.655,1.655,1.605], $
 color=0
;
;Outer divertor
;
polyfill,s*[.625,.625,1.205,1.205,.625], $
 [1.63,1.655,1.427,1.402,1.63],color=0
polyfill,s*[.625,.625,1.205,1.205,.625], $
 -1.0*[1.63,1.657,1.429,1.402,1.63], $
 color=0
;
;Secondary passive plates
;
polyfill,s*[1.12,1.337,1.317,1.10,1.12], $
 [1.36,1.05,1.034,1.344,1.36], $
 color=0
polyfill,s*[1.12,1.337,1.317,1.10,1.12], $
 -1.0*[1.36,1.05,1.034,1.344,1.36], $
 color=0
;Primary passive plates
polyfill,s*[1.361,1.51,1.486,1.337,1.361], $
 [1.008,.553,.545,1.0,1.008], $
 color=0
polyfill,s*[1.361,1.51,1.486,1.337,1.361], $
 -1.0*[1.008,.553,.545,1.0,1.008], $
 color=0
end

;--------------------------------------------------------------------------------------
PRO plot_nstx_vessel_day0,print=print

title	= '!5NSTX DAY 0 Conducting Elements'

xs	= 6.5
ys	= 10.0

if keyword_set(print) then begin
 SET_PLOT,'ps'
 xo      = 0.5*(8.5-xs)
 yo      = 0.5*(11.-ys)
 DEVICE,/portrait,/inches,xoffset=xo,xsize=xs,yoffset=yo,ysize=ys
 DEVICE,/color
endif

RESTORE,'nstx_day0_conductors.dat'

plot_nstx_vessel,s0,title

if keyword_set(print) then begin
 DEVICE,/close
 SET_PLOT,'x'
 PRINT,'PS image stored in idl.ps'
endif

RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_nstx_vessel,s,title

xmin	= 0
xmax	= 2.4
ymin	= -2.4
ymax	= +2.4
size	= 0.7

plot_element,s,xmin,xmax,ymin,ymax,size,title

RETURN
END

 ;--------------------------------------------------------------------------------------
PRO shade_elements_resist,s,boxcolor
r	= s.full.cedata.r
z	= s.full.cedata.z
w	= s.full.cedata.w/2.
h	= s.full.cedata.h/2.
rv	= [[r],[r+w],[r-w],[r-w],[r+w]]
zv	= [[z],[z+h],[z+h],[z-h],[z-h]]
FOR i=0, N_ELEMENTS(r)-1 DO BEGIN
  rvec	= TRANSPOSE([[rv(i,1:4)],[rv(i,1)]])
  zvec	= TRANSPOSE([[zv(i,1:4)],[zv(i,1)]])
  color = s.full.cedata.eta.color[s.full.cedata.Iresistance(i)]
  POLYFILL,rvec,zvec,color=color,NOCLIP=0
  OPLOT,rvec,zvec,color=boxcolor
ENDFOR
RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_element,s,xmin,xmax,ymin,ymax,size,title

;!P.FONT = 0
;DEVICE,/HELVETICA,/BOLD, FONT_INDEX=5
;!p.background=255
plot,[xmin,xmax],[ymin,ymax],psym=3,xstyle=1,ystyle=1,charsize=1.3,$
	xtitle='!5R(m)',ytitle='!5Z(m)',title=title
shade_elements_resist,s,0

RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_elements,file,xmin,xmax,ymin,ymax,size
title	= '!5Outline of Conducting Elements'
s0 	= read_conductor_data(file)
plot_element,s0,xmin,xmax,ymin,ymax,size,title
RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_nstx_vessel_day0,print=print

title	= '!5NSTX DAY 0 Conducting Elements'

xs	= 6.5
ys	= 10.0

if keyword_set(print) then begin
 SET_PLOT,'ps'
 xo      = 0.5*(8.5-xs)
 yo      = 0.5*(11.-ys)
 DEVICE,/portrait,/inches,xoffset=xo,xsize=xs,yoffset=yo,ysize=ys
 DEVICE,/color
endif

RESTORE,'nstx_day0_conductors.dat'

plot_nstx_vessel,s0,title

if keyword_set(print) then begin
 DEVICE,/close
 SET_PLOT,'x'
 PRINT,'PS image stored in idl.ps'
endif

RETURN
END

;--------------------------------------------------------------------------------------
FUNCTION nstx_conductor_data_structure,shot
file	= '/u/jmenard/idl/efit/nstx/nstx_day1_conductors.dat'
IF shot GE 109078 THEN file = '/u/jmenard/idl/efit/nstx/nstx_chiabsv2_conductors.dat'
IF shot GE 115000 THEN file = '/u/jmenard/idl/efit/nstx/nstx_newpf1a_conductors.dat'
RESTORE,file
RETURN, s0
END

;--------------------------------------------------------------------------------------
PRO plot_nstx_vessel,s,title

xmin	= 0
xmax	= 2.4
ymin	= -2.4
ymax	= +2.4
size	= 0.7

plot_element,s,xmin,xmax,ymin,ymax,size,title

RETURN
END

;--------------------------------------------------------------------------------------
PRO shade_elements_resist,s,boxcolor
r	= s.full.cedata.r
z	= s.full.cedata.z
w	= s.full.cedata.w/2.
h	= s.full.cedata.h/2.
rv	= [[r],[r+w],[r-w],[r-w],[r+w]]
zv	= [[z],[z+h],[z+h],[z-h],[z-h]]
vessel_colors = [9,10,11,15,1,16,10,16,11]
FOR i=0, N_ELEMENTS(r)-1 DO BEGIN
  rvec	= TRANSPOSE([[rv(i,1:4)],[rv(i,1)]])
  zvec	= TRANSPOSE([[zv(i,1:4)],[zv(i,1)]])
  color = s.full.cedata.eta.color(s.full.cedata.Iresistance[i])
  color = vessel_colors[s.full.cedata.Iresistance[i]]
;  color = s.full.cedata.eta.color[s.full.cedata.Iresistance[i]]
  POLYFILL,rvec,zvec,color=color,NOCLIP=0
  OPLOT,rvec,zvec,color=boxcolor
ENDFOR
RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_element,s,xmin,xmax,ymin,ymax,size,title

;!P.FONT = 0
;DEVICE,/HELVETICA,/BOLD, FONT_INDEX=5
;!p.background=255
plot,[xmin,xmax],[ymin,ymax],psym=3,xstyle=1,ystyle=1,charsize=1.3,$
	xtitle='!5R(m)',ytitle='!5Z(m)',title=title
shade_elements_resist,s,0

RETURN
END

;--------------------------------------------------------------------------------------
PRO plot_elements,file,xmin,xmax,ymin,ymax,size
title	= '!5Outline of Conducting Elements'
s0 	= read_conductor_data(file)
plot_element,s0,xmin,xmax,ymin,ymax,size,title
RETURN
END

;--------------------------------------------------------------------------------------
PRO oplot_efit_contours,shot=shot,time=time,index=index,posvec=posvec,double=double
COMMON plot_data, data
dposvec = [.13,.1,.93,.92]
IF KEYWORD_SET(posvec) THEN dposvec = posvec
plot_efit,shot=shot,time=time,index=index,data=data,posvec=dposvec,double=double,/nologo
END

;--------------------------------------------------------------------------------------
pro plot_efit,shot=shot,time=time,animate=animate,tree=tree,	$
              plotwall=plotwall,jsodata=jsodata,gfile=gfile,	$
              lstwall=lstwall,colwall=colwall,fillwall=fillwall,$
              gstruc=gstruc,index=index,nologo=nologo,		$
              beam=beam,comment=comment,double=double,save=save,$
              ps=ps,eps=eps,path=path,extension=extension,	$
              num_vfs=num_vfs,					$
              vfs_space_cm=vfs_space_cm,ptrz=ptrz,		$
              xpos=xpos,ypos=ypos,autoupdate=autoupdate,	$
              posvec=posvec,notitle=notitle,			$
              data=data,quiet=quiet,winsize=winsize,		$
              poloidal_mirnovs=poloidal_mirnovs,		$
              oplot_bnd=oplot_bnd,				$
              modify_oplot_bnd=modify_oplot_bnd,		$
              notiles=notiles,cds=cds,nohhfw=nohhfw,		$
              plot_range=plot_range,limdat=limdat,tfdat=tfdat,	$
              bndthk=bndthk,oplot_limiter=oplot_limiter,	$
              outg=outg

;--- ignore extraneous info?
dwinsize	= 800L
dwinsize_auto	= 600L

IF KEYWORD_SET(winsize)    THEN dwinsize = winsize
IF KEYWORD_SET(autoupdate) THEN dwinsize = dwinsize_auto

;--- Get tree name
dtree = 'EFIT'
IF KEYWORD_SET(tree) THEN dtree = STRUPCASE(STRTRIM(tree,2))

;--- Get index of EFIT
dindex = 1
;
IF KEYWORD_SET(index) THEN BEGIN
    dindex = index
    IF index LT 0 THEN BEGIN
        dindex = ABS(index)
        dtree = 'LRDFIT'
    ENDIF
ENDIF
index_ = strcompress(dindex,/remove_all)

;--- Event to wait for for autoupdate
efit_done_event = 'phoenixDone_P'+index_

;--- ignore extraneous info?
iquiet	= KEYWORD_SET(quiet) OR KEYWORD_SET(autoupdate)

IF KEYWORD_SET(autoupdate) THEN animate = 1L

IF KEYWORD_SET(gfile) OR KEYWORD_SET(gstruc) THEN BEGIN
    animate 	= 0L
    autoupdate 	= 0L
    GOTO, noauto
ENDIF

;--- Assume this is a plasma shot
GOTO, plasmashot

;--- If on autoupdate mode, and a vacuum shot, wait for the next shot
vacuumshot:
PRINT,'Waiting for EFIT to complete for next shot...'
mdswfevent,efit_done_event

;--- A shot or geqdsk with plasma data
plasmashot:

shotin = lastshot()
IF KEYWORD_SET(shot) THEN shotin = shot
oldshot = shotin
newshot = oldshot

;--- Setup for animation
IF KEYWORD_SET(animate) THEN pixmap=1 ELSE pixmap=0

;--- Skip automation stuff
noauto:

;--- Set x,y position of plot window
dxpos	= 730
dypos	= 478
IF KEYWORD_SET(xpos) THEN dxpos = xpos
IF KEYWORD_SET(ypos) THEN dypos = ypos

;--- Postscript file or plot?
print	= 0
IF KEYWORD_SET(ps ) THEN PRINT = 1
IF KEYWORD_SET(eps) THEN PRINT = 1
;
; Read input data (efit and Jon's vessel)
;

; Put in D.G. patented colortable
;
;

IF ((NOT KEYWORD_SET(gfile)) AND (NOT KEYWORD_SET(gstruc))) THEN BEGIN

    IF KEYWORD_SET(data) THEN BEGIN
        IF (data.shot EQ shotin) AND (data.index EQ index) THEN GOTO, usedata
    ENDIF
    
    mdsopen,dtree,shotin

    psirzt 	= mdsvalue('\'+DTREE+'0'+index_+'::PSIRZ')
    times 	= mdsunits('\'+DTREE+'0'+index_+'::PSIRZ',/time)
    zt 		= mdsvalue('\'+dtree+'0'+index_+'::z')
    rt 		= mdsvalue('\'+dtree+'0'+index_+'::r')
    ssimagt 	= mdsvalue('\'+dtree+'0'+index_+'::ssimag')
    ssibryt 	= mdsvalue('\'+dtree+'0'+index_+'::ssibry')

;    rmidint 	= mdsvalue('\'+dtree+'0'+index_+'::rmidin')
;    rmidoutt 	= mdsvalue('\'+dtree+'0'+index_+'::rmidout')

    zbdryt 	= mdsvalue('\'+dtree+'0'+index_+'::zbdry')
    rbdryt 	= mdsvalue('\'+dtree+'0'+index_+'::rbdry')
    nbdryt 	= mdsvalue('\'+dtree+'0'+index_+'::nbdry')

    zlimt 	= REFORM(mdsvalue('\'+dtree+'0'+index_+'::zlim'))
    rlimt 	= REFORM(mdsvalue('\'+dtree+'0'+index_+'::rlim'))
    nlimt 	= mdsvalue('\'+dtree+'0'+index_+'::limitr')

    wsize 	= MAX(ROUND(mdsvalue('\'+dtree+'0'+index_+'::mw')))
    hsize 	= MAX(ROUND(mdsvalue('\'+dtree+'0'+index_+'::mh')))
    codename	= dtree+'_'

    rmidint	= MIN(rbdryt,/dim)
    rmidoutt	= MAX(rbdryt,/dim)

    mdsclose,dtree,shotin

;stop

    data = {	shot:shotin,		$
                index:dindex,		$
                psirzt:psirzt,		$
                times:times, 		$
                zt:zt,			$
                rt:rt,			$
                ssimagt:ssimagt,	$
                ssibryt:ssibryt,	$
                rmidint:rmidint,	$
                rmidoutt:rmidoutt,	$
                zbdryt:zbdryt,		$
                rbdryt:rbdryt,		$
                nbdryt:nbdryt,		$
                zlimt:zlimt,		$
                rlimt:rlimt,		$
                nlimt:nlimt,		$
                wsize:wsize,		$
                hsize:hsize		}

    usedata:

                                ;--- Use only time-slices with closed flux
    iplasma 	= WHERE(data.nbdryt GT 0)
    stext 	= STRTRIM(STRING(data.shot),2)
    ntext 	= STRTRIM(STRING(N_ELEMENTS(iplasma)),2)
    IF (iplasma[0] NE -1) THEN BEGIN
        PRINT,'There are '+ntext+' time-slices with closed boundaries for shot '+stext
    ENDIF
    IF (iplasma[0] EQ -1) THEN BEGIN
        PRINT,'There are no EFIT time-slices with closed boundaries for shot '+stext
        IF KEYWORD_SET(animate) AND KEYWORD_SET(autoupdate) THEN GOTO, vacuumshot
        GOTO, finished
    ENDIF
    
    k          = iplasma
    psirzt 	= data.psirzt[*,*,k]
    times 	= data.times[k]
    zt 	= data.zt[*,k]
    rt 	= data.rt[*,k]
    ssimagt 	= data.ssimagt[k]
    ssibryt 	= data.ssibryt[k]
    rmidint 	= data.rmidint[k]
    rmidoutt 	= data.rmidoutt[k]
    zbdryt 	= data.zbdryt[*,k]
    rbdryt 	= data.rbdryt[*,k]
    nbdryt 	= data.nbdryt[k]
    zlimt 	= data.zlimt[*,k]
    rlimt 	= data.rlimt[*,k]
    nlimt 	= data.nlimt[k]
    wsize	= data.wsize
    hsize	= data.hsize

    dtime = [MIN(times),MAX(times)]
    IF KEYWORD_SET(time) THEN dtime = time
    time = dtime

ENDIF

IF KEYWORD_SET(gfile) OR KEYWORD_SET(gstruc) THEN BEGIN
    IF KEYWORD_SET(gfile)  THEN g = readg(gfile)
    IF KEYWORD_SET(gstruc) THEN g = gstruc

    outg	= g

    wsize	= g.mw
    hsize	= g.mh
    psirzt	= dblarr(wsize,hsize,1)
    psirzt[*,*,0]	= g.psirz 
    rt 		= g.r[0:wsize-1]
    zt 		= g.z[0:hsize-1]
    ssimagt	= g.ssimag
    ssibryt 	= g.ssibry

    ;-------------------------------------
    nbnd 	= g.nbdry
    rbnd	= g.bdry[0,0:nbnd-1]
    zbnd	= g.bdry[1,0:nbnd-1]

    rmidin	= MIN(rbnd)
    rmidout	= MAX(rbnd)

    rmidint 	= rmidin
    rmidoutt 	= rmidout

    rbdryt	= DBLARR(nbnd,1)
    zbdryt	= DBLARR(nbnd,1)
    rbdryt[*,0] 	= rbnd
    zbdryt[*,0] 	= zbnd
    nbdryt 	= g.nbdry

    ;-------------------------------------
    nlim 	= g.limitr
    rlim	= g.lim[0,0:nlim-1]
    zlim	= g.lim[1,0:nlim-1]

    rlimt	= DBLARR(nlim,1)
    zlimt	= DBLARR(nlim,1)
    rlimt[*,0] 	= rlim
    zlimt[*,0] 	= zlim
    nlimt 	= g.limitr
    ;-------------------------------------

    time	= g.time/1000.0
    times 	= time[0]
    shotin	= g.shot
    cpasma	= g.cpasma
    ecase	= g.ecase
    codename	= STRTRIM(ecase[0]+STRMID(ecase[1],0,2),2)+'_'
ENDIF

;--- Get the data for the vessel     
IF KEYWORD_SET(cds) THEN s0 = cds
IF NOT KEYWORD_SET(cds) THEN s0	= nstx_conductor_data_structure(shotin)

ntime = n_elements(time)
if ntime eq 1 then begin
    tind = max(where(times le time[0]))
endif else begin
    tind = WHERE(times GE time[0] AND times LE time[1])
endelse

;--- Use these times during plots
;
tuse	= times[tind]
ntimes = N_ELEMENTS(tind)
ngood 	= ntimes
win 	= intarr(ntimes)

fmt 	= '("Plotting time range:  t(min) = ",f8.4,"s, t(max) = ",f8.4,"s")'  	
PRINT,MIN(tuse),MAX(tuse),format=fmt

FOR i=0,ntimes-1 DO BEGIN

   loadct,2,file='/u/dgates/public/mine.tbl'
   red_ = [0,159,186,214,241,255,255,255,255,0,135,255,255,217,220,190]
   green_ = [0,0,17,52,88,124,160,196,232,72,255,174,255,217,12,190]
   blue_ = [0,0,0,0,0,0,58,133,207,255,0,0,255,217,103,190]
   tvlct,red_,green_,blue_

   black = 0
   red   = 1
   orange= 6
   blue	 = 9
   white = 12
   
   !p.background = white
   
   !p.color = black
   
   psirz 	= psirzt[*,*,tind[i]]
   ssibry  	= ssibryt[tind[i]]
   ssimag  	= ssimagt[tind[i]]
   rmidin  	= rmidint[tind[i]]
   rmidout 	= rmidoutt[tind[i]]
   ztop 	= max(zbdryt[*,tind[i]])
   zbot 	= min(zbdryt[*,tind[i]])
   rbdry 	= rbdryt[0:nbdryt[tind[i]]-1,tind[i]]
   zbdry 	= zbdryt[0:nbdryt[tind[i]]-1,tind[i]]

   rlim 	= rlimt[0:nlimt[tind[i]]-1,tind[i]]
   zlim 	= zlimt[0:nlimt[tind[i]]-1,tind[i]]
   
   r 	= rt[*,tind[i]] & z = zt[*,tind[i]]
   rind 	= WHERE(r gt .18)
   zind 	= WHERE(z gt -1.5 and z lt 1.5)
   rimax 	= MAX(rind) 
   rimin 	= MIN(rind) 
   zimax 	= MAX(zind)
   zimin 	= MIN(zind)

;
; Print or plot choices
;

    IF KEYWORD_SET(posvec) THEN GOTO, skip_windows

    if keyword_set(print) then begin
        if i eq 0 then begin
            dpath = '.'
            IF KEYWORD_SET(path) THEN dpath = path
            IF NOT KEYWORD_SET(eps) THEN deps = 0
            IF     KEYWORD_SET(eps) THEN deps = 1
            extn = '.ps'
            IF deps THEN extn = '.eps'
            IF KEYWORD_SET(extension) THEN extn = '_'+extension+extn
            set_plot,'ps'
            filename= dpath+'/'+codename+strcompress(shotin,/remove_all)+'_'+ $
              strcompress(fix(tuse[i]*1000),/remove_all)+extn
            if not keyword_set(double) then begin
                device,filename=filename,/color,bits_per_pixel=8,	$
                  xsize=14,ysize=28*.833*1.05, 			$
                  xoff=3,yoff=2,encapsulated=deps
            endif else begin
                device,filename=filename,/color,bits_per_pixel=8,	$
                  xsize=26,ysize=26*.833*1.05, 			$
                  xoff=-3,yoff=2,encapsulated=deps
            endelse
        endif
    endif else begin

;        shades(where(shades eq !d.n_colors)) = 0

        if keyword_set(animate) then begin
            if not keyword_set(double) then begin
                window,/free,xsize=dwinsize/2,ysize=dwinsize,/pixmap,xpos=dxpos,ypos=dypos
                win[i] = !d.window
            endif else begin
                window,/free,xsize=dwinsize,ysize=dwinsize,/pixmap,xpos=dxpos,ypos=dypos
                win[i] = !d.window
            endelse

            fmt = '("Finished window = ",i6," for t=",f8.4,"s")'  	
            IF NOT iquiet THEN PRINT,win[i],tuse[i],format=fmt
;;            PRINT,win[i],tuse[i],format=fmt
            
        endif else begin
            print,'Defining window ',i
            if not keyword_set(double) then $
              window,i,xsize=dwinsize/2,ysize=dwinsize,xpos=dxpos,ypos=dypos $
            else $
              window,i,xsize=dwinsize,ysize=dwinsize,xpos=dxpos,ypos=dypos
        endelse
    endelse
;
; Put in my patented colortable
;

skip_windows:

;--- set the shades	
    shades 	= LONG(!d.n_colors*float(bytscl(psirz(rimin:rimax, zimin:zimax)<0))/255.) > 1
    IF KEYWORD_SET(print) THEN shades = LONG(smooth(bytscl(psirz(rimin:rimax,zimin:zimax)<0),2)) > 1

;
; Make shaded surface
;
    dplot_range = 2.4

    IF KEYWORD_SET(plot_range) THEN dplot_range = plot_range

    if not keyword_set(double) then begin
        xmin	= -dplot_range * 0.05
;        xmin	= -0
        xmax	= +dplot_range
        ymin	= -dplot_range
        ymax	= +dplot_range
        size	= 0.7
    endif else begin
        xmin	= -dplot_range
        xmax	=  dplot_range
        ymin	= -dplot_range
        ymax	= +dplot_range
        size	= 0.7
    endelse
    
    if keyword_set(print) then begin
        pos_vec = [.2,.1,.9,.94]
        ch_siz = 1.4
        th=2
    endif else begin
        ch_siz=2
        pos_vec = [.13,.1,.93,.92]
        th=4
    endelse
    
;--- Set the position for plotting
    IF KEYWORD_SET(posvec) THEN pos_vec = posvec
    
;--- Shade the surfaces and show the title

    title_text 	= ''+DTREE+'0'+index_+' shot='+strcompress(shotin,/remove_all)+', time=' + $
      strcompress(fix(1000*tuse[i]),/remove_all)+'ms'
    IF KEYWORD_SET(notitle) THEN title_text = ''
    
    title_text1	= ''
    IF KEYWORD_SET(posvec)  THEN title_text1 = title_text

    shade_surf,psirz(MIN(rind):MAX(rind),zind) < 0,r(rind),z(zind),			$
      shades= shades,ax=90,az=0,xr=[xmin,xmax],yr=[ymin,ymax],xstyle=1,ystyle=1, 	$
      charsize=ch_siz,xtitle='!5R(m)',ytitle='!5Z(m)',zs=4,clip= 			$
      [xmin,ymin,xmax,ymax],position=pos_vec,/nod,font=0,back=white,color=black, 	$
      thick=th,title=title_text1
;
    polyfill,rbdry,zbdry,color=1

    if keyword_set(double) then polyfill,-1.0*rbdry,zbdry,color=1
    colors = [2+indgen(8),white+intarr(11)]

;--- Set spacing between vacuum flux surfaces
;
    inlevs  	= 8     ; number of flux surfaces inside plasma
    xnlevs	= 6     ; number of flux surfaces external to plasma
    dxcm	= 1.0 	; spatial separation between flux surfaces in cm

    IF KEYWORD_SET(vfs_space_cm) THEN dxcm   = DOUBLE(vfs_space_cm)
    IF KEYWORD_SET(num_vfs) 	 THEN xnlevs = num_vfs

    deltpsi	= (ssimag-ssibry) / DOUBLE(inlevs-1) 
    levels	= ssibry + deltpsi * DINDGEN(inlevs)

    rbnd_rmax = MAX(rbdry)
    ibnd_rmax = (WHERE(rbdry EQ rbnd_rmax))[0]
    zbnd_rmax = zbdry[ibnd_rmax]
    zdiff	  = ABS(z-REPLICATE(zbnd_rmax,N_ELEMENTS(z)))
    izrbmax   = WHERE(zdiff EQ MIN(zdiff))
    psirzline = psirz[0:wsize-1,izrbmax]
    dpsidr 	  = DERIV(r,psirzline)
    dpsidrob  = SPLINE(r,dpsidr,rbnd_rmax)
    deltpsicm = -deltpsi/dpsidrob*100.
    rbnd_txt  = STRTRIM(STRING(rbnd_rmax),2)+'m'
    delp_txt  = STRTRIM(STRING(dxcm),2)+'cm'
    xlevfac	  = dxcm / deltpsicm * deltpsi
    xlevels	  = ssibry - xlevfac * DINDGEN(xnlevs)

    IF NOT iquiet THEN PRINT, 'Vacuum flux surface separation at R(max) = '+ $
      rbnd_txt+' on boundary = '+delp_txt

; Calculate contour bounds

    rmin_ind = max(where(r le 0.98*rmidin))
    if rmin_ind EQ -1 then rmin_ind = 0
    rmax_ind = min(where(r ge 1.02*rmidout))
    zmin_ind = max(where(z le 1.02*zbot))
    zmax_ind =  min(where(z ge 1.02*ztop))

; Calculate contour indices

    IF (zmax_ind EQ zmin_ind) AND (zmax_ind NE 0) THEN zmin_ind = zmax_ind - 1L
    IF (rmax_ind EQ rmin_ind) AND (rmax_ind NE 0) THEN rmin_ind = rmax_ind - 1L

; Plot the contour

    contour,psirz(rmin_ind:rmax_ind,zmin_ind:zmax_ind),r(rmin_ind:rmax_ind),  $
      z(zmin_ind:zmax_ind), levels = (levels(sort(levels))),/over,path_xy=xy, $
      path_info=info,/path_data_coords

skip_contour:
    keep=-1
    
;
; Throw out extra contours
;
    for k=0,n_elements(info)-1 do begin
        s = indgen(info(k).n)
        xgood = where(xy(0,info(k).offset+s) ge $
                      1.01*min(rbdry) and xy(0,info(k).offset+s) le 0.99*max(rbdry))
        zgood = where(xy(1,info(k).offset+s) ge $
                      0.999*min(zbdry) and xy(1,info(k).offset+s) le 0.999*max(zbdry))
        if n_elements(zgood) eq n_elements(xy(1,info(k).offset+s)) and $
          n_elements(xgood) eq n_elements(xy(0,info(k).offset+s)) then begin
            if n_elements(keep) eq 1 then begin
                if keep eq -1 then keep=k else keep=[keep,k]
            endif else begin
                keep=[keep,k]
            endelse
        endif
    endfor
    
;
; Make shaded plasma plots
;
;    
;IF NOT KEYWORD_SET(jsodata) THEN BEGIN
    
    IF keep[0] EQ -1 THEN GOTO, skip_polyfill
    
    for k = 0,n_elements(keep)-1 do begin
        s = [indgen(info(keep(k)).n),0]
        polyfill,xy(0,info(keep(k)).offset+s),xy(1,info(keep(k)).offset+s), $
          color=k+2
        if keyword_set(double) then begin
            polyfill,-1.*xy(0,info(keep(k)).offset+s),xy(1,info(keep(k)).offset+s), $
              color=k+2
        endif
    endfor

skip_polyfill:

;ENDIF

; Plot limiter surface from equilibrium data?
;
    IF KEYWORD_SET(oplot_limiter) THEN BEGIN
        OPLOT,rlim,zlim,COLOR=blue,THICK=6
        IF KEYWORD_SET(double) THEN OPLOT,-rlim,zlim,COLOR=blue,THICK=6
    ENDIF

; Plot a user-specified limiter surface?
;
    IF KEYWORD_SET(limdat) THEN BEGIN
        rlim	= limdat.rlim
        zlim	= limdat.zlim
        OPLOT,rlim,zlim,COLOR=blue,THICK=6
        IF KEYWORD_SET(double) THEN OPLOT,-rlim,zlim,COLOR=blue,THICK=6
    ENDIF

; Plot a user-specified TF coil?
;
    IF KEYWORD_SET(tfdat) THEN BEGIN

        tfrout	= tfdat.rout
        tfzout	= tfdat.zout
        tfrin	= tfdat.rin
        tfzin	= tfdat.zin

        polyfill,[tfrout,REVERSE(tfrin)],$
          [tfzout,REVERSE(tfzin)],color=orange
        OPLOT,tfrout,tfzout,color=black,THICK=1
        OPLOT,tfrin ,tfzin ,color=black,THICK=1

    ENDIF
;
; Make the vessel a nice color
;
    if not keyword_set(notiles) then begin
        plot_nstx_tiles
        if keyword_set(double) then plot_nstx_tiles,/reverse
    endif

    contour,psirz[0:wsize-1,0:hsize-1],r,z,levels = levels(sort(levels)),/over
    contour,psirz[0:wsize-1,0:hsize-1],r,z,levels = xlevels(sort(xlevels)),/over

    ;--- plot boundary flux contour red
    dbndthk	= 1
    IF KEYWORD_SET(bndthk) THEN dbndthk = bndthk
    contour,psirz[0:wsize-1,0:hsize-1],r,z,levels = xlevels[0],/over,thick=dbndthk,color=red

    rarr	= dblarr(wsize,hsize)
    zarr	= dblarr(wsize,hsize)
    for k=0,hsize-1 do rarr[*,k] = r
    for k=0,wsize-1 do zarr[k,*] = z

    IF KEYWORD_SET(ptrz) THEN BEGIN
        plevels	= DINDGEN(inlevs)/DOUBLE(inlevs) * MAX(ptrz)
        print,plevels
        contour,ptrz[0:wsize-1,0:hsize-1],rarr,zarr,nlevels = 51,/overplot,color=1
        help,ptrz
        oplot,[0,1],[0,1]
    ENDIF

    if keyword_set(double) then BEGIN
        contour,psirz,-1*r,z, levels = levels(sort(levels)),/over
        contour,psirz,-1*r,z, levels = xlevels(sort(xlevels)),/over
        contour,psirz,-1*r,z, levels = xlevels[0],/over,thick=2
    ENDIF

    IF NOT KEYWORD_SET(nohhfw) THEN plot_nstx_hhfw_antenna

    shade_elements_resist,s0,0
    if keyword_set(double) then begin
        if i eq 0 then begin
            s1 = s0
            s1.full.cedata.r = -1.0*s1.full.cedata.r
        endif
        shade_elements_resist,s1,0
    endif

;
; Info
;

    IF NOT KEYWORD_SET(posvec) THEN BEGIN
        if not keyword_set(double) then begin
            xyouts,0.50,.95,title_text,/norm,charsize=ch_siz,font=1,alignment=0.5
        endif else begin
            xyouts,0.50,.97,title_text,/norm,charsize=ch_siz,font=1,alignment=0.5
        endelse
    ENDIF

    if keyword_set(comment) then xyouts,.25,1,'!17'+comment,Charsize=3.0,/norm

;
; Logo and Header 
;

    if keyword_set(beam) then begin
        color1 = 0 & color2 = white
        oplot,[-.05,.05,.05,-.05,-.05]+.5,[-.2,-.2,.2,.2,-.2],thick=4, $
          color=color1
        oplot,[-.05,.05,.05,-.05,-.05]+.6,[-.2,-.2,.2,.2,-.2],thick=4, $
          color=color1
        oplot,[-.05,.05,.05,-.05,-.05]+.7,[-.2,-.2,.2,.2,-.2],thick=4, $
          color=color1
        polyfill,[-.04,.04,.04,-.04,-.04]+.5,[-.19,-.19,.19,.19,-.19],thick=4,$
          color=color2
        polyfill,[-.04,.04,.04,-.04,-.04]+.6,[-.19,-.19,.19,.19,-.19],thick=4,$
          color=color2
        polyfill,[-.04,.04,.04,-.04,-.04]+.7,[-.19,-.19,.19,.19,-.19],thick=4,$
          color=color2
    endif

    IF KEYWORD_SET(plotwall) then begin
        v = name2vec(plotwall)
        loadct,39
        dwlst = 0
        dwcol = 250
        IF KEYWORD_SET(lstwall) THEN dwlst = lstwall
        IF KEYWORD_SET(colwall) THEN dwcol = colwall
        print,lstwall,colwall
        print,dwlst,dwcol
        plotwall,6,v,linestyle=dwlst,color=dwcol,fill=fillwall
    ENDIF

    IF KEYWORD_SET(oplot_bnd) then begin
        FOR bindx=0,oplot_bnd.nwbps-1 DO BEGIN
            opbr = oplot_bnd.rrbps[*,bindx]
            opbz = oplot_bnd.rzbps[*,bindx]
            IF KEYWORD_SET(modify_oplot_bnd) THEN BEGIN
                opbr = opbr + modify_oplot_bnd.dr
                opbz = opbz + modify_oplot_bnd.dz
                opbr = opbr * modify_oplot_bnd.scale
                opbz = opbz * modify_oplot_bnd.scale
            ENDIF
            oplot,opbr,opbz,thick=15,color=09;+bindx
        ENDFOR
    ENDIF

    IF KEYWORD_SET(jsodata) THEN BEGIN
        x_j	= jsodata.x
        z_j	= jsodata.z
        nt	= n_elements(x_j[*,0])
        nr	= n_elements(x_j[0,*])
        ndivs	= 16
        nsurf	= ndivs+1
        ind	= (LINDGEN(ndivs)+1)*(nr/ndivs)
        xjs	= x_j[*,ind]
        zjs	= z_j[*,ind]
        press	= DBLARR(nt,ndivs)
        FOR i=0,nt-1 do press[i,*] = jsodata.p[ind]
        press = press/max(press)
        cscale	= 0.4
        col0	= 256.*DOUBLE(ndivs-1)/DOUBLE(ndivs)
        cols	= (DINDGEN(ndivs)/DOUBLE(ndivs-1)*cscale+(1.0-cscale))*col0
        levs	= REFORM(REVERSE(press[0,*],2))

        loadct,3
        contour,press,xjs,zjs,levels=levs,/fill,c_colors=cols,/overplot
        polyfill,xjs[*,0],zjs[*,0],color=col0
        for i=0,ndivs-1 do oplot,xjs[*,i],zjs[*,i],color=0
        oplot,x_j[*,nr-1],z_j[*,nr-1],color=0

    ENDIF

    if not keyword_set(nologo) then begin
        if keyword_set(print) then begin
            a = read_bmp('/p/nstxusr1/util/fas_cvs/hires_logo.bmp',red,green,blue)
            a(where(a eq 255)) = 0
            tvlct,red,green,blue
            if not keyword_set(double) then $
              tv,a,10,0,xsize=5,ysize=1,/centimeters $
            else $
              tv,a,14,-1,xsize=5,ysize=1,/centimeters
        endif else begin
            a = read_bmp('/p/nstxusr1/util/fas_cvs/lores_logo.bmp',red,green,blue)
            a(where(a eq 216)) = 100
            a(where(a eq 192)) = 45
            if not keyword_set(double) then $
              tv,a,300,10 $
            else $
              tv,a,500,0
        endelse
    endif

    IF KEYWORD_SET(poloidal_mirnovs) THEN plot_nstx_poloidal_mirnovs
    
endfor


;
; Clean up
;
if keyword_set(print) then begin
    device,/close
    set_plot,'x'
    IF NOT deps THEN spawn,'ps2pdf '+filename
endif

if keyword_set(animate) then begin
    done=0
    if not keyword_set(double) then $
      xsize=dwinsize/2 $
    else $
      xsize=dwinsize
    window,0,xsize=xsize,ysize=dwinsize,xpos=dxpos,ypos=dypos

    stext = STRTRIM(STRING(newshot),2)
    PRINT,'---> Plotting animation sequence for shot '+stext+' <---'

    WHILE NOT done DO BEGIN     ;    * Animate *
        FOR i = 0,ngood-1 DO BEGIN
            DEVICE,copy=[0,0,xsize,dwinsize,0,0,win[i]]
            WAIT,0.3

                                ; Check for new shot number every mgk plots
            mgk = 10L
            ngk = i MOD mgk
            IF KEYWORD_SET(autoupdate) AND (ngk EQ 0) THEN BEGIN
                newshot = lastshot()
                IF newshot NE oldshot THEN BEGIN

                    FOR i=0,ngood-1 DO BEGIN
;;;                        PRINT,'Cleaning up video memory for window ',win[i]
                        WDELETE, win[i] ; * When done, free video memory *
                    ENDFOR

                    stext = STRTRIM(STRING(newshot),2)
                    PRINT,'------------------------------------------------------'
                    PRINT,'New shot number detected = '+stext
                    PRINT,'Waiting for EFIT to complete for shot '+stext+'...'
                    mdswfevent,efit_done_event
                    PRINT,'EFIT autocode completed for shot '+stext
                    GOTO, plasmashot						
                ENDIF
            ENDIF

        ENDFOR
        WAIT,3.0

; 	res = get_kbrd(0)
;	if res ne '' then done = 1

    ENDWHILE

    if keyword_set(save) then begin
        for i=0,ngood-1 do begin
;   wset,win(i)
;   a = tvrd()
;   a2 = rebin(a,xsize/2,dwinsize/2)
            if i eq 0 then mpegID = MPEG_OPEN([xsize,dwinsize],file=codename+ $
                                              strcompress(shotin,/remove_all)+'_'+ $
                                              strcompress(fix(tuse[i]*1000),/remove_all)+'.mpg')
            print,'Now saving frame number: ',strcompress(i+1),' of ', $
              strcompress(ngood),' frames'
            MPEG_PUT,mpegID,window=win(i),/COLOR,/order,frame=i
        endfor
        print,'Now writing file :'+codename+ $
          strcompress(shotin,/remove_all)+'_'+ $
          strcompress(fix(tuse[0]*1000),/remove_all)+'.mpg'
        MPEG_SAVE,mpegID,filename=codename+strcompress(shotin,/remove_all)+'_'+ $
          strcompress(fix(tuse[0]*1000),/remove_all)+'.mpg'
        MPEG_CLOSE,mpegID
    endif

    for i=0,ngood-1 do BEGIN
        print,'Cleaning up video memory for window ',win[i]
        wdelete,win[i]          ; * When done, free video memory *
    endfor

endif
;
; Done
;
finished:
;
end