; $Id: tas2d.pro,v 1.2 2002/01/14 20:59:26 amunter Exp $
PRO put_new_label, id, label, delimeter
;replaces everything in the widget label after the delimeter by "label"
;Tries to keep the same overall string length
  WIDGET_CONTROL, id, GET_VALUE=oldlabel
  loc=STRPOS(oldlabel,delimeter)+1
  IF (loc NE -1) THEN BEGIN
   newlabel=STRMID(oldlabel,0,loc)+label
   WIDGET_CONTROL, id, SET_VALUE=newlabel
  END
  RETURN                
END



FUNCTION plotit, hardcopy=hardcopy, which=which, lastcall=lastcall

COMMON stuff, d, info
  
  IF keyword_set(hardcopy) THEN hardcopy=1 ELSE hardcopy=0
  IF keyword_set(which) THEN which=which ELSE which = 0
  ;0=plot both, 1=plot contour, 2=plot surface, 3=plot both
  ;****** Get the indices for the plot variables for the TASCOM data
  xind=tasid(info.x,d)
  IF (xind eq -1) THEN xind=0
  zind=tasid(info.z,d)
  IF (strupcase(info.z) eq 'I') then sind=tasid('SIGI',d) ELSE $
      IF (strupcase(info.z) eq 'MON') THEN sind=tasid('SIGM',d) ELSE sind=-1
  yind=tasid(info.y,d)
  tind=tasid('TIM',d)
  if (tind ge 0) then sd=stdev(d.data(tind,*),mean) else mean=0

  ;***** Set up min and max x *****
  minr=(0.001 > 0.01*max(d.data(xind,*)))
  dif=((minr) > (0.1*(max(d.data(xind,*))-min(d.data(xind,*)))))/2
  if (info.new) then info.xmin=min(d.data(xind,*))-dif
  if (info.new) then info.xmax=max(d.data(xind,*))+dif

  minr=(0.001 > 0.01*max(d.data(yind,*)))
  dif=((minr) > (0.1*(max(d.data(yind,*))-min(d.data(yind,*)))))/2
  if (info.new) then info.ymin=min(d.data(yind,*))-dif
  if (info.new) then info.ymax=max(d.data(yind,*))+dif

  ;***** Set up interpolated surface *****
  triangulate,d.data(xind,*),d.data(yind,*),tr,b       
  rangex = (info.xmax-info.xmin)/info.gridsize
  rangey = (info.ymax-info.ymin)/info.gridsize
  inter=trigrid(d.data(xind,*),d.data(yind,*),d.data(zind,*),$
        tr,[rangex,rangey],[info.xmin,info.ymin,info.xmax,info.ymax])
  

  intersize=size(inter)
  xvals=rangex*FINDGEN(intersize(1))
  yvals=rangey*FINDGEN(intersize(2))
  xvals=xvals+info.xmin
  yvals=yvals+info.ymin
  
  ;****** Calculate the volume under the surface in a primitive fashion *****
  vol=0
  i=INT(0) & j=INT(0)
  FOR i=0,intersize(1)-1,1 DO FOR j=0,intersize(2)-1,1 DO vol = vol + inter(i,j)
  vol = vol*rangex*rangey/intersize(1)/intersize(2)
  print, 'Volume = ',vol
  
  IF keyword_set(lastcall) THEN return, {x:xvals, y:yvals, z:inter}

  ;****** Plot it, baby! ******
  IF ((which EQ 0) OR (which EQ 1) OR (which EQ 3)) THEN BEGIN
  ;First the contour
   IF NOT keyword_set(hardcopy) THEN BEGIN
        WIDGET_CONTROL, GET_VALUE=who, info.canvas1
        WSET, who
   ENDIF


   
  ;****** Figure out a nice set of levels ******
  zmax = MAX(inter)	;Get the max value
  del=ABS(FLOAT(zmax/info.nlevels)) ;Nominal interval between levels
  IF (del NE 0) THEN BEGIN
  	jim=10.0^(FIX(ALOG10(del)))
  	del=FIX(del/jim)*FIX(jim)
  END ELSE del=1
  IF del EQ 0 THEN del=1
  levels = FINDGEN(info.nlevels)*del
  levels = levels + del
  clabels = FINDGEN(info.nlevels)
  clabels(*) = 1
  ;inter = min_curve_surf(inter,xvals,yvals)
  contour, inter,xvals,yvals,/FOLLOW,LEVELS=levels,c_labels=clabels,$
  xtitle='!6'+info.x,ytitle='!6'+info.y,$
  yrange=[info.ymin,info.ymax],xrange=[info.xmin,info.xmax],title=info.title,$
  xstyle=1,ystyle=1,xminor=5,yminor=5,$
  xtickformat='(F7.3)',ytickformat='(F7.3)'
  END
  
 IF ((which EQ 0) OR (which EQ 2) OR (which EQ 3)) THEN BEGIN
   ;Now the surface
   IF NOT keyword_set(hardcopy) THEN BEGIN
        WIDGET_CONTROL, GET_VALUE=who, info.canvas2
        WSET, who
   END

 
