;+
;  Name:
;       BATSE_SPEC_MODEL
;	
; PURPOSE: 
;	calculate the low energy cross section for 
;	BATSE spec detectors from materials and geometry
;	Only considers the front surface.  Validity limited
;	to energies less than 14.0 keV.  The detector is
;	a right cylinder with a depth of 7.195 cm and 6.515 cm radius
; CATEGORY:
;	BATSE
; CALLING SEQUENCE:
;	batse_spec_model, theta, edges, drm 
; INPUTS:
;
;	theta - angle in degrees le 90 degrees
; OPTIONAL INPUTS:
;	
; OUTPUTS:
;	drm
;	edges - energy in keV for input and output edges for drm
; OPTIONAL OUTPUTS:
;
; PROCEDURE:
;	Sets up to run analytic model in RESP_CALC.
;
; CALLS:
;	xsec, other_filters, resp_calc, edge_products
; COMMON BLOCKS:
;
; RESTRICTIONS:
;
; MODIFICATION HISTORY:
;	ras, 2-mar-95	
;-
;
pro batse_spec_model, theta, edges,drm

;common esidecom, eside, xyz_rad, theta_rad, cos_rad

;.....................
;Need to convert thicknesses in cm of window materials to gm/cm2 for resp_calc.pro
;Specify the element z and use the structure coef to obtain the conversion factors

  coefdata, coef
  conv_cm_2_gmcm2 = coef.at_den * coef.conv ;multiply cm * this factor to obtain gm/cm2
  elem_windows= [4,13, 14]
  welem = elem_windows* 0
  for i=0,n_elements(welem)-1 do welem(i)=(where(elem_windows(i) eq coef.z))(0)
  window_cfactors = conv_cm_2_gmcm2(welem)
;.....................

thetar=theta /!radeg


nflux   = 100 ;100 channels in model

barea   = !pi*3.81^2

nseg = 20.
xyz_rad= [sin( (findgen(nseg)+.5)/nseg *90./!radeg), $
	cos( (findgen(nseg)+.5)/nseg *90./!radeg), fltarr(nseg)]
xyz_rad = transpose(reform(xyz_rad,nseg,3))


rsource = [sin(theta/!radeg),0.0,cos(theta/!radeg)]
cos_rad = (rsource#xyz_rad)(*)
theta_rad = acos(cos_rad) *!radeg

;.............
;% CHKARG: recalling resp_calc.pro from memory
;---- Module: resp_calc.pro
;---- From:   /usr/users/schwartz/soft/spex/new/
;---> Call: 
;pro resp_calc, detector,area,func,func_par,d,z,gmcm,nflux,elo,ehi, $
;eloss_mat, pls_ht_mat, edges, smatrix,                     $
;nosigma=nosigma, theta=theta

pfwhm = [.5,.5,0,0]
;First calculate for the interior region on the top of the detector
resp_calc,/nosigma, theta=theta, $
 'NAI', barea, 'fwhm', pfwhm, 7.195, [4], [0.127*window_cfactors(0)], $
	nflux, 3., 100., eloss0, pls_ht, edges

edge_products, edges, mean=em, edges_2=edges

att_multi = other_filters(em,/multi)^(4./cos(thetar))
for i=0,nflux-1 do eloss0(*,i) = eloss0(*,i) * att_multi
drm = eloss0 

;Now calculate for the outer annulus
;The silicon compound acts like silicon with a density of 1gm/cm3
;with the pe xsection of 1/2 that much silicon.

ann_area = !pi*(6.515^2-3.81^2)
resp_calc,/nosigma, theta=theta, $
 'NAI', ann_area, 'fwhm', pfwhm, 7.195, $
	[13,14],[0.068*window_cfactors(1),0.15/2.], nflux, 3.,100., eloss0

drm = drm + eloss0
;Now loop over each slice of the detector on the side

;eside = fltarr(nflux,nseg)
;ww = findgen(nflux)*(nflux+1)

for i=0,nseg-1 do begin

   slice_area = 7.195*!pi*6.515/nseg

   resp_calc,/nosigma, theta=theta_rad(i), $
    'NAI', slice_area, 'fwhm', pfwhm, 6.515*2*xyz_rad(0,i), $
	[13,14],[0.13*window_cfactors(1),0.15/2.], nflux, 3., 100., eloss0
   ;eside(0,i) = eloss0(ww)
   drm = drm + eloss0
endfor 
end