!  ********************************************************
      subroutine CANONICAL(pxc,pyc)
!
! Correct Px and Py for vector potential of solenoid field
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      pxc = pp(1) - 0.5d0 * pcharge * xp(2) * bfld(3)
      pyc = pp(2) + 0.5d0 * pcharge * xp(1) * bfld(3)
!
      return
      end
!  ********************************************************
      subroutine CLEAR
c
c  Clear statistics arrays
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      do i=1,20
      do j=1,5        
      st(j,i) = 0.
      end do
      st(6,i) = 1.e6
      st(7,i) = -1.e6
      end do
c
      do i=1,20
      do j=1,14        
      stsp(j,i) = 0.
      end do
      stsp(6,i) = 1.e6
      stsp(7,i) = -1.e6
      stsp(12,i) = 1.e6
      stsp(13,i) = -1.e6
      end do
!
      return
      end

c  *********************************************************
      subroutine COVARIANCE
!
!  Compute covariance matrix of particle distribution
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      dimension q(6),r(6,6),s(6,6,6),cor2(6,6),cor3(6,6,6)
      dimension cov(6,6),sd(6),v(6)
      real*8 n
      integer indx(6)
      character*2 lab(6)
!
      data lab/'x','y','ct','Px','Py','Pz'/
!
      if( ncovar .eq. 0 ) go to 900
!
      do 600 iz=1,ncovar
! Test is this is a desired covariance plane
      if( icvdes(iz) .ne. j7rec ) go to 600
!
      do 50 i=1,6
      q(i) = 0.
      do 50 j=1,6
      r(i,j) = 0.
      do 50 k=1,6
      s(i,j,k) = 0.
50    continue
      n = 0.
!
      do 100 ip=1,npart
      call READ_EVT(ip)
!
      if( ipflg .ne. 0 ) go to 100
!
      v(1) = xp(1)
      v(2) = xp(2)
      v(3) = tp * cc
      v(4) = pp(1)
      v(5) = pp(2)
      v(6) = pp(3)
!      if( refpar .ne. 0 ) then
!        v(3) = (tp -tpref) * cc
!        v(6) = pp(3) - ppref(3)
!      end if
!
!  Accumulate sums
      n = n + evtwt
!
      do i=1,6
      q(i) = q(i) + v(i) * evtwt
      do j=1,6
      r(i,j) = r(i,j) + v(i)*v(j) * evtwt     
      do k=1,6
      s(i,j,k) = s(i,j,k) + v(i)*v(j)*v(k) * evtwt
      end do
      end do  
      end do  
!
100   continue
!
      if( n .eq. 0. ) go to 900
!
!  Get expectation values
      do i=1,6
      q(i) = q(i) / n          ! 1st moment (mean)
      do j=1,6
      r(i,j) = r(i,j) / n      ! 2nd moment
      do k=1,6
      s(i,j,k) = s(i,j,k) / n  ! 3rd moment
      end do
      end do
      end do
!
!  Get standard deviations
      do i=1,6
      arg = r(i,i) - q(i)**2
      if( arg .gt. 0. ) then
        sd(i) = SQRT( arg )
      else
        sd(i) = -998.
      end if
      end do
!
!  Convert moments to covariances and correlations
      do i=1,6
      do j=1,6
      cov(i,j) = r(i,j) - q(i)*q(j)
      fac = sd(i) * sd(j)
      if( fac .gt. 1e-10 ) then
        cor2(i,j) = cov(i,j) / fac
      else
        cor2(i,j) = -999.
      end if
      do k=1,6
      fac = sd(i) * sd(j) * sd(k)
      if( fac .gt. 1e-10 ) then
        cor3(i,j,k) = s(i,j,k) - q(i)*r(j,k) - q(j)*r(i,k)
     &    - q(k)*r(i,j) + 2.*q(i)*q(j)*q(k)
        cor3(i,j,k) = cor3(i,j,k) / fac 
      else
        cor3(i,j,k) = -999.
      end if
      end do
      end do
      end do
!
! Compute dispersions
      if( sd(6) .gt. 0. ) then
        disx = cor2(1,6) * sd(1) * q(6) / sd(6)
        disy = cor2(2,6) * sd(2) * q(6) / sd(6)
        disxp = cor2(4,6) * sd(4) / sd(6)
        disyp = cor2(5,6) * sd(5) / sd(6)
      else
        disx = -999.
        disy = -999.
        disxp = -999.
        disyp = -999.
      end if
!
      write(2,*)
      write(2,110) j7rec
110   format(15('*'),t18,' PARTICLE COVARIANCE MATRIX, PLANE ',i5,
     &   t63,15('*'))
      write(2,115)
115   format(t12,'Mean values')
      do i=1,6
      if( i .eq. 3 ) then
        write(2,118) lab(i),q(i),'c*tpref',cc*tpref 
      else if( i .eq. 6 ) then
        write(2,118) lab(i),q(i),'pzref',ppref(3) 
      else
        write(2,118) lab(i),q(i)
118     format(t5,a2,t10,e13.6,t30,a7,t38,e13.6)
      end if
      end do
!
      write(2,*)                    
      write(2,120)
120   format(t10,'SD',t30,'Second order correlations')
      write(2,130) (lab(i),i=1,5)
130   format(t25,a2,t35,a2,t45,a2,t55,a2,t65,a2)
!
      do i=1,6
      write(2,150) lab(i),sd(i),(cor2(i,j),j=1,i-1)
150   format(a3,e12.4,2x,5f10.3)
      end do
!
      write(2,*)
      write(2,'(a)')' Third order correlations'
      write(2,160) (lab(i),i=1,6)