;   info.title = '100/010 Reflection'
;   shade_surf, inter,xvals,yvals,$
;        ztitle='!6'+'Intensity'+'(<Time>='+ string(mean,format='(i4)') + ' seconds)',$
;        xtitle='!6'+'H (r.l.u.)',ytitle='!6'+'K (r.l.u.)',title=info.title,$
;        yrange=[info.ymin,info.ymax],xrange=[info.xmin,info.xmax],$
;        ;zrange=[0,max(d.data(zind,*))],$
;        xstyle=1,ystyle=1,NOERASE=info.noerase,$
 ;       ytickformat='(F6.2)',xtickformat='(F4.2)'

   shade_surf, inter,xvals,yvals,$
        ztitle='!6'+info.z+'(<Time>='+ string(mean,format='(i4)') + ' seconds)',$
        xtitle='!6'+info.x,ytitle='!6'+info.y,title=info.title,$
        yrange=[info.ymin,info.ymax],xrange=[info.xmin,info.xmax],$
        zrange=[0,max(d.data(zind,*))],$
        xstyle=1,ystyle=1,NOERASE=info.noerase,charsize=1.5
        ;xtickformat='(G6.2)',ytickformat='(G6.2)'
  END
  jims = 1
  RETURN, jims
END


PRO level_event, ev
COMMON stuff, d, info
 info.nlevels=ev.value
 jim=plotit(which=1)
RETURN
END 
PRO grid_event, ev
COMMON stuff, d, info
 info.gridsize=ev.value
 jim=plotit(which=3)
RETURN
END

FUNCTION tas2d,inputdata,xaxis=x,yaxis=y,zaxis=z,title=title,noerase=noerase,$
	xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax,$
        gridsize=gridsize, nlevels=nlevels, nokill=nokill, color=color, twopage=twopage, $
        logscale=logscale
;+
; NAME:         TAS2D
; PURPOSE:
;               Plot a CONTOUR map and a SURFACE of a TASCOM file read in with TASRDM
; CATEGORY:
; CALLING SEQUENCE:
;               output=tas2d(tasdata)
; INPUTS:
;               tasdata = a tascom data structure read in with TASRDM
;               
; OPTIONAL INPUT PARAMTERS: none
; KEYWORD PARAMETERS:
;               xaxis,yaxis,zaxis = which axes to plot. Must be a string contained
;			in the tascom data structure: tasdata.PRP.
;			Default: xaxis='H',yaxis='K',zaxis='I', if they exist.
;			  Otherwise xaxis=tasdata.PRP(0), yaxis=tasdata.PRP(1), 
;				zaxis='I'
;		title=title of the plot. Default: title='firstfile--lastfile'
;		xmin,xmax,ymin,ymax = plot limits. Default: Extent of the data
;		gridsize = number of grid points. Default: gridsize=50
;		nlevels = number of levels for the contour plot. 
;			Default: nlevels=10
;		/nokill	If set, then the plot remains after quitting
;		/color	If set, then color postscript code will be used.
;			Note that a regular postscript printer can usually
;			interpret color postscript code. On the laserjet 4,
;			color code results in a nicer image, eventhough it
;			isn't a color printer.
;		/twopage If set, then the hardcopy will put the contour and 
;			surface on seperate pages
;		
;                  
; OUTPUTS:
;               a contour and surface plot of the data.
; OPTIONAL OUTPUT PARAMETERS:
; COMMON BLOCKS: stuff, d, info	
;		This is a local common block used to pass info to the
;		event handlers.
;		d is the tascom data, and info is a plotinfo structure
;		for surface data
;
; SIDE EFFECTS: no idea
; RESTRICTIONS:
;               
; PROCEDURE:
; MODIFICATION HISTORY:
;               Paul Schleger, Riso, August 1993
;-

COMMON stuff, d, info

d = inputdata

;****** Process keywords ******
IF NOT keyword_set(slices) THEN slices=0
IF NOT keyword_set(NOERASE) THEN  noerase=0
IF NOT keyword_set(twopage) THEN twopage=0
IF NOT keyword_set(nlevels) THEN nlevels=10
IF NOT keyword_set(title) THEN title=d.text(1)
new=((not keyword_set(xmin)) and (not keyword_set(xmax)))
new=(new and (not keyword_set(ymin)) and (not keyword_set(ymax)))


