      program ECALC9F
!
! Modification of ECALC9 using an input command file
!
! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
c
c  Program Ecalc9, version 2
c
c  LBL emittance calculations, revised 6/10/2001
c  reads particle data file in for009.dat format
c  now has same capabilities as built-in icool diagnostics
c     sigma cuts enabled, and can turn off energy-amplitude correlations
c     also stores processing information in header of output file
c  extra diagnostics beyond what is built in to icool
c     correlations of energy and time with transverse amplitude
c     dispersion calc's:  x, y, radial, also r^2
c
c  canonical momentum calc only works when reference particle is used
c     else assumes Bz=0
c  calculates # particles within transverse emittance cut
c     for two different emittances at once
c  option of adjusting particle time modulus RF period, 
c     centered at reference particle
c  reads input file for009.dat
c  writes to output file ecalc9.out
c
c
      implicit double precision (a-h,o-z)
      parameter ( npmax = 150000 )
      parameter ( version = 2.09 )
c
      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(8)
      real*8 cc,ee,pi,twopi,avagn,relcl,mass
c
      common/CONTROL/sigma_cut,rffreq,transcuta,transcutb,
     &  longcut,bzaxis,refcharge,refmass,zloc,npart,ncurrent,
     &  pzcorr,wrnew9
      real*8 sigma_cut,rffreq,transcuta,transcutb,longcut
      real*8 bzaxis,refcharge,refmass,zloc
      integer npart,ncurrent
      logical pzcorr
      logical wrnew9
!
      common/PART/xpar(3,npmax),ppar(3,npmax),tpar(npmax),
     &  epar(npmax),evtwtpar(npmax),bfldpar(3,npmax),
     &  efldpar(3,npmax),sarcpar(npmax),polpar(3,npmax),
     &  iptypar(npmax),usepar(npmax),ievtpar(npmax),
     &  ipnumpar(npmax),ipflgpar(npmax),jsrgpar(npmax)
      real*8 xpar,ppar,tpar,epar,evtwtpar,bfldpar,efldpar
      real*8 sarcpar,polpar
      integer iptypar,usepar,ievtpar,ipnumpar
      integer ipflgpar,jsrgpar
c
      integer useptype,ioc,i,ip
      real*8 pzmax,pzmin,tp,pmass
      real*8 Bp(3),Ep(3),phasep,polp(3),reftime,refenergy,energy
      integer t_ievt,t_ipnum,t_iptyp,t_ipflg
      real*8 t_evtwt,t_xp(3),t_pp(3),count,t_tp
      logical loopflag, flag
!      character*1 letter
      integer*4 nregion
      character*80 title
      character*20 response
!
! const
      data cc/2.9979246d8/,ee/0.29979246d0/,pi/3.1415926d0/,
     +  avagn/6.0221367d23/,relcl/2.8179409d-13/,
     &  twopi/6.2831853d0/
!      data mass/0.51099906d-3, 0.10565839d0, 0.13956995d0, 0.493677d0,
!     +  0.93827231d0/  ! GeV/c^2
      data mass/0.51099892d-3, 0.105658367d0, 0.13957018d0,   ! GeV/c^2
     &  0.493677d0, 0.93827201d0,  
     &  1.8756128d0,  2.809413d0,  6.535366/    ! d,He3,Li7 
!
!  obtain from user the following (includes default option):
!    type of particle (0 = all types)
!    Bz on axis
!    Pzmax
!    Pzmin
!    transverse emittance cutoff 1
!    transverse emittance cutoff 2
!    longitudinal emittance cutoff
!    RF frequency
!    sigma cut
!    subtract out correlations?
!    
c   default values
      useptype=0
      pzmin = 0.0d0
      pzmax = 0.0d0
      transcuta=0.0d0
      transcutb=0.0d0
      longcut=0.0d0
      rffreq=0.0d0
      sigma_cut = 0.0d0
      pzcorr = .true.
      response = ' '
!
      write(*,'(a,f6.2)')' Ecalc9f version ',version
! vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
!
      open(unit=1,file='ecalc9f.inp',status='old',iostat=ioc)
	if( ioc .ne. 0 ) then
	  write(*,*) 'Error opening ECALC9F.INP'
	  stop
	end if
!
      read(1,'(a80)') title
      read(1,*,iostat=ioc) useptype
      read(1,*,iostat=ioc) pzmin,pzmax
      read(1,*,iostat=ioc) transcuta,transcutb
      read(1,*,iostat=ioc) longcut
      read(1,*,iostat=ioc) rffreq
!      read(1,*) fper
      read(1,*,iostat=ioc) sigma_cut
      read(1,*,iostat=ioc) pzcorr
      read(1,*,iostat=ioc) wrnew9
      if( ioc .lt. 0 ) then         ! EOF
        wrnew9 = .false.
      end if
!
      write(*,'(a80)') title
      write(*,*)'   particle type: ',useptype
      write(*,*)'   pzmin/pzmax: ',pzmin,'/',pzmax
      write(*,*)'   transverse cuts: ',transcuta,transcutb
      write(*,*)'   longitudinal cut: ',longcut
      write(*,*)'   RF frequency (MHz): ',rffreq
!      write(*,*)'   RF offset: ',fper
      write(*,*)'   sigma cut: ',sigma_cut
      if (pzcorr) then
         write(*,*)'   subtract out amplitude correlation'
      else
         write(*,*)'   do not subtract out amplitude correlation'
      end if
      write(*,*)'   wrnew9 = ',wrnew9
c
! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      open(unit=3,file='for009.dat',status='unknown',iostat=ioc)
      read(3,'(a80)') title     ! title card on beam output file
      read(3,*)
      read(3,*)
!
      if( wrnew9 ) then
        open(unit=9,file='ecalc9f.cut',status='unknown')
      end if
c
      open(unit=2,file='ecalc9f.dat',status='unknown',iostat=ioc)
      write(2,*)title
      write(2,'(a,f6.2)')' output from Ecalc9f version ',version
      write(2,*)' settings: '
      write(2,*)'   particle type: ',useptype
      write(2,*)'   pzmin/pzmax: ',pzmin,'/',pzmax
      write(2,*)'   transverse cuts: ',transcuta,transcutb
      write(2,*)'   longitudinal cut: ',longcut
      write(2,*)'   RF frequency (MHz): ',rffreq
!      write(2,*)'   RF offset: ',fper
      write(2,*)'   sigma cut: ',sigma_cut
      if (pzcorr) then
         write(2,*)'   subtract out amplitude correlation'
      else
         write(2,*)'   do not subtract out amplitude correlation'
      end if
      write(2,*)'   wrnew9 = ',wrnew9
      write(2,*)
      write(2,30)
30    format(t1,'regn #',t10,'Z',t22,'Bz',
     &  t34,'eperp',t46,'elong',t58,'e6D',t70,'Ldim',
     &  t82,'Pzavg',t94,'beta',t106,'alpha',
     &  t118,'betaL',t130,'alphaL',t142,'n0',t154,'n1',t166,'n2',
     &  t178,'Lcan(eVs)',t190,'sigmaE',t202,'sigmaT',
     &  t214,'corrE',t226,'corrT',t238,'sigmaE_c',
     &  t250,'xavg',t262,'yavg',t274,'Dx',t286,'Dy',
     &  t298,'Dr',t310,'Dr2')
c
      loopflag = .true.
      read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,nregion,
     &     t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp
!      BP(3) = 0.
      if (ioc .ne. 0) then
         loopflag = .false.
      end if
      do while (loopflag)   ! loop over regions
!
         flag = .true.
         reftime = 0.d0
         refenergy = 0.d0
         if (useptype .eq. 0) then
            refcharge = ee
            refmass = mass(2)
         else
            refcharge = SIGN( ee, DBLE(useptype) )
            refmass = mass( ABS(useptype) )
         end if
         zloc = 0.0d0
         bzaxis = 0.0d0
         ip = 0
         ncurrent = nregion
         if (.not. loopflag) flag = .false.
