; $Id: annsca_procs.pro,v 1.7 2005/02/03 18:05:03 azuah Exp $
;###############################################################################
; LICENSE:
;  The software in this file is written by an employee of 
;  National Institute of Standards and Technology 
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;************************************************************************************************
;
; These procedures are used by the program "annsca" (q.v.)
;
;************************************************************************************************
pro annsca_checkinput,checkstring,vals,ok,messages
;************************************************************************************************
;
compile_opt strictarr
;
;
ok=1
;
messages=''
case checkstring of
	"angles": begin
		philo=vals[0]
		phihi=vals[1]
		dphi=vals[2]
		if (philo ge phihi) then begin
			ok=0
			messages=[messages,'Start angle must be less than stop angle!']
		endif
		if (dphi le 0.d0) then begin
			ok=0
			messages=[messages,'Step in angle must be positive!']
		endif
	end
;
	"sigt_ins": begin
		nann=n_elements(vals)
		for k = 0,nann-1 do begin
		  if (vals[k] lt 0.d) then begin
				ok=0
    		messages=[messages,'SigT_in cannot be negative!']
  		endif
		endfor
		if (vals[nann-1] eq 0.0d) then begin
			ok=0
			messages=[messages,'SigT_in must be positive for the outermost annulus!']
		endif
	end
;
	"sigt_outs": begin
		nann=n_elements(vals)
		for k = 0,nann-1 do begin
		  if (vals[k] lt 0.d) then begin
    		ok=0
    		messages=[messages,'SigT_out cannot be negative!']
  		endif
		endfor
		if (vals[nann-1] eq 0.0d) then begin
			ok=0
			messages=[messages,'SigT_out must be positive for the outermost annulus!']
		endif
	end
;
	"radii": begin
		rprev = 0.d0
		nann=n_elements(vals)
		for k = 0,nann-1 do begin
		  if (vals[k] le rprev) then begin
	  		if (k eq 0) then begin
  				ok=0
	  	  	messages=[messages,'First radius must be positive!']
	  		endif else begin
		  		ok=0
		    	messages=[messages,'Radii must increase monotonically!']
				endelse
			endif
    	rprev=vals[k]
		endfor
	end
;
	"ybeamsize": begin
		beam_wid=vals[0]
		if (beam_wid lt 0.0) then begin
			ok=0
			messages=[messages,'Beam width cannot be negative!']
		endif
		beam_ydispl=vals[1]
		max_sam_radius=vals[2]
		y_min=beam_ydispl-beam_wid*0.5
		y_max=beam_ydispl+beam_wid*0.5
		if (y_min ge max_sam_radius or y_max le -max_sam_radius) then begin
			ok=0
  		messages=[messages,'Beam completely misses sample!']
		endif
	end
;
	"zbeamsize": begin
		beam_hgt=vals[0]
		if (beam_hgt lt 0.0) then begin
			ok=0
			messages=[messages,'Beam height cannot be negative!']
		endif
		beam_zdispl=vals[1]
		sam_halfheight=vals[2]
		z_min=beam_zdispl-beam_hgt*0.5
		z_max=beam_zdispl+beam_hgt*0.5
		if (z_min ge sam_halfheight or z_max le -sam_halfheight) then begin
			ok=0
  		messages=[messages,'Beam completely misses sample!']
		endif
	end
;
	else:
endcase
;
end


;************************************************************************************************
pro annsca_ShowGeometry,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,get_uvalue = pState
;
(*pState).overplot=0
widget_control,(*pState).butt_oplot,set_value="Overplot"
;
if ((*pState).plotting_geometry eq 0) then begin
	(*pState).plotting_geometry=1
	(*pState).autoscale=1
endif
;
; Get the radii of the circles.
radius = (*pState).radius
ncircles=(*pState).nann
if (ncircles eq 0) then return
;
; Number of points on each circle.
np = 100
theta = dindgen(np)/(np-1.)*360.
xc = fltarr(ncircles,np)
yc = fltarr(ncircles,np)
for i = 0,ncircles-1 do begin
  xc[i,0:np-1] = radius[i]*cos(theta[0:np-1] * !dtor)
  yc[i,0:np-1] = radius[i]*sin(theta[0:np-1] * !dtor)
