;+
;  NAME: spherical_draw_field
; 
;  PURPOSE: 
;    Given vector fieldline data on a spherical grid, this procedure renders
;    these fieldlines using an orthographic projection.  The result is a
;    true-color image.
; 
;  CALLING SEQUENCE: 
;    spherical_draw_field,sph_data,xsize=xsize,ysize=ysize,bcent=bcent,
;      lcent=lcent,width=width,im_data=im_data,imsc=imsc,thick=thick,
;      outim=outim,/onscreen,/movie,/drawopen,/drawclosed,/nolines,/noimage,
;      /for_ps,/indexed_colors,/quiet,extra_objects=extra_objects
; 
;  INPUTS: 
;    sph_data = a structure of type spherical_field_data (see
;               spherical_field_data__define.pro) with the following fields
;               defined on input: ptr,ptth,ptph,nstep.  Basically, one needs
;               the trajectories of each fieldline in (r,theta,phi)-space.
;    xsize = number of pixels in the width of the rendered image (default=512)
;    ysize = number of pixels in the height of the rendered image (default=512)
;    bcent = central latitude (in degrees) of projection centroid (default=0)
;    lcent = central longitude (in degrees) of projection centroid (default=0)
;    roll = roll angle of projection counterclockwise from north (default=0)
;    width = width of horizontal extent of final image relative to the
;            diameter of the outer radial boundary (default=1)
;    im_data = a structure of type spherical_image_data (see
;                 spherical_image_data__define.pro), with all fields defined,
;                 to be texture-mapped, as a grayscale image, onto an opaque
;                 sphere located at the specified radius, assuming /noimage is
;                 not set (default = radial component of the vector field data
;                 at inner radius)
;    imsc = data value(s) to which to scale the image contained in the
;           im_data structure; if a scalar, then image is scaled to
;           [-imsc,imsc], and if a two-element vector then image is scaled to
;           [imsc(0),imsc(1)]  (default = image scale, centered around 0)
;    thick = thickness of field lines
;    onscreen = if set, then display image(s) onscreen
;    movie = if set, creates movie sequence of field-line data, with each
;            image rotated 1 degree from the last
;    drawopen = if set, only open field lines are drawn
;    drawclosed = if set, only closed field lines are drawn
;    nolines = if set, no field lines are drawn, supersedes both the drawopen
;              and drawclosed keywords
;    noimage = if set, no central image is inserted and the central region is
;              left transparent
;    for_ps = if set, interchanges white and black colors (for .ps printouts)
;    indexed_colors = set this flag to return a 8-bit color image (default is
;                    to return a true color image)
;    quiet = set to inhibit screen output
;    extra_objects = extra objects to be added to the views
; 
;  OUTPUTS:
;    outim = the rendered image in true-color format (unless /indexed_colors
;            flag is set)
;
;  NOTES:
;    1.  Ever confusing color issues: In reality, only 256 colors are used.
;        In this code, however, 24-bit colors are used in the object graphics
;        routines, and then converted at the end to 8-bit color images in the
;        end (if the /indexed_colors flag is set).
;    2.  Occasionally, the draw method from the IDLgrBuffer destination object
;        sometimes generates underflow, overflow, and illegal operand errors
;        (at least in IDL v.6.1), but these tend to happen only
;        intermittently, and I don't see any problems with the output, so I
;        usually ignore the errors.  I believe these are exceptions returned
;        from the renderer, and IDL passes them on.
;    3.  If both /movie and /onscreen are set, then routine only displays the
;        first frame of the movie, but will return full movie to outim.
;     
;  MODIFICATION HISTORY:
;    M.DeRosa - 11 Jan 2006 - created
;               12 Jan 2006 - fixed clipping problem when width<1
;               12 Jan 2006 - also added im_data keyword
;               18 Jan 2006 - fixed problem with non-square [xsize,ysize]
;               27 Jan 2006 - takes into account invalid lines
;               26 Mar 2007 - fixed minor bug in cleanup at end
;               26 Mar 2007 - now deals with background images in im_data 
;                             that do not span 360 degrees of longitude 
;               21 Apr 2010 - added extra_objects keyword
;               25 Apr 2011 - added an ambient and a fixed light (at the
;                             observer's position) so that surfaces
;                             provided through the extra_objects keyword are
;                             shaded properly
;               25 Apr 2011 - now accommodates the extra_objects tag in the
;                             sph_data structure
;               18 May 2011 - added roll keyword
;               22 Nov 2011 - if im_data set, now renders image at radius
;                             stored in the im_data structure
;               22 Nov 2011 - added outlines for bounded images (for full
;                             spherical images, only the latitudinal
;                             boundaries of the lower boundary image are
;                             outlined)
;               25 Jul 2013 - fixed bug in texture mapping for nonuniform grids
;
;-

pro spherical_draw_field,sph_data,xsize=xsize,ysize=ysize,bcent=bcent, $
  lcent=lcent,roll=roll,width=width,im_data=im_data,imsc=imsc, $
  thick=thick,outim=outim,onscreen=onscreen,movie=movie,nolines=nolines, $
  noimage=noimage,drawopen=drawopen,drawclosed=drawclosed,for_ps=for_ps, $
  indexed_colors=indexed_colors,quiet=quiet,extra_objects=extra_objects

;  usage message
if n_elements(sph_data) eq 0 then begin
  print,'  ERROR in spherical_draw_field: no input data provided'
  print,'  Calling sequence: spherical_draw_field,sph_data,xsize=xsize,'+ $
    'ysize=ysize,bcent=bcent,lcent=lcent,roll=roll,width=width,'+ $
    'im_data=im_data,'+ $
    'imsc=imsc,thick=thick,outim=outim,/onscreen,/movie,/drawopen,'+ $
    '/drawclosed,/nolines,/noimage,/for_ps,/indexed_colors,/quiet,'+ $
    'extra_objects=extra_objects'
  return
endif

;  set defaults
if n_elements(xsize) eq 0 then xsize=512 else xsize=round(xsize(0))
if n_elements(ysize) eq 0 then ysize=512 else ysize=round(ysize(0))
if n_elements(bcent) eq 0 then bc=0.0 else bc=bcent(0)
if n_elements(lcent) eq 0 then lc=0.0 else lc=lcent(0)
if n_elements(width) eq 0 then width=1 else width=width(0)
if n_elements(thick) eq 0 then thick=1 else thick=thick(0)
if n_elements(roll) eq 0 then rang=0.0 else rang=roll(0)

if (n_elements(drawclosed) eq 0) and (n_elements(drawopen) eq 0) then begin
  drawclosed=1b
  drawopen=1b
endif
if n_elements(drawclosed) ne 0 $
  then drawclosed=drawclosed gt 0 else drawclosed=0b
if n_elements(drawopen) ne 0 $
  then drawopen=drawopen gt 0 else drawopen=0b

;  set RGB triplets of line colors
bla=byte([  0,  0,  0])  ;  0 in palette
yel=byte([255,255,  0])  ;  not in palette
gre=byte([  0,255,  0])  ;  250/251 in palette
red=byte([255,  0,255])  ;  252/253 in palette
whi=byte([255,255,255])  ;  254/255 in palette

;  get ranges and scaling of x,y,z
rmax=max(*sph_data.rix,min=rmin)
thmin=min(*sph_data.theta,max=thmax)
phmin=min(*sph_data.phi,max=phmax)
xra=[-rmax,rmax]*width
yra=[-rmax,rmax]*width
zra=[-rmax,rmax]*(1>width)  ;  prevents clipping of items having z>width
zs=(ys=(xs=[0.,1.]))  ;  no scaling

;  create the view object
aspect=float(ysize)/float(xsize)
viewrect=[-width,-width*aspect,2*width,2*width*aspect]*rmax
if keyword_set(for_ps) then backcol=whi else backcol=bla
oview=obj_new('IDLgrView',color=backcol,viewplane=viewrect,proj=1, $
  dim=[xsize,ysize],zclip=[0.1,zra(0)-zra(1)-0.1],eye=2*rmax+1)

;  create the outer (top) model object
omodeltop=obj_new('IDLgrModel')
oview->add,omodeltop  ;  add the model to the view

;  move top model so that the front coincides with the viewplane at z=0
omodeltop->translate,0,0,-rmax*(1>width)  ;  prevents clipping of z>width

;  add an ambient light and a directional light
olight1=obj_new('IDLgrLight',type=1,loc=[0,0,2*rmax+1],intensity=0.8)
olight2=obj_new('IDLgrLight',type=0,intensity=0.5)
omodeltop->add,olight1
omodeltop->add,olight2

;  put an inner model inside the outer model (keeps the light with the viewer)
omodelin=obj_new('IDLgrModel')
omodeltop->add,omodelin

;  tilt the inner model so that (lon,lat)=(lc,bc) is at image center
omodelin->rotate,[1,0,0],-90  ;  point rotation axis up
omodelin->rotate,[0,1,0],-90-lc  ;  get correct lc
omodelin->rotate,[1,0,0],bc  ;  get correct bc
omodelin->rotate,[0,0,1],rang  ;  get correct roll angle

;  put another model inside the inner model (for moviemaking below)
omodelmov=obj_new('IDLgrModel')
omodelin->add,omodelmov  ;  add the new model to the inner model

;  create the texture map for the spherical surface at the lower radius
if not keyword_set(noimage) then begin

  ;  create the spherical_image_data structure
  if n_elements(im_data) eq 0 then $
    im_data=spherical_image_create((*sph_data.br)(*,*,0),$
      *sph_data.lon,*sph_data.lat)

  ;  check to see if image completely spans all longitudes
  dphi=min(abs((*im_data.phi)(1:2)-(*im_data.phi)(0:1)))  ;  get grid spacing
  if round((2*!dpi)/dphi) ne im_data.nlon then imbound=1 else imbound=0

  ;  replicate last longitude so that no seam in the texture map appears
  imd2=im_data
  if imbound eq 0 then begin
    imd2.image=ptr_new([*im_data.image,(*im_data.image)(0,*)])
    imd2.lon=ptr_new([*im_data.lon,360+(*im_data.lon)(0)])
    imd2.phi=ptr_new([*im_data.phi,2*!dpi+(*im_data.phi)(0)])
    imd2.nlon=im_data.nlon+1
  endif

  ;  get byte image to texture map
  if n_elements(imsc) eq 0 then $
    imsc=max([-min(*imd2.image),max(*imd2.image)])
  scim,*imd2.image,outim=byteimage,sc=imsc,/nowin,/quiet,top=249
  otexmap=obj_new('IDLgrImage',data=byteimage,palette=imd2.palette)

  ;  create vertex list for polygonal spherical surface
  phgrid=(*imd2.phi)#replicate(1,imd2.nlat)
  thgrid=!dpi/2-replicate(1,imd2.nlon)#(*imd2.theta)
  radgrid=replicate(imd2.rad,imd2.nlon,imd2.nlat)
  sph_vert=transpose([[[phgrid]],[[thgrid]],[[radgrid]]],[2,0,1])
  sph_vert=reform(sph_vert,3,imd2.nlon*imd2.nlat,/overwrite)
  rect_vert=cv_coord(from_sph=sph_vert,/to_rect)

  ;  create the connectivity array
  connect=lonarr(5,imd2.nlon-1,imd2.nlat-1)
  for i=0,(imd2.nlat-1)-1 do for j=0,(imd2.nlon-1)-1 do begin
    stpt=i*imd2.nlon+j  ;  index of starting point in vertex list
    connect(*,j,i)=[4,stpt,stpt+1,stpt+1+imd2.nlon,stpt+imd2.nlon]
  endfor

  ;  determine the texture map coords
  xg=linrange(imd2.nlon,0,1)#replicate(1,imd2.nlat)
  yg=replicate(1,imd2.nlon)#linrange(imd2.nlat,0,1)
  texture_coord=reform(transpose([[[xg]],[[yg]]],[2,0,1]),2,imd2.nlon*imd2.nlat)
  ;; ;  old stuff (wrong) below
  ;; if imbound then begin
  ;;   phirange=(*imd2.phi)(imd2.nlon-1)-(*imd2.phi)(0)
  ;;   if phirange lt 0 then phirange=phirange+2*!dpi
  ;;   texture_coord=[(sph_vert(0,*)-(*imd2.phi)(0))/phirange, $
  ;;     (sph_vert(1,*)-min(thgrid))/(max(thgrid)-min(thgrid))]
  ;; endif else begin
  ;;   texture_coord=[(sph_vert(0,*)-min(sph_vert(0,*)))/(2*!dpi), $
  ;;     (sph_vert(1,*)-min(thgrid))/(max(thgrid)-min(thgrid))]
  ;; endelse

  ;  add the texture mapped image to the model
  osphere=obj_new('IDLgrPolygon',data=rect_vert,polygons=connect,col=whi,$
    style=2,texture_map=otexmap,texture_coord=texture_coord, $
    /texture_interp,/reject)
  omodelmov->add,osphere

  ;  create outlines for image border
  if imbound then begin
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph1,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(imd2.rad,sph_data.nlat)]]))
    blobj=obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph2,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(imd2.rad,sph_data.nlat)]]))
    blobj=[blobj,obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
    omodelmov->add,blobj
  endif

  ;  create outlines for image_border
  latcoords=cv_coord(/to_rect,from_sph=transpose( $
    [[*sph_data.phi],[replicate(!dpi/2-thmin,sph_data.nlon)], $
    [replicate(imd2.rad,sph_data.nlon)]]))
  blobj=obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
   xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
  latcoords=cv_coord(/to_rect,from_sph=transpose( $
    [[*sph_data.phi],[replicate(!dpi/2-thmax,sph_data.nlon)], $
    [replicate(imd2.rad,sph_data.nlon)]]))
  blobj=[blobj,obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
    xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
  omodelmov->add,blobj

  ;  create outlines for longitudinal boundaries at bottom
  if imbound then begin
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph1,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(rmin,sph_data.nlat)]]))
    blobj=obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph2,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(rmin,sph_data.nlat)]]))
    blobj=[blobj,obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
    omodelmov->add,blobj
  endif

  ;  create outlines for latitudinal boundaries at bottom
  latcoords=cv_coord(/to_rect,from_sph=transpose( $
    [[*sph_data.phi],[replicate(!dpi/2-thmin,sph_data.nlon)], $
    [replicate(rmin,sph_data.nlon)]]))
  blobj=obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
   xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
  latcoords=cv_coord(/to_rect,from_sph=transpose( $
    [[*sph_data.phi],[replicate(!dpi/2-thmax,sph_data.nlon)], $
    [replicate(rmin,sph_data.nlon)]]))
  blobj=[blobj,obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
    xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
  omodelmov->add,blobj

  ;  create outlines for longitudinal boundaries at top
  if imbound then begin
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph1,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(rmax,sph_data.nlat)]]))
    blobj=obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    loncoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(ph2,sph_data.nlat)],[!dpi/2-(*sph_data.theta)], $
      [replicate(rmax,sph_data.nlat)]]))
    blobj=[blobj,obj_new('IDLgrPolyline',loncoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
    omodelmov->add,blobj
  endif

  ;  create outlines for latitudinal boundaries at top
  if imbound then begin
    latcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[*sph_data.phi],[replicate(!dpi/2-thmin,sph_data.nlon)], $
      [replicate(rmax,sph_data.nlon)]]))
    blobj=obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    latcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[*sph_data.phi],[replicate(!dpi/2-thmax,sph_data.nlon)], $
      [replicate(rmax,sph_data.nlon)]]))
    blobj=[blobj,obj_new('IDLgrPolyline',latcoords,color=whi,thick=thick, $
      xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)]
    omodelmov->add,blobj
  endif

  ;  create outlines for radial boundaries
  if imbound then begin
    radcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(phmin,sph_data.nr)],[replicate(!dpi/2-thmin,sph_data.nr)], $
      [*sph_data.rix]]))
    blobj=obj_new('IDLgrPolyline',radcoords,color=whi,thick=thick, $
     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    omodelmov->add,blobj
    radcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(phmax,sph_data.nr)],[replicate(!dpi/2-thmin,sph_data.nr)], $
      [*sph_data.rix]]))
    blobj=obj_new('IDLgrPolyline',radcoords,color=whi,thick=thick, $
     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    omodelmov->add,blobj
    radcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(phmin,sph_data.nr)],[replicate(!dpi/2-thmax,sph_data.nr)], $
      [*sph_data.rix]]))
    blobj=obj_new('IDLgrPolyline',radcoords,color=whi,thick=thick, $
     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    omodelmov->add,blobj
    radcoords=cv_coord(/to_rect,from_sph=transpose( $
      [[replicate(phmax,sph_data.nr)],[replicate(!dpi/2-thmax,sph_data.nr)], $
      [*sph_data.rix]]))
    blobj=obj_new('IDLgrPolyline',radcoords,color=whi,thick=thick, $
     xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
    omodelmov->add,blobj
  endif

endif

;  next create fieldline objects
if ptr_valid(sph_data.nstep) and (not keyword_set(nolines)) then begin

  nlines=n_elements(*sph_data.nstep)
  open=intarr(nlines)
  for i=0l,nlines-1 do begin

    ;  only draw lines that have line data
    ns=(*sph_data.nstep)(i)
    if ns gt 0 then begin

      ;  determine whether line is open or closed 
      if (max((*sph_data.ptr)(0:ns-1,i))-rmin)/(rmax-rmin) gt 0.99 then begin
        irc=get_interpolation_index(*sph_data.rix,(*sph_data.ptr)(0,i))
        ithc=get_interpolation_index( $
          *sph_data.lat,90-(*sph_data.ptth)(0,i)*!radeg)
        iphc=get_interpolation_index( $
          *sph_data.lon,((*sph_data.ptph)(0,i)*!radeg+360) mod 360)
        brc=interpolate(*sph_data.br,iphc,ithc,irc)
        if brc gt 0 then open(i)=1 else open(i)=-1
      endif  ;  else open(i)=0, which has already been done

      ;  flag those lines that go higher than the first radial gridpoint
      heightflag=max((*sph_data.ptr)(0:ns-1,i)) gt (*sph_data.rix)(1)

      ;  flag those lines that need to be drawn
      drawflag=(drawopen and (open(i) ne 0)) or (drawclosed and (open(i) eq 0))

      ;  create object for each valid line
      if (heightflag and drawflag) then begin

        ;  first determine color
        case open(i) of
          -1: col=red
           0: if keyword_set(for_ps) then col=bla else col=whi
           1: col=gre
        endcase

        ;  transform from spherical to cartesian coordinates
        linecoords=cv_coord(/to_rect,from_sph=transpose( $
          [[(*sph_data.ptph)(0:ns-1,i)],[!dpi/2-(*sph_data.ptth)(0:ns-1,i)], $
          [(*sph_data.ptr)(0:ns-1,i)]]))

        ;  create object
        flobj=obj_new('IDLgrPolyline',linecoords,color=col,thick=thick, $
          xcoord_conv=xs,ycoord_conv=ys,zcoord_conv=zs)
        if n_elements(ofieldlines) eq 0 then ofieldlines=flobj $
          else ofieldlines=[ofieldlines,flobj]

      endif

    endif
  
  endfor
endif

;  add fieldlines to the model
if n_elements(ofieldlines) gt 0 then omodelmov->add,ofieldlines

;  add extra objects to the model
nextra1=n_elements(extra_objects)
if nextra1 gt 0 then begin
  for k=0,nextra1-1 do begin
    extra_objects(k)->setproperty,xcoord_conv=xs,ycoord_conv=ys,$
      zcoord_conv=zs
    omodelmov->add,extra_objects(k),/alias
  endfor
endif
if ptr_valid(sph_data.extra_objects) then begin
  nextra2=n_elements(*sph_data.extra_objects)
  if nextra2 gt 0 then begin
    for k=0,nextra2-1 do begin
      (*sph_data.extra_objects)(k)->setproperty,xcoord_conv=xs,$
        ycoord_conv=ys,zcoord_conv=zs
      omodelmov->add,(*sph_data.extra_objects)(k),/alias
    endfor
  endif
endif

;  render the image and read it out to an image array
obuffer=obj_new('IDLgrBuffer',dim=[xsize,ysize])
obuffer->draw,oview
oimage=obuffer->read()
oimage->getproperty,data=outim

;  do the movie
if keyword_set(movie) then begin

  ;  allocate the array space
  outimarray=bytarr(3,xsize,ysize,360,/noz)
  outimarray(*,*,*,0)=outim

  print_time,'  ',1,360,slen,tst
  for i=1,359 do begin

    print_time,'  ',i+1,360,slen,tst

    ;  rotate the inner model by 1 degree
    omodelmov->rotate,[0,0,1],1

    ;  redraw the image and save
    obuffer->draw,oview
    oimage=obuffer->read()
    oimage->getproperty,data=outim2
    outimarray(*,*,*,i)=outim2

  endfor

  outim=temporary(outimarray)

endif

;  convert from 24-bit color to 8-bit color
if keyword_set(indexed_colors) then begin

  ;  color palette
  re=[bindgen(250),0b,0b,255b,255b,255b,255b]
  gr=[bindgen(250),255b,255b,0b,0b,255b,255b]
  bl=[bindgen(250),0b,0b,255b,255b,255b,255b]

  ;  do the conversion
  palindex=re+256l*gr+65536l*bl
  palsort=sort(palindex)
  paldupl=where(palindex(palsort(1:255)) eq palindex(palsort(0:255)),nwh)
  if nwh gt 0 then palindex(palsort(paldupl))=-1
  outindex=reform(outim(0,*,*)+256l*outim(1,*,*)+65536l*outim(2,*,*))
  if keyword_set(movie) then outim=make_array(dim=[xsize,ysize,360],/byte) $
    else outim=make_array(dim=[xsize,ysize],/byte)
  for i=0b,255 do outim=temporary(outim)+i*(outindex eq palindex(i))

endif

;  display onscreen
if keyword_set(onscreen) then begin
  if keyword_set(indexed_colors) then begin
    tvlct,re,gr,bl
    if keyword_set(movie) then scim,outim(*,*,0),/quiet else scim,outim,/quiet
  endif else begin
    if keyword_set(movie) then scim,outim(*,*,*,0),/true,/quiet $
      else scim,outim,/true,/quiet
  endelse
endif

;  cleanup
if ptr_valid(sph_data.nstep) and (not keyword_set(nolines)) then $
  obj_destroy,ofieldlines
obj_destroy,oimage
obj_destroy,obuffer
obj_destroy,oview
obj_destroy,omodeltop
obj_destroy,omodelin
obj_destroy,omodelmov
if not keyword_set(noimage) then obj_destroy,otexmap
if not keyword_set(noimage) then obj_destroy,osphere

end