plot3d_new.pro (Documentation for /home/davin/idl/socware/)

;+
;NAME:
;  PLOT3D_NEW
;PROCEDURE: plot3d_new, dat , latitude, longitude
;PURPOSE:
;  Draws a series of 3d color plots, one plot per energy step.
;INPUT:
;  dat:  3d data structure.
;  latitude:  latitude of center of plot
;  longitude:  longitude of center of plot
;KEYWORDS:
;   EBINS:     Specifies energy bins to plot.
;   SUM_EBINS: Specifies how many bins to sum, starting with EBINS.  If
;              SUM_EBINS is a scalar, that number of bins is summed for
;              each bin in EBINS.  If SUM_EBINS is an array, then each
;              array element will hold the number of bins to sum starting
;              at each bin in EBINS.  In the array case, EBINS and SUM_EBINS
;              must have the same number of elements.
;   BNCENTER:  if > 0 then lat,lon set so that B direction is in center
;              if < 0 then lat,lon set so that -B direction is in center
;              ('magf' element must exist in dat structure. See "ADD_MAGF")
;   BINS:    bins to use for autoscaling.
;   SMOOTH:  smoothing parameter.  0=none,  2 is typical
;   TITLE:   overrides default plot title.
;   NOTITLE: if set, suppresses plot title.
;   NOCOLORBAR: if set, suppresses colorbar.
;   NOBORDER: if set, suppresses plot border.
;SEE ALSO:  "PLOT3D_OPTIONS" to see how to set other options.
;
;;NOTE: plot3d_new has been created to eliminate a naming conflict between plot3d and a library routine that was added to IDL in version 8.1.  
;      The original plot3d.pro remains to preserve backwards compatibility for users using IDL 7.1 or earlier.
;
;CREATED BY:	Davin Larson
;LAST MODIFICATION:	@(#)plot3d.pro	1.43 02/04/17
;-
pro plot3d_new,tbindata, latitude,longitude,rotation, $
   BINS=bins,            $
   TBINS=tbins,          $
   BNCENTER=bncenter,    $
   LABELS=labs,          $
   SMOOTH= smoth,        $
   EBINS=eb,             $
   SUM_EBINS=seb,        $
   PTLIMIT=plot_limits,  $
   MAP= ptmap,           $
   UNITS = units,        $
   SETLIM = setlim,      $
   PROC3D = proc3d,      $
   ZRANGE = zrange,      $
   ZERO = zero,          $
   STACK = stack,        $     ;!p.multi[1:2]
   ROW_MAJOR = row_major,$     ;!p.multi[4]
   NOERASE = noerase,    $     ;!p.multi[0]
   TITLE = title,	 $
   NOTITLE = notitle,    $
   NOCOLORBAR = nocolorbar,	$
   NOBORDER = noborder,$
   log=log

@plot3d_com.pro
plot3d_options          ; Set Defaults

if size(/type,tbindata) ne 8 || tbindata.valid eq 0 then begin
  dprint,'Invalid data'
  return
endif

if not keyword_set(units) then units='Eflux'
bindata = conv_units(tbindata,units)

str_element,bindata,'magf',value=magf

if keyword_set(magf) then begin
   xyz_to_polar,magf,theta=th,phi=phi
   bdir = [th,phi]
endif

str_element,bindata,'vsw',value=vsw

if n_elements(latitude)  eq 0 then latitude = latitude_3d
if n_elements(longitude) eq 0 then longitude = longitude_3d
if n_elements(rotation) eq 0 then rotation = 0d
if n_elements(smoth)     eq 0 then smoth = smooth_3d
if n_elements(log) eq 0 then log = logscale_3d

if keyword_set(bncenter) and n_elements(bdir) eq 2 then begin
   latitude = bdir(0)
   longitude = bdir(1)
   if bncenter lt 0 then begin
      latitude = - latitude
      longitude = longitude + 180
   endif
endif

n_e = bindata.nenergy
nbins = bindata.nbins
if keyword_set(bins) eq 0 then str_element,bindata,'bins',bins

if n_elements(bins) eq (nbins*n_e) then begin
	if n_elements(tbins) ne nbins then tbins = total(bins,1) gt 0
	bins = total(bins,2) gt 0
	binsindx = where(tbins eq 0,binscnt)
	if binscnt eq 0 then tbins = 0
	binsindx = where(bins eq 0,binscnt)
	if binscnt eq 0 then bins = 0
endif

if n_elements(bins) ne nbins then bins = replicate(1,bindata.nbins)
goodbins = where(bins,ngood)
badbins = where(bins eq 0,nbad)
if nbad ne 0 then bindata.data(*,badbins) = !values.f_nan

triang = keyword_set(triang_3d)

str_element,bindata,'map',value=ptmap
if triang eq 0 then begin
 if keyword_set(ptmap) then begin
   np = dimen1(ptmap)
   ptmap = rotate(ptmap,7)  ; this line and the next should be removed when
   ptmap = rotate(ptmap,5)  ; the c-code is fixed
   ptmap = reverse(ptmap,2) ; assumes !order = 0
   dprint,'Using supplied mapping'
 endif else begin
   ares = 16
   ptmap = make_3dmap(bindata,64*ares,32*ares)
   ptmap = shift(ptmap,32*ares,0)  ; map routine expects phi: -180 to 180
 endelse
endif



; handle ebins
if n_elements(eb) eq 0 then                $  ; undefined ebins
    start=indgen(n_e)                      $
else                                       $  ; scalar or 1-D ebins
    start=eb

; handle sum_ebins
if keyword_set(seb) eq 0 then              $  ; SUM_EBINS keyword not set
    seb=1                                  $
else begin
    if n_elements(seb) gt 1 and (n_elements(seb) ne n_elements(eb)) then begin
        message,'if SUM_EBINS is an array, it must have same dimensions'+ $
        ' as EBINS.'
    endif
endelse

; get stop bins
    stop = start + seb - 1

n_p = dimen1(start)  ; number of plots

stacking= [ [1,1],   [1,1],   [1,2],   [2,2],   [2,2],   [2,3],   [2,3],  $
   [3,3],   [3,3],   [3,3],   [3,4],   [3,4],   [3,4],   [3,5],   [3,5],  $
   [3,5],   [4,4],   [4,5],   [4,5],   [4,5],   [4,5],   [5,5],   [5,5],  $
   [5,5],   [5,5],   [5,5],   [5,6],   [5,6],   [5,6],   [5,6],   [5,6],  $
   [6,6],   [6,6],   [6,6],   [6,6],   [6,6],   [6,6],   [6,7],   [6,7],  $
   [6,7],   [6,7],   [6,7],   [6,7],   [7,7],   [7,7],   [7,7],   [7,7],  $
   [7,7],   [7,7],   [7,7],   [7,8],   [7,8],   [7,8],   [7,8],   [7,8],  $
   [7,8],   [7,8],   [8,8],   [8,8],   [8,8],   [8,8],   [8,8],   [8,8],  $
   [8,8],   [8,8]]

if n_elements(stack) ne 2 then stack=stacking(*,n_p)

if not keyword_set(noerase)   then noerase   = 0 else noerase   = 1
if not keyword_set(row_major) then row_major = 0 else row_major = 1
!p.multi=[noerase,stack(0),stack(1),0,row_major]

if !p.charsize ne 0 then charsize = !p.charsize else charsize = 1
!y.omargin = [1,4]
!x.omargin = [0,6]
if !p.multi(1) gt 2 or !p.multi(2) gt 2 then begin
;   charsize = !p.charsize *2
   !y.omargin = !y.omargin * 2
   !x.omargin = !x.omargin * 2
endif

xposmax = 0.
yposmax = 0.
yposmin = 1.



for plot_num=0,n_p-1 do begin
   e1 = start(plot_num)
   en1 = total(bindata.energy(e1,*))/n_elements(bindata.energy(e1,*))

   e2 = stop(plot_num)
   en2 = total(bindata.energy(e2,*))/n_elements(bindata.energy(e2,*))

   if e1 ne e2 then begin
     ttl = string(en1)+' - '+string(en2)+' eV'
   endif else begin
     ttl = strcompress(string(en1,format='(F12.1)')+' eV', $
           /remove_all)
   endelse
   ttl=strcompress(ttl)

nsteps = e2-e1+1

mapname = strlowcase(strmid(mapname_3d,0,3))
if keyword_set(setlim) then begin
   map_set, latitude,longitude,0., /WHOLE_MAP, $
     /noerase, $
     MOLL  = mapname eq 'mol' ,     $
     CYLIN = mapname eq 'cyl' ,     $
     ORTHO = mapname eq 'ort' ,     $
     AITOFF= mapname eq 'ait' ,     $
     LAMBERT=mapname eq 'lam' ,     $
     GNOMIC =mapname eq 'gno' ,     $
     MERCATO=mapname eq 'mer' ,     $
     LIMIT=plot_limits        , $
     TITLE=ttl, $
     charsize = charsize, $
     XMARGIN=[1,1], YMARGIN=[1,1], noborder=noborder
     return
endif


data = total(bindata.data(e1:e2,*),1)  ; sum over energies

if keyword_set(zrange) then range = zrange[*]  $
else  begin
   range = minmax(data(goodbins),/pos)
endelse

if keyword_set(triang) then begin
;   data = total(bindata.data(e1:e2,*),1)  ; sum over energies
;   range = minmax(data(goodbins),/pos)

   theta = total(bindata.theta(e1:e2,*),1)/nsteps
   phi = total(bindata.phi(e1:e2,*),1)/nsteps

   if keyword_set(tbins) then  ind = where(tbins ne 0)  $
   else ind = where(finite(theta) and finite(phi))

if ind(0) eq -1 then goto, contu
   data  = data(ind)
   theta = theta(ind)
   phi   = phi(ind)
   triangulate,phi,theta,fvalue=data,sphere=sphere,/degrees
   image = trigrid(data,sphere=sphere,[2.,2.],[-180.,-90.,178.,90.],/deg)
   image = bytescale(image,log=log, range=range,zero=zero)

endif else begin
;   data = total(bindata.data(e1:e2,*),1)  ; sum over energies
;   range = minmax(data(goodbins))
   data = [!values.f_nan,data]
;   data = [0.,data]
   data = data(ptmap+1)                     ; make 2-dimensional
   if( keyword_set(smoth)) then begin
      dim = dimen(data) * 1*smoth   ; increase the dimensions
      data = rebin(data,dim(0),dim(1))
      data = smooth_periodic(data,smoth)
   endif

   image = bytescale(data,log=log, range=range,zero=zero)

   image = congrid(image,360,180)  ; assumes  !order eq 0
endelse

   mapname = strlowcase(strmid(mapname_3d,0,3))
   map_set, latitude,longitude,rotation, /advance, /WHOLE_MAP, $
     MOLL  = mapname eq 'mol' ,     $
     CYLIN = mapname eq 'cyl' ,     $
     ORTHO = mapname eq 'ort' ,     $
     AITOFF= mapname eq 'ait' ,     $
     LAMBERT=mapname eq 'lam' ,     $
     GNOMIC =mapname eq 'gno' ,     $
     MERCATO=mapname eq 'mer' ,     $
     LIMIT=plot_limits        , $
     TITLE=ttl, $
     charsize = charsize, $
     XMARGIN=[1,1], YMARGIN=[1,1], noborder = noborder

   new = map_image(image,sx,sy,xsize,ysize, $
      COMPRESS=compress_3d, BILIN=0)

;   new = map_patch(image,xstart=sx,ystart=sy,xsize=xsize,ysize=ysize)

   tv,new,sx,sy,xsize=xsize,ysize=ysize,order=0

   if n_elements(bdir) eq 2 then begin     ; draw mag field
      plots,bdir(1),bdir(0),psym=1
      plots,bdir(1)+180,-bdir(0),psym=4
   endif
   if n_elements(vsw) then begin
      xyz_to_polar,vsw,theta=th,phi=ph
      plots,ph,th,psym=2
   endif else  plots,179.99,0.,psym=2   ; draw the sun.


   if keyword_set(grid_3d) then map_grid,londel=grid_3d(0),latdel=grid_3d(1)

   str_element,opts,'plot3d_proc',value=proc3d
   str_element,bindata,'plot3d_proc',value=proc3d
   if keyword_set(proc3d) then begin
       proc = ''
       if ndimen(proc3d) eq 0 then proc=proc3d else begin
           if plot_num lt n_elements(proc3d) then proc = proc3d(plot_num)
       endelse
       if keyword_set(proc) then call_procedure,proc,bindata,en2,plot_num
   endif

   if keyword_set(labs) then begin
      lab=strcompress(indgen(nbins),/rem)
      colorcode = [!p.background,!p.color]
      xyouts,phi,theta,lab,align=.5,COLOR=colorcode(bins)
   endif

   ;map_grid,/label,latdel =45,londel=45,CHARSIZE=charsize
   ;map_continents
   if xposmax lt !x.window(1) then xposmax=!x.window(1)
   if yposmax lt !y.window(1) then yposmax=!y.window(1)
   if yposmin gt !y.window(1) then yposmin=!y.window(0)
   contu:
endfor

if keyword_set(notitle) then ttl = '' else $
	if not keyword_set(title) then begin
        str_element,bindata,'trange',trange
;        if not keyword_set(trange) then $
           trange=[bindata.time,bindata.end_time]
		ttl = bindata.project_name+'  '+bindata.data_name
	;	ttl =ttl+'!C'+ trange_str(bindata.time,bindata.end_time)
		ttl =ttl+'!C'+ trange_str(trange)
	endif else ttl = title

tcharsize = 1.25 * charsize
titleheight = tcharsize * float(!d.y_ch_size)/!d.y_size
xyouts, 0.5,1.-titleheight,ttl,/normal,alignment=.5,charsize=tcharsize

time_stamp

!p.multi = 0
!y.omargin=0
!x.omargin=0

space = charsize * float(!d.x_ch_size)/!d.x_size
pos =[xposmax+space,yposmin,xposmax+3*space,yposmax]

if not keyword_set(nocolorbar) then $
	draw_color_scale,range = [0.,1.],pos=pos,chars=charsize

return
end