!spin with constant pitch program spin_check use bmad use spin_test_parameters ! use muon_mod ! use muon_interface ! use parameters_bmad ! use quad_scrape_parameters type (lat_struct),target:: lat type (branch_struct),pointer:: branch type (ele_struct) ele type (coord_struct), allocatable :: co(:) type (coord_struct) orb_at_s procedure(em_field_custom_def) :: em_field_custom real(rp) w,s1,s2, half real(rp) offset real(rp) dist real(rp) f_rev real(rp) delta_e, chrom_x, chrom_y real(rp) delta_B/0./ real(rp) pangle real(rp) s real(rp) px,py,pz real(rp) betahat(1:3) real(rp) beta_dot_s real(rp) deltap real(rp) xoff real(rp) voff real(rp) xp real(rp) theta, beta(1:3),beta_cart(1:3), spin_cart(1:3), sCrossBeta(1:3), sCrossB(1:3) real(rp) B(1:3)/0.,1.,0./ real(rp) gamma real(rp) radius/7.112/ real(rp) beta0,betaz,energy,mass, gamma0 real(rp) gammaz real(rp) betaxy real(rp) betazz real(rp) DELTA real(rp) eta real(rp) symp_err integer ix, i, j, m, ii, n_turns, index, ixx, k integer nargs,iargc integer datetime_values(8) integer n integer ib/0/, track_state integer lun character*16 date, time, zone character*100 lat_file, pitch_angle/'0.5'/ character*100 delta_p/'0.118034'/ character*100 x_offset/'0.0'/ character*100 v_offset/'0.0'/ character*100 x_angle/'0.0'/ character*100 blong/'0.0'/ character*100 brad/'0.0'/ character*100 harmonic/'0'/ character*100 k_quad_volts/'0'/ character*100 edm_eta/'0'/ character*100 Bradial_frac/'1.'/ logical err_flag logical itexists namelist /spin_check_input/lat_file, x_offset, x_angle , v_offset, pitch_angle, delta_p, blong, brad, harmonic, k_quad_volts, edm_eta, Bradial_frac call date_and_time(date,time,zone,datetime_values) print '(a10,a10,a1,a10)',' date = ',date,' time = ', time ! directory = trim(date)//'_'//trim(time(1:6)) ! status= system('mkdir ' // directory) ! if(status == 0)print '(a20,a)',' create directory = ', trim(directory) ! quad_fringe_energy_change = .false. ! quad_fringe_energy_change = .true. em_field_custom_ptr => em_field_custom bmad_com%spin_tracking_on = .true. bmad_com%aperture_limit_on = .false. bmad_com%max_aperture_limit = 10000. nargs = command_argument_count() ! filename = 'test' ! ix_branch = 0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(nargs <=0)then inquire (file='spin_check_input.dat', exist=itexists) if(.not. itexists)print *,'spin_check_input.dat not found' if (itexists)then OPEN (UNIT=5, FILE='spin_check_input.dat', STATUS='old', IOSTAT=ios) READ (5, NML=spin_check_input, IOSTAT=ios) print *, 'ios=', ios rewind(unit=5) CLOSE(5) write(6,nml = spin_check_input) else print *, 'spin_check_input.dat not found' print *, 'Input' print *, 'lattice (bmad.)' print *,'x offset ' print *,'x angle ' print *,'v offset' print *,'pitch angle' print *,' momentum offset ' print *, 'B_longitudinal/Bz ' print *, 'B_radial/Bz ' print *,' longitudinal field harmonic (integer) ' print *,' k_quad (~ -150000)' print *,' edm = 0.' print *,' radial_frac = 1.' stop endif endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(nargs >= 1)then call get_command_argument(1,lat_file) print *, 'Using lattice ',trim(lat_file) if(nargs >= 2)then call get_command_argument(2,x_offset) print *,'Using x offset ', trim(x_offset) endif if(nargs >= 3)then call get_command_argument(3,x_angle) print *,'Using x angle ', trim(x_angle) endif if(nargs >= 4)then call get_command_argument(4,v_offset) print *,'Using v offset', trim(v_offset) endif if(nargs >= 5)then call get_command_argument(5,pitch_angle) print *,'Using pitch angle', trim(pitch_angle) endif if(nargs >= 6)then call get_command_argument(6,delta_p) print *,'Using momentum offset ', trim(delta_p) endif if(nargs >= 7)then call get_command_argument(7,blong) print *, 'B_longitudinal/Bz ', trim(blong) endif if(nargs >= 8)then call get_command_argument(8,brad) print *, 'B_radial/Bz ', trim(brad) endif if(nargs >= 9)then call get_command_argument(9,harmonic) print *, 'harmonic ', trim(harmonic) endif if(nargs >= 10)then call get_command_argument(10,k_quad_volts) print *, 'k_quad ', trim(k_quad_volts) endif if(nargs >= 11)then call get_command_argument(11,edm_eta) print *, 'electric dipole moment: eta=', trim(edm_eta) endif if(nargs >= 12)then call get_command_argument(12,Bradial_frac) print *, 'Fraction with radial B: Bradial_frac=', trim(Bradial_frac) endif else lat_file = 'bmad.' read(x_offset,*)xoff read(pitch_angle,*)betaz read(x_angle,*)xp read(v_offset,*)voff read(delta_p,*)deltap endif read(x_offset,*)xoff read(x_angle,*)xp read(v_offset,*)voff read(pitch_angle,*)betaz read(delta_p,*)deltap read(blong,*)B_long_over_Bz read(brad,*)B_rad_over_Bz read(harmonic,*)nn read(k_quad_volts,*)k_quad read(edm_eta,*)eta read(Bradial_frac,*)radial_frac print '(a,es12.4)', trim(k_quad_volts), k_quad print '(a,es12.4)', trim(edm_eta), eta print '(a,es12.4)' , trim(Bradial_frac), radial_frac print '(2a12,7(a12,es12.4),a12,i10, 3(a12,es12.4))','lat_file = ',lat_file,' x_offset = ',xoff,' x_angle = ',xp,' v_offset = ', voff,' betaz = ',betaz,' delta_p = ',deltap, ' B_long_over Bz =', B_long_over_Bz,' B_rad_over_Bz =', B_rad_over_Bz,' harmonic =',nn,' k_quad =', k_quad,' eta_edm=',eta,' Brad_frac =',radial_frac bmad_com%electric_dipole_moment = eta print '(2a)', ' lat_file = ', lat_file call bmad_parser (lat_file, lat) bmad_com%electric_dipole_moment = eta ! print *,' element offset = ', offset call reallocate_coord (co, lat%n_ele_track) ! do ie=1,lat%branch(ix_branch)%n_ele_track ! print *,ie,lat%branch(ix_branch)%ele(ie)%name ! end do call init_coord(co(0), particle=antimuon$) print *,' mass=',mass_of(co(0)%species) print *,' edm=', eta*.197e-4/4./mass_of(co(0)%species) !e-cm call closed_orbit_calc(lat, co) print '(a,6es12.4)',' Closed orbit ',co(0)%vec call twiss_from_tracking(lat, co(0), symp_err, err_flag) !, status=status, ix_branch=ix_branch) print '(a,es12.4)', 'Energy = ',lat%ele(0)%value(E_TOT$) print '(8a12)',' ','Q_x','Q_y','beta_x','beta_y','alpha_x','alpha_y','eta_x' print '(a,7f12.4)','Closed ring:',lat%a%tune/twopi, lat%b%tune/twopi, lat%ele(0)%a%beta, lat%ele(0)%b%beta, & & lat%ele(0)%a%alpha,lat%ele(0)%b%alpha,lat%ele(0)%a%eta print '(a,es16.8)',' Spin tune = ',lat%param%spin_tune/twopi f_rev=c_light*co(0)%beta/(lat%ele(lat%n_ele_track)%s-co(lat%n_ele_track)%vec(5)) print '(a,es16.8, a, es16.8)',' f_rev = v/l ', f_rev,' f_rev = 1/t ',1./co(lat%n_ele_track)%t print '(a,es16.8,a,es16.8)',' f_a = ', f_rev*lat%param%spin_tune/twopi,' or ', lat%param%spin_tune/twopi/co(lat%n_ele_track)%t co(0)%vec=0 ! print '(a,es12.4)',' Pitch angle = ', pangle co(0)%vec(1) = xoff co(0)%vec(2) = xp mass=mass_of(co(0)%species) ! for fixed transverse momentum increase the total energy gammaz = 1./sqrt(1-betaz**2) energy = lat%ele(0)%value(E_TOT$) gamma0=energy/mass beta0 = (1. - 1./gamma0**2/gammaz**2)**.5 print '(4(a,e12.4))',' mass= ', mass, ' energy=',energy,' gamma =', gamma0,' beta=', beta0 co(0)%vec(4) = betaz/beta0 * gammaz co(0)%vec(6) = deltap+sqrt(1+co(0)%vec(4)**2)-1 co(0)%vec(3) = voff print '(a,6es12.4)',' vec = ', co(0)%vec n=lat%n_ele_track print '(a10, 6a12)','element','s','x','y','z','theta' ! write(11, '(a10, 6a12)')'element','s','x','y','z','theta' do i=1,n print '(i10,6es12.4)',i,lat%ele(i)%s, lat%ele(i)%floor%r, lat%ele(i)%floor%theta ! write(11, '(i10,6es12.4)')i,lat%ele(i)%s, lat%ele(i)%floor%r, lat%ele(i)%floor%theta end do print *,' lat%n_ele_track = ', lat%n_ele_track, n lun=lunget() open(unit=lun,file = 'tbt_'//trim(pitch_angle)//'_bl_bz_'//trim(blong)//'_br_bz_'//trim(brad)//'_'//trim(harmonic)//'_'//trim(k_quad_volts)//'_eta_'//trim(edm_eta)//'_radial_frac_'//trim(Bradial_frac)) print '(a,a)',' Write to file :','tbt_'//trim(pitch_angle)//'_bl_bz_'//trim(blong)//'_br_bz_'//trim(brad)//'_'//trim(harmonic)//'_'//trim(k_quad_volts)//'_eta_'//trim(edm_eta)//'_radial_frac_'//trim(Bradial_frac) write(lun, '(2a12,7(a12,es12.4),a12,i10, a12,es12.4,a12,es12.4)')'lat_file = ',lat_file,' x_offset = ',xoff,' x_angle = ',xp,' v_offset = ', voff, & ' betaz = ',betaz,' delta_p = ',deltap,' B_long_over Bz =', B_long_over_Bz,' B_rad_over_Bz =', B_rad_over_Bz,' harmonic =',nn,' k_quad =', k_quad,' edm eta =',eta write(lun, '(2a10, 7a14,3a14, a14, a15,3a15, a15,9a14,3a14)') 'turn,','ele','s', 'x', 'xp', 'y', 'yp','z','zp', 's_x','s_y','s_z','beta_dot_s','time','beta_x','beta_y','beta_z','theta', & 'beta_x cart','beta_y cart','beta_z cart','spin_x cart','spin_y cart','spin_z cart','sCrossBeta_x','sCrossBeta_y','sCrossBeta_z', & 'spin X B x','spin X B y','spin X B z' co(0)%spin(:) = 0 co(0)%spin(3) = -1. print '(6es14.6,3es14.6)',co(0)%vec, co(0)%spin s=0 do i=1,3000 track_state = 0 call track_all(lat,co, ib, track_state, err_flag) if(co(lat%n_ele_track)%state /= alive$)then print *,'Particle lost at element ', track_state print *,'Lost at ', coord_state_name(co(track_state)%state) print *,' co(track_state)%vec ', co(track_state)%vec exit else do j=1, lat%n_ele_track s = (i-1)*co(lat%n_ele_track)%s + co(j)%s px = co(j)%vec(2) py = co(j)%vec(4) pz = sqrt((co(j)%vec(6)+1)**2 - co(j)%vec(2)**2-co(j)%vec(4)**2) betahat(1:3)=(/px,py,pz/)/sqrt(px**2+py**2+pz**2) beta(1:3) = betahat(1:3) * co(j)%beta theta = co(j)%s/lat%ele(lat%n_ele_track)%s * twopi beta_cart = beta beta_cart(1) = beta(1)*cos(theta) - beta(3)*sin(theta) beta_cart(3) = beta(1)*sin(theta) +beta(3)*cos(theta) beta_dot_s = dot_product(betahat,co(j)%spin) spin_cart = co(j)%spin spin_cart(1) = co(j)%spin(1)*cos(theta) - co(j)%spin(3)*sin(theta) spin_cart(3) = co(j)%spin(1)*sin(theta) + co(j)%spin(3)*cos(theta) sCrossBeta(1) = spin_cart(2)*beta_cart(3) - spin_cart(3)*beta_cart(2) sCrossBeta(2) = spin_cart(3)*beta_cart(1) - spin_cart(1)*beta_cart(3) sCrossBeta(3) = spin_cart(1)*beta_cart(2) - spin_cart(2)*beta_cart(1) sCrossB(1) = -spin_cart(3)*B(2) sCrossB(2) = 0 sCrossB(3) = spin_cart(1)*B(2) gamma = 1./sqrt(1-dot_product(beta_cart,beta_cart)) betaxy=sqrt(beta(1)**2+beta(3)**2) betazz = beta(2) write(lun, '(2i10,7es14.6,3es14.6,es14.6,es14.6,3es14.6, es14.6,9es14.6,3es14.6)')i,j,s,co(j)%vec, co(j)%spin, beta_dot_s, co(j)%t, betahat(1:3), theta, beta_cart, spin_cart, sCrossBeta, & sCrossB write(11, '(i10,3es12.4)')j,co(j)%s,radius*cos(co(j)%s *twopi/lat%ele(lat%n_ele_track)%s),radius*sin(co(j)%s *twopi/lat%ele(lat%n_ele_track)%s) end do co(0) = co(lat%n_ele_track) endif end do close(lun) print '(5(a12,es12.4))', ' The length is ',lat%ele(lat%n_ele_track)%s, ' gamma = ', gamma,' betaxy=', betaxy, ' omega = ',betaxy*c_light/radius,' betaz= ', betazz end !+ ! Subroutine em_field_custom (ele, param, s_rel, orbit, local_ref_frame, field, calc_dfield, err_flag, & ! calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles) ! ! Routine for handling custom (user supplied) EM fields. ! This routine is called when ele%field_calc = custom$ or when ele is a custom element (ele%key = custom$) ! In order to be used, this stub file must be modified appropriately. See the Bmad manual for more details. ! ! Note: Unlike all other elements, "s_rel" and "here" areguments for a patch element are with respect to ! the exit reference frame of the element. See the Bmad manual for more details. ! ! Note: Fields should not have any unphysical discontinuities. ! Discontinuities may cause Runge-Kutta integration to fail resulting in particles getting marked as "lost". ! The mode of failure here is that RK will try smaller and smaller steps to integrate through the ! discontinuity until the step size gets lower than bmad_com%min_ds_adaptive_tracking. At this ! point the particle gets marked as lost. ! ! General rule: Your code may NOT modify any argument that is not listed as ! an output agument below. ! ! Input: ! ele -- Ele_struct: Custom element. ! param -- lat_param_struct: Lattice parameters. ! s_rel -- Real(rp): Longitudinal position relative to the start of the element. ! orbit -- Coord_struct: Coords with respect to the reference particle. ! local_ref_frame -- Logical, If True then take the ! input coordinates and output fields as being with ! respect to the frame of referene of the element. ! calc_dfield -- Logical, optional: If present and True then the field ! derivative matrix is wanted by the calling program. ! use_overlap -- logical, optional: Add in overlap fields from other elements? Default is True. ! calc_potential -- logical, optional: Calc electric and magnetic potentials? Default is false. ! grid_allow_s_out_of_bounds ! -- logical, optional: For grids, allow s-coordinate to be grossly out of bounds ! and return zero instead of an error? Default: False. Used internally for overlapping fields. ! rf_time -- real(rp), optional: Set the time relative to the RF clock. Normally this time is calculated using ! orbit%t or orbit%vec(5) but sometimes it is convenient to be able to override this. ! For example, time_runge_kutta uses this. ! used_eles(:) -- ele_pointer_struct, allocatable, optional: For use if this routine is called recursively. ! This argument should be passed back to em_field_calc if em_field calc is called by this routine. ! ! Output: ! field -- Em_field_struct: Structure hoding the field values. ! err_flag -- Logical, optional: Set true if there is an error. False otherwise. !- recursive subroutine em_field_custom (ele, param, s_rel, orbit, local_ref_frame, field, calc_dfield, err_flag, & calc_potential, use_overlap, grid_allow_s_out_of_bounds, rf_time, used_eles) use em_field_mod !, except_dummy => em_field_custom use spin_test_parameters implicit none type (ele_struct), target :: ele type (lat_param_struct) param type (coord_struct), intent(in) :: orbit type (em_field_struct) :: field type (ele_pointer_struct), allocatable, optional :: used_eles(:) real(rp), intent(in) :: s_rel real(rp), optional :: rf_time real(rp) k real(rp) x,z real(rp) Bl, Bradial, Bz real(rp) radius/7.112/ real(rp) zzlength logical local_ref_frame logical, optional :: calc_dfield, err_flag, calc_potential, grid_allow_s_out_of_bounds, use_overlap character(*), parameter :: r_name = 'em_field_custom' ! Bz=1.4513007 zzlength=length bradial = B_rad_over_bz * Bz !* sin( twopi* orbit%s/zzlength) Bl = B_long_over_Bz * Bz * cos( twopi*nn* orbit%s/zzlength) field%B(1)=0. !if(orbit%s/zzlength<=radial_frac)then field%B(1)= Bradial !endif !field%B(1)= 0. field%B(2) = Bz field%B(3)= Bl !field%B(3)= 0. !print '(2(a,es12.4), a,3es12.4)', 'orbit%s=',orbit%s,' radial_frac = ', radial_frac,' field%B = ', field%B !write(12,'(2es12.4)')orbit%s,Bl k= k_quad !k= -12500000. !(10 times nominal) !k= -36000. !omega_a \sim omega_p !k= -18000. !omega_a \sim 2*omega_p !k=0 field%E(1) = -k * orbit%vec(1) field%E(2) = k * orbit%vec(3) !print '(a7,3es12.4,a2,es12.4,a4,6es12.4)','field%E',field%E,' s', orbit%s,' vec', orbit%vec !if (present(err_flag)) err_flag = .true. end subroutine