subroutine airgap(en0, x0, y0, z0, thx0, thy0, txin1, txin2,
     &     tyin1, tyin2, tzin, enp, x, y, z, thx, thy, scatflg)

      implicit none

      include 'const.inc'

      double precision en0, x0, y0, z0, thx0, thy0, rx0, ry0
	double precision txin1, txin2, tyin1, tyin2, tzin
      double precision x, y, z, thx, thy, ep, enp
	double precision thp, php, sinthp, sinphp
      double precision xdis, ydis
	double precision rfx, rfy, rfz, rdist(3), dreact
	double precision thsc, phsc, thsclab, th0, ph0, tharg, pharg
	double precision sig, thxp, thyp, rand
	double precision, external :: rsqrt,ratan
	real vec(3)
	logical scatflg
c     x0, etc. are current location from main
c     txin1 etc. are constraints of this air gap  from main
c      enp,x,y etc are results of passing through this air gap 



c	write(6,*) ' in air gap'
	    scatflg = .false.

      if((x0.ge.txin1.and.x0.le.txin2).and.
     &   (y0.ge.tyin1.and.y0.le.tyin2)) then   ! we're inside air gap

        call ranlux(vec,3)
c	write(6,*) ' rand vector 1, 2',vec(1),vec(2)
c	write(6,*) ' rand vector 3',vec(3)

c       Get polar LAB angles from thx0, thy0
         if(thx0.ge.0..and.thy0.ge.0.) then
            ph0 = ratan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.ge.0.) then
            ph0 = pi + ratan(dtan(thy0)/dtan(thx0))
         elseif(thx0.lt.0..and.thy0.lt.0.) then
            ph0 = pi + ratan(dtan(thy0)/dtan(thx0))
         elseif(thx0.ge.0..and.thy0.lt.0.) then
            ph0 = 2.*pi + ratan(dtan(thy0)/dtan(thx0))
         endif

         th0      = ratan(rsqrt(dtan(thx0)**2 + dtan(thy0)**2))

c       choose scattering angle -- isotropic in CM - z-axis along
c            neutron momentum
	
		thsc  = dacos(1.d0-2.0*vec(1)) ! isotropic (0,pi) bec
		phsc  = 2*pi*vec(2)

c	 New kinetic energy from LAB scattering angle but still in CS 
c      with z along k

        thsclab = ratan(dsin(thsc)/(dcos(thsc)+(m_n/m_air)))
	if (thsclab.lt.0.) thsclab = pi + thsclab

	  ep = en0*(m_n/(m_n+m_air))**2*( dcos(thsclab)+
     &       rsqrt( (m_air/m_n)**2-(dsin(thsclab))**2 ) )**2
	if (ep.le.0.01d0) then 
		write(*,*) ' ep < 0.01:',ep
	    x = 1.d-6; y = 1.d-6
	   goto 100
	end if


c      Transform angles back from CS with z along k, to z along beam axis

	tharg = -dsin(th0)*dsin(thsclab)*dcos(phsc)+dcos(th0)*dcos(thsclab)

	if (tharg.gt.1.0d0) then
		write(*,*) ' theta_p error :  cos(th_p) =',tharg
		tharg=1.0d0
	else if (tharg.lt.-1.0d0) then
		write(*,*) ' theta_p error :  cos(th_p) =',tharg
		tharg=-1.0d0
	endif
	thp = dacos(tharg)

	if (thp.lt.1.0d-9) then
	 pharg = 1.0
	else
	 pharg =  (dsin(thsclab)*dcos(phsc)*dcos(ph0)*dcos(th0) -
     &      dsin(thsclab)*dsin(phsc)*dsin(ph0) + 
     &	  dcos(thsclab)*dcos(ph0)*dsin(th0) )/dsin(thp)		
	end if

	if (pharg.gt.1.0d0) then
c		write(*,*) ' phi_p error :  cos(ph_p) =',pharg
		pharg=1.0d0
	else if (pharg.lt.-1.0d0) then
