!two different sets of constants are used for E < 300 and E > 300eV
	data (a0(i),a1(i),a2(i),a3(i),e_zero(i),e_region(i),i=1,2 )/
     .    2.069989d1,-7.07605d0,4.83547d-1,0.0,5.6914686d1,3.0d2,
     .   3.00207076d-1,4.4915014d-2,-1.554915014d-1,9.50318d-4,0.0,2.d3/

	energy_in=(vx0**2+vy0**2+vz0**2)/clight**2*emass*0.5

      if(energy_in.le.e_region(1))then
      i=1      !E < 300
      else if(energy_in.le.e_region(2))then
      i=2 !300 < E < 3000
      else
      i=3 !E > 1000
      endif

	!f_ref = (reflected electrons)/(total generated secondaries)
      if(i.eq.3)then
          f_ref=0.0
      else
          eln=log(energy_in+e_zero(i))
          f_ref=a0(i)+a1(i)*eln+a2(i)*eln**2+a3(i)*eln**3
          f_ref=exp(f_ref)
      endif
!	to turn off reflective electron, keep the following line
!		f_ref=0.0
        rfel_vs_se=rn(seed2) !rn = random number between 0 and 1

	if(y.eq.0.0d0.and.x.gt.0.0d0)then
		theta_rot=pi/2.0d0
	else if(y.eq.0.0d0.and.x.lt.0.0d0)then
		theta_rot=-pi/2.0d0
	else if(y.gt.0.0d0)then
		theta_rot=atan(-x/y)+pi
	else
		theta_rot=atan(-x/y)+twopi
	endif
!	MAP
	vx=vx0*dcos(theta_rot)+vy0*dsin(theta_rot)
	vy=-vx0*dsin(theta_rot)+vy0*dcos(theta_rot) 
	!y direction is normal to the wall 
	if(vy.eq.0.0d0)then
	theta_in=0.0
	else
	theta_in=atan(dsqrt(vx**2+vz**2)/dabs(vy)) !angle to normal
	endif

	!Main yield calculation, Emax = energy of SEY peak, deltamax = max SEY
	xi=energy_in/Emax*(1.0d0+0.7d0*(1.0d0-cos(theta_in)))
	yield=deltamax*1.11*xi**(-0.35)*(1.-   
     .      exp(-2.3*xi**1.35))*exp(0.5d0*(1.0d0-dcos(theta_in)))
      yield=yield/(1.d0-f_ref)
!     for ture secondary emission
!	secondary  e- energy e_pe
1	e_pe=e20+sige2*tgauss(iseed) !gaussian energy distribution
	if(e_pe.lt.0.) goto 1
	absv=dsqrt(e_pe*2./emass)*clight !incoming speed
	theta2=ctheta2(x,y) 
	theta=asin(rn(seed2)) !cos(theta) distribution
	puci=rn(seed2)*2.0*pi !phi

	vx=absv*sin(theta)*sin(puci)
	vy=absv*cos(theta)
	vz=absv*sin(theta)*cos(puci)
	vxy=dsqrt(vx**2+vy**2)
	theta1=acos(vx/vxy)
	thetavxy=theta1+theta2+pi/2.
	vx=vxy*cos(thetavxy)
	vy=vxy*sin(thetavxy)



!     reflective option
      if(rfel_vs_se.le.f_ref)goto 2


      goto 999

!     reflective e-
	!reflected electrons keep their original energy
2	  trss11=dcos(theta_rot)**2-dsin(theta_rot)**2
      trss22=-trss11
      trss12= 2.0d0*dcos(theta_rot)*dsin(theta_rot)
      trss21=trss12
  	  vx=vx0*trss11+vy0*trss12
	  vy=vx0*trss21+vy0*trss22
      vz=vz0     
      erremit=x*vx+y*vy
!			output emission angle with tangle in pi
	   if(erremit.gt.0.)then
			   write(*,*)'Reflective e- emission err'
			   write(*,*)x,y,vx,vy
			   write(9,*)'Refelective e- emission err'
			   write(9,'(6(2x,1pe12.5))')x,y,z,vx,vy,vz
			   pause
	   endif

999	return  
	end