160   format(t14,a2,t24,a2,t34,a2,t44,a2,t54,a2,t64,a2)
!
      do i=1,6
      do j=i,6
      if( i.eq. 1) then
        write(2,'(2a3,6f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      else if( i .eq. 2 ) then
        write(2,'(2a3,10x,5f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      else if( i .eq. 3 ) then
        write(2,'(2a3,20x,4f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      else if( i .eq. 4 ) then
        write(2,'(2a3,30x,3f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      else if( i .eq. 5 ) then
        write(2,'(2a3,40x,2f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      else if( i .eq. 6 ) then
        write(2,'(2a3,50x,f10.3)') lab(i),lab(j),(cor3(i,j,k),k=i,j)
      end if
      end do
      write(2,*)
      end do
!
      write(2,'(a)')' Dispersions  [m]'
      write(2,165) 'D(X)= ',disx,'D(Y)= ',disy,"D'(X)= ",
     &  disxp,"D'(Y)= ",disyp  
165   format(a, e12.5, t20, a, e12.5, t40, a, e12.5, t60, a, e12.5)
!
! Get LBNL value of 6-D emittance
      call LUDCMP(cov,6,6,indx,d)
!
      if( d .eq. 0. ) then        ! flag singular matrix
        d = -998.
        go to 580
      end if
!
      do i=1,6
      d = d * cov(i,i)
      end do
      if( d .gt. 0. ) then
        d = SQRT( d )
      else
        d = -998.
      end if
!
580   continue
! Normalization = 1 / (mc)^3
      d = d / pmass**3
!      emit6(iz) = d
      write(2,*)
      write(2,'(a,e13.6)') '6-D emittance [m^3] = ',d
      write(2,*)
!
      call LBNL_EMIT(iz)
!
600   continue        ! end of covariance plane loop
!
900   continue
      return
      end
c  *********************************************************
      subroutine DISP_CORR(n,emx,emy,xarg,yarg)
!
! Compute 2D emittances in presence of dispersion
! Uses correction scheme of D. Neuffer
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      dimension emx(5),emy(5),emd(6),evx(5),evy(5)
      real*8 n
!
      do i=1,6
      emd(i) = 0.
      end do
!
      do 100 ip=1,npart
!
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 100
!
      pzc = pp(3)
      x = xp(1)
      y = xp(2)
      if( pxycorr ) then 
        call CANONICAL(pxc,pyc)  
      else
        pxc = pp(1)
        pyc = pp(2)
      end if
!
!      cpx = pxc / pzc         ! x'
      cpx = pxc
!      cpy = pyc / pzc         ! y'
      cpy = pyc
!
      emd(1) = emd(1) + evtwt * pzc
      emd(2) = emd(2) + evtwt * pzc * pzc
      emd(3) = emd(3) + evtwt * x * pzc
      emd(4) = emd(4) + evtwt * cpx * pzc
      emd(5) = emd(5) + evtwt * y * pzc
      emd(6) = emd(6) + evtwt * cpy * pzc
!
100   continue
!
! Get expectation values
      do i=1,5
      evx(i) = emx(i) / n
      evy(i) = emy(i) / n
      end do
      do i=1,6
      emd(i) = emd(i) / n
      end do
!
      vp = emd(2) - emd(1)*emd(1)
      if( vp .gt. 0. ) then       ! get dispersions
        fac = emd(1) / vp
        dx = (emd(3)-evx(1)*emd(1)) * fac
        dxp = (emd(4)-evx(3)*emd(1)) * fac
        dy = (emd(5)-evy(1)*emd(1)) * fac
        dyp = (emd(6)-evy(3)*emd(1)) * fac
      else
        xarg = -1.        ! flag failure
        yarg = -1. 
        go to 900       
      end if
!
! Correct the x and y distributions
      fac = vp / emd(1) / emd(1)
      evx(2) = evx(2) - dx*dx*fac
      evx(4) = evx(4) - dxp*dxp*fac
      evx(5) = evx(5) - dx*dxp*fac
      evy(2) = evy(2) - dy*dy*fac
      evy(4) = evy(4) - dyp*dyp*fac
      evy(5) = evy(5) - dy*dyp*fac
!
      xarg = evx(2)*evx(4) - evx(5)*evx(5)
      yarg = evy(2)*evy(4) - evy(5)*evy(5)
!
900   continue
      end
c  *********************************************************
      subroutine EMITTANCE
!
!  Compute emittance of beam at this plane
!   including possible tail cuts and Bz corrections
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      dimension emx(5),emy(5),emz(5),stem(2,6)
      dimension xmean(6),xstand(6),ystand(6)
      real*8 n
      integer ncut(6)
!
      data iteste/20/,sig_thre/1e-4/
!
      if( nemit .eq. 0 ) go to 900
! ---------------------------------------------------------
      do 600 iz=1,nemit
! Test if this is one of the desired emittance planes
      if( iemdes(iz) .ne. j7rec ) go to 600
!
      iflg = 0
! .........................................................
! Tail cutting logic
      if( sigmacut ) then
        knt = 0       ! count iterations
        do i=1,6
        ystand(i) = 1e6
        ncut(i) = 10
        end do
!
100     continue      ! enter here for iterations on tail cut
!
        knt = knt + 1
        if( knt .gt. iteste ) then
           iflg = 1
          go to 300
        end if
!
        do j=1,6      ! zero array for {x...Pz} statistics
        do i=1,2
        stem(i,j) = 0.
        end do
        end do
!
        stev = 0.
!..........................................................
!  Read back particle information for this plane
        do 200 ip=1,npart
!
        call READ_EVT(ip)
        if( ipflg .ne. 0 ) go to 200
!
        if( knt .eq. 1 ) emevt(iz) = emevt(iz) + evtwt

!
! Correct Px and Py for vector potential
        if( pxycorr ) then 
          call CANONICAL(pxc,pyc)  
        else
          pxc = pp(1)
          pyc = pp(2)
         end if
!
        energy = SQRT( pxc**2 + pyc**2 + pp(3)**2 + pmass**2 )
        z = pp(3) / energy * cc * tp
!
        stev = stev + evtwt
        stem(1,1) = stem(1,1) +evtwt* xp(1)
        stem(2,1) = stem(2,1) +evtwt* xp(1)**2
        stem(1,2) = stem(1,2) +evtwt* xp(2)
        stem(2,2) = stem(2,2) +evtwt* xp(2)**2
        stem(1,3) = stem(1,3) +evtwt* z
        stem(2,3) = stem(2,3) +evtwt* z*z
        stem(1,4) = stem(1,4) +evtwt* pxc
        stem(2,4) = stem(2,4) +evtwt* pxc**2
        stem(1,5) = stem(1,5) +evtwt* pyc
        stem(2,5) = stem(2,5) +evtwt* pyc**2
        stem(1,6) = stem(1,6) +evtwt* pp(3)
        stem(2,6) = stem(2,6) +evtwt* pp(3)**2
!
200     continue        ! end of particle loop
! .........................................................
        if( stev .eq. 0. ) then
          iflg = 2
          go to 300
        end if  
!
        call MEANSD(stev,stem,xmean,xstand)
!
! Test if tail cutting process has converged
! Stop if fractional change in sigma < sig_thre and ncuts=0
        do 250 i=1,6
        denom = xstand(i)
        if( denom .le. 0. ) go to 250
!
        arg = ABS(ystand(i) - xstand(i)) / denom
        if( arg.gt.sig_thre .or. ncut(i).ne.0 ) then
! If any of the 6 variables requires a tail cut, load and cut all 6
          do j=1,6
          ystand(j) = xstand(j)
          end do
!
          call SIGMA_CUT(xmean,xstand,ncut) ! cut tails on {x...Pz}
!
          go to 100
        end if
!
250     continue
      end if      ! end of test on sigmacut
! ---------------------------------------------------------
300   continue
!
!  Accumulate final statistics for emittance calculation
!
!  We do x and y first because we may need them for corrections
!   to the z emittance
!
      do i=1,5
      emx(i) = 0.d0
      emy(i) = 0.d0
      end do
!
      embg = 0.d0
! .........................................................
      do 350 ip=1,npart
!
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 350
      emevcut(iz) = emevcut(iz) + evtwt
      if( .not. sigmacut ) emevt(iz) = emevt(iz) + evtwt
!
      if( pxycorr ) then
          call CANONICAL(pxc,pyc)  
      else
        pxc = pp(1)
        pyc = pp(2)
      end if
!
      pzc = pp(3)
      x = xp(1)
      y = xp(2)
!      cpx = pxc / pzc
      cpx = pxc
!      cpy = pyc / pzc
      cpy = pyc
      pm = SQRT( pxc**2 + pyc**2 + pzc**2 )
!
      embg = embg + evtwt*  pm / pmass
!
      emx(1) = emx(1) + evtwt* x
      emx(2) = emx(2) + evtwt* x * x
      emx(3) = emx(3) + evtwt* cpx
      emx(4) = emx(4) + evtwt* cpx * cpx
      emx(5) = emx(5) + evtwt* x * cpx
!
      emy(1) = emy(1) + evtwt* y
      emy(2) = emy(2) + evtwt* y * y
      emy(3) = emy(3) + evtwt* cpy
      emy(4) = emy(4) + evtwt* cpy * cpy
      emy(5) = emy(5) + evtwt* y * cpy
!
350   continue
! .........................................................
! Calculate x and y normalized emittances
! For statistical definition, cf. J. Buon, CERN 91-04, p. 48
!
      emitx(iz) = -999.
      emity(iz) = -999.
      n = emevcut(iz)
      if( n .lt. 2 ) go to 550
      betgam = embg / n    ! mean value of beta * gamma
!
      if( discorr ) then
        call DISP_CORR(n,emx,emy,xarg,yarg)
        if( xarg .gt. 0. )then
!          emitx(iz) = SQRT( xarg ) * betgam
           emitx(iz) = 1.d0/pmass * SQRT( xarg )
        end if
        if( yarg .gt. 0. )then
!          emity(iz) = SQRT( yarg ) * betgam
          emity(iz) = 1.d0/pmass * SQRT( yarg )
        end if
!
      else            ! no dispersion correction
        qx = emx(1) / n
        vx = emx(2)/n - qx*qx
        if(vx .lt. 0.) go to 410
        qpx = emx(3) / n
        vpx = emx(4)/n - qpx*qpx
        if(vpx .lt. 0.) go to 410
        cx = emx(5) / n - qx*qpx
        arg = vx * vpx - cx * cx
!
        if( arg .gt. 0. )then
!          emitx(iz) = SQRT( arg ) * betgam
          emitx(iz) = 1.d0/pmass * SQRT( arg )
        end if
c
410     continue
        qx = emy(1) / n
        vx = emy(2)/n - qx*qx
        if(vx .lt. 0.) go to 420
        qpy = emy(3) / n
        vpx = emy(4)/n - qpy*qpy
        if(vpx .lt. 0.) go to 420
        cx = emy(5) / n  - qx*qpy
        arg = vx * vpx - cx * cx
!
        if( arg .gt. 0. )then
!          emity(iz) = SQRT( arg ) * betgam
          emity(iz) = 1.d0/pmass * SQRT( arg )
        end if
      end if      ! end of test on discorr
! ---------------------------------------------------------
420   continue
!  Now accumulate final statistics for z emittance calculation
!  
      if( ipzcor .ne. 0 ) then   ! get Pz-amp correlation value
        if( emitx(iz).eq.0. .or. emity(iz).eq.0. ) go to 550
        betax = betgam * emx(2)/n  / emitx(iz)
        betay = betgam * emy(2)/n  / emity(iz)
        if( betax.gt.0. .and. betay.gt.0. ) then
!
          call PZ_AMP_COR(betax,betay,slp)
!
        end if  
      end if
!
      emitz(iz) = -999.
      do i=1,5
      emz(i) = 0.d0
      end do
! .........................................................
      do 450 ip=1,npart
!
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 450
!
      if( pxycorr ) then 
        call CANONICAL(pxc,pyc)  
      else
        pxc = pp(1)
        pyc = pp(2)
      end if
!
      if( ipzcor .eq. 1 ) then    ! use corrected Pz for emittance
        amp = SQRT( (pxc/pp(3))**2 + (pyc/pp(3))**2 + 
     &   (xp(1)/betax)**2 + (xp(2)/betay)**2 )
        pzc = pp(3) - slp*amp
!
      else if( ipzcor .eq. 2 ) then
        r2 = xp(1)**2 + xp(2)**2
        pt2 = pxc**2 + pyc**2
        pm2 = pmass**2
        amp = pt2/pm2 + ((pcharge*bfld(3))**2*r2)/(4.*pm2)
        amp = SQRT(amp)
        pzc = pp(3) - slp*amp
!
      else
        pzc = pp(3)
      end if
!
      energy = SQRT( pxc**2 + pyc**2 + pzc**2 + pmass**2 )
!
      if( refpar.ne.0  .and. diagref) then            ! with reference particle
        cpz = pzc - ppref(3)
        z = pp(3)/energy*cc * (tp - tpref)
      else                         ! no reference particle
        cpz = pzc
        z = pp(3)/energy*cc * tp
      end  if
!
      emz(1) = emz(1) + evtwt* z
      emz(2) = emz(2) + evtwt* z * z
      emz(3) = emz(3) + evtwt* cpz
      emz(4) = emz(4) + evtwt* cpz * cpz
      emz(5) = emz(5) + evtwt* z * cpz
!
450   continue
! .........................................................
! Calculate normalized z emittance
!
      emitz(iz) = -999.
      n = emevcut(iz)
      if( n .lt. 2 ) go to 550
!
      qx = emz(1) / n                 ! qx = <betaz*c*t> 
      vx = emz(2)/n - qx*qx
      if(vx .lt. 0.) go to 550
      qpz = emz(3) / n
      vpx = emz(4)/n - qpz*qpz
      if(vpx .lt. 0.) go to 550
      cx = emz(5) / n  - qx*qpz
      arg = vx * vpx - cx * cx
!
      if( arg .gt. 0.)then
        emitz(iz) = 1.d0/pmass * SQRT( arg )
      end if
c
! ---------------------------------------------------------
550   continue
      if( sigmacut ) then
        if( iflg .ne. 0 ) then      ! flag problem cases
          if(emitx(iz) .gt. 0.) emitx(iz) = -emitx(iz)
          if(emity(iz) .gt. 0.) emity(iz) = -emity(iz)
          if(emitz(iz) .gt. 0.) emitz(iz) = -emitz(iz)
        end if
!
!  Reset particle flag for particles in tails
        do ip=1,npart
!
        call READ_EVT(ip)
!
        if( ipflg .eq. 10 ) then    ! tail particle
          ipflg = 0
          call WRITE_EVT(ip)
        end if
!
        end do
      end if
!
600   continue        ! end of emittance plane loop
!
900   continue
      return
      end
c  ********************************************************
      subroutine FILL_HIST(izp)
C
c  Fills histogram arrays
C
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables: var
c
      if( nhist .eq. 0 ) go to 900
      call LOADVAR( izp )
c
      DO 110 I=1,NHISt
      if(ihdes(i) .ne. izp) go to 110
      val = var( ihvar(i) )
      call STATH(val,i)        ! accumulate statistics
      vhi = hxmin(i)  + nhbins(i)*hdx(i)
!
      if(val .lt. hxmin(i))then
        ihund(i) = ihund(i) + 1
        go to 110
      else if(val .gt. vhi)then
        ihovf(i) = ihovf(i) + 1
        go to 110
      else
        ihcon(i) = ihcon(i) + 1
        ip = lhpnt(i) + int( (val-hxmin(i))/hdx(i) )
        hval(ip) = hval(ip) + 1
      end if
!
110   continue
c
900   continue
      return
      end
c  *********************************************************
      subroutine FILL_POL(izp)
!
!  Fill polarization arrays at plane  IZP
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      do 100 i=1,nemit
      if( iemdes(i) .ne. izp ) go to 100
!
      if( izp .ge. 0 ) then       ! normal plane
        if( jpar.gt.0 .and. ABS(iptyp).eq.2 .and. spintrk.gt.0 ) then
          px = pol(1)
          py = pol(2)
          pz = pol(3)
          hel = HELICITY( pol )
        else
          px = 0.
          py = 0.
          pz = 0.
          hel = pol(3)    ! DECAY sets this when spintrk=F
        end if
!
      else                        ! production
        if( jpar.gt.0 .and. ABS(iptyp).eq.2 ) then
          px = polinit(1)
          py = polinit(2)
          pz = polinit(3)
          hel = HELICITY( polinit )
        end if
      end if
!
      polx(i) = polx(i) + px*evtwt
      poly(i) = poly(i) + py*evtwt
      polz(i) = polz(i) + pz*evtwt
      helic(i) = helic(i) + hel*evtwt
100   continue
!
      return
      end

c  ********************************************************
      SUBROUTINE FILL_SCAT(izp)
C
c  Fills scatterplot arrays
C
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables: var
c
!      dimension QST(20)
c
      IF(nscat .EQ. 0) GO TO 900
c
      call LOADVAR( izp )
c
      DO 200 I=1,Nscat
      if(isxdes(i).ne.izp .and. isydes(i).ne.izp) go to 200
!      idx = isxdes(i)
!      idy = isydes(i)
c
c  check if one of the plot quantities is downstream
c   if so, save in qst array
!      if(idx .gt. izp)then    ! x variable coming later
!        qst(i) = var( isyvar(i) )    ! save y variable from this place
!        go to 200
!      end if
!
!      if(idy .gt. izp)then    ! y variable coming later
!        qst(i) = var( isxvar(i) )  ! save x variable from this place
!        go to 200
!      end if
c
c  check if one of the plot quantities was upstream
!      if(idx .lt. izp)then    ! previously saved the x variable
!        vx = qst(i)
!        vy = var( isyvar(i) )
!      else if(idy .lt. izp)then  ! previously saved the y variable
!        vx = var( isxvar(i) )
!        vy = qst(i)
!      else                    ! both quantites are given at this plane
        vx = var( isxvar(i) )
        vy = var( isyvar(i) )
!      end if
c
      call STATSP(vx,vy,i)        ! accumulate statistics
!
c  check limits
      vxhi = sxmin(i) + nsxbin(i)*sdx(i)
      if(vx .lt. sxmin(i))then
        isxund(i) = isxund(i) + 1
        go to 200
      else if(vx .gt. vxhi)then
        isxovf(i) = isxovf(i) + 1
        go to 200
      end if
c
      vyhi = symin(i) + nsybin(i)*sdy(i)
      if(vy .lt. symin(i))then
        isyund(i) = isyund(i) + 1
        go to 200
      else if(vy .gt. vyhi)then
        isyovf(i) = isyovf(i) + 1
        go to 200
      end if
c
      iscon(i) = iscon(i) + 1
      ip = int( (vy-symin(i))/sdy(i) ) * nsxbin(i)
      ip = ip + int( (vx-sxmin(i))/sdx(i) ) + lspnt(i)
      sval(ip) = sval(ip) + 1
200   continue
c
900   continue
      return
      end

c  ********************************************************
      subroutine FILL_ZHIST
!
!  Fill z-history arrays
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables: var
!  Keep zhxval a function of particle # for print purposes
!
      call LOADVAR( 1 )
      zp = xp(3)
!
      do i=1,nzhist
      if( jpar .le. ABS(nzhpar(i)) ) then
        if( zp.ge.zhxcur(jpar,i) .and. zp.lt.zhxmax(i) ) then
!
            izhknt(jpar,i) = izhknt(jpar,i) + 1
            knt = izhknt(jpar,i)
            if( knt .gt. 70 ) knt = 70 
            zhxval( knt,jpar,i ) = zp
            zhyval( knt,jpar,i ) = var( izhyvar(i) )
! update current pointer
          zhxcur(jpar,i) = zhxcur(jpar,i) + zhdx(i)
        end if
!   
      end if
      end do
!
      return
      end
!  *********************************************************
      subroutine LBNL_EMIT(iz)
!
! Compute emittances, covariances, and related quantities
!  following the wisdom of the LBNL group
!   
! Adapted from G. Penn, 3 Nov 99
!         modified 21 Nov 99
!         modified 20 April 00
!         modified 12 Sept 00
!
!  if reference particle used, 
!     only particles of same type used in calculations
! Sigma cuts implemented, but different from EMITTANCE     
! "pzcorr" only affects longitudinal emittance calculation 
! "pxycorr" has no effect
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      dimension xtemp(3)
      dimension w(6),moments(6,6),mcovar(6,6),mtrans(4,4),avg(6)
      dimension xpc(3),ppc(3),mcovar2(6,6),sdev(6)
      dimension avglong(3),mlong(3,3),mclong(3,3),mc2(3,3)
      real*8 kappa,moments,mcovar,mtrans,mcovar2
      real*8 Lmech,Lcanon,Ldim,mlong,mclong,mc2
      character*8 axes(6)
      integer indperp(4),ishift(4),indx(6),indlong(3)
      integer usepar(mxpar),knt,iteste,ncut
      logical fset,longcorr,loopflag
!      LOGICAL FIRST
!
      data axes/'x','y','tp-tpref','Px','Py','Ep-Epref'/
      data iteste/20/
!      DATA FIRST/.TRUE./
!
!      IF( FIRST ) THEN
!        OPEN(UNIT=18,FILE='MTRANS.DAT',STATUS='UNKNOWN')
!        FIRST = .FALSE.
!      END IF
      fset = .false.
!
      if( refpar .ne. 0 ) then
        call READ_EVT(0)
        refenergy = DSQRT(ppref(3)**2 + pmass**2)
        reftime = tpref
        zloc = xp(3)
        call FIELD(xp,tp)
        bzaxis = bfld(3)
        charge = pcharge
        fset = .true.
      else
        refenergy = 0.
        reftime = 0.
        bzaxis = 0.
        charge = 0.
        zloc = 0.
      end if
!
      ishift(1) = 1
      ishift(2) = 2
      ishift(3) = 4
      ishift(4) = 5
!
      dtrans = 0.
      dlong = 0.
!
!      longcorr = pzcorr
      if( ipzcor .ne. 0 ) then
        longcorr = .true.
      else
        longcorr = .false.
      end if
!
      do ip=1,npart
         usepar(ip)=0
      end do
!
      loopflag = .true.  ! test for iterations
      knt = 0       ! count iterations
!
! ---------------------------------------------------------
!
      do while (loopflag)
!
! Particle cutting logic
        ncut = 0
        do 200 ip=1,npart
        if (usepar(ip) .lt. 0) go to 200
!
        call READ_EVT(ip)
!  energy in GeV
        energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2)
        if(ipflg .ne. 0) usepar(ip)=-1
        if((refpar .ne. 0) .and. (iptyp .ne. iptypref)) usepar(ip)=-2 
        if (usepar(ip) .lt. 0) go to 200
!
        if (.not. fset) then  ! find z and Bz
           zloc = xp(3)
           xtemp(1)=0.
           xtemp(2)=0.
           xtemp(3)=xp(3)
           call FIELD(xtemp,tp)
           bzaxis = bfld(3)
           charge = pcharge
           fset = .true.
        end if
!
        if (sigmacut) then  ! check if amplitudes too large
           Atrans = 0.
           if (dtrans .gt. 0.) then
              kappa = 0.5 * charge * bzaxis / avgpz
              do i=1,3
                 xpc(i)=xp(i)-avg(i)
                 ppc(i)=pp(i)-avg(i+3)
              end do
      Atrans = (1/pmass)*( (csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2)  +
     & (xpc(1)**2 + xpc(2)**2)* avgpz * csgamma +
     & 2*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) +
     & 2*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) )
              if (Atrans .gt. 2 * sig_cut**2 * dtrans) then
!  factor of two to make up for combining 4D at once 
                 usepar(ip)=-3 
                 ncut = ncut + 1
                 go to 200
              end if
           end if 
           if (dlong .gt. 0.) then
              ctime = tp - reftime - avglong(1) 
              cenergy = energy - refenergy - avglong(2) 
              if (longcorr) then
                 ctime = ctime - corr_t * (Atrans-avglong(3)) 
                 cenergy = cenergy - corr_E * (Atrans-avglong(3)) 
              end if
      Along = (cc/pmass) * ( ctime**2/delta +
     &   delta * (cenergy - alphalong * ctime / delta)**2 )
              if (Along .gt. sig_cut**2 * dlong) then
                 usepar(ip)=-4 
                 ncut = ncut + 1
                 go to 200
              end if
           end if 
        end if
200     continue        ! end of particle loop
!
!  calculate beam statistics
!     initialize
        avgpz = 0.
        do 250 i=1,6
        avg(i) = 0.
        sdev(i) = 0.
        do 250 j=1,6
        moments(i,j) = 0.
250     continue
        stev = 0.
!     first loop
        do 300 ip=1,npart
        if (usepar(ip) .lt. 0) go to 300
!
      call READ_EVT(ip)
      w(1) = xp(1)
      w(2) = xp(2)
      w(3) = tp - reftime
      w(4) = pp(1)
      w(5) = pp(2)
!  energy in GeV
      energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2)
      w(6) = energy - refenergy
!
!  Accumulate sums
      stev = stev + evtwt
      avgpz = avgpz + pp(3) * evtwt
      do i=1,6
      avg(i) = avg(i) + w(i) * evtwt
      do j=1,6
      moments(i,j) = moments(i,j) + w(i)*w(j) * evtwt
      end do
      end do
!
300     continue      ! end of first particle summation
!
      if( stev .eq. 0. ) go to 900
!
!  Get expectation values
      avgpz = avgpz / stev
      kappa = 0.5 * charge * bzaxis / avgpz
      do i=1,6
      avg(i) = avg(i) / stev          ! 1st moment (mean)
      do j=1,6
      moments(i,j) = moments(i,j) / stev      ! 2nd moment
      end do
      end do
!
!  Convert moments to covariances and standard deviations
      do i=1,6
      do j=1,6
         mcovar(i,j) = moments(i,j) - avg(i)*avg(j)
         mcovar2(i,j) = mcovar(i,j)
      end do
      if( mcovar(i,i) .lt. 0. ) go to 900     ! RCF 22 Nov 99
      sdev(i)=DSQRT(mcovar(i,i))
      end do
!
!   fill mtrans with x, y, Px, Py data (corr to index 1,2,4,5)
      do i=1,4
      do j=1,4
         mtrans(i,j) = mcovar(ishift(i),ishift(j))
      end do
      end do
!
! get 6-D emittance
      call LUDCMP(mcovar,6,6,indx,dfull)
!
      if( dfull .eq. 0. ) then        ! flag singular matrix
        dfull = -998.
        go to 680
      end if
!
      do i=1,6
      dfull = dfull * mcovar(i,i)
      end do
      if( dfull .gt. 0. ) then
        dfull = DSQRT( dfull )
      else
        dfull = -998.
        go to 680
      end if
! Normalization = c / ( (mc)^2 (mc^2))
      dfull = dfull * cc / pmass**3
! third pmass is actually mc^2 in MeV
680   continue
!
! get transverse emittance (square root of 4D emittance)
!      IF( J7REC.EQ.3 .OR. J7REC.EQ.52 ) THEN
!        WRITE(18,'(I5)')J7REC
!        DO I=1,4
!        WRITE(18,'(4E14.6)') (MTRANS(I,J),J=1,4)
!        END DO
!      END IF
      call LUDCMP(mtrans,4,4,indperp,dtrans)
!
      if( dtrans .eq. 0. ) then        ! flag singular matrix
        dtrans = -998.
        go to 780
      end if
!
!      IF( J7REC.EQ.3 .OR. J7REC.EQ.52 ) THEN
!        WRITE(18,'(I5)')J7REC
!        DO I=1,4
!        WRITE(18,'(4E14.6)') (MTRANS(I,J),J=1,4)
!        END DO
!      END IF
      do i=1,4
      dtrans = dtrans * mtrans(i,i)
      end do
      if( dtrans .gt. 0. ) then
!   start with 4D transverse emittance
        dtrans = DSQRT( dtrans )
!   take another square root to get pseudo 2D emittance
        dtrans = DSQRT( dtrans )
      else
        dtrans = -998.
        go to 780
      end if
!
! Normalization = 1 / (mc)
      dtrans = dtrans / pmass
780   continue
!
!
      Lmech = moments(1,5) - moments(2,4)
      Lcanon = Lmech + 0.5 * charge * bzaxis
     &   * (moments(1,1) + moments(2,2))
      if( dtrans.gt.0. ) then
        Ldim = (mcovar2(1,5) - mcovar2(2,4) + 
     &  0.5 * charge * bzaxis * (mcovar2(1,1) + mcovar2(2,2)))
     &  / ( 2 * pmass * dtrans)
        csalpha = - (mcovar2(1,4) + mcovar2(2,5)) 
     &     / (2 * pmass * dtrans)
        csbeta = (mcovar2(1,1) + mcovar2(2,2)) * avgpz 
     &     / (2 * pmass * dtrans)
        csgamma = (mcovar2(4,4) + mcovar2(5,5)) 
     &     / (2 * pmass * dtrans * avgpz)
        altgamma = (1 / csbeta) * (1 + csalpha**2 
     &   + (csbeta * kappa - Ldim)**2) 
      else
        longcorr = .false.
        Ldim = 0.
        csbeta = 0.
        csalpha = 0.
        csgamma = 0.
        altgamma = 0.
       end if
!
!  longitudinal beam parameters
!  initialize
      do 350 i=1,3
      avglong(i) = 0.
      do 350 j=1,3
      mlong(i,j) = 0.
 350     continue
      stev = 0.
!
!   second loop
        do 400 ip=1,npart
        if (usepar(ip) .lt. 0) go to 400
!
      call READ_EVT(ip)
      w(1) = tp - reftime
      energy = DSQRT(pp(1)**2 + pp(2)**2 + pp(3)**2 + pmass**2)
      w(2) = energy - refenergy
      if (longcorr) then
              do i=1,3
                 xpc(i)=xp(i)-avg(i)
                 ppc(i)=pp(i)-avg(i+3)
              end do
      Atrans = (1/pmass)*( (csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2)  +
     & (xpc(1)**2 + xpc(2)**2)* avgpz * csgamma +
     & 2*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) +
     & 2*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) )
         w(3) = Atrans
      else
         w(3) = 0.
      end if
