program beamline_simulation ! to do: ! after distribution update: ! - check mb_track_beam subroutine & update, if necessary use bmad use beam_mod use spin_mod ! use wake_mod implicit none complex(rp), parameter :: i_imag = (0._rp, 1._rp) integer int1, int2, NconsistWarn, bunch, particle, linenr, elenr, date_time(8), sums3, sums3target, startele, endele real(rp) phi, s3 logical, target :: tracking_error, calc_bunch_error character*8 date character*10 time character*5 zone type (lat_struct) beamline ! holds the lattice info type (beam_struct), target :: beam type (beam_init_struct) beam_init type (bunch_params_struct), target :: bunch_params type init_params_type character(200) latfile, beamfile real(rp) polarization, charge, sigma0, sx_offset, sy_offset, ss_offset, sx_pitch, sy_pitch, stilt, sroll logical human_readable integer verbiage, spin_model, elestart, elemax, FFstart, FFend integer, dimension(1) :: randomseed end type init_params_type type (init_params_type) init_params ! ************************************************************************************************* bmad_com%spin_tracking_on = .true. bmad_com%auto_bookkeeper = .false. ! *************************************** ! *** read in initial beam parameters *** ! *************************************** namelist / init_params_namelist / init_params namelist / initial_values / beam_init open (11, file = "/scratch/mbeckman/bmad/bmad_dist/bmad/Init_beam.txt", iostat = int1, status = "old") if (int1 == 0) then read (11, nml = init_params_namelist) read (11, nml = initial_values) else stop 'STOP: Could not open file with initial parameters: /scratch/mbeckman/bmad/bmad_dist/bmad/Init_beam.txt' endif close (11) ! *********************************** ! *** initialize random generator *** ! *********************************** call date_and_time(date, time, zone, date_time) if( init_params%randomseed(1)<0 ) then init_params%randomseed(1)=60000*24*date_time(7)+1000*24*date_time(6)+24*date_time(8)+date_time(5) endif call random_seed( put=init_params%randomseed ) call ran_seed_put( init_params%randomseed(1) ) ! *********************** ! *** read in lattice *** (calls xsif_parser if latfile ending is .xsif or .XSIF; otherwise calls bmad_parser) ! *********************** int1 = 200 do while (init_params%latfile(int1:int1) == ' ') int1 = int1 - 1 enddo if( init_params%latfile(int1-4:int1) == '.xsif' .or. init_params%latfile(int1-3:int1) == '.XSIF' ) then call xsif_parser (init_params%latfile, beamline) else call bmad_parser (init_params%latfile, beamline) endif startele = init_params%elestart if(init_params%elemax == -1) then endele = beamline%n_ele_track else endele = init_params%elemax endif if(init_params%elestart < 0 .or. init_params%elemax < -1) then do elenr = 0, beamline%n_ele_track if(beamline%ele(elenr)%key==marker$ .and. beamline%ele(elenr)%name=="IP") then if (init_params%elestart==-1) then startele = elenr endif if (init_params%elemax==-3) then endele = elenr endif endif if(beamline%ele(elenr)%key==marker$ .and. beamline%ele(elenr)%name=="MEXFOC") then if (init_params%elemax==-2) then endele = elenr endif exit endif enddo endif ! ************************************* ! replace key == x by key == elementx$ ! *** generate magnet misalignments *** ! does not account for split elements ! ************************************* do elenr = 0, beamline%n_ele_track ! drift$ sbend$ quadrupole$ group$ sextupole$ overlay$ custom$ taylor$ ! rfcavity$ elseparator$ beambeam$ wiggler$ sol_quad$ marker$ kicker$ hybrid$ ! octupole$ rbend$ multipole$ accel_sol$ def_beam$ ab_multipole$ solenoid$ ! patch$ lcavity$ def_parameter$ null_ele$ init_ele$ hom$ match$ monitor$ ! instrument$ hkicker$ vkicker$ rcollimator$ ecollimator$ ! girder$ bend_sol_quad$ def_beam_start$ photon_branch$ branch$ mirror$ crystal$ pipe$ ! if (beamline%ele(elenr)%key == sbend$ .or. beamline%ele(elenr)%key == quadrupole$ .or. ...) then if (elenrinit_params%FFend) then beamline%ele(elenr)%value(x_offset$) = Gauss(0._rp, init_params%sx_offset) beamline%ele(elenr)%value(y_offset$) = Gauss(0._rp, init_params%sy_offset) beamline%ele(elenr)%value(s_offset$) = Gauss(0._rp, init_params%ss_offset) beamline%ele(elenr)%value(x_pitch$) = Gauss(0._rp, init_params%sx_pitch) beamline%ele(elenr)%value(y_pitch$) = Gauss(0._rp, init_params%sy_pitch) beamline%ele(elenr)%value(tilt$) = Gauss(0._rp, init_params%stilt) beamline%ele(elenr)%value(roll$) = Gauss(0._rp, init_params%sroll) else beamline%ele(elenr)%value(x_offset$) = Gauss(0._rp, 0.001*init_params%sx_offset) beamline%ele(elenr)%value(y_offset$) = Gauss(0._rp, 0.001*init_params%sy_offset) beamline%ele(elenr)%value(s_offset$) = Gauss(0._rp, 0.001*init_params%ss_offset) beamline%ele(elenr)%value(x_pitch$) = Gauss(0._rp, 0.001*init_params%sx_pitch) beamline%ele(elenr)%value(y_pitch$) = Gauss(0._rp, 0.001*init_params%sy_pitch) beamline%ele(elenr)%value(tilt$) = Gauss(0._rp, 0.001*init_params%stilt) beamline%ele(elenr)%value(roll$) = Gauss(0._rp, 0.001*init_params%sroll) endif ! endif enddo call lattice_bookkeeper (beamline) ! *********************************************** ! *** propagate Twiss parameters thru lattice *** ! *********************************************** call twiss_propagate_all (beamline) if(init_params%human_readable == .false.) then int1 = beam_init%n_bunch*beam_init%n_particle print '(2A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A)', & 'element/I:x_offset/D:y_offset/D:s_offset/D:x_pitch/D:y_pitch/D:tilt/D:roll/D:charge/D:', & 't_center[', beam_init%n_bunch, ']/D:z_center[', beam_init%n_bunch, ']/D:lost[', int1, & ']/I:x[', int1, ']/D:px[', int1, ']/D:y[', int1, ']/D:py[', int1, ']/D:z[', int1, ']/D:pz[', int1, & ']/D:xhel[', int1, ']/D:yhel[', int1, ']/D:zhel[', int1, ']/D # Branch Descriptor for TTree:ReadFile(...)' ! print '(2A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A)', & ! 'element/I:x_offset/D:y_offset/D:s_offset/D:x_pitch/D:y_pitch/D:tilt/D:roll/D:', & ! 'bunch[', int1, ']/I:particle[', int1, ']/I:lost[', int1, ']/I:x[', int1, ']/D:px[', int1, ']/D:y[', int1, & ! ']/D:py[', int1, ']/D:z[', int1, ']/D:pz[', int1, ']/D:xhel[', int1, ']/D:yhel[', int1, & ! ']/D:zhel[', int1, ']/D # Branch Descriptor for TTree:ReadFile(...)' ! ! 'betaX/D:alphaX/D:gammaX/D:phiX/D:etaX/D:etapX/D:sigmaX/D:sigmpaX/D:emitX/D:normemitX/D:', & ! ! 'betaY/D:alphaY/D:gammaY/D:phiY/D:etaY/D:etapY/D:sigmaY/D:sigmpaY/D:emitY/D:normemitY/D:', & ! ! 'betaS/D:alphaS/D:gammaS/D:phiS/D:etaS/D:etapS/D:sigmaS/D:sigmpaS/D:emitS/D:normemitS/D:', & endif write (*, '(12A)') '# generated (YYYY-MM-DD hh:mm): ',date(1:4),'-',date(5:6),'-',date(7:8),' ',time(1:2),':',time(3:4),' UTC',zone write (*, '(2A)') '# lattice file: ', trim(init_params%latfile) if( init_params%beamfile /= '' ) then write (*, '(2A)') '# beam file: ', trim(init_params%beamfile) endif ! *************************************************** ! *** initialize/read beam (produces text output) *** ! *************************************************** if( init_params%beamfile == '' ) then beam_init%bunch_charge = beam_init%bunch_charge * e_charge ! beam_init%bunch_charge = abs(beam_init%bunch_charge) * e_charge ! before neg. sign was allowed call init_beam_distribution (beamline%ele(startele), beamline%param, beam_init, beam) do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle if(abs(beam%bunch(bunch)%particle(particle)%charge*beam_init%n_particle-beam_init%bunch_charge)>1.e-7) then stop 'STOP: Particle charges initialized unexpectedly' endif enddo enddo else call read_beam (beam, init_params%beamfile, beamline%ele(startele), beam_init) endif write (*, '(A,I13)') '# verbiage: ', init_params%verbiage write (*, '(A,I13)') '# random seed: ', init_params%randomseed write (*, '(A,ES13.6)') '# particle charge (init_beam) [e]: ', init_params%charge write (*, '(A,A15)') '# lattice parameters: particle: ', particle_name(beamline%param%particle) write (*, '(A,ES13.6)') '# mass [eV]: ', mass_of(beamline%param%particle) write (*, '(A,I13)') '# charge [e]: ', charge_of(beamline%param%particle) ! write (*, '(A,ES13.6,2ES15.6)') '# emittance (a,b,z) []: ', beamline%a%emit, beamline%b%emit, & ! beamline%z%emit ! write (*, '(A,ES13.6,2ES15.6)') '# chromaticity (a,b,z) []: ', beamline%a%chrom, beamline%b%chrom, & ! beamline%z%chrom ! write (*, '(A,ES13.6,2ES15.6)') '# beam size (a,b,z) []: ', beamline%a%sigma, beamline%b%sigma, & ! beamline%z%sigma ! write (*, '(A,ES13.6,2ES15.6)') '# beam divergence (a,b,z) []: ', beamline%a%sigmap, beamline%b%sigmap, & ! beamline%z%sigmap ! write (*, '(A,ES13.6)') '# random (Gaussian) magnet misalignment [m]: ', init_params%sigmaalign write (*, '(A)') '# random (Gaussian) magnet misalignments [m]: ' write (*, '(A,ES13.6,2ES15.6)') '# x_offset, y_offset, s_offset: ', init_params%sx_offset, & init_params%sy_offset, init_params%ss_offset write (*, '(A,ES13.6,3ES15.6)') '# x_pitch, y_pitch, tilt, roll: ', init_params%sx_pitch, & init_params%sy_pitch, init_params%stilt, init_params%sroll if (init_params%FFstart /= -1 .or. init_params%FFend /= -1) then write (*, '(A,I13)') '# final focus region start: ', init_params%FFstart write (*, '(A,I13)') '# (misalignments*0.001) end: ', init_params%FFend endif write (*, '(A)') '#' ! ************************************ ! *** manipulate spin of particles *** ! ************************************ if( init_params%beamfile == '' ) then select case(init_params%spin_model) case (0) ! all particles with same spin do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle beam%bunch(bunch)%particle(particle)%r%spin = (/ sqrt( (1._rp+init_params%polarization)*0.5_rp )*(1._rp, 0._rp), & sqrt( (1._rp-init_params%polarization)*0.5_rp )*(1._rp, 0._rp) /) ! (/ (1., 0.), (0., 0.) /) enddo enddo case (1) ! spin Gaussian-distributed around mean polarization with width sigma0, transversal component uniform distributed in phi do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle s3 = Gauss(init_params%polarization, init_params%sigma0) call random_number(phi) phi = 2._rp*pi*phi beam%bunch(bunch)%particle(particle)%r%spin = (/ (1._rp, 0._rp) * sqrt( (1._rp+s3)/2._rp ), & exp( i_imag*phi ) * sqrt( (1._rp-s3)/2._rp ) /) enddo enddo case (2) ! spin up and down, exact value for polarization, last spin chosen to match pol. value perfectly sums3 = 0 sums3target = idnint(init_params%polarization*beam_init%n_particle) do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle-1 call random_number(s3) if ( (sums3-1 < sums3target-(beam_init%n_particle-particle)) .or. & ((s3 >= (1._rp-init_params%polarization)/2._rp) .and. & .not.(sums3+1 > sums3target+(beam_init%n_particle-particle))) ) then beam%bunch(bunch)%particle(particle)%r%spin = (/ (1._rp, 0._rp), (0._rp, 0._rp) /) ! hel = +1 sums3 = sums3 + 1 else beam%bunch(bunch)%particle(particle)%r%spin = (/ (0._rp, 0._rp), (1._rp, 0._rp) /) ! hel = -1 sums3 = sums3 - 1 endif enddo ! adjustment of last spin per bunch to match polarization value exactly s3 = init_params%polarization*beam_init%n_particle-sums3 if(s3>1. .or. s3<-1.) then stop 'STOP: Spin manipulation: (s3>1 or s3<-1)' endif call random_number(phi) beam%bunch(bunch)%particle(beam_init%n_particle)%r%spin = (/ (1._rp, 0._rp) * sqrt( (1._rp+s3)/2._rp ), & exp( i_imag*phi ) * sqrt( (1._rp-s3)/2._rp ) /) enddo case (3) ! spin up and down, non-exact value for polarization sums3 = 0 sums3target = idnint(init_params%polarization*beam_init%n_particle) do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle call random_number(s3) if ( (sums3-1 < sums3target-(beam_init%n_particle-particle)) .or. & ((s3 >= (1._rp-init_params%polarization)/2._rp) .and. & .not.(sums3+1 > sums3target+(beam_init%n_particle-particle))) ) then beam%bunch(bunch)%particle(particle)%r%spin = (/ (1._rp, 0._rp), (0._rp, 0._rp) /) ! hel = +1 sums3 = sums3 + 1 else beam%bunch(bunch)%particle(particle)%r%spin = (/ (0._rp, 0._rp), (1._rp, 0._rp) /) ! hel = -1 sums3 = sums3 - 1 endif enddo enddo case default stop 'STOP: init_params%spin_model: Value out of range' end select endif write (*, '(A,I13)') '# beam parameters: # bunches: ', beam_init%n_bunch write (*, '(A,I13)') '# # sim. particles per bunch: ', beam_init%n_particle if( init_params%beamfile == '' ) then write (*, '(A,ES13.6)') '# a-mode emittance [m] (normalized): ', beam_init%a_norm_emitt write (*, '(A,ES13.6)') '# a bunch emittance rms jitter (norm.): ', beam_init%emitt_jitter(1) write (*, '(A,ES13.6)') '# b-mode emittance [m] (normalized): ', beam_init%b_norm_emitt write (*, '(A,ES13.6)') '# b bunch emittance rms jitter (norm.): ', beam_init%emitt_jitter(2) write (*, '(A,ES13.6)') '# dPz_dz: ', beam_init%dPz_dz write (*, '(A,ES13.6)') '# time between bunches [s]: ', beam_init%dt_bunch write (*, '(A,ES13.6)') '# z sigma [m]: ', beam_init%sig_z write (*, '(A,ES13.6)') '# e_sigma in dE/E: ', beam_init%sig_e endif write (*, '(A,ES13.6)') '# bunch charge [C]: ', beam_init%bunch_charge write (*, '(A,ES13.6)') '# bunch charge [e]: ', & (beam_init%bunch_charge/e_charge) if( init_params%beamfile == '' ) then write (*, '(A,ES13.6)') '# bunch length RMS jitter: ', beam_init%sig_z_jitter write (*, '(A,ES13.6)') '# energy spread RMS jitter: ', beam_init%sig_e_jitter write (*, '(A,ES13.6,2ES15.6)') '# bunch center offset rel. to reference: ', beam_init%center(1:3) write (*, '(A,ES13.6,2ES15.6)') '# ', beam_init%center(4:6) write (*, '(A,ES13.6,2ES15.6)') '# bunch center rms jitter: ', & beam_init%center_jitter(1:3) write (*, '(A,ES13.6,2ES15.6)') '# ', & beam_init%center_jitter(4:6) write (*, '(A)') '#' write (*, '(A,I9)') '# misc settings: spin model: ', init_params%spin_model write (*, '(A,F9.6)') '# beam polarization: ', init_params%polarization if ( init_params%spin_model == 1 ) then write(*, '(A,F9.6)') '# sigma: ', init_params%sigma0 endif endif write (*, '(A)') '##' write (*, '(A)') '## end of header' write (*, '(A)') '##' NconsistWarn = 0 linenr = 100 ! ******************************* ! *** track beam thru lattice *** ! ******************************* !!! call track_beam (beamline, beam, 0, beamline%n_ele_track, tracking_error) ! zero_lr_wakes_in_lat (beamline) calc_bunch_error = 0 call mb_track_beam (beamline, beam, beamline%ele(startele), beamline%ele(endele), err=tracking_error) if (tracking_error) then write (*, '(A)') '! Beam tracking error occured' endif if (calc_bunch_error) then write (*, '(A)') '! Calc_bunch_params: error occured' endif if ( NconsistWarn /= 0 ) then write (*, '(A,I,A)') '# Consistency check failed in ', NconsistWarn, ' cases' endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ******************************************************************************************* ! ! ************************************* SUBROUTINES ************************************* ! ! ******************************************************************************************* ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! contains ! ***************************************** ! *** mb_track_beam: track thru lattice *** ! *** (modified version of track_beam) *** ! ***************************************** ! Subroutine track_beam (lat, beam, ele1, ele2, err) ! ! Subroutine to track a beam of particles from the end of ! ele1 Through to the end of ele2. Both must be in the same lattice branch. ! ! Note: zero_lr_wakes_in_lat needs to be called to initial wakes before tracking. ! ! Modules needed: ! use beam_mod ! ! Input: ! lat -- lat_struct: Lattice to track through. ! beam -- Beam_struct: Beam at end of element ix1. ! ele1 -- Ele_struct, optional: Starting element (this element ! is NOT tracked through). Default is lat%ele(0). ! ele2 -- Ele_struct, optional: Ending element. ! Default is lat%ele(lat%n_ele_track). ! ! Output: ! beam -- beam_struct: Beam at end of element ix2. ! err -- Logical: Set true if there is an error. ! EG: Too many particles lost for a CSR calc. !- subroutine mb_track_beam (lat, beam, ele1, ele2, err) implicit none type (lat_struct), target :: lat type (beam_struct) :: beam type (branch_struct), pointer :: branch type (ele_struct), optional, target :: ele1, ele2 type (ele_struct), pointer :: e1, e2 integer i, j logical err ! Init e1 => lat%ele(0) if (present(ele1)) e1 => ele1 e2 => lat%ele(lat%n_ele_track) if (present(ele2)) e2 => ele2 ! Loop over all elements in the lattice branch => lat%branch(e1%ix_branch) call writeout (e1%ix_ele) do i = e1%ix_ele+1, e2%ix_ele call track1_beam (beam, lat, branch%ele(i), beam, err) if (err) return call writeout (i) if (calc_bunch_error) return enddo end subroutine mb_track_beam ! ********************************************************* ! *** writeout: write output for each beam-line element *** ! ********************************************************* subroutine writeout (i) integer i, bunch, particle, component ! ************************************************************ ! *** decides wheter to write output at all *** ! *** see init_beam.txt for explanation of verbiage level *** ! ************************************************************ select case(init_params%verbiage) case (0) if(beamline%ele(i)%key/=marker$ .and. beamline%ele(i)%key/=init_ele$) then return else if(beamline%ele(i)%key==marker$ .and. beamline%ele(i)%name/="POL" .and. beamline%ele(i)%name/="IP" .and. beamline%ele(i)%name/="MEXFOC") then return endif endif case (1) if(beamline%ele(i)%value(l$)==0. .and. (beamline%ele(i)%key==monitor$ .or. beamline%ele(i)%key==instrument$ .or. & (beamline%ele(i)%key==marker$ .and. beamline%ele(i)%name/="POL" .and. beamline%ele(i)%name/="IP" .and. beamline%ele(i)%name/="MEXFOC"))) then return endif case (3) case default stop 'STOP: init_params%verbiage: Value out of range' end select ! ******************** ! *** write output *** ! ******************** ! element info if(init_params%human_readable) then print '(A)', '#' print '(2A)', '# # key type name s[m] E[GeV] p[GeV] Xoff[mu] Yoff[mu] Soff[mu] Xpitch Ypitch', & ' tilt roll' ! print '(A,I10)', '# ', i print '(A,2I4,x,a12,x,a12,f8.2,2f10.5,x,f8.3,x,f8.3,x,f8.3,x,4f7.2)', '# ', i, beamline%ele(i)%key, & key_name(beamline%ele(i)%key), beamline%ele(i)%name, beamline%ele(i)%s, beamline%ele(i)%value(e_tot$)/1.e9, & beamline%ele(i)%value(p0c$)/1.e9, beamline%ele(i)%value(x_offset$)/1.e-6, beamline%ele(i)%value(y_offset$)/1.e-6, & beamline%ele(i)%value(s_offset$)/1.e-6, beamline%ele(i)%value(x_pitch$), beamline%ele(i)%value(y_pitch$), & beamline%ele(i)%value(tilt$), beamline%ele(i)%value(roll$) print '(A)', '#' print '(A)', "# x: beta alpha gamma phi eta eta'" print '(A,6ES12.4)', '# ', beamline%ele(i)%a%beta, beamline%ele(i)%a%alpha, beamline%ele(i)%a%gamma, beamline%ele(i)%a%phi, & beamline%ele(i)%a%eta, beamline%ele(i)%a%etap print '(A)', '#' print '(A)', "# y: beta alpha gamma phi eta eta'" print '(A,6ES12.4)', '# ', beamline%ele(i)%b%beta, beamline%ele(i)%b%alpha, beamline%ele(i)%b%gamma, beamline%ele(i)%b%phi, & beamline%ele(i)%b%eta, beamline%ele(i)%b%etap print '(A)', '#' print '(A)', "# z: beta alpha gamma phi eta eta'" print '(A,6ES12.4)', '# ', beamline%ele(i)%z%beta, beamline%ele(i)%z%alpha, beamline%ele(i)%z%gamma, beamline%ele(i)%z%phi, & beamline%ele(i)%z%eta, beamline%ele(i)%z%etap print '(A)', '#' linenr = 100 else write (*, '(I4,7(x,ES23.16))', advance="no"), i, beamline%ele(i)%value(x_offset$), beamline%ele(i)%value(y_offset$), beamline%ele(i)%value(s_offset$), & beamline%ele(i)%value(x_pitch$), beamline%ele(i)%value(y_pitch$), beamline%ele(i)%value(tilt$), beamline%ele(i)%value(roll$) !, & endif ! bunch info if(init_params%human_readable) then do bunch = 1, beam_init%n_bunch call calc_bunch_params (beam%bunch(bunch), beamline%ele(i), beamline%param, bunch_params, calc_bunch_error, .false.) if (calc_bunch_error) return print '(A)', '#' print '(A)', '# calc_bunch_params: TWISS:' print '(A)', "# x: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%x%beta, bunch_params%x%alpha, bunch_params%x%gamma, bunch_params%x%phi, & bunch_params%x%eta, bunch_params%x%etap, bunch_params%x%norm_emit print '(A)', '#' print '(A)', "# y: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%y%beta, bunch_params%y%alpha, bunch_params%y%gamma, bunch_params%y%phi, & bunch_params%y%eta, bunch_params%y%etap, bunch_params%y%norm_emit print '(A)', '#' print '(A)', "# z: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%z%beta, bunch_params%z%alpha, bunch_params%z%gamma, bunch_params%z%phi, & bunch_params%z%eta, bunch_params%z%etap, bunch_params%z%norm_emit print '(A)', '#' print '(A)', "# a: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%a%beta, bunch_params%a%alpha, bunch_params%a%gamma, bunch_params%a%phi, & bunch_params%a%eta, bunch_params%a%etap, bunch_params%a%norm_emit print '(A)', '#' print '(A)', "# b: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%b%beta, bunch_params%b%alpha, bunch_params%b%gamma, bunch_params%b%phi, & bunch_params%b%eta, bunch_params%b%etap, bunch_params%b%norm_emit print '(A)', '#' print '(A)', "# c: beta alpha gamma phi eta eta' norm. emit" print '(A,7ES12.4)', '# ', bunch_params%c%beta, bunch_params%c%alpha, bunch_params%c%gamma, bunch_params%c%phi, & bunch_params%c%eta, bunch_params%c%etap, bunch_params%c%norm_emit ! bunch_params%c%sigma, bunch_params%c%sigma_p, bunch_params%c%emit print '(A)', '#' print '(A)', "# SIGMA: 11 12 13 14 15 16 22 23 24 25 26" print '(A,11ES12.4)', '# ', bunch_params%sigma(1), bunch_params%sigma(2), bunch_params%sigma(3), bunch_params%sigma(4), bunch_params%sigma(5), & bunch_params%sigma(6), bunch_params%sigma(7), bunch_params%sigma(8), bunch_params%sigma(9), bunch_params%sigma(10), & bunch_params%sigma(11) print '(A)', "# 33 34 35 36 44 45 46 55 56 66" print '(A,10ES12.4)', '# ', bunch_params%sigma(12), bunch_params%sigma(13), bunch_params%sigma(14), bunch_params%sigma(15), & bunch_params%sigma(16), bunch_params%sigma(17), bunch_params%sigma(18), bunch_params%sigma(19), & bunch_params%sigma(20), bunch_params%sigma(21) print '(A)', '#' print '(A)', "# centroid: x px y py z pz s" print '(A,7ES12.4)', '# ', bunch_params%centroid%vec, bunch_params%s print '(A)', '#' enddo endif ! particle info if(init_params%human_readable) then do bunch = 1, beam_init%n_bunch print '(2A)', '# bunch ix_bunch t_center z_center ix_ele charge' write (*, '(A,I9,x,I9,x,ES10.3,x,ES10.3,x,I7,x,ES10.3)') '# ', bunch, beam%bunch(bunch)%ix_bunch, & beam%bunch(bunch)%t_center, beam%bunch(bunch)%z_center, beam%bunch(bunch)%ix_ele, beam%bunch(bunch)%charge linenr = 100 do particle = 1, beam_init%n_particle if (linenr == 100) then print '(2A)', '# particle lost x[mm] px*1E3 y[mm] py*1E3 z[mm] pz*1E3 ', & ' xhel yhel zhel' linenr = 0 endif write (*, '(A,I12,x,I5,x,6f9.4,3f7.3)') '# ', particle, beam%bunch(bunch)%particle(particle)%ix_lost, & (1000._rp* beam%bunch(bunch)%particle(particle)%r%vec), spinvec(beam%bunch(bunch)%particle(particle)%r%spin) linenr = linenr + 1 enddo enddo else write (*, '(x,ES23.16)', advance="no"), beam_init%bunch_charge !beam%bunch(1)%charge do bunch = 1, beam_init%n_bunch write (*, '(x,ES23.16)', advance="no"), beam%bunch(bunch)%t_center enddo do bunch = 1, beam_init%n_bunch write (*, '(x,ES23.16)', advance="no"), beam%bunch(bunch)%z_center enddo do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle write (*, '(x,I0)', advance="no"), beam%bunch(bunch)%particle(particle)%ix_lost enddo enddo do component = 1, 6 do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle write (*, '(x,ES23.16)', advance="no"), beam%bunch(bunch)%particle(particle)%r%vec(component) enddo enddo enddo do component = 1, 3 do bunch = 1, beam_init%n_bunch do particle = 1, beam_init%n_particle write (*, '(x,ES23.16)', advance="no"), spinvecn(beam%bunch(bunch)%particle(particle)%r%spin, component) enddo enddo enddo write (*, '(A)'), '' ! new line endif end subroutine writeout ! **************************************** ! *** Gauss: Gaussian random generator *** ! **************************************** function Gauss (mean0, sigma0) use random_mod real(rp), optional, intent(in) :: mean0 real(rp), optional, intent(in) :: sigma0 real(rp) mean, sigma, Gauss mean = 0._rp if (present(mean0)) mean = mean0 sigma = 1._rp if (present(sigma0)) sigma = sigma0 call ran_gauss(Gauss) Gauss = sigma * Gauss + mean end function Gauss ! ******************************************************** ! *** spinvec: computes 3-dim. spin-vector from spinor *** ! ******************************************************** function spinvec (spinor) complex (rp) :: spinor(2) real (rp) :: spinvec(3) ! real (rp) sv2 ! xhel[j] = 2.*(re_plus*re_minus+im_plus*im_minus); ! yhel[j] = 2.*(re_plus*im_minus-im_plus*re_minus); ! zhel[j] = re_plus*re_plus+im_plus*im_plus-re_minus*re_minus-im_minus*im_minus; spinvec(1) = 2._rp*(real(spinor(1))*real(spinor(2))+aimag(spinor(1))*aimag(spinor(2))) spinvec(2) = 2._rp*(real(spinor(1))*aimag(spinor(2))-aimag(spinor(1))*real(spinor(2))) spinvec(3) = real(spinor(1))**2+aimag(spinor(1))**2-real(spinor(2))**2-aimag(spinor(2))**2 ! sv2 = abs(spinvec) ! if( abs(abs(spinvec)-1._rp)>1.e-16 ) then ! print '(A,F)', '# Consistency check failed: |\vec{spin}| = ', abs(spinvec) ! stop ! endif end function spinvec ! ***************************************************************** ! *** spinvecn: computes n-th spin-vector component from spinor *** ! ***************************************************************** function spinvecn (spinor, n) complex (rp) :: spinor(2) real (rp) :: spinvecn integer n select case(n) case (1) spinvecn = 2._rp*(real(spinor(1))*real(spinor(2))+aimag(spinor(1))*aimag(spinor(2))) case (2) spinvecn = 2._rp*(real(spinor(1))*aimag(spinor(2))-aimag(spinor(1))*real(spinor(2))) case (3) spinvecn = real(spinor(1))**2+aimag(spinor(1))**2-real(spinor(2))**2-aimag(spinor(2))**2 case default spinvecn = log(-1.) ! NaN end select end function spinvecn ! ********************************************************************* ! *** read_beam: read beam_struct from text file (format see below) *** ! ********************************************************************* subroutine read_beam (beam, beamfile, ele, beaminit) ! Example beam file: ! some lines of comments ! not longer than 200 characters each ! and concluded by an empty line (even if there are no comments!) ! ! 2, 2, -3.2E-9 ! 0., 0. ! -1, 0., 0., 0., 0., 0., 0., 0., 0., -1. ! -1, 1.23, 2.34, -3.45, -4.56, 5.67, 6.7E1, 1., 0., 0. ! 1.E-6, 2.9E+2 ! -1, 0., 0., 0., 0., 0., 0., 0., 0., -1. ! +5, 1.23, 2.34, -3.45, -4.56, 5.67, 6.7E1, 1., 0., 0. ! use beam_utils use spin_mod implicit none type (ele_struct) ele type (beam_struct), target :: beam type (beam_init_struct) beaminit character(200) beamfile character(200) :: str type (bunch_struct), pointer :: bunch type (particle_struct), pointer :: p real(rp), allocatable :: ri_t_center(:), ri_z_center(:), ri_vec(:,:,:), ri_spin(:,:,:) real(rp) ri_charge integer, allocatable :: ri_lost(:,:) integer i_bunch, i_particle, n_bunch, n_particle, stat open (12, file = beamfile, iostat = stat, status = "old") if (stat == 0) then write(*,'(A)') '# *** description from beam file ***' do read(12, '(A)', iostat=stat) str if (stat /= 0) exit if ( len(trim(str))==0 ) exit write(*,'(2A)') '# ', trim(str) end do if (stat /= 0) stop 'STOP: Error while reading file with initial beam information (comments)' write(*,'(A)') '# *** end of description ***' write(*,'(A)') '#' read(12, *, iostat=stat) n_bunch, n_particle, ri_charge if (stat /= 0) stop 'STOP: Error while reading file with initial beam information (first line)' ! write (*, '(A,I,X,I,X,ES)') '##### ', n_bunch, n_particle, ri_charge allocate(ri_t_center(n_bunch)) allocate(ri_z_center(n_bunch)) allocate(ri_vec(6,n_bunch,n_particle)) allocate(ri_spin(3,n_bunch,n_particle)) allocate(ri_lost(n_bunch,n_particle)) do i_bunch = 1, n_bunch read(12, *, iostat=stat) ri_t_center(i_bunch), ri_z_center(i_bunch) if (stat /= 0) exit do i_particle = 1, n_particle read(12, *, iostat=stat) ri_lost(i_bunch,i_particle), ri_vec(:,i_bunch,i_particle), & ri_spin(:,i_bunch,i_particle) if (stat /= 0) exit ! print '(2A)', '# bunch t_center z_center charge' ! write (*, '(A,I9,x,ES10.3,x,ES10.3,x,ES10.3)') '# ', i_bunch, & ! ri_t_center(i_bunch), ri_z_center(i_bunch), ri_charge ! print '(2A)', '# lost x[m] px*1E0 y[m] py*1E0 z[m] pz*1E0 ', & ! ' xhel yhel zhel' ! write (*, '(A,I5,9(x,E10.3))') '# ', ri_lost(i_bunch,i_particle), & ! ri_vec(:,i_bunch,i_particle), ri_spin(:,i_bunch,i_particle) ! stop enddo enddo ! do i_bunch = 1, n_bunch ! print '(2A)', '# bunch t_center z_center charge' ! write (*, '(A,I9,x,ES10.3,x,ES10.3,x,ES10.3)') '# ', i_bunch, & ! ri_t_center(i_bunch), ri_z_center(i_bunch), ri_charge ! linenr = 100 ! ! do i_particle = 1, n_particle ! if (linenr == 100) then ! print '(2A)', '# particle lost x[m] px*1E0 y[m] py*1E0 z[m] pz*1E0 ', & ! ' xhel yhel zhel' ! linenr = 0 ! endif ! write (*, '(A,I12,x,I5,x,6f9.4,3f7.3)') '# ', i_particle, ri_lost(i_bunch,i_particle), & ! ri_vec(:,i_bunch,i_particle), ri_spin(:,i_bunch,i_particle) ! linenr = linenr + 1 ! enddo ! enddo call reallocate_beam (beam, n_bunch, n_particle) do i_bunch = 1, n_bunch bunch => beam%bunch(i_bunch) bunch%charge = ri_charge bunch%z_center = ri_z_center(i_bunch) bunch%t_center = ri_t_center(i_bunch) bunch%ix_ele = ele%ix_ele bunch%ix_bunch = i_bunch do i_particle = 1, n_particle p => bunch%particle(i_particle) p%r%vec = ri_vec(:,i_bunch,i_particle) call vec_to_spinor (ri_spin(:,i_bunch,i_particle), p%r) ! write (*, '(A,3(x,F))') '# vec:', vec ! write (*, '(A,3(x,F))') '# polar:', polar%theta, polar%phi, polar%xi p%charge = bunch%charge / n_particle p%ix_lost = ri_lost(i_bunch,i_particle) enddo enddo deallocate(ri_t_center, ri_z_center, ri_vec, ri_spin, ri_lost) if (stat /= 0) stop 'STOP: Error while reading file with initial beam information' else stop 'STOP: Could not open file with initial beam information' endif close (12) beam_init%n_bunch = n_bunch beam_init%n_particle = n_particle beam_init%bunch_charge = ri_charge end subroutine read_beam ! *** backup *** ! from bmad_dist_2010_0713_d ! ! ! ! subroutine mb_track_beam (err, ix1, ix2) ! ! ! ! ! ! implicit none ! ! ! ! ! ! integer, optional, intent(in) :: ix1, ix2 ! ! ! integer i, i1, i2, bunch, particle ! ! ! ! real :: last_s = -1.e10 ! for verbiage == 1 ! ! ! ! complex :: last_spin(2,beam_init%n_bunch,beam_init%n_particle) ! ! ! ! ! ! logical err ! ! ! ! ! ! ! Init ! ! ! ! ! ! i = 0 ! ! ! i1 = 0 ! ! ! if (present(ix1)) i1 = ix1 ! ! ! i2 = beamline%n_ele_track ! ! ! if (present(ix2)) i2 = ix2 ! ! ! ! ! ! ! ! initialize last_spin ! ! ! ! do bunch = 1, beam_init%n_bunch ! ! ! ! do particle = 1, beam_init%n_particle ! ! ! ! last_spin(:,bunch,particle) = (/ (2., 2.), (2., 2.) /) ! ! ! ! enddo ! ! ! ! enddo ! ! ! ! ! ! ! call writeout (i, last_s, last_spin) ! ! ! call writeout (i) ! ! ! ! ! ! ! Loop over all elements in the lattice ! ! ! do i = i1+1, i2 ! ! ! call track1_beam (beam, beamline, beamline%ele(i), beam, err) ! ! ! if (err) return ! ! ! ! call writeout (i, last_s, last_spin) ! ! ! call writeout (i) ! ! ! enddo ! ! ! ! ! ! end subroutine mb_track_beam end program