      program EMITCALC

c

c  Program EMITCALC, version 1.0

c  dated 3/16/2004

c

c  calculates many beam parameters and amplitude cuts

c  main new feature:  the 3 invariant emittances

c     note that the order in which they are listed may not

c     be as expected, especially if coupling is strong

c

c  inputs control parameters from 'emitcalc.inp'

c  reads particle data file in ICOOL for009.dat format

c  calculations are performed on blocks of particles 

c     having identical region # and Z

c  uses particle weights

c  most calculations assume axisymmetry

c  uses transverse position, canonical transverse momentum, energy and time

c

c  options:

c     broad particle cuts (particle type, pzmin, pzmax)

c     several models for calculating vector potential in axisymmetric field

c        set up to be able to add your own analytic model for vector potential

c        0:  assume constant B0, taken from input

c        1:  use Bz at particle plus reference particle

c            if reference particle not available, use model 2

c        2:  calculate using particle data, assuming axisymmetric 

c        3:  algebraic expansion from input file

c     may remove linear dependence of transverse correlations on E, t

c        (though implemented in a very poor way)

c     may remove nonlinear dependence of E, t, on transverse amplitude

c     may cut tails of distribution

c     may overlap bunches based on RF period (2 options for doing this)

c        0:  no overlap

c        1:  center about reference particle, if available (else model 2)

c        2:  use weighting to find proper centered phase

c     

c  outputs some beam averages, plus:

c     number of particles within hard amplitude cuts

c     3 invariant emittances, no axisymmetry assumption

c        note:  it is possible for emittances to be out of order

c     transverse, longitudinal emittance

c     transverse and longitudinal beta, alpha parameters

c     x and y component of energy dispersion

c     nonlinear correlation parameters, E and t with transverse amplitude

c     sigma_E

c     asymmetry test (if axisymmetric, should be very small number)

c

c  units:  meters, GeV, GeV/c, Tesla 

c     note <Lcanon> given in GeV-m/c, multiply by 3.33 to get in eV-s

c     phase space area is (particle rest mass) * c * emittance

c     RF frequencies are in MHz

c

c  note that e_6D=e_1*e_2*e_3 exactly, same as LBL definition except

c     for option to remove some nonlinear correlations

c

c  ************************************************************

c

      IMPLICIT NONE

c

      integer npmax

      parameter ( npmax = 100000 )

c

      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5)

      real*8 cc,ee,pi,twopi,avagn,relcl,mass

c

      common/CONTROL/sigma_cut,rffreq,transcut,longcut,

     &  bzaxis,refcharge,refmass,reftime,refenergy,zloc,

     &  pzmax,pzmin,curr_region,vp_model,rf_model,useptype,

     &  ampl_corr,disp_corr

      real*8 sigma_cut,rffreq,transcut,longcut

      real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc

      real*8 pzmax,pzmin

      integer curr_region,vp_model,rf_model,useptype

      logical ampl_corr,disp_corr

c

      common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax),

     &  epar(npmax),bzpar(npmax),evtwtpar(npmax),

     &  iptypar(npmax),usepar(npmax),npart

      real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar

      integer iptypar,usepar,npart

c

      integer i,ip,ioc

      real*8 Bp(3),Ep(3),phasep,polp(3),energy

      integer t_ievt,t_ipnum,t_iptyp,t_ipflg

      real*8 t_evtwt,t_xp(3),t_pp(3),count,t_tp,store_bz

      logical loopflag, flag, found_ref

      integer iregion, store_rf_model

!

! 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

!

      call SETUP

!

      store_bz = bzaxis

      store_rf_model = rf_model

      loopflag = .true.

      read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,iregion,

     &     t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp

      if (ioc .ne. 0) then

         loopflag = .false.

      end if

      do while (loopflag)       ! loop over regions

!     

         flag = .true.

         found_ref = .false.

         reftime = 0.

         refenergy = 0.

         if (useptype .eq. 0) then

            refcharge = ee

            refmass = mass(2)

         else

            refcharge = SIGN( ee, REAL(useptype) )

            refmass = mass( ABS(useptype) )

         end if

         bzaxis = store_bz

         rf_model = store_rf_model

         zloc = 0.0

         ip = 0

         curr_region = iregion

         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)  

            end do

!

            tpar(ip) = t_tp

            bzpar(ip) = Bp(3)

            evtwtpar(ip) = t_evtwt

            iptypar(ip) = t_iptyp

            if(t_ipflg .eq. 0) then

               usepar(ip) = 0

               energy = DSQRT(t_pp(1)**2 + t_pp(2)**2 + 

     &              t_pp(3)**2 + mass(ABS(iptypar(ip)))**2)  

            else

               usepar(ip) = -1

               energy = 0.

            end if

            epar(ip) = energy

            if((useptype .ne. 0) .and. (t_iptyp .ne. useptype)) 

     &           usepar(ip) = -2

            if((pzmax .gt. 0.0) .and. (t_pp(3) .gt. pzmax))

     &           usepar(ip) = -3

            if(t_pp(3) .lt. pzmin)

     &           usepar(ip) = -4

            if(t_ievt .eq. 0) then

               found_ref = .true.

               usepar(ip) = -5

               if (vp_model .ne. 0) then

                  bzaxis = Bp(3)

               end if

               zloc = t_xp(3)

               reftime = t_tp

               refenergy = energy

               if (useptype .eq. 0) then

                  refcharge = SIGN( ee, REAL(t_iptyp) )

                  refmass = mass( ABS(t_iptyp) )

               end if   

            end if

            if ((.not. found_ref) .and. (usepar(ip) .eq. 0)) then

               zloc = t_xp(3)

               reftime = t_tp

               if (useptype .eq. 0) then

                  refcharge = SIGN( ee, REAL(t_iptyp) )

                  refmass = mass( ABS(t_iptyp) )

               end if   

            end if

c     

            read(3,*,iostat=ioc)t_ievt,t_ipnum,t_iptyp,t_ipflg,iregion,

     &           t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp

            if (ioc .ne. 0) loopflag = .false.

            if (iregion .ne. curr_region) flag = .false.

            if (.not. loopflag) flag = .false.

!

         end do                 ! end loop over particles 

         npart = ip

         if ((.not. found_ref) .and. (rf_model .eq. 1)) rf_model = 2 

c     

c Compute emittances and related quantities for this region

         call EMITTANCE

c

      end do                    ! loop over regions

c

      close(3)

      close(2)

c

      end

c

c  *********************************************************

      subroutine EMITTANCE

c

      IMPLICIT NONE

c

      integer npmax,nloopmax

      real*8 min_conf,rf_scale

      parameter ( npmax = 100000 )

      parameter ( nloopmax = 20 )

      parameter ( min_conf = 0.3)

      parameter ( rf_scale = 1.E6)