;****** First rename all H0, K0, L0 in d.prp string to H K L ******
length=(size(d.prp))(1)
FOR I=0,length-1,1 DO IF (STRPOS(d.prp(I),'H0') NE -1) THEN d.prp(I)='H'
FOR I=0,length-1,1 DO IF (STRPOS(d.prp(I),'K0') NE -1) THEN d.prp(I)='K'
FOR I=0,length-1,1 DO IF (STRPOS(d.prp(I),'L0') NE -1) THEN d.prp(I)='L'

;****** Figure out what to plot ******
;default values:
Hid = tasid('H',d)
Kid = tasid('K',d)
Lid = tasid('L',d)
IF ((Hid NE -1) AND (Kid NE -1) AND (Lid NE -1)) THEN BEGIN
  IF (min(d.data(Hid,*)) NE max(d.data(Hid,*))) THEN hasH=1 ELSE hasH=0
  IF (min(d.data(Kid,*)) NE max(d.data(Kid,*))) THEN hasK=1 ELSE hasK=0
  IF (min(d.data(Lid,*)) NE max(d.data(Lid,*))) THEN hasL=1 ELSE hasL=0

  IF NOT keyword_set(zaxis) THEN z='I'
  IF ((NOT keyword_set(xaxis)) AND (NOT keyword_set(yaxis))) THEN BEGIN
    IF (hasH AND hasK) THEN BEGIN
        x='H'
        y='K'
    ENDIF ELSE IF (hasH AND hasL) THEN BEGIN
        x='H'
        y='L'
    ENDIF ELSE BEGIN
        x='K'
        y='L'
    ENDELSE
  ENDIF
ENDIF ELSE BEGIN
	IF (NOT keyword_set(xaxis)) THEN x=d.prp(0) ELSE x=x
	IF (NOT keyword_set(yaxis)) THEN y=d.prp(1) ELSE y=y
	IF (NOT keyword_set(zaxis)) THEN z='I'      ELSE z=z
ENDELSE

;Figure out a reasonable grid
IF NOT keyword_set(gridsize) THEN BEGIN
  	IF (d.data(tasid(x,d),0) EQ d.data(tasid(x,d),1)) THEN scanvar=d.data(tasid(x,d),*) $
   	ELSE	scanvar=d.data(tasid(y,d),*)
   	nxvals=where(scanvar EQ scanvar(2),count)
   	if (count NE 0) AND (n_elements(nxvals) GT 1) THEN gridsize=min([nxvals(1),n_elements(nxvals)]) ELSE gridsize=50
   	print, gridsize
END 
IF (gridsize LE 1) THEN gridsize=50   

base = WIDGET_BASE(TITLE='2D TASCOM Plot',/COLUMN)
lcol = WIDGET_BASE(base, /FRAME, /ROW)
rcol = WIDGET_BASE(base, /FRAME, /ROW)
junk = { CW_PDMENU_S, flags:0, name:'' }
plotaxes = REPLICATE(junk,3*(length+1)+7)
plotaxes(0) = {CW_PDMENU_S, 1, 'x-axis:'+x+'         '}
plotaxes(length+1) = {CW_PDMENU_S, 1, 'y-axis:'+y+'        '}
plotaxes(2*length+2) = {CW_PDMENU_S, 1, 'z-axis:'+z+'        '}
 FOR I=1,length-1,1 DO BEGIN
        plotaxes(I) = {CW_PDMENU_S, 0, d.prp(I-1)}
        plotaxes(I+length+1) = {CW_PDMENU_S, 0, d.prp(I-1)}
        plotaxes(I+2*length+2) = {CW_PDMENU_S, 0, d.prp(I-1)}
END
plotaxes(length) = {CW_PDMENU_S, 2, d.prp(length-1)}
plotaxes(2*length+1) = {CW_PDMENU_S, 2, d.prp(length-1)}
plotaxes(3*length+2) = {CW_PDMENU_S, 2, d.prp(length-1)}
plotaxes(3*length+3) = {CW_PDMENU_S, 0, 'Plot it'}
plotaxes(3*length+4) = {CW_PDMENU_S, 0, 'Zoom in'}
plotaxes(3*length+5) = {CW_PDMENU_S, 0, 'Zoom out'}
plotaxes(3*length+6) = {CW_PDMENU_S, 0, 'Hardcopy'}
plotaxes(3*length+7) = {CW_PDMENU_S, 2, 'Quit'}

menuaxes = cw_pdmenu(lcol,plotaxes,/RETURN_FULL_NAME,IDS=axesID)
gridslider = WIDGET_SLIDER(lcol,MINIMUM=1,MAXIMUM=500,VALUE=gridsize,TITLE='grid elements',EVENT_PRO='grid_event')
levelslider = WIDGET_SLIDER(lcol,MINIMUM=1,MAXIMUM=30,VALUE=nlevels,TITLE='no. of levels',EVENT_PRO='level_event')

