; $Id: annsca_calc_monteCarlo_new.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.)
;
;************************************************************************************************
function annsca_calc_path0,radsq,nann,x,y
;************************************************************************************************
;
; Calculates path lengths through a set of concentric circles centered at the origin.
; The number of circles is nann and their squared radii, in increasing order, are in radsq.
;	The neutron starts at (x,y) and is travelling in the +x direction.
; The quantity returned is a set of nann*2 path lengths.
;
compile_opt strictarr
;
;
dp=sqrt((radsq-y^2) > 0.0)
dmp=[-rotate(dp,5),0.0,dp]
dis=(dmp-x) > 0.0
path=[dis,0.0]-[0.0,dis]
return,path[1:nann*2]
end


;************************************************************************************************
function annsca_calc_path1,radsq,nann,x0,y0,u,v
;************************************************************************************************
;
; Calculates path lengths through a set of concentric circles centered at the origin.
; The number of circles is nann and their squared radii, in increasing order, are in radsq.
;	The neutron starts at (x0,y0) and is travelling in the (u,v) direction.
; The quantity returned is a set of nann*2 path lengths.
;
compile_opt strictarr
;
;
x=(x0*u+y0*v)
y=(x0*v-y0*u)
dp=sqrt((radsq-y^2) > 0.0)
dmp=[-rotate(dp,5),0.0,dp]
dis=(dmp-x) > 0.0
path=[dis,0.0]-[0.0,dis]
return,path[1:nann*2]
end