endfor
;
!x=(*pState).bangx
!y=(*pState).bangy
!p=(*pState).bangp
;
; Change the margins here (so there are no margins) in order
; to get the proper aspect ratio for the circles.  Store the
; existing margins and restore them after rendering the circles.
xmargin = !x.margin
ymargin = !y.margin
!x.margin = [0.0,0.0]
!y.margin = [0.0,0.0]
wset,(*pState).winPix
;
factor=5
if ((*pState).autoscale eq 1) then begin
  (*pState).ylo = -radius[ncircles-1]*factor
  (*pState).yhi =  radius[ncircles-1]*factor
  (*pState).xlo = -radius[ncircles-1]*factor
  (*pState).xhi =  radius[ncircles-1]*factor
endif else begin
	xmid=0.5*((*pState).xlo+(*pState).xhi)
	ymid=0.5*((*pState).ylo+(*pState).yhi)
	xhwid=0.5*((*pState).xhi-(*pState).xlo)
	yhwid=0.5*((*pState).yhi-(*pState).ylo)
	if (xhwid lt yhwid) then begin
		(*pState).xlo=xmid-yhwid
		(*pState).xhi=xmid+yhwid
	endif else begin
		(*pState).ylo=ymid-xhwid
		(*pState).yhi=ymid+xhwid
	endelse
endelse
;
xmid=0.5*((*pState).xlo+(*pState).xhi)
ymid=0.5*((*pState).ylo+(*pState).yhi)
hwid=0.5*((*pState).xhi-(*pState).xlo)
aspectratio=float((*pState).winysize)/float((*pState).winxsize)
if (aspectratio lt 1.0) then begin
	(*pState).xlo=xmid-hwid/aspectratio
	(*pState).xhi=xmid+hwid/aspectratio
endif else begin
	(*pState).ylo=ymid-hwid*aspectratio
	(*pState).yhi=ymid+hwid*aspectratio
endelse
;
; Plot the outermost circle.
plot,xc[ncircles-1,*],yc[ncircles-1,*],psym = 0,xstyle = 5,ystyle = 5, color=1, background=0,$
	xrange = [(*pState).xlo,(*pState).xhi],yrange = [(*pState).ylo,(*pState).yhi]
polyfill,xc[ncircles-1,*],yc[ncircles-1,*],/line_fill,orientation=45+(ncircles-1)*90,$
	color=1+ncircles
oplot,xc[ncircles-1,*],yc[ncircles-1,*], color=1, thick=2
;
; Plot additional circles if any.
if (ncircles gt 1) then begin
	for i = ncircles-2,0,-1 do begin
		polyfill,xc[i,*],yc[i,*],color=0
		if ((*pState).sigt_in[i] gt 0.0) then $
			polyfill,xc[i,*],yc[i,*],/line_fill,orientation=45+i*90,$
			color=2+i
  	oplot,xc[i,*],yc[i,*],psym = 0, color=1, thick=2
	endfor
endif
;
;	Plot the incident beam.
xmin=!x.crange[0]
xmax=!x.crange[1]
ymin=!y.crange[0]
ymax=!y.crange[1]
oplot,[xmin,xmax],((*pState).beam_ydispl+0.5*(*pState).beam_wid)*[1,1],psym = 0, color=1, thick=1
oplot,[xmin,xmax],(*pState).beam_ydispl*[1,1],psym = 0, color=1, thick=1
oplot,[xmin,xmax],((*pState).beam_ydispl-0.5*(*pState).beam_wid)*[1,1],psym = 0, color=1, thick=1
;
; Plot directions to the detectors.
nphi=n_elements(*(*pState).phiPtr)
dis=sqrt((xmax-xmin)^2+(ymax-ymin)^2)*0.5
for k=0,nphi-1 do begin
	x=dis*cos((*(*pState).phiPtr)[k]*!dtor)
	y=dis*sin((*(*pState).phiPtr)[k]*!dtor)
	oplot,[0.0,x],[0.0,y],psym = 0, color=1, thick=1
endfor
;
wset,(*pState).winVis
device, copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0, (*pState).winPix]
;
!x.margin = xmargin
!y.margin = ymargin
(*pState).bangx=!x
(*pState).bangy=!y
(*pState).bangp=!p
;
widget_control,event.top,set_uvalue = pState
;
return
end