c

      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5)

      real*8 cc,ee,pi,twopi,avagn,relcl,mass

c

      common/CONTROL/sigma_cut,rffreq,transcut,longcut,

     &  bzaxis,refcharge,refmass,reftime,refenergy,zloc,

     &  pzmax,pzmin,curr_region,vp_model,rf_model,useptype,

     &  ampl_corr,disp_corr

      real*8 sigma_cut,rffreq,transcut,longcut

      real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc

      real*8 pzmax,pzmin

      integer curr_region,vp_model,rf_model,useptype

      logical ampl_corr,disp_corr

c

      common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax),

     &  epar(npmax),bzpar(npmax),evtwtpar(npmax),

     &  iptypar(npmax),usepar(npmax),npart

      real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar

      integer iptypar,usepar,npart

c

      integer nscut,i,j,ip,iptyp,knt

      real*8 phasespace(6,npmax),transverse(4,npmax)

      real*8 phaseavg(6),phasecorr(6,6),transavg(4),transcorr(4,4)

      real*8 vector3(3),vector4(4),longavg(3),longcorr(3,3)

      real*8 vpx,vpy,bzp,xp(3),pp(3),pcharge,Lcanon,avgpz

      real*8 ampl_trans(npmax),ampl_long(npmax)

      real*8 wtot,weight,save_wtot,disp_t,disp_E,avg_E

      real*8 e_perp,beta_perp,alpha_perp,Ldim,gamma1

      real*8 e_long,beta_long,alpha_long,delta_long

      real*8 invar_emit(3),sigma_Aperp,avg_Aperp,corr_t,corr_E

      real*8 w_aperture,disp_x,disp_y,check_axisym

      real*8 store_corr_t,store_corr_E

      real*8 temp_t,temp_E,conf

      real*8 cut_perp,cut_long

      logical cut_flag,store_ampl_corr

      real*8 TIMESHIFT

!

      cut_flag = .true.         ! test for iterations

      knt = 1                   ! count iterations

      store_ampl_corr = ampl_corr

!

      do while (cut_flag)

!

c     initializations

         ampl_corr = store_ampl_corr

         if (vp_model .eq. 2) then

            call CALC_VP(bzaxis)

         end if

         if ((rf_model .gt. 0) .and. (rffreq .ne. 0.)) then

            call CALC_TC(rffreq*rf_scale,reftime,conf)

            if (conf .lt. min_conf) then

               write(*,*)' difficulty calculating bunch center'

            end if

         end if

c     

c     calculate averages, canonical momentum

c     shift time, energy by reference values

c     store new values in phasespace

c     

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               do i=1,3

                  xp(i) = xpar(i,ip)

                  pp(i) = ppar(i,ip)

               end do

               bzp = bzpar(ip)

!     

               call GET_VP(bzp,xp(1),xp(2),vpx,vpy)

!     

               iptyp = iptypar(ip)

               pcharge = SIGN( ee, REAL(iptyp) )

c     

c   *****   use canonical momentum, subtract reference time and energy *****

c

               phasespace(1,ip) = xp(1)

               phasespace(2,ip) = pp(1) + pcharge * vpx 

               phasespace(3,ip) = xp(2)

               phasespace(4,ip) = pp(2) + pcharge * vpy 

               if (rf_model .eq. 0) then

                  phasespace(5,ip) = tpar(ip) - reftime

               else

                  phasespace(5,ip) = 

     &                 TIMESHIFT(tpar(ip)-reftime,rffreq*rf_scale)

               end if

               phasespace(6,ip) = epar(ip) - refenergy

            end if

         end do

!     

!     particle weights and some averages

         wtot = 0.

         avgpz = 0.

         Lcanon = 0.

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               wtot = wtot + weight

               avgpz = avgpz + weight * ppar(3,ip)

               Lcanon = Lcanon + weight * 

     &              (phasespace(1,ip)*phasespace(4,ip)-

     &              phasespace(3,ip)*phasespace(2,ip))

            end if

         end do

         if (knt.eq.1) save_wtot=wtot ! total particle weight before cuts

         if (wtot>0.) then

            avgpz = avgpz / wtot

            Lcanon = Lcanon / wtot

         end if

!

!     6D correlations matrix 

         do i=1,6

            phaseavg(i)=0.

            do j=1,6

               phasecorr(i,j)=0.

            end do

         end do

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               do i=1,6

                  phaseavg(i)=phaseavg(i)+weight*phasespace(i,ip)

                  do j=1,6

                     phasecorr(i,j)=phasecorr(i,j)+

     &                    weight*phasespace(i,ip)*phasespace(j,ip)

                  end do

               end do

            end if

         end do

         if (wtot>0.) then

            do i=1,6

               phaseavg(i)=phaseavg(i)/wtot

            end do

            do i=1,6

               do j=1,6

                  phasecorr(i,j)=phasecorr(i,j)/wtot -

     &                 phaseavg(i)*phaseavg(j)

               end do

            end do

         end if

         avg_E = phaseavg(6) + refenergy

!     

!     *****   for transverse parameters,

!     *****   subtract out linear dispersion with energy, time

         do i=1,4

            if (phasecorr(5,5) .gt. 0.) then

               disp_t = phasecorr(i,5)/phasecorr(5,5)

            else

               disp_t = 0.

            end if

            if (phasecorr(6,6) .gt. 0.) then

               disp_E = phasecorr(i,6)/phasecorr(6,6)

            else

               disp_E = 0.

            end if

c     save some parameters for output

            if (i .eq. 1) disp_x = disp_E * avg_E

            if (i .eq. 3) disp_y = disp_E * avg_E

            if (.not. disp_corr) then

               disp_t = 0.

               disp_E = 0.

            end if

            do ip=1,npart

               if (usepar(ip) .eq. 0) then

                  transverse(i,ip)=phasespace(i,ip)

     &                 - disp_t*(phasespace(5,ip)-phaseavg(5))

     &                 - disp_E*(phasespace(6,ip)-phaseavg(6))

               end if

            end do

         end do

!     

!     averages, correlation matrix

         do i=1,4

            transavg(i)=0.

            do j=1,4

               transcorr(i,j)=0.

            end do

         end do

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               do i=1,4

                  transavg(i)=transavg(i)+weight*transverse(i,ip)

                  do j=1,4

                     transcorr(i,j)=transcorr(i,j)+

     &                    weight*transverse(i,ip)*transverse(j,ip)

                  end do

               end do

            end if

         end do

         if (wtot>0.) then

            do i=1,4

               transavg(i)=transavg(i)/wtot

            end do

            do i=1,4

               do j=1,4

                  transcorr(i,j)=transcorr(i,j)/wtot -

     &                 transavg(i)*transavg(j)

               end do

            end do

         end if

!     

!     calculate transverse emittance, beta, alpha, assuming axisymmetry

         call AXISYM(refmass,avgpz,transcorr,e_perp,

     &        beta_perp,alpha_perp,Ldim,check_axisym)