!
         do while (flag)        ! loop over particles in region
          ip = ip + 1
          do i=1,3
            xpar(i,ip) = t_xp(i)
            ppar(i,ip) = t_pp(i)
            bfldpar(i,ip) = Bp(i)
            efldpar(i,ip) = Ep(i)
            polpar(i,ip) = polp(i)  
          end do
!
          tpar(ip) = t_tp
          iptypar(ip) = t_iptyp
          evtwtpar(ip) = t_evtwt
          ievtpar(ip) = t_ievt
          ipnumpar(ip) = t_ipnum
          ipflgpar(ip) = t_ipflg
          jsrgpar(ip) = nregion
          sarcpar(ip) = phasep
!
          usepar(ip) = 0
          if(t_ipflg .ne. 0) usepar(ip) = -1
          if((useptype .ne. 0) .and. (t_iptyp .ne. useptype)) 
     &           usepar(ip) = -2
          if((pzmax .gt. 0.0d0) .and. (t_pp(3) .gt. pzmax))
     &           usepar(ip) = -3
          if(t_pp(3) .lt. pzmin)
     &           usepar(ip) = -4
          if(t_ievt .eq. 0) then
             bzaxis = Bp(3)
             usepar(ip) = -5
             reftime = t_tp
             refenergy = DSQRT(t_pp(3)**2 + 
     &            mass(ABS(t_iptyp))**2)
             zloc = t_xp(3)
             if (useptype .eq. 0) then
                refcharge = SIGN( ee, DBLE(t_iptyp) )
                refmass = mass( ABS(t_iptyp) )
             end if   
          end if
c
          read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,nregion,
     &         t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp
!      BP(3) = 0.
          if (ioc .ne. 0) loopflag = .false.
          if (nregion .ne. ncurrent) flag = .false.
          if (.not. loopflag) flag = .false.
         end do   ! end loop over particles 
         npart = ip
c
c loop over particles to adjust time and energy if necessary
         do ip=1,npart
            if (usepar(ip) .eq. 0) then
               if (refenergy .eq. 0.) zloc = xpar(3,ip)
               pmass = mass( ABS(iptypar(ip)) )
               energy = DSQRT(ppar(1,ip)**2 + ppar(2,ip)**2 + 
     &              ppar(3,ip)**2 + pmass**2)  
               epar(ip) = energy 
               tp = tpar(ip) - reftime
               if (rffreq .ne. 0.d0) then
                  per = 1.0d-6 / rffreq
                  if( refcharge .lt. 0.d0 ) then
                    toff = 0.5d0 * per
                  else
                    toff = 0.d0
                  end if
!                  toff = toff + fper * per 
                  count = 0.5d0 + (tp+toff) * rffreq * 1.0d+6
                  if (count .lt. 0.d0) count = count - 1.
                  tp = tp - INT(count) * per
!
               end if
               tpar(ip) = tp
            else
               epar(ip)=0.d0
            end if
         end do
c
c Compute emittances and related quantities for this region
         call EMIT
c
      end do ! loop over regions
c
      close(3)
      close(2)
c
      end
c
c  *********************************************************
      subroutine EMIT
c
      implicit double precision (a-h,o-z)
      parameter ( npmax = 150000 )
c
      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(8)
      real*8 cc,ee,pi,twopi,avagn,relcl,mass
c
      common/CONTROL/sigma_cut,rffreq,transcuta,transcutb,
     &  longcut,bzaxis,refcharge,refmass,zloc,npart,ncurrent,
     &  pzcorr,wrnew9
      real*8 sigma_cut,rffreq,transcuta,transcutb,longcut
      real*8 bzaxis,refcharge,refmass,zloc
      integer npart,ncurrent
      logical pzcorr,wrnew9
!
      common/PART/xpar(3,npmax),ppar(3,npmax),tpar(npmax),
     &  epar(npmax),evtwtpar(npmax),bfldpar(3,npmax),
     &  efldpar(3,npmax),sarcpar(npmax),polpar(3,npmax),
     &  iptypar(npmax),usepar(npmax),ievtpar(npmax),
     &  ipnumpar(npmax),ipflgpar(npmax),jsrgpar(npmax)
      real*8 xpar,ppar,tpar,epar,evtwtpar,bfldpar,efldpar
      real*8 sarcpar,polpar
      integer iptypar,usepar,ievtpar,ipnumpar
      integer ipflgpar,jsrgpar