!
!  Accumulate sums
      stev = stev + evtwt
      do i=1,3
      avglong(i) = avglong(i) + w(i) * evtwt
      do j=1,3
      mlong(i,j) = mlong(i,j) + w(i)*w(j) * evtwt
      end do
      end do
!
400   continue
!
      if( stev .eq. 0. ) go to 900
!
!  Get expectation values
      do i=1,3
      avglong(i) = avglong(i) / stev          ! 1st moment (mean)
      do j=1,3
      mlong(i,j) = mlong(i,j) / stev      ! 2nd moment
      end do
      end do
!
!  Convert moments to covariances and standard deviations
      do i=1,3
      do j=1,3
         mclong(i,j) = mlong(i,j) - avglong(i)*avglong(j)
         mc2(i,j) = mclong(i,j)
      end do
        if( mclong(i,i) .lt. 0. ) go to 900
      end do
!
! check some amplitude spread;  otherwise, no correlation
      sigmaA = DSQRT(mclong(3,3))
      if (sigmaA .eq. 0.) longcorr = .false. 
      if (longcorr) then
! get 3-D emittance
      call LUDCMP(mclong,3,3,indlong,dlong)
!
      if( dlong .eq. 0. ) then        ! flag singular matrix
        dlong = -998.
        go to 880
      end if
