!+ ! Program dark_current_tracker ! ! Program to track multiple particles through an arbitrary lattice ! (backwards or forwards through elements) using a custom time- ! based tracking routine. ! ! Note: Partcles are assumed to be electrons ! ! Modules Needed: ! use bmad ! use beam_def_struct ! ! Input (namelist) ! dark_current_tracker.in ! ! Output ! [particle file name].out !- program dark_current_tracker use bmad use beam_def_struct implicit none type (lat_struct) :: lat type (ele_struct), pointer :: ele type (track_struct) :: single_track type (particle_struct), allocatable :: start_particles(:) type (coord_struct) :: start, middle type (ele_pointer_struct), allocatable :: monitor_eles(:) character(100) :: in_file, lat_name, particle_file_name, monitor_file_name, outfile_name character(400) :: element_monitor_list real(rp) :: e_tot, max_charge, min_charge_ratio, dt_save integer :: outfile, namelist_file, monitor_file integer, parameter :: max_iteration = 100000 integer :: ele_id, particle_id, iteration, point_id integer :: ix_ele, n_monitor_eles logical :: wall_flag, exit_flag logical :: save_tracks, verbose, err namelist / dark_current_tracker_params / lat_name, particle_file_name, save_tracks, & verbose, min_charge_ratio, dt_save, element_monitor_list !------ print '(a100)', '' print '(a100)', ' __ __ __ __ __ ___ ___ ___ __ __ ___ __ ' print '(a100)', '| \ /\ |__) |__/ / ` | | |__) |__) |__ |\ | | | |__) /\ / ` |__/ |__ |__)' print '(a100)', '|__/ /~~\ | \ | \ \__, \__/ | \ | \ |___ | \| | | | \ /~~\ \__, | \ |___ | \' print '(a100)', '' !------------------------------------------ !Defaults for namelist lat_name = 'lat.bmad' particle_file_name = 'ReferenceParticles.dat' save_tracks = .false. verbose = .false. min_charge_ratio = 1e-6 dt_save = 0.001/c_light element_monitor_list = '' !Read namelist in_file = 'dark_current_tracker.in' if (command_argument_count() > 0) call get_command_argument(1, in_file) namelist_file = lunget() print *, 'Opening: ', trim(in_file) open (namelist_file, file = in_file, status = "old") read (namelist_file, nml = dark_current_tracker_params) close (namelist_file) if (save_tracks) then print *, "save_tracks ON; all particle tracks will be saved" else print *, "save_tracks OFF; only final coordinates will be saved" endif !Prepare output file name call file_suffixer (particle_file_name, outfile_name, '.out', .true.) !Parse Lattices call bmad_parser (lat_name, lat) !------------------------------------------ !Track through multiple elements !Import BMAD-T style particles call import_time_distribution(particle_file_name, & start_particles) !Switch all elements to custom tracking do ele_id = 1, lat%n_ele_track lat%ele(ele_id)%tracking_method = custom$ !Also set logic default for monitor use lat%ele(ele_id)%logic = .false. end do !Monitor elements if (element_monitor_list /= '') then call lat_ele_locator (element_monitor_list, lat, monitor_eles, n_monitor_eles, err) if (err) then print *, 'Problem with monitor eles' call err_exit end if print *, 'Monitoring: ' do ix_ele = 1, n_monitor_eles monitor_eles(ix_ele)%ele%logic = .true. print *, ' ', trim( monitor_eles(ix_ele)%ele%name ), ' at s = ', monitor_eles(ix_ele)%ele%s, ' m' end do monitor_file = lunget() call file_suffixer (particle_file_name, monitor_file_name, '.monitor', .true.) open(monitor_file, file = monitor_file_name) write (monitor_file, '(10a19)'), 't', 'x', 'cp_x', 'y', 'cp_y', 's', 'cp_s', 'hit_angle', 'charge', 's_origin' write (monitor_file, '(10a19)'), 's', 'm', 'eV', 'm', 'eV', 'm', 'eV', 'rad', 'C', 'm' !Coordinates will be monitored at the exit end of elements in the loop below end if !Open output file outfile = lunget() open(outfile, file = outfile_name) !Write header write (outfile, '(9a19)'), 't', 'x', 'cp_x', 'y', 'cp_y', 's', 'cp_s', 'hit_angle', 'charge' write (outfile, '(9a19)'), 's', 'm', 'eV', 'm', 'eV', 'm', 'eV', 'rad', 'C' !-------------------- !Iterate through particles to find maximum charge max_charge = maxval(start_particles%charge) !Track electrons lat%param%particle = electron$ !-------------------- !Track list of particles do particle_id = 1, size(start_particles) if (verbose) print *, "Tracking particle: ", particle_id !If particle has charge less than max_charge * min_charge_ratio, move on if (start_particles(particle_id)%charge < max_charge * min_charge_ratio) & cycle !Set start coords and starting element; reset flags start = start_particles(particle_id)%r !if s < 0, drift particle forward to s = 0 if (start%s < 0) call drift_orbit_time(start, mass_of(lat%param%particle), -start%s) call ele_at_s(lat, start%s, ix_ele) wall_flag = .false. exit_flag = .false. !Write to file: header write (outfile, '(a, i8)'), 'particle ', particle_id !Track particle as it traverses up to max_iteration elements do iteration = 1, max_iteration ele => lat%ele(ix_ele) if (verbose) print *, " Tracking element:", ix_ele !Track particle until it hits something (beginning, end, wall) if (save_tracks) then single_track%ds_save = dt_save*c_light call track1( start, ele, lat%param, middle, track = single_track) else call track1( start, ele, lat%param, middle) end if select case (lat%param%particle_at) case (entrance_end$) ix_ele = ix_ele - 1 !Write monitor coordinates if going 'backwards' and the next element is a monitor element if (lat%ele(ix_ele)%logic) write (monitor_file, '(10es19.10E3)') & middle%t, middle%vec, middle%phase_x, & start_particles(particle_id)%charge, start_particles(particle_id)%r%s if (ix_ele .eq. 0) then if (verbose) print *, " Particle exited beginning of lattice" exit_flag = .true. endif case (exit_end$) !Write monitor coordinates at exit if (ele%logic) write (monitor_file, '(10es19.10E3)') & middle%t, middle%vec, middle%phase_x, & start_particles(particle_id)%charge, start_particles(particle_id)%r%s ix_ele = ix_ele + 1 if (ix_ele .eq. lat%n_ele_track + 1) then if (verbose) print *, " Particle exited end of lattice" exit_flag = .true. endif case (no_end$) if (verbose) print *, " Particle lost in open area or & hit wall at element index:", lat%param%ix_lost wall_flag = .true. end select if (verbose) print *, " Number of track points: ", single_track%n_pt !Write to file: track data, if it exists, for particles that didn't ! get lost in an open area if (save_tracks .and. & (.not. (lat%param%particle_at == no_end$ .and. & lat%param%ix_lost == -1))) then do point_id = 0, single_track%n_pt !Change coordinates to lab frame single_track%orb(point_id)%t = & single_track%orb(point_id)%t + & (ele%ref_time - ele%value(delta_ref_time$)) single_track%orb(point_id)%s = & single_track%orb(point_id)%s + (ele%s - ele%value(l$)) single_track%orb(point_id)%vec(5) = single_track%orb(point_id)%s write (outfile, '(9es19.10E3)') & single_track%orb(point_id)%t, & single_track%orb(point_id)%vec,& single_track%orb(point_id)%phase_x, & start_particles(particle_id)%charge end do end if if (save_tracks .and. wall_flag .and. single_track%n_pt <= 2) & print *, "WARNING: particle hit wall immediately: ", particle_id !if particle hit wall or exited lattice, move onto next particle if ((wall_flag) .or. (exit_flag)) then !Write to file: final coordinates, if there is no save_tracks, ! for particles that didn't get lost in an open area if ((.not. save_tracks) .and. & (.not. (lat%param%particle_at == no_end$ .and. & lat%param%ix_lost == -1))) then !Write to file: header !write (outfile, '(a, i8)'), 'particle ', particle_id !Write to file: final point write (outfile, '(9es19.10E3)') & middle%t, middle%vec, middle%phase_x, start_particles(particle_id)%charge end if exit end if !New start coords for next lattice element start = middle end do end do !Close output file close(outfile) !Close monitor file if used if (element_monitor_list /= '') then print *, "Written monitor file: ", monitor_file_name close(monitor_file) end if print *, "Written: ", outfile_name contains !+ ! Subroutine import_time_distribution(dist_file, mc2, orb) ! ! Simple routine to import time-based distribution particles from a file. ! ! Modules Needed: ! use bmad ! use beam_def_struct ! ! Input: ! dist_file -- character(*): File name for time distibution. ! The time format is: ! nparticles(integer) ! x/m m*c^2*\beta_x*\gamma_x/eV y/m m*c^2*\beta_y*\gamma_y z/m m*c^2*\beta_z*\gamma_z/eV time/s charge/C ! . ! . ! ! Output: ! orb(nparticles) -- coord_struct: array containing all of these particles ! !- subroutine import_time_distribution(dist_file, particles) use bmad use beam_def_struct implicit none character(*), intent(in) :: dist_file type (particle_struct), allocatable, intent(out) :: particles(:) real(rp) timevec(8) integer nparticles, i integer, parameter :: dist = 999 !Temporary - better number? open(dist, file = dist_file) read(dist, '(i8)') nparticles if (.not. allocated (particles)) then allocate(particles(1:nparticles)) !write (outfile, '(a)' ) "Start coordinates" !write (outfile, '(6es18.10)') particle_beg%vec !write (outfile, '(a)' ) "End coordinates" !write (outfile, '(6es18.10)') particle_end%vec endif do i = 1, nparticles read(dist, *) timevec !Bmad-T variables use c*px, c*py, c*ps in eV particles(i)%r%vec = [ timevec(1), timevec(2), timevec(3), timevec(4), timevec(5), timevec(6) ] particles(i)%r%t = timevec(7) particles(i)%r%s = timevec(5) particles(i)%charge = timevec(8) end do close(dist) end subroutine import_time_distribution !+ ! Subroutine drift_orbit_time(orbit, mc2, delta_s) ! ! Simple routine to drift a particle orbit in time-based coordinates by a distance delta_s ! ! Modules Needed: ! use bmad_struct ! ! Input: ! orbit -- coord_struct: particle orbit in time-based coordinates ! mc2 -- real(rp): particle mass in eV ! delta_s -- real(rp): s-coordinate distance to drift particle ! . ! ! Output: ! orbit -- coord_struct: particle orbit in time-based coordinates ! !- subroutine drift_orbit_time(orbit, mc2, delta_s) use bmad_struct implicit none type (coord_struct) :: orbit real(rp) :: mc2, delta_s, delta_t, v_s, e_tot, vel(3) e_tot = sqrt( orbit%vec(2)**2 + orbit%vec(4)**2 + orbit%vec(6)**2 + mc2**2) !Get e_tot from momentum vel(1:3) = c_light*[ orbit%vec(2), orbit%vec(4), orbit%vec(6) ]/ e_tot ! velocities v_x, v_y, v_s: c*[c*p_x, c*p_y, c*p_s]/e_tot delta_t = delta_s / vel(3) !Drift x, y, s orbit%vec(1) = orbit%vec(1) + vel(1)*delta_t !x orbit%vec(3) = orbit%vec(3) + vel(2)*delta_t !y orbit%vec(5) = orbit%vec(5) + vel(3)*delta_t !s orbit%s = orbit%s + delta_s orbit%t = orbit%t + delta_t end subroutine drift_orbit_time end program