! Subroutine create_phase_space(lat, nmuons, create_new_distribution, new_file, muon_file, & ! epsdistr, epsx, epsy, tdistr, tlength, tsigma, pzdistr, pz, pzsigma, twiss2, muons, twiss_ref, mat_inv) ! ! Routine to create distribution of muons. ! ! Input: ! lat: lat_struct ! nmuons -- Integer: number of muons in the distribution ! create_new_distribution -- Logical: If true create a new distribution. If false, read an existing distribution from a file ! new_file -- Character: Name of the file with the new distribution. If empty, that is if new_file = '', then do not write the distribution to a file ! muon_file -- Character: Name of the file of the distribution to read ! epsdistr -- Character: If 'gaus' then new distribution will be gaussian in transverse phase space ! If 'flat' then new distribution will be flat in transverse phase space ! epsx,epsy -- Real: Emittance of gaussian distribution ! tdistr -- Character: If tdistr = 'e821' then time distribution is gaussian with sigma =25ns ! If tdistr = 'e989' then time distribution is FNAL W ! If tdistr = 'gaus' then time distribution is gaussian with sigma = tsigma with tails cut off beyond +- tlength/2 ! tlength -- Real: Maximum time in temporal distribution if the distribution is gaussian (tdistr='gaus') ! tsigma -- Real: Width of gaussian temporal distribution if the distribution is gaussian (tdistr='gaus') ! pzdistr -- Character: If pzdistr = 'gaus' then momentum (vec(6)) distribution is gaussian, with width pzsigma and tails cut off beyond +- pz ! If pzdistr = 'flat' then momentum distribution is flat, with width +- pz/2 ! pz -- Real: Maximum momentum offset of momentum distribution if gaussian (pzdistr='gaus') ! pzsigma -- Real: Width of momentum distribution (vec(6)) if pzdistr = 'gaus' ! twiss2 -- G2_TWISS_STRUCT: Twiss parameters in the inflector. The distribution will be generated so that in the inflector it will have these twiss parameters ! muons -- Muon_struct: Phase space coordinates of the muon distribution generated (or read from a file) ! twiss_ref -- character: 'center' put twiss beam parameters at inflector midpoint, 'end' put twiss beam parameters at end ! mat_inv -- Real,optional :: dimension(6,6): Matrix to propagate phase space defined by TWISS2 backwards to the start of the injection line, (branch = 0) ! nmuon_first - when reading from a list start at nmuon_first subroutine create_phase_space(lat, nmuons, create_new_distribution, new_file, muon_file, & epsdistr, epsx, epsy, tdistr, tlength, tsigma, pzdistr, pz, pzsigma, twiss2, muons, twiss_ref, mat,mat_inv, nmuon_first) !use bmad USE nr use parameters_bmad use muon_mod use muon_interface, dummy => create_phase_space use random_mod use bmad_interface use materials_mod !use lat_ele_loc_mod IMPLICIT NONE type(lat_struct) lat type (averages_struct) averages type (ele_struct), pointer:: ele type (ele_pointer_struct), allocatable :: eles(:) real(rp) Jx_0, Jy_0, beta_x, alpha_x, beta_y, alpha_y, eta, etap, tlength, tsigma, pz, pzsigma, gamma_x, gamma_y real(rp) G(2,2), G_inv(2,2), vec_rot(2), phi real(rp) xvec(2), yvec(2) real(rp) x2_average/0/, y2_average/0/,px2_average/0/, py2_average/0/,xpx_average/0/, ypy_average/0/ real(rp) epsx, epsy real(rp), save :: epsx0, epsy0 real(rp) pz2_average/0/, xpz_average/0/, pxpz_average/0/, ypz_average/0/, pypz_average/0/ real(rp) deltapz real(rp) sigma_x/0.001/, sigma_xp/0.0/, sigma_y/0.001/, sigma_yp/0.0/ real(rp) T(6,6) real(rp) vec(6), ent_inf_vec(6) real(rp), optional :: mat_inv(6,6), mat(6,6) real(rp) L_drift ! drift length for Diktys's distribution real(rp) s0 real(rp) theta real(rp) dtheta real(rp) svec(1:3), pvec(1:3) integer nmuons, nmuon_first integer i integer unit integer lost_at_inflector integer tot integer lun integer pulse, ix integer n_loc integer turn/0/ integer seed integer j character*16 epsdistr, tdistr, pzdistr, inf_aperture, twiss_ref character*120 new_file, muon_file character*100 string logical scatter_inf_end_us, scatter_inf_end_ds, withinInflectorAperture logical create_new_distribution logical energy_loss logical err logical first/.true./ real(rp) s_target/-1000./ ! distance to target (m) real(rp) betagamma real(rp) sigma(6,6) real(rp) Energy, gamma, mc2 type (muon_struct), allocatable :: muons(:), muons_raw(:) type (g2twiss_struct) twiss1, twiss2, twiss_inf call initializeMaterials() call AssignMaterialPointer(material_ptr,'Al') allocate(muons(nmuons)) !do i=1, nmuons ! call init_coord(muons(i)%coord,ele = lat%ele(0), element_end=upstream_end$, particle=antimuon$) !end do lost_at_inflector = 0 ! initialize spins forall(i=1:nmuons)muons(i)%coord%spin = [0, 0, -1] forall(i=1:nmuons)muons(i)%coord%r=1.e7 !lifetime !initial sum BetaCrossE forall(i=1:nmuons)muons(i)%betacrossE(1:3)=0 forall(i=1:nmuons)muons(i)%betacrossE_time=0 forall(i=1:nmuons)muons(i)%path=0 forall(i=1:nmuons)muons(i)%pathx=0 If(create_new_distribution) then if(first)then epsx0=epsx epsy0=epsy first = .false. endif !================================================ ! CREATE 6D PHASE-SPACE DISTRIBUTION (AT TARGET) !================================================ ! Assume for simplicity that transverse positions and momenta are uncorrelated at target, i.e. alpha_{x,y}=0 sigma_xp = epsx0/sigma_x sigma_yp = epsy0/sigma_y IF (epsdistr=='gaus') THEN ! Gaussian distribution call ran_gauss(muons(:)%gas) muons(:)%coord%vec(1) = muons(:)%gas * sigma_x call ran_gauss(muons(:)%gas) muons(:)%coord%vec(2) = muons(:)%gas * sigma_xp call ran_gauss(muons(:)%gas) muons(:)%coord%vec(3) = muons(:)%gas * sigma_y call ran_gauss(muons(:)%gas) muons(:)%coord%vec(4) = muons(:)%gas * sigma_yp ELSEIF (epsdistr=='flat') THEN ! Uniform distribution (default). The width of the flat distribution ! is chosen so as to recover the user's input emittances numerically call ran_uniform(muons(:)%flat) muons(:)%coord%vec(1) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_x call ran_uniform(muons(:)%flat) muons(:)%coord%vec(2) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_xp call ran_uniform(muons(:)%flat) muons(:)%coord%vec(3) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_y call ran_uniform(muons(:)%flat) muons(:)%coord%vec(4) = ((muons(:)%flat)-0.5) * 2*1.73421323168796*sigma_yp ELSEIF (epsdistr == 'round')THEN dtheta = twopi/size(muons(:)) theta=0. i=0 print *,' epsx0 = ',epsx0,' epsy0= ',epsy0 do while (theta <= twopi .and. i < nmuons) theta = theta + dtheta i=i+1 muons(i)%coord%vec(1) = epsx0 * cos(theta) muons(i)%coord%vec(2) = 0. muons(i)%coord%vec(3) = epsy0 * sin(theta) muons(i)%coord%vec(4) = 0. end do ! print *,'line 144' ! do i=1,size(muons(:)) ! print '(i10,1x,6es12.4)',i,muons(i)%coord%vec ! end do ELSE ! delta-function sigma_x = 0. sigma_y = 0. sigma_xp = 0. sigma_yp = 0. muons(:)%coord%vec(1) = 0. muons(:)%coord%vec(2) = 0. muons(:)%coord%vec(3) = 0. muons(:)%coord%vec(4) = 0. ENDIF ! Generate a momentum at the target. Assume that position and momentum are uncorrelated. IF (pzdistr=="gaus" .or. pzdistr=='Gaus' .or. pzdistr=='GAUS') THEN DO i=1, nmuons ! Generate a random momentum call ran_gauss(muons(i)%gas) muons(i)%coord%vec(6) = muons(i)%gas * pzsigma ! Make sure the random momentum lies within [-pz +pz], i.e. "cut the tails off the Gaussian" DO WHILE ( abs(muons(i)%coord%vec(6))>pz ) call ran_gauss(muons(i)%gas) muons(i)%coord%vec(6) = muons(i)%gas * pzsigma ENDDO ENDDO ELSE ! Uniform distribution (default) call ran_uniform(muons(:)%flat) muons(:)%coord%vec(6) = ((muons(:)%flat)-0.5) * pz ENDIF print *,'s_target/betaMagic/c_light =',s_target/betaMagic/c_light endif !only if instructed to create new distribution ! If so instructed read muons from a file if(.not. create_new_distribution .and. len(trim(muon_file)) > 2 .and. index(muon_file,'none') == 0)then lun = lunget() open(unit = lun, file = trim(muon_file)) print '(a,a)',' Read muons from ',trim(muon_file) if(trim(muon_file) == "VDstop_DS_436_12000.dat" .or. trim(muon_file) == "particles_endm4m5_100.txt" & .or. trim(muon_file) == "particles_M4M5End_400_mod.txt" .or. index(muon_file,'particles')/=0)then print '(a,a)',' Read muons from ',trim(muon_file) if(trim(muon_file) == "VDstop_DS_436_12000.dat") & call read_Vmuons(lun,nmuons, muons, s_target/betaMagic/c_light, "Volodya's distribution", tot) ! if(trim(muon_file) == "particles_endm4m5_100.txt") & !Kim 8Apr2016 print '(a,a)',' Read muons from ',trim(muon_file) if(index(muon_file,'particles')/=0) then !Kim 8Apr2016 print '(a,a)',' Read muons from ',trim(muon_file) call read_Dmuons(lun,nmuons, muons, s_target/betaMagic/c_light, "Diktys's distribution", tot, nmuon_first) !Kim 8Apr2016 endif else call read_phase_space(lun, nmuons, muons,'Distribution off target', tot, nmuon_first) endif close(unit=lun) if(tot < nmuons)then print '(a,a,i10,a)',' The number of muons in file ',trim(muon_file), tot,' is less than the number required ' stop endif endif ! how far is the 'BACKLEG_START' element from the beginning of the line s0=0. call lat_ele_locator('BACKLEG_START',lat,eles,n_loc,err) ele => eles(1)%ele s0=ele%s if(err)then print *,' Cannot find BACKLEG_START in create_phase_space ' stop endif print '(a,es12.4)', 'BACKLEG_START at s = ',s0 if(index(twiss_ref,'use_as_is') == 0)call generate_time_long_dist(tdistr, muons, nmuons, tlength, tsigma) ! If so instructed write the distribution of muons to a file print '(3a,i10)', 'new_file =' ,new_file, ' len(new_file) =', len(trim(new_file)) if(len(trim(new_file)) > 2 .and. index(new_file,'none')==0)then muons(:)%coord%state = alive$ call write_phase_space(nmuons, muons,new_file, tot, turn) endif !DO i=1, nmuons ! betagamma = ( 1+muons(i)%coord%vec(6) )*pMagic/mmu ! muons(i)%coord%t = muons(i)%coord%t - s_target/( 1+muons(i)%coord%vec(6)/cosh(asinh(betagamma))**2 )/betaMagic/c_light ! muons(i)%coord%vec(5) = tanh(asinh(betagamma)) * c_light * muons(i)%coord%t !ENDDO !!!!!!!!!!! ! compute twiss parameters at target muons(:)%coord%state = alive$ if (nmuons==1) then muons(1)%Jx = epsx muons(1)%Jy = epsy muons(1)%coord%vec = 0. muons(1)%coord%s = 0. muons(1)%coord%t = 0. if(index(tdistr,'e-')/=0)read(tdistr,*)muons(1)%coord%t ! muons(1)%coord%spin = [0,0,1] muons(1)%coord%spin = [1,0,0] return endif !if(twiss2%betax <0 .or. twiss2%betay<0)then !!!!!!!!!!!!!!!!!! ! Diktys's distribution is produced at 1589 mm after the last final focus Q025 ! End of Q023 is located 9.8773 upstream from inflector exit according to Brynn's plot 'IBMS_correct_positions.pdf' ! End of Q025 is 4.113 m downstream from Q023 ! Our simulation starting point (INJ line) is 9.8773-4.113m - 4.4m (4.3 + 0.1) = 1.366m downstream from Q025 ! ! So the location of Diktys's distribution is L_D=1.589 downstream from Q025 ! and the location of our start point is L_start = 1.366m downstream from Q025 ! So Diktys's distribution must be drifted by (L_start-L_D)=-0.223m if(trim(muon_file) == "VDstop_DS_436_12000.dat" .or. trim(muon_file) == "particles_endm4m5_100.txt" & .or. trim(muon_file) == "particles_M4M5End_400_mod.txt" .or. index(muon_file,'particles')/=0)then L_drift = -0.223 T(:,:) = 0 forall(i=1:6)T(i,i)=1. T(1,2) = L_drift T(3,4) = L_drift ! Propagate muons to the simulation starting point do i=1,nmuons muons(i)%coord%vec = matmul(T,muons(i)%coord%vec) end do endif !muon_file call compute_emittance_beta(nmuons,muons, twiss1, epsx,epsy,averages) !of the raw distribution (Volodya's muons for example) call compute_beam_params(muons, 1, 'RawDistributionAtStart') twiss1%phix = 0. twiss1%phiy = 0. allocate(muons_raw(0:nmuons)) muons_raw = muons if(index(twiss_ref,'use_as_is')/=0)then WRITE(*,'(a)') 'At start of tracking' else WRITE(*,'(a)') 'AT TARGET or raw distribution:' endif WRITE(*,'(5(a15,es12.5))') 'epsx =', epsx, 'sigma_x =', sqrt(averages%x2_average), 'sigma_px =', sqrt(averages%px2_average), 'sigma_xpx =', averages%xpx_average WRITE(*,'(4(a15,es12.5))') 'epsy =', epsy, 'sigma_y =', sqrt(averages%y2_average), 'sigma_py =', sqrt(averages%py2_average), 'sigma_ypy =', averages%ypy_average WRITE(*,'(3(a15,es12.5))') 'sigma_pz =', averages%deltapz,'xpz_average =',averages%xpz_average,'pxpz_average =',averages%pxpz_average WRITE(*,*) WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'beta_x', 'alpha_x', twiss1%betax, twiss1%alphax WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'alpha_x','gamma_x', twiss1%alphax, twiss1%gammax WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'eta_x','eta_px', twiss1%etax, twiss1%etapx WRITE(*,*) WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'beta_y', 'alpha_y', twiss1%betay, twiss1%alphay WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'alpha_y','gamma_y', twiss1%alphay, twiss1%gammay WRITE(*,*) if(index(twiss_ref,'end')/=0)WRITE(*,'(a)') 'AT INFLECTOR EXIT' if(index(twiss_ref,'center')/=0)WRITE(*,'(a)') 'Midway through inflector' if(index(twiss_ref,'brynn')/=0)WRITE(*,'(a)') '10 cm upstream from IBMS1 (9.7cm upstream from start of injection channel)' if(index(twiss_ref,'use_as_is')/=0)WRITE(*,'(a)') 'Use distribution as is.' WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'beta_x', 'alpha_x', twiss2%betax, twiss2%alphax WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'alpha_x','gamma_x', twiss2%alphax, twiss2%gammax WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'eta_x','eta_px', twiss2%etax, twiss2%etapx WRITE(*,*) WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'beta_y', 'alpha_y', twiss2%betay, twiss2%alphay WRITE(*,'(7x,2a9,7x,es12.5,1x,es12.5)') 'alpha_y','gamma_y', twiss2%alphay, twiss2%gammay call write_phase_space( nmuons, muons,'test1', tot, turn) !Raw distribution inflector exit assuming no scattering or apertures !construct transfer matrix from target (or wherever is the origin of the distribution) to final twiss parameters and dispersion halfway through inflector (or at entrance to iron) if(index(twiss_ref,'use_as_is')==0)then if(twiss2%betax > 0 .and. twiss2%betay > 0)then call make_transfer(twiss1,twiss2,T) !transfer matrix from twiss parameters of the raw distribution to the twiss parameters specified in input file, namely midway through inflector or at start else T = mat !matrix that propagates from twiss1 (at target) to twiss2 endif print '(a)', 'Matrix that propagates from twiss1 to twiss2' print '(4(a,es12.4))', 'twiss1: betax =',twiss1%betax,' alphax =',twiss1%alphax,' betay=',twiss1%betay,' alphay=',twiss1%alphay print '(4(a,es12.4))', 'twiss2: betax =',twiss2%betax,' alphax =',twiss2%alphax,' betay=',twiss2%betay,' alphay=',twiss2%alphay print '(4es12.4)',((T(i,j),i=1,4),j=1,4) !propagate muons from production target to halfway through inflector, (or entrance to iron) , from twiss1 to twiss2 do i=1,nmuons muons(i)%Jx = epsx muons(i)%Jy = epsy vec = muons(i)%coord%vec muons(i)%coord%vec = matmul(T,vec) end do !construct transfer matrix from midpoint of inflector to exit (forwards) !assuming downstream half of inflector is a drift which is a pretty good approximation T(:,:) = 0 forall(i=1:6)T(i,i)=1. if(index(twiss_ref,'end')/=0)then T(1,2) = 0. ! dlr 20141009 inf_length/2 T(3,4) = 0. !dlr 20141009 inf_length/2 print '(a,1x,a)', ' Twiss_ref = ',twiss_ref elseif(index(twiss_ref,'center')/=0)then T(1,2) = inf_length/2 T(3,4) = inf_length/2 print '(a,1x,a)', ' Twiss_ref = ',twiss_ref elseif(index(twiss_ref,'brynn')/=0)then T(1,2) = 0.097 !Brynn's twiss parameters are evaluated 10 cm upstream of IBMS1 which is 3mm from the start of the injection line T(3,4) = 0.097 !and the twiss parameters are read in from the lattice file print '(a,1x,a)', ' Twiss_ref = ',twiss_ref else print '(a)',' You must specify twiss parameters reference point, (twiss_ref) in index file' stop endif call write_phase_space( nmuons, muons,'test2', tot, turn) !Raw distribution inflector exit assuming no scattering or apertures ! Propagate muons from midpoint of inflector to exit, or from ref to start of line if(twiss2%betax >0 .and. twiss2%betay > 0)then do i=1,nmuons muons(i)%coord%vec = matmul(T,muons(i)%coord%vec) end do endif else ! skip to here is twiss_end ='use_as_is' T(:,:) = 0 forall(i=1:6)T(i,i)=1. endif if(trim(twiss_ref)=='end' .or. trim(twiss_ref)=='center')string = 'InfExitNoScatter' if(trim(twiss_ref)=='brynn')string = 'AtStartInjectionChannel' if(trim(twiss_ref)=='use_as_is')string='UseAsIs' call compute_beam_params(muons, 1, string) call write_phase_space( nmuons, muons,string, tot, turn) !Raw distribution inflector exit assuming no scattering or apertures !propagate twiss parameters from midpoint of inflector to exit !This is not used anywhere but note that if twiss2 is at start of injection line, this calculation is not correct if(twiss2%betax>0)call prop_phase_space(T, twiss2, twiss_inf) !IF (epsdistr=='flat' .or. epsdistr=='gaus') THEN if(index(twiss_ref,'center')/=0 .or. index(twiss_ref,'end')/=0)then WRITE(*,*) WRITE(*,*) 'AT INFLECTOR MIDPOINT:' if(index(twiss_ref,'center')/=0)WRITE(*,'(a)') 'AT INFLECTOR MIDPOINT' if(index(twiss_ref,'end')/=0)WRITE(*,'(a)') 'AT INFLECTOR END' WRITE(*,'(4(a15,es12.5))') 'betax =', twiss2%betax, 'alphax =', twiss2%alphax, 'etax =', twiss2%etax, 'etapx =', twiss2%etapx WRITE(*,'(4(a15,es12.5))') 'betay =', twiss2%betay, 'alphay =', twiss2%alphay, 'etay =', twiss2%etay, 'etapy =', twiss2%etapy WRITE(*,*) WRITE(*,*) 'AT INFLECTOR EXIT:' WRITE(*,'(4(a15,es12.5))') 'betax =', twiss_inf%betax, 'alphax =', twiss_inf%alphax, 'etax =', twiss_inf%etax, 'etapx =', twiss_inf%etapx WRITE(*,'(4(a15,es12.5))') 'betay =', twiss_inf%betay, 'alphay =', twiss_inf%alphay, 'etay =', twiss_inf%etay, 'etapy =', twiss_inf%etapy endif !ENDIF ! if mat_inv exits then use it to propagate distribution to start of injection line if(present(mat_inv) .and. trim(twiss_ref) /= 'use_as_is' .and. trim(twiss_ref)/= 'brynn')then do i=1,nmuons muons(i)%coord%vec = matmul(mat_inv,muons(i)%coord%vec) end do endif ! set polarization so a)svec/|svec| = -pvec/|pvec| or b) energy dependence from DR if(set_polarization)then Energy=lat%branch(1)%ele(0)%value(e_tot$) ! p=lat%branch(1)%ele(0)%value(p0c$) mc2=mass_of(lat%branch(1)%param%particle) gamma= Energy/mc2 call manage_polarization(nmuons,muons, gamma) endif call write_phase_space(nmuons,muons,'DistStartInjLine', tot, turn) ! Raw distribution propagated forward to end of inflector then back using mat_inv to start of injection line call compute_beam_params(muons,1,'DistStartInjLine') return END SUBROUTINE create_phase_space subroutine makeGR(beta,alpha,phi,G,Ginv,R) use precision_def implicit none real(rp) beta,alpha, phi, G(2,2), Ginv(2,2), R(2,2) G = 0. Ginv = 0. G(1,1) = 1./sqrt(beta) G(2,1) = -alpha/sqrt(beta) G(2,2) = -sqrt(beta) Ginv(1,1) = sqrt(beta) Ginv(2,1) = -alpha/sqrt(beta) Ginv(2,2) = -1./sqrt(beta) R(1,1) = cos(phi) R(1,2) = sin(phi) R(2,1) = -R(1,2) R(2,2) = R(1,1) return end subroutine make_m(eta1,eta0,etap1,etap0,MM,m) use precision_def implicit none real(rp) eta1, eta0, etap1, etap0 real(rp) eta_b(2,2), MM(2,2), m(2,2), eta_e(2,2) eta_e = 0 eta_e(1,2) = eta1 eta_e(2,2) = etap1 eta_b = 0 eta_b(1,2) = eta0 eta_b(2,2) = etap0 m=0 m = eta_e + matmul(MM,eta_b) return end subroutine make_transfer(twiss0,twiss1,T) use precision_def use muon_mod implicit none real(rp) T(6,6), G(2,2), G_inv(2,2), R(2,2), m(2,2),Mx(2,2),Ny(2,2) type (g2twiss_struct) twiss0, twiss1 integer i T=0 do i=1,6 T(i,i)=1. end do call makeGR(twiss0%betax,twiss0%alphax,twiss1%phix-twiss0%phix,G, G_inv,R) Mx = matmul(R,G) call makeGR(twiss1%betax,twiss1%alphax,twiss1%phix-twiss0%phix,G, G_inv,R) T(1:2,1:2) = matmul(G_inv, Mx) call makeGR(twiss0%betay,twiss0%alphay,twiss1%phiy - twiss0%phiy,G, G_inv,R) Ny = matmul(R,G) call makeGR(twiss1%betay,twiss1%alphay,twiss1%phiy - twiss0%phiy,G, G_inv,R) T(3:4,3:4) = matmul(G_inv, Ny) Mx = T(1:2,1:2) call make_m(twiss1%etax,twiss0%etax,twiss1%etapx,twiss0%etapx,Mx,m) T(1:2,5:6) = m call make_m(twiss1%etay,twiss0%etay,twiss1%etapy,twiss0%etapy,T(3:4,3:4),m) T(3:4,5:6) = m return end subroutine write_phase_space(nmuons, muons, string, tot, turn,BetaCrossE_start_time) use muon_mod !use spin_mod use parameters_bmad implicit none integer nmuons integer i, tot integer j integer lun integer turn character*(*) string character(72) plotting_script(20)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu', & 'hist_off_target.gnu',9*' ','hist_start_injection_line.gnu','hist_after_inflector_scatter',4*' '/ type (muon_struct), allocatable :: muons(:) type (muon_struct), allocatable, save :: muons_first_turn(:) real(rp), allocatable, save:: freq(:) real(rp), optional :: BetaCrossE_start_time real(rp) BetaCrossE_st real(rp) pvec(3),svec(3),pdots_norm, smag, phi real(rp) delta_t logical first/.true./ tot=0 lun = lunget() BetaCrossE_st = 0. if(present(BetaCrossE_start_time))BetaCrossE_st = BetaCrossE_start_time if(first)then allocate(freq(nmuons)) freq(:)=0 first = .false. endif if(turn == 1 .and. index(string,'END')/= 0)then if(.not. allocated(muons_first_turn))allocate(muons_first_turn(nmuons)) do i=1,nmuons muons_first_turn(i)%coord%t =muons(i)%coord%t end do endif open(unit = lun, file = trim(directory)//'/'//trim(string)//'_phase_space.dat') !write(lun,'(a,a10,1x,10a12)')'!', 'muon', 'Jx', 'Jy', 'x', 'xp','y','yp','vec(5)','pz','t' write(lun,'(a,a10,1x,14a16,13x,a16,19x,a16,16x,a16,16x,a16,6a16,a16,16x,3a16, a16)')'!', 'muon', 'Jx', 'Jy', 'x', 'xp','y','yp','vec(5)','pz','t','s','s_x','s_y','s_z','tau', & 'BetaCrossE', 'BCrossETime','BetaDotB','frequency','','','','','','','lifetime','B_field','state','decay','acos(s . p)' if(allocated(muons_first_turn) .and. turn > 1)then do i=1,nmuons if(muons(i)%coord%state /= alive$ .and. .not. muons(i)%decay)cycle ! if(muons(i)%coord%state /= alive$ .and. .not. muons(i)%decay)cycle delta_t = muons(i)%coord%t - muons_first_turn(i)% coord%t if(delta_t /= 0)freq(i) = (muons(i)%path_1) / delta_t end do endif do i = 1,nmuons if(muons(i)%coord%state == alive$ .or. (muons(i)%decay .and. muons(i)%coord%t>BetaCrossE_st))then 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 smag = sqrt(dot_product(svec,svec)) phi=0 if(smag > 0)then pdots_norm = dot_product(pvec,svec)/smag/sqrt(dot_product(pvec,pvec)) phi = acos(pdots_norm) endif write(lun,'(i10,1x, 29es16.8,3es16.8,1x,i10,1x,l, 1x,es16.8)') i, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin, muons(i)%coord%r,& muons(i)%BetaCrossE, muons(i)%BetaCrossE_time, muons(i)%BetaDotB, freq(i), muons(i)%ave_vec, muons(i)%coord%r, muons(i)%B_field,muons(i)%coord%state, muons(i)%decay, phi tot = tot+1 endif end do !do i = 1,nmuons ! if(muons(i)%coord%state == alive$)then ! print '(i10,1x,14es12.4)', i, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & ! muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& ! muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin, muons(i)%coord%r ! tot = tot+1 ! endif !end do close(unit=lun) print '(/,a,1x,i10,1x,a,1x,a)','Phase space vector for ',tot,'particles written to:',trim(string)//'_phase_space.dat' print '(6a)',' Use to plot' end subroutine !---- subroutine write_phase_space2(nmuons, muons, string, tot, turn) !by Kim use muon_mod !use spin_mod use parameters_bmad implicit none integer nmuons integer i, tot integer lun integer turn ! turn number character(100) str_turn character*(*) string !character(72) plotting_script(20)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu', & !'hist_off_target.gnu',9*' ','hist_start_injection_line.gnu','hist_after_inflector_scatter',4*' '/ type (muon_struct), allocatable :: muons(:) tot=0 lun = lunget() print *, 'writing~', turn, 'turn' write(str_turn,fmt=*) turn !open(unit = lun, file = trim(string)//'_phase_space.dat') open(unit = lun, file = trim(directory)//'/'//trim(string)//'_'//trim(adjustl(str_turn))//'_phase_space.dat') !write(lun,'(a,a10,1x,9a12)')'!', 'muon', 'Jx', 'Jy', 'x', 'xp','y','yp','vec(5)','pz','t' write(lun,'(a,a10,1x,13a12)')'!', 'muon', 'Jx', 'Jy', 'x', 'xp','y','yp','vec(5)','pz','t','s','s_x','s_y','s_z' !by Kim do i = 1,nmuons if(muons(i)%coord%state == alive$)then write(lun,'(i10,1x,13es13.5)') i, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin !by Kim tot = tot+1 endif end do close(unit=lun) print '(/,a,1x,i10,1x,a,1x,a)','Phase space vector for ',tot,'particles written to:',trim(string)//'_'//trim(adjustl(str_turn))//'_phase_space.dat' !print '(6a)',' Use to plot' return end !read format generated by this program which also works for Valetov data file subroutine read_phase_space(unit, nmuons, muons, string, tot, nmuon_first) use muon_mod implicit none integer unit, nmuons integer nmuon_first integer i, tot,j, ios integer nread character*(*) string character*290 line character(72) plotting_script(4)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu'/ logical include_spin/.false./ logical include_spin_valetov/.false./ type (muon_struct), allocatable :: muons(:) real(rp)initx, inity, initz tot=1 nread=0 do while(tot <= nmuons) read(unit,'(a)', end=99) line if((index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0) .and. (index(line,'s_x')/=0))include_spin=.true. if((index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0) .and. (index(line,'sx')/=0))include_spin_valetov=.true. if(index(line(1:3),'!')/= 0 .or. index(line(1:3),'#')/= 0)cycle if(line(1:10)== ' ')cycle if(nread < nmuon_first) then nread = nread+1 cycle endif if(nread < nmuon_first)cycle i=tot if(include_spin_valetov)then read(line,*, end=99) j, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin, initx, inity, initz elseif(include_spin)then read(line,'(i10,1x,13es16.8)') j, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& muons(i)%coord%t, muons(i)%coord%s, muons(i)%coord%spin else read(line,'(i10,1x,9es16.8)') i, muons(i)%Jx, muons(i)%Jy, muons(i)%coord%vec(1), & muons(i)%coord%vec(2), muons(i)%coord%vec(3), muons(i)%coord%vec(4), muons(i)%coord%vec(5), muons(i)%coord%vec(6),& muons(i)%coord%t endif tot = tot+1 end do 99 print '(a,i10,a)',' Skip ', nread,' muons' print '(/,a,1x,i10,1x,a,1x,i10)','Phase space vector for ',tot,' particles read from unit ',unit return end subroutine read_phase_space subroutine read_Vmuons(unit, nmuons, muons,toff, string, tot) use muon_mod implicit none integer unit, nmuons, nmuon_first ! read nmuons, start at nmuon_first integer i, tot character*(*) string character*120 line character(72) plotting_script(4)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu'/ type (muon_struct), allocatable :: muons(:) real(rp) z,px,py,pz,t, deltap real(rp) vec(1:9), time real(rp) toff time =0. tot=0 i=0 do while(tot < nmuons) read(unit,'(a)', end=99) line if(index(line(1:3),'!')/= 0)cycle i=i+1 read(line,*)vec(1:9) muons(i)%coord%vec(1:5) = vec(1:5)/1000. muons(i)%coord%t = vec(9) * 1.e-9 px=vec(6) py=vec(7) pz=vec(8) deltap = (sqrt(px**2+py**2+pz**2) -3094.35) muons(i)%coord%vec(6) = deltap/3094.35 time = time + vec(9) tot = tot+1 end do print *,' time/tot = ', time/tot muons(1:tot)%coord%t = (muons(1:tot)%coord%t - time/tot * 1.e-9) ! + toff 99 print '(/,a,1x,i10,1x,a,1x,i10)','Phase space vector for ',tot,' particles read from Volodyas file unit ',unit return end subroutine read_Dmuons(unit, nmuons, muons,toff, string, tot, nmuon_first) ! read distributions from Diktys' M4M5 simulation, 6th Apr 2016 ! Under construction use muon_mod !use spin_mod ! by SCK 23May2016 use parameters_bmad implicit none integer unit, nmuons, lun integer nmuon_first integer i, tot, ii integer j character*(*) string character*300 line character(72) plotting_script(4)/'hist_initial_dist.gnu','hist_after_inflector.gnu','hist_start_tracking.gnu','hist_out.gnu'/ type (muon_struct), allocatable :: muons(:) real(rp) z,px,py,pz,t, deltap,ptot real(rp) vec(1:28), time real(rp) svec(1:3) ! spin polarization in cartesian coordinate real(rp) svec_test(1:3) ! spin polarization in cartesian coordinate, debug real(rp) toff real(rp) s s=0; time =0. tot=0 i=0 j=0 do while(i < nmuons) read(unit,'(a)', end=99) line if(index(line(1:3),'#')/= 0 .or. line(1:10) ==' ')cycle read(line,*)vec(1:28) if(vec(7)>9889)cycle if(j < nmuon_first) then j=j+1 cycle endif i=i+1 !muons(i)%coord%vec(1:5) = vec(1:5)/1000. px=vec(4) py=vec(5) pz=vec(6) ptot = sqrt(px**2+py**2+pz**2) muons(i)%coord%vec(1) = vec(1)/1000. !muons(i)%coord%vec(2) = vec(4)/ptot !1st version debug muons(i)%coord%vec(2) = vec(4)/pz muons(i)%coord%vec(3) = vec(2)/1000. !muons(i)%coord%vec(4) = vec(5)/ptot !1st version debug muons(i)%coord%vec(4) = vec(5)/pz muons(i)%coord%vec(5) = vec(3)/1000. muons(i)%coord%t = vec(7) * 1.e-9 deltap = (sqrt(px**2+py**2+pz**2) -3094.35) muons(i)%coord%vec(6) = deltap/3094.35 time = time + vec(7) svec(1) = vec(21) svec(2) = vec(22) svec(3) = vec(23) muons(i)%coord%spin=svec !debug !print '(a,3es12.4)',' muon spin 0',svec !debug !print '(a,3es12.4)',' muon spin 1',svec_test !debug tot = tot+1 end do print *,' time/tot = ', time/tot muons(1:tot)%coord%t = (muons(1:tot)%coord%t - time/tot * 1.e-9) ! + toff lun= lunget() open(unit = lun, file=trim(directory)//'/'//'Diktys_phase_space.dat') write(lun,'(a,a10,1x,10a12,3a12)')'!', 'muon', 'Jx','Jy','x','xp','y','yp','vec(5)','pz','t','s','sx','sy','sz' do ii = 1, nmuons ! min(nmuons,100) write(lun,'(i10,1x,10es12.4,3es12.4)') ii, muons(ii)%jx,muons(ii)%jy,muons(ii)%coord%vec(1), & muons(ii)%coord%vec(2), muons(ii)%coord%vec(3), muons(ii)%coord%vec(4),muons(ii)%coord%vec(5), muons(ii)%coord%vec(6),& muons(ii)%coord%t,s,muons(ii)%coord%spin end do close(unit=lun) 99 print '(/,a,i10,a)', 'Skip first ', j, ' muons' print '(a,1x,i10,1x,a,1x,i10)','Phase space vector for ',tot,' particles read from Diktys file unit ',unit return end subroutine inflector_end_scatter(nmuons, muons) use muon_mod use parameters_bmad use nr implicit none integer unit, nmuons integer i, tot real(rp) dev1, dev2, f, sum real(rp) theta0 type (muon_struct), allocatable :: muons(:) ! apply a scattering angle: do i=1, nmuons f = (sqrt(muons(i)%coord%vec(2)**2+muons(i)%coord%vec(4)**2) + 1.) sum = (f*thickness_al/radlength_al)*(1+0.038*log(f*thickness_al/radlength_al))**2 + & (f*thickness_cu/radlength_cu)*(1+0.038*log(f*thickness_cu/radlength_cu))**2 + & (f*thickness_nbti*0.5/radlength_nb)*(1+0.038*log(f*thickness_nbti*0.5/radlength_nb))**2 + & (f*thickness_nbti*0.5/radlength_ti)*(1+0.038*log(f*thickness_nbti*0.5/radlength_ti))**2 theta0 = ((13.6/momentumunit_MeVperc)/magicmomentum) * sqrt(sum) !print '(i10,1x,2es12.4)',i,f,theta0 call ran_gauss(dev1) call ran_gauss(dev2) muons(i)%coord%vec(2) = muons(i)%coord%vec(2)+theta0*dev1 muons(i)%coord%vec(4) = muons(i)%coord%vec(4)+theta0*dev2 enddo return end ! approximation: leave tangential momentum unchanged, no energy loss subroutine prop_phase_space(T, twiss1, twiss2) use parameters_bmad use muon_mod use muon_interface implicit none type (g2twiss_struct) twiss1, twiss2 real(rp) T(6,6), N(2,2), A(2,2), temp(2,2), Ttrans(2,2), M_inv(2,2) N(1,1:2) = [(1+twiss1%alphax**2)/twiss1%betax, twiss1%alphax] N(2,1:2) = [twiss1%alphax, twiss1%betax] M_inv(1,1:2) = [T(2,2),-T(1,2)] M_inv(2,1:2) = [-T(2,1),T(1,1)] temp = matmul(N,M_inv) Ttrans(1,1:2) = [M_inv(1,1),M_inv(2,1)] Ttrans(2,1:2) = [M_inv(1,2),M_inv(2,2)] A = matmul(Ttrans,temp) twiss2%betax = A(2,2) twiss2%alphax = A(1,2) twiss2%etax= T(1,1)*twiss1%etax + T(1,2)*twiss1%etapx twiss2%etapx= T(2,1)*twiss1%etax + T(2,2)*twiss1%etapx N(1,1:2) = [(1+twiss1%alphay**2)/twiss1%betay, twiss1%alphay] N(2,1:2) = [twiss1%alphay, twiss1%betay] M_inv(1,1:2) = [T(4,4),-T(3,4)] M_inv(2,1:2) = [-T(4,3),T(3,3)] temp = matmul(N,M_inv) Ttrans(1,1:2) = [M_inv(1,1),M_inv(2,1)] Ttrans(2,1:2) = [M_inv(1,2),M_inv(2,2)] A = matmul(Ttrans,temp) twiss2%betay = A(2,2) twiss2%alphay = A(1,2) twiss2%etax= T(1,1)*twiss1%etax + T(1,2)*twiss1%etapx twiss2%etapx= T(2,1)*twiss1%etax + T(2,2)*twiss1%etapx return end SUBROUTINE E989InflectorAperture(x,y,within,us,ds) USE parameters_bmad IMPLICIT NONE REAL(rp), INTENT(IN) :: x,y logical us,ds LOGICAL, INTENT(INOUT) :: within REAL(rp) :: ax,ay ! aperture dimensions real(rp) :: ax_min, ax_max REAL(rp) dx, xtemp ! Inflector aperture half-widths ay = inflector_height ! same as E821 ax = inflector_width xtemp = x - 0.009 !inner aperture does not move. if(us)then dx = 1.7 *tilt xtemp = xtemp+dx endif IF ((xtemp**2/ax**2 + y**2/ay**2) < 1.) THEN ! IF (xtemp > ax_min .and. xtemp < ax_max .and. abs(y) < ay) THEN ! IF (abs(xtemp) < ax .and. abs(y) < ay) THEN within = .true. ENDIF RETURN END SUBROUTINE E989InflectorAperture SUBROUTINE E821InflectorAperture(x,y,within,us,ds) USE parameters_bmad IMPLICIT NONE REAL(rp), INTENT(IN) :: x,y LOGICAL, INTENT(INOUT) :: within logical us,ds REAL(rp) :: m,b ! slope and intercept of slanted part of E821 inflector aperture real(rp) dx, xtemp m = (0.016-0.028)/(0.009-0.002) ! slope b = 0.028 - m*0.002 ! intercept xtemp = x if(us)then dx = 1.7 *tilt xtemp = xtemp+dx ! print '(a,4es12.4)',' x,xtemp,tilt,dx ',x,xtemp,tilt,dx endif IF (abs(xtemp)>0.009 .or. abs(y)>0.028) THEN within = .false. ELSE IF (xtemp>0.002 .and. abs(y)>(m*xtemp+b)) THEN within = .false. ELSE within = .true. ENDIF !if(.not. within)print *, 'not within' RETURN END SUBROUTINE E821InflectorAperture SUBROUTINE RectangularInflectorAperture(x,y,within) USE parameters_bmad IMPLICIT NONE REAL(rp), INTENT(IN) :: x,y LOGICAL, INTENT(INOUT) :: within within = (abs(x)inflector_scatter use materials_mod IMPLICIT NONE type (muon_struct), allocatable :: muons(:) integer i,nmuons logical us, ds, energy_loss ! Scattering and energy loss at upstream inflector end energy_loss = eloss IF (us) THEN DO i=1,nmuons call inflector_scatter1(muons(i)%coord,us,ds,energy_loss) ENDDO ENDIF ! Scattering and energy loss at the downstream end of the inflector IF (ds) THEN DO i=1,nmuons call inflector_scatter1(muons(i)%coord,us,ds,energy_loss) ENDDO ENDIF return end subroutine !inflector_scatter !********************************************* !*********************************************** !*********************************************** subroutine check_inflector_aperture(nmuons,muons, inf_aperture, lost_at_inflector) use parameters_bmad use muon_mod use muon_interface, dummy => check_inflector_aperture use materials_mod IMPLICIT NONE type (muon_struct), allocatable :: muons(:) integer i,nmuons, lost_at_inflector logical withinInflectorAperture, us,ds character*16 inf_aperture real(rp) vec(6) us = .true. ds = .true. ! Check to see if muons clear the inflector aperture DO i=1, nmuons vec = muons(i)%coord%vec ! upstream end of inflector withinInflectorAperture = .false. IF (inf_aperture=='e989' .or. inf_aperture=='E989') THEN call E989InflectorAperture(vec(1),vec(3),withinInflectorAperture,us,ds) ELSEIF (inf_aperture=='e821' .or. inf_aperture=='E821') THEN call E821InflectorAperture(vec(1),vec(3),withinInflectorAperture, us, ds) ELSEIF (inf_aperture=='rect' .or. inf_aperture=='Rect' .or. inf_aperture=='RECT') THEN call RectangularInflectorAperture(vec(1),vec(3),withinInflectorAperture) ELSE withinInflectorAperture = .true. ENDIF ! Kill the track if necessary IF (muons(i)%coord%state==alive$ .and. .not.withinInflectorAperture) THEN muons(i)%coord%state = lost$ lost_at_inflector = lost_at_inflector + 1 ENDIF ENDDO return end subroutine !check_inflector_aperture subroutine compute_emittance_beta(nmuons,muons, twiss1, epsx,epsy,averages) use bmad use muon_mod implicit none type (muon_struct), allocatable :: muons(:) type (averages_struct) averages type (g2twiss_struct) twiss1 real(rp) x2_average/0/, y2_average/0/,px2_average/0/, py2_average/0/,xpx_average/0/, ypy_average/0/ real(rp) epsx, epsy real(rp) pz2_average/0/, xpz_average/0/, pxpz_average/0/, ypz_average/0/, pypz_average/0/ real(rp) deltapz real(rp) sigma_x/0.001/, sigma_xp/0.0/, sigma_y/0.001/, sigma_yp/0.0/ real(rp) sigma(6,6) integer nmuons, i x2_average = 0. y2_average = 0. px2_average = 0. py2_average = 0. xpx_average = 0. ypy_average = 0. pz2_average = 0. xpz_average = 0. pxpz_average = 0. ypz_average = 0. pypz_average = 0. do i=1,nmuons x2_average = x2_average + muons(i)%coord%vec(1)**2 y2_average = y2_average + muons(i)%coord%vec(3)**2 px2_average = px2_average + muons(i)%coord%vec(2)**2 py2_average = py2_average + muons(i)%coord%vec(4)**2 xpx_average = xpx_average + muons(i)%coord%vec(1)*muons(i)%coord%vec(2) ypy_average = ypy_average + muons(i)%coord%vec(3)*muons(i)%coord%vec(4) pz2_average = pz2_average + muons(i)%coord%vec(6)**2 xpz_average = xpz_average + muons(i)%coord%vec(1) * muons(i)%coord%vec(6) pxpz_average = pxpz_average + muons(i)%coord%vec(2)*muons(i)%coord%vec(6) ypz_average = ypz_average + muons(i)%coord%vec(3)*muons(i)%coord%vec(6) pypz_average = pypz_average + muons(i)%coord%vec(4)*muons(i)%coord%vec(6) end do x2_average = x2_average /nmuons y2_average = y2_average /nmuons px2_average = px2_average/nmuons py2_average = py2_average/nmuons xpx_average = xpx_average/nmuons ypy_average = ypy_average/nmuons pz2_average = pz2_average/nmuons xpz_average = xpz_average/nmuons pxpz_average = pxpz_average/nmuons ypz_average = ypz_average/nmuons pypz_average = pypz_average/nmuons if(abs(pz2_average)>0)then sigma(1,1) = x2_average - xpz_average**2/pz2_average sigma(1,2) = xpx_average - xpz_average*pxpz_average/pz2_average sigma(2,2) = px2_average - pxpz_average**2/pz2_average sigma(3,3) = y2_average - ypz_average**2/pz2_average sigma(3,4) = ypy_average - ypz_average*pypz_average/pz2_average sigma(4,4) = py2_average - pypz_average**2/pz2_average else sigma(1,1) = x2_average sigma(1,2) = xpx_average sigma(2,2) = px2_average sigma(3,3) = y2_average sigma(3,4) = ypy_average sigma(4,4) = py2_average endif epsx = sqrt(x2_average*px2_average - xpx_average**2) epsy = sqrt(y2_average*py2_average - ypy_average**2) epsx = sqrt(sigma(1,1)*sigma(2,2)-sigma(1,2)**2) epsy = sqrt(sigma(3,3)*sigma(4,4)-sigma(3,4)**2) deltapz = sqrt(pz2_average) averages%x2_average = x2_average /nmuons averages%y2_average = y2_average /nmuons averages%px2_average = px2_average/nmuons averages%py2_average = py2_average/nmuons averages%xpx_average = xpx_average/nmuons averages%ypy_average = ypy_average/nmuons averages%pz2_average = pz2_average/nmuons averages%xpz_average = xpz_average/nmuons averages%pxpz_average = pxpz_average/nmuons averages%ypz_average = ypz_average/nmuons averages%pypz_average = pypz_average/nmuons twiss1%betax = averages%x2_average/epsx twiss1%alphax = averages%xpx_average/epsx twiss1%gammax = averages%px2_average/epsx twiss1%betay = averages%y2_average/epsy twiss1%alphay = averages%ypy_average/epsy twiss1%gammay = averages%py2_average/epsy if(abs(averages%pz2_average)>0)then twiss1%etax = (averages%xpz_average)/pz2_average twiss1%etapx = (averages%pxpz_average)/pz2_average endif twiss1%betax = sigma(1,1)/epsx twiss1%alphax = -sigma(1,2)/epsx twiss1%gammax = sigma(2,2)/epsx twiss1%betay = sigma(3,3)/epsy twiss1%alphay = -sigma(3,4)/epsy twiss1%gammay = sigma(4,4)/epsy if(abs(averages%pz2_average)>0)then twiss1%etay = (averages%ypz_average)/averages%pz2_average twiss1%etapy = (averages%pypz_average)/averages%pz2_average endif return end subroutine inflector_scatter1(orb,us,ds,energy_loss) use bmad use materials_mod use parameters_bmad implicit none type (coord_struct) orb logical us, ds, energy_loss REAL(rp) surface_normal(3)/0.,0.,1./ !print '(a,6es12.4)','inflector scatter 1 in:',orb%vec(1:6) if((us .and. us_scatter) .or. (ds .and. ds_scatter))then call scatter( Al, 0.0015_rp, orb,surface_normal) ! flange window = 1.5 mm Al !print '(a,6es12.4)','inflector scatter 1 Al:',orb%vec(1:6) call scatter( InflectorCoilE821, 0.0033_rp, orb,surface_normal) ! coil = 3.3 mm kapton-wrapped, aluminum-stabilized NbTi/Cu superconducting wires !print '(a,6es12.4)','inflector scatter 1 InfCoil:',orb%vec(1:6) call scatter( Al, 0.0015_rp, orb,surface_normal) ! mandrel window = 1.5 mm Al IF (eloss) THEN call energyLoss( Al, 0.0015_rp, orb) ! flange window = 1.5 mm Al call energyLoss( InflectorCoilE821, 0.0033_rp, orb) ! coil = 3.3 mm kapton-wrapped, aluminum-stabilized NbTi/Cu superconducting wires call energyLoss( Al, 0.0015_rp, orb) ! flange window = 1.5 mm Al ENDIF endif !print '(a,6es12.4)','inflector scatter 1 out:',orb%vec(1:6) return end !generate temporal distribution centered at 0, and translate to vec(5) subroutine generate_time_long_dist(tdistr,muons, nmuons, tlength, tsigma) use muon_mod use parameters_bmad use muon_interface, dummy => generate_time_long_dist implicit none type (muon_struct), allocatable :: muons(:) character*16 tdistr integer nmuons integer ix, pulse, i real(rp) tlength, tsigma ! Generate a time and longitudinal position at the target. By default, the target is assumed to be at s=-1000 meters (CDR). call str_upcase(tdistr,tdistr) IF (tdistr=='e989' .or. tdistr=='E989') THEN ! Generate a random time with W distribution DO i=1, nmuons muons(i)%coord%t = muons(i)%coord%t + fnalw(0.0_rp,tlength) ENDDO ELSEIF (index(tdistr, 'E989-P')/=0 )then ix = index(tdistr,'P') print *,' tdistr = ',tdistr, ' ix = ', ix, ' len(tdistr) = ', len(trim(tdistr)) if(len(trim(tdistr))== 7)then read(tdistr(ix+1:ix+1),'(i1)')pulse elseif(len(trim(tdistr))== 8) then read(tdistr(ix+1:ix+2),'(i2)')pulse else print *, 'generate_time_long_dist Cannot translate tdistr to 7 or 8 characters ' stop endif if(pulse < 1 .or. pulse > 10)then !pulse = 1-8 are the 8 pulses. pulse=9 is the average of all 8 Chris's average is 10 print *,' error reading pulse type in create_phase_space' stop endif DO i=1,nmuons call muon_pulse_shape(pulse, muons(i)%coord%t) END DO ELSEIF (tdistr=='e821' .or. tdistr=='E821') THEN ! Generate a random time (Guassian with sigma=25ns) call ran_gauss(muons(:)%gas) muons(:)%coord%t = muons(:)%gas * 25.e-9 ELSEIF (tdistr=='gaus' .or. tdistr=='Gaus' .or. tdistr=='GAUS') THEN DO i=1, nmuons ! Generate a random time call ran_gauss(muons(i)%gas) muons(i)%coord%t = muons(i)%gas * tsigma ! Ensure the random time lies within [-tlength/2, +tlength/2], i.e. "cut the tails off the Gaussian" DO WHILE (abs(muons(i)%coord%t)>tlength/2.) call ran_gauss(muons(i)%gas) muons(i)%coord%t = muons(i)%gas * tsigma ENDDO ENDDO ELSE ! Uniform distribution (default) call ran_uniform(muons(:)%flat) muons(:)%coord%t = ((muons(:)%flat)-0.5) * tlength ENDIF ! vec(5) is distance from center of bunch muons(:)%coord%vec(5) = betaMagic*c_light*muons(:)%coord%t return end subroutine FUNCTION fnalw(center,width) USE parameters_bmad ! fnalw_integral table IMPLICIT NONE REAL(rp) :: fnalw ! quantity to determine REAL(rp), INTENT(IN) :: center, width ! center and width of distribution [a.u.] INTEGER, PARAMETER :: nbins=100 ! number of bins in fnalw_integral table REAL(rp) rand, di, frac_low, frac_high, i_interp ! helper variables INTEGER i ! iterator integer seed/1/ logical first/.true./ ! Generate a random number in [0,1) ! if(first)call ran_seed_put(seed) ! first=.false. call ran_uniform(rand) ! Find the random number in the fnalw_integral table, and convert to a number based on the user's center/width arguments DO i=1, 99 IF (fnalw_integral(i)<=rand .and. rand