!
      do i=1,3
      dlong = dlong * mclong(i,i)
      end do
      if( dlong .gt. 0. ) then
        dlong = DSQRT( dlong / mc2(3,3))
      else
        dlong = -998.
        go to 880
      end if
      else
         dlong = mclong(1,1)*mclong(2,2)-mclong(1,2)**2
         if( dlong .gt. 0. ) then
           dlong = DSQRT(dlong)
         else
           dlong = -998.
           go to 880
         end if
      end if
!   normalize to c/(mc^2)
        dlong = dlong * cc / pmass
! pmass is actually mc^2 in MeV
 880       continue
!
      if(dlong .gt. 0.) then
      if (longcorr) then
        alphalong = (mc2(1,2) - mc2(1,3)*mc2(2,3)/mc2(3,3))
     &    / (dlong * pmass / cc)
        delta = (mc2(1,1) - mc2(1,3)**2/mc2(3,3))
     &    / (dlong * pmass / cc)
        corr_t = mc2(1,3)/mc2(3,3)
        corr_E = mc2(2,3)/mc2(3,3)
      else
        alphalong = mc2(1,2) / (dlong * pmass / cc)
        delta = mc2(1,1) / (dlong * pmass / cc)
        corr_t = 0.
        corr_E = 0.
      end if
      else
        alphalong = 0.
        delta = 0.
        corr_t = 0.
        corr_E = 0.
      end if
!
!  check whether should repeat loop
      knt = knt + 1
      if ((.not. sigmacut) .or. (knt .gt. iteste)) then
         loopflag=.false.
      else
! Test if tail cutting process has converged
! Stop if ncut = 0  (need to test changes in parameters?)
         if ((ncut .eq. 0) .and. (knt .gt. 1)) loopflag = .false.
      end if
!
      end do
!
!
! ---------------------------------------------------------
!  OUTPUTS
      write(2,*)
      write(2,210) j7rec
210   format(15('*'),t18,' PARTICLE MOMENTS, PLANE ',i5,
     &   t63,15('*'))        
      write(2,212)
212   format(t17,'Mean values',t35,'Standard Deviation')
      do i=1,6
        write(2,214) axes(i),avg(i),sdev(i)
214   format(t5,a8,t16,e13.6,t34,e13.6)
      end do
      write(2,216)
216   format(t2,' Reference Particle:',t30,'time',t45,'Energy',t60,'Pz')
      write(2,218) reftime,refenergy,ppref(3)
218   format(t24,e15.8,t39,e15.8,t54,e15.8)
      write(2,*)
!
      write(2,310)
310   format(t2,'Emittances: ',t22,'6D',t30,'Transverse',
     &   t42,'Corrected Longitudinal')
      write(2,320) dfull,dtrans,dlong
320   format(t16,3e12.4)
      write(2,'(a,e12.4)')'    eL*eT*eT = ',dtrans**2 * dlong
      write(2,330)
