; $Id: slab_ss_calc.pro,v 1.6 2005/09/13 20:48:52 jrdcopley 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.
;
;###############################################################################
;************************************************************************************************
; Program to calculate SAFs and SSFs for a plane slab sample in a not necessarily "thin"
; container, by John Copley
;
;	slab_ss_calc,sigmatS,sigmatC,psi0,tth,signconv,ssf=ssf,fcsc=fcsc,fcc=fcc,saf=saf,$
;		transm_s=transm_s,transm_c=transm_c,ts0=ts0,tc0=tc0
; sigmatS:				the dimensionless product sigmaT_S*t_S
; sigmatC:				the dimensionless product sigmaT_C*t_C
;	psi0, psi:	are angles in degrees
; tth:			  is the array that contains the diffraction angles where SSF etc are to be calculated
;     								(note that tth is a double precission array)
; signconv: is the sign convention for angles phi0 and tth
;
; ! Default values are created if parameters were not defined on input
;************************************************************************************************
pro slab_ss_calc,sigmatS,sigmatC,psi0,tth,signconv,ssf=ssf,fcsc=fcsc,fcc=fcc,saf=saf,$
	transm_s=transm_s,transm_c=transm_c,ts0=ts0,tc0=tc0
;************************************************************************************************
;
compile_opt strictarr
;
;
; sigmatS is the dimensionless product sigmaT_S*t_S
; sigmatC is the dimensionless product sigmaT_C*t_C
; the sample thickness is t_S and the thickness of each wall of the container is t_C so the total
; thickness is t_S+2*t_C and the transmission normal to the slab is exp(-sigmatS-2*sigmatC)
;	psi0 and psi are angles in degrees
;
; Create default values if parameters were not defined on input
if (n_elements(sigmatS) eq 0) then sigmatS=0.05
if (n_elements(sigmatC) eq 0) then sigmatC=0.00
if (n_elements(psi0) eq 0) then psi0=60.0
if (n_elements(tth) eq 0) then begin
	start=28
	stop=32
	step=0.005
	n=(stop-start)/step+1.01
	tth=(start+findgen(n)*step)
endif
;
sigmatS=double(sigmatS)
sigmatC=double(sigmatC)
psi0=double(psi0)
;
ntth=n_elements(tth)
ssf=dblarr(ntth)
fcsc=dblarr(ntth)
fcc=dblarr(ntth)
saf=dblarr(ntth)
transm_s=dblarr(ntth)
transm_c=dblarr(ntth)
;
; If sigmatS < 0 or sigmatC < 0 set ssf etc = -1 and return.
if (sigmatS lt 0.0 or sigmatC lt 0.0) then begin
	ssf=ssf-1.0
	fcsc=fcsc-1.0
	fcc=fcc-1.0
	saf=saf-1.0
	transm_s=transm_s-1.0
	transm_c=transm_c-1.0
	ts0=-1.0
	tc0=-1.0
	return
endif
;
; If sigmatS = 0 and sigmatC = 0 set ssf etc = 1 and return.
if (sigmatS eq 0.0 and sigmatC eq 0.0) then begin
	ssf=ssf+1.0
	fcsc=fcsc+1.0
	fcc=fcc+1.0
	saf=saf+1.0
	transm_s=transm_s+1.0
	transm_c=transm_c+1.0
	ts0=1.0
	tc0=1.0
	return
endif
;
; Require that -!pi/2 < psi0 (in radians) <= !pi/2
psi0_in=psi0*!dtor
while (psi0_in gt !pi/2) do psi0_in=psi0_in-!pi
while (psi0_in le -!pi/2) do psi0_in=psi0_in+!pi
;
cos_psi0=cos(psi0_in)
if (cos_psi0 eq 0.0) then return
sec_psi0=1.0/cos_psi0
TS0=exp(-sigmatS*sec_psi0)
TC0=exp(-sigmatC*sec_psi0)
;
;print,'sigmatS,sigmatC,psi0'
;print,sigmatS,sigmatC,psi0
; Construct psi (in radians).
if (signconv[0] eq signconv[1]) then psi=psi0_in-tth*!dtor else psi=psi0_in+tth*!dtor
;
; For all values of psi require that -!pi < psi <= !pi
repeat begin
	ind=where(psi gt !pi, count)
	if (count gt 0) then psi[ind]=psi[ind]-2*!pi
endrep until (count eq 0)
;
repeat begin
	ind=where(psi le -!pi, count)
	if (count gt 0) then psi[ind]=psi[ind]+2*!pi
endrep until (count eq 0)
;
cos_psi=cos(psi)
;
; Treat cos(psi) > 0 (transmission), with cos(psi) <> cos(psi0)
ind=where(cos_psi gt 0.0 and abs(cos_psi-cos_psi0) gt 0.00001,count)
;print,'Transmission count: ',count
if (count gt 0) then begin
	sec_psi=1.0/cos_psi[ind]
	transm_s[ind]=exp(-sigmatS*sec_psi)
	transm_c[ind]=exp(-sigmatC*sec_psi)
	delta_sec=sec_psi0-sec_psi
	ssf[ind]=TC0*transm_c[ind]*(transm_s[ind]-TS0)/(sigmatS*delta_sec)
	if (sigmatC gt 0.0) then begin
		fcsc[ind]=0.5*(TS0*TC0+transm_s[ind]*transm_c[ind])*(transm_c[ind]-TC0)/(sigmatC*delta_sec)
		fcc[ind]=0.5*(TC0+transm_c[ind])*(transm_c[ind]-TC0)/(sigmatC*delta_sec)
		saf[ind]=(TS0*TC0+transm_s[ind]*transm_c[ind])/(TC0+transm_c[ind])
	endif else begin
		fcsc[ind]=0.5*(TS0+transm_s[ind])
		fcc[ind]=1.0
		saf[ind]=fcsc[ind]
	endelse
endif
;
; Treat cos(psi) < 0 (reflection)
ind=where(cos_psi lt 0.0,count)
;print,'Reflection count: ',count
if (count gt 0) then begin
	sec_psi=1.0/cos_psi[ind]
	transm_s[ind]=exp(sigmatS*sec_psi)
	transm_c[ind]=exp(sigmatC*sec_psi)
	delta_sec=sec_psi0-sec_psi
	ssf[ind]=TC0*transm_c[ind]*(1.0-TS0*transm_s[ind])/(sigmatS*delta_sec)
	if (sigmatC gt 0.0) then begin
		fcsc[ind]=0.5*(1.0+TS0*TC0*transm_s[ind]*transm_c[ind])*(1.0-TC0*transm_c[ind])/(sigmatC*delta_sec)
		fcc[ind]=0.5*(1.0+TC0*transm_c[ind])*(1.0-TC0*transm_c[ind])/(sigmatC*delta_sec)
		saf[ind]=(1.0+TS0*TC0*transm_s[ind]*transm_c[ind])/(1.0+TC0*transm_c[ind])
	endif else begin
		fcsc[ind]=0.5*(1.0+TS0*transm_s[ind])
		fcc[ind]=1.0
		saf[ind]=fcsc[ind]
	endelse
endif
;
; Treat cos(psi) > 0 (transmission), with cos(psi) = cos(psi0)
ind=where(cos_psi gt 0.0 and abs(cos_psi-cos_psi0) le 0.00001,count)
;print,'Transmission parallel count: ',count
if (count gt 0) then begin
	ssf[ind]=TC0*TS0*TC0
	fcsc[ind]=TC0*TS0*TC0
	fcc[ind]=TC0*TC0
	saf[ind]=TS0
	transm_s[ind]=TS0
	transm_c[ind]=TC0
endif
;
end