; $Id: annsca_calc.pro,v 1.6 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_ninsect,nann,rann,x0,y0,nphi,u,v,dism,disp,nint
;************************************************************************************************
;
compile_opt strictarr
;
; Calculates the intersections of nphi directed lines with a set of
; concentric circles centered at the origin.  The number of circles
; is nann and their radii, in increasing order, are (rann(i),i=1,nann).
; The lines go through the point (x0,y0) and their directions are (u,v).
; The calculated distances to intersection with circle k are dism(k),
; disp(k), such that dism(k).le.disp(k).  If there is no intersection
; both quantities are set equal to the same value.
; The numbers of intersections are nint.
;
front = (-x0*u-y0*v)##(intarr(nann)+1)
disc0 = front^2-x0^2-y0^2
disc = (disc0+(intarr(nphi)+1)##(rann^2)) > 0.0
root = sqrt(disc)
dism = front-root
disp = front+root
ndim = size(disc,/n_dimensions)
nint = 2*round(total(disc gt 0.0,ndim))
return
end


;************************************************************************************************
pro annsca_nout_calc,nann,dism,disp,sigt_out,nphi,arg_out
;************************************************************************************************
;
compile_opt strictarr
;
; Calculates the attenuation from scattering point to all detectors.
;
dm=dism>0
dp=disp>0
sg=sigt_out[0:nann-1]#(intarr(nphi)+1)
s1=total((dp-dm)*sg,1)
ssg=shift(sg,-1)
s2=total((dp-dm)*ssg,1)
s3=(dp[nann-1,*]-dm[nann-1,*])*ssg[nann-1]
arg_out=s1-s2+s3
return
end


;************************************************************************************************
pro annsca_orc_calc,x,y,cphi,sphi,bmax,tmax,del,orcfac
;************************************************************************************************
;
compile_opt strictarr
;
; Calculates the visibility of detectors looking through ORC.
;
if (bmax eq 0.0) then orcfac=1.0 else begin
	b = abs(cphi*y - sphi*x)
	orcfac = $
	(b lt del)*tmax*(1.0-del/bmax) + (b gt del)*(b lt bmax)*tmax*(1.0-b/bmax)
endelse
return
end


;************************************************************************************************
pro annsca_calc_area,rann,d,area
;************************************************************************************************
;
compile_opt strictarr
;
;	Calculates illuminated areas between y=0 and y=d.
;
nann=n_elements(rann)
area=fltarr(nann)
if (d eq 0.0) then return
sd=round(d/abs(d))
ad=abs(d)
for k=0,nann-1 do begin
	if (ad gt rann[k]) then begin
		area[k]=0.5*!pi*rann[k]^2
	endif else begin
		cth=ad/rann[k]
		sth=sqrt(1.0-cth^2)
		theta=acos(cth)
		area[k]=(0.5*!pi-theta+cth*sth)*rann[k]^2
	endelse
	if (k gt 0) then area[k]=area[k]-total(area[0:k-1])
endfor
area=sd*area
end


;************************************************************************************************
pro annsca_Calc_numerical,$
	event,$
	y_wid,y_max,y_min,$
	nann,radius,$
	sigs_in,sigs_out,$
	sigt_in,sigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans
;************************************************************************************************
;
compile_opt strictarr
;
;
; Within this routine the arrays radius,sigs_in,sigs_out,sigt_in and sigt_out are dimensioned nann.
;
widget_control,event.top,get_uvalue = pState
;
trans=0.0
path=fltarr(nann)
tot0=fltarr(nann)
*(*pState).ttotalPtr=fltarr(nann,nphi)
*(*pState).etotalPtr=fltarr(nann,nphi)
*(*pState).ssfPtr=fltarr(nann,nphi)
*(*pState).esfPtr=fltarr(nann,nphi)
;
ptr_free,(*pState).ttotal2Ptr
ptr_free,(*pState).etotal2Ptr
(*pState).ttotal2Ptr=ptr_new(/allocate_heap)
(*pState).etotal2Ptr=ptr_new(/allocate_heap)
;
ny=(*pState).ny
y = fltarr(ny)
dy = y_wid/ny
y = y_min+0.5*dy+dy*dindgen(ny)
;
; Main loop over y-coordinate
t1=systime(1)
divy=20
nyd=ny/divy
if (nyd eq 0) then nyd=1
for iy = 0,ny-1 do begin
	if (iy/nyd eq float(iy)/nyd) then begin
		message="Numerical calculation is"+strcompress(fix(100.0*iy/ny))+"% complete"
		widget_control,(*pState).message_area,set_value=message
	endif
;
; Calculate distances from the origin to inter-region interfaces
  annsca_ninsect,nann,radius,0.d,y[iy],1,1.d,0.d,dism_in,disp_in,nint
;
; Calculate contribution to transmission for this value of y
 	arg = 0.d
  for k = 0,nann-1 do begin
    if (k eq 0) then begin
      path[0] = disp_in[0]-dism_in[0]
      cumpath = path[0]
    endif else begin
      path[k] = disp_in[k] - dism_in[k] - cumpath
      cumpath = cumpath+path[k]
    endelse
    arg = arg + sigt_in[k]*path[k]
  endfor
  trans = trans + exp(-arg)
;
; Calculation for scattering points within each region
  for k = 0,nann-1 do begin
;
; Only do calculation if region has sigt_in>0 and neutron actually enters region
    if (sigt_in[k]*path[k] gt 0.0) then begin
;
; Decide step in x
 	    nx = path[k]/dy
      if (nx/2*2 ne nx) then nx = nx + 1
      if (nx lt 40) then nx = 40
      dx = path[k]/nx
      sum = 0.0*dindgen(nphi)
;
; First half of path (x<0)
      xmin = dism_in[k]
 	    xmax = 0.d
   	  if (k gt 0) then xmax = dism_in[k-1]
     	nx = fix((xmax-xmin)/dx)
      x = fltarr(nx)
      x = xmin+0.5*dx+dx*dindgen(nx)
; Sum over x
      for ix = 0,nx-1 do begin
;
; Calculate arg_in, i.e. sum of sigt_in*distance over regions, before scattering
 	      arg_in = 0.d
        for kp = k,nann-1 do begin
          if (kp eq k) then dist_in = x[ix] - dism_in[kp]
          if (kp gt k) then dist_in = 0.5*path[kp]
          arg_in = arg_in+sigt_in[kp]*dist_in
        endfor
;
; For each scattering angle calculate contribution to SSF
        annsca_ninsect,nann,radius[0:nann-1],x[ix],y[iy],nphi,cphi,sphi,dism_out,disp_out,nint
				annsca_nout_calc,nann,dism_out,disp_out,sigt_out,nphi,arg_out
				annsca_orc_calc,x[ix]-(*pState).orc_X,y[iy]-(*pState).orc_Y,cphi,sphi,bmax,tmax,del,orcfac
				sum=sum+exp(-arg_in-arg_out)*orcfac
   	  endfor
;
; Second half of path (x>0)
     	xmin = 0.0d
      if (k gt 0) then xmin = disp_in[k-1]
      xmax = disp_in[k]
      nx = fix((xmax-xmin)/dx)
      x = fltarr(nx)
      x = xmin+0.5*dx+dx*dindgen(nx)
;
; Sum over x
      for ix = 0,nx-1 do begin
;
; Calculate arg_in, i.e. sum of sigt_in*distance over regions, before scattering
        arg_in = 0.d
        for kp = 0,nann-1 do begin
          if (kp lt k) then dist_in = path[kp]
          if (kp eq k) then dist_in = 0.5*path[kp]+x[ix]-xmin
          if (kp gt k) then dist_in = 0.5*path[kp]
          arg_in = arg_in+sigt_in[kp]*dist_in
        endfor
;
; For each scattering angle calculate contribution to SSF
        annsca_ninsect,nann,radius[0:nann-1],x[ix],y[iy],nphi,cphi,sphi,dism_out,disp_out,nint
				annsca_nout_calc,nann,dism_out,disp_out,sigt_out,nphi,arg_out
				annsca_orc_calc,x[ix]-(*pState).orc_X,y[iy]-(*pState).orc_Y,cphi,sphi,bmax,tmax,del,orcfac
				sum=sum+exp(-arg_in-arg_out)*orcfac
      endfor
;
; Update totals
 	    for j = 0,nphi-1 do begin
        (*(*pState).ttotalPtr)[k,j]=(*(*pState).ttotalPtr)[k,j]+sum[j]*dx
      endfor
    endif
 	  tot0[k] = tot0[k] + path[k]
 	endfor
endfor
;
t2=systime(1)
;
(*pState).message='Computation time '+string(t2-t1,format='(f10.3)')+' seconds'
;
trans = trans/ny
;
annsca_calc_area,radius[0:nann-1],y_max,area1
annsca_calc_area,radius[0:nann-1],y_min,area2
area=area1-area2
;
*(*pState).ttotalPtr=$
	*(*pState).ttotalPtr*(dy/y_wid)*(sigt_in#(fltarr(nphi)+1))
;
tot0=tot0*dy
;
ind=where(sigt_in gt 0.d0)
(*(*pState).ttotalPtr)[ind,*]=(*(*pState).ttotalPtr)[ind,*]*((sigs_in[ind]/sigt_in[ind])#(fltarr(nphi)+1))
;
ind=where(sigs_in gt 0.d0)
sig_teff=sigs_in*(tot0/y_wid)#(fltarr(nphi)+1)
(*(*pState).ssfPtr)[ind,*]=(*(*pState).ttotalPtr)[ind,*]/sig_teff[ind,*]
;
end


;************************************************************************************************
pro annsca_Calc_MonteCarlo,$
	event,$
	y_wid,y_max,y_min,$
	nann,radius,$
	sigs_in,sigs_out,$
	sigt_in,sigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans
;************************************************************************************************
;
compile_opt strictarr
;
;
; Within this routine the arrays radius,sigs_in,sigs_out,sigt_in and sigt_out are dimensioned nann.
;
widget_control,event.top,get_uvalue = pState
;
trans=0.0
path=fltarr(nann)
*(*pState).ttotalPtr=fltarr(nann,nphi)
*(*pState).etotalPtr=fltarr(nann,nphi)
*(*pState).ssfPtr=fltarr(nann,nphi)
*(*pState).esfPtr=fltarr(nann,nphi)
;
ptr_free,(*pState).ttotal2Ptr
ptr_free,(*pState).etotal2Ptr
(*pState).ttotal2Ptr=ptr_new(/allocate_heap)
(*pState).etotal2Ptr=ptr_new(/allocate_heap)
;
attenh=fltarr(nann)
area=fltarr(nann)
prob1=fltarr(nann)
prob2=fltarr(nann)
;
nmc=(*pState).nmc
seed=(*pState).seed
if (strmatch(seed,'[0-9]*') eq 0) then seed=0 else seed=fix(seed)
if (seed eq 0) then begin
	prand=ptr_new(randomu(unknown,nmc,1+2*nann))
endif else begin
	prand=ptr_new(randomu(seed,nmc,1+2*nann))
endelse
;
t1=systime(1)
divmc=20
nmcd=nmc/divmc
if (nmcd eq 0) then nmcd=1
for imc = 0,nmc-1 do begin
	if (imc/nmcd eq float(imc)/nmcd) then begin
		message="Monte Carlo Calculation is"+strcompress(fix(100.0*imc/nmc))+"% complete"
		widget_control,(*pState).message_area,set_value=message
	endif
;
	y=y_min+y_wid*(*prand)[imc,0]
;
; Calculate distances from the origin to inter-region interfaces
 	annsca_ninsect,nann,radius,0.d,y,1,1.d,0.d,dism_in,disp_in,nint
;
; Calculate contribution to transmission for this value of y
  for k = 0,nann-1 do begin
 	  if (k eq 0) then begin
   	  path[0] = disp_in[0]-dism_in[0]
      cumpath = path[0]
    endif else begin
      path[k] = disp_in[k] - dism_in[k] - cumpath
      cumpath = cumpath+path[k]
    endelse
    attenh[k] = exp(-sigt_in[k]*path[k]/2.d)
  endfor
;
	prob_thru=1.d0
;
	for k=nann-1,0,-1 do begin
		prob1[k]=prob_thru*(1.d0-attenh[k])
		prob_thru=prob_thru*attenh[k]
	endfor
;
	for k=0,nann-1 do begin
		prob2[k]=prob_thru*(1.d0-attenh[k])
		prob_thru=prob_thru*attenh[k]
	endfor
;
	trans=trans+prob_thru
;
; Calculation for scattering points within each region
  for k = 0,nann-1 do begin
;
; Only do calculation if region has sigt_in>0 and neutron actually enters region
    if (sigt_in[k]*path[k] gt 0.0) then begin
;
; Pick a distance between 0 and 0.5*path[k]
			dis=-alog(1.d0-(*prand)[imc,k+1]*(1.d0-attenh[k]))/sigt_in[k]
;
; Decide between the two parts of the trajectory in region k. hence x.
			prob12=prob1[k]+prob2[k]
			if ((*prand)[imc,k+nann+1] lt prob1[k]/prob12) then begin
				x=dism_in[k]+dis
			endif else begin
				if (k eq 0) then x=dis else x=disp_in[k-1]+dis
			endelse
;
; For each scattering angle calculate contribution to SSF
      annsca_ninsect,nann,radius[0:nann-1],x,y,nphi,cphi,sphi,dism_out,disp_out,nint
  		annsca_nout_calc,nann,dism_out,disp_out,sigt_out,nphi,arg_out
			annsca_orc_calc,x-(*pState).orc_X,y-(*pState).orc_Y,cphi,sphi,bmax,tmax,del,orcfac
			term=prob12*exp(-arg_out)*orcfac
			(*(*pState).ttotalPtr)[k,*]=(*(*pState).ttotalPtr)[k,*]+term
			(*(*pState).etotalPtr)[k,*]=(*(*pState).etotalPtr)[k,*]+term^2
   	endif
  endfor
endfor
;
ptr_free,prand
;
t2=systime(1)
;
(*pState).message=string((t2-t1)/nmc*1000,format='(f10.3)')+' secs per 1000 neutrons'
;
trans = trans/nmc
;
*(*pState).ttotalPtr=*(*pState).ttotalPtr/nmc; "<x>"
*(*pState).etotalPtr=*(*pState).etotalPtr/nmc; "<x^2>"
*(*pState).etotalPtr=sqrt((*(*pState).etotalPtr-*(*pState).ttotalPtr^2)/(nmc-1)); "std dev=sqrt((<x^2>-<x>^2)/(nmc-1))"
;
annsca_calc_area,radius[0:nann-1],y_max,area1
annsca_calc_area,radius[0:nann-1],y_min,area2
area=area1-area2
;
ind=where(sigt_in gt 0.d0)
(*(*pState).ttotalPtr)[ind,*]=(*(*pState).ttotalPtr)[ind,*]*((sigs_in[ind]/sigt_in[ind])#(fltarr(nphi)+1))
(*(*pState).etotalPtr)[ind,*]=(*(*pState).etotalPtr)[ind,*]*((sigs_in[ind]/sigt_in[ind])#(fltarr(nphi)+1))
;
ind=where(sigs_in gt 0.d0)
sig_teff=sigs_in*(area/y_wid)#(fltarr(nphi)+1)
(*(*pState).ssfPtr)[ind,*]=(*(*pState).ttotalPtr)[ind,*]/sig_teff[ind,*]
(*(*pState).esfPtr)[ind,*]=(*(*pState).etotalPtr)[ind,*]/sig_teff[ind,*]
;
end


;************************************************************************************************
pro annsca_Calc,event
;************************************************************************************************
;
compile_opt strictarr
;
;
widget_control,event.top,get_uvalue = pState
;
numerical=(*pState).calcType eq 0
montecarlo=(*pState).calcType eq 1
mc_methodA=(*pState).calcType eq 2
mc_methodB=(*pState).calcType eq 3
;
!except = 0
;
nann = (*pState).nann
;
;orcsel = 0 (1)	if the ORC is excluded from (included in) the calculation
orcsel = (*pState).orcsel
;
; scatsel = 0 (1) for elastic (inelastic) scattering
scatsel = (*pState).scatsel
;
nradius = (*pState).radius[0:nann-1]
height=(*pState).height
;
; sigt_out = vector of outgoing sigma values
nsigs_in = (*pState).sigs_in[0:nann-1]
nsigt_in = (*pState).sigt_in[0:nann-1]
if (scatsel eq 0) then begin
  nsigs_out = (*pState).sigs_in[0:nann-1]
  nsigt_out = (*pState).sigt_in[0:nann-1]
endif else begin
  nsigs_out = (*pState).sigs_out[0:nann-1]
  nsigt_out = (*pState).sigt_out[0:nann-1]
endelse
;
; Information about incident beam.
beam_wid  = (*pState).beam_wid
beam_ydispl  = (*pState).beam_ydispl
beam_hgt  = (*pState).beam_hgt
beam_zdispl  = (*pState).beam_zdispl
;
; Information about ORC location.
orc_X = (*pState).orc_X
orc_Y = (*pState).orc_Y
;
; Information about ORC dimensions.
orc_bs = (*pState).orc_bs
orc_ir = (*pState).orc_ir
orc_or = (*pState).orc_or
orc_bt = (*pState).orc_bt
;
bmax=0.0
tmax=0.0
twodel=0.0
if (orcsel) then begin
	radc=(orc_ir*orc_or)/(orc_or-orc_ir)
	two_alpha_radc=orc_bs*!dtor*radc
	bmax=two_alpha_radc-(orc_ir+orc_or)/(orc_ir*orc_or)*radc*0.5*twodel
	tmax=bmax/two_alpha_radc
endif
del=0.5*twodel
;
; Get ready for the calculation.
sam_wid = nradius[nann-1]*2.0
sam_hgt = height
phi = *(*pState).phiPtr
phi1 = min(phi)
phi2 = max(phi)
dphi = phi[1]-phi[0]
nphi = n_elements(phi)
;
cphi 		= fltarr(nphi)
sphi 		= fltarr(nphi)
dis_out 	= fltarr(nphi,3)
sum 		= fltarr(nphi)
ssf 		= fltarr(nann,nphi)
esf 		= fltarr(nann,nphi); (only used for MC)
;
dism_in 	= fltarr(nann)
disp_in 	= fltarr(nann)
dism_out 	= fltarr(nann)
disp_out 	= fltarr(nann)
attenh 		= fltarr(nann)
area		= fltarr(nann)
prob1		= fltarr(nann); (only used for MC)
prob2		= fltarr(nann); (only used for MC)
;
phird = fltarr(nphi)
phird = phi*!dtor
cphi = cos(phird)
sphi = sin(phird)
;
y_min = (beam_ydispl - beam_wid*0.5) > (-0.5*sam_wid)
y_max = (beam_ydispl + beam_wid*0.5) < (+0.5*sam_wid)
y_wid = y_max - y_min
;
z_min = (beam_zdispl - beam_hgt*0.5) > (-0.5*sam_hgt)
z_max = (beam_zdispl + beam_hgt*0.5) < (+0.5*sam_hgt)
z_hgt = z_max - z_min
;
message="Starting the calculation....Please wait."
widget_control,(*pState).message_area,set_value=message
(*pState).message=''
;
; Numerical calculation starts here.
if (numerical) then annsca_Calc_numerical,$
	event,$
	y_wid,y_max,y_min,$
	nann,nradius,$
	nsigs_in,nsigs_out,$
	nsigt_in,nsigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans
;
; Monte Carlo calculation starts here.
if (montecarlo) then annsca_Calc_MonteCarlo,$
	event,$
	y_wid,y_max,y_min,$
	nann,nradius,$
	nsigs_in,nsigs_out,$
	nsigt_in,nsigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans
;
if (mc_methodA) then annsca_Calc_MonteCarlo_new,$
	event,$
	y_wid,y_max,y_min,$
	z_hgt,z_max,z_min,$
	nann,nradius,height,$
	nsigs_in,nsigs_out,$
	nsigt_in,nsigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans,1,0
;
if (mc_methodB) then annsca_Calc_MonteCarlo_new,$
	event,$
	y_wid,y_max,y_min,$
	z_hgt,z_max,z_min,$
	nann,nradius,height,$
	nsigs_in,nsigs_out,$
	nsigt_in,nsigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	trans,0,1
;
s1=total(nradius*nsigt_in)         ; e.g. R0*S0 + R1*S1 + R2*S2 + R3*S3
ssg=shift(nsigt_in,-1)
s2=total(nradius*ssg)              ; e.g. R0*S1 + R1*S2 + R2*S3 + R3*S0
s3=(nradius[nann-1])*ssg[nann-1]   ; e.g. R3*S0
arg_out=s1-s2+s3                   ; e.g. R0*S0 + (R1-R0)*S1 + (R2-R1)*S2 + (R3-R2)*S3
transc=exp(-2.0*arg_out)
;
widget_control,(*pState).meantrans,set_value=$
	'T_average ='+strcompress(string(trans,format='(f8.4)'))
widget_control,(*pState).mintrans,set_value =$
	'T_diameter ='+strcompress(string(transc,format='(f8.4)'))
widget_control,(*pState).message_area,set_value=(*pState).message
;
ymax = 1.2*max(ssf)
ymin = 0.0
(*pState).ylo = 0.0
(*pState).yhi = ymax
(*pState).xlo = min(phi)
(*pState).xhi = max(phi)
;
*(*pState).phiPtr = phi
;
if ((*pState).n_plots lt (*pState).max_n_plots) then begin
	if (numerical) then calc="Num."
	if (montecarlo) then calc="MC"
	if (mc_methodA) then calc="MC-A"
	if (mc_methodB) then calc="MC-B"
	(*pState).explanation[(*pState).n_plots]=calc
endif
;
annsca_2plotter,event,redraw=0
;
widget_control,(*pState).outbase2,/sensitive
;
widget_control,event.top,set_uvalue = pState
;
return
end
