program main_WUp
c	same as main_slim (main_cwn with no 1000s or 2000s histograms) except targets on opposite side
c		West up (x<0) and East dn (x>0)
c        calculating path length and rotation from entire flight path, not just while in target
      implicit none

      include 'const.inc'


	integer loop
	double precision ep, wp, kp, k1, thp,php,qlast
      integer i, j, counter,lim,Event,cin,cot,trig  !bec
	integer Nenscat1,Ntar1,Nenscat2,Ntar2,Ntar12,NEdet
	integer encnt, Nsdet,Nthdet, Ntar1out, Ntar2out
	integer Nalscat1, Nalscat2, Nthscat1, Nthscat2
	integer Nalout1, Nalout2, Nalin1, Nairin1, Nairout1, Nairscat1
	integer NPSMout, Nincoilout, NOutcoilout, NTargetout   ! bec
	integer numtar1flg, numtar2flg, numendif1, numendif2
c

      real xlam, HRNDM1, vect(5), lm(3), const, xinit, yinit, plr,xtest
      character*28 outfname
      double precision th0, ph0, thindex, enstart, entar1, entar2
      double precision en0, x0, y0, z0, thx0, thy0
      double precision en1, x1, y1, z1, thx1, thy1
      double precision time1
	logical tar2flg,enflg1,enflg2,thflg1,thflg2,alflg1,alflg2
	logical airflg1, histflg, histflg2, empflg1, empflg2
      integer ranno
      integer seedc
      integer n1, n2, n3, Nempty1, Nempty2
C      parameter(ranno=1)
      double precision r2d, lambda, sig, intlen, thscat, sigsom, sigsq
      double precision pl1_ns, th1_ns, pl2_ns, th2_ns, pl_ns, th_ns
      double precision pl1a_sc, pl1b_sc, th1a_sc, th1b_sc
      double precision pl2a_sc, pl2b_sc, th2a_sc, th2b_sc
      double precision pl1_sc, th1_sc, pl2_sc, th2_sc, pl3, th3, phi
	double precision pltot,thtot,pl_del
	double precision lsrot,lsrotsq, lsrotsc,lsrotscsq
	double precision lsrotscle,lsrotsclesq,lsde,lsdesq
	double precision rsrot,rsrotsq, rsrotsc,rsrotscsq
	double precision rsrotscle,rsrotsclesq,rsde,rsdesq
	integer nr,nrs,nrsl,nl,nls,nlsl,nder,nders,ndel,ndels
      CHARACTER day*2,month*2,hour*2,minute*2, filename*20
c-----------------------------------
      integer qmin, qmax
      double precision k0, qdummy, q
c-----------------------------------
c
c     -------------------------
c     fill the ntuple 'spin.rz' //// ntuple variables ....
      integer max_paw, iquest
      integer istat, icycle, datetime(8)
      parameter(max_paw = 1000000) ! PAW common block size
      real quest, p
	character date*12, time*12, zone*12
c      common/quest/iquest(100)
      common/pawc/p(max_paw)

c      time1=time()
c      ranno=mod(time1,2000000)
c      ranno=mod(time1,1089926518)

       CALL Date_and_time(date,time,zone,datetime)  !bec 

      day= CHAR(INT(datetime(3)/10)+48)//CHAR(mod(datetime(3),10)+48)
      month= CHAR(int(datetime(2)/10)+48)//CHAR(MOD(datetime(2),10)+48)
      hour= CHAR(INT(datetime(5)/10)+48)//CHAR(mod(datetime(5),10)+48)
      minute= CHAR(INT(datetime(6)/10)+48)//CHAR(mod(datetime(6),10)+48)

	filename ='xportWUp-'//month//day//'-'//hour//minute//'.o'
	open(unit=3,FILE=filename,Status='new')

	write(*,*) '  Nevents: ', nevts        ! max_paw  !bec
	write(*,*) '  Start time: ', date,time			!bec
	write(3,*) '  Nevents: ', nevts        ! max_paw  !bec
	write(3,*) '  Start time: ', date,time			!bec
	write(*,*) '  Upstream target on left, x<0.  Dn on right '
	write(3,*) '  Upstream target on left, x<0.  Dn on right '

c      print *,"asking for seed"
c      read(*,*) ranno;
c      ranno=seedc()
	ranno=datetime(6)*datetime(7)*datetime(8)+datetime(7)+datetime(8)	!bec
      write(*,*) 'ranNo is ', ranno         
	write(3,*) 'ranNo is ', ranno   
c      stop
c      stop/N/B/msarsour/spin/new_ih/
      call hlimit(max_paw)
c      iquest(10) = 6500000
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c
      call hrget(100,'lam.hbook','')
c      call hrget(102,'xyspot.hbook','')
      call hrget(151,'gradx.hbook','')
      call hrget(152,'gradx.hbook','')
      call hrget(301,'sq.hbook','')
      call hcopyr(301,201,"",1,100,0.,0.,"")

c-----------------------------------------------------------------------
c      read(*,*) ranno
      call rluxgo(3,ranno,0,0)     ! RANLUX luxury level 2