330   format(t2,'Angular Momentum: ',t21,'Mechanical',t33,'Canonical',
     &   t45,'Dimensionless = Lcanon/2mceT')
      write(2,340) Lmech,Lcanon,Ldim
340   format(t20,3e12.4)
      write(2,*)
      write(2,420) 
420   format(t10,'Modified Courant-Snyder Parameters:  ',
     &   'ASSUMES A CYLINDRICAL BEAM')
      write(2,432)
432   format(t5,'beta_perp',t17,'alpha_perp',t29,'gamma_perp')
      write(2,434) csbeta,csalpha,csgamma
434   format(t3,3e12.4)
      write(2,436) altgamma
436   format(t5,'consistency check:',t27,e12.4)
      write(2,*)
      write(2,438)
438   format('Long:',t7,'alphaL',t22,'delta',t37,'corr_t',t52,'corr_E',
     &   t67,'sigma Atrans')
      write(2,440) alphalong,delta,corr_t,corr_E,sigmaA
440   format(5e15.4)
      write(2,*)
      write(2,442)
442   format('MISC:',t7,'Z',t22,'N (weighted)',t37,'Bz on axis',
     &   t52,'kappa',t67,'<Pz>')
      write(2,444) zloc,stev,bzaxis,kappa,avgpz
444   format(5e15.4)
      write(2,*)
!
! Save some quantities for a summary table
        emit6(iz) = dfull
        emit4(iz) = dtrans
        emitl(iz) = dlong
        canang(iz) = Lcanon
        betper(iz) = csbeta
!
900   continue
      return
      end
!  *********************************************************
      subroutine LOADVAR(izp)
c
c  Load kinematic variables for diagnostics
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      if( izp .le. -1 ) then
        do i=1,3
        var(i) = xpint(i)
        var(i+3) = ppint(i)
        end do
        var(7) = tpint * cc
      else
        do i=1,3
        var(i) = xp(i)
        var(i+3) = pp(i)
        end do
        var(7) = tp * cc
        if( refpar.ne.0 .and. diagref ) then      ! redefine pz and t
          var(6) = var(6) - ppref(3)
          var(7) = var(7) - tpref*cc
        end if
      end if
!
      call VMAG(pp(1),var(8))
      var(9) = SQRT(var(8)**2 + pmass**2)      ! E
      gam = var(9) / pmass
      var(10) = pmass * (gam - 1.)             ! KE
      var(11) = pp(1) / pp(3)                  ! x'
      var(12) = pp(2) / pp(3)                  ! y'
      var(13) = ACOS( pp(3) / var(8) )         ! space angle
      var(14) = SQRT( xp(1)**2 + xp(2)**2 )    ! r
!
      if( .not.(xp(1).eq.0. .and. xp(2).eq.0.) ) then
        phi = ATAN2( xp(2),xp(1) )   ! phi
      else
        phi = 0.
      end if
!
      var(15) = phi  
      cph = COS( phi )
      sph = SIN( phi )
      var(16) = pp(1)*cph + pp(2)*sph
      var(17) = pp(2)*cph - pp(1)*sph
!
      x = xp(1)
      vlx = xp(2)*pp(3) - xp(3)*pp(2) 
      vly = xp(3)*pp(1) - x*pp(3) 
      vlz = x*pp(2) - xp(2)*pp(1) 
      vl2 = vlx*vlx + vly*vly + vlz*vlz 
      var(18) = vlz                 ! Lz
      var(19) = vl2                 ! L2
      var(20) = sarc
!
      var(21) = SQRT(pp(1)**2 + pp(2)**2)       ! pT
!      do i=1,3
!      var(20+i) = xpint(i)
!      var(23+i) = ppint(i)
!      end do
!      var(27) = tpint * cc
!
      if( betaperp .ne. 0. ) then                                   ! Palmer A2
        var(28) = (var(14)/betaperp)**2 + var(11)**2 + var(12)**2 
      else
        var(28) = 0.
      end if
c
      var(29) = var(14)**2                    ! r2
!
      hlab = 0.
      if( spin .and. jpar.gt.0 .and. ABS(iptyp).eq.2 .and.
     &  spintrk.gt.0 ) then
         hlab = HELICITY( pol )
      else
         hlab = pol(3)
      end if
      var(30) = hlab
!
      do i=1,3
      var(30+i) = bfld(i)
      var(33+i) = efld(i)
      var(36+i) = pol(i)
      end do
!
      var(40) = phaserf
!
      return
      end
!  *********************************************************
      subroutine MEANSD(stev,stem,xmean,xstand)
!
!  Compute mean and standard deviation of stored array
!
      implicit double precision (a-h,o-z)
      dimension stem(2,6),xmean(6),xstand(6)
!
      do i=1,6
      qx = stem(1,i) / stev
      xmean(i) = qx
      vx = stem(2,i)/stev - qx*qx
      if( vx .ge. 0. ) then
        xstand(i) = SQRT( vx )
      else
        xstand(i) = -999.
      end if
      end do
!
      return
      end
c  *********************************************************
      subroutine OUT_EMIT
c
c  Print out normalized emittances
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      real*8 n0
c
      if( nemit .eq. 0 ) go to 800
!
      write(2,*)
      write(2,40)
40    format(18('*'),t24,'Weighted NORMALIZED 2-D EMITTANCES',
     &   t61,18('*'))
      write(2,*)
      write(2,70)
70    format(t3,'des',t10,'No',t17,'N',t26,'emitX',t38,'emitY',
     +  t50,'emitZ',t61,'Tr',t69,'fET',t77,'fEL')
      write(2,75)
75    format(t27,'[m]',t39,'[m]',t51,'[m]')
!
      do 200 i=1,nemit
      if( i .eq. 1 ) then
        n0 = emevt(i)
        ex0 = emitx(i)
        ey0 = emity(i)
        ez0 = emitz(i)
      end if
!
      if( n0*ex0*ey0*ez0 .eq. 0. ) go to 800
!
      ftr = emevcut(i) / n0
      fet = 0.5 * (emitx(i)/ex0 + emity(i)/ey0) 
      fez = emitz(i) / ez0
      write(2,150)iemdes(i),emevt(i),emevcut(i),emitx(i),emity(i),
     &   emitz(i),ftr,fet,fez
150   format(i5,2f7.0,1x,3e12.4,3f8.3)
200   continue
!
800   continue
      if( ncovar .eq. 0 ) go to 900
      write(2,*)
      write(2,805)
805   format(15('*'),t20,'LBNL NORMALIZED DET(COV) EMITTANCES',
     &   t60,19('*'))
      write(2,*)
      write(2,810)
810   format(t2,'des',t12,'emitT',t25,'emitL',t38,'emit6D',
     &  t51,'betaT',t64,'Lcanon') 
      write(2,815)
815   format(t12,'[m]',t25,'[m]',t38,'[m^3]',t51,'[m]',t64,
     &  '[m GeV/c]')
      write(2,*)
!
      do 850 i=1,ncovar
      write(2,'(i5,2x,5e13.5)') icvdes(i),emit4(i),emitl(i),emit6(i),
     &  betper(i),canang(i) 
850   continue
!
900   continue
      return
      end

c  *********************************************************
      subroutine OUT_HIST
c
c  Generate histograms in character format
c
c  maximum histogram size is 50 bins
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables:
!             /dvars/ labvar 
c
      integer lowb,highb,xrange
      character*80 scr(25)
      character*10 cstr*10,fname
      character*6 csti
c
      if(nhist .le. 0) go to 900
c
      if( hcprn.ge.20 .and. hcprn.le.99 ) then    ! write out histogram contents
        write(fname,20) hcprn
20      format('for0',i2,'.dat')
        open(unit=hcprn,file=fname,status='unknown')
        write(hcprn,'(a)')' HISTOGRAM CONTENTS'
        do ih=1,nhist
        ip = lhpnt(ih) - 1
        nhb = nhbins(ih)
        write(hcprn,'(a,i5)')' Histogram # ',ih
        write(hcprn,'(2e12.4,2x,i4)') hxmin(ih),hdx(ih),nhbins(ih)
        write(hcprn,'(10i7)') (hval(ip+i),i=1,nhb)
        end do
        close(hcprn)
      end if
c  for each histogram, prepare virtual page in array SCR
c
      do 700 ih=1,nhist
      if( ihcon(ih) .eq. 0 ) go to 700
      ip = lhpnt(ih) - 1
      nhb = nhbins(ih)
c  clear out SCR
      do i=1,25
      do j=1,80
      scr(i)(j:j) = ' '
      end do
      end do
c       
c  find range of occupied x bins
      lowb = 0
      do i = 1,nhb
      ival = hval(ip+i)
      if(ival .gt. 0)then
      lowb = i
      go to 30
      end if
      end do
c     
30    continue
      highb = 51
      do i=nhb,1,-1
      ival = hval(ip+i)
      if(ival .gt. 0)then
      highb = i
      go to 40
      end if
      end do
c
40    continue
      xrange = highb - lowb + 1
c     
c  left and right vertical lines
      do i=2,24
      scr(i)(10:10) = '|'
      k = nhb + 11
      scr(i)(k:k) = '|'
      end do
c
c  top and bottom horizontal lines
      do i=1,nhb
      k = i + 10
      scr(1)(k:k) = '-'
      scr(25)(k:k) = '-'
      end do
c
c  histogram parameters
      scr(1)(63:69) = 'ICOOL  '
!      write(cstr,'(f10.2)')version
!      scr(1)(71:80) = cstr
      scr(1)(71:80) = version
      do i=62,80
      scr(2)(i:i) = '-'
      end do
      scr(3)(63:73) = 'histogram ='
      write(csti,'(i6)') ih
      scr(3)(74:79) = csti
      scr(4)(63:73) = 'variable  ='
      scr(4)(77:79) = labvar( ihvar(ih) )
      scr(5)(63:73) = 'z plane   ='
      ival = ihdes(ih)
      write(csti,'(i6)') ival
      scr(5)(74:79) = csti