!     calculate single particle amplitudes

         if (beta_perp .gt. 0.) then

            gamma1 = (1.+alpha_perp**2+Ldim**2)/beta_perp

         else 

            gamma1 = 0.

         end if

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               do i=1,4

                  vector4(i)=transverse(i,ip)-transavg(i)

               end do

               ampl_trans(ip) = (1/refmass)*((beta_perp/avgpz)*

     &              (vector4(2)**2+vector4(4)**2) + 

     &              gamma1*avgpz*(vector4(1)**2+vector4(3)**2) + 

     &              2.*alpha_perp*(vector4(1)*vector4(2) + 

     &              vector4(3)*vector4(4)) - 2.*Ldim*

     &              (vector4(1)*vector4(4)-vector4(2)*vector4(3)))

            end if

         end do

!

!     *****   now for longitudinal parameters,

!     *****   subtract out nonlinear correlation with transverse amplitude

         do i=1,3

            longavg(i)=0.

            do j=1,3

               longcorr(i,j)=0.

            end do

         end do

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               vector3(1) = phasespace(5,ip)

               vector3(2) = phasespace(6,ip)

               vector3(3) = ampl_trans(ip)

               do i=1,3

                  longavg(i)=longavg(i)+weight*vector3(i)

                  do j=1,3

                     longcorr(i,j)=longcorr(i,j)+

     &                    weight*vector3(i)*vector3(j)

                  end do

               end do

            end if

         end do

         if (wtot>0.) then

            do i=1,3

               longavg(i)=longavg(i)/wtot

            end do

            do i=1,3

               do j=1,3

                  longcorr(i,j)=longcorr(i,j)/wtot -

     &                 longavg(i)*longavg(j)

               end do

            end do

         end if

         avg_Aperp = longavg(3)

         sigma_Aperp = longcorr(3,3)

         if (sigma_Aperp .le. 0.) ampl_corr = .false.

!

         if (sigma_Aperp .gt. 0.) then

            corr_t = longcorr(1,3)/sigma_Aperp

            corr_E = longcorr(2,3)/sigma_Aperp

         else

            corr_t = 0.

            corr_E = 0.

         end if

c     save correlation parameters for output

         store_corr_t = corr_t

         store_corr_E = corr_E

         if (.not. ampl_corr) then

            corr_t = 0.

            corr_E = 0.

         end if

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               phasespace(5,ip)=phasespace(5,ip)

     &              - corr_t*(ampl_trans(ip)-longavg(3))

               phasespace(6,ip)=phasespace(6,ip)

     &              - corr_E*(ampl_trans(ip)-longavg(3))

            end if

         end do

!     

!     create revised 6D correlation matrix after amplitude corr removed

         do i=1,6

            phaseavg(i)=0.

            do j=1,6

               phasecorr(i,j)=0.

            end do

         end do

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               do i=1,6

                  phaseavg(i)=phaseavg(i)+weight*phasespace(i,ip)

                  do j=1,6

                     phasecorr(i,j)=phasecorr(i,j)+

     &                    weight*phasespace(i,ip)*phasespace(j,ip)

                  end do

               end do

            end if

         end do

         if (wtot>0.) then

            do i=1,6

               phaseavg(i)=phaseavg(i)/wtot

            end do

            do i=1,6

               do j=1,6

                  phasecorr(i,j)=phasecorr(i,j)/wtot -

     &                 phaseavg(i)*phaseavg(j)

               end do

            end do

         end if

!

!     calculate longitudinal emittance, beta, alpha

         e_long = phasecorr(5,5)*phasecorr(6,6)-phasecorr(5,6)**2

         if (e_long .gt. 0.) then

            e_long = DSQRT(e_long)

c        normalize to c/(mc^2)

c        refmass is actually mc^2 in GeV

            e_long = e_long * cc / refmass

            delta_long = phasecorr(5,5) / (e_long * refmass / cc)

            alpha_long = phasecorr(5,6) / (e_long * refmass / cc)

         else 

            if (e_long .lt. 0.) 

     &           write(*,*)' error in e_long= ',e_long,

     &           ' region #',curr_region,' npart= ',npart

            delta_long = phasecorr(5,5) 

            alpha_long = phasecorr(5,6) 

         end if

c  define longitudinal beta roughly like transverse, <Pz> sigma_z/sigma_pz

         beta_long = delta_long * cc * avgpz**3 / 

     &        (refmass**2 + avgpz**2)

!     calculate single particle amplitudes

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               if ((e_long .le. 0.) .or. (delta_long .le. 0.)) then

                  ampl_long(ip) = 0.

               else

                  temp_t = phasespace(5,ip) - phaseavg(5)

                  temp_E = phasespace(6,ip) - phaseavg(6)

                  ampl_long(ip) = (cc/refmass) * 

     &                 ( temp_t**2/delta_long + delta_long*

     &                 (temp_E-alpha_long*temp_t/delta_long)**2 )

               end if

            end if

         end do

c

c     calculate invariant emittances for modified correlation matrix

c

         call EMIT3D(refmass,phasecorr,invar_emit)

!

!     calculate number within amplitude cuts, and apply sigma cuts

         w_aperture = 0.

         nscut = 0      ! count # particles removed by sigma cuts this round

         do ip=1,npart

            if (usepar(ip) .eq. 0) then

               weight = evtwtpar(ip)

               if ((transcut .le. 0.) .or. 

     &              (ampl_trans(ip) .le. transcut)) then

                  if ((longcut .le. 0.) .or. 

     &                 (ampl_long(ip) .le. longcut)) then

                     w_aperture = w_aperture + weight

                  end if

               end if

               if (sigma_cut .gt. 0.) then

                  cut_perp = 2*sigma_cut**2 * e_perp

!                 factor of two to make up for combining 4D at once

                  cut_long = sigma_cut**2 * e_long

                  if ((cut_perp .gt. 0.) .and. 

     &                 (ampl_trans(ip) .gt. cut_perp)) then

                     usepar(ip) = knt

                     nscut = nscut + 1

                  else if ((cut_long .gt. 0.) .and.

     &                 (ampl_long(ip) .gt. cut_long)) then

                     usepar(ip) = knt

                     nscut = nscut + 1

                  end if

               end if

            end if

         end do

!

!  check whether should repeat loop

         knt = knt + 1

         if ((sigma_cut .le. 0.) .or. (knt .gt. nloopmax)) then

            cut_flag=.false.

         else

!     Test if tail cutting process has converged, stop if nscut=0

            if (nscut .eq. 0) cut_flag = .false.

         end if

!     

      end do                    ! loop to test cut_flag

!

!  OUTPUTS

      write(2,150)curr_region,zloc,bzaxis,avgpz,Lcanon,

     &     transavg(1),transavg(3),phaseavg(5)+reftime,save_wtot,

     &     w_aperture,invar_emit(1),invar_emit(2),invar_emit(3),

     &     e_perp,e_long,beta_perp,alpha_perp,Ldim,beta_long,alpha_long,

     &     disp_x,disp_y,store_corr_E,store_corr_t,

     &     DSQRT(phasecorr(6,6)),check_axisym

150   format(i6,25e12.4)

!

c

      return

      end

c  *********************************************************

      subroutine LUDCMP(a,n,np,indx,d)

c

      IMPLICIT NONE

c

! 

! Make LU decomposition of matrix

!

      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

c  *********************************************************

      subroutine AXISYM(refmass,avgpz,transcorr,e_perp,

     &     beta_perp,alpha_perp,Ldim,check_axisym)

c

      IMPLICIT NONE

c

      real*8 refmass,avgpz,transcorr(4,4)

      real*8 e_perp,beta_perp,alpha_perp

      real*8 Ldim,check_axisym

c

      integer i,j,indperp(4)

      real*8 matrix(4,4),denom,asymmetry



      do i=1,4

         do j=1,4

            matrix(i,j) = transcorr(i,j)

         end do

      end do

!

! get transverse emittance (square root of 4D emittance)

      call LUDCMP(matrix,4,4,indperp,e_perp)

!

      if (e_perp .ne. 0.) then  ! else singular

         do i=1,4

            e_perp = e_perp * matrix(i,i)

         end do

      end if

!

      if (e_perp .gt. 0.) then

         e_perp = DSQRT(e_perp)

c     take another square root to get a pseudo 2D emittance

         e_perp = DSQRT(e_perp)

c        normalize to 1/(mc)

         e_perp = e_perp / refmass

         beta_perp = (transcorr(1,1) + transcorr(3,3)) * avgpz

     &        / (2. * refmass * e_perp)

         alpha_perp = - (transcorr(1,2) + transcorr(3,4))

     &        / (2. * refmass * e_perp)

         Ldim = (transcorr(1,4) - transcorr(3,2)) 

     &        / ( 2. * refmass * e_perp)

      else 

         if (e_perp .lt. 0.) 

     &        write(*,*)' error in e_perp= ',e_perp

         beta_perp = (transcorr(1,1) + transcorr(3,3)) * avgpz / 2.

         alpha_perp = - (transcorr(1,2) + transcorr(3,4)) / 2.

         Ldim = (transcorr(1,4) - transcorr(3,2)) / 2.

      end if

!

      check_axisym = 0.

      denom = transcorr(1,1) + transcorr(3,3)

      denom = denom**2

      if (denom .gt. 0.) then

         asymmetry = transcorr(1,2)**2 + 

     &        (transcorr(1,1)-transcorr(3,3))**2

         check_axisym = check_axisym + asymmetry/denom

      end if

      denom = transcorr(2,2) + transcorr(4,4)

      denom = denom**2

      if (denom .gt. 0.) then

         asymmetry = transcorr(2,4)**2 + 

     &        (transcorr(2,2)-transcorr(4,4))**2

         check_axisym = check_axisym + asymmetry/denom

      end if

      denom = transcorr(1,1)*transcorr(2,2) + 

     &     transcorr(3,3)*transcorr(4,4) + transcorr(1,1)*transcorr(4,4) 

     &     + transcorr(3,3)*transcorr(2,2)

      if (denom .gt. 0.) then

         asymmetry = (transcorr(1,4)+transcorr(3,2))**2 +

     &        (transcorr(1,2)-transcorr(3,4))**2

         check_axisym = check_axisym + asymmetry/denom

      end if



      return

      end

c  *********************************************************

      subroutine SETUP

c

      IMPLICIT NONE

c

      common/CONTROL/sigma_cut,rffreq,transcut,longcut,

     &  bzaxis,refcharge,refmass,reftime,refenergy,zloc,

     &  pzmax,pzmin,curr_region,vp_model,rf_model,useptype,

     &  ampl_corr,disp_corr

      real*8 sigma_cut,rffreq,transcut,longcut

      real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc

      real*8 pzmax,pzmin

      integer curr_region,vp_model,rf_model,useptype

      logical ampl_corr,disp_corr

c

      integer ioc

      character*80 title

      character*40 file_name_in, file_name_out, file_name_vp

!

!  input file structure:

!    description

!    input file name

!    output file name

!    type of particle (0 = all types)

!    Pzmin, Pzmax (if Pzmax<=0. then no cut)

!    vector potential model, Bz on axis

!    transverse, longitudinal emittance cutoffs

!    subtract out linear correlations with energy and time?

!    subtract out nonlinear correlations with transverse amplitude?

!    time overlap model, RF frequency (in units of MHz)

!    sigma cut (no cut if <=0.)

!    

c   default values

      useptype=0

      pzmin=0.0

      pzmax=0.0

      vp_model=1

      bzaxis=0.0

      transcut=0.0

      longcut=0.0

      disp_corr = .true.

      ampl_corr = .true.

      rf_model = 0

      rffreq=0.0

      sigma_cut=0.0

!

      open(unit=1,file='emitcalc.inp',status='old',iostat=ioc)

	if( ioc .ne. 0 ) then

	  write(*,*) 'Error opening emitcalc.inp'

	  stop

	end if

!

      read(1,'(a80)',iostat=ioc) title

      read(1,*,iostat=ioc) file_name_in

      read(1,*,iostat=ioc) file_name_out

      read(1,*,iostat=ioc) useptype

      read(1,*,iostat=ioc) pzmin,pzmax

      read(1,*,iostat=ioc) vp_model,bzaxis,file_name_vp

      read(1,*,iostat=ioc) transcut,longcut

      read(1,*,iostat=ioc) disp_corr

      read(1,*,iostat=ioc) ampl_corr

      read(1,*,iostat=ioc) rf_model,rffreq

      read(1,*,iostat=ioc) sigma_cut

!     

      close(1)

!

      if (vp_model .ge. 3) then

         open(unit=4,file=file_name_vp,status='unknown',iostat=ioc)

         if( ioc .ne. 0 ) then

            write(*,*) 'Error opening vector potential file ',

     &           file_name_vp

            write(*,*) 'switching to vp model 2, numerical fit'

            vp_model = 2

         else

            call INPUT_VP

         end if

         close(4)

      end if

      if (rffreq .eq. 0.) rf_model = 0 

!

      open(unit=3,file=file_name_in,status='old',iostat=ioc)

      if( ioc .ne. 0 ) then

         write(*,*) 'Error opening data file ',file_name_in

         stop

      end if

      read(3,'(a80)') title     ! title card on beam output file

      read(3,*)

      read(3,*)

c     

      open(unit=2,file=file_name_out,status='unknown',iostat=ioc)

      write(2,*)title

      write(*,*)title

      write(2,*)' output from Emitcalc version 1.0'

      write(*,*)' Emitcalc version 1.0'

      write(2,*)' settings: '

      write(*,*)' settings: '

      write(2,*)'   input file: ',file_name_in

      write(*,*)'   input file: ',file_name_in

      write(2,*)'   output file: ',file_name_out

      write(*,*)'   output file: ',file_name_out

      write(2,*)'   particle type: ',useptype

      write(*,*)'   particle type: ',useptype

      write(2,*)'   pzmin/pzmax: ',pzmin,'/',pzmax

      write(*,*)'   pzmin/pzmax: ',pzmin,'/',pzmax

      if (vp_model .eq. 0) then

         write(2,*)'   vp model: ',vp_model,'   Bz on axis: ',bzaxis

         write(*,*)'   vp model: ',vp_model,'   Bz on axis: ',bzaxis

      else if (vp_model .le. 2) then

         write(2,*)'   vp model: ',vp_model

         write(*,*)'   vp model: ',vp_model

      else 

         write(2,*)'   vp model: ',vp_model,

     &        '   input file: ',file_name_vp

         write(*,*)'   vp model: ',vp_model,

     &        '   input file: ',file_name_vp

      end if

      write(2,*)'   trans/long cuts: ',transcut,longcut

      write(*,*)'   trans/long cuts: ',transcut,longcut

      if (disp_corr) then

         write(2,*)'   remove linear dispersions'

         write(*,*)'   remove linear dispersions'

      else

         write(2,*)'   do not remove linear dispersions'

         write(*,*)'   do not remove linear dispersions'

      end if

      if (ampl_corr) then

         write(2,*)'   subtract out amplitude correlation'

         write(*,*)'   subtract out amplitude correlation'

      else

         write(2,*)'   do not subtract out amplitude correlation'

         write(*,*)'   do not subtract out amplitude correlation'

      end if

      if (rf_model .eq. 0) then

         write(2,*)'   not periodic in time' 

         write(*,*)'   not periodic in time' 

      else 

         write(2,*)'   RF frequency (MHz): ',rffreq,

     &        '   RF model: ',rf_model

         write(*,*)'   RF frequency (MHz): ',rffreq,

     &        '   RF model: ',rf_model

      end if

      write(2,*)'   sigma cut: ',sigma_cut

      write(*,*)'   sigma cut: ',sigma_cut

      write(2,*)

      write(2,30)

 30   format(t1,'regn #',t10,'Z',t22,'Bz0',

     &     t34,'<Pz>',t46,'<Lcanon>',t58,'<x>',t70,'<y>',

     &     t82,'<t>',t94,'total #',t106,'# in cuts',

     &     t118,'emit_1',t130,'emit_2',t142,'emit_3',

     &     t154,'e_perp',t166,'e_long',

     &     t178,'beta_perp',t190,'alpha_perp',t202,'Ldim',

     &     t214,'beta_long',t226,'alpha_long',

     &     t238,'D_x',t250,'D_y',t262,'Corr_E',t274,'Corr_t',

     &     t286,'sigma_E',t298,'asymmetry')

c

      return

      end

c  *********************************************************

      real*8 function TIMESHIFT(time,freq)

c

c     returns time shifted by integer number of periods so that 

c     new time is between +/- 1/(2*freq)

      IMPLICIT NONE

c

      real*8 time, freq

c

      integer count

c

      if (freq .ne. 0.) then

         count = 0.5 + time * freq

         if (count .lt. 0.) count = count - 1.

         TIMESHIFT = time - INT(count) / freq

      else

         TIMESHIFT = time

      end if

c

      return

      end

c  *********************************************************

      subroutine INPUT_VP

c

      IMPLICIT NONE

c

      common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6)

      real*8 vp_axi,vp_x,vp_y

c

      integer ioc

c

      read(4,*,iostat=ioc)vp_x(1,1),vp_x(1,2)

      read(4,*,iostat=ioc)vp_x(2,1),vp_x(2,2),vp_x(2,3)

      read(4,*,iostat=ioc)vp_x(3,1),vp_x(3,2),vp_x(3,3),vp_x(3,4)

      read(4,*,iostat=ioc)vp_x(4,1),vp_x(4,2),vp_x(4,3),vp_x(4,4),

     &     vp_x(4,5)

      read(4,*,iostat=ioc)vp_x(5,1),vp_x(5,2),vp_x(5,3),vp_x(5,4),

     &     vp_x(5,5),vp_x(5,6)

      read(4,*,iostat=ioc)

      read(4,*,iostat=ioc)vp_y(1,1),vp_y(1,2)

      read(4,*,iostat=ioc)vp_y(2,1),vp_y(2,2),vp_y(2,3)

      read(4,*,iostat=ioc)vp_y(3,1),vp_y(3,2),vp_y(3,3),vp_y(3,4)

      read(4,*,iostat=ioc)vp_y(4,1),vp_y(4,2),vp_y(4,3),vp_y(4,4),

     &     vp_y(4,5)

      read(4,*,iostat=ioc)vp_y(5,1),vp_y(5,2),vp_y(5,3),vp_y(5,4),

     &     vp_y(5,5),vp_y(5,6)

c

      return

      end

c  *********************************************************

      subroutine CALC_VP(bzaxis)

c

c     numerical fit to Bz = B0 + B1 r^2 + B2 r^4

c

      IMPLICIT NONE

c

      integer npmax

      parameter ( npmax = 100000 )

c

      common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax),

     &  epar(npmax),bzpar(npmax),evtwtpar(npmax),

     &  iptypar(npmax),usepar(npmax),npart

      real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar

      integer iptypar,usepar,npart

c

      common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6)

      real*8 vp_axi,vp_x,vp_y

c

      real*8 bzaxis

c

      integer i,j,ip

      real*8 corrmatrix(3,3),avg(3),bfit(3),vector3(3),rsq,weight,wtot

      real*8 denom

c

      wtot = 0.

      do i=1,3

         avg(i)=0.

         do j=1,3

            corrmatrix(i,j)=0.

         end do

      end do

      do ip=1,npart

         if (usepar(ip) .eq. 0) then

            weight = evtwtpar(ip)

            wtot = wtot + weight

            rsq = xpar(1,ip)**2 + xpar(2,ip)**2

            vector3(1) = bzpar(ip)

            vector3(2) = rsq

            vector3(3) = rsq*rsq

            do i=1,3

               avg(i)=avg(i)+weight*vector3(i)

               do j=1,3

                  corrmatrix(i,j)=corrmatrix(i,j)+

     &                 weight*vector3(i)*vector3(j)

               end do

            end do

         end if

      end do

      if (wtot>0.) then

         do i=1,3

            avg(i)=avg(i)/wtot

         end do

         do i=1,3

            do j=1,3

               corrmatrix(i,j)=corrmatrix(i,j)/wtot -

     &              avg(i)*avg(j)

            end do

         end do

      end if

c

      denom = corrmatrix(2,2)*corrmatrix(3,3) - corrmatrix(2,3)**2

      if (corrmatrix(2,2) .le. 0.) then

c     no variation

         bfit(3) = 0.

         bfit(2) = 0.

         bfit(1) = avg(1)

      else if (denom .le. 0.) then

c     use linear fit

         bfit(3) = 0.

         bfit(2) = corrmatrix(1,2) / corrmatrix(2,2)

         bfit(1) = avg(1) - bfit(2) * avg(2)

      else

         bfit(3) = (corrmatrix(1,3)*corrmatrix(2,2) - 

     &        corrmatrix(2,3)*corrmatrix(1,2)) / denom

         bfit(2) = (corrmatrix(1,2) - corrmatrix(2,3) * bfit(3)) / 

     &        corrmatrix(2,2)

         bfit(1) = avg(1) - avg(2) * bfit(2) - avg(3) * bfit(3)

      end if

c

      vp_axi(1) = bfit(1) / 2.

      vp_axi(2) = bfit(2) / 4.

      vp_axi(3) = bfit(3) / 6.

c

      bzaxis = bfit(1)

c

      return

      end

c  *********************************************************

      subroutine CALC_TC(freq,centertime,conf)

c

c     for systems periodic in time, find central time

c

      IMPLICIT NONE

c

      integer npmax

      parameter ( npmax = 100000 )

c

      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5)

      real*8 cc,ee,pi,twopi,avagn,relcl,mass

c

      common/REGION/xpar(3,npmax),ppar(3,npmax),tpar(npmax),

     &  epar(npmax),bzpar(npmax),evtwtpar(npmax),

     &  iptypar(npmax),usepar(npmax),npart

      real*8 xpar,ppar,tpar,epar,bzpar,evtwtpar

      integer iptypar,usepar,npart

c

      real*8 freq,centertime,conf

c

      real*8 tcos(npmax),tsin(npmax),xbar,ybar,wtot,weight,avgtime

      integer ip

c

      wtot = 0.

      avgtime = 0.

      xbar = 0.

      ybar = 0.

      do ip=1,npart

         if (usepar(ip) .eq. 0) then

            weight = evtwtpar(ip)

            wtot = wtot + weight

            if (freq .eq. 0.) then

               avgtime = avgtime + weight*tpar(ip)

            else

               xbar = xbar + weight*DCOS(2.*pi*tpar(ip)*freq)

               ybar = ybar + weight*DSIN(2.*pi*tpar(ip)*freq)

            end if

         end if

      end do

      if (wtot>0.) then

            avgtime = avgtime/wtot

            xbar = xbar/wtot

            ybar = ybar/wtot

      end if

c

      if (freq .eq. 0.) then

         if (wtot .eq. 0.) then

            conf = 0.

         else

            conf = 1.

         end if

      else

         conf = DSQRT(xbar**2 + ybar**2)

      end if

      if (conf .eq. 0.) then

         centertime = 0.

      else

         centertime = DATAN2(ybar,xbar) / (2.*pi*freq)

      end if

c

      return

      end

c  *********************************************************

      subroutine EMITXYZ(esq,f,g,h)

c

c     associate emittances with x, y, z planes, and reorder them

c     actually uses esq = (m c emit)^2

c     f = sum of 2x2 determinants across a row

c     g = sum of codeterminants across a row

c     h = projected emittances (just in case other method fails)

c

      IMPLICIT NONE

c

      real*8 SMALL

      parameter (SMALL = 1.d-8)

c

      real*8 esq(3),f(3),g(3),h(3)

c

      integer i,isingle,ntuple,iflip

      real*8 detxy,detxz,detyz,ex,ey,ez,best

      real*8 invx,invy,esingle,epair

      real*8 factor1,factor2,factor3,factor4

      real*8 savex,savey,savez,temp

c

      ntuple=1

      if (dabs(esq(1)-esq(2)).le.small*(esq(1)+esq(2))) then

         ntuple=ntuple+1

         isingle = 3

      end if

      if (dabs(esq(1)-esq(3)).le.small*(esq(1)+esq(3))) then

         ntuple=ntuple+1

         isingle = 2

      end if

      if (dabs(esq(2)-esq(3)).le.small*(esq(2)+esq(3))) then

         ntuple=ntuple+1

         isingle = 1

      end if

      if (ntuple.ge.3) then    ! three emittances basically equal

c     rearrange to be in same order as projected emittances

         call REORDER(3,h,esq)

         return

      else if (ntuple.eq.2) then

c     define solitary and double emittance

         esingle=esq(isingle)

         epair=1.d0

         do i=1,3

            if (i.ne.isingle) epair=epair*esq(i)

         end do

         epair=dsqrt(epair)

         invx=(esingle-f(1))/(esingle-epair)

         invy=(esingle-f(2))/(esingle-epair)

         if (invx.ge.0.5d0) then ! either ex=ey or ex=ez

            if (invy.ge.1.d0-0.5d0*invx) then ! ex=ey

               esq(isingle)=esq(3)

               esq(3)=esingle

            else                ! ex=ez

               esq(isingle)=esq(2)

               esq(2)=esingle

            end if

         else                   ! should have ey=ez

            if ((invy.ge.0.5d0).and.(invy.ge.2.d0*(1.d0-invx))) 

     &           then

               esq(isingle)=esq(1)

               esq(1)=esingle

            else                ! nothing seems to work (large coupling?)

               write(*,*)'difficulty ordering eigenvalues'

               call REORDER(3,h,esq)

            end if

         end if

         return

      end if

c     

c     if made it this far, ntuple=1, all three emittances are distinct

c     don't have to worry about zeroes in denominators

c     first use detxz,detyz to choose ez;  independent of order, ex and ey

c     then use detxy to choose ordering between ex and ey

      best=4.d0

      do i=1,3

         ez=esq(i)

         ex=esq(1)

         ey=esq(2)

         if (i.eq.1) ex=esq(3)

         if (i.eq.2) ey=esq(3)

         factor1=g(1)-ex*ey+ez*(f(1)-ez)

         factor2=g(2)-ex*ez+ez*(f(2)-ey)

         if ((factor1.ne.0.).and.((factor1+factor2).ne.0)) then

            detxz = - (ez-ex)*(ez-ey)/factor1

            detyz = factor1/(factor1+factor2)

c     detxz-1. represents 'amount' of coupling in some sense, 

c     magnitude must be <=1 else corresponds to flipping the beam

            if ((dabs(detxz-1.).le.1).and.(dabs(detyz-1.).le.1))

     &           then           ! good values, check if detxy okay

               factor3=g(1)-ex*ey+ex*(f(1)-ez)

               if (factor3.ne.0.) then

                  detxy = (ex-ey)*(ex-ez)/(detxz*factor3)

                  if (dabs(detxy-1.).le.1) then ! valid option

                     temp=dabs(detxy-1.)+dabs(detxz-1.)

     &                    +dabs(detyz-1.)

                     if(temp.lt.best) then ! best so far

                        if (best.le.3.d0) then

                           write(*,*)'multiple orderings found'

                        end if

                        best=temp

                        savez=ez

                        savex=ex

                        savey=ey

                     end if

                  end if

               end if

c     try exchanging ex and ey

               factor4=g(1)-ey*ex+ey*(f(1)-ez)

               if (factor4.ne.0.) then

                  detxy = (ey-ex)*(ey-ez)/(detxz*factor4)

                  if (dabs(detxy-1.).le.1) then ! valid option

                     temp=dabs(detxy-1.)+dabs(detxz-1.)

     &                    +dabs(detyz-1.)

                     if(temp.lt.best) then ! best so far

                        if (best.le.3.d0) then

                           write(*,*)'multiple orderings found'

                        end if

                        best=temp

                        savez=ez

                        savex=ey

                        savey=ex

                     end if

                  end if

               end if

            end if              ! end of checking detxy

         end if                 ! end of test on detxz and detyz

      end do                    ! loop over choosing ez plane

c

      if (best.gt.3.d0) then    ! no good solution found

         write(*,*)'difficulty ordering eigenvalues'

         call REORDER(3,h,esq)

      else ! did find a valid choice

         esq(1)=savex

         esq(2)=savey

         esq(3)=savez

      end if

c

      return

      end

c  *********************************************************

      subroutine REORDER(ndim,x,y)

c

c     reorder vector y, not sorted but in same order as x

c

      IMPLICIT NONE

c

      integer ndim

      real*8 x(ndim),y(ndim)

c

      integer NMAX

      parameter (NMAX=500)

      integer ishift(NMAX),i,j,itemp

      real*8 temp,new(NMAX)

c

      do i=1,ndim

         ishift(i)=i

      end do

      do i=1,ndim-1

         do j=i+1,ndim

            if (x(i) .gt. x(j)) then

               temp = x(i)

               x(i)=x(j)

               x(j)=temp

               itemp=ishift(i)

               ishift(i)=ishift(j)

               ishift(j)=itemp

            end if

         end do

      end do

      do i=1,ndim-1

         do j=i+1,ndim

            if (y(i) .gt. y(j)) then

               temp = y(i)

               y(i)=y(j)

               y(j)=temp

            end if

         end do

      end do

      do i=1,ndim

         new(i) = y(ishift(i))

      end do

      do i=1,ndim

         y(i) = new(i)

      end do

c

      return

      end

c  *********************************************************

      subroutine GET_VP(bz,x,y,vpx,vpy)

c

c     get transverse components of vector potential for particle

c

      IMPLICIT NONE

c

      common/CONTROL/sigma_cut,rffreq,transcut,longcut,

     &  bzaxis,refcharge,refmass,reftime,refenergy,zloc,

     &  pzmax,pzmin,curr_region,vp_model,rf_model,useptype,

     &  ampl_corr,disp_corr

      real*8 sigma_cut,rffreq,transcut,longcut

      real*8 bzaxis,refcharge,refmass,reftime,refenergy,zloc

      real*8 pzmax,pzmin

      integer curr_region,vp_model,rf_model,useptype

      logical ampl_corr,disp_corr

c

      common/VP/vp_axi(3),vp_x(5,6),vp_y(5,6)

      real*8 vp_axi,vp_x,vp_y

c

      integer i,j

      real*8 bz,x,y,vpx,vpy

      real*8 Aphi,rsq

c

      rsq = x**2 + y**2

      if (vp_model .eq. 0) then

         Aphi = bzaxis / 2.

         vpx = - Aphi * y 

         vpy = Aphi * x

      else if (vp_model .eq. 1) then

         Aphi = (bzaxis + bz)/4.

         vpx = - Aphi * y 

         vpy = Aphi * x

      else if (vp_model .eq. 2) then

         Aphi = vp_axi(1) + vp_axi(2)*rsq + vp_axi(3)*(rsq**2)

         vpx = - Aphi * y 

         vpy = Aphi * x

      else if (vp_model .ge. 3) then

         vpx = 0.

         vpy = 0.

         do i=1,5

            do j=1,i+1

               vpx = vpx + vp_x(i,j)*(x**(i+1-j))*(y**(j-1))

               vpy = vpy + vp_y(i,j)*(x**(i+1-j))*(y**(j-1))

            end do

         end do

      end if

c

      return

      end

c  *********************************************************

      subroutine EMIT3D(refmass,matrix,emittances)

c

c     calculate the 3 separate invariant emittances

c

      IMPLICIT NONE

c

      common/CONST/cc,ee,pi,twopi,avagn,relcl,mass(5)

      real*8 cc,ee,pi,twopi,avagn,relcl,mass

c

      real*8 refmass,matrix(6,6),emittances(3)

      real*8 projected(3),sum,product,det,newmatrix(6,6)

      real*8 factor1,factor2,drow(3),crow(3)

      integer i,j,indx(6)

      real*8 SMALLDET,CODET

c

c     need to rescale t to (ct), so dimensions work out properly

c        including in cross-terms

      do i=1,6

         if (i .eq. 5) then

            factor1 = cc

         else

            factor1 = 1.

         end if

         do j=1,6

            if (j .eq. 5) then

               factor2 = cc

            else

               factor2 = 1.

            end if

            newmatrix(i,j) = matrix(i,j) * factor1 * factor2

         end do

      end do

      do i=1,3

         projected(i) = SMALLDET(i,i,6,newmatrix)

      end do

      sum = 0.

      product = 0.

      do i=1,3

         drow(i)=0.

         crow(i)=0.

         do j=1,3

            drow(i) = drow(i) + SMALLDET(i,j,6,newmatrix)

            crow(i) = crow(i) + CODET(i,j,6,newmatrix)

         end do

         sum = sum + drow(i)

         product = product + crow(i)

      end do

      call LUDCMP(newmatrix,6,6,indx,det)

      if (det .ne. 0.) then  ! else singular

         do i=1,6

            det = det * newmatrix(i,i)

         end do

      end if

c

c     solve x^3 - sum x^2 + product x - det = 0, 

c     solutions are x=(m c emittance)^2

c

      call SOLVE_CUBIC(-sum,product,-det,emittances)

c

c     emittances need to be associated with x, y, z planes:  

c        this is a rough method, which tries to minimize the 

c        necessary decoupling matrix

c

      call EMITXYZ(emittances,drow,crow,projected)

c

c     take square root, and divide by (mc) to get emittance in meters

c

      do i=1,3

         if (emittances(i) .gt. 0.) then

            emittances(i) = DSQRT(emittances(i))/refmass

         end if

      end do

c

      return

      end

c  *********************************************************

      real*8 function SMALLDET(i,j,ndim,matrix)

c

c     2x2 determinant of submatrix using conjugate coordinate pair

c     ndim = # dimensions should be even

c

      IMPLICIT NONE

c

      integer i,j,ndim

      real*8 matrix(ndim,ndim)

c

      SMALLDET = matrix(2*i-1,2*j-1)*matrix(2*i,2*j) -

     &     matrix(2*i,2*j-1)*matrix(2*i-1,2*j)

c

      return

      end

c  *********************************************************

      real*8 function CODET(i,j,ndim,matrix)

c

c     codeterminant of submatrix, from conjugate coordinate pair

c

      IMPLICIT NONE

c

      integer i,j,ndim

      real*8 matrix(ndim,ndim)

      integer NMAX

      parameter (NMAX = 500)

      integer m,n,tm,tn,indx(NMAX)

      real*8 reduced(NMAX,NMAX),determ

c

      do m=1,ndim

         if ((m .ne. 2*i-1) .and. (m .ne. 2*i)) then

            if (m .lt. 2*i-1) then

               tm = m

            else

               tm = m - 2

            end if

            do n=1,ndim

               if ((n .ne. 2*j-1) .and. (n .ne. 2*j)) then

                  if (n .lt. 2*j-1) then

                     tn = n

                  else

                     tn = n - 2

                  end if

                  reduced(tm,tn)=matrix(m,n)

               end if

            end do

         end if

      end do

c

      call LUDCMP(reduced,ndim-2,NMAX,indx,determ)

!

      if (determ .ne. 0.) then  ! else singular

         do m=1,ndim-2

            determ = determ * reduced(m,m)

         end do

      end if

      CODET = determ

c

      return

      end

c  *********************************************************

      subroutine SOLVE_CUBIC(a,b,c,solutions)

c

c     get nonnegative solutions to cubic y=x^3 + a x^2 + b x + c

c     there must be a better way of doing this, hard to get accurate though

c

      IMPLICIT NONE

c

      real*8 SMALL

      parameter (SMALL = 1.d-12)

c

      integer i,j

      real*8 a,b,c,solutions(3)

      real*8 x0,y0,y1,sigma,quadr

      real*8 u,us,f,df,err,close,xc,yc,d,ds,xcs

c

      x0 = -a/3.d0

      y0 = x0**3 + a * x0**2 + b * x0 + c

      quadr = a**2-3.d0*b

      if (quadr.lt.0.d0) then

         if (dabs(quadr).gt.SMALL) then

            write(*,*)'complex roots '

         else

            quadr = 0.d0

         end if

      else

         quadr = DSQRT(quadr)

      end if

      y1 = quadr/3.d0

      sigma = SIGN(1.,y0)

      if (c.eq.0.d0) then

         if ((a.gt.0.d0).or.(b.lt.0.d0)) then

            write(*,*)'negative roots'

         end if

         solutions(1)=0.d0

         quadr=a**2-4.d0*b

         if (quadr.lt.0.d0) then

            if (dabs(quadr).gt.SMALL*a**2) then

               write(*,*)'complex roots '

            else

               quadr = 0.d0

            end if

         else

            quadr = DSQRT(quadr)

         end if

         solutions(2) = 0.5d0*(-a + quadr)

         if (solutions(2).eq.0.d0) then

            solutions(3) = 0.5d0*(-a - quadr)

         else

            solutions(3)=2.d0*b/solutions(2)

         end if

      else

         if ((a.ge.0.d0).or.(b.le.0.d0)) then

            write(*,*)'negative roots '

         end if

         if (y0.le.0.d0) then

            u=2.d0*y1+x0

         else

            u=-2.d0*y1+x0

         end if

         do i=1,20

            f = u**3 + a * u**2 + b * u + c

            df = 3.d0*u**2 + 2.d0* a * u + b

            us = u

            u = u - f/df

         end do

         if (u.eq.0.d0) then

            err = (u-us)/a

         else

            err = (u-us)/u

         end if

         solutions(1)=u

         xc=x0+sigma*y1

         do i=1,20

            xcs = xc

            xc = xc - (3.d0*xc**2+2.d0*a*xc+b)/(6.d0*xc+2.d0*a)

         end do

         if (xc.eq.0.d0) then

            err = (xc-xcs)/y1

         else

            err = (xc-xcs)/xc

         end if

         yc = xc**3 + a * xc**2 + b * xc + c

         close=-sigma*yc/(3.d0*y1)

         if (close.lt.0.d0) then

            if (dabs(close).gt.SMALL*y1**2) then

               write(*,*)'complex roots '

            else

               close = 0.d0

            end if

         end if

         if (close.lt.SMALL*y1**2) then

            d=dsqrt(close)

            do i=1,20

               ds = d

               d= dsqrt(close-sigma*d**3/(3.d0*y1))

            end do

            if (d.eq.0.d0) then

               err = (d-ds)/y1

            else

               err = (d-ds)/d

            end if

            solutions(2)=xc+d

            d=-dsqrt(close)

            do i=1,20

               ds = d

               d= -dsqrt(close-sigma*d**3/(3.d0*y1))

            end do

            if (d.eq.0.d0) then

               err = (d-ds)/y1

            else

               err = (d-ds)/d

            end if

            solutions(3)=xc+d

         else

            if (u.eq.0.d0) then

               quadr=a**2-4.d0*b

            else

               quadr=(a+u)**2+4.d0*c/u

            end if

            if (quadr.lt.0.d0) then

               if (dabs(quadr).gt.SMALL) then

                  write(*,*)'complex roots '

               else

                  quadr = 0.d0

               end if

            else

               quadr = DSQRT(quadr)

            end if

            solutions(2) = 0.5d0*(-a-u+quadr)

            if (u.eq.0.d0) then

               solutions(3) = 0.5d0*(-a-u-quadr)

            else

               solutions(3) = -c/(u*solutions(2))

            end if

            if (quadr.lt.SMALL) then

               if (u.eq.0.d0) then

                  solutions(2)=dsqrt(b)

                  solutions(3)=solutions(2)

               else

                  solutions(2)=x0+sigma*y1

                  solutions(3)=solutions(2)

               end if

            else

               do j=2,3

                  u=solutions(j)

                  do i=1,20

                     f = u**3 + a * u**2 + b * u + c

                     df = 3.d0*u**2 + 2.d0* a * u + b

                     us = u

                     u = u - f/df

                  end do

                  if (u.eq.0.d0) then

                     err = (u-us)/a

                  else

                     err = (u-us)/u

                  end if

                  solutions(j)=u

               end do

            end if

         end if                 ! end of test on near-double root

      end if                    ! end of test on c=0

c

      return

      end

c  *********************************************************