c
      integer nscut,iteste,ip,n,i,j,knt,iptyp
      real*8 evtwt,tp,energy,pmass
      real*8 xp(3),pp(3),xpc(3),ppc(3)
      real*8 w(6),moments(6,6),mcovar(6,6),mtrans(4,4)
      real*8 m6(6,6),m3(3,3),mlong(3,3),avglong(3)
      real*8 avg(6),dfull,dtrans,Lmech,Lcanon,Ldim
      real*8 stev,csbeta,csalpha,csgamma,avgpz
      real*8 Along,Atrans,sigmaA,sigmaEc,altgamma
      real*8 ampsq, kappa,dlong,corr_t,corr_E
      real*8 alphalong,etrans,elong,ntotal,ncut0,ncut1,ncut2
      integer*4 indperp(4),ishift(4),indx(6),indlong(3)
      real*8 delta,betalong,ctime,cenergy,mclong(3,3)
      real*8 Dx,Dy,Dr,Dr2,avgR,avgER,avgER2
      logical longcorr,loopflag
!
      data iteste/20/
!
      ishift(1) = 1
      ishift(2) = 2
      ishift(3) = 4
      ishift(4) = 5
!
      dtrans = 0.d0
      dlong = 0.d0
!
      longcorr = pzcorr
!
      loopflag = .true.  ! test for iterations
      knt = 0       ! count iterations
      nrec1 = 0
!
! ---------------------------------------------------------
!
      do while (loopflag)
!
        if( wrnew9 .and. nrec1.ne.0 ) then
          do i=1,nrec1
          backspace(9)
          end do
          write(9,'(3i7)') 0,1,jsrgpar(1)      ! reference particle
          nrec1 = 1
        end if
!
! Particle cutting logic
        nscut = 0     ! # particles removed by sigma cuts this round
        do 200 ip=1,npart
        if (usepar(ip) .ne. 0) go to 200
!
        do i=1,3
          xp(i) = xpar(i,ip)
          pp(i) = ppar(i,ip)
        end do
!
        tp = tpar(ip)
        energy = epar(ip)
        iptyp = iptypar(ip)
        evtwt = evtwtpar(ip)
        pmass = mass( ABS(iptyp) )
!        pcharge = SIGN( ee, REAL(iptyp) )
        if (sigma_cut .gt. 0.d0) then  ! check if amplitudes too large
           Atrans = 0.d0
           if (dtrans .gt. 0.d0) then
              kappa = 0.5d0 * refcharge * bzaxis / avgpz
              do i=1,3
                 xpc(i)=xp(i)-avg(i)
                 ppc(i)=pp(i)-avg(i+3)
              end do
      Atrans = (1d0/pmass)*( (csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2)  +
     & (xpc(1)**2 + xpc(2)**2)* avgpz * csgamma +
     & 2d0*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) +
     & 2d0*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) )
              if (Atrans .gt. 2d0 * sigma_cut**2 * dtrans) then
!  factor of two to make up for combining 4D at once
                 usepar(ip) = knt
                 nscut = nscut + 1
                 go to 200
              end if
           end if
           if (dlong .gt. 0.) then
              ctime = tp - avglong(1)
              cenergy = energy - 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. sigma_cut**2 * dlong) then
                 usepar(ip) = knt 
                 nscut = nscut + 1
                 go to 200
              end if
           end if
        end if
 200         continue        ! end of particle loop
!
!  calculate beam statistics
!  initialize
      avgpz = 0.d0
      avgR = 0.d0
      avgER = 0.d0
      avgER2 = 0.d0
      do 250 i=1,6
         avg(i) = 0.
      do 250 j=1,6
         moments(i,j) = 0.
250   continue
      stev = 0.d0
      n = 0
!
!   first loop
      do 300 ip=1,npart
!
        if(usepar(ip) .ne. 0) go to 300
!
        do i=1,3
          xp(i) = xpar(i,ip)
          pp(i) = ppar(i,ip)
        end do