c
      scr(7)(63:70) = 'lo lim ='
      write(cstr,'(e10.3)') hxmin(ih)
      scr(7)(71:80) = cstr
      scr(8)(63:70) = 'hi lim ='
      xhi = hxmin(ih) + real(nhb)*hdx(ih)
      write(cstr,'(e10.3)') xhi
      scr(8)(71:80) = cstr
      scr(9)(63:70) = 'step   ='
      write(cstr,'(e10.3)') hdx(ih)
      scr(9)(71:80) = cstr
c
      scr(11)(63:73) = 'contents  ='
      ival = ihcon(ih)
      write(csti,'(i6)') ival
      scr(11)(74:79) = csti
      scr(12)(63:73) = 'underflow ='
      ival = ihund(ih)
      write(csti,'(i6)') ival
      scr(12)(74:79) = csti
      scr(13)(63:73) = 'overflow  ='
      ival = ihovf(ih)
      write(csti,'(i6)') ival
      scr(13)(74:79) = csti
c
      scr(15)(63:73) = 'x range =  '
      write(csti,'(i6)') xrange
      scr(15)(74:79) = csti
c
      do i=62,80
      scr(16)(i:i) = '-'
      end do
      do i=1,16
      scr(i)(61:61) = '|'
      end do
c
c  y axis values
c
c   find maximum value
      mxval = 0
      do i=1,nhb
      ival = hval(ip+i)
      if( ival .gt. mxval ) mxval = ival
      end do
      if(hauto .and. mxval.gt.23)then
c  scale bin contents to maximum value
c  label y axis values
        do i=1,23
        val = real(mxval)*real(i)/23.
        write(cstr,'(f9.1)')val
        scr(25-i)(1:9) = cstr
        end do
      else
c  label y axis from 1 to 23
        do i=1,23
        write(csti,'(i6)') i
        scr(25-i)(3:8) = csti
        end do
        mxval = 23
      end if
c
      step = real(mxval) / 23.
c
c  fill in the histogram data
c   loop over columns
      do 200 i=1,nhb
      val = hval(ip+i)
      k = 10 + i
c  loop over rows
      do 150 j=1,23
      height = (j-1)*step
      if(val .gt. height)then
        scr(25-j)(k:k) = 'X'
      else
        go to 170
      end if
150   continue
c
170   continue
c  mark contents that are off scale when not using autoscaling
      if(int(val) .gt. mxval)scr(2)(k:k) = '^'
200   continue
c
!  Fit two plots per printed page
      if( ffcr .and. mod(ih,2).eq.1 ) then
        call FORMFEED
      else if( ffcr .and. mod(ih,2).eq.0 ) then
        write(2,*)
      end if
!
      do i=1,25
      write(2,'(a80)')scr(i)
      end do
700   continue
c
900   continue
      return
      end

c  *********************************************************
      subroutine OUT_POL
!
!  Print summary of polarization values at "emittance" planes
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 n
!
      write(2,*)
      write(2,40)
40    format(20('*'),t28,'Weighted POLARIZATION',t58,20('*'))
      write(2,50)
50    format(t3,'des',t12,'N',t25,'polX',t37,'polY',
     &  t49,'polZ',t60,'helicity')
!
      do 200 i=1,nemit
      n = emevt(i)
      if( n .lt. 1 ) go to 200
      px = polx(i) / n
      py = poly(i) / n
      pz = polz(i) / n
      hel = helic(i) / n
      write(2,150) iemdes(i),n,px,py,pz,hel
150   format(i5,f8.0,5x,4f12.4)
200   continue
!
      return
      end
!  *********************************************************
      subroutine OUT_RHIST
!
!  Generate R-history plots in character format
!
c  plots have 70 horizontal bins
c  plots have 23 vertical bins
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      integer*2 r,c
      character*80 scr(25)
      character*10 cstr
      character*6 csti
c
      if(nrhist .le. 0) go to 900
c  
      if( rhprin ) then    ! print out R-history contents
        write(2,'(a)')' R-HISTORY CONTENTS'
        do ih=1,nrhist
        write(2,'(a,i5)')' R-History # ',ih
        do iz=1,70
        izp = irhzmin(ih) + iz*irhdz(ih) - 1
        write(2,'(i4,4e12.4)') izp,(rhyval(iz,j,ih),j=1,4)
        end do
        end do
      end if
!
c  for each R-history, prepare virtual page in array SCR
c
      do 700 ih=1,nrhist
!
c  clear out SCR
      do i=1,25
      do j=1,80
      scr(i)(j:j) = ' '
      end do
      end do
c       
c  Left vertical line
      do i=1,24
      scr(i)(10:10) = '|'
      end do
      scr(12)(10:10) = '+'
      scr(12)(3:5) = labvar( irhyvar(ih) )
c
!  Horizontal line
      do i=10,80
      scr(24)(i:i) = '-'
      end do
      scr(24)(27:27) = '+'
      scr(24)(45:45) = '+'
      scr(24)(62:62) = '+'
      scr(25)(42:48) = 'Plane #'
!
c  Plot parameters
      scr(5)(3:7) = 'ICOOL'
!      write(cstr,'(f10.2)')version
!      scr(6)(1:7) = cstr(4:10)
      scr(6)(1:7) = version
      scr(8)(1:9) = 'R history'
      write(csti,'(i6)') ih
      scr(9)(1:5) = csti(2:6)
      do i=1,9
      scr(4)(i:i) = '-'
      scr(10)(i:i) = '-'
      end do
!
c  if rhauto = true => force all points to fill plot
c
      if(rhauto)then
        ymin = 1e9
        ymax = -1e9
c  Find maximum and minimum values
!
        do ip=1,4
        do 80 iz=1,70
        y = rhyval(iz,ip,ih)
        if( y .eq. 0. ) go to 75
        if( y .gt. ymax) ymax = y
        if( y .lt. ymin) ymin = y
75      continue
80      end do
        end do
!
      else                    ! use input file limits
        ymin = rhymin(ih)
        ymax = rhymax(ih)
      end if
!
      dy = (ymax-ymin) / 23.
!
!  Label minimum and maximum values on plot
      izp = irhzmin(ih)
      write(csti,'(i6)') izp
      scr(25)(5:10) = csti
      write(csti,'(i6)') izp+69*irhdz(ih)
      scr(25)(75:80) = csti
      write(cstr,'(e10.3)') ymin
      scr(24)(1:10) = cstr
      write(cstr,'(e10.3)') ymax
      scr(1)(1:10) = cstr
!
c  Fill in the R-history data
!
      if( dy .gt. 0. ) then
        do ip=1,4
!
        do 180 iz=1,70
        c = 10 + iz
        y = rhyval(iz,ip,ih)
        if( y .eq. 0. ) go to 175
        r = 24. - (y-ymin) / dy
        if( r .lt. 1 ) go to 175
        if( r .gt. 24 ) go to 175
!
        if( ip .eq. 1 ) then
          scr(r)(c:c) = '-'
        else if( ip .eq. 2 ) then
          scr(r)(c:c) = '+'
        else if( ip .eq. 3 ) then
          scr(r)(c:c) = 'm'
        else
          scr(r)(c:c) = 's'
        end if
!
175     continue
180     end do
        end do
!
      end if
!
!  Fit two plots per printed page
      if( ffcr .and. mod(ih,2).eq.1 ) then
        call FORMFEED
      else if( ffcr .and. mod(ih,2).eq.0 ) then
        write(2,*)
      end if
!
!  Write out the virtual page
!
      do i=1,25
      write(2,'(a80)')scr(i)
      end do
!
700   continue                    ! end of loop on R-histories
c
900   continue
      return
      end
c  *********************************************************
      subroutine OUT_SCAT
c
c  Generate scatterplots in character format
c
c  plots can have a maximum of 50 horizontal bins
c  plots can have a maximum of 23 vertical bins
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables:
!                /dvars/ labvar 
c
      integer*2 r,c
      character*80 scr(25)
      character*10 cstr
      character*6 csti
c
      if(nscat .le. 0) go to 900
c  
c  for each scatterplot, prepare virtual page in array SCR
c
      do 700 ih=1,nscat
      if( iscon(ih) .eq. 0 ) go to 700
      ip = lspnt(ih) - 1
      nxb = nsxbin(ih)
      nyb = nsybin(ih)
c  clear out SCR
      do i=1,25
      do j=1,80
      scr(i)(j:j) = ' '
      end do
      end do
c       
c  left and right vertical lines
      do i=1,nyb
      scr(25-i)(10:10) = '|'
      k = nxb + 11
      scr(25-i)(k:k) = '|'
      end do
c
c  top and bottom horizontal lines
      do i=1,nxb
      k = i + 10
      scr(25)(k:k) = '-'
      scr(24-nyb)(k:k) = '-'
      end do
c
c  scatterplot parameters
      scr(1)(63:69) = 'ICOOL  '
!      write(cstr,'(f10.2)')version
!      scr(1)(71:80) = cstr
      scr(1)(71:80) = version
      do i=62,80
      scr(2)(i:i) = '-'
      end do
      scr(3)(63:73) = 'plot      ='
      write(csti,'(i6)') ih
      scr(3)(74:79) = csti
      scr(4)(63:73) = 'var(x)    ='
      scr(4)(77:79) = labvar( isxvar(ih) )
      scr(5)(63:73) = 'dest(x)   ='
      ival = isxdes(ih)
      write(csti,'(i6)') ival
      scr(5)(74:79) = csti
      scr(6)(63:73) = 'var(y)    ='
      scr(6)(77:79) = labvar( isyvar(ih) )
      scr(7)(63:73) = 'dest(y)   ='
      ival = isydes(ih)
      write(csti,'(i6)') ival
      scr(7)(74:79) = csti