;************************************************************************************************
function annsca_calc_pathn,radsq,nann,x0,y0,u,v,nphi,unann,unann2,unphi1,unphit
;************************************************************************************************
;
; Calculates path lengths through a set of concentric circles centered at the origin.
; The number of circles is nann and their squared radii, in increasing order, are in radsq.
;	The neutron starts at (x0,y0) and calculations are performed for nphi pairs of (u,v) directions.
; The quantity returned is a 2-d array of path lengths, nphi rows, nann*2 columns.
;
compile_opt strictarr
;
;
x=(x0*u+y0*v)##unann
y=(x0*v-y0*u)##unann
dp=sqrt((unphi1##radsq-y^2) > 0.0)
dmp=[-rotate(dp,5),unphit,dp]
x2=(x0*u+y0*v)##unann2
dis=(dmp-x2) > 0.0
path=[dis,unphit]-[unphit,dis]
return,path[1:nann*2,*]
end



;************************************************************************************************
pro annsca_Calc_MonteCarlo_new,$
	event,$
	y_wid,y_max,y_min,$
	z_hgt,z_max,z_min,$
	nann,radius,height,$
	sigs_in,sigs_out,$
	sigt_in,sigt_out,$
	nphi,cphi,sphi,$
	bmax,tmax,del,orcfac,$
	tot_trans,$
	methodA,methodB
;************************************************************************************************
;
; Monte Carlo calculation of single and double scattering.
;
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
;
; Ratio of scattering to total cross sections.
scatrat=fltarr(nann)
ind=where(sigt_in gt 0.0)
scatrat[ind]=sigs_in[ind]/sigt_in[ind]
;
; Number of Monte Carlo neutrons.
nmc=(*pState).nmc
;
; Random number seed.
seed=(*pState).seed
if (strmatch(seed,'[0-9]*') eq 0) then seed=0 else seed=fix(seed)
;
; Generate some random numbers, for y and distance to scattering.
; For double scattering generate some more random numbers, for z,
; polar and azimuthal angles and distance to second scattering.
if (seed eq 0) then begin
	prand=ptr_new(randomu(unknown,nmc,2))
	prand2=ptr_new(randomu(unknown,nmc,4))
endif else begin
	prand=ptr_new(randomu(seed,nmc,2))
	prand2=ptr_new(randomu(seed,nmc,4))
endelse
;
; [Code for progress messages.]
divmc=20
nmcd=nmc/divmc
if (nmcd eq 0) then nmcd=1
;
; Construct mirrored vectors, e.g. if sigma is [sigma_1,sigma_2,sigma_3],
; then sigma_mirrored=[sigma_3,sigma_2,sigma_1,sigma_1,sigma_2,sigma_3]
sigt_in_mirrored=[rotate(sigt_in,5),sigt_in]
sigt_out_mirrored=[rotate(sigt_out,5),sigt_out]
scatrat_mirrored=[rotate(scatrat,5),scatrat]
;
radsq=radius^2
nann2=nann*2
;
tpi=2.0*!pi
;
; Unit vectors needed later.
unann=intarr(nann)+1
unann2=intarr(2*nann+1)+1
unphi1=intarr(nphi)+1
unphit=transpose(intarr(nphi))
;
; Subregion vector, e.g. [3,2,1,0,0,1,2,3]
subregion=[rotate(indgen(nann),5),indgen(nann)]
;
; Starting value of x (which must be on or outside the outermost annulus).
xstart=-radius[nann-1]
;
; Initalize various quantities.
tot_trans=0.0
*(*pState).ttotalPtr=fltarr(nann,nphi)
*(*pState).etotalPtr=fltarr(nann,nphi)
*(*pState).ssfPtr=fltarr(nann,nphi)
*(*pState).esfPtr=fltarr(nann,nphi)
;
*(*pState).ttotal2Ptr=fltarr(nann,nann,nphi)
*(*pState).etotal2Ptr=fltarr(nann,nann,nphi)
;
; Start a timer.
t1=systime(1)
;
; Main Monte Carlo loop.
for imc = 0,nmc-1 do begin
;
; [More code for progress messages.]
	if (imc/nmcd eq float(imc)/nmcd) then begin
		message="Calculation is"+strcompress(fix(100.0*imc/nmc))+"% complete"
		widget_control,(*pState).message_area,set_value=message
	endif
;
; Pick a value of y.
	y=y_min+y_wid*(*prand)[imc,0]
;
; Calculate distances to inter-region interfaces, hence path lengths ("t").
	tee=annsca_calc_path0(radsq,nann,xstart,y)
;
;	Cumulative path lengths ("L").
	ell=[0,total(tee,/cumulative)]
;
;	Attenuation lengths ("b").
 	bee=tee*sigt_in_mirrored
;
;	Cumulative attenuation lengths ("A").
 	uc_A2=total(bee,/cumulative)
	uc_A1=[0.0,uc_A2[0:nann2-2]]
;
;	Probability of escape.
 	trans=exp(-uc_A2[nann2-1])
;
	if (methodA) then begin
;
; Pick an attenuation distance between 0 and uc_A2[2*nann].
		q=-alog(1.d0-(*prand)[imc,1]*(1.d0-trans))
;
; Determine the subregion.
		k=where(q gt uc_A1 and q lt uc_A2, count)
		if (count ne 1) then begin
			res=dialog_message('PROBLEM - single, method A, imc = '+strcompress(imc))
			if (q le uc_A1[0]) then k=0
			if (q ge uc_A2[nann2-1]) then k=nann2-1
		endif
;
;	Distance traveled.
		k=k[0]
		dis=ell[k]+(q-uc_A1[k])/sigt_in_mirrored[k]
		x=-ell[nann]+dis
;
;	Statistical weight.
		wgt=(1.0-trans)*scatrat_mirrored[k]
;
	endif
;
	if (methodB) then begin
;
; Generate scattering probabilities ("P").
		pee=exp(-uc_A1)*(1.0-exp(-bee))*scatrat_mirrored
;
; Generate cumuative scattering probabilities ("c").
		cee2=total(pee,/cumulative)
		cee1=[0.0,cee2[0:nann2-2]]
;
; Determine the subregion.
		zee=(*prand)[imc,1]*cee2[nann2-1]
		k=where(zee gt cee1 and zee lt cee2, count)
		if (count ne 1) then begin
			res=dialog_message('PROBLEM - single, method B, imc = '+strcompress(imc))
			if (zee le cee1[0]) then k=0
			if (zee ge cee2[nann2-1]) then k=nann2-1
		endif
;
;	Distance traveled.
		k=k[0]
		dis=ell[k]-alog(1.d0-(zee-cee1[k])*exp(uc_A1[k])/scatrat_mirrored[k])/sigt_in_mirrored[k]
		x=-ell[nann]+dis
;
;	Statistical weight.
		wgt=cee2[nann2-1]
;
	endif
;
; For each scattering angle calculate contribution to SSF
	path=annsca_calc_pathn(radsq,nann,x,y,cphi,sphi,nphi,unann,unann2,unphi1,unphit)
;
	thing=(intarr(nphi)+1)##sigt_out_mirrored
	arg_out=total(thing*path,1)
;
	if (bmax eq 0.0) then begin
		orcfac=1.0
	endif else begin
		annsca_orc_calc,x-(*pState).orc_X,y-(*pState).orc_Y,cphi,sphi,bmax,tmax,del,orcfac
	endelse
;
;	Contribution to T0
 	tot_trans=tot_trans+trans
;
	term=wgt*exp(-arg_out)*orcfac
;
; Given the subregion k, determine the region j.
	j=subregion[k]
;
	(*(*pState).ttotalPtr)[j,*]=(*(*pState).ttotalPtr)[j,*]+term
	(*(*pState).etotalPtr)[j,*]=(*(*pState).etotalPtr)[j,*]+term^2
;
; START OF DOUBLE SCATTERING CALCULATION.
; Pick a value of z.
	z=z_min+z_hgt*(*prand2)[imc,0]
;
;	Now pick a new direction.
; First pick the cosine of a polar angle relative to the cylindrical axis.
	cospolar=1.0-2.0*(*prand2)[imc,1]
	sinpolar=sqrt(1.0-cospolar^2)
;
;	Next pick the sine and cosine of an azimuthal angle.
	azimu=(*prand2)[imc,2]*tpi
	cosazimu=cos(azimu)
	sinazimu=sin(azimu)
;
;	Hence direction cosines.
	u=sinpolar*cosazimu
	v=sinpolar*sinazimu
	w=cospolar
;
; Calculate distances to inter-region interfaces, hence path lengths.
	tee=annsca_calc_path1(radsq,nann,x,y,cosazimu,sinazimu)/sinpolar
;
;	Cumulative path lengths.
	ell=[0,total(tee,/cumulative)]
;
;	Determine the maximum cumulative path length.
	ellmax=(cospolar lt 0.0) ? ((-height*0.5-z)/cospolar) : $
		((cospolar gt 0.0) ? ((height*0.5-z)/cospolar) : (1000.0))
;
;	Truncate the cumulative path lengths and then recalculate path lengths.
	ell=ell < ellmax
	tee=[ell,0.0]-[0.0,ell]
	tee=tee[1:nann2]
;
;	Attenuation lengths.
 	bee=tee*sigt_in_mirrored
;
;	Cumulative attenuation lengths.
 	uc_A2=total(bee,/cumulative)
	uc_A1=[0.0,uc_A2[0:nann2-2]]
;
;	Probability of escape.
 	trans2=exp(-uc_A2[nann2-1])
;
	if (methodA) then begin
;
; Pick an attenuation distance between 0 and uc_A2[2*nann].
		q=-alog(1.d0-(*prand2)[imc,3]*(1.d0-trans2))
;
; Determine the subregion.
		k=where(q gt uc_A1 and q lt uc_A2, count)
		if (count ne 1) then begin
			res=dialog_message('PROBLEM - double, method A, imc = '+strcompress(imc))
			if (q le uc_A1[0]) then k=0
			if (q ge uc_A2[nann2-1]) then k=nann2-1
		endif
;
;	Distance traveled.
		k=k[0]
		dis=ell[k]+(q-uc_A1[k])/sigt_in_mirrored[k]
		x=x+u*dis
		y=y+v*dis
		z=z+w*dis
;
;	Statistical weight.
		wgt=wgt*(1.0-trans2)*scatrat_mirrored[k]
;
	endif
;
	if (methodB) then begin
;
; Generate scattering probabilities ("P").
		pee=exp(-uc_A1)*(1.0-exp(-bee))*scatrat_mirrored
;
; Generate cumuative scattering probabilities ("c").
		cee2=total(pee,/cumulative)
		cee1=[0.0,cee2[0:nann2-2]]
;
; Determine the subregion.
		zee=(*prand)[imc,1]*cee2[nann2-1]
		k=where(zee gt cee1 and zee lt cee2, count)
		if (count ne 1) then begin
			res=dialog_message('PROBLEM - single, method B, imc = '+strcompress(imc))
			if (zee le cee1[0]) then k=0
			if (zee ge cee2[nann2-1]) then k=nann2-1
		endif
;
;	Distance traveled.
		k=k[0]
		dis=ell[k]-alog(1.d0-(zee-cee1[k])*exp(uc_A1[k])/scatrat_mirrored[k])/sigt_in_mirrored[k]
		x=x+u*dis
		y=y+v*dis
		z=z+w*dis
;
;	Statistical weight.
		wgt=wgt*cee2[nann2-1]
;
	endif
;
; For each scattering angle calculate contribution to SSF
	path=annsca_calc_pathn(radsq,nann,x,y,cphi,sphi,nphi,unann,unann2,unphi1,unphit)
;
	thing=(intarr(nphi)+1)##sigt_out_mirrored
	arg_out=total(thing*path,1)
;
	if (bmax eq 0.0) then begin
		orcfac=1.0
	endif else begin
		annsca_orc_calc,x-(*pState).orc_X,y-(*pState).orc_Y,cphi,sphi,bmax,tmax,del,orcfac
	endelse
;
	term=wgt*exp(-arg_out)*orcfac
;
; Given the subregion k, determine the region j.
	j2=subregion[k]
;
	(*(*pState).ttotal2Ptr)[j,j2,*]=(*(*pState).ttotal2Ptr)[j,j2,*]+term
	(*(*pState).etotal2Ptr)[j,j2,*]=(*(*pState).etotal2Ptr)[j,j2,*]+term^2
endfor
;
t2=systime(1)
;
(*pState).message=string((t2-t1)/nmc*1000,format='(f10.3)')+' secs per 1000 neutrons'
;
ptr_free,prand
ptr_free,prand2
;
tot_trans=tot_trans/nmc
;
; Compute "<x>", "<x^2>", and "sqrt((<x^2>-<x>^2)/(nmc-1))" (single scattering)
*(*pState).ttotalPtr=*(*pState).ttotalPtr/nmc
*(*pState).etotalPtr=*(*pState).etotalPtr/nmc
*(*pState).etotalPtr=sqrt((*(*pState).etotalPtr-*(*pState).ttotalPtr^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(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,*]
;
; Compute "<x>", "<x^2>", and "sqrt((<x^2>-<x>^2)/(nmc-1))" (double scattering)
*(*pState).ttotal2Ptr=*(*pState).ttotal2Ptr/nmc
*(*pState).etotal2Ptr=*(*pState).etotal2Ptr/nmc
*(*pState).etotal2Ptr=sqrt((*(*pState).etotal2Ptr-*(*pState).ttotal2Ptr^2)/(nmc-1))
;
end


