pro g_c,energy,light,dlight, light2=light2
;+
;  Name:
;       G_C
;c*************************************************************************
;     pro g_c,energy,light,dlight
;
;c    This is a subroutine to find the light output from the NaI crystal
;c    for a given energy. The light output is  determined by a spline fit
;c    to a series of points off a "best eye line"  through the energy 
;c    loss-light out graph (L/E vs. log10(E)) of the crystal for energies 
;c    less than 720 keV. For larger energies, L/E is set equal to the 
;c    value of the spline at 720 keV. The spline data is stored in scoeff.dat.
;c    Input:
;c       energy--energy incident on the crystal, assumed to be an
;c	           ascending vector in keV
;c    Output:
;c       light--Light emitted/energy normalized to 511 keV
;c       dlight--d(light)/d(energy)
;     Optional Output:
;	light2 - light output with severe dip above 34 keV, 
;		 The severe dip was removed from light and dlight so the
;		 expression could be inverted.
;c    written 10/16/91            last modified 12/6/91      LAF
;
;c***************************************************************************
;               RAS, 7-Jul-1997, changed PERM_DATA to SSWDB_BATSE
;-
      common spline_com, break,c0,c1,c2,c3
  
;read in ASCII file of SCOEFF.PRO, parameterization of light per unit energy
;loss (L/E) as a function of photon energy
 
if n_elements(break) eq 0 then begin
	sc = dblarr(5,58)
	openr,lu,/get,concat_dir(getenv('SSWDB_BATSE'),'scoeff.dat')
	readf,lu,sc
	free_lun,lu
	sc = transpose(sc)
	break = sc(*,4)
	c0    = sc(*,0)
	c1    = sc(*,1)
	c2    = sc(*,2)
	c3    = sc(*,3)
endif

      e=exp(1.)
      dloge_de=alog10(e)/energy

;c    find where the data belongs in the spline fit
      en=alog10(energy)
	n_en = n_elements(en) ;number of input energies
	in = intarr(n_en) ;
		
	test = sort([en, break]) ;order of combined energies
	
	in_current= 0 + n_en
	j= 0

	for i=0, 58 + n_en-1 do begin ;sort through the combined energies

      		if (test(i) lt n_en) then begin
			in(j) = in_current ;index referenced to break energies
			j = j+1 
		endif else in_current = test(i)
        endfor
	in = ((in -n_en) > 0) < 57

;	print,in

;c    Calculate L and dL/dE

      x=en-break(in)
      light=c0(in) + x*(c1(in) + x*0.5*(c2(in) + x*c3(in)/3.))
      dlight=(c1(in) + x*(c2(in) + (c3(in)*x*0.5)))*dloge_de
      light2 = light
      w2 = where( energy ge 33.17 and energy lt 38.0, nw2)
      if nw2 ge 1 then begin
	moredip = 1.-((1. -tanh( (energy(w2)-33.17 )/2. ))*.015)
        light2(w2) = light2(w2) * moredip
      endif
      return
      end