c     
      scr(9)(63:70) = 'lo x   ='
      write(cstr,'(e10.3)') sxmin(ih)
      scr(9)(71:80) = cstr
      scr(10)(63:70) = 'hi x   ='
      xhi = sxmin(ih) + real( nxb*sdx(ih) )
      write(cstr,'(e10.3)') xhi
      scr(10)(71:80) = cstr
      scr(11)(63:70) = 'step x ='
      write(cstr,'(e10.3)') sdx(ih)
      scr(11)(71:80) = cstr
      scr(12)(63:70) = 'lo y   ='
      write(cstr,'(e10.3)') symin(ih)
      scr(12)(71:80) = cstr
      scr(13)(63:70) = 'hi y   ='
      yhi = symin(ih) + real( nyb*sdy(ih) )
      write(cstr,'(e10.3)') yhi
      scr(13)(71:80) = cstr
      scr(14)(63:70) = 'step y ='
      write(cstr,'(e10.3)') sdy(ih)
      scr(14)(71:80) = cstr
c
      scr(16)(63:73) = 'contents  ='
      ival = iscon(ih)
      write(csti,'(i6)') ival
      scr(16)(74:79) = csti
      scr(17)(63:73) = 'under x   ='
      ival = isxund(ih)
      write(csti,'(i6)') ival
      scr(17)(74:79) = csti
      scr(18)(63:73) = 'over x    ='
      ival = isxovf(ih)
      write(csti,'(i6)') ival
      scr(18)(74:79) = csti
      scr(19)(63:73) = 'under y   ='
      ival = isyund(ih)
      write(csti,'(i6)') ival
      scr(19)(74:79) = csti
      scr(20)(63:73) = 'over y    ='
      ival = isyovf(ih)
      write(csti,'(i6)') ival
      scr(20)(74:79) = csti
c
      do i=62,80
      scr(21)(i:i) = '-'
      end do
      do i=1,21
      scr(i)(61:61) = '|'
      end do
c     
!  if sauto = true and maxval > 61,
!   then scale all contents to a maximum of 61
c
      if(sauto)then
        mxval = 0
        k = 0
c   find maximum value
        do i=1,nxb
        do j=1,nyb
        k = k + 1
        ival = sval(ip+k)
        if( ival .gt. mxval ) mxval = ival
        end do
        end do
        if(mxval .gt. 61)then
          k = 0
          do i=1,nxb
          do j=1,nyb
          k = k + 1
          sval(ip+k) = sval(ip+k)*61/mxval
          end do
          end do
        end if
      end if
c
c  fill in the scatterplot data
      k = 0
      do 200 j=1,nyb
      r = 25 - j
      do 150 i=1,nxb
      c = 10 + i
      k = k + 1
      ival = sval(ip+k)
!
      if(ival .eq. 0)then   ! use bin content to determine plot symbol
        scr(r)(c:c) = ' '
      else if(ival .eq. 1)then
        scr(r)(c:c) = '+'
      else if(ival.gt.1 .and. ival.lt.10)then
        scr(r)(c:c) = CHAR(ival + 48)
      else if(ival.gt.9 .and. ival.lt.36)then
        scr(r)(c:c) = CHAR(ival + 87)
      else if(ival.gt.35 .and. ival.lt.62)then
        scr(r)(c:c) = CHAR(ival + 29)
      else if(ival .gt. 61)then
        scr(r)(c:c) = '*'
      end if
c
150   continue
200   continue
c
!  Fit two plots per printed page
      if( ffcr .and. MOD(ih,2).eq.1 ) then
        call FORMFEED
      else if( ffcr .and. MOD(ih,2).eq.0 ) then
        write(2,*)
      end if
!
      do i=1,25         ! write out the array
      write(2,'(a80)')scr(i)
      end do
700   continue
c
900   continue
      return
      end

c  *********************************************************
      subroutine OUT_STATS
c
c  Compute and list statistics for histogram and scatter variables
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      real*8 kurtx,kurty
      real*8 m1x(20),m2x(20),m3x(20),m4x(20)
      real*8 m1y(20),m2y(20),m3y(20),m4y(20)
      real*8 n,sd_x(20),sd_y(20)
!
! ---------------------------------------------------------
      if( nhist .eq. 0 ) go to 232
!
      write(2,*)
      write(2,40)
40    format(20('*'),t28,'HISTOGRAM STATISTICS',t58,20('*'))
      write(2,100)
100   format('  ih',t7,'var',t13,'des',t22,'N',t31,'mean',t43,'sd',
     +  t55,'lo',t69,'hi')
c
      do 200 i=1,nhist
      n = st(1,i)
      if(n .ge. 1.)then
        m1x(i) = st(2,i) / n
        m2x(i) = st(3,i) / n
        m3x(i) = st(4,i) / n
        m4x(i) = st(5,i) / n
        vr = m2x(i) - m1x(i)*m1x(i)
        if(vr .ge. 0.)then
          sd_x(i) = sqrt(vr)
        else
          sd_x(i) = -998.
        end if
      else
        m1x(i) = -999.
        sd_x(i) = -999.
      end if
      write(2,150)i,ihvar(i),ihdes(i),n,m1x(i),sd_x(i),st(6,i),st(7,i)
150   format(2i4,i6,f11.1,4e12.4)
200   continue
! ---------------------------------------------------------
      write(2,*)
      write(2,205)
205   format(20('*'),t30,'HISTOGRAM MOMENTS',t58,20('*'))
      write(2,210)
210   format('  ih',t7,'var',t13,'des',t19,'Skewness',t31,'Kurtosis')
c
      do 230 i=1,nhist
      n = st(1,i)
      if( n .eq. 0 ) go to 230
      if( sd_x(i) .le. 0. ) go to 230
      skewx = m3x(i) - 3.*m1x(i)*m2x(i) + 2.*m1x(i)**3
      skewx = skewx / sd_x(i)**3
      kurtx = m4x(i) - 4.*m1x(i)*m3x(i) + 6.*m1x(i)**2*m2x(i) 
     &  -3.*m1x(i)**4
      kurtx = kurtx / sd_x(i)**4 - 3.
      write(2,220)i,ihvar(i),ihdes(i),skewx,kurtx
220   format(2i4,i6,2f12.4)
230   continue
! ---------------------------------------------------------
232   continue
      if( nscat .eq. 0 ) go to 900
!
      write(2,*)
      write(2,235)
235   format(20('*'),t28,'SCATTERPLOT STATISTICS',t58,20('*'))
      write(2,240)
240   format('  is',t6,'var',t12,'des',t16,'ax',t27,'N',t35,'mean'
     &  ,t47,'sd',t59,'lo',t71,'hi')
c
      do 300 i=1,nscat
      n = stsp(1,i)
      if(n .ge. 1.)then
        m1x(i) = stsp(2,i) / n
        m2x(i) = stsp(3,i) / n
        m3x(i) = stsp(4,i) / n
        m4x(i) = stsp(5,i) / n
        vr = m2x(i) - m1x(i)*m1x(i)
        if(vr .ge. 0.)then
          sd_x(i) = sqrt(vr)
        else
          sd_x(i) = -998.
        end if
      else
        m1x(i) = -999.
        sd_x(i) = -999.
      end if
      write(2,250)i,isxvar(i),isxdes(i),n,m1x(i),sd_x(i),
     +  stsp(6,i),stsp(7,i)
250   format(2i4,i6,2x,'x',2x,f12.1,4e12.4)
!  ........................................................
      if(n .ge. 1.)then
        m1y(i) = stsp(8,i) / n
        m2y(i) = stsp(9,i) / n
        m3y(i) = stsp(10,i) / n
        m4y(i) = stsp(11,i) / n
        vr = m2y(i) - m1y(i)*m1y(i)
        if(vr .ge. 0.)then
          sd_y(i) = sqrt(vr)
        else
          sd_y(i) = -998.
        end if
      else
        m1y(i) = -999.
        sd_y(i) = -999.
      end if
      write(2,270)i,isyvar(i),isydes(i),n,m1y(i),sd_y(i),
     +  stsp(12,i),stsp(13,i)
270   format(2i4,i6,2x,'y',2x,f12.1,4e12.4)
300   continue
! ---------------------------------------------------------
!
      write(2,*)
      write(2,305)
305   format(20('*'),t29,'SCATTERPLOT MOMENTS',t58,20('*'))
      write(2,310)
310   format('  is',t7,'vx',t13,'dx',t17,'vy',t23,'dy',
     &  t33,'r',t42,'skewx',t52,'kurtx',t62,'skewy',t72,'kurty')
c
      do 330 i=1,nscat
      n = stsp(1,i)
      if( n .eq. 0 ) go to 330
      if( sd_x(i).le.0. .or. sd_y(i).le.0. ) go to 330
!
      skewx = m3x(i) - 3.*m1x(i)*m2x(i) + 2.*m1x(i)**3
      skewx = skewx / sd_x(i)**3
      kurtx = m4x(i) - 4.*m1x(i)*m3x(i) + 6.*m1x(i)**2*m2x(i) 
     &  -3.*m1x(i)**4
      kurtx = kurtx / sd_x(i)**4 - 3.
!
      skewy = m3y(i) - 3.*m1y(i)*m2y(i) + 2.*m1y(i)**3
      skewy = skewy / sd_y(i)**3
      kurty = m4y(i) - 4.*m1y(i)*m3y(i) + 6.*m1y(i)**2*m2y(i) 
     &  -3.*m1y(i)**4
      kurty = kurty / sd_y(i)**4 - 3.
!
      fac = sd_x(i) * sd_y(i)
      if( fac .le. 0. ) go to 330
      r = (stsp(14,i) / n - m1x(i) * m1y(i)) / fac
!
      write(2,320)i,isxvar(i),isxdes(i),isyvar(i),isydes(i),
     &  r,skewx,kurtx,skewy,kurty
320   format(2i4,i6,i4,i6,f12.3,4f10.4)
330   continue
      write(2,*)
! ---------------------------------------------------------
900   continue
      return
      end
c  *********************************************************
      subroutine OUT_ZHIST
c
!  Generate Z-history plots in character format
c
c  plots have 70 horizontal bins
c  plots have 23 vertical bins
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!  Common variables:
!               /dvars/ labvar 
c
      integer*2 r,c
      character*80 scr(25)
      character*10 cstr
      character*6 csti
c
      if(nzhist .le. 0) go to 900