;************************************************************************************************
pro annsca_Legends_handler,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,get_uvalue=state
wdelete,state.winVis
widget_control,event.top,/destroy
end


;************************************************************************************************
pro annsca_Legends,event
;************************************************************************************************
;
compile_opt strictarr
;
bangx=!x
bangy=!y
bangp=!p
!x.margin = [0.0,0.0]
!y.margin = [0.0,0.0]
;
widget_control,event.top,get_uvalue = pState
;
tlb=widget_base(title="Explanation of curves",/col,group_leader=event.top,/modal)
		win=widget_draw(tlb,xsize=400,ysize=200)
		base=widget_base(tlb)
			done=widget_button(base,value="Dismiss")
centertlb,tlb
widget_control,tlb,/realize
;
widget_control,win,get_value=winVis
wset,winVis
;
xrange=[-0.5,2.5]
yrange=[-0.5,2.5]
x0=[0.25,0.4]
y0=[-0.5,-0.5]
dx=0.5
dy=0.25
plot,x0,y0,color=1,xrange=xrange,yrange=yrange,/nodata,xstyle=5,ystyle=5
nplots=(*pState).n_plots < (*pState).max_n_plots
for j=0,nplots-1 do begin
	ckn=(*(*pState).storeplotPtr[j]).ckn
	nann=ckn[0]
	labels="ssf"+strcompress((indgen(nann)+1),/remove_all)
	labels=[labels,"S"+strcompress((indgen(nann)+1),/remove_all)]
	if (ckn[2] gt ckn[1]) then begin
		for m=0,nann-1 do begin
			number=strcompress((m+1),/remove_all)
			labels=[labels,"S"+number+strcompress((indgen(nann)+1),/remove_all)]
		endfor
		labels=[labels,"S.S.","D.S."]
	endif else begin
		labels=[labels,"S.S."]
	endelse
	kv=annsca_get_kvals((*pState).plotwhat,ckn)
	for k=kv[0],kv[1] do begin
		xs=xrange[0]+j*dx
		ys=yrange[1]-(k-kv[0])*dy
		oplot,xs+x0,ys+y0,color=2+k,psym=-dave_selsym(j+5)
		xyouts,xs+x0[1]+0.1,ys+y0,labels[k],color=1
	endfor
	xyouts,xs+x0[0],yrange[1]-0.25,(*pState).explanation[j],color=1
endfor
;
state={winVis:winVis}
widget_control,tlb,set_uvalue=state
;
xmanager,'annsca_Legends',tlb,event_handler='annsca_Legends_handler'
;
!x=bangx
!y=bangy
!p=bangp
;
end


;************************************************************************************************
pro annsca_WinDraw,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,get_uvalue = pState
if (not (*pState).plotting_geometry and n_elements(*(*pState).ssfPtr) eq 0) then return
;
!x=(*pState).bangx
!y=(*pState).bangy
!p=(*pState).bangp
;
wset,(*pState).winVis
;
case event.type of
0: begin  ; button press
	(*pState).mouse = event.press
	if (*pState).mouse eq 4 then begin
	  (*pState).autoscale = 1
			if ((*pState).plotting_geometry) then annsca_ShowGeometry,event
			if (not (*pState).plotting_geometry) then annsca_2plotter,event,/redraw
   	  (*pState).mouse = 0
   	  widget_control,event.top,set_uvalue = pState
	  return
	endif
	(*pState).xp1 = event.x
	(*pState).yp1 = event.y
	(*pState).autoscale = 0
   end
;
1: begin	; button release
	if ((*pState).mouse eq 1) then begin
		x = (*pState).xp1 < (*pState).xp2
		y = (*pState).yp1 < (*pState).yp2
		w = abs((*pState).xp1 - (*pState).xp2)
		h = abs((*pState).yp1 - (*pState).yp2)
		x2 = x+w
		y2 = y+h
		lc = convert_coord(x, y, /device, /to_data)
		uc = convert_coord(x2, y2, /device, /to_data)
		wset,(*pState).winVis
		(*pState).xlo = lc[0] & (*pState).xhi = uc[0]
		(*pState).ylo = lc[1] & (*pState).yhi = uc[1]
		if ((*pState).plotting_geometry) then annsca_ShowGeometry,event
		if (not (*pState).plotting_geometry) then annsca_2plotter,event,/redraw
	endif
	(*pState).mouse = 0
   end