!
        tp = tpar(ip)
        energy = epar(ip)
        iptyp = iptypar(ip)
        evtwt = evtwtpar(ip)
        pmass = mass( ABS(iptyp) )
!        pcharge = SIGN( ee, REAL(iptyp) )
!
      n = n + 1
      stev = stev + evtwt
      w(1) = xp(1)
      w(2) = xp(2)
      w(3) = tp 
      w(4) = pp(1)
      w(5) = pp(2)
      w(6) = energy
!
!  Accumulate sums
      avgpz = avgpz + pp(3) * evtwt
      avgR = avgR + DSQRT(xp(1)**2 + xp(2)**2) * evtwt
      avgER = avgER + energy * DSQRT(xp(1)**2 + xp(2)**2) * evtwt
      avgER2 = avgER2 + energy * (xp(1)**2 + xp(2)**2) * 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
!
      if (stev .eq. 0.d0) go to 900
!
!  Get expectation values
      ntotal = stev
      avgpz = avgpz / stev
      kappa = 0.5d0 * refcharge * bzaxis / avgpz
      avgR = avgR / stev
      avgER = avgER / stev
      avgER2 = avgER2 / stev
      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)
            m6(i,j) = mcovar(i,j)
         end do
         if( mcovar(i,i) .lt. 0.d0 ) go to 900    
      end do
!
! Energy dispersions
      Dx = 0.d0
      Dy = 0.d0
      Dr = 0.d0
      Dr2 = 0.d0
      if (mcovar(6,6) .gt. 0.) then
         Dx = mcovar(1,6) * avg(6) / mcovar(6,6)
         Dy = mcovar(2,6) * avg(6) / mcovar(6,6)
         Dr = (avgER - avgR*avg(6)) * avg(6) / mcovar(6,6)
         Dr2 = (avgER2 - (mcovar(1,1)+mcovar(2,2))*avg(6)) 
     &        * avg(6) / mcovar(6,6)
      end if
!
! get 6-D emittance
      call LUDCMP(m6,6,6,indx,dfull)
!
      if( dfull .eq. 0.d0 ) then        ! flag singular matrix
         dfull = -998.
         go to 680
      end if
!
      do i=1,6
         dfull = dfull * m6(i,i)
      end do
      if( dfull .ge. 0. ) then
         dfull = DSQRT( dfull )
      else
         dfull = -998.
      end if
! Normalization = c / ( (mc)^2 (mc^2)) 
      dfull = dfull * cc / refmass**3
! third refmass is actually mc^2 in MeV
680   continue
!
!
!   fill mtrans with x, y, Px, Py data (corr to index 1,2,4,5)
      ishift(1) = 1
      ishift(2) = 2
      ishift(3) = 4
      ishift(4) = 5
      do i=1,4
      do j=1,4
         mtrans(i,j) = mcovar(ishift(i),ishift(j))
      end do
      end do
!
! get transverse emittance (square root of 4D emittance)
      call LUDCMP(mtrans,4,4,indperp,dtrans)
!
      if( dtrans .eq. 0.d0 ) then        ! flag singular matrix
        dtrans = -998.
        go to 780
      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 / refmass
780   continue
!
!
      kappa = 0.5d0 * refcharge * bzaxis / avgpz
      Lmech = moments(1,5) - moments(2,4)
      Lcanon = Lmech + 0.5d0 * refcharge * bzaxis
     &   * (moments(1,1) + moments(2,2))
      if( dtrans.gt.0. ) then
        Ldim = (mcovar(1,5) - mcovar(2,4) +
     &     0.5d0 * refcharge * bzaxis * (mcovar(1,1) + mcovar(2,2)))
     &     / ( 2d0 * refmass * dtrans)
        csalpha = - (mcovar(1,4) + mcovar(2,5))
     &     / (2d0 * refmass * dtrans)
        csbeta = (mcovar(1,1) + mcovar(2,2)) * avgpz
     &     / (2d0 * refmass * dtrans)
        csgamma = (mcovar(4,4) + mcovar(5,5))
     &     / (2d0 * refmass * dtrans * avgpz)
        altgamma = (1d0 / csbeta) * (1d0 + csalpha**2
     &     + (csbeta * kappa - Ldim)**2)
      else
        Ldim = 0.d0
        csbeta = 0.d0
        csalpha = 0.d0
        csgamma = 0.d0
        altgamma = 0.d0
       end if