c  
      if( zhprin ) then    ! print out z-history contents
        write(2,'(a)')' Z-HISTORY CONTENTS'
        do ih=1,nzhist
        write(2,'(a,i5)')' Z-History # ',ih
        do ip=0,nzhpar(ih)
        write(2,'(a,i5)')' Particle # ',ip
        do iz=1,izhknt(ip,ih)
        write(2,'(i4,2e12.4)') iz,zhxval(iz,ip,ih),zhyval(iz,ip,ih)
        end do
        end do
        end do
      end if
!
c  for each z-history, prepare virtual page in array SCR
c
      do 700 ih=1,nzhist
!
c  clear out SCR
      do i=1,25
      do j=1,80
      scr(i)(j:j) = ' '
      end do
      end do
c       
c  Left vertical line
      do i=1,24
      scr(i)(10:10) = '|'
      end do
      scr(12)(10:10) = '+'
      scr(12)(3:5) = labvar( izhyvar(ih) )
c
!  Horizontal line
      do i=10,80
      scr(24)(i:i) = '-'
      end do
      scr(24)(27:27) = '+'
      scr(24)(45:45) = '+'
      scr(24)(62:62) = '+'
      scr(25)(45:45) = 'Z'
!
c  Scatterplot parameters
      scr(5)(3:7) = 'ICOOL'
!      write(cstr,'(f10.2)')version
!      scr(6)(1:7) = cstr(4:10)
      scr(6)(1:7) = version
      scr(8)(1:9) = 'Z history'
      write(csti,'(i6)') ih
      scr(9)(1:5) = csti(2:6)
      do i=1,9
      scr(4)(i:i) = '-'
      scr(10)(i:i) = '-'
      end do
!
      xmin = zhxmin(ih)
      xmax = zhxmax(ih)
!
c  if zhauto = true => force all points to fill plot
c
      if(zhauto)then
        dx = (xmax - xmin) / 70.
        xmaxp = xmax
!
        ymin = 1e9
        ymax = -1e9
c  Find maximum and minimum values
!
        do ip=0,nzhpar(ih)
        do 100 iz=1,izhknt(ip,ih)
        if( iz .gt. 70 ) go to 100
        y = zhyval(iz,ip,ih)
        if( y .gt. ymax) ymax = y
        if( y .lt. ymin) ymin = y
100     continue
        end do
!
      else                    ! use input file limits
        dx = zhdx(ih)
        xmaxp = xmin + 70 * dx
        ymin = zhymin(ih)
        ymax = zhymax(ih)
      end if
!
      dy = (ymax-ymin) / 23.
!
!  Label minimum and maximum values on plot
      write(cstr,'(e10.3)') xmin
      scr(25)(10:19) = cstr
      write(cstr,'(e10.3)') xmaxp
      scr(25)(71:80) = cstr
      write(cstr,'(e10.3)') ymin
      scr(24)(1:10) = cstr
      write(cstr,'(e10.3)') ymax
      scr(1)(1:10) = cstr
!
c  Fill in the z-history data
!
      if( dx .gt. 0. .and. dy .gt. 0. ) then
        do 250 ip=0,nzhpar(ih)
!
        do 200 iz=1,izhknt(ip,ih)
        if( iz .gt. 70 ) go to 200
        x = zhxval(iz,ip,ih)
        c = (x-xmin) / dx + 10
        if( c .lt. 11 ) c = 11
        if( c .gt. 80 ) c = 80
        y = zhyval(iz,ip,ih)
        r = 24. - (y-ymin) / dy
        if( r .lt. 1 ) go to 200
        if( r .gt. 24 ) go to 200
!
        if( ip .lt. 10 ) then
          scr(r)(c:c) = CHAR(ip + 48)
        else
          scr(r)(c:c) = CHAR(ip + 55)
        end if
!
200     continue
!
250     continue
c
      end if
!
!  Fit two plots per printed page
      if( ffcr .and. mod(ih,2).eq.1 ) then
        call FORMFEED
      else if( ffcr .and. mod(ih,2).eq.0 ) then
        write(2,*)
      end if
!
!  Write out the virtual page
!
      do i=1,25
      write(2,'(a80)')scr(i)
      end do
!
700   continue                    ! end of loop on z-histories
c
900   continue
      return
      end
!  *********************************************************
      subroutine PZ_AMP_COR(betax,betay,slp)
!
!  Determine Pz - amplitude correlation
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 v(5)
!
      do i=1,5
      v(i) = 0.
      end do
!
      do 200 ip=1,npart
!
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 200
!
      if( pxycorr ) then 
        call CANONICAL(pxc,pyc)  
      else
        pxc = pp(1)
        pyc = pp(2)
      end if
!
      if( ipzcor .eq. 1 ) then
        xpr = pxc / pp(3)
        ypr = pyc / pp(3)
        amp = SQRT( xpr**2 + ypr**2 + (xp(1)/betax)**2
     &   + (xp(2)/betay)**2 )
!
      else if( ipzcor .eq. 2 ) then
        r2 = xp(1)**2 + xp(2)**2
        pt2 = pxc**2 + pyc**2
        pm2 = pmass**2
        amp = pt2/pm2 + ((pcharge*bfld(3))**2*r2)/(4.*pm2)
        amp = SQRT(amp)
!
      else
        amp = 0.
      end if
!
      v(1) = v(1) + evtwt * pp(3)
      v(2) = v(2) + evtwt * amp
      v(3) = v(3) + evtwt * amp*amp
      v(4) = v(4) + evtwt * amp*pp(3)
      v(5) = v(5) + evtwt     
200   continue
!
! Get parameters of least square fit to straight line (cf. NR p. 656)
      s = v(5)
      sx = v(2)
      sy = v(1)
      sxx = v(3)
      sxy = v(4)
      del = s*sxx - sx*sx
      if( del .ne. 0. ) then
!        off = (sxx*sy - sx*sxy) / del
        slp = (s*sxy - sx*sy) / del
      end if
!
      return
      end
!  *********************************************************
      subroutine RHISTORY 
!
!  Compute and load R-history data into array
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      do 600 ih=1,nrhist
      if( j7rec.lt.irhzmin(ih) .or. j7rec.gt.irhzmax(ih) ) go to 600
      izp = (j7rec - irhzmin(ih))/irhdz(ih) + 1
!
      stev = 0.
      stem = 0.
      stem2 = 0.
      stemlo = 1e20
      stemhi = -1e20
!
!  Accumulate statistics; get min and max values
      do 100 ip=1,npart
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 100
      call LOADVAR(j7rec)
      v = var( irhyvar(ih) )
      stev = stev + 1.
      stem = stem + v
      stem2 = stem2 + v*v
      if( v .lt. stemlo ) stemlo = v
      if( v .gt. stemhi ) stemhi = v
100   continue
!
!  Calculate mean, standard deviation; load array
      if( stev .gt. 0. ) then
        qx = stem / stev
        vx = stem2/stev - qx*qx
      else
        qx = 0.
        vx = 0.
      end if    
      if( vx .gt. 0. ) then
        sd = SQRT( vx )
      else
        sd = 0.
      end if
!
      rhyval(izp,1,ih) = stemlo
      rhyval(izp,2,ih) = stemhi
      rhyval(izp,3,ih) = qx
      rhyval(izp,4,ih) = sd
!
600   continue
!            
      return
      end
!  *********************************************************
      subroutine SIGMA_CUT(xmean,xstand,ncut)
!
!  Flag particle with {x...Pz} in tails of distributions
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 xmean(6),xstand(6),v(6)
      integer ncut(6)
!
      do i=1,6
      ncut(i) = 0
      end do
!
      do 200 ip=1,npart
!
      call READ_EVT(ip)
      if( ipflg .ne. 0 ) go to 200
!
! Correct Px and Py for vector potential
      if( pxycorr ) then   
        call CANONICAL(pxc,pyc)  
      else
        pxc = pp(1)
        pyc = pp(2)
        end if
!
      energy = SQRT( pxc**2 + pyc**2 + pp(3)**2 + pmass**2 )
      z = pp(3) / energy * cc * tp
!
      v(1) = xp(1)
      v(2) = xp(2)
      v(3) = z
      v(4) = pxc
      v(5) = pyc
      v(6) = pp(3)
!
      do i=1,6
      arg = v(i) - xmean(i)     ! difference from mean
      if( ABS(arg).gt.sig_cut*xstand(i)
     &   .and. xstand(i).gt.0. ) then
        ncut(i) = ncut(i) + 1
        ipflg = 10
        call WRITE_EVT(ip)
        go to 200
      end if
      end do
!
200   continue
!
      return
      end

c  *********************************************************
      subroutine STATH(v,ih)
c
c  Accumulate statistics on the histogram variables
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      st(1,ih) = st(1,ih) + 1.
      st(2,ih) = st(2,ih) + v
      st(3,ih) = st(3,ih) + v*v
      st(4,ih) = st(4,ih) + v*v*v
      st(5,ih) = st(5,ih) + v*v*v*v
      if(v .lt. st(6,ih)) st(6,ih) = v
      if(v .gt. st(7,ih)) st(7,ih) = v
c
      return
      end
c  *********************************************************
      subroutine STATSP(vx,vy,is)
c
c  Accumulate statistics on the scatterplot variables
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      stsp(1,is) = stsp(1,is) + 1.
c
      stsp(2,is) = stsp(2,is) + vx
      stsp(3,is) = stsp(3,is) + vx*vx
      stsp(4,is) = stsp(4,is) + vx*vx*vx
      stsp(5,is) = stsp(5,is) + vx*vx*vx*vx
      if(vx .lt. stsp(6,is)) stsp(6,is) = vx
      if(vx .gt. stsp(7,is)) stsp(7,is) = vx
c
      stsp(8,is) = stsp(8,is)   + vy
      stsp(9,is) = stsp(9,is)   + vy*vy
      stsp(10,is) = stsp(10,is) + vy*vy*vy
      stsp(11,is) = stsp(11,is) + vy*vy*vy*vy
      if(vy .lt. stsp(12,is)) stsp(12,is) = vy
      if(vy .gt. stsp(13,is)) stsp(13,is) = vy
c
      stsp(14,is) = stsp(14,is) + vx * vy
c
      return
      end
