Merge pull request #683 from bmad-sim/devel/69

Corrected code for using natural emit for initial beam emit.

bmad/low_level/track_a_foil.f90

@@ -58,7 +58,7 @@ subroutine track_a_foil (orbit, ele, param, mat6, make_matrix)
 
 I_excite = mean_excitation_energy_over_z(z_material) * z_material
 mass_material = mass_of(material) / atomic_mass_unit
-elec_area_density = N_avogadro * ele%value(area_density_used$) * z_material / (1.0e-3_rp * mass_material) ! number_electrons / m^2
+elec_area_density = (1.0e3_rp * N_avogadro) * ele%value(area_density_used$) * z_material / mass_material ! number_electrons / m^2
 
 n_step = nint(ele%value(num_steps$))
 do i = 1, n_step

bmad/modules/bmad_routine_interface.f90

@@ -338,14 +338,6 @@ subroutine c_to_cbar (ele, cbar_mat)
  real(rp) cbar_mat(2,2)
 end subroutine
 
-subroutine calc_emit_from_beam_init (beam_init, ele, species)
- import
- implicit none
- type (beam_init_struct) beam_init
- type (ele_struct) ele
- integer species
-end subroutine
-
 subroutine calc_next_fringe_edge (track_ele, s_edge_body, fringe_info, orbit, init_needed, time_tracking)
  import
  type (ele_struct), target :: track_ele
@@ -2136,6 +2128,16 @@ recursive subroutine set_ele_status_stale (ele, status_group, set_slaves)
  logical, optional :: set_slaves
 end subroutine
 
+function set_emit_from_beam_init (beam_init_in, ele, species, modes, err_flag) result (beam_init_set)
+ import
+ implicit none
+ type (beam_init_struct), target :: beam_init_set, beam_init_in
+ type (ele_struct) ele
+ type (normal_modes_struct), optional :: modes
+ integer species
+ logical, optional :: err_flag
+end function
+
 subroutine set_fringe_on_off (fringe_at, ele_end, on_or_off) 
  import
  implicit none

bmad/multiparticle/set_emit_from_beam_init.f90

@@ -0,0 +1,143 @@
+!+
+! Function set_emit_from_beam_init (beam_init_in, ele, species, modes, err_flag) result (beam_init_set)
+!
+! Routine to put in the beam_init_set structure values that should be used to initialize a beam.
+! For example, if beam_init_in%sig_z is set negative, beam_init_set%sig_z will be set to the value in modes.
+! Also, if emit is set, norm_emit will be computed and vice versa.
+!
+! Input:
+! beam_init_in -- beam_init_struct: Input parameters
+! ele -- ele_struct:
+! species -- integer: Beam particle species.
+! modes -- normal_modes_struct, optional: Normal mode parameters. Used if stuff like beam_init_in%a_emit set negative.
+!
+! Ouput:
+! err_flag -- logical, optional: Set true if there is an error. False otherwise.
+! beam_init_set -- beam_init_struct: See above.
+!-
+
+function set_emit_from_beam_init (beam_init_in, ele, species, modes, err_flag) result (beam_init_set)
+
+use bmad_routine_interface, dummy => set_emit_from_beam_init
+
+implicit none
+
+type (beam_init_struct), target :: beam_init_in, beam_init_set
+type (beam_init_struct), pointer :: bis
+type (ele_struct) ele
+type (normal_modes_struct), optional :: modes
+
+real(rp) ran_g(2), beta_gamma
+integer species
+logical, optional :: err_flag
+
+character(*), parameter :: r_name = 'set_emit_from_beam_init'
+
+! Check
+
+if (present(err_flag)) err_flag = .true.
+bis => beam_init_set
+bis = beam_init_in
+if (ele%value(p0c$) == 0) return
+
+! Sanity checks.
+
+if (bis%a_emit /= 0 .and. bis%a_norm_emit /= 0 .and. .not. (bis%a_emit < 0 .and. bis%a_norm_emit < 0)) then
+ call out_io (s_error$, r_name, 'I am confused! Both a_emit and a_norm_emit are set non-zero in the beam_init_struct structure.', &
+ 'Please set one or the other to zero.')
+ return
+endif
+
+if (bis%b_emit /= 0 .and. bis%b_norm_emit /= 0 .and. .not. (bis%b_emit < 0 .and. bis%b_norm_emit < 0)) then
+ call out_io (s_error$, r_name, 'I am confused! Both b_emit and b_norm_emit are set non-zero in the beam_init_struct structure.', &
+ 'Please set one or the other to zero.')
+ return
+endif
+
+! a-emit
+
+beta_gamma = ele%value(p0c$) / mass_of(species)
+
+if (bis%a_emit < 0 .or. bis%a_norm_emit < 0) then
+ if (present(modes)) then
+ bis%a_emit = modes%a%emittance
+ bis%a_norm_emit = modes%a%emittance * beta_gamma
+ else
+ call out_io (s_warn$, r_name, &
+ 'a_emit and/or a_norm_emit is set negative but the calling program has not provided lattice emittance numbers!', &
+ 'Will use a value of zero...')
+ bis%a_emit = 0
+ bis%a_norm_emit = 0
+ endif
+elseif (bis%a_norm_emit /= 0) then
+ bis%a_emit = bis%a_norm_emit / beta_gamma
+else
+ bis%a_norm_emit = bis%a_emit * beta_gamma
+endif
+
+! b-emit
+
+if (bis%b_emit < 0 .or. bis%b_norm_emit < 0) then
+ if (present(modes)) then
+ bis%b_emit = modes%b%emittance
+ bis%b_norm_emit = modes%b%emittance * beta_gamma
+ else
+ call out_io (s_warn$, r_name, &
+ 'b_emit and/or b_norm_emit is set negative but the calling program has not provided lattice emittance numbers!', &
+ 'Will use a value of zero...')
+ bis%b_emit = 0
+ bis%b_norm_emit = 0
+ endif
+elseif (bis%b_norm_emit /= 0) then
+ bis%b_emit = bis%b_norm_emit / beta_gamma
+else
+ bis%b_norm_emit = bis%b_emit * beta_gamma
+endif
+
+! sig_z
+
+if (bis%sig_z < 0) then
+ if (present(modes)) then
+ bis%sig_z = modes%sig_z
+ else
+ call out_io (s_warn$, r_name, &
+ 'sig_z is set negative but the calling program has not provided lattice emittance numbers!', &
+ 'Will use a value of zero...')
+ bis%sig_z = 0
+ endif
+endif
+
+! sig_pz
+
+if (bis%sig_pz < 0) then
+ if (present(modes)) then
+ bis%sig_pz = modes%sigE_E
+ else
+ call out_io (s_warn$, r_name, &
+ 'sig_pz is set negative but the calling program has not provided lattice emittance numbers!', &
+ 'Will use a value of zero...')
+ bis%sig_pz = 0
+ endif
+endif
+
+! Checking that |bis%dpz_dz| < mode%sigE_E / mode%sig_z
+
+if (abs(bis%dPz_dz * bis%sig_z) > bis%sig_pz) then
+ call out_io (s_error$, r_name, "|dpz_dz| MUST be < mode%sigE_E / mode%sig_z")
+ return
+endif
+
+! Add jitter if needed
+
+if (any(bis%emit_jitter /= 0)) then
+ call ran_gauss(ran_g) ! ran(3:4) for z and e jitter used below
+ bis%a_emit = bis%a_emit * (1 + bis%emit_jitter(1) * ran_g(1))
+ bis%b_emit = bis%b_emit * (1 + bis%emit_jitter(2) * ran_g(2))
+ bis%a_norm_emit = bis%a_norm_emit * (1 + bis%emit_jitter(1) * ran_g(1))
+ bis%b_norm_emit = bis%b_norm_emit * (1 + bis%emit_jitter(2) * ran_g(2))
+endif
+
+if (present(err_flag)) err_flag = .false.
+
+end function set_emit_from_beam_init 
+