;
2: begin	; mouse motion
    if ((*pState).mouse eq 1) then begin
    	(*pState).xp2 = event.x
    	(*pState).yp2 = event.y
    	xc = [(*pState).xp1, event.x, event.x, (*pState).xp1, (*pState).xp1]
    	yc = [(*pState).yp1, (*pState).yp1, event.y, event.y, (*pState).yp1]
    	wset, (*pState).winVis
	    device, copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0, (*pState).winPix]
    	plots, xc, yc, /device, color=1
    endif
    xy=convert_coord(event.x,event.y,/device,/to_data)
    widget_control,(*pState).butt_disp,set_value=$
    	'2theta ='+string(xy[0],format='(f6.1)')+'. Ordinate ='+string(xy[1],format='(f7.4)')
   end
else:
endcase
widget_control,event.top,set_uvalue = pState
return
end


;************************************************************************************************
pro annsca_WriteFactors,event
;************************************************************************************************
;
compile_opt strictarr
;
widget_control,event.top,get_uvalue = pState
;
if (n_elements(*(*pState).ssfPtr) eq 0) then begin
	res=dialog_message("Nothing found to save to a file.",/information)
	return
endif
;
cd,(*!dave_defaults).workDir
filename=dialog_pickfile(dialog_parent=event.top,filter='*.txt',/write,title='Output file')
if (strlowcase(strmid(filename,strlen(filename)-4)) ne '.txt') then filename=filename+'.txt'
openw,lun,filename,/get_lun
;
printf,lun,(*pState).prob_comment
printf,lun,"Number of annuli"
printf,lun,(*pState).nann
nann=(*pState).nann
printf,lun,"Radii and height"
printf,lun,((*pState).radius)[0:nann-1]
printf,lun,(*pState).height
printf,lun,"Cross sections S(in), A(in), S(out), A(out)"
printf,lun,((*pState).sigs_in)[0:nann-1]
printf,lun,((*pState).siga_in)[0:nann-1]
printf,lun,((*pState).sigs_out)[0:nann-1]
printf,lun,((*pState).siga_out)[0:nann-1]
printf,lun,"Beam width, height, y and z displacements"
printf,lun,(*pState).beam_wid,(*pState).beam_hgt,(*pState).beam_ydispl,(*pState).beam_zdispl
printf,lun,"Numbers of steps and neutrons, random number seed"
printf,lun,(*pState).ny,(*pState).nmc
printf,lun,(*pState).seed
printf,lun,"ORC flag, ORC x and y displacements, and ORC parameters"
printf,lun,(*pState).orcsel,(*pState).orc_X,(*pState).orc_Y
printf,lun,(*pState).orc_bs,(*pState).orc_ir,(*pState).orc_or,(*pState).orc_bt
printf,lun,"Calculation type"
printf,lun,(*pState).calctype
;
if ((*pState).calctype le 1) then begin
	thisFormat = '('+strtrim(string(1+4*nann),2)+'f10.2)'
endif else begin
	thisFormat = '('+strtrim(string(1+4*nann+2*nann^2),2)+'f10.2)'
endelse
;
nphi=n_elements(*(*pState).phiPtr)
printf,lun,"Angle, SSFs and errors, Scattered intensities and errors"
;
for i=0,nphi-1 do begin
  dataOut=(*(*pState).phiPtr)[i]
  dataOut=[dataOut,(*(*pState).ssfPtr)[0:nann-1,i]]
  dataOut=[dataOut,(*(*pState).esfPtr)[0:nann-1,i]]
  dataOut=[dataOut,(*(*pState).ttotalPtr)[0:nann-1,i]]
  dataOut=[dataOut,(*(*pState).etotalPtr)[0:nann-1,i]]
  if ((*pState).calctype gt 1) then begin
    dataOut=[dataOut,reform((*(*pState).ttotal2Ptr)[0:nann-1,0:nann-1,i],nann^2)]
    dataOut=[dataOut,reform((*(*pState).etotal2Ptr)[0:nann-1,0:nann-1,i],nann^2)]
	endif
  printf,lun,dataOut,format=thisFormat
endfor
;
free_lun,lun
;
return
end


;************************************************************************************************
pro annsca_procs
;************************************************************************************************
;
compile_opt strictarr
;
end