c		write(*,*) ' phi_p error :  cos(ph_p) =',pharg
		pharg=-1.0d0
	end if

	php = dacos(pharg)
	sinphp =     (dsin(thsclab)*dcos(phsc)*dsin(ph0)*dcos(th0) +
     &      dsin(thsclab)*dsin(phsc)*dcos(ph0) + 
     &      dcos(thsclab)*dsin(ph0)*dsin(th0) )/dsin(thp) 
      if (sinphp.lt.0) php = 2.0d0*pi - php

c
c	write(6,*) ' thsc ', thsc,',  thsclab ',thsclab
c	write(6,*)  ' en0 ', en0,', ep ',ep

c       Get sigma(en0)  en0 in (meV) and find reaction length
c      !!This sigma(E) for air gap is fit of intepreted data from ENDF --
c            for n-14N  , 14N cross section much bigger than 16O  
c            but not much actual data!!!!!

	if (en0.lt.2.2) then
	   sig = 1.0742/(en0/1000.0)**(0.47)    
	 else
	   sig = 6.7799/(en0/1000.0)**(0.162)    
	endif	   

c         sig = 2.6
c
c	write(6,*) 'En ',en0,' sig ',sig

        dreact = -1.0*log(1-vec(3))/(sig*ndenair)  ! n*sig in 1/cm if sig in barn

c	write(6,*) ' nsig ',sig*ndenair, ' dreact ',dreact  

            thxp = ratan(dtan(thp)*dcos(php))
            thyp = ratan(dtan(thp)*dsin(php))
            rfx = dreact*dsin(th0)*dcos(ph0)  ! before scat dist.
            rfy = dreact*dsin(th0)*dsin(ph0)
            rfz = dreact*dcos(th0)
c
c	write(6,*) 'x,y,z',x0,y0,z0
c	write(6,*) 'dx,dy,dz',rfx,rfy,rfz

            rdist(1) = x0 + rfx   ! distance to scatter point
            rdist(2) = y0 + rfy
            rdist(3) = z0 + rfz

c      ****Inner Loop  to determine if scattered, if in target, etc..***	 
       if(  (rdist(1).ge.txin1.and.rdist(1).le.txin2) .and.
     &      (rdist(2).ge.tyin1.and.rdist(2).le.tyin2) .and.
     &      (rdist(3).ge.z0.and.rdist(3).le.(z0+tzin)) ) then  !scat in tar
c       if(  (rdist(3).ge.z0.and.rdist(3).le.(z0+tzin)) ) then  !scat in tar
             
		   xdis  = rdist(1) + (tzin-rfz)*dtan(thxp)		! x loc at end of tar
             ydis  = rdist(2) + (tzin-rfz)*dtan(thyp)		! y loc at end of tar
		   scatflg=.true.
         if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &        (ydis.ge.tyin1.and.ydis.le.tyin2) ) then   ! still in air gap
            x = xdis
            y = ydis
            z = z0 + tzin
            thx = thxp
            thy = thyp
            enp = ep
         else			!  starts in air gap but leaves
c			write(6,*) ' air gap scatter - leaves target '
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6
         endif
	 else     ! doesn't scatter

           xdis = x0 + tzin*dtan(thx0)
           ydis = y0 + tzin*dtan(thy0)

         if( (xdis.ge.txin1.and.xdis.le.txin2) .and.
     &        (ydis.ge.tyin1.and.ydis.le.tyin2) ) then   ! still in air gap
c		write(6,*) ' air gap NO scatter - still in air gap '
            x = xdis
            y = ydis
            z = z0 + tzin
            thx = thx0
            thy = thy0
            enp = en0
         else			!  starts in air gap but leaves
c			write(6,*) ' air gap NO scatter - leaves air gap '
            x = -1.d6
            y = -1.d6
            z = -1.d6
            thx = -1.d6
            thy = -1.d6
            enp = -1.d6
         endif
	 endif
      else				!  never made it into air gap 
c			write(6,*) ' air gap Never in air gap '
         x = -1.d6
         y = -1.d6
         z = -1.d6
         thx = -1.d6
         thy = -1.d6
         enp = -1.d6
      endif
100	continue
      return
      end