bsim/modules/lt_tracking_mod.f90

@@ -121,6 +121,7 @@ module lt_tracking_mod
  type (ltt_section_struct), allocatable :: sec(:) ! Array of sections indexed from 0. The first one marks the beginning.
  type (ele_pointer_struct), allocatable :: ramper(:) ! Ramper element locations.
  type (beam_init_struct) beam_init
+ type (beam_init_struct) beam_init_used
  type (random_state_struct) ramper_ran_state
  integer :: n_ramper_loc = 0
  integer, allocatable :: ix_wake_ele(:) ! List of element indexes where a wake is applied.
@@ -1161,8 +1162,8 @@ subroutine ltt_init_beam_distribution (lttp, ltt_com, beam)
 
 !
 
-call init_beam_distribution (ele_start, lat%param, ltt_com%beam_init, beam, err_flag, &
-  modes = ltt_com%modes, print_p0c_shift_warning = .false.)
+call init_beam_distribution (ele_start, lat%param, ltt_com%beam_init, beam, err_flag, ltt_com%modes, &
+ ltt_com%beam_init_used, print_p0c_shift_warning = .false.)
 if (err_flag) stop
 print '(a, i8)', 'n_particle: ', size(beam%bunch(1)%particle)
 ltt_com%n_particle = size(beam%bunch(1)%particle)
@@ -1475,7 +1476,10 @@ subroutine ltt_setup_high_energy_space_charge(lttp, ltt_com, branch)
 type (branch_struct) branch
 type (ltt_com_struct), target :: ltt_com
 type (normal_modes_struct) modes
+type (ele_struct), pointer :: ele_start, ele_stop
+
 real(rp) n_particle, gamma
+logical err
 
 !
 
@@ -1495,10 +1499,12 @@ subroutine ltt_setup_high_energy_space_charge(lttp, ltt_com, branch)
 gamma = branch%ele(0)%value(e_tot$) / mass_of(branch%ele(0)%ref_species)
 modes = ltt_com%modes
 
-if (ltt_com%beam_init%a_emit > 0) modes%a%emittance = ltt_com%beam_init%a_emit
-if (ltt_com%beam_init%b_emit > 0) modes%b%emittance = ltt_com%beam_init%b_emit
-if (ltt_com%beam_init%a_norm_emit > 0) modes%a%emittance = ltt_com%beam_init%a_norm_emit / gamma
-if (ltt_com%beam_init%b_norm_emit > 0) modes%b%emittance = ltt_com%beam_init%b_norm_emit / gamma
+call ltt_pointer_to_map_ends(lttp, ltt_com%tracking_lat, ele_start, ele_stop)
+ltt_com%beam_init_used = set_emit_from_beam_init(ltt_com%beam_init, ele_start, ele_start%ref_species, ltt_com%modes, err)
+if (err) stop
+
+modes%a%emittance = ltt_com%beam_init_used%a_emit
+modes%b%emittance = ltt_com%beam_init_used%b_emit
 
 if (modes%a%emittance == 0 .or. modes%b%emittance == 0) then
  print *, 'WARNING! No a-mode or b-mode emittance set in beam_init structrue. Cannot compute high energy space charge kick.'