!
!  Now accumulate final statistics for z emittance calculation
!      subtracts out amplitude correlation, neglecting
!         any correlation with canonical angular momentum
!  
      do 350 i=1,3
      avglong(i) = 0.
      do 350 j=1,3
      mlong(i,j) = 0.
350   continue
      n = 0.
      stev = 0.d0
!
!  second loop
      do 400 ip=1,npart
!
        if(usepar(ip) .ne. 0) go to 400
!
        do i=1,3
          xp(i) = xpar(i,ip)
		xpc(i) = xp(i) - avg(i)
          pp(i) = ppar(i,ip)
		ppc(i) = pp(i) - avg(i+3)
        end do
!
        tp = tpar(ip)
        energy = epar(ip)
        iptyp = iptypar(ip)
        evtwt = evtwtpar(ip)
        pmass = mass( ABS(iptyp) )
!        pcharge = SIGN( ee, REAL(iptyp) )
!
      n = n + 1.
      stev = stev + evtwt
      w(1) = tp
      w(2) = energy
      ampsq = (1d0/pmass)*((csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2)  +
     & (xpc(1)**2 + xpc(2)**2) * csgamma * avgpz +
     & 2d0*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) +
     & 2d0*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) )
      w(3) = ampsq
!
!  Accumulate sums
      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
!
!  Get expectation values
      ntotal = stev
      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)
         m3(i,j) = mclong(i,j) 
      end do
	if( mclong(i,i) .lt. 0.d0 ) go to 900    
      end do
!
!     check some amplitude spread;  otherwise, no correlation
      sigmaA = DSQRT(mclong(3,3))
      if (sigmaA .eq. 0.d0) longcorr = .false.
      if (longcorr) then
