!+ !subroutine compute_spins(muons,turn,nmuons,unit,ele_name,ix) ! !compute the (average) spin for alive muons and the components of polorization, Qx,Qy,Qz. !Record the information in the specified unit. ! !Modules Needed: ! ! use bmad ! use muon_mod ! use spin_mod ! !Input: ! muons -- muon_struct:the muon array ! turn -- real(rp):number of turns the muons have gone through ! nmuons --integer: number of muons in the array ! unit --integer: unit to record the information ! ele_name --character*(*):name of the element the muons are in ! !Output: ! None !- subroutine compute_spins(muons,turn,nmuons,unit,ele_name,ix,ele_s) USE bmad USE muon_mod USE spin_mod IMPLICIT NONE integer unit integer n,i,j,nmuons,ix integer nmu real(rp) turn,ele_s logical first38/.true./,first39/.true./ real(rp),dimension(3):: Q,tvec !polarization !complex(rp) spin_ave(2) type (muon_struct), dimension(nmuons) :: muons !defined this way to avoid some exceptions 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 if(nmu <= 0)then return endif if (nmu /= 0) then coord%spin = coord%spin/nmu end if ! obtain polarization do j=1,3 Q(j) = real(dot_product(coord%spin,matmul(coord%spin,transpose(pauli(j)%sigma)))) end do if (first38 .and. unit == 38) then write (unit, '(5a10,a19,a19)') 'Element','turn','Qx','Qy','Qz','Remaining','Element Name' first38 = .false. end if if (first39 .and. unit == 39) then write (unit, '(5a10,a19,2a19)') 'Element','turn','Qx','Qy','Qz','Remaining','Element Name','s' first39 = .false. end if write(unit,'(i10,4es12.4,i10,a7,a19,es12.4)')ix,turn,Q(1),Q(2),Q(3),nmu,' ',adjustl(ele_name),ele_s 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 real(rp),dimension(3):: Q !polarization integer n !which muon integer j type (coord_struct) :: coord,ecoord ! for muon and electron logical first/.true./ if (first) write(unit,'(2a10,a14,a11,a6,5a13,3a11,3a13,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','Pol_x','Pol_y','Pol_z' first = .false. ! obtain polarization do j=1,3 Q(j) = real(dot_product(coord%spin,matmul(coord%spin,transpose(pauli(j)%sigma)))) end do write(unit,'(i10,9es12.4,6es12.4,1x,3es12.4)') n,turn,coord%t,coord%path_len,ecoord%vec,coord%vec, Q 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 end do exit end if end do return end subroutine electron_info_ele subroutine electron_info_s(electrons_s,n_row,amount,unit,n,tracker1,unit2) 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,n,unit2 logical sample/.true./,tracker1 real(rp) sum type (electron_struct),dimension(0:n_row,0:amount):: electrons_s if (sample .and. tracker1) then !do j=1,n_row j=1 !if (electrons_s(j,0)%index_mu == n) then do k = 0,amount write(unit2,'(3es12.4,3i10)') electrons_s(j,k)%coord%t,electrons_s(j,k)%coord%vec(1),electrons_s(j,k)%coord%s,electrons_s(j,0)%ix_start,electrons_s(j,k)%coord%state,n end do !exit !end if !end do sample = .false. end if if (.not. tracker1) then !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= 150,n_row !if (electrons_s(j,0)%exist) then !write(101,'(3es12.4,2i10)') electrons_s(j,1)%coord%s,electrons_s(j,2)%coord%s,electrons_s(j,3)%coord%s,electrons_s(j,0)%index_mu,j j = 1 do k = 1,amount write(unit,'(3es12.4,2i10)') electrons_s(j,k)%coord%t,electrons_s(j,k)%coord%vec(1),electrons_s(j,k)%coord%s,electrons_s(j,0)%ix_start,electrons_s(j,k)%coord%state !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 end if !not tracker1 return end subroutine electron_info_s !+ ! Subroutine Compute_moments(muons, turn, moment, nmuons, nmu, unit,ele_name,ix) ! Compute (,,,,,) for all "alive" muons ! Compute moments (6 diagnonal terms >2 ... , ! and 3 cross terms ()(), ()(), ()() ! And write to a file (fort.16 or fort.26) ! Modules Needed: ! use muon_mod ! ! Input: ! muons -- Muon_struct: phase space coordinates of all muons ! turn -- Integer: Turn around ring ! nmuon -- Integer: number of muons in muon structure ! unit -- Integer: unit where output is written ! ele_name -- Character*(*): name of element where moments are calculated ! ix -- Integer: Element number ! ! Output: ! moment -- Moment_struct ! nmu -- Integer: number of "alive" muons that are included in the moments !- SUBROUTINE compute_moments(muons, turn, moment, nmuons, nmu, lun,ele_name,ix) USE muon_mod IMPLICIT NONE INTEGER n,i,nmuons,ix INTEGER nmu integer lun,lun2 REAL(rp) turn,moment_ave_sqr !average of momentum error square logical first/.true./ TYPE (muon_struct), ALLOCATABLE :: muons(:) TYPE (g2moment_struct) moment character*(*) ele_name If (ix == -1) then ! get a time slice of the x and x' distributions (along with other distributions) if (first) then lun2 = lunget() open(unit=lun2,file="180degree_distribution.dat",status='replace') do n=1, nmuons if (muons(n)%coord%state /= alive$)cycle ! .or. muons(n)%lost) cycle write(lun2,'(6es12.4)') muons(n)%coord%vec(1:6) end do close(lun2) first = .false. end if return end if ! Compute average position and sigma for each muon at each turn moment%ave = 0. moment%sigma = 0. moment_ave_sqr = 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) moment_ave_sqr = moment_ave_sqr + (muons(n)%coord%vec(6))**2 nmu = nmu + 1 end do if(nmu <= 0)then print '(a,i10,a)', 'Number of surviving muons = ', nmu,' in COMPUTE_MOMENTS' return endif moment%ave = moment%ave/nmu moment_ave_sqr = moment_ave_sqr/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 (index(ele_name,'BEGINNING') /=0) write(lun,'(a7,a13,12a12,a13,a10,1x,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 sqr_ave','remaining','element name ' write(lun,'(i7,es13.5,12es12.4,es13.4,i10,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_sqr, nmu, ele_name RETURN END SUBROUTINE compute_moments !+ ! Subroutine Collect_time(muons, nturns, NN,bin_max) ! Histogram of muon time coordinate at each turn. Suppose that we measure the arrival times of muons at some detector. ! The arrival times will tend to spread out due to the spread of energies. We want to count the number of muons that ! arrive at the detector in time bins. The histogram is dN/dt ! ! Modules Needed: ! use parameters_bmad ! sue muon_mod ! ! Input: ! muons -- Muon_struct: phase space coordinates of all muons ! nturns -- Integer: Turn number ! ! Output: ! NN -- Real, allocatable :: NN(i) is the temporal distribution of muons at ! bin_max -- Integer: Number of time bins !- 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(-100:nbins)) endif do n=1, size(muons) if(muons(n)%coord%state /= alive$)cycle t = muons(n)%coord%t * 1.e9 t0 = muons(n)%coord%t * 1.e9 bin = max(t/bin_width + 0.5, 0.) bin = t/bin_width + 0.5 if(abs(bin)<=100)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) ! compute maxima and minima of moments of the phase space distribution ! ! Modules needed: ! muon_mod ! ! Input: ! moment -- Moment_struct: moments already computed with subroutine compute_moments ! ! Output: ! moment -- Moment_struct: max_sigma, min_sigma, max_ave, min_ave ! 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) ! compute sigma matrix for the distribution and corresponding twiss parameters ! (Same as subroutine compute_moments with a bit more) ! ! Modules Needed ! muon_mod ! ! Input: ! muons -- Muon_struct: distribution ! verbosity -- Integer, optional (typeout) ! where -- character*(*), optional: string identifying location ! ! Output: ! output, twiss parameters and sigma matrix written to file="BeamParametersInjectionLine.dat" if(where) exists ! and to terminal depending on and whether or not ! 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 logical first/.true./ real(rp) moment_ave_sqr ! average of momentum error square ! This variable controls how much information is printed ! If verbosity is negative, its absolute value is the number of turns verbosity_=1 IF (PRESENT(verbosity)) verbosity_ = verbosity ! Initialize n_alive = 0 moment%ave(:) = 0. moment%sigma(:,:) = 0. moment_ave_sqr = 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) moment_ave_sqr = moment_ave_sqr+ (muons(n)%coord%vec(6))**2 n_alive = n_alive+1 ENDDO IF (n_alive==0) RETURN moment%ave = moment%ave/n_alive moment_ave_sqr= moment_ave_sqr/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) if(present(where))then WRITE(*,'(2x,a,a)') 'CALCULATED BEAM PARAMETERS at:', trim(where) else WRITE(*,'(2x,a)') 'CALCULATED BEAM PARAMETERS at:' endif 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 If (verbosity < 0) then ! now verbosity = - (the number of turns) if (first) then lun = lunget() open(unit=lun,file='params_in_ring.dat', status= 'replace') write(lun,'(a6,4a10,a14)') "turn","emitx","emity",'etax','betax','dp/p sqr_ave' first = .false. close(lun) end if lun = lunget() open(unit=lun,file='params_in_ring.dat', access= 'append') write(lun,'(i4,5es12.4)') -verbosity,emitx,emity,etax,betax,moment_ave_sqr close(lun) end if RETURN END SUBROUTINE compute_beam_params subroutine compute_spin_single_muon(coord,Q) USE bmad USE spin_mod IMPLICIT NONE integer n,i,j,nmuons,ix real(rp),dimension(3):: Q,tvec !polarization !complex(rp) spin_ave(2) type (coord_struct) coord ! obtain polarization do j=1,3 Q(j) = real(dot_product(coord%spin,matmul(coord%spin,transpose(pauli(j)%sigma)))) end do return end subroutine compute_spin_single_muon