; $Id: asc.pro,v 1.11 2005/02/08 18:24:44 azuah Exp $
; asc.pro
;
; Calculates absorption (self-shielding) correction factors for annular sample
; geometries.  The algorithm used here is from John Copley's program
; originally written in Fortran called RINGS.F.  Either numerical integration
; or Monte-Carlo integration over the area of the sample can be performed.
;
; Written by R.M.Dimeo and John Copley (8/6/01)
;
;###############################################################################
; 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 or if the code in this file is
;  included in another product.
;
;###############################################################################
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro insect_old,nann,rann,x0,y0,u,v,dism,disp,nint
; Calculates the intersections of a directed line 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 line goes through the point (x0,y0) and its direction is (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 zero.  The number of intersections is
; nint.
; This is the old version.
;
dism = dblarr(nann) & disp = dblarr(nann)
front = -x0*u-y0*v
disc0 = front^2-x0^2-y0^2
dism = 0.0*dindgen(nann)
disp = 0.0*dindgen(nann)
nint = 0
for i = nann-1,0,-1 do begin
  disc = disc0+rann[i]^2
  if (disc lt 0.0) then return
  root = sqrt(disc)
  dism[i] = front-root
  disp[i] = front+root
  nint = nint+2
endfor
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro insect,nann,rann,x0,y0,u,v,dism,disp,nint
; Calculates the intersections of a directed line 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 line goes through the point (x0,y0) and its direction is (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 zero.  The number of intersections is
; nint.
;
front = -x0*u-y0*v
disc0 = front^2-x0^2-y0^2
disc = (disc0+rann[0:nann-1]^2) > 0.0
root = sqrt(disc)
dism = front-root
disp = front+root
ind = where(disc gt 0.0,nint)
nint = nint*2
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ninsect,nann,rann,x0,y0,nphi,u,v,dism,disp,nint
; 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
nint = 2*round(total(disc gt 0.0,2))
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro out_calc0,nann,dism,disp,sigma_out,arg_out,k
; calculates the attenuation from scattering point to detector
; This is the old version.
;
nneg = 0
npos = 0
whereneg1 = where(dism lt 0.0d,nneg1)
whereneg2 = where(disp lt 0.0d,nneg2)
nneg = nneg1 + nneg2
wherepos1 = where(dism gt 0.0d,npos1)
wherepos2 = where(disp gt 0.0d,npos2)
npos = npos1 + npos2

arg_out = 0.0d
if (npos le nneg) then begin
; if npos<=nneg, path from scattering point to detector goes outward
  for kp = k,nann-1 do begin
    if (kp eq k) then dist_out = disp[kp]
    if (kp gt k) then dist_out = disp[kp]-disp[kp-1]
    arg_out = arg_out+sigma_out[kp]*dist_out
  endfor
endif else begin
; if npos>nneg, path from scattering point to detector first goes inward
  for kp = 0,nann-1 do begin
    if kp lt k then begin
      if (kp eq 0) then dist_out = disp[kp]-dism[kp]
      ;if (kp gt 0) then dist_out = disp[kp]-dism[kp]-dist_out
      if (kp gt 0) then dist_out = disp[kp]-dism[kp]-(disp[kp-1]-dism[kp-1])
    endif
    if (kp eq k) then dist_out = dism[kp-1]+disp[kp]-disp[kp-1]
    if kp gt k then dist_out = disp[kp] - disp[kp-1]
    arg_out = arg_out + sigma_out[kp]*dist_out
  endfor
endelse
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro out_calc,nann,dism,disp,sigma_out,arg_out,k
; calculates the attenuation from scattering point to detector
dm=dism>0
dp=disp>0
s1=total((dp-dm)*sigma_out)
ssigma=shift(sigma_out,-1)
s2=total((dp-dm)*ssigma)
s3=(dp[nann-1]-dm[nann-1])*ssigma[nann-1]
arg_out=s1-s2+s3
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro nout_calc,nann,dism,disp,sigma_out,nphi,arg_out
; calculates the attenuation from scattering point to all detectors
dm=dism>0
dp=disp>0
sg=sigma_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 orc_calc,x,y,cphi,sphi,bmax,tmax,del,orcfac
; 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 calc_area,rann,d,area
nann=n_elements(rann)
area=fltarr(nann)
if (d eq 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 RMDAnnulusPlot,event
widget_control,event.top,get_uvalue = pState,/no_copy
widget_control,(*pState).selSlider,get_value = slideValue
slideValue = fix(slideValue[0])-1
widget_control,(*pState).r,set_value = ((*pState).rvals)[slideValue]
widget_control,(*pState).sigma,set_value = ((*pState).sigvals)[slideValue]
widget_control,(*pState).sigma_out,set_value = ((*pState).sigoutvals)[slideValue]

if not ptr_valid((*pState).ssfPtr) then begin
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif
widget_control,(*pState).selSlider,get_value = selSlider
selSlider = fix(selSlider[0])-1

; Define the plotting symbol..a filled circle
theta = findgen(30)/29.*360.
scale = 1.0
xsym = cos(theta * !dtor) * scale
ysym = sin(theta * !dtor) * scale
usersym, xsym,ysym,/fill

ssf = *(*pState).ssfPtr
phi = *(*pState).phiPtr
ssfSize = size(ssf)
nann = ssfSize[1]
if selSlider gt nann-1 then begin
  selSlider = nann-1
  widget_control,(*pState).selSlider,set_value = nann
endif
nphi = ssfSize[2]
if (*pState).autoscale eq 1 then begin
  (*pState).ylo = 0.0
  (*pState).yhi = 1.2*max(ssf[selSlider,*])
  (*pState).xlo = min(phi)
  (*pState).xhi = max(phi)
endif

if (*pState).calcType eq 0 then begin
  thisyTitle = '!3SSF ['+strtrim(string(selSlider+1),2)+']: Analytic'
endif else begin
  thisyTitle = '!3SSF ['+strtrim(string(selSlider+1),2)+']: Monte-Carlo'
endelse

wset,(*pState).winPix
plot,phi,ssf[selSlider,*],psym = -8,yrange = [(*pState).ylo,(*pState).yhi],$
     ystyle = 1,xrange = [(*pState).xlo,(*pState).xhi],xstyle = 1,$
     xtitle = '!32!4h!3 (degrees)',ytitle = thisyTitle,$
     title = '!3Self-Shielding Attenuation Factor'
if (*pState).calcType eq 1 then begin
  esf = *(*pState).esfPtr
  errplot,phi,ssf[selSlider,*]-esf[selSlider,*],ssf[selSlider,*]+esf[selSlider,*],width = 0.0
endif

wset,(*pState).winVis
device,copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0,(*pState).winPix]

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusPrint,event
widget_control,event.top,get_uvalue = pState,/no_copy
if not ptr_valid((*pState).ssfPtr) then begin
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif
widget_control,(*pState).selSlider,get_value = selSlider
selSlider = fix(selSlider[0])-1

; Define the plotting symbol..a filled circle
theta = findgen(30)/29.*360.
scale = 1.0
xsym = cos(theta * !dtor) * scale
ysym = sin(theta * !dtor) * scale
usersym, xsym,ysym,/fill

ssf = *(*pState).ssfPtr
phi = *(*pState).phiPtr
ssfSize = size(ssf)
nann = ssfSize[1]
if selSlider gt nann-1 then begin
  selSlider = nann-1
  widget_control,(*pState).selSlider,set_value = nann
endif
nphi = ssfSize[2]
if (*pState).autoscale eq 1 then begin
  (*pState).ylo = 0.0
  (*pState).yhi = 1.2*max(ssf[selSlider,*])
  (*pState).xlo = min(phi)
  (*pState).xhi = max(phi)
endif

plotxsize = 5.5 & plotysize = 4.5

; Send to the printer
filename = DIALOG_PICKFILE(FILTER = '*.ps')
filename = filename+'.ps'
SET_PLOT, 'PS'
DEVICE, FILENAME=filename,/inches,xsize=5.5,ysize=4.5

plot,phi,ssf[selSlider,*],psym = -8,yrange = [(*pState).ylo,(*pState).yhi],$
     ystyle = 1,xrange = [(*pState).xlo,(*pState).xhi],xstyle = 1,$
     xtitle = '!62!4h!6 (degrees)',$
     title = 'Self-Shielding Attenuation Factor'
DEVICE, /CLOSE_FILE
case !version.os of
  'Win32':  begin
         set_plot,'win'
            end
  'IRIX':  begin
         set_plot,'X'
            end
  'linux':  begin
         set_plot,'X'
            end
else:
endcase
!p.multi = 0
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro RMDCalcannulusmc,event
widget_control,event.top,get_uvalue = pState,/no_copy
; sigma_out = vector of outgoing sigma values

!except = 0
nann = (*pState).nann
widget_control,(*pState).nmc,get_value = nmc
nmc = fix(nmc[0])

(*pState).calcType = 1

widget_control,(*pState).orcsel,get_value = orcsel
orcsel = 1 - fix(orcsel[0])
;orcsel = 0 ; do not use orc in calculation
;orcsel = 1 ; use orc in calculation

widget_control,(*pState).scatsel,get_value = scatsel
scatsel = fix(scatsel[0])
; scatsel = 0:    elastic scattering
; scatsel = 1:    inelastic scattering

sigma = (*pState).sigvals
sigma_out = dblarr((*pState).num)

if scatsel eq 0 then begin
  sigma_out = (*pState).sigvals
endif else begin
  sigma_out = (*pState).sigoutvals
endelse

radius = (*pState).rvals
widget_control,(*pState).beam_wid,get_value = beam_wid
beam_wid = 2.0*double(beam_wid[0])
widget_control,(*pState).beam_displ,get_value = beam_displ
beam_displ = double(beam_displ[0])
widget_control,(*pState).orc_bs,get_value = orc_bs
orc_bs = double(orc_bs[0])
widget_control,(*pState).orc_ir,get_value = orc_ir
orc_ir = double(orc_ir[0])
widget_control,(*pState).orc_or,get_value = orc_or
orc_or = double(orc_or[0])
widget_control,(*pState).orc_bt,get_value = orc_bt
orc_bt = double(orc_bt[0])

bmax=0.0
tmax=0.0
twodel=0.0
if (orcsel) then begin
   if (orc_or le orc_ir) then begin
    void = dialog_message('Invalid ORC radii!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
   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

rprev = 0.d
for k = 0,nann-1 do begin
  if (sigma[k] lt 0.d) then begin
    void = dialog_message('Invalid sigma values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
  if (sigma_out[k] lt 0.d) then begin
    void = dialog_message('Invalid sigma_out values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
  if (radius[k] le rprev) then begin
    void = dialog_message('Invalid radius values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
endfor

wid = radius[nann-1]*2.0
hwid = wid*0.5
phi = *(*pState).phiPtr
phi1 = min(phi)
phi2 = max(phi)
dphi = phi[1]-phi[0]
nphi = n_elements(phi)

cphi     = dblarr(nphi)
sphi     = dblarr(nphi)
dis_out  = dblarr(nphi,3)
sum    = dblarr(nphi)
ttotal    = dblarr(nann,nphi)
etotal    = dblarr(nann,nphi)
ssf    = dblarr(nann,nphi)
esf    = dblarr(nann,nphi)

dism_in  = dblarr(nann)
disp_in  = dblarr(nann)
dism_out    = dblarr(nann)
disp_out    = dblarr(nann)
tot0     = dblarr(nann)
etot0      = dblarr(nann)
path     = dblarr(nann)
attenh    = dblarr(nann)
area   = dblarr(nann)
prob1    = dblarr(nann)
prob2    = dblarr(nann)

phird = dblarr(nphi)
phird = phi*!dtor
cphi = cos(phird)
sphi = sin(phird)

trans = 0.0

cumpath = 0.d
wset,(*pState).winPix
erase
!p.charsize = 1.0
xyouts,0.15,0.5,'Thinking....please wait',/normal
!p.charsize = 1.0
wset,(*pState).winVis
device,copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0,(*pState).winPix]

y_min = (beam_displ - beam_wid*0.5) > (-hwid)
y_max = (beam_displ + beam_wid*0.5) < (+hwid)
y_wid = y_max - y_min
if (y_min ge radius[nann-1] or y_max le -radius[nann-1]) then begin
  void = dialog_message('Beam completely misses sample!',/information)
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif

;prand=ptr_new(randomu(100l,nmc,1+2*nann))
prand=ptr_new(randomu(unknown,nmc,1+2*nann))

; Monte Carlo loop
for imc = 0,nmc-1 do begin
   y=y_min+y_wid*(*prand)[imc,0]

; calculate distances from the origin to inter-region interfaces
  insect,nann,radius,0.d,y,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(-sigma[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 sigma>0 and neutron actually enters region
    if sigma[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]))/sigma[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
      ninsect,nann,radius[0:nann-1],x,y,nphi,cphi,sphi,dism_out,disp_out,nint
      nout_calc,nann,dism_out,disp_out,sigma_out,nphi,arg_out
      orc_calc,x,y,cphi,sphi,bmax,tmax,del,orcfac
      term=prob12*exp(-arg_out)*orcfac
      ttotal[k,*] = ttotal[k,*] + term
      etotal[k,*] = etotal[k,*] + term^2
;      for j = 0,nphi-1 do begin
;        insect,nann,radius,x,y,cphi[j],sphi[j],dism_out,disp_out,nint
;        out_calc,nann,dism_out,disp_out,sigma_out,arg_out,k
;        ttotal[k,j] = ttotal[k,j] + prob12*exp(-arg_out)
;        etotal[k,j] = etotal[k,j] + (prob12*exp(-arg_out))^2
;      endfor

    endif
    tot0[k] = tot0[k] + path[k]
    etot0[k] = etot0[k] + path[k]^2
  endfor
endfor

ptr_free,prand

trans = trans/nmc
s1=total(radius*sigma)
ssg=shift(sigma,-1)
s2=total(radius*ssg)
s3=(radius[nann-1])*ssg[nann-1]
arg_out=s1-s2+s3
transc=exp(-2.0*arg_out)
;print,'Transmission = ',trans,transc
widget_control,(*pState).meantrans,set_value = trans
widget_control,(*pState).mintrans,set_value = transc

ttotal=ttotal/nmc; "<x>"
etotal=etotal/nmc; "<x^2>"
etotal=sqrt((etotal-ttotal^2)/(nmc-1)); "std dev=sqrt((<x^2>-<x>^2)/(nmc-1))"

for k = 0,nann-1 do begin
  for j = 0,nphi-1 do begin
    ttotal[k,j] = ttotal[k,j]*y_wid/sigma[k]
    etotal[k,j] = etotal[k,j]*y_wid/sigma[k]
  endfor
endfor

tot0=tot0/nmc; "<x>"
etot0=etot0/nmc; "<x^2>"
etot0=sqrt((etot0-tot0^2)/(nmc-1)); "std dev=sqrt((<x^2>-<x>^2)/(nmc-1))"
tot0=tot0*y_wid
etot0=etot0*y_wid

calc_area,radius[0:nann-1],y_max,area1
calc_area,radius[0:nann-1],y_min,area2
area=area1-area2
for k=0,nann-1 do begin
;  if (k eq 0) then area[k]=!pi*radius[k]^2 else begin
;   area[k]=!pi*(radius[k]^2-radius[k-1]^2)
;  endelse
   if (sigma[k]*path[k] gt 0.0) then print,radius[k],area[k],tot0[k],etot0[k]
endfor

for j = 0,nphi-1 do begin
  for k = 0,nann-1 do begin
    ssf[k,j] = 0.0
    esf[k,j] = 0.0
    if sigma[k] gt 0.d then begin
      ssf[k,j] = ttotal[k,j]/area[k]
      esf[k,j] = etotal[k,j]/area[k]
    endif
  endfor
endfor

ymax = 1.2*max(ssf)
ymin = 0.0
(*pState).ylo = 0.0
(*pState).yhi = ymax
(*pState).xlo = min(phi)
(*pState).xhi = max(phi)

(*pState).ssfPtr = ptr_new(ssf,/no_copy)
(*pState).phiPtr = ptr_new(phi,/no_copy)
(*pState).esfPtr = ptr_new(esf,/no_copy)

widget_control,event.top,set_uvalue = pState,/no_copy
RMDannulusPlot,event
widget_control,event.top,get_uvalue = pState,/no_copy

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro RMDCalcAnnulus,event
widget_control,event.top,get_uvalue = pState,/no_copy
; sigma_out = vector of outgoing sigma values

widget_control,(*pState).orcsel,get_value = orcsel
orcsel = 1 - fix(orcsel[0])
;orcsel = 0 ; do not use orc in calculation
;orcsel = 1 ; use orc in calculation

widget_control,(*pState).scatsel,get_value = scatsel
scatsel = fix(scatsel[0])
; scatsel = 0:    elastic scattering
; scatsel = 1:    inelastic scattering

nann = (*pState).nann
widget_control,(*pState).ny,get_value = ny
ny = fix(ny[0])

(*pState).calcType = 0

sigma = dblarr((*pState).num)
sigma = (*pState).sigvals
sigma_out = dblarr((*pState).num)
if scatsel eq 0 then begin
  sigma_out = (*pState).sigvals
endif else begin
  sigma_out = (*pState).sigoutvals
endelse
radius = dblarr((*pState).num)
radius = (*pState).rvals

widget_control,(*pState).beam_wid,get_value = beam_wid
beam_wid = 2.0*double(beam_wid[0])
widget_control,(*pState).beam_displ,get_value = beam_displ
beam_displ = double(beam_displ[0])
widget_control,(*pState).orc_bs,get_value = orc_bs
orc_bs = double(orc_bs[0])
widget_control,(*pState).orc_ir,get_value = orc_ir
orc_ir = double(orc_ir[0])
widget_control,(*pState).orc_or,get_value = orc_or
orc_or = double(orc_or[0])
widget_control,(*pState).orc_bt,get_value = orc_bt
orc_bt = double(orc_bt[0])

bmax=0.0
tmax=0.0
twodel=0.0
if (orcsel) then begin
   if (orc_or eq orc_ir) then begin
    void = dialog_message('Invalid ORC radii!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
   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

rprev = 0.d
for k = 0,nann-1 do begin
  if (sigma[k] lt 0.d) then begin
    void = dialog_message('Invalid sigma values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
  if (sigma_out[k] lt 0.d) then begin
    void = dialog_message('Invalid sigma_out values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
  if (radius[k] le rprev) then begin
    void = dialog_message('Invalid radius values!',/information)
    widget_control,event.top,set_uvalue = pState,/no_copy
    return
  endif
endfor

wid = radius[nann-1]*2.0
hwid = wid*0.5
phi = *(*pState).phiPtr
phi1 = min(phi)
phi2 = max(phi)
dphi = phi[1]-phi[0]
nphi = n_elements(phi)

cphi     = dblarr(nphi)
sphi     = dblarr(nphi)
dis_out  = dblarr(nphi,3)
sum    = dblarr(nphi)
ttotal    = dblarr(nann,nphi)
ssf    = dblarr(nann,nphi)

dism_in  = dblarr(nann)
disp_in  = dblarr(nann)
dism_out    = dblarr(nann)
disp_out    = dblarr(nann)
tot0     = dblarr(nann)
path     = dblarr(nann)
attenh    = dblarr(nann)
area   = dblarr(nann)

phird = dblarr(nphi)
phird = phi*!dtor
cphi = cos(phird)
sphi = sin(phird)

trans = 0.0

y_min = (beam_displ - beam_wid*0.5) > (-hwid)
y_max = (beam_displ + beam_wid*0.5) < (+hwid)
y_wid = y_max - y_min
if (y_min ge radius[nann-1] or y_max le -radius[nann-1]) then begin
  void = dialog_message('Beam completely misses sample!',/information)
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif

y = dblarr(ny)
dy = y_wid/ny
y = y_min+0.5*dy+dy*dindgen(ny)

cumpath = 0.d
wset,(*pState).winPix
erase
!p.charsize = 1.0
xyouts,0.15,0.5,'!3Thinking....please wait',/normal
!p.charsize = 1.0
wset,(*pState).winVis
device,copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0,(*pState).winPix]


; Main loop over y-coordinate
for iy = 0,ny-1 do begin

; calculate distances from the origin to inter-region interfaces
  insect,nann,radius,0.d,y[iy],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 + sigma[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 sigma>0 and neutron actually enters region
    if sigma[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 = dblarr(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 sigma*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+sigma[kp]*dist_in
        endfor

; for each scattering angle calculate contribution to SSF
        ninsect,nann,radius[0:nann-1],x[ix],y[iy],nphi,cphi,sphi,dism_out,disp_out,nint
         nout_calc,nann,dism_out,disp_out,sigma_out,nphi,arg_out
         orc_calc,x[ix],y[iy],cphi,sphi,bmax,tmax,del,orcfac
         sum=sum+exp(-arg_in-arg_out)*orcfac
;        for j = 0,nphi-1 do begin
;          insect,nann,radius,x[ix],y[iy],cphi[j],sphi[j],dism_out,disp_out,nint
;          out_calc,nann,dism_out[*,j],disp_out[*,j],sigma_out[0:nann-1],arg_out,k
;          sum[j] = sum[j] + exp(-arg_in-arg_out)
;        endfor
      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 = dblarr(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 sigma*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+sigma[kp]*dist_in
        endfor

; for each scattering angle calculate contribution to SSF
        ninsect,nann,radius[0:nann-1],x[ix],y[iy],nphi,cphi,sphi,dism_out,disp_out,nint
         nout_calc,nann,dism_out,disp_out,sigma_out,nphi,arg_out
         orc_calc,x[ix],y[iy],cphi,sphi,bmax,tmax,del,orcfac
         sum=sum+exp(-arg_in-arg_out)*orcfac
;        for j = 0,nphi-1 do begin
;          insect,nann,radius,x[ix],y[iy],cphi[j],sphi[j],dism_out,disp_out,nint
;          out_calc,nann,dism_out[*,j],disp_out[*,j],sigma_out[0:nann-1],arg_out,k
;          sum[j] = sum[j] + exp(-arg_in-arg_out)
;        endfor
      endfor

; update totals
      for j = 0,nphi-1 do begin
        ttotal[k,j] = ttotal[k,j]+sum[j]*dx
      endfor
    endif
    tot0[k] = tot0[k] + path[k]
  endfor
endfor

trans = trans/ny
s1=total(radius*sigma)
ssg=shift(sigma,-1)
s2=total(radius*ssg)
s3=(radius[nann-1])*ssg[nann-1]
arg_out=s1-s2+s3
transc=exp(-2.0*arg_out)
;print,'Transmission = ',trans,transc
widget_control,(*pState).meantrans,set_value = trans
widget_control,(*pState).mintrans,set_value = transc
for k = 0,nann-1 do begin
  for j = 0,nphi-1 do begin
    ttotal[k,j] = ttotal[k,j]*dy
  endfor
  tot0[k] = tot0[k]*dy
endfor

for j = 0,nphi-1 do begin
  for k = 0,nann-1 do begin
    ssf[k,j] = 0.0
    if sigma[k] gt 0.d then ssf[k,j] = ttotal[k,j]/tot0[k]
  endfor
endfor

ymax = 1.2*max(ssf)
ymin = 0.0
(*pState).ylo = 0.0
(*pState).yhi = ymax
(*pState).xlo = min(phi)
(*pState).xhi = max(phi)

(*pState).ssfPtr = ptr_new(ssf,/no_copy)
(*pState).phiPtr = ptr_new(phi,/no_copy)

widget_control,event.top,set_uvalue = pState,/no_copy
RMDAnnulusPlot,event
widget_control,event.top,get_uvalue = pState,/no_copy

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusWinDraw,event
widget_control,event.top,get_uvalue = pState,/no_copy
if not ptr_valid((*pState).ssfPtr) then begin
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif

case event.type of
0: begin  ; button press
   (*pState).mouse = event.press
   if (*pState).mouse eq 4 then begin
     (*pState).autoscale = 1
     widget_control,event.top,set_uvalue = pState,/no_copy
      RMDAnnulusPlot,event
     widget_control,event.top,get_uvalue = pState,/no_copy
        (*pState).mouse = 0
        widget_control,event.top,set_uvalue = pState,/no_copy
     return
   endif
   (*pState).xp1 = event.x
   (*pState).yp1 = event.y
   (*pState).autoscale = 0
   wset, (*pState).winVis
   device, copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0, (*pState).winPix]
   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]

     widget_control,event.top,set_uvalue = pState,/no_copy
     RMDAnnulusPlot,event
     widget_control,event.top,get_uvalue = pState,/no_copy
   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
    endif
   end
else:
endcase
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ascCleanup,tlb
widget_control,tlb,get_uvalue = pState,/no_copy
s = size((*pState).notifyIDs)
if s[0] eq 1 then count = 0 else count = s[2]-1
for j = 0,count do begin
  ascInfo = {ascEvent,$
                  ID:(*pState).notifyIDs[0,j],$
                  Top:(*pState).notifyIDs[1,j],$
                  Handler:0l,$
                  inPtr:ptr_new(*(*pState).inPtr)}
   if widget_info((*pState).notifyIDs[0,j],/valid_id) then begin $
      widget_control,(*pState).notifyIDs[0,j],send_event = ascInfo
   endif
endfor
wdelete,(*pState).winPix
free_ptr,(*pState).xPtr
free_ptr,(*pState).yPtr
free_ptr,(*pState).origPhiPtr
free_ptr,(*pState).esfPtr
free_ptr,(*pState).phiPtr
free_ptr,(*pState).ssfPtr
free_ptr,(*pState).x0Ptr
free_ptr,(*pState).y0Ptr
free_ptr,(*pState).x1Ptr
free_ptr,(*pState).y1Ptr
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusQuit,event
widget_control,event.top,get_uvalue = pState,/no_copy
s = size((*pState).notifyIDs)
if s[0] eq 1 then count = 0 else count = s[2]-1
for j = 0,count do begin
  ascInfo = {ascEvent,$
                  ID:(*pState).notifyIDs[0,j],$
                  Top:(*pState).notifyIDs[1,j],$
                  Handler:0l,$
                  inPtr:ptr_new(*(*pState).inPtr)}
   if widget_info((*pState).notifyIDs[0,j],/valid_id) then begin $
      widget_control,(*pState).notifyIDs[0,j],send_event = ascInfo
   endif
endfor
widget_control,event.top,set_uvalue = pState,/no_copy
widget_control,event.top,/destroy
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusHFBSAngles,event
widget_control,event.top,get_uvalue = pState,/no_copy
detangles = dblarr(16)
detangles = [14.46,20.98,27.08,32.31,36.0,43.75,51.5,59.25,67.0,$
             74.75,82.5,90.25,98.0,105.75,113.5,121.5]

(*pState).origPhiPtr = ptr_new(detangles)
widget_control,(*pState).thlo,set_value = min(detangles)
widget_control,(*pState).thhi,set_value = max(detangles)
widget_control,(*pState).dth,set_value = detangles[1]-detangles[0]
(*pState).phiPtr = ptr_new(detangles,/no_copy)
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusDCSAngles,event
widget_control,event.top,get_uvalue = pState,/no_copy
file = (*pState).dcsFile
ndets = 913
detangles = dblarr(ndets)

catch,errno
if (errno ne 0) then begin
  catch,/cancel
  void = dialog_message(/error,["Unknown error detected",!err_string])
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif
; display the hourglass symbol to denote we're doing some work here
widget_control,/hourglass
angle = 0.0
groupnum = 0
openr, lun, file, /get_lun
x1    = 0.0
y1    = 0.0
instr = ''
for i = 0,4 do begin
  readf,lun,instr
endfor
for i = 0,ndets-1 do begin
  readf,lun,groupnum,y1,y1,y1,y1,y1,y1,y1,y1,angle
  detangles[i] = angle
endfor
free_lun,lun
(*pState).origPhiPtr = ptr_new(detangles)
y = detangles[uniq(detangles,sort(detangles))]
widget_control,(*pState).thlo,set_value = min(y)
widget_control,(*pState).thhi,set_value = max(y)
widget_control,(*pState).dth,set_value = y[1]-y[0]
(*pState).phiPtr = ptr_new(y,/no_copy)
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro calcHydrogenSigma,event
widget_control,event.top,get_uvalue = pState,/no_copy
rho = 0.088    ; g/cc
sigma_s = 82.02   ; barns
sigma_a = 0.3326  ; barns
m = 2.0      ; g/mol
widget_control,(*pState).lambda,get_value = lam & lam = double(lam[0])
sigmatot = sigma_s+(1.798/lam)*sigma_a
cs = rho*sigmatot*0.602/m

widget_control,(*pState).selSlider,get_value = sliderValue
sliderValue = fix(sliderValue[0])-1
(*pState).sigvals[sliderValue] = cs
widget_control,(*pState).sigma,set_value = cs
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro calcCH3ISigma,event
widget_control,event.top,get_uvalue = pState,/no_copy
rho = 2.27              ; g/cc
sigma_s = 5.551+3.0*82.02+3.81     ; barns
sigma_a = 0.0035+3.0*0.3326+6.15 ; barns
m = 142.0               ; g/mol
widget_control,(*pState).lambda,get_value = lam & lam = double(lam[0])
sigmatot = sigma_s+(1.798/lam)*sigma_a
cs = rho*sigmatot*0.602/m
widget_control,(*pState).selSlider,get_value = sliderValue
sliderValue = fix(sliderValue[0])-1
(*pState).sigvals[sliderValue] = cs
widget_control,(*pState).sigma,set_value = cs
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro calcAlSigma,event
widget_control,event.top,get_uvalue = pState,/no_copy
rho = 2.70             ; g/cc
sigma_s = 1.503           ; barns
m = 27.0            ; g/mol

widget_control,(*pState).lambda,get_value = lam & lam = double(lam[0])
sigma_a = 0.231
sigmatot = sigma_s+(1.798/lam)*sigma_a
cs = rho*sigmatot*0.602/m
widget_control,(*pState).selSlider,get_value = sliderValue
sliderValue = fix(sliderValue[0])-1
(*pState).sigvals[sliderValue] = cs
widget_control,(*pState).sigma,set_value = cs

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro calcHe4Sigma,event
widget_control,event.top,get_uvalue = pState,/no_copy
rho = 1.45             ; g/cc
sigma_s = 1.34          ; barns
m = 4.0              ; g/mol

widget_control,(*pState).lambda,get_value = lam & lam = double(lam[0])
sigma_a = 0.0
sigmatot = sigma_s+(1.798/lam)*sigma_a
cs = rho*sigmatot*0.602/m
widget_control,(*pState).selSlider,get_value = sliderValue
sliderValue = fix(sliderValue[0])-1
(*pState).sigvals[sliderValue] = cs
widget_control,(*pState).sigma,set_value = cs
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusWriteFactors,event
widget_control,event.top,get_uvalue = pState,/no_copy
if not ptr_valid((*pState).ssfPtr) then begin
  widget_control,event.top,set_uvalue = pState,/no_copy
  return
endif

origangles = *(*pState).origPhiPtr
angles = *(*pState).phiPtr
ssf = *(*pState).ssfPtr
sfSize = size(ssf)
nphi = sfSize[2]
nann = sfSize[1]

norigangles = n_elements(origangles)
;ssfOut = dblarr(nann,norigangles)
nangles = n_elements(angles)
ssfOut = dblarr(nann,nangles)
for i = 0,nangles-1 do begin
;  whereok = where(origangles eq angles[i],count)
;  if count gt 0 then begin
;    for j = 0,count-1 do begin
;      ssfOut[0:nann-1,whereok[j]] = ssf[0:nann-1,i]
      ssfOut[0:nann-1,i] = ssf[0:nann-1,i]
;    endfor
;  endif
endfor
cd,(*!dave_defaults).workDir
filename = DIALOG_PICKFILE(dialog_parent = event.top,$
                FILTER = '*.txt',/write,title = 'Output file')
filename = filename+'.txt'
openw,lun,filename,/get_lun
thisFormat = '('+strtrim(string(2+nann),2)+'f15.5)'
;dataOut = fltarr(2+nann)
dataOut = fltarr(1+nann)

;for i = 0,norigangles-1 do begin
for i = 0,nPhi-1 do begin
;  dataOut[0:2+nann-1] = [i+1,origangles[i],ssfOut[0:nann-1,i]]
  dataOut[0] = angles[i]
  dataOut[1:nann] = ssfOut[0:nann-1,i]
;  dataOut[0:2+nann-1] = [angles[i],ssfOut[0:nann-1,i]]
;  printf,lun,dataOut[0:2+nann-1],format = thisFormat
  printf,lun,dataOut[0:nann],format = thisFormat
;  print,dataOut[*],format = thisFormat
endfor

free_lun,lun


widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusDefineAngles,event
widget_control,event.top,get_uvalue = pState,/no_copy
widget_control,(*pState).dth,get_value = dth & dth = double(dth[0])
widget_control,(*pState).thhi,get_value = thhi & thhi = double(thhi[0])
widget_control,(*pState).thlo,get_value = thlo & thlo = double(thlo[0])
nangles = 1+fix((thhi-thlo)/dth)
phi = dblarr(nangles)
phi = thlo+dth*dindgen(nangles)
(*pState).phiPtr = ptr_new(phi,/no_copy)
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro annulusShowGeometry,event
widget_control,event.top,get_uvalue = pState,/no_copy

; Define the circle
ncircles = 5
radius = (*pState).rvals;dblarr(ncircles)
;for i = 0,ncircles-1 do begin
;  widget_control,(*pState).r[i],get_value = rval
;  radius[i] = double(rval[0])
;endfor

;radius = dblarr(ncircles)
;radius = 1.0+dindgen(ncircles)
np = 100
theta = dindgen(np)/(np-1.)*360.
xc = dblarr(ncircles,np)
yc = dblarr(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

; 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
plot,xc[ncircles-1,*],yc[ncircles-1,*],psym = 0,xstyle = 4,ystyle = 4
xyouts,radius[ncircles-1],0.0,strtrim(string(ncircles),2),/data
for i = 0,ncircles-2 do begin
  oplot,xc[i,*],yc[i,*],psym = 0
  xyouts,radius[i],0.0,strtrim(string(i+1),2),/data
endfor
wset,(*pState).winVis
device, copy = [0,0,(*pState).winxsize,(*pState).winysize,0,0, (*pState).winPix]

!x.margin = xmargin
!y.margin = ymargin

widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro ascHelp,event
widget_control,event.top,get_uvalue = pState
;pdf_file = file_which('asc.pdf',/include_current_dir)
pdf_file = !DAVE_AUXILIARY_DIR+'asc.pdf'
void = launch_help(pdf_file,tlb = event.top)
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro asc_event,event
widget_control,event.top,get_uvalue = pState,/no_copy
case event.id of
(*pState).dth: begin
         widget_control,event.top,set_uvalue = pState,/no_copy
         annulusDefineAngles,event
         widget_control,event.top,get_uvalue = pState,/no_copy
         end
(*pState).thlo:   begin
         widget_control,event.top,set_uvalue = pState,/no_copy
         annulusDefineAngles,event
         widget_control,event.top,get_uvalue = pState,/no_copy
         end
(*pState).thhi:   begin
         widget_control,event.top,set_uvalue = pState,/no_copy
         annulusDefineAngles,event
         widget_control,event.top,get_uvalue = pState,/no_copy
         end
(*pState).r:   begin
           widget_control,(*pState).selSlider,get_value = slideValue
           slideValue = fix(slideValue[0])-1
           widget_control,(*pState).r,get_value = rnew
           (*pState).rvals[slideValue] = double(rnew[0])
         end
(*pState).sigma:begin
           widget_control,(*pState).selSlider,get_value = slideValue
           slideValue = fix(slideValue[0])-1
           widget_control,(*pState).sigma,get_value = signew
           (*pState).sigvals[slideValue] = double(signew[0])
         end
(*pState).sigma_out:begin
           widget_control,(*pState).selSlider,get_value = slideValue
           slideValue = fix(slideValue[0])-1
           widget_control,(*pState).sigma_out,get_value = signew
           (*pState).sigoutvals[slideValue] = double(signew[0])
         end
else:
endcase
widget_control,event.top,set_uvalue = pState,/no_copy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro asc,notifyIDs,inPtr,group_leader = group_leader
; Annular Self-shielding Corrections program
dcsFile = !DAVE_AUXILIARY_DIR+'det_pos_jcc2.dat'
;dcsFile = file_which('det_pos_jcc2.dat',/include_current_dir)

loadct,0,/silent

tlb = widget_base(/row,title = $
      'Self-Shielding Calculator for Annular and Cylindrical Geometries',$
      mbar = bar,tlb_frame_attr = 1,$
      group_leader = group_leader)
filemenu = WIDGET_BUTTON(bar, VALUE='File', /MENU)
xs = 225l
xs1 = xs - 10l
butbase = widget_base(tlb,/col,/base_align_right,xsize = xs)

butbase2 = widget_base(tlb,/col,/base_align_right,xsize = xs)
plotbase = widget_base(tlb,/col,/base_align_center)
buttonSize = 100
void = widget_button(filemenu,value = 'Write Correction factors out to a file',$
       event_pro = 'annulusWriteFactors',xsize = buttonSize)

changebase = widget_base(butbase,/col,/base_align_right,/frame,xsize = xs1)
void = widget_label(changebase,value = 'Scattering Type',frame = 0,/align_center)
scattype = ['Elastic', 'Inelastic']
scatsel = cw_bgroup(changebase, scattype, /row, /exclusive,$
         /no_release,$
         set_value=0,/return_index)
selSlider = widget_slider(changebase,title='Select Region',value = 1,$
             minimum = 1,maximum = 5,$
             event_pro='RMDAnnulusPlot')
num = 5
sigvals = dblarr(num)
sigvals = [1.0,2.0,0.0,0.0,0.0]
sigma = cw_field(changebase,title='sigma (cm-1)',value=sigvals[0],$
                    xsize=8,/string,/return_events)
sigoutvals = [1.0,2.0,0.0,0.0,0.0]
sigma_out = cw_field(changebase,title='sigma out (cm-1)',value=sigoutvals[0],$
                    xsize=8,/string,/return_events)

rvals = [1.4,1.5,1.6,1.7,1.8]
r = cw_field(changebase,title='radius (cm)',value=rvals[0],$
                    xsize=8,/string,/return_events)
anglebase = widget_base(butbase,/col,/base_align_right,/frame,xsize = xs1)
thlo = cw_field(anglebase,title='Lowest angle (deg)',value=0.0,$
                    xsize=8,/string,/return_events)
thhi = cw_field(anglebase,title='Highest angle (deg)',value=180.0,$
                    xsize=8,/string,/return_events)
dth = cw_field(anglebase,title='Angle step (deg)',value=10.0,$
                    xsize=8,/string,/return_events)

nangles = 1+fix((180.0)/10.0)
th = dblarr(nangles)
th = 0.0 + 10.0*dindgen(nangles)

void = widget_button(filemenu,value = 'Show geometry',$
       event_pro = 'annulusShowGeometry',xsize = buttonSize)
void = widget_button(filemenu,value = 'HFBS Detector Angles',$
       event_pro = 'annulusHFBSAngles',xsize = buttonSize)
void = widget_button(filemenu,value = 'DCS Detector Angles',$
       event_pro = 'annulusDCSAngles',xsize = buttonSize,$
       sensitive = 1)

if !version.release eq 5.4 then begin
  void = widget_button(filemenu,value = 'Print to PS file',$
         event_pro = 'annulusPrint',xsize = buttonSize,sensitive = 1)
endif else begin
  void = widget_button(filemenu,value = 'Print current factors',$
         event_pro = 'annulusPrint',xsize = buttonSize)

endelse
void = widget_button(filemenu,value = 'Help',$
       event_pro = 'ascHelp')
void = widget_button(filemenu,value = 'Quit',$
       event_pro = 'annulusQuit')

beambase = widget_base(butbase,/col,/base_align_right,/frame,xsize = xs1)
beam_wid = cw_field(beambase,title='Beam halfwidth (cm)',value=5.0,$
                    xsize=8,/string,/return_events)
beam_displ = cw_field(beambase,title='Beam displacement (cm)',value=0.0,$
                    xsize=8,/string,/return_events)


orcbase = widget_base(butbase2,/col,/base_align_right,/frame,xsize = xs1)
void = widget_label(orcbase,value = 'Oscillating Radial Collimator',frame = 0,/align_center)
orctype = ['ORC', 'No ORC']
orcsel = cw_bgroup(orcbase, orctype, /row, /exclusive,$
         /no_release,$
         frame=1,set_value=1,/return_index)
orc_bs = cw_field(orcbase,title='Blade separation (deg)',value=2.0,$
                    xsize=8,/string,/return_events)
orc_ir = cw_field(orcbase,title='Inner radius (cm)',value=20.0,$
                    xsize=8,/string,/return_events)
orc_or = cw_field(orcbase,title='Outer radius (cm)',value=30.0,$
                    xsize=8,/string,/return_events)
orc_bt = cw_field(orcbase,title='Blade thickness (cm)',value=0.036,$
                    xsize=8,/string,/return_events)

newbutbase = widget_base(butbase2,/col,/frame,xsize = xs1)
void = widget_label(newbutbase,value = 'Calculation Type',frame = 0,/align_center)
butbase4 = widget_base(newbutbase,/row)
butbase3 = widget_base(newbutbase,/row)

void = widget_button(butbase3,value = 'Analytic',$
       event_pro = 'RMDCalcAnnulus')
ny = cw_field(butbase4,/col,title='Int steps',value=30,$
                    xsize=5,/string,/return_events)


nmc = cw_field(butbase4,/col,title='MC steps',value=700,$
                    xsize=5,/string,/return_events)
void = widget_button(butbase3,value = 'Monte-Carlo',$
       event_pro = 'RMDCalcAnnulusmc')

winxsize = 400
winysize = 400
win = widget_draw(plotbase,xsize = winxsize,ysize = winysize,$
      /button_events,/motion_events,event_pro='annulusWinDraw')

transbase = widget_base(butbase2,/row,/base_align_right)
meantrans = cw_field(transbase,/col,title='T_average',value=0.0,$
                    xsize=8,/string,/noedit)
mintrans = cw_field(transbase,/col,title='T_diameter',value=0.0,$
                    xsize=8,/string,/noedit)

centertlb,tlb
widget_control,tlb,/realize

!p.charsize = 1.0
window,/free,/pixmap,xsize = winxsize,ysize = winysize
winPix = !d.window
widget_control,win,get_value = winVis

state = {winxsize:winxsize,$
      winysize:winysize,$
      meantrans:meantrans,$
      mintrans:mintrans,$
      orcsel:orcsel,$
      scatsel:scatsel,$
      sigoutvals:sigoutvals,$
      sigma_out:sigma_out,$
      xPtr:ptr_new(),$
      yPtr:ptr_new(),$
      beam_wid:beam_wid,$
      beam_displ:beam_displ,$
      orc_bs:orc_bs,$
      orc_ir:orc_ir,$
      orc_or:orc_or,$
      orc_bt:orc_bt,$
      ny:ny,$
      nmc:nmc,$
      nann:5,$
      sigma:sigma,$
      sigvals:sigvals,$
      dcsFile:dcsFile,$
      r:r,$
      rvals:rvals,$
      winVis:winVis,$
      xlo:0.0,$
      ylo:0.0,$
      xhi:1.0,$
      yhi:1.0,$
      mouse:0,$
      xp1:0.0,$
      xp2:0.0,$
      yp1:0.0,$
      yp2:0.0,$
      num:num,$
      thlo:thlo,$
      thhi:thhi,$
      autoscale:1,$
      calcType:0,$   ; 0 for analytical and 1 for Monte-Carlo
      dth:dth,$
      selSlider:selSlider,$
      origPhiPtr:ptr_new(th),$
      esfPtr:ptr_new(),$
      phiPtr:ptr_new(th),$
      ssfPtr:ptr_new(),$
      x0Ptr:ptr_new(),$
      y0Ptr:ptr_new(),$
      x1Ptr:ptr_new(),$
      y1Ptr:ptr_new(),$
      winPix:winPix,$
      inPtr:inPtr,$
      notifyIDs:notifyIDs}

widget_control,tlb,set_uvalue = ptr_new(state),/no_copy

xmanager,'asc',tlb,cleanup = 'ascCleanup'

return
end