! get 3-D emittance
         call LUDCMP(m3,3,3,indlong,dlong)
         if( dlong .eq. 0.d0 ) then ! flag singular matrix
            dlong = -998.
            go to 880
         end if
         do i=1,3
            dlong = dlong * m3(i,i)
         end do
         if( dlong .gt. 0. ) then
            dlong = DSQRT(dlong / mclong(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.d0 ) then
           dlong = DSQRT(dlong)
         else
           dlong = -998.
           go to 880
         end if
      end if
!
!   normalize to c/(mc^2)
        dlong = dlong * cc / refmass
! refmass is actually mc^2 in MeV
880   continue
!
      if (dlong .gt. 0.) then
         if (longcorr) then
            alphalong = (mclong(1,2) - mclong(1,3)*mclong(2,3)/
     &                   mclong(3,3)) / (dlong * refmass / cc)
            delta = (mclong(1,1) - mclong(1,3)**2/mclong(3,3))
     &           / (dlong * refmass / cc)
            corr_t = mclong(1,3)/mclong(3,3)
            corr_E = mclong(2,3)/mclong(3,3)
         else
            alphalong = mclong(1,2) / (dlong * refmass / cc)
            delta = mclong(1,1) / (dlong * refmass / cc)
            corr_t = 0.d0
            corr_E = 0.d0
         end if
      else
        alphalong = 0.d0
        delta = 0.d0
        corr_t = 0.d0
        corr_E = 0.d0
      end if
c  define longitudinal beta roughly like transverse, <Pz> sigma_z/sigma_pz
      betalong = delta * cc * avgpz**3 / (refmass**2 + avgpz**2)
c
c  hard cuts in emittance
      ncut0 = 0.d0
      ncut1 = 0.d0
      ncut2 = 0.d0
      do 600 ip=1,npart
!
	  if(usepar(ip) .ne. 0) go to 600
!
        do i=1,3
          xp(i) = xpar(i,ip)
		xpc(i) = xp(i) - avg(i)
          pp(i) = ppar(i,ip)
		ppc(i) = pp(i) - avg(i+3)
        end do
!
        tp = tpar(ip)
        energy = epar(ip)
        iptyp = iptypar(ip)
        evtwt = evtwtpar(ip)
        pmass = mass( ABS(iptyp) )
!        pcharge = SIGN( ee, REAL(iptyp) )
!
	if (dtrans .gt. 0.d0) then
       ampsq = (1d0/pmass)*((csbeta/avgpz)*(ppc(1)**2 + ppc(2)**2)  +
     & (xpc(1)**2 + xpc(2)**2) * csgamma * avgpz +
     & 2d0*csalpha*(xpc(1)*ppc(1) + xpc(2)*ppc(2)) +
     & 2d0*(csbeta*kappa-Ldim)*(xpc(1)*ppc(2) - xpc(2)*ppc(1)) )
       etrans = ampsq
      else
	 ampsq = 0.
	 etrans = 0.
	end if 
      ctime = tp - avglong(1) - corr_t * (ampsq - avglong(3))
      cenergy = energy - avglong(2) - corr_E * (ampsq - avglong(3))
!
      if (dlong .gt. 0.d0) then
	 elong = (cc/pmass) * (ctime**2/delta + 
     &   delta * (cenergy - alphalong * ctime / delta)**2)
	else
	 elong = 0.d0
	end if
      ncut0 = ncut0 + evtwt
!
      if((etrans .lt. transcuta) .or. (transcuta .eq. 0.d0)) then
        if ((elong .lt. longcut) .or. (longcut .eq. 0.)) then 
          ncut1 = ncut1 + evtwt
          nrec1 = nrec1 + 1
          if( wrnew9 ) then
            write(9,'(3i7)')ievtpar(ip),ipnumpar(ip),jsrgpar(ip)
          end if
        end if
      end if
! 
      if((etrans .lt. transcutb) .or. (transcutb .eq. 0.d0)) then
       if ((elong .lt. longcut) .or. (longcut .eq. 0.)) 
     &   ncut2 = ncut2 + evtwt
      end if 
!
600   continue
!
!  check whether should repeat loop
      knt = knt + 1
      if ((sigma_cut .le. 0.) .or. (knt .gt. iteste)) then
         loopflag=.false.
      else
! Test if tail cutting process has converged
! Stop if nscut = 0  (need to test changes in parameters?)
         if ((nscut .eq. 0) .and. (knt .gt. 1)) loopflag = .false.
      end if
!
      end do            ! end of loop on loopflag?
!
!
! ---------------------------------------------------------
!555   continue
!
      sigmaEc = DSQRT(mclong(2,2) - corr_E**2 * mclong(3,3))
!  OUTPUTS
      write(2,150)ncurrent,zloc,bzaxis,dtrans,dlong,dfull,
     &   ldim,avgpz,csbeta,csalpha,betalong,alphalong,ncut0,
     &   ncut1,ncut2,Lcanon*10.E9/cc,DSQRT(mcovar(6,6)),
     &   DSQRT(mcovar(3,3)),corr_E,corr_t,sigmaEc,
     &   avg(1),avg(2),Dx,Dy,Dr,Dr2
150   format(i6,26e12.4)
!
900   continue
c
      end
c
c  *********************************************************
      subroutine LUDCMP(a,n,np,indx,d)
! 
! Make LU decomposition of matrix
!
      implicit double precision (a-h,o-z)
      INTEGER n,np,indx(n),NMAX
! ---------------------------------------------------------
! RCF modified
      REAL*8 d,a(np,np),TINY
! ---------------------------------------------------------
      PARAMETER (NMAX=500,TINY=1.0e-20)
      INTEGER i,imax,j,k
! ---------------------------------------------------------
! RCF modified
      REAL*8 aamax,dum,sum,vv(NMAX)
! ---------------------------------------------------------
      d=1.
      do 12 i=1,n
        aamax=0.
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.) then ! singular matrix in ludcmp
          d = 0.
          go to 900
        end if
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.
        do 16 i=j,n
          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.)a(j,j)=TINY
        if(j.ne.n)then
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
!
900   continue
      return
      END