sx = 400
sy = 400
canvas1=WIDGET_DRAW(rcol,xsize=sx,ysize=sy)
canvas2=WIDGET_DRAW(rcol,xsize=sx,ysize=sy)

;****** fill up the plotinfo structure ******
info = {plotinfo,x:x,y:y,z:z,noerase:noerase,gridsize:gridsize,nlevels:nlevels,new:new,$
        xmin:0.,xmax:0.,ymin:0.,ymax:0.,canvas1:canvas1,canvas2:canvas2,title:title}
IF keyword_set(xmin) THEN info.xmin=xmin
IF keyword_set(xmax) THEN info.xmax=xmax
IF keyword_set(ymin) THEN info.ymin=ymin
IF keyword_set(ymax) THEN info.ymax=ymax

IF keyword_set(logscale) THEN BEGIN
	jim=abs(d.data(tasid(z,d),*))
	nocounts=where(jim EQ 0,counts)
	IF counts NE 0 THEN jim(nocounts)=1
	jim=ALOG10(jim)-ALOG10(min(jim))
	d.data(tasid(z,d),*)=jim
END

WIDGET_CONTROL, /REALIZE, base
REPEAT BEGIN
        ev = WIDGET_EVENT(base)
        IF (STRPOS(ev.value,'x-axis') NE -1) THEN BEGIN
                x=STRTRIM(STRMID(ev.value,STRPOS(ev.value,'.')+1,5),2)
                put_new_label, axesID(0),x,':'
                info.x=x
        END
        IF (STRPOS(ev.value,'y-axis') NE -1) THEN BEGIN
                y=STRTRIM(STRMID(ev.value,STRPOS(ev.value,'.')+1,5),2)
                put_new_label,axesID(length+1),y,':'
                info.y=y
        END
        IF (STRPOS(ev.value,'z-axis') NE -1) THEN BEGIN
                z=STRTRIM(STRMID(ev.value,STRPOS(ev.value,'.')+1,5),2)
                put_new_label,axesID(2*length+2),z,':'
                info.z=z
        END
        IF (ev.value EQ 'Plot it') THEN jims = plotit(which=3)
        IF ev.value EQ 'Hardcopy' THEN BEGIN
                ;name=choose_name()
                IF NOT keyword_set(twopage) THEN BEGIN
                	psinit, FILE='tas2d.ps', /VECTOR, WINDOW=[1,4.54,1,3.56] 
                	pmultisave = !P.MULTI
        		!P.MULTI = [0,1,2]
        		device,YSIZE=25.0,XSIZE=16.0,$
        			XOFFSET=1.5,YOFFSET=1.5
        	END ELSE psinit, /LANDSCAPE, FILE='tas2d.ps', /VECTOR, WINDOW=[1,4.54,1,3.56] 
                IF keyword_set(color) THEN device,/COLOR $
                ELSE device,BITS_PER_PIXEL=8
                DEVICEPS, FREQUENCY=60, DOTS_PER_INCH=600
                jims=plotit(/hardcopy,which=3)
                IF NOT keyword_set(twopage) THEN !P.MULTI=pmultisave
                psterm,/nostamp
        END
        IF ev.value EQ 'Zoom in' THEN BEGIN
                WIDGET_CONTROL, GET_VALUE=who, info.canvas1
                WSET, who
                r=recroi(sx,sy,xroi,yroi)
                loadct, 0
                rng = convert_coord([info.xmin,info.xmax],[info.ymin,info.ymax],/DATA, /TO_DEVICE)
                info.xmin = (min(xroi)-rng(0,0))*(info.xmax-info.xmin)/(rng(0,1)-rng(0,0))+info.xmin
                info.xmax = (max(xroi)-rng(0,0))*(info.xmax-info.xmin)/(rng(0,1)-rng(0,0))+info.xmin
                info.ymin = (min(yroi)-rng(1,0))*(info.ymax-info.ymin)/(rng(1,1)-rng(1,0))+info.ymin
                info.ymax = (max(yroi)-rng(1,0))*(info.ymax-info.ymin)/(rng(1,1)-rng(1,0))+info.ymin
                info.new=0
                jims=plotit(which=3)
        END
        IF ev.value EQ 'Zoom out' THEN BEGIN
                info.xmin = 0;min(d.data(tasid(x,d),*))
                info.xmax = 0;max(d.data(tasid(x,d),*))
                info.ymin = 0;min(d.data(tasid(y,d),*))
                info.ymax = 0;max(d.data(tasid(y,d),*))
                info.new=1
                jims=plotit(which=3)
        END             
END UNTIL ev.value EQ 'Quit'
jims=plotit(/lastcall)
IF NOT keyword_set(nokill) THEN WIDGET_CONTROL, /DESTROY, base 
RETURN, jims
END
