!+ !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 use parameters_bmad IMPLICIT NONE integer unit integer n,i,j,nmuons,ix integer nmu integer, save :: lun(0:1) real(rp) turn,ele_s, ave_spin(3) real(rp),dimension(3):: tvec !polarization !complex(rp) spin_ave(2) type (muon_struct), dimension(nmuons) :: muons !defined this way to avoid some exceptions character*(*) ele_name logical first/.true./ if(unit < 0)then close(lun(0)) close(lun(1)) return endif if(first)then lun(0)= lunget() open(unit=lun(0), file = trim(directory)//'/'//'average_spin_tbt.dat') lun(1) = lunget() open(unit=lun(1), file = trim(directory)//'/'//'average_spin_by_element.dat') first=.false. endif ! compute average spin for alive muons nmu = 0 ave_spin = 0 do i =1, nmuons if (muons(i)%coord%state /= alive$) cycle ! .or. muons(n)%lost) cycle !ave_spin = muons(i)%ave_spin ave_spin = ave_spin + muons(i)%coord%spin nmu = nmu + 1 end do if(nmu <= 0) return ave_spin = ave_spin/nmu ! obtain polarization if (index(ele_name,'BEGINNING') /=0) then ! on first call, write header to file write (lun(unit),'(a7,a13,3a12,a10,1x,a16,a12)') & 'Element','turn ','Qx ','Qy ','Qz ','Remaining','Element Name ','s ' else ! otherwise write average spin info to file write (lun(unit),'(i7,es13.5,3es12.4,i10,1x,a16,es12.4)') & ix,turn, ave_spin,nmu,adjustl(ele_name),ele_s endif end subroutine compute_spins !---- subroutine decay_info(turn,n,ecoord,coord, muons) use bmad use parameters_bmad use muon_mod IMPLICIT NONE integer, save :: lun41 real(rp) turn real(rp) t,t0 !the actually obtained decay time vs the assigned lifetime integer n !which muon integer j type (coord_struct) :: coord,ecoord ! for muon and electron type (muon_struct), optional :: muons logical first/.true./ if (first)then lun41 = lunget() open(unit=lun41,file=trim(directory)//'/'//'muon_decay.dat') write(lun41,'(a10,3a12,6a12,6a12,4a12,12x,a12,12x,a16,12x,a12)') & '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 ','frac_turn ','BetaCrossE','BetaCrossE_time','BetaDotB' first = .false. endif ! if(present(muons))then write(lun41,'(i10,3es12.4,6es12.4,6es12.4,3es12.4, es12.4,3es12.4,es16.4,3es12.4)') & n,turn,coord%t,coord%r,ecoord%vec,coord%vec, coord%spin, turn-int(turn),& muons%betaCrossE,muons%betaCrossE_time, muons%betaDotB else write(lun41,'(i10,3es12.4,6es12.4,6es12.4,3es12.4, es12.4)') & n,turn,coord%t,coord%r,ecoord%vec,coord%vec, coord%spin, turn-int(turn) endif end subroutine decay_info !---- subroutine electron_info_ele(electrons_ele,nmuons,n_ele,unit) use bmad use muon_mod use parameters_bmad 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,n,tracker1) use bmad use muon_mod use parameters_bmad 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 n integer, save :: lun, unit2 logical sample/.true./,tracker1 logical first/.true./ real(rp) sum type (electron_struct),dimension(0:n_row,0:amount):: electrons_s if(first)then lun = lunget() open(unit = lun,file = trim(directory)//'/'//'file42.dat') unit2 = lunget() open(unit=unit2, file = trim(directory)//'/'//'file43.dat') first=.false. endif 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(lun,'(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 ('Beam_moments_tbt.dat' or 'Beam_moments_by_element.dat') ! 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, unit,ele_name,ix) USE muon_mod use parameters_bmad use muon_interface, dummy => compute_moments IMPLICIT NONE INTEGER n,i,nmuons,ix INTEGER nmu integer, save :: lun(0:2), lun2, lun3 integer unit REAL(rp) turn,moment_ave_sqr !average of momentum error square logical first/.true./ logical first_time/.true./ TYPE (muon_struct), ALLOCATABLE :: muons(:) TYPE (g2moment_struct) moment character*(*) ele_name logical, allocatable, save :: already_lost(:) if(first_time)then allocate(already_lost(size(muons))) already_lost = .false. lun(0)=lunget() open(unit=lun(0), file = trim(directory)//'/'//'Beam_moments_tbt.dat') lun(1) = lunget() open(unit=lun(1), file = trim(directory)//'/'//'Beam_moments_by_element.dat') lun(2) = lunget() lun3=lunget() open(unit = lun3, file = trim(directory)//'/'//'Lost_muons_by_element.dat') first_time = .false. do n=1,nmuons muons(n)%n=0 muons(n)%ave_vec(:)=0 end do endif 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=trim(directory)//'/'//'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 muons(n)%n = muons(n)%n+1 muons(n)%ave_vec(:) = (muons(n)%ave_vec(:)*(muons(n)%n-1) + muons(n)%coord%vec(:))/muons(n)%n end do ! this is now done by subroutine 'time_bins' which is called from 'tracking_master' ! if(bin_width > 0 .and. pbin_width >0 & ! .and. trim(time_binning) == 'momentum_slice' & ! .and. index(ele_name,'BEGINNING')==0 )call bin_spin_angle_by_time(turn,unit, nmuons,muons) 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(unit),'(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','remaining','element name ' write(lun(unit),'(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(6), nmu, ele_name if(unit == 1)then if( index(ele_name,'BEGINNING')/=0)write(lun3,'(a10,1x,a7,1x,a16,1x,a5,1x,a20,1x,a,1x,a)') & 'muon','element','name','state','state name','vec','time' do n=1,nmuons if(muons(n)%coord%t > 10.e-6 .and. muons(n)%coord%state /= alive$ .and. .not. already_lost(n) )then write(lun3,'(i10,1x,i7,1x,a16,1x,i5,1x,a20,1x,6es12.4,1x,es12.4)') & n, ix,ele_name,muons(n)%coord%state, coord_state_name(muons(n)%coord%state), & muons(n)%coord%vec,muons(n)%coord%t already_lost(n) = .true. endif end do endif 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) local_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/local_bin_width + 0.5, 0.) bin = t/local_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 use parameters_bmad 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 use parameters_bmad 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 if(abs(moment%sigma(6,6))>0)then sigma_beta(i,j) = moment%sigma(i,j) - moment%sigma(i,6)*moment%sigma(j,6)/moment%sigma(6,6) ! Eq. (14.28) else sigma_beta(i,j) = moment%sigma(i,j) endif 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 ! ... if(abs(moment%sigma(6,6))>0)then 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) ! ... endif ! 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=trim(directory)//'/'//'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)')'Number of muons = ', 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=trim(directory)//'/'//'params_in_ring_180.dat', status= 'replace') write(lun,'(a6,4a10,a14)') "turn","emitx","emity",'etax','betax','dp/p sqr_ave' first = .false. close(lun) lun = lunget() open(unit=lun,file=trim(directory)//'/'//'params_in_ring_270.dat', status= 'replace') write(lun,'(a6,4a10,a14)') "turn","emitx","emity",'etax','betax','dp/p sqr_ave' close(lun) end if if (-verbosity >=1000000) then lun = lunget() open(unit=lun,file=trim(directory)//'/'//'params_in_ring_270.dat', access= 'append') write(lun,'(i4,5es12.4)') -verbosity-1000000,emitx,emity,etax,betax,moment_ave_sqr close(lun) else lun = lunget() open(unit=lun,file=trim(directory)//'/'//'params_in_ring_180.dat', access= 'append') write(lun,'(i4,5es12.4)') -verbosity,emitx,emity,etax,betax,moment_ave_sqr close(lun) end if end if RETURN END SUBROUTINE compute_beam_params subroutine bin_spin_angle_by_time(turn, unit, nmuons,muons, dphi_dp_ext, dphi_dp_err_ext) use nr use muon_mod use parameters_bmad implicit none type (muon_struct), allocatable :: muons(:) real(rp), save :: p_min, p_max, p_over real(rp) pvec(3), svec(3),sdotp,sdotp_norm, spin_phase real(rp), save, allocatable :: spin_phase_sum(:,:),spin_phase_prod(:,:) real(rp), allocatable :: avg_phase(:), p(:), rms(:), err_avg(:) real(rp), allocatable :: dphi_dp(:),dphi_dp_err(:) real(rp), optional, allocatable :: dphi_dp_ext(:),dphi_dp_err_ext(:) real(rp) siga, sigb, a,b, chi2,q real(rp) turn real(rp) one/1./ integer, save, allocatable :: npart(:,:) integer unit, nmuons integer i integer, save :: n_time_bins integer m integer, save :: npbins integer time_bin, p_bin logical first/.true./ if(first .and. unit ==1)then n_time_bins = number_turns*149.e-9/bin_width + 100 p_min = -0.005 p_max= 0.005 ! npbins = 50 npbins = (p_max-p_min)/pbin_width + 2 p_over=0. print *,' number of time bins = ',n_time_bins allocate(spin_phase_sum(-npbins/2:npbins/2,0:n_time_bins), spin_phase_prod(-npbins/2:npbins/2,0:n_time_bins)) allocate(npart(-npbins/2:npbins/2,0:n_time_bins)) do i=1, n_time_bins spin_phase_sum(:,i)=0 spin_phase_prod(:,i)=0 npart(:,i)=0 end do first=.false. endif if(turn > 0 .and. unit ==1)then do i=1, nmuons if (muons(i)%coord%state /= alive$)cycle ! .or. muons(n)%lost) cycle time_bin = muons(i)%coord%t/bin_width +0.5 ! associate integer with each bin_width time slice if(time_bin > n_time_bins .or. time_bin <0)cycle pvec(1) = muons(i)%coord%vec(2) pvec(2) = muons(i)%coord%vec(4) pvec(3) = sqrt((1+muons(i)%coord%vec(6))**2-pvec(1)**2-pvec(2)**2) svec = muons(i)%coord%spin if(dot_product(svec,svec) == 0)cycle sdotp = dot_product(svec, pvec) sdotp_norm= sdotp/sqrt(dot_product(svec,svec)*dot_product(pvec,pvec)) if(muons(i)%coord%vec(6) >= p_min .and. muons(i)%coord%vec(6) <= p_max)then p_bin = muons(i)%coord%vec(6)/pbin_width + 0.5 spin_phase = acos(sdotp_norm) * sign(one, (pvec(3)*svec(1)-pvec(1)*svec(3))) spin_phase_sum(p_bin,time_bin) = spin_phase + spin_phase_sum(p_bin,time_bin) spin_phase_prod(p_bin,time_bin) = spin_phase**2 + spin_phase_prod(p_bin,time_bin) npart(p_bin,time_bin) = 1+npart(p_bin,time_bin) else p_over=p_over+1 endif ! end do elseif(turn == -1 .or. unit == 0)then if(.not. allocated(avg_phase))allocate(avg_phase(1:npbins)) if(.not. allocated(rms))allocate( rms(1:npbins), p(1:npbins), err_avg(1:npbins)) if(.not. allocated(dphi_dp))allocate(dphi_dp(0:n_time_bins), dphi_dp_err(0:n_time_bins)) if(.not. allocated(dphi_dp).and. present(dphi_dp_ext))allocate(dphi_dp_ext(0:n_time_bins), dphi_dp_err_ext(0:n_time_bins)) dphi_dp(0:n_time_bins)=0; dphi_dp_err(0:n_time_bins)=0 if(pbin_width <=0)return !nothing to calculate, this routine just sets dphi_dp and dphi_dp_err to zero open(unit = 345,file = trim(directory)//'/'//'spin_vs_momentum_at_time.dat') open(unit = 346,file = trim(directory)//'/'//'domega_ddelta_vs_time.dat') do i=1, n_time_bins write(345,'(/,a1,/)')'#' m=0 do p_bin = -npbins/2, npbins/2 if(npart(p_bin,i) < 3)cycle m=m+1 avg_phase(m) = spin_phase_sum(p_bin,i)/npart(p_bin,i) rms(m) = sqrt(spin_phase_prod(p_bin,i)/npart(p_bin,i) - avg_phase(m)**2) ! err_avg(m) = sqrt(rms(m)**2*(1 + 1./4/npart(p_bin,i))+avg_phase(m)**2/npart(p_bin,i)) err_avg(m) = 1./sqrt(float(npart(p_bin,i))) p(m) = p_bin write(345,'(2es12.4,i10,4es12.4)')i*bin_width,p_bin*pbin_width,m,p(m),avg_phase(m),rms(m), err_avg(m) end do if(m>=3)then call fit(p(1:m),avg_phase(1:m),a,b,siga,sigb,chi2,q,err_avg(1:m)) write(345,'(a1,i10,a12,es12.4,a12,es12.4,a12,es12.4)')'#',i,' slope= ',b,' offset=',a,' err=',sigb dphi_dp(i) = b dphi_dp_err(i) = sigb write(346,'(3es12.4)')i*bin_width,dphi_dp(i), dphi_dp_err(i) else write(345,'(a,i10)')'Number of points m = ',m endif end do if(present(dphi_dp_ext))then dphi_dp_ext = dphi_dp dphi_dp_err_ext = dphi_dp_err endif close(unit=346) close(unit=345) endif return end subroutine bin_spin_angle_by_time