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 = 100000 ) parameter ( version = 2.05 ) 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,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 ! ! 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 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, REAL(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, 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,nregion, & t_tp,t_xp,t_pp,Bp,t_evtwt,Ep,phasep,polp 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 per = 1.0d-6 / rffreq if (rffreq .ne. 0.d0) then 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 = 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,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,0,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, 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