cc      call rluxgo(3,710,0,0)     ! RANLUX luxury level 2
c-----------------------------------------------------------------------
c      Print *, "setting up Hbooks" 
      call hbook1(305,'West q sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(306,'West w sim(1/psec)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(307,'West Eo sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(308,'West Ep sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(309,'West ko sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(310,'West kp sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(205,'East q sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(206,'East w sim(1/psec)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(207,'East Eo sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(208,'East Ep sim(meV)',3000,0.,40.0,0.)      ! bec test q choices
      call hbook1(209,'East ko sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
      call hbook1(210,'East kp sim(1/Ang)',100,0,4.0,0.)      ! bec test q choices
	call hbook1(1001,'x (cm) ',100,-3.5,3.5,0.)
		call hbook1(1,'x (cm) start',100,-3.5,3.5,0.)
		call hbook1(2,'x (cm) guide out',100,-3.5,3.5,0.)
		call hbook1(3,'x (cm) incoil in',100,-3.5,3.5,0.)
		call hbook1(4,'x (cm) incoil out',100,-3.5,3.5,0.)
		call hbook1(5,'x (cm) outcoil out',100,-3.5,3.5,0.)
		call hbook1(6,'x (cm) ASM out',100,-3.5,3.5,0.)
      call hbook1(20,'E (meV) start Both',2000,0.,40.,0.)

      call hbook1(1002,'y (cm) UpWtar',100,-3.,3.,0.)
c     call hbook1(1003,'thx (mrad) Up',160,-1.6,1.6,0.)
c      call hbook1(1004,'thy (mrad) Up',160,-1.6,1.6,0.)
      call hbook1(1005,'E (meV)/ non-scat WE',2000,0.,40.,0.)
c      call hbook1(1006,'Path len (cm)/ non-scat UpB',100,35.0,100.0,0.)
c      call hbook1(1007,'th_rot (mrad)/ non-scat UpB',100,0.,10000.,0.)
c      call hbook1(1008,'Path len (cm)/ non-scat UpR',100,35.0,100.,0.)
c      call hbook1(1009,'th_rot (mrad)/ non-scat UpR',100,0.,10000.,0.)
c      call hbook1(1010,'Path len (cm)/ non-scat UL',100,35.0,100.0,0.)
c      call hbook1(1011,'th_rot (mrad)/ non-scat UL',100,0.,10000.,0.)
      call hbook1(1012,'E (meV)/ scat UpWE',2000,0.,40.,0.)
c      call hbook1(1013,'Path len (cm)/ scat UpB',100,35.0,100.0,0.)
c      call hbook1(1014,'th_rot (mrad)/ scat UpB',100,0.,10000.,0.)
c      call hbook1(1015,'Path len (cm)/ scat UpR',100,35.0,100.0,0.)
c      call hbook1(1016,'th_rot (mrad)/ scat UpR',100,0.,10000.,0.)
c      call hbook1(1017,'Path len (cm)/ scat UpL',100,35.0,100.0,0.)
c      call hbook1(1018,'th_rot (mrad)/ scat UpL',100,0.,10000.,0.)
*
      call hbook1(2001,'x (cm) Dn',100,-3.5,3.5,0.)
      call hbook1(2002,'y (cm) Dn',100,-3.,3.,0.)
c      call hbook1(2003,'thx (mrad) Dn',160,-1.6,1.6,0.)
c      call hbook1(2004,'thy (mrad) Dn',160,-1.6,1.6,0.)
      call hbook1(2005,'E (meV)/ non-scat DnWE',2000,0.,40.,0.)
c      call hbook1(2006,'Path len (cm)/ non-scat DnB',100,41.5,41.7,0.)
c      call hbook1(2007,'th_rot (mrad)/ non-scat DnB',100,0.,1000.,0.)
c      call hbook1(2008,'Path len (cm)/ non-scat DnR',100,41.5,41.7,0.)
c      call hbook1(2009,'th_rot (mrad)/ non-scat DnR',100,0.,1000.,0.)
c      call hbook1(2010,'Path len (cm)/ non-scat DnL',100,41.5,41.7,0.)
c      call hbook1(2011,'th_rot (mrad)/ non-scat DnL',100,0.,1000.,0.)
      call hbook1(2012,'E (meV)/ scat DnWE',2000,0.,40.,0.)
c      call hbook1(2013,'Path len (cm)/ scat DnB',100,41.,45.,0.)
c      call hbook1(2014,'th_rot (mrad)/ scat DnB',100,0.,1000.,0.)
c      call hbook1(2015,'Path len (cm)/ scat DnR',100,41.,45.,0.)
c      call hbook1(2016,'th_rot (mrad)/ scat DnR',100,0.,1000.,0.)
c      call hbook1(2017,'Path len (cm)/ scat DnL',100,41.,45.,0.)
c      call hbook1(2018,'th_rot (mrad)/ scat DnL',100,0.,1000.,0.)
*
      call hbook1(3001,'x (cm) Det',100,-3.5,3.5,0.)
      call hbook1(3002,'y (cm) Det',100,-3.,3.,0.)
      call hbook1(3003,'thx (mrad) Det',200,0.,0.1,0.)
      call hbook1(3004,'thy (mrad) Det',200,0.,0.1,0.)
      call hbook1(3005,'E (meV) DetWE',2000,0.,40.,0.)
      call hbook1(4005,'Lambda (Angs) DetWE',1000,0.,50.,0.)
      call hbook1(3006,'Path length (cm) DeWE',200,114.58,114.7,0.)
      call hbook1(3007,'th_rotation (mrad) DetWE',200,-5.,5.0,0.)
      call hbook1(3008,'Path length (cm) DetEast',200,114.58,114.7,0.)
      call hbook1(3009,'th_rotation (mrad) DetEast',200,-5.,5.0,0.)
      call hbook1(3010,'Path length (cm) DetWest',200,114.58,114.7,0.)
      call hbook1(3011,'th_rotation (mrad) DetWest',200,-5.,5.0,0.)
      call hbook1(3012,'E (meV) DetEast',2000,0.,40.,0.)
      call hbook1(3013,'E (meV) DetWest',2000,0.,40.,0.)
      call hbook1(3014,'E (meV) scat DetWE',2000,0.,40.,0.)
      call hbook1(3020,'E_loss (meV) scat DetWE',2000,0.,10.,0.)
c
c   adding more spectra below for diagnostics   bec
	call hbook2(8010,'sig tot (b/atom) vs. q (1/Ang)',100,0.,4.,
     &      200,0.,1.,0.)
      call hbook1(4010,'West sig tot (b/atom)',200,0.,1.,0.)
c	call hbook1(4011,'West interaction len (cm)',200,0.,2000.,0.)
      call hbook1(4020,'East sig tot (b/atom)',200,0.,1.,0.)
c	call hbook1(4021,'East interaction len (cm)',200,0.,2000.,0.)
	call hbook1(4012,'theta scattered (deg)',200,0.,180.,0.)
	call hbook1(4013,'West sig Sommers (b/atom)',200,0.,1.,0.)
	call hbook1(4014,'West sig Sq (b/atom)',200,0.,1.,0.)
	call hbook1(4023,'East sig Sommers (b/atom)',200,0.,1.,0.)
	call hbook1(4024,'East sig Sq (b/atom)',200,0.,1.,0.)
	call hbook1(4015,'theta-p (deg) U',200,0.,180.,0.)
	call hbook1(4016,'phi-p (deg) U',200,0.,360.,0.)
	call hbook1(4017,'th-0 (deg)',200,0.,10.,0.)
	call hbook1(4018,'phi-0 (deg) ',200,0.,360.,0.)
	call hbook1(5012,'theta scattered (deg) D',200,0.,180.,0.)
	call hbook1(5015,'theta-p (deg) D',200,0.,180.,0.)
	call hbook1(5016,'phi-p (deg) D',200,0.,360.,0.)
	call hbook1(5017,'theta-x (deg) U',200,-2.,2.,0.)
	call hbook1(5018,'theta-y (deg) U',200,-2.,2.,0.)
	call hbook1(5019,'theta-x (deg) D',200,-2.,2.,0.)
	call hbook1(5020,'theta-y (deg) D',200,-2.,2.,0.)
	call hbook1(5021,'phi (deg) af-guide',200,0.,360.,0.)
	call hbook1(5022,'phi (deg) af-target 1',200,0.,360.,0.)	
	call hbook1(5023,'phi (deg) af-target 2',200,0.,360.,0.)

      call hbook1(6001,'Path length (cm) enter Up WE',200,8.0,
     &      12.0,0.)
      call hbook1(6002,'th_rotation (mrad) enter Up WE',200,-1.,
     &      3.,0.)
      call hbook1(6003,'Path length (cm) Det scat WE',200,114.58,
     &      114.7,0.)
      call hbook1(6004,'th_rotation (mrad) Det scat WE',200,-5.,
     &      5.0,0.)
      call hbook1(6005,'Path length (cm) Det scat West',200,114.58,
     &      114.7,0.)
      call hbook1(6006,'th_rotation (mrad) Det scat West',200,-5.,
     &      5.0,0.)
      call hbook1(6007,'Path length (cm) Det scat East',200,114.58,
     &      114.7,0.)
      call hbook1(6008,'th_rotation (mrad) Det scat East',200,-5.,
     &      5.0,0.)
      call hbook1(6009,'Path length (cm) Det no scat WE',200,114.58,
     &      114.7,0.)
      call hbook1(6010,'th_rotation (mrad) Det no scat WE',200,-5.,
     &      5.0,0.)
      call hbook1(6011,'Path length (cm) Det no scat West',200,114.58,
     &      114.7,0.)
      call hbook1(6012,'th_rotation (mrad) Det no scat West',200,-5.,
     &      5.0,0.)
      call hbook1(6013,'Path length (cm) Det no scat East',200,114.58,
     &      114.7,0.)
      call hbook1(6014,'th_rotation (mrad) Det no scat East',
     &      200,-20.,20.,0.)
	call hbook2(6020,'th_rot (mrad) vs. E (meV) det WE',1000,0.,40.,
     &	  200,-5.,5.,0.)
	call hbook2(6021,'th rot (mrad) vs. Path len (cm) detWE',200,
     &      114.58,114.7,200,-20.0,5.0,0.)
c	call hbook2(6022,'th rot (mrad) vs. E (meV) detR',1000,0.,40.,
c     &			200,-10.,50.,0.)
c	call hbook2(6023,'th rot (mrad) vs. Path len (cm) detR',200,480.6,477.34,
c     &			200,-10.,50.,0.)
c	call hbook2(6024,'th rot (mrad) vs. E (meV) detL',1000,0.,40.,
c     &			200,-10.,50.,0.)
c	call hbook2(6025,'th rot (mrad) vs. Path len (cm) detL',200,480.6,477.34,
c     &			200,-10.,50.,0.)
c	call hbook2(6026,'E (meV) vs. Path len (cm) detR',200,480.6,477.34,
c     &			1000,0.,40.,0.)
c	call hbook2(6027,'E (meV) vs. Path len (cm) detL',200,480.6,477.34,
c     &			1000,0.,40.,0.)
	call hbook2(8011,'Thx vs. E (meV) scat DetB',1000,0.,40.,
     &       200,0.,0.1,0.)
	call hbook2(8012,'Thx vs. E-loss (meV) scat DetB',1000,0.,40.,200,
     &	0.,0.1,0.)

c
c------------Input file -------------------------
c      open(39,file='inp.dat',status='old')
c      read(39,*) gl, gin, T
c------------------------------------------------
      r2d=180./pi
c------------
      counter = 0
      Event = 0
      cin=0
      cot=0
	Nenscat1=0		!bec
	Ntar1=0		!bec
	Nthscat1=0
	Nenscat2=0		!bec
	Ntar2=0		!bec
	Nthscat2=0
	Nsdet=0		!bec
	NEdet=0		!bec
	NTargetout=0; Ntar1out=0; Ntar2out=0    !bec
	Nthdet=0
	NPSMout=0 ; Nincoilout=0 ; NOutcoilout = 0   !bec
	lsrot=0;lsrotsq=0;lsrotsc=0;lsrotscsq=0;lsrotscle=0;lsrotsclesq=0
	lsde=0;lsdesq=0;rsrot=0;rsrotsq=0;rsrotsc=0;rsrotscsq=0
	rsrotscle=0
	rsrotsclesq=0;rsde=0;rsdesq=0;nr=0;nrs=0;nrsl=0;nl=0;nls=0;nlsl=0
	nder=0;nders=0;ndel=0;ndels=0
	Nalscat1=0;Nalscat2=0;Nalout1=0;Nalout2=0; Nalin1=0
	tar2flg = .false.  ! bec
	Nairin1 = 0; Nairout1 = 0; Nairscat1 = 0
	histflg = .false. ; histflg2 = .false.
	numtar1flg = 0; numendif1 = 0
	numtar2flg = 0; numendif2 = 0
c
c

      do 101 i = 1, nevts
c
         Event = 1
		pltot = 0.0
		thtot = 0.0
		enflg1=.false.
		thflg1=.false.
		enflg2=.false.
		thflg2=.false.

c	print *, " Envent number ", i
c
c     Generate initial neutron wavelength
c
         xlam = HRNDM1(100)
c
c     The energy calculated from the neutron's wavelength
c
       if(xlam.lt.1.) goto 100
c		call ranlux(vect,4)  !  temporary
       en0 = (1/(2*m_n))*(2*pi*h_par/xlam)**2
c		en0 = 50.*vect(1)	!  temporary
	 enstart = en0	
            call hf1(20,real(en0),1.)
	   if (en0.le.1.0) encnt =encnt+1
c
c     Generates the initial trajectory of the neutron.
c
 321        call ranlux(vect,4)

         thindex = ng6index*xlam  ! NG6 about 1.7mrad/ang

	if (thindex.lt.1E-7) thindex = 1E-7
	if (abs(1.d0-(1.d0-dcos(thindex))*vect(1)).gt.1.0d0) then
		write(6,*) ' Arccos(theta_start)>1 ' 
		goto 321
	end if
		th0  = dacos(1.d0-(1.d0-dcos(thindex))*vect(1)) ! isotropic (0,thindex) bec

	if (th0.lt.1.E-8) th0=1.E-8

         ph0  = 2.0d0*pi*vect(2) ! azimuthal angle (radians)
         call hf1(4017,real(th0*180./3.14159),1.)
         call hf1(4018,real(ph0*180./3.14159),1.)

c         call HRNDM2(102,xinit,yinit)
c         x0=xinit/10.           ! cm
c         y0=yinit/10.           ! cm
c         write(18,*) x0, y0; goto 100
         xtest = xin*(2.*vect(3)-1.)
c
		if (i.eq.1) write(*,*) ' ! Running with No x gradient !'
		if (i.eq.1) write(3,*) ' ! Running with No x gradient !'
		x0 = xtest
c
c		if (i.eq.1) write(*,*) ' ! Running with R to L x gradient !'
c		if (i.eq.1) write(3,*) ' ! Running with R to L x gradient !'
c		x0 = 3.0d0 - 2.d0*HRNDM1(151)  !bec R to L gradient
c
c		if (i.eq.1) write(*,*) ' ! Running with L to R x gradient !'
c		if (i.eq.1) write(3,*) ' ! Running with L to R x gradient !'
c		x0 = 2.d0*HRNDM1(151) - 3.0d0  !bec L to R gradient
c
c   this gives double gradient -- gradient across each half...
c         if(xtest.gt.0.) then
c            x0 = HRNDM1(151)  ! gives gradient only acrros half
cc            x0 = 3.d0 - HRNDM1(151)
c         elseif(xtest.lt.0) then
c            x0 = HRNDM1(151) - 3.
cc            x0 = -HRNDM1(151)
c         endif
c
c         y0 = 2.25*(2.*vect(4)-1.)
c      square cross section
	  y0 = yin*(2.*vect(4)-1.)    
         z0= 0.0d0
c
c     Calculate the incident angles in the xz-plane and yz-plane 
c       note: these are not "direction cosines"= angle between momentum and x (or y)
c
         thx0 = datan(dtan(th0)*dcos(ph0))
         thy0 = datan(dtan(th0)*dsin(ph0))
c

	   NPSMout = NPSMout + 1   !  bec count exiting PSM
         call hf1(1,real(x0),1.)

********************************************************************
*     Neutrons passing thruogh the gap (3.0 cm) after the PSM before input guides
*   this is being added in to old code, so x0, y0, z0, etc just redefined as result
*    of this 3 cm transport
*
c         x1 = x0 + tzair1*dtan(thx0)
c         y1 = y0 + tairz1*dtan(thy0)
c         z1 = z0 + tairz1

	 if((x0.ge.tairx1.and.x0.le.tairx2).and.
     &   (y0.ge.tairy1.and.y0.le.tairy2)) Nairin1 = Nairin1 + 1
c
	   call  airgap(en0, x0, y0, z0, thx0, thy0, 
     &         tairx1,tairx2,tairy1,tairy2,tairz1, en1, x1, y1, z1, 
     &         thx1,thy1,airflg1)
	if (airflg1) Nairscat1 = Nairscat1+1  ! count scatters with energy change  bec
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100

cbec		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot   !  mrad
	if (pltot.gt.500) write(18,*) '0th air gap',pltot
	Nairout1 = Nairout1 + 1
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if
********************************************************************
*     Neutrons passing through FIRST Al Window
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	 if((x0.ge.talx1.and.x0.le.talx2).and.
     &   (y0.ge.taly1.and.y0.le.taly2)) Nalin1 = Nalin1 + 1
c
	   call  alwindow(en0, x0, y0, z0, thx0, thy0, 
     &         talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, 
     &         thx1,thy1,alflg1)
	if (alflg1) Nalscat1 = Nalscat1+1  ! count scatters with energy change  bec
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
cbec		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0 
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	Nalout1 = Nalout1 + 1
	if (thtot.gt.1e9) then
		write(6,*) ' Al thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if


********************************************************************
c     Neutrons passing thruogh the guide"
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
         call guide(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,
     &      thx1,thy1)
         if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
         call hf1(2,real(x1),1.)
cbec		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' guide thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

*
********************************************************************
********************************************************************
*     Neutrons passing thruogh the input coil
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
         call incoil(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1)
         if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		Nincoilout = Nincoilout + 1  !  bec   counttransmission through guide and Input coil
	         call hf1(3,real(x1),1.)
cbec		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'incoil ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if
*
********************************************************************
********************************************************************
*     Neutrons passing thruogh the gap (10.0 cm) between the input
*     coil and the target.
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  airgap(en0, x0, y0, z0, thx0, thy0, 
     &         tairx1,tairx2,tairy1,tairy2,tairz2, en1, x1, y1, z1, 
     &         thx1,thy1,airflg1)

       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100

	         call hf1(4,real(x1),1.)
	call hf1(5017,real(180./3.14159*thx1),1.)
	call hf1(5018,real(180./3.14159*thy1),1.)
	   if(thx1.ge.0..and.thy1.ge.0.) then
            phi = datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.ge.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.lt.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.ge.0..and.thy1.lt.0.) then
            phi = 2.*pi + datan(dtan(thy1)/dtan(thx1))
         endif
	call hf1(5021,real(180./3.14159*phi),1.)
		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
            call hf1(6001,real(pltot),1.)
            call hf1(6002,real(thtot),1.)
	if (pltot.gt.500) write(18,*) 'gap 1 ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

*
********************************************************************
*     Neutrons passing through Al window
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  alwindow(en0, x0, y0, z0, thx0, thy0, 
     &         talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, 
     &         thx1,thy1,alflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' Al thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if


********************************************************************
********************************************************************
c	print *, " first half target"
*     Neutrons passing thruogh the 1st half of the target (41.6 cm).
C  Note: below checks only where the neutron is in x direction.  Once it enters target
c   target.f checks the y-direction too, so might have neutrons that get into target.f
c   that don't run target.f really, since the y-values are out of range, hence sig or thscat values
c   get nothing put in them 
*
************Defining x<0 to imply West which in main_WUp.exe has upstream target
C
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1

c      Given current energy determine q from S(q) and delE from omega(q)
         k0 = sqrt(2*m_n*en0)/h_par

         qmin = 0
c         qmax = int(k0*25.)     ! for sq_1.hbook
         qmax = int(2.0*k0*25.)     ! bec  qmax = 2*k  from p-consv. 
         if((2.0*k0).gt.3.6) qmax=int(3.6*25.)
c        Check to see if we've been here before and thus created HIST 302
c         histflg set false before loop, then always true
         if(histflg)  call hdelet(302)
		histflg = .true.
         call hcopyr(301,302,"",qmin,qmax,0.,0.,"")
c
c         qdummy=HRNDM1(301)  ! bec 
c   HRNDM1 replaces spec with integrated spec!, so next pass through the loop
c   hcopyr is copying the integrated spec into 302, not the original spec!!
c   without the qdummy=HRNDM1(301) line, the simulated q-spec (305) looks
c   much more reasonable!  bec
		loop=0
 123        q=HRNDM1(302)  
	   if (loop.gt.1000) then 
		write(6,*) ' 1 q-loop ',loop,en0,q
		write(6,*) k0,kp,ep
	   end if
		loop=loop+1
       	call wq(q,wp,en0,k0)  !  wq modified to prevent E,p consv. probs
		ep  = en0- (h_par*wp)
		if (ep.le.0.) goto 123

		kp=sqrt(2*m_n*ep)/h_par
		if (dabs((k0**2+kp**2-q**2)/(2.*k0*kp)).le.1.d0) goto 124
      goto 123
 124	continue


         call hf1(305,real(q),1.)
         call hf1(306,real(wp),1.)
         call hf1(307,real(en0),1.)
         call hf1(308,real(ep),1.)
         call hf1(309,real(k0),1.)
         call hf1(310,real(kp),1.)

c
	enflg1 = .false.
	thflg1 = .false.
	empflg1 = .false.

         if( (x0.ge.tarx3.and.x0.le.tarx4) ) then  ! West=pos. x, target UPstream
	sig = 0.0 ; sigsom = 0.0 ; sigsq =0.0 ; thscat = 0.0
	thp = 0.0; php = 0.0
			call target(q, en0, x0, y0, z0, thx0, thy0, tarx3, tarx4,
     &           en1, x1, y1, z1, thx1, thy1, sig,
     &           thscat,sigsom,sigsq,enflg1,thflg1,thp,php,Ntar1,
     &           pl1a_sc,pl1b_sc)  ! bec

c		make sure really went in target -- we know x-position is ok, check y-position
		if (y0.ge.tary1.and.y0.le.tary2) then
			Ntar1=Ntar1+1
			if (enflg1) Nenscat1 = Nenscat1+1  ! count scatters with energy change  bec
			if (thflg1) Nthscat1 = Nthscat1+1  ! count scatters with theta change  bec
c      capture diagnostics for target interaction stuff   bec
			call hf1(4010,real(sig),1.) 
		call hf2(8010,real(q),real(sig),1.)  !temporary
c			call hf1(4011,real(intlen),1.) 
			call hf1(4012,real(180./3.14159*thscat),1.)
			call hf1(4013,real(sigsom),1.)  
			call hf1(4014,real(sigsq),1.)  
			call hf1(4015,real(180./3.14159*thp),1.)  
			call hf1(4016,real(180./3.14159*php),1.)  
		end if

          if ((x1.ge.tarx3.and.x1.le.tarx4).and.
     &        (y1.ge.tary1.and.y1.le.tary2)) 
     &          Ntar1out = Ntar1out+1     ! out of target west side only

         elseif( (x0.ge.tarx1.and.x0.le.tarx2) ) then  ! East=neg. x, empty UP
			call empty(en0, x0, y0, z0, thx0, thy0, tarx1, tarx2,
     &           en1, x1, y1, z1, thx1, thy1,empflg1)

	if (empflg1) Nempty1 = Nempty1+1  ! count scatters with energy change  bec
	
         else
			x1 = -1.d6
			y1 = -1.d6
			z1 = -1.d6
			thx1 = -1.d6
			thy1 = -1.d6
			en1 = -1.d6
         endif
*
        if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
*
		entar1 = en1


         call hf1(1001,real(x1),1.)
         call hf1(1002,real(y1),1.)
c         call hf1(1003,real(thx1),1.)
c         call hf1(1004,real(thy1),1.)
         if(en0.eq.en1) then
            pl1_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
            th1_ns = gy*mag_B*pl1_ns*0.1/sqrt(2.*en1/m_n)
c            write(16,*) th1_ns
            pl1_sc = 0.
            th1_sc = 0.
 		  pltot = pl1_ns + pltot
		  thtot = th1_ns + thtot
            call hf1(1005,real(en1),1.)
c            call hf1(1006,real(pl1_ns),1.)
c            call hf1(1007,real(th1_ns),1.)
            if(x0.gt.0..and.x1.gt.0.) then
c               call hf1(1008,real(pl1_ns),1.)
c               call hf1(1009,real(th1_ns),1.)
            elseif(x0.lt.0..and.x1.lt.0.) then
c               call hf1(1010,real(pl1_ns),1.)
c               call hf1(1011,real(th1_ns),1.)
            endif
         else
            th1a_sc = gy*mag_B*pl1a_sc*0.1/sqrt(2.*en0/m_n)
            th1b_sc = gy*mag_B*pl1b_sc*0.1/sqrt(2.*en1/m_n)
            pl1_ns = 0.
            th1_ns = 0.
            pl1_sc  = pl1a_sc + pl1b_sc
            th1_sc  = th1a_sc + th1b_sc
 		  pltot = pl1_sc + pltot
		  thtot = th1_sc + thtot
            call hf1(1012,real(en1),1.)
c            call hf1(1013,real(pl1_sc),1.)
c            call hf1(1014,real(th1_sc),1.)
            if(x0.gt.0..and.x1.gt.0.) then
c               call hf1(1015,real(pl1_sc),1.)
c               call hf1(1016,real(th1_sc),1.)
            elseif(x0.lt.0..and.x1.lt.0.) then
c               call hf1(1017,real(pl1_sc),1.)
c               call hf1(1018,real(th1_sc),1.)
            endif
         endif
	if (pltot.gt.500) write(18,*) 'target 1 ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl1_ns ',pl1_ns
		write(6,*) ' pl1_sc ',pl1_sc,',  en ',en1
		write(6,*) ' en0 ',en0, pltot
	end if

*
********************************************************************
*     Neutrons passing through Al window
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  alwindow(en0, x0, y0, z0, thx0, thy0, 
     &         talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, 
     &         thx1,thy1,alflg1)
      if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' Al thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if


********************************************************************
********************************************************************
c	print *,"   Gap btw targets" 
*     Neutrons passing thruogh the gap (pi-coil = 7.3 cm) between the
*     two tagets = 8.0cm.
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  airgap(en0, x0, y0, z0, thx0, thy0, 
     &         tairx1,tairx2,tairy1,tairy2,tairz3, en1, x1, y1, z1, 
     &         thx1,thy1,airflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100

c       do NOT allow neutrons from one target cell into opposite side cell
         if( ((x0.ge.tarx1).and.(x0.le.tarx2)).and.((x1.lt.tarx1).or.
     &        (x1.gt.tarx2)) ) 
     &		goto 100
         if( ((x0.ge.tarx3).and.(x0.le.tarx4)).and.((x1.lt.tarx3).or.
     &        (x1.gt.tarx4)) ) 
     &		goto 100

	call hf1(5019,real(180./3.14159*thx1),1.)
	call hf1(5020,real(180./3.14159*thy1),1.)
	   if(thx1.ge.0..and.thy1.ge.0.) then
            phi = datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.ge.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.lt.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.ge.0..and.thy1.lt.0.) then
            phi = 2.*pi + datan(dtan(thy1)/dtan(thx1))
         endif
	call hf1(5022,real(180./3.14159*phi),1.)


c*************PI Coil reversal of rotation angle
c	using 1/2 gap length -- rotate to halfway point, reverse, rotate the other half
	    pl_del = 0.5*sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot

	if (thtot.gt.1e4) then
		write(6,*) ' thtot ', thtot, ',  pl1_ns ',pl1_ns
		write(6,*) ' pl1_sc ',pl1_sc,',  en ',en1
		write(6,*) ' en0 ',en0, pltot
	    write(6,*) thx0,thy0, thx1, thy1
	end if

c     reverse thtot in center of gap and rotate the other half of gap length
		thtot = - thtot
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot

	if (pltot.gt.500) write(18,*) 'gap 2 ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if
	 
*
********************************************************************
*     Neutrons passing through Al window
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  alwindow(en0, x0, y0, z0, thx0, thy0, 
     &         talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, 
     &         thx1,thy1,alflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' Al thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if


********************************************************************
********************************************************************
c	print *, "        second half target"
*     Neutrons passing thruogh the 2nd part of the target (41.6 cm).
*
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
c         call hf1(203,real(x0),1.)
	enflg2 = .false.
	thflg2 = .false.
	empflg2 = .false.

         if( (x0.gt.tarx3.and.x0.lt.tarx4) ) then  ! West = pos. x, empty down
			call empty(en1, x0, y0, z0, thx0, thy0, tarx3, tarx4,
     &           en1, x1, y1, z1, thx1, thy1, empflg2)

	if (empflg2) Nempty2 = Nempty2+1  ! count scatters with energy change  bec

         elseif( (x0.gt.tarx1.and.x0.lt.tarx2) ) then  ! East = neg. x, target DOWN
			if (enflg1.or.thflg1) then   ! if scat in t1, now scat in both targets
				Ntar12=Ntar12+1
			endif
	    sig = 0.0 ; sigsom = 0.0 ; sigsq =0.0 ; thscat = 0.0 
	    thp = 0.0; php = 0.0
c
c         Given current energy determine q from S(q) and delE from omega(q)
c          Doing this only if in east side where target is located
c
           k0 = sqrt(2*m_n*en0)/h_par

           qmin = 0
c            qmax = int(k0*25.)     ! for sq_1.hbook
           qmax = int(2.0*k0*25.)     ! bec  qmax = 2*k  from p-consv. 
           if((2.0*k0).gt.3.6) qmax=int(3.6*25.)
           if(histflg2)  call hdelet(202)
	 	 histflg2 = .true.

           call hcopyr(201,202,"",qmin,qmax,0.,0.,"")
c
		 loop=0
 223        q=HRNDM1(202)  
	       if (loop.gt.1000) then 
		     write(6,*) ' 2 q-loop ',loop,en0,q
		     write(6,*) k0,kp,ep
	       end if
		 loop=loop+1
       	 call wq(q,wp,en0,k0)  !  wq modified to prevent E,p consv. probs
		 ep  = en0- (h_par*wp)
		 if (ep.le.0.) goto 223

		 kp=sqrt(2*m_n*ep)/h_par

		 if (dabs((k0**2+kp**2-q**2)/(2.*k0*kp)).le.1.d0) goto 224
           goto 223
 224	     continue

         call hf1(205,real(q),1.)
         call hf1(206,real(wp),1.)
         call hf1(207,real(en0),1.)
         call hf1(208,real(ep),1.)
         call hf1(209,real(k0),1.)
         call hf1(210,real(kp),1.)

			call target(q, en0, x0, y0, z0, thx0, thy0, tarx1, tarx2,
     &            en1, x1, y1, z1, thx1, thy1, sig,
     &            thscat,sigsom,sigsq,enflg2,thflg2,thp,php,Ntar2,
     &            pl2a_sc,pl2b_sc)  ! bec

c		make sure really went in target -- we know x-position is ok, check y-position
		if (y0.ge.tary1.and.y0.le.tary2) then
			Ntar2=Ntar2+1
			if (enflg2) Nenscat2 = Nenscat2+1  ! count scatters with ef.ne.ei  bec
			if (thflg2) Nthscat2 = Nthscat2+1  ! count scatters with thsc>0  bec
c      capture diagnostics for target interaction stuff   bec
			call hf1(4020,real(sig),1.) 
c			call hf1(4021,real(intlen),1.) 
			call hf1(5012,real(180./3.14159*thscat),1.)
			call hf1(4023,real(sigsom),1.)  
			call hf1(4024,real(sigsq),1.)  
			call hf1(5015,real(180./3.14159*thp),1.)  
			call hf1(5016,real(180./3.14159*php),1.)  
		end if

           if ((x1.ge.tarx1.and.x1.le.tarx2).and.
     &         (y1.ge.tary1.and.y1.le.tary2)) Ntar2out = Ntar2out+1     ! out of target east side only

         else
			x1 = -1.d6
			y1 = -1.d6
			z1 = -1.d6
         endif
*
         if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
*
		entar2 = en1

         call hf1(2001,real(x1),1.)
         call hf1(2002,real(y1),1.)
c         call hf1(2003,real(thx1),1.)
c         call hf1(2004,real(thy1),1.)
         if(.not.enflg2) then
            pl2_ns = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
            th2_ns = gy*mag_B*pl2_ns*0.1/sqrt(2.*en1/m_n)
            pl2_sc = 0.
            th2_sc = 0.
	      pltot = pl2_ns + pltot
		  thtot = th2_ns + thtot
            call hf1(2005,real(en1),1.)
c            call hf1(2006,real(pl2_ns),1.)
c            call hf1(2007,real(th2_ns),1.)
            if(x0.gt.0..and.x1.gt.0.) then
c               call hf1(2008,real(pl2_ns),1.)
c               call hf1(2009,real(th2_ns),1.)
            elseif(x0.lt.0..and.x1.lt.0.) then
c               call hf1(2010,real(pl2_ns),1.)
c               call hf1(2011,real(th2_ns),1.)
            endif
         else
            th2a_sc = gy*mag_B*pl2a_sc*0.1/sqrt(2.*en0/m_n)
            th2b_sc = gy*mag_B*pl2b_sc*0.1/sqrt(2.*en1/m_n)
            pl2_ns = 0.
            th2_ns = 0.
            pl2_sc  = pl2a_sc + pl2b_sc
            th2_sc  = th2a_sc + th2b_sc
	      pltot = pl2_sc + pltot
		  thtot = th2_sc + thtot
            call hf1(2012,real(en1),1.)
c            call hf1(2013,real(pl2_sc),1.)
c            call hf1(2014,real(th2_sc),1.)
            if(x0.gt.0..and.x1.gt.0.) then
c               call hf1(2015,real(pl2_sc),1.)
c               call hf1(2016,real(th2_sc),1.)
            elseif(x0.lt.0..and.x1.lt.0.) then
c               call hf1(2017,real(pl2_sc),1.)
c               call hf1(2018,real(th2_sc),1.)
            endif
         endif
	if (pltot.gt.500) write(18,*) 'target 2 ',pltot

	if (thtot.gt.1e4) then
		write(6,*) ' thtot ', thtot, ',  pl2_ns ',pl2_ns
		write(6,*) ' pl2_sc ',pl2_sc,',  en ',en1
		write(6,*) ' en0 ',en0, pltot
	    write(6,*) thx0,thy0, thx1, thy1
	end if

*
********************************************************************
		NTargetout = NTargetout + 1
********************************************************************
********************************************************************
*     Neutrons passing through Al window
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  alwindow(en0, x0, y0, z0, thx0, thy0, 
     &         talx1,talx2,taly1,taly2,talz1, en1, x1, y1, z1, 
     &         thx1,thy1,alflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'guide ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' Al thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

***********************************************************************
*     Neutrons passing thruogh the gap (13.0 cm) between the taget and
*     the output coil.
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  airgap(en0, x0, y0, z0, thx0, thy0, 
     &         tairx1,tairx2,tairy1,tairy2,tairz4, en1, x1, y1, z1, 
     &         thx1,thy1,airflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100

	   if(thx1.ge.0..and.thy1.ge.0.) then
            phi = datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.ge.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.lt.0..and.thy1.lt.0.) then
            phi = pi + datan(dtan(thy1)/dtan(thx1))
         elseif(thx1.ge.0..and.thy1.lt.0.) then
            phi = 2.*pi + datan(dtan(thy1)/dtan(thx1))
         endif
	call hf1(5023,real(180./3.14159*phi),1.)

	    pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'gap 3 ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

********************************************************************
********************************************************************
*     Neutrons passing thruogh the output coil.
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call outcoil(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1)
         if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
		         call hf1(5,real(x1),1.)
cbec	    pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'out coil ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

*
********************************************************************
********************************************************************
*     Neutrons passing through the gap (10.0 cm) between the output coil
*	and the analyzing super mirror (ASM).
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
	   call  airgap(en0, x0, y0, z0, thx0, thy0, 
     &         tairx1,tairx2,tairy1,tairy2,tairz5, en1, x1, y1, z1, 
     &         thx1,thy1,airflg1)
       if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100

cbec	    pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'gap 5 ',pltot
	NOutcoilout = NOutcoilout + 1
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

********************************************************************
********************************************************************
*     Neutrons passing thruogh the supper mirror."
c	print *,"   Supper Mirror" 
*
c      Update x,y,z,thx,thy,en values to new values
	x0=x1;y0=y1;z0=z1;thx0=thx1;thy0=thy1;en0=en1
c
         call sm(en0,x0,y0,z0,thx0,thy0,x1,y1,z1,thx1,thy1)
         if(x1.lt.-1d5.or.y1.lt.-1d5) goto 100
         cot=cot+1
		         call hf1(6,real(x1),1.)
cbec	    pl_del = sqrt((x1-x0)**2+(y1-y0)**2+(z1-z0)**2)
		pl_del = 0.d0
		pltot =  pl_del + pltot
          thtot = gy*mag_B*pl_del*0.1/sqrt(2.*en1/m_n) + thtot
	if (pltot.gt.500) write(18,*) 'ASM ',pltot
	if (thtot.gt.1e9) then
		write(6,*) ' thtot ', thtot, ',  pl_del ',pl_del
		write(6,*) ' en ',en1, ',  pltot ',pltot
	end if

*
********************************************************************
********************************************************************
*
         call hf1(3001,real(x1),1.)
         call hf1(3002,real(y1),1.)
         call hf1(3003,real(thx1),1.)
         call hf1(3004,real(thy1),1.)
         pl3 = pl1_ns + pl2_ns + pl1_sc + pl2_sc
         th3 = th1_ns + th2_ns + th1_sc + th2_sc

         call hf1(3005,real(en1),1.)
         lambda=(2.*pi*h_par)/sqrt(2*m_n*en1)
         call hf1(4005,real(lambda),1.)
         call hf1(3006,real(pltot),1.)
         call hf1(3007,real(thtot),1.)
         if(x1.ge.0.) then   ! x>0 = West (=old right)
            call hf1(3008,real(pltot),1.)
            call hf1(3009,real(thtot),1.)
            call hf1(3012,real(en1),1.)
			rsrot = rsrot + dsin(thtot/1000.)
			rsrotsq = rsrotsq + dsin(thtot/1000.)**2
			nr = nr + 1
         else				   ! x<0 = East (=old left)
            call hf1(3010,real(pltot),1.)
            call hf1(3011,real(thtot),1.)
            call hf1(3013,real(en1),1.)
			lsrot = lsrot + dsin(thtot/1000.)
			lsrotsq = lsrotsq + dsin(thtot/1000.)**2
			nl = nl + 1
         endif

	if (entar1.ne.enstart)  numendif1 = numendif1 + 1
	if (enflg1)  numtar1flg = numtar1flg + 1
	if (entar2.ne.entar1)  numendif2 = numendif2 + 1
	if (enflg2)  numtar2flg = numtar2flg + 1

c	if (entar1 .ne. enstart) write(17,*) numendif1, enstart, entar1, 
c     &    enstart-entar1
c	if ((entar1.ne.enstart).and.(.not.enflg1)) 
c     &    write(6,*) ' flag problem'


		if(enflg1.or.enflg2) then   ! scattered and changed energy in either target
			 NEdet=NEdet+1
			call hf1(6003, real(pltot),1.)
			call hf1(6004, real(thtot),1.)
				if (x1.lt.0.) then			! x<0 = East (=old left)
				  call hf1(6005, real(pltot),1.)
				  call hf1(6006, real(thtot),1.)
c				  call hf2(6024,real(en1),real(thtot),1.)
c				  call hf2(6025,real(pltot),real(thtot),1.)
c				  call hf2(6027,real(pltot),real(en1),1.)
					lsrotsc = lsrotsc + dsin(thtot/1000.)
					lsrotscsq = lsrotscsq + dsin(thtot/1000.)**2
					nls = nls + 1
					if ((entar1.le.0.35).or.(entar2.le.0.35)) then
						lsrotscle = lsrotscle + dsin(thtot/1000.)
						lsrotsclesq=lsrotsclesq+dsin(thtot/1000.)**2
						nlsl = nlsl + 1
					end if
					if(enflg1)  then
						call hf1(3020,real(enstart-entar1),1.)
						call hf2(8011,real(entar2),real(thx1),1.)
						call hf2(8012,real(enstart-entar1),real(thx1),1.)
	if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', 
     &         thx1,entar2,enstart-entar1
						lsde = lsde + enstart-entar1
						lsdesq = lsdesq + (enstart-entar1)**2
						ndel = ndel + 1
c			write(6,*) ' 1W eloss tar1 tar2 ', entar1-entar2 
					endif
					if(enflg2)  then
						call hf1(3020,real(entar1-entar2),1.)
						call hf2(8011,real(entar2),real(thx1),1.)
						call hf2(8012,real(entar1-entar2),real(thx1),1.)
	if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', thx1,entar2,
     &         entar1-entar2
						lsde = lsde + entar1-entar2
						lsdesq = lsdesq + (entar1-entar2)**2
						ndel = ndel + 1
c			write(6,*) ' 2W eloss tar1 tar2 ', entar1-entar2 
					endif
				else						! x>0 = West (=old right)
				  call hf1(6007, real(pltot),1.)
				  call hf1(6008, real(thtot),1.)
c				  call hf2(6022,real(en1),real(thtot),1.)
c				  call hf2(6023,real(pltot),real(thtot),1.)
c				  call hf2(6026,real(pltot),real(en1),1.)
					rsrotsc = rsrotsc + dsin(thtot/1000.)
					rsrotscsq = rsrotscsq + dsin(thtot/1000.)**2
					nrs = nrs + 1
					if ((entar1.le.0.35).or.(entar2.le.0.35)) then
						rsrotscle = rsrotscle + dsin(thtot/1000.)
						rsrotsclesq=rsrotsclesq+dsin(thtot/1000.)**2
						nrsl = nrsl + 1
c			write(6,*) ' 1E eloss tar1 tar2 ', entar1-entar2 
					end if
					if(enflg1)  then
						call hf1(3020,real(enstart-entar1),1.)
						call hf2(8011,real(entar2),real(thx1),1.)
						call hf2(8012,real(enstart-entar1),real(thx1),1.)
	if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', 
     &         thx1,entar2,enstart-entar1
						rsde = rsde + enstart-entar1
						rsdesq = rsdesq + (enstart-entar1)**2
						nder = nder + 1
c			write(6,*) ' 2E eloss tar1 tar2 ', entar1-entar2 
					endif
					if(enflg2)  then
						call hf1(3020,real(entar1-entar2),1.)
						call hf2(8011,real(entar2),real(thx1),1.)
						call hf2(8012,real(entar1-entar2),real(thx1),1.)
	if (thx1.gt.0.1) write(6,*) ' thx, e, eloss', 
     &         thx1,entar2,entar1-entar2
						rsde = rsde + entar1-entar2
						rsdesq = rsdesq + (entar1-entar2)**2
						nder = nder + 1
					endif
				endif
		else
			call hf1(6009, real(pltot),1.)
			call hf1(6010, real(thtot),1.)
				if (x1.lt.0.) then
				  call hf1(6011, real(pltot),1.)
				  call hf1(6012, real(thtot),1.)
c				  call hf2(6024,real(entar2),real(thtot),1.)
c				  call hf2(6025,real(pltot),real(thtot),1.)
c				  call hf2(6027,real(pltot),real(en1),1.)
				else
				  call hf1(6013, real(pltot),1.)
				  call hf1(6014, real(thtot),1.)
c				  call hf2(6022,real(entar2),real(thtot),1.)
c				  call hf2(6023,real(pltot),real(thtot),1.)
c				  call hf2(6026,real(pltot),real(en1),1.)
				endif
		endif
c	if (en1.lt.0.36) then
c		write(16,888) ' Edet<0.36 -- E0, E1, E2 ', enstart, entar1, entar2
c		write(16,889) ' E0-1, E1-2 ', enstart-entar1, entar1-entar2
c		write(16,888) ' thx, thy', thx1, thy1
c		write(16,889) ' q, thscat ', q, thscat
c		write(16,889) ' Th-rot, pltot ',thtot, pltot
c		write(16,*)
c	endif
c 888  format(2x,a30,1x,1pe11.4,1x,1pe11.4,1x,1pe11.4)
c 889  format(2x,a30,1x,1pe11.4,1x,1pe11.4)
		if(thflg1.or.thflg2) Nthdet=Nthdet+1
		if(enflg1.or.enflg2.or.thflg1.or.thflg2) then  ! bec
            call hf1(3014,real(en1),1.)
		  Nsdet=Nsdet+1
         endif
	call hf2(6020,real(en1),real(thtot),1.)
c	call hf2(6021,real(pltot),real(thtot),1.)
c	if (i.lt.100) write(17,*) ' end ', pltot,thtot
********************************************************************
*



 34      format(i12,2x,a6)
 100        do j=1,10
            lim = float(j)*float(nevts)/10.
            if(float(i).eq.lim) write(*,34) i, 'events'
         enddo

 101	continue


	print *
 	print *, "   through loop, writing to ", filename


	write(6,*) ' num of enflg1 ',numtar1flg
	write(6,*) ' num of endif 1',numendif1
	write(3,*) ' num of enflg1 ',numtar1flg
	write(3,*) ' num of endif 1',numendif1
	write(6,*) ' num of enflg2 ',numtar2flg
	write(6,*) ' num of endif 2',numendif2
	write(3,*) ' num of enflg1 ',numtar2flg
	write(3,*) ' num of endif 2',numendif2


	print *
		if (hegastar) then
			write(6,*) ' 4He gas target '
			write(3,*) ' 4He gas target '
		elseif (airtar) then
			write(6,*) ' Air target '
			write(3,*) ' Air target '
		else
			write(6,*) ' Liquid He target; Temp= ',T,' K '
			write(3,*) ' Liquid He target; Temp= ',T,' K '
		endif


	print *
      write(*,*) 'input                   ',nevts
      write(*,*) 'PSM output              ',NPSMout
	write(*,*) 'Air gap 1 out           ',Nairout1
	write(*,*) 'Al window 1 in          ',Nalin1
	write(*,*) 'Al window 1 out         ',Nalout1
      write(*,*) 'Input Coil output       ',Nincoilout
	write(*,*) 'West target IN          ',Ntar1
	write(*,*) 'West target OUT         ',Ntar1out
	write(*,*) 'East target IN          ',Ntar2 
	write(*,*) 'East target OUT         ',Ntar2out 
	write(*,*) 'E-W Target Outputs      ',NTargetout 
      write(*,*) 'Output Coil output      ',NOutcoilout
	write(*,*) 'ASM output (geom. only) ', cot     
	write(*,*)

      write(3,*) 'input                   ',nevts
      write(3,*) 'PSM output              ',NPSMout
	write(3,*) 'Air gap 1 out           ',Nairout1
	write(3,*) 'Al window 1 out         ',Nalout1
      write(3,*) 'Input Coil output       ',Nincoilout
	write(3,*) 'West target IN          ',Ntar1
	write(3,*) 'West target OUT         ',Ntar1out
	write(3,*) 'East target IN          ',Ntar2 
	write(3,*) 'East target OUT         ',Ntar2out 
	write(3,*) 'E-W Target Outputs      ',NTargetout 
      write(3,*) 'Output Coil output      ',NOutcoilout
	write(3,*) 'ASM output (geom. only) ', cot     
	write(3,*)

	write(*,*) 'Number scattered in Al window 1 ',Nalscat1
	write(3,*) 'Number scattered in Al window 1 ',Nalscat1
	write(*,*) 'Number scattered in Air gap   1 ',Nairscat1
	write(3,*) 'Number scattered in Air gap   1 ',Nairscat1
	write(*,*) 'Number scattered in UpTar (En0.ne.En1)',Nenscat1
	write(3,*) 'Number scattered in UpTar (En0.ne.En1)',Nenscat1
	write(*,*) 'Number scattered in DnTar (En0.ne.En1)',Nenscat2
	write(3,*) 'Number scattered in DnTar (En0.ne.En1)',Nenscat2
	write(*,*) 'Number scattered in EmptyTar East ',Nempty1
	write(*,*) 'Number scattered in EmptyTar West ',Nempty2
	write(3,*) 'Number scattered in EmptyTar East ',Nempty1
	write(3,*) 'Number scattered in EmptyTar West ',Nempty2
c	write(*,*) 'Number Theta scattered in Target 1 ',Nthscat1
c	write(3,*) 'Number Theta scattered in Target 1 ',Nthscat1
c	write(*,*) 'Number Theta scattered in Target 2 ',Nthscat2
c	write(3,*) 'Number Theta scattered in Target 2 ',Nthscat2
c
c     -------------------------

      if(ranno<1e1) then
         Write(outfname,'(1Hh,I1,6H.hbook)')ranno
      elseif(ranno<1e2) then
         Write(outfname,'(1Hh,I2,6H.hbook)')ranno
      elseif(ranno<1e3) then
         Write(outfname,'(1Hh,I3,6H.hbook)')ranno
      elseif(ranno<1e4) then
         Write(outfname,'(1Hh,I4,6H.hbook)')ranno
      elseif(ranno<1e5) then
         Write(outfname,'(1Hh,I5,6H.hbook)')ranno
      elseif(ranno<1e6) then
         Write(outfname,'(1Hh,I6,6H.hbook)')ranno
      elseif(ranno<1e7) then
         Write(outfname,'(1Hh,I7,6H.hbook)')ranno
      elseif(ranno<1e8) then
         Write(outfname,'(1Hh,I8,6H.hbook)')ranno
      endif
      call hrput(0,outfname,'N')
c     -------------------------
c     Adding some printouts at end of run  -- bec
	print *
	Call Date_and_time(date,time,zone,datetime)	

	write(6,*) '  End time: ', date, time			
	write(6,*) '      HBOOK filename: ',outfname
	write(6,*) 

	write(6,*) '  Frac # of scatters with Ef.ne.Ei tar1 ', 
     &		float(Nenscat1)/float(Ntar1)	
	write(6,*) '  Frac # of scatters with Ef.ne.Ei tar2 ', 
     &		float(Nenscat2)/float(Ntar2)	
	write(6,*) '  Frac num of tar1 entering tar2 ',
     &				float(Ntar12)/float(Ntar1)
	write(6,*) '  Frac of scat of total reaching det ',
     &				float(Nsdet)/float(Nevts)
	write(6,*) '   Frac of det hits that were scattered ',
     &				float(Nsdet)/float(cot)
	write(6,*) '   Frac of det hits that were scattered with thsc>0 ',
     &				float(Nthdet)/float(cot)
	write(6,*) '  Frac of det hits that were scattered with Ef.ne.Ei ',
     &				float(NEdet)/float(cot)

	write(6,713) float(nr), float(nl)
	write(6,714) float(nr)/float(nl)
	write(6,717) float(nr)-float(nl),Sqrt(float(nl+nr))


	write(6,*)
	write(6,700) mag_B

	write(6,704) 1000*dasin(rsrot/float(nr)),
     &		1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))
	write(6,703) 1000*dasin(lsrot/float(nl)),
     &		1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))
	write(6,711) 1000*(dasin(rsrot/float(nr))-dasin(lsrot/float(nl))),
     &		1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 
     &                + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2)
	write(6,*)

	if (nls.gt.0) then

c	write(6,*) rsde, nder, lsde, ndel
	write(6,702) rsde/float(nder),
     &             sqrt((rsdesq-rsde**2/nder)/nder/(nder-1))
	write(6,701) lsde/float(ndel),
     &             sqrt((lsdesq-lsde**2/ndel)/ndel/(ndel-1))
	write(6,706) 1000*dasin(rsrotsc/float(nrs)),
     &	   1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))
	write(6,705) 1000*dasin(lsrotsc/float(nls)),
     &	   1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))
	write(6,712) 1000*(dasin(rsrotsc/float(nrs))
     & - dasin(lsrotsc/float(nls))),
     & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))**2 
     &        + dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))**2)
	
	write(6,*)
	if (nrsl.gt.1) then
		write(6,708) 1000*dasin(rsrotscle/float(nrsl)),
     &  1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1)))
	else
		write(6,708) 0.,0.
	endif
	if(nlsl.gt.1) then
		write(6,707) 1000*dasin(lsrotscle/float(nlsl)),
     &  1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1)))
	else
		write(6,707) 0.,0.
	end if
	
	write(6,710) rsrotsc/rsrot,rsrotscle/rsrot
	write(6,709) lsrotsc/lsrot,lsrotscle/lsrot
	write(6,*)
	write(6,715) float(nrs), float(nls)
	write(6,716) float(nrs)/float(nls)
	write(6,718) float(nrs)-float(nls),Sqrt(float(nls+nrs))

	else
	write(6,*) '  NO scattering so limited output written '
	endif


 700  format(3x,'   Rotations for Incoil Ouput to Outcoil input  114.6~cm with 
     & B(mGauss)=',1f9.4)
 701	format(3x,' Avg. Energy loss at Det East Dn ',1pe11.4,' +- '
     &           ,1pe11.4)
 702	format(3x,' Avg. Energy loss at Det West Up ',1pe11.4,' +- '
     &           ,1pe11.4)
 703	format(3x,' Avg. Rot East Dn              (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 704	format(3x,' Avg. Rot West Up               (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 705	format(3x,' Avg. Rot Scat East Dn          (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 706	format(3x,' Avg. Rot Scat West Up        (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 707	format(3x,' Avg. Rot Scat East Dn E<0.35meV (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 708	format(3x,' Avg. Rot Scat West Up E<0.35meV  (mrad) ',1pe11.4,
     &              ' +- ',1pe11.4)
 709	format(3x,' Frac of rotation from scat East Dn   ',1pe11.4,
     &			', oflow En ',1pe11.4)
 710	format(3x,' Frac of rotation from scat West Up    ',1pe11.4,
     &			', oflow En ',1pe11.4)
 711	format(3x,' Rotation UpW-DownE           (mrad) ',1pe11.4,
     &            ' +- ',1pe11.4)
 712	format(3x,' Rotation UpqW-DownE scat      (mrad) ',1pe11.4,
     &            ' +- ',1pe11.4)
 713  format(3x,' # in DetW Up',1pe13.6,2x,', # in Det E Dn ',1pe13.6)
 714  format(3x,' Ratio # DetWUp/DetEDn ,'1pe13.6)
 717	format(3x,' Det Counts UpW-DownE   ',1pe11.4,' +- ',1pe11.4)
 715  format(3x,' # scat in DetW Up ',1pe11.4,2x,
     &       ', # scat in Det E Dn ',1pe11.4)
 716  format(3x,' Ratio # scat DetWUp/DetEDn ,'1pe11.4)
 718	format(3x,' Scat Det Counts UpW-DownE   ',1pe11.4,' +- ',1pe11.4)


c
	write(3,*) '  End time: ', date, time			
	write(3,*) '      HBOOK filename: ',outfname
	write(3,*) 

	write(3,*) '  Frac # of scatters with Ef.ne.Ei tar1 ', 
     &		float(Nenscat1)/float(Ntar1)	
	write(3,*) '  Frac # of scatters with Ef.ne.Ei tar2 ', 
     &		float(Nenscat2)/float(Ntar2)	
	write(3,*) '  Frac num of tar1 entering tar2 ',
     &				float(Ntar12)/float(Ntar1)
	write(3,*) '  Frac of scat of total reaching det ',
     &				float(Nsdet)/float(Nevts)
	write(3,*) '   Frac of det hits that were scattered ',
     &				float(Nsdet)/float(cot)
	write(3,*) '   Frac of det hits that were scattered with thsc>0 ',
     &				float(Nthdet)/float(cot)
	write(3,*) '  Frac of det hits that were scattered with Ef.ne.Ei ',
     &				float(NEdet)/float(cot)

	write(3,713) float(nr), float(nl)
	write(3,714) float(nr)/float(nl)
	write(3,717) float(nr)-float(nl),Sqrt(float(nl+nr))


	write(3,*)
	write(3,700) mag_B

	write(3,704) 1000*dasin(rsrot/float(nr)),
     &		1000*dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))
	write(3,703) 1000*dasin(lsrot/float(nl)),
     &		1000*dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))
	write(3,711) 1000*(dasin(rsrot/float(nr))-dasin(lsrot/float(nl))),
     &		1000*sqrt(dasin(sqrt((lsrotsq-lsrot**2/nl)/nl/(nl-1)))**2 
     &                + dasin(sqrt((rsrotsq-rsrot**2/nr)/nr/(nr-1)))**2)
	write(3,*)

	if (nls.gt.0) then

c	write(3,*) rsde, nder, lsde, ndel
	write(3,702) rsde/float(nder),
     &             sqrt((rsdesq-rsde**2/nder)/nder/(nder-1))
	write(3,701) lsde/float(ndel),
     &             sqrt((lsdesq-lsde**2/ndel)/ndel/(ndel-1))
	write(3,706) 1000*dasin(rsrotsc/float(nrs)),
     &	   1000*dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))
	write(3,705) 1000*dasin(lsrotsc/float(nls)),
     &	   1000*dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))
	write(3,712) 1000*(dasin(rsrotsc/float(nrs))
     & - dasin(lsrotsc/float(nls))),
     & 1000*sqrt(dasin(sqrt((lsrotscsq-lsrotsc**2/nls)/nls/(nls-1)))**2 
     &        + dasin(sqrt((rsrotscsq-rsrotsc**2/nrs)/nrs/(nrs-1)))**2)

	write(3,*)
	if (nrsl.gt.1) then
		write(3,708) 1000*dasin(rsrotscle/float(nrsl)),
     &  1000*dasin(sqrt((rsrotsclesq-rsrotscle**2/nrsl)/nrsl/(nrsl-1)))
	else
		write(3,708) 0.,0.
	endif
	if(nlsl.gt.1) then
		write(3,707) 1000*dasin(lsrotscle/float(nlsl)),
     &  1000*dasin(sqrt((lsrotsclesq-lsrotscle**2/nlsl)/nlsl/(nlsl-1)))
	else
		write(3,707) 0.,0.
	end if
	
	write(3,*)
	write(3,710) rsrotsc/rsrot,rsrotscle/rsrot
	write(3,709) lsrotsc/lsrot,lsrotscle/lsrot
	write(3,715) float(nrs), float(nls)
	write(3,716) float(nrs)/float(nls)
	write(3,718) float(nrs)-float(nls),Sqrt(float(nls+nrs))

	else
	write(3,*) '  NO scattering so limited output written '
	endif
c

	close(unit=3)
c bec      return
      end