subroutine time_bins(co, n_ele_track, nmuons, time_binning_start, ele) use bmad use nr use bmad_struct use bmad_interface use parameters_bmad use muon_interface, dummy=>time_bins implicit none type (coord_struct), allocatable :: co(:) type (ele_struct),allocatable :: ele(:) 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), save, allocatable :: spin_phase_save(:) real(rp), allocatable :: avg_phase(:), p(:), rms(:), err_avg(:) real(rp), allocatable :: dphi_dp(:),dphi_dp_err(:) real(rp), allocatable :: dphi_dp_ext(:),dphi_dp_err_ext(:) real(rp) time_min, time_max real(rp) write_time/20.e-6/ real(rp) one/1./ real(rp) bin_center_time real(rp) siga, sigb, a,b, chi2,q real(rp) time_binning_start real(rp) sxp(3) !svec X pvec rotated so that in +- y direction real(rp) spin_phase1 integer i integer, save :: n_time_bins integer m integer, save :: npbins integer p_bin integer bin,j,n_ele_track integer bin_max, bin_min integer, save, allocatable :: npart(:,:) integer nmuons integer, save :: start_time_bin integer index integer ncycles integer, allocatable, save :: write_bins(:) logical first/.true./ integer, allocatable, save :: already_used(:,:) integer n_writes/20/ logical, allocatable,save :: mask(:) logical close_files/.false./ character*4 cindex ! time bins, bin_width wide, starting from zero !call reallocate(co, n_ele_track) if(.not. allocated(write_bins))allocate(write_bins(1:n_writes)) if(n_ele_track > 0)then if(first)then start_time_bin = time_binning_start/bin_width n_time_bins = number_turns*149.e-9/bin_width + 100 - start_time_bin forall(i=1:n_writes) write_bins(i) = float(n_time_bins-100)*i/n_writes p_min = -0.015 p_max= 0.015 ! 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)) allocate(mask(1:n_ele_track)) allocate(already_used(1:nmuons,0:n_time_bins)) allocate(spin_phase_save(1:nmuons)) spin_phase_save(:)=0 do i=1, n_time_bins spin_phase_sum(:,i)=0 spin_phase_prod(:,i)=0 npart(:,i)=0 already_used(:,i)=0 end do open(unit = 301,file = trim(directory)//'/'//'spin_phase_vs_time.dat') write(301,'(a12,a12,a12)')'muon_number', 'time_bin', 'spin_phase' first=.false. endif ! hardwired so that if any element has CALO in name, only those elements with CALO in name will be used if(any(index(ele(:)%name,'CALO')/=0))then mask(:)=.false. do i=1,n_ele_track if(index(ele(i)%name,'CALO')/=0)mask(i)=.true. end do ! print *,'CALO YES' else mask(1:n_ele_track)=.true. ! print *,'NO CALO' endif ! do i=1,n_ele_track ! print '(i10,a,l)',i,' mask =', mask(i) ! end do time_min = co(1)%t time_max= co(n_ele_track)%t ! if(n_ele_track>0)print '(a,es12.4,a,es12.4)','subroutine time_bin: time_min =',time_min, ' time_max = ', time_max bin_min = time_min/bin_width - start_time_bin bin_max = time_max/bin_width + 1 - start_time_bin do bin=bin_min,bin_max bin_center_time = (bin+start_time_bin) *bin_width -bin_width/2. j=minloc(abs(bin_center_time - co(1:n_ele_track)%t),1,mask) if(abs(bin_center_time - co(j)%t) > bin_width/2)cycle if(already_used(muon_number,bin)>0)cycle mask(j)=.false. already_used(muon_number,bin)= already_used(muon_number,bin)+1 pvec(1) = co(j)%vec(2) pvec(2) = co(j)%vec(4) pvec(3) = sqrt((1+co(j)%vec(6))**2-pvec(1)**2-pvec(2)**2) svec = co(j)%spin if(dot_product(svec,svec) == 0)cycle if(co(j)%vec(6) >= p_min .and. co(j)%vec(6) <= p_max)then p_bin = co(j)%vec(6)/pbin_width + 0.5 call RightOrLeft(pvec,svec,sxp) ! finds part of svec in the x-p plane and rotates pvec,svec so that cross product (sxp) is in the +/- y direction ! print '(a,9es12.4)','pvec,svec,sxp =',pvec,svec,sxp sdotp = dot_product(svec, pvec) sdotp_norm= sdotp/sqrt(dot_product(svec,svec)*dot_product(pvec,pvec)) spin_phase = atan2(sxp(2),sdotp) ! spin_phase1 = acos(sdotp_norm) * sign(one, (pvec(3)*svec(1)-pvec(1)*svec(3))) ! print '(a,2es12.4)','spin_phase_method 1 and2', spin_phase, spin_phase1 if(spin_phase < 0)spin_phase = spin_phase + twopi if(spin_phase < spin_phase_save(muon_number))then ncycles = 1 + (spin_phase_save(muon_number)-spin_phase)/twopi spin_phase = spin_phase + ncycles * twopi endif spin_phase_save(muon_number) = spin_phase write(301,'(i12,i12,es12.4)')muon_number, bin, spin_phase if(npart(p_bin,bin)>0)then if((spin_phase - spin_phase_sum(p_bin,bin)/npart(p_bin,bin))>twopi/2.) spin_phase=spin_phase-twopi if((spin_phase - spin_phase_sum(p_bin,bin)/npart(p_bin,bin))<-twopi/2.) spin_phase=spin_phase+twopi endif spin_phase_sum(p_bin,bin) = spin_phase + spin_phase_sum(p_bin,bin) spin_phase_prod(p_bin,bin) = spin_phase**2 + spin_phase_prod(p_bin,bin) npart(p_bin,bin) = 1+npart(p_bin,bin) else p_over=p_over+1 endif if(any(write_bins(:)==bin))then write_time = (bin+start_time_bin)*bin_width - bin_width/2 call write_phase_space_at_time(co(j),spin_phase, write_time, bin, write_bins, .false.) endif end do ! time bins elseif(n_ele_track <0)then close_files = .true. call write_phase_space_at_time(co(1),spin_phase, write_time, bin, write_bins, close_files) 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)) 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') write(345,'(2a12,5a12)')'time','m','p(m)','avg_phase','rms', 'err_avg','a+b*p' open(unit = 346,file = trim(directory)//'/'//'domega_ddelta_vs_time.dat') write(346,'(4a12,a10)')'time','dphi_dp', 'dphi_dp_err','chi2','index' ! open(unit = 347,file = trim(directory)//'/'//'already_used.dat') ! write(347,'(3a12)')'muon','time_bin',' uses' 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*pbin_width end do if(m>=3)then do j=1,m if(j>1)then if(avg_phase(j)-avg_phase(1) > twopi/2.)avg_phase(j)=avg_phase(j)-twopi if(avg_phase(j)-avg_phase(1) < -twopi/2.)avg_phase(j)=avg_phase(j)+twopi if(avg_phase(j)-avg_phase(1) > twopi/2.)avg_phase(j)=avg_phase(j)-twopi if(avg_phase(j)-avg_phase(1) < -twopi/2.)avg_phase(j)=avg_phase(j)+twopi endif end do a=0. b=0. 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 do j=1,m write(345,'(es12.4,i12,5es12.4)')(i+start_time_bin)*bin_width,j,p(j),avg_phase(j),rms(j), err_avg(j), a+b*p(j) end do dphi_dp(i) = b dphi_dp_err(i) = sigb write(346,'(4es12.4,i10)')(i+start_time_bin)*bin_width,dphi_dp(i), dphi_dp_err(i),chi2,i else write(345,'(a,i10)')'Number of points m = ',m endif end do ! do j=1, nmuons ! do i=1,n_time_bins ! if(already_used(j,i)>0)write(347,'(3i12)')j,i, already_used(j,i) ! end do ! enddo ! 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 !muon_number < 0 return end subroutine time_bins ! rotate pvec and svec so pvec_p is along z-axis. ! rotate svec_p so that it is in the x-z plane ! compute cross product cross = svec X pvec ! rotate cross to cross'' - it will be in the +- ydirection subroutine RightOrLeft(pvec,svec,sxp_pp) use sim_utils implicit none real(rp) svec(3), pvec(3) real(rp) svec_p(3), pvec_p(3), pvec_pp(3), svec_pp(3), pvec_ppp(3) real(rp) Rz(3,3), Ry(3,3), Rzp(3,3) real(rp) sxp(3) real(rp)sxp_p(3), sxp_pp(3) real(rp) sdir(3), khat(3) real(rp) zero/0./ real(rp) theta_z, theta_y, theta_zp !Find the part of the polarization in the x-p plane ! The unit vector normal to the x-p plane is khat = (/zero,-pvec(3),pvec(2)/) khat = khat /sqrt(dot_product(khat,khat)) !the polarization in the x-p plane is sdir = svec-dot_product(svec,khat)*khat svec = sdir /sqrt(dot_product(sdir,sdir)) theta_z = atan2(pvec(2),pvec(1)) Rz = 0. Rz(1,1) = cos(theta_z) Rz(1,2) = sin(theta_z) Rz(2,1) = -Rz(1,2) Rz(2,2) = Rz(1,1) Rz(3,3) = 1 pvec_p = matmul(Rz,pvec) theta_y= atan2(-pvec_p(1),pvec_p(3)) Ry = 0. Ry(1,1) = cos(theta_y) Ry(1,3) = sin(theta_y) Ry(3,1) = -Ry(1,2) Ry(3,3) = Ry(1,1) Ry(2,2) = 1 pvec_pp=matmul(Ry,pvec_p) svec_p = matmul(Ry,matmul(Rz,svec)) theta_zp = atan2(svec_p(2),svec_p(1)) Rzp = 0 Rzp(1,1) = cos(theta_zp) Rzp(1,2) = sin(theta_zp) Rzp(2,1) = -Rzp(1,2) Rzp(2,2) = Rzp(1,1) Rzp(3,3) = 1 pvec_ppp = matmul(Rzp,pvec_pp) svec_pp = matmul(Rzp,svec_p) sxp = (/svec(2)*pvec(3)-svec(3)*pvec(2),svec(3)*pvec(1)-svec(1)*pvec(3),svec(1)*pvec(2)-svec(2)*pvec(1)/) sxp_p = matmul(Rzp,matmul(Ry,matmul(Rz,sxp))) sxp_pp = (/svec_pp(2)*pvec_ppp(3)-svec_pp(3)*pvec_ppp(2), & svec_pp(3)*pvec_ppp(1)-svec_pp(1)*pvec_ppp(3),svec_pp(1)*pvec_ppp(2)-svec_pp(2)*pvec_ppp(1)/) !print '(a,9es12.4)',' RightOrLeft, pvec_ppp, svec_pp,sxp_pp= ',pvec_ppp,svec_pp, sxp_pp sxp_pp(2) = sign(sxp_pp(2),sxp(2)) return end subroutine RightOrLeft subroutine write_phase_space_at_time(co, spin_phase,write_time, bin, write_bins, close_files) use bmad use bmad_struct use bmad_interface use parameters_bmad implicit none type (coord_struct) co real(rp) spin_phase real(rp) write_time real(rp) jx/0./,jy/0./ integer time_nanosec integer, allocatable, save :: lunn(:) integer, allocatable :: write_bins(:) integer bin integer ix integer n_write integer i character*5 time_word character*29, allocatable, save :: file_name(:) logical, allocatable, save :: first(:) logical close_files n_write = size(write_bins) if(close_files)then do i=1,n_write close(unit = lunn(i)) print '(a,i10,a,a)', 'close unit ',lunn(i),' file ', trim(directory)//'/'//file_name(i) end do return endif if(.not. allocated(first))then allocate(first(1:n_write)) first(:) = .true. allocate(lunn(1:n_write)) allocate(file_name(1:n_write)) endif ix=minloc(abs(write_bins(:)-bin),1) if(first(ix)) then time_nanosec = write_time*1.e9 time_word = '00000' if(time_nanosec >=10 .and. time_nanosec<100) write(time_word(4:5),'(i2)')time_nanosec if(time_nanosec >=100 .and. time_nanosec<1000) write(time_word(3:5),'(i3)')time_nanosec if(time_nanosec >=1000 .and. time_nanosec<10000) write(time_word(2:5),'(i4)')time_nanosec if(time_nanosec >=10000) write(time_word,'(i5)')time_nanosec file_name(ix) = 'phase_space_at_'//time_word//'_ns.dat' lunn(ix)=lunget() open(lunn(ix), file = trim(directory)//'/'//file_name(ix)) print '(a,i10,a,a)', 'open unit ',lunn(ix),' file ', trim(directory)//'/'//file_name(ix) write(lunn(ix),'(a,a10,1x,14a16)')'!', 'muon', 'jx','jy', 'x', 'xp','y','yp','vec(5)','pz','t','s','s_x','s_y','s_z','acos(s . p)' first(ix)=.false. endif write(lunn(ix),'(i10,1x, 14es16.8)') muon_number,jx,jy, co%vec, co%t, co%s, co%spin, spin_phase return end subroutine write_phase_space_at_time