subroutine compute_spins(muons,turn,nmuons,unit,ele_name,ix) USE bmad USE muon_mod USE spin_mod IMPLICIT NONE integer unit integer n,i,j,nmuons,ix integer nmu real(rp) turn logical first38/.true./,first39/.true./ real(rp),dimension(3):: Q !polarization !complex(rp) spin_ave(2) type (muon_struct), dimension(nmuons) :: muons !defined this way to avoid some execptions type (pauli_struct),dimension(3):: sig !pauli matrices type (coord_struct) coord character*(*) ele_name ! compute average spin for alive muons nmu = 0 coord%spin = (/(0.,0.),(0.,0.)/) do i =1, nmuons if (muons(i)%coord%state /= alive$) cycle ! .or. muons(n)%lost) cycle coord%spin = muons(i)%coord%spin !coord%spin = coord%spin + muons(i)%coord%spin !nmu = nmu + 1 end do !coord%spin = coord%spin/nmu ! obtain polarization sig(1)%sigma = reshape((/(0.,0.),(1.,0.),(1.,0.),(0.,0.)/),(/2,2/)) sig(2)%sigma = reshape((/(0.,0.),(0.,1.),(0.,-1.),(0.,0.)/),(/2,2/)) sig(3)%sigma = reshape((/(1.,0.),(0.,0.),(0.,0.),(-1.,0.)/),(/2,2/)) do j=1,3 Q(j) = real(dot_product(coord%spin,matmul(coord%spin,transpose(sig(j)%sigma)))) end do if (first38 .and. unit == 38) then write (unit, '(5a10,2a19)') 'Element','turn','Qx','Qy','Qz','Remaining','Element Name' first38 = .false. end if if (first39 .and. unit == 39) then write (unit, '(5a10,2a19)') 'Element','turn','Qx','Qy','Qz','Remaining','Element Name' first39 = .false. end if write(unit,'(i10,4es12.4,i10,a19)')ix,turn,Q(1),Q(2),Q(3),nmu,adjustl(ele_name) return end subroutine compute_spins subroutine decay_info(unit,turn,n,ecoord,coord) use bmad IMPLICIT NONE integer unit real(rp) turn real(rp) t,t0 !the actually obtained decay time vs the assigned lifetime integer n !which muon type (coord_struct) :: coord,ecoord ! for muon and electron logical first/.true./ if (first) write(unit,'(2a10,a14,a11,a6,5a13,3a11,3a13)') 'n','turn','decay time', 'lifetime','e_x','e_px','e_y','e_py','e_z','e_pz', & 'mu_x','mu_px','mu_y','mu_py','mu_z','mu_pz' first = .false. write(unit,'(i10,9es12.4,6es12.4)') n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec return end subroutine Decay_info subroutine electron_info_ele(electrons_ele,nmuons,n_ele,unit) use bmad use muon_mod IMPLICIT NONE integer nmuons integer n_ele !number of elements integer j,k,count/0./,ix_start integer unit real(rp) sum type (electron_struct),dimension(0:nmuons,0:n_ele):: electrons_ele do j = 1,nmuons sum = 0. if (electrons_ele(j,0)%exist) then do k = 1,n_ele sum = sum + electrons_ele(j,k)%coord%vec(1)**2 end do if (sum == 0) then count = count+1 ix_start = electrons_ele(j,0)%ix_start print *, 'state',electrons_ele(j,ix_start)%coord%state end if end if end do print *,'count',count do j= 1,nmuons if (electrons_ele(j,0)%exist) then do k = 1,n_ele write(unit,'(2es12.4,i10)') electrons_ele(j,k)%coord%t,electrons_ele(j,k)%coord%vec(1),electrons_ele(j,0)%ix_start print *,'test',electrons_ele(j,k)%coord%t,electrons_ele(j,k)%coord%vec(1),electrons_ele(j,0)%ix_start,'n',j end do exit end if end do return end subroutine electron_info_ele subroutine electron_info_s(electrons_s,n_row,amount,unit) use bmad use muon_mod IMPLICIT NONE integer nmuons integer amount !number of columns in electrons_s integer n_row !number of rows in electrons_s integer j,k,count/0./,ix_start integer unit real(rp) sum type (electron_struct),dimension(0:n_row,0:amount):: electrons_s do j = 1,n_row !sum = 0 if (electrons_s(j,0)%exist) then !do k = 1,amount !sum = sum + electrons_s(j,k)%coord%vec(1)**2 !end do !if (sum == 0.) then count = count+1 !ix_start = electrons_s(j,0)%ix_start !print *, 'state',electrons_s(j,ix_start)%coord%state !end if end if end do print *,'count',count do j= 240,n_row if (electrons_s(j,0)%exist) then do k = 1,amount write(unit,'(2es12.4,i10)') electrons_s(j,k)%coord%t,electrons_s(j,k)%coord%vec(1),electrons_s(j,0)%ix_start print *,'test',electrons_s(j,k)%coord%t,electrons_s(j,k)%coord%vec(1),electrons_s(j,0)%ix_start,'n',j end do exit end if end do return end subroutine electron_info_s SUBROUTINE compute_moments(muons, turn, moment, nmuons, nmu, unit,ele_name,ix) USE muon_mod IMPLICIT NONE INTEGER unit INTEGER n,i,nmuons,ix INTEGER nmu REAL(rp) turn LOGICAL first16/.true./, first26/.true./, last TYPE (muon_struct), ALLOCATABLE :: muons(:) TYPE (g2moment_struct) moment character*(*) ele_name ! Compute average position and sigma for each muon at each turn moment%ave = 0. moment%sigma = 0. ! Compute (,,,,,) for alive muons nmu=0 do n=1, nmuons if (muons(n)%coord%state /= alive$)cycle ! .or. muons(n)%lost) cycle moment%ave(1) = moment%ave(1) + muons(n)%coord%vec(1) moment%ave(2) = moment%ave(2) + muons(n)%coord%vec(2) moment%ave(3) = moment%ave(3) + muons(n)%coord%vec(3) moment%ave(4) = moment%ave(4) + muons(n)%coord%vec(4) moment%ave(5) = moment%ave(5) + muons(n)%coord%t moment%ave(6) = moment%ave(6) + muons(n)%coord%vec(6) nmu = nmu + 1 end do moment%ave = moment%ave/nmu ! Comment do n=1, nmuons if (muons(n)%coord%state /= alive$)cycle ! .or. muons(n)%lost) cycle ! Compute terms on the diagonal moment%sigma(1,1) = moment%sigma(1,1) + ( muons(n)%coord%vec(1)-moment%ave(1) )**2 moment%sigma(2,2) = moment%sigma(2,2) + ( muons(n)%coord%vec(2)-moment%ave(2) )**2 moment%sigma(3,3) = moment%sigma(3,3) + ( muons(n)%coord%vec(3)-moment%ave(3) )**2 moment%sigma(4,4) = moment%sigma(4,4) + ( muons(n)%coord%vec(4)-moment%ave(4) )**2 moment%sigma(5,5) = moment%sigma(5,5) + ( muons(n)%coord%t -moment%ave(5) )**2 moment%sigma(6,6) = moment%sigma(6,6) + ( muons(n)%coord%vec(6)-moment%ave(6) )**2 ! Compute off-diagonal terms of 2x2 sub-blocks moment%sigma(1,2) = moment%sigma(1,2) + ( muons(n)%coord%vec(1)-moment%ave(1) )*( muons(n)%coord%vec(2)-moment%ave(2) ) moment%sigma(3,4) = moment%sigma(3,4) + ( muons(n)%coord%vec(3)-moment%ave(3) )*( muons(n)%coord%vec(4)-moment%ave(4) ) moment%sigma(5,6) = moment%sigma(5,6) + ( muons(n)%coord%t -moment%ave(5) )*( muons(n)%coord%vec(6)-moment%ave(6) ) end do moment%sigma(1:6,1:6) = moment%sigma(1:6,1:6)/nmu if (first26 .or. first16) write(unit,'(2a10,1x,13a12,1x,a8,a16)')'Element','Turn','x_ave','y_ave','t_ave','sigma(1,1)', & 'sigma(1,2)','sigma(2,2)','sigma(3,3)','sigma(3,4)','sigma(4,4)','sigma(5,5)','sigma(5,6)','sigma(6,6)','dp/p ave', 'remaining', & 'element name' if(unit == 16)first16=.false. if(unit == 26)first26 = .false. write(unit,'(i10,es12.4,1x,13es12.4,1x,i7,1x,a16)')ix,turn,moment%ave(1),moment%ave(3), moment%ave(5), & moment%sigma(1,1), moment%sigma(1,2), moment%sigma(2,2), & moment%sigma(3,3), moment%sigma(3,4), moment%sigma(4,4), & moment%sigma(5,5), moment%sigma(5,6), & moment%sigma(6,6), moment%ave(6), nmu, adjustl(ele_name) RETURN END SUBROUTINE compute_moments SUBROUTINE collect_times(muons, nturns, NN, bin_max) USE parameters_bmad USE muon_mod IMPLICIT NONE REAL(rp) bin_width/5./ ! ns REAL(rp) t, t0 INTEGER, SAVE :: nbins INTEGER bin, n REAL(rp), ALLOCATABLE :: NN(:) INTEGER nturns INTEGER bin_max TYPE (muon_struct), ALLOCATABLE :: muons(:) if (.not. allocated(NN)) then nbins = 150 * (nturns+2) allocate(NN(0:nbins)) endif do n=1, size(muons) t = muons(n)%coord%t * 1.e9 t0 = muons(n)%coord%t * 1.e9 bin = max(t/bin_width + 0.5, 0.) NN(bin) = NN(bin) + 1*(1.+0.1*cos((omega_a/1.e9)*(t-t0)))/2.*exp(-(t-t0)/(lifetime*1.e9)) if (bin>nbins) then print *,' nturns, bin, nbin, size(NN) = ',nturns, bin, nbins, size(NN) print *,' t = ', t endif bin_max = max(bin, bin_max) end do RETURN END SUBROUTINE collect_times SUBROUTINE init_moment(moment) USE muon_mod IMPLICIT NONE TYPE (g2moment_struct) moment moment%max_sigma = -100. moment%min_sigma = 100. moment%max_ave = -100. moment%min_ave = 100. RETURN END SUBROUTINE init_moment SUBROUTINE moment_range(moment) USE muon_mod IMPLICIT NONE TYPE (g2moment_struct) moment moment%max_sigma = max(moment%max_sigma, moment%sigma) moment%min_sigma = min(moment%min_sigma, moment%sigma) moment%max_ave = max(moment%max_ave, moment%ave) moment%min_ave = min(moment%min_ave, moment%ave) RETURN END SUBROUTINE moment_range SUBROUTINE compute_beam_params(muons,verbosity,where) USE muon_mod IMPLICIT NONE ! Subroutine arguments TYPE (muon_struct), ALLOCATABLE, INTENT(IN) :: muons(:) INTEGER, OPTIONAL :: verbosity character*(*), optional :: where ! Subroutine internal variables TYPE(g2moment_struct) :: moment ! see "muon_mod.f90" REAL(RP) :: sigma_beta(6,6) ! betatron-only part of covariance matrix REAL(RP) :: emitx, betax, alphax, gammax, etax, etapx ! parameters to calculate REAL(RP) :: emity, betay, alphay, gammay, etay, etapy ! parameters to calculate INTEGER :: i, j, n, n_alive ! counters & loop indices INTEGER :: verbosity_ ! silly fortran integer lun ! This variable controls how much information is printed verbosity_=1 IF (PRESENT(verbosity)) verbosity_ = verbosity ! Initialize n_alive = 0 moment%ave(:) = 0. moment%sigma(:,:) = 0. ! Compute {,,,,,} for alive muons DO n=1, size(muons) IF (muons(n)%coord%state/=alive$) CYCLE moment%ave(1) = moment%ave(1) + muons(n)%coord%vec(1) moment%ave(2) = moment%ave(2) + muons(n)%coord%vec(2) moment%ave(3) = moment%ave(3) + muons(n)%coord%vec(3) moment%ave(4) = moment%ave(4) + muons(n)%coord%vec(4) moment%ave(5) = moment%ave(5) + muons(n)%coord%vec(5) moment%ave(6) = moment%ave(6) + muons(n)%coord%vec(6) n_alive = n_alive+1 ENDDO IF (n_alive==0) RETURN moment%ave = moment%ave/n_alive ! Compute covariance matrix DO n=1, size(muons) IF (muons(n)%coord%state/=alive$) CYCLE DO i=1, 6 DO j=i, 6 moment%sigma(i,j) = moment%sigma(i,j) + ( muons(n)%coord%vec(i)-moment%ave(i) )*( muons(n)%coord%vec(j)-moment%ave(j) ) ENDDO ENDDO ENDDO DO i=1, 6 DO j=i, 6 moment%sigma(i,j) = moment%sigma(i,j)/n_alive moment%sigma(j,i) = moment%sigma(i,j) ENDDO ENDDO ! Compute betatron-only part of covariance matrix (see e.g. BMAD Manual v22.0, Section 14.2) DO i=1, 6 DO j=1, 6 sigma_beta(i,j) = moment%sigma(i,j) - moment%sigma(i,6)*moment%sigma(j,6)/moment%sigma(6,6) ! Eq. (14.28) ENDDO ENDDO ! Compute beam parameters emitx = sqrt( sigma_beta(1,1)*sigma_beta(2,2)-sigma_beta(1,2)**2 ) ! RMS emittance (@Gaussian) emity = sqrt( sigma_beta(3,3)*sigma_beta(4,4)-sigma_beta(3,4)**2 ) ! ... betax = sigma_beta(1,1)/emitx ! Eqs. (14.32) betay = sigma_beta(3,3)/emity ! ... alphax = -sigma_beta(1,2)/emitx ! ... alphay = -sigma_beta(3,4)/emity ! ... gammax = sigma_beta(2,2)/emitx ! ... gammay = sigma_beta(4,4)/emity ! ... etax = moment%sigma(1,6)/moment%sigma(6,6) ! Eqs. (14.25) etay = moment%sigma(3,6)/moment%sigma(6,6) ! ... etapx = moment%sigma(2,6)/moment%sigma(6,6) ! ... etapy = moment%sigma(4,6)/moment%sigma(6,6) ! ... ! FIXME: Aren't these formulas more accurate? ! etax = (moment%sigma(1,6)-(-alphax/gammax)*moment%sigma(2,6))/moment%sigma(6,6) ! etay = (moment%sigma(3,6)-(-alphay/gammay)*moment%sigma(4,6))/moment%sigma(6,6) IF (verbosity_>1) THEN ! Print covariance matrix WRITE(*,*) WRITE(*,'(a)') REPEAT('=',150) WRITE(*,'(2x,a)') 'COVARIANCE MATRIX (BETATRON + DISPERSION PARTS):' WRITE(*,'(a)') REPEAT('-',150) DO i=1,6 WRITE(*,'(6es25.15)') moment%sigma(i,1:6) ENDDO WRITE(*,'(a)') REPEAT('=',150) ! Print betatron-only part of covariance matrix WRITE(*,*) WRITE(*,'(a)') REPEAT('=',150) WRITE(*,'(2x,a)') 'COVARIANCE MATRIX (BETATRON PART ONLY):' WRITE(*,'(a)') REPEAT('-',150) DO i=1,6 WRITE(*,'(6es25.15)') sigma_beta(i,1:6) ENDDO WRITE(*,'(a)') REPEAT('=',150) ENDIF ! Print calculated beam parameters if(present(where))then lun = lunget() open(unit=lun,file='BeamParamsInjectionLine.dat', access= 'append') print '(a26,1x,a)','Beam Parameters written to', 'BeamParamsInjectionLine.dat' WRITE(lun,*) WRITE(lun,'(a,i5)') 'Number of muons = ',n_alive WRITE(lun,'(a)') REPEAT('=',150) WRITE(lun,'(2x,2a)') 'CALCULATED BEAM PARAMETERS at:', trim(where) WRITE(lun,'(a)') REPEAT('=',150) WRITE(lun,'(a4,7a20)') '', 'emit(m-rad)', '4*emit(m-rad)', 'beta(m)', 'alpha()', 'gamma(1/m)', 'eta(m)', 'etap()' WRITE(lun,'(a)') REPEAT('-',150) WRITE(lun,'(a4,7es20.10)') ' x:', emitx, 4.*emitx, betax, alphax, gammax, etax, etapx WRITE(lun,'(a4,7es20.10)') ' y:', emity, 4.*emity, betay, alphay, gammay, etay, etapy WRITE(lun,'(a)') REPEAT('=',150) WRITE(lun,*) close(unit=lun) endif IF (verbosity_>0) THEN WRITE(*,*) WRITE(*,'(a,i5)') n_alive WRITE(*,'(a)') REPEAT('=',150) WRITE(*,'(2x,a)') 'CALCULATED BEAM PARAMETERS:' WRITE(*,'(a)') REPEAT('=',150) WRITE(*,'(a4,7a20)') '', 'emit(m-rad)', '4*emit(m-rad)', 'beta(m)', 'alpha()', 'gamma(1/m)', 'eta(m)', 'etap()' WRITE(*,'(a)') REPEAT('-',150) WRITE(*,'(a4,7es20.10)') ' x:', emitx, 4.*emitx, betax, alphax, gammax, etax, etapx WRITE(*,'(a4,7es20.10)') ' y:', emity, 4.*emity, betay, alphay, gammay, etay, etapy WRITE(*,'(a)') REPEAT('=',150) WRITE(*,*) ENDIF RETURN END SUBROUTINE compute_beam_params