Merge pull request #572 from bmad-sim/devel/step28

Start devel of ASCII::4 beam storage format.

bmad/code/lat_sanity_check.f90

@@ -267,6 +267,15 @@ subroutine lat_sanity_check (lat, err_flag)
       endif
     endif
 
+    if (has_attribute(ele, 'X_PITCH_TOT')) then
+      if (abs(ele%value(x_pitch_tot$)) + abs(ele%value(y_pitch_tot$)) > 0.5*pi) then
+        call out_io (s_fatal$, r_name, &
+              'HAVING |X_PITCH_TOT| + |Y_PITCH_TOT| > PI/2 FOR AN ELEMENT (' // trim(ele%name) // ') DOES NOT MAKE SENSE ', &
+              'SINCE PARTICLES WILL NOT BE MOVING IN THE RIGHT DIRECTION WITHIN THE ELEMENT.')
+        err_flag = .true.
+      endif
+    endif
+
     ! Do not check the extra elements temporarily inserted by bmad_parser2.
 
     select case (ele%key)

bmad/doc/beam-init.tex

@@ -322,7 +322,7 @@ \section{File Based Beam Initialization}
   [bunch loop: ib = 1 to n_bunch]
     BEGIN_BUNCH    ! Marker to mark the beginning of a bunch specification block.
     <species_name> ! Species of particle
-    <bunch_charge> ! Charge of bunch. 0 => Use <macro_charge>.
+    <charge_tot>   ! Total charge of bunch (alive + dead). 0 => Use <macro_charge>.
     <z_center>     ! z position at center of bunch.
     <t_center>     ! t position at center of bunch.
     [particle loop: Stop when END_BUNCH marker found]
@@ -338,7 +338,7 @@ \section{File Based Beam Initialization}
   25000   ! n_particle
   BEGIN_BUNCH
     POSITRON
-    3.2E-9   ! bunch_charge
+    3.2E-9   ! charge_tot
     0.0      ! z_center
     0.0      ! t_center
    -6.5E-3  9.6E-3 -1.9E-2  8.8E-3  2.2E-2 -2.4E-2  1.2E-13  1 1.0 0.0 0.0
@@ -366,7 +366,7 @@ \section{File Based Beam Initialization}
 After this, there are \vn{<n_bunch>} blocks of data, one for each bunch. Each one of these blocks
 starts with a \vn{BEGIN_BUNCH} line to mark the beginning of the block and ends with a
 \vn{END_BUNCH} marker line. In between, the first four lines give the \vn{<species_name>} name,
-\vn{<bunch_charge>}, \vn{<z_center>}, and \vn{<t_center>} values. The \vn{<species_name>} name may
+\vn{<charge_tot>}, \vn{<z_center>}, and \vn{<t_center>} values. The \vn{<species_name>} name may
 be one of:
 \begin{example}
   positron  ! default
@@ -381,7 +381,7 @@ \section{File Based Beam Initialization}
 The lines following the \vn{<t_center>} line specify particle coordinates. One line for each
 particle.  Only the first six numbers, which are the phase space coordinates, need to be specified
 for each particle. If the \vn{<macro_charge>} column is not present, or is zero, it defaults to
-\vn{<bunch_charge>/<n_particle>}.
+\vn{<charge_tot>/<n_particle>}.
 
 The \vn{<state>} parameter indicates whether a particle is alive or dead. Values are
 \begin{example}
@@ -394,9 +394,9 @@ \section{File Based Beam Initialization}
 
 The particle spin is specified by $x$, $y$ and $z$ components.
 
-Each particle has an associated \vn{<macro_charge>}. If \vn{<bunch_charge>} is set to a non-zero
+Each particle has an associated \vn{<macro_charge>}. If \vn{<charge_tot>} is set to a non-zero
 value, the charge of all the particles will be scaled by a factor to make the total macro charge
-equal to \vn{<bunch_charge>}. The macro charge is ignored in tracking. The charge of the particle
+equal to \vn{<charge_tot>}. The macro charge is ignored in tracking. The charge of the particle
 used in tracking is the charge as calculated for the particle species. On the other hand, the macro
 charge is used to calculate such things as the total charge in a particular region or the field
 produced by a particle. That is, the macro charge acts as a weighting factor for a particle when the

bmad/modules/bmad_routine_interface.f90

@@ -1875,11 +1875,12 @@ subroutine reallocate_beam (beam, n_bunch, n_particle, save)
   logical, optional :: save
 end subroutine
 
-subroutine reallocate_bunch (bunch, n_particle)
+subroutine reallocate_bunch (bunch, n_particle, save)
   import
   implicit none
   type (bunch_struct) bunch
   integer n_particle
+  logical, optional :: save
 end subroutine
 
 function relative_mode_flip (ele1, ele2)

bmad/modules/bmad_struct.f90

@@ -495,8 +495,13 @@ module bmad_struct
 integer, parameter :: lost$ = 2
 integer, parameter :: lost_neg_x_aperture$ = 3, lost_pos_x_aperture$ = 4 
 integer, parameter :: lost_neg_y_aperture$ = 5, lost_pos_y_aperture$ = 6
-integer, parameter :: lost_pz_aperture$ = 7  ! Particle "turned around" when not tracking with time_runge_kutta.
-integer, parameter :: lost_z_aperture$ = 9
+integer, parameter :: lost_z_aperture$ = 7
+integer, parameter :: lost_pz_aperture$ = 8  ! Particle "turned around" when not tracking with time_runge_kutta.
+
+! State_name is not the full list of coord%state possible settings! Missing is not_set$
+character(12), parameter ::state_name(0:8) = [character(12):: 'Pre_Born', 'Alive', 'Lost', 'Hit_Neg_X', &
+                                                 'Hit_Pos_X', 'Hit_Neg_Y', 'Hit_Pos_Y', 'Hit_Pz_Aper', 'Hit_Z_Aper']
+
 
 real(rp), parameter :: vec0$(6) = 0
 
@@ -2414,39 +2419,59 @@ end function next_in_branch
 !-------------------------------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------------------------------
 !+
-! Function coord_state_name (coord_state) result (state_str)
+! Function coord_state_name (coord_state, one_word) result (state_str)
 !
 ! Routine to return the string representation of a coord%state state.
 !
 ! Input:
 !   coord_state -- integer: coord%state value
+!   one_word    -- logical, optional. Default False. If True then output string will be one word (no blanks).
+!   
 !
 ! Output:
 !   state_str   -- character(16): String representation.
 !-
 
-function coord_state_name (coord_state) result (state_str)
+function coord_state_name (coord_state, one_word) result (state_str)
 
 implicit none
 
 integer coord_state
 character(16) state_str
+logical, optional :: one_word
 
 !
 
-select case (coord_state)
-case (pre_born$);              state_str = 'Pre_Born'
-case (alive$);                 state_str = 'Alive'
-case (lost$);                  state_str = 'Lost'
-case (not_set$);               state_str = 'Not_Set'
-case (lost_neg_x_aperture$);   state_str = 'Hit -X Side'
-case (lost_pos_x_aperture$);   state_str = 'Hit +X Side'
-case (lost_neg_y_aperture$);   state_str = 'Hit -Y Side'
-case (lost_pos_y_aperture$);   state_str = 'Hit +Y Side'
-case (lost_pz_aperture$);      state_str = 'Hit Energy Aper'
-case (lost_z_aperture$);       state_str = 'Hit Z Side'
-case default;                  state_str = 'UNKNOWN!'
-end select
+if (logic_option(.false., one_word)) then
+  select case (coord_state)
+  case (not_set$);               state_str = 'Not_Set'
+  case (pre_born$);              state_str = 'Pre_Born'
+  case (alive$);                 state_str = 'Alive'
+  case (lost$);                  state_str = 'Lost'
+  case (lost_neg_x_aperture$);   state_str = 'Hit_Neg_X'
+  case (lost_pos_x_aperture$);   state_str = 'Hit_Pos_X'
+  case (lost_neg_y_aperture$);   state_str = 'Hit_Neg_Y'
+  case (lost_pos_y_aperture$);   state_str = 'Hit_Pos_Y'
+  case (lost_pz_aperture$);      state_str = 'Hit_Pz_Aper'
+  case (lost_z_aperture$);       state_str = 'Hit_Z_Aper'
+  case default;                  state_str = 'UNKNOWN!'
+  end select
+
+else
+  select case (coord_state)
+  case (not_set$);               state_str = 'Not_Set'
+  case (pre_born$);              state_str = 'Pre_Born'
+  case (alive$);                 state_str = 'Alive'
+  case (lost$);                  state_str = 'Lost'
+  case (lost_neg_x_aperture$);   state_str = 'Hit -X Side'
+  case (lost_pos_x_aperture$);   state_str = 'Hit +X Side'
+  case (lost_neg_y_aperture$);   state_str = 'Hit -Y Side'
+  case (lost_pos_y_aperture$);   state_str = 'Hit +Y Side'
+  case (lost_z_aperture$);       state_str = 'Hit Z Side'
+  case (lost_pz_aperture$);      state_str = 'Hit Energy Aper'
+  case default;                  state_str = 'UNKNOWN!'
+  end select
+endif
 
 end function coord_state_name
 

bmad/multiparticle/beam_file_io.f90

@@ -72,7 +72,7 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat)
   bunch => beam%bunch(ib)
   write (iu, *) 'BEGIN_BUNCH'
   write (iu, *) '  ', trim(species_name(bunch%particle(1)%species))
-  write (iu, *) bunch%charge_tot, '  ! bunch_charge_tot'
+  write (iu, *) bunch%charge_tot, '  ! charge_tot'
   write (iu, *) bunch%z_center,   '  ! z_center'
   write (iu, *) bunch%t_center,   '  ! t_center'
   do ip = 1, size(bunch%particle)
@@ -97,7 +97,7 @@ end subroutine write_beam_file
 ! If non_zero, the following components of beam_init are used to rescale the beam:
 !     %n_bunch
 !     %n_particle
-!     %bunch_charge
+!     %charge_tot
 !
 ! If the beam file has '.h5' or '.hdf5' suffix then the file is taken to be an HDF5 file.
 ! Otherwise the file is assumed to be ASCII.
@@ -128,7 +128,7 @@ subroutine read_beam_file (file_name, beam, beam_init, err_flag, ele, print_mom_
 integer i, j, k, n, np, ix, iu, ix_word, ios, ix_ele, species
 integer n_bunch, n_particle, n_particle_lines, ix_lost
 
-real(rp) vec(6), sum_charge, bunch_charge
+real(rp) vec(6), sum_charge, charge_tot
 complex(rp) spinor(2)
 
 character(*) file_name
@@ -554,7 +554,7 @@ end subroutine read_beam_file
 ! If non_zero, the following components of beam_init are used to rescale the beam:
 !     %n_bunch
 !     %n_particle
-!     %bunch_charge
+!     %charge_tot
 !
 ! If the beam file has '.h5' or '.hdf5' suffix then the file is taken to be an HDF5 file.
 ! Otherwise the file is assumed to be ASCII.
@@ -579,59 +579,133 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin
 type (ele_struct), optional :: ele
 type (bunch_struct), pointer :: bunch
 type (coord_struct), pointer :: p
+type (coord_struct) p0
 
-real(rp) bunch_charge, z_center, t_center
+real(rp) charge_tot, z_center, t_center
 
 integer iu, ip, ix, ios, n_particle, n_bunch
 
 character(*) file_name
 character(*), parameter :: r_name = 'read_beam_ascii4'
-character(200) cols, line
+character(200) cols, line, str
 
-logical err_flag
+logical err_flag, err
 logical, optional :: print_mom_shift_warning, conserve_momentum
 
 !
 
-
-cols = ''
 n_bunch = 0
+err_flag = .true.
+err = .false.
+
+! bunch loop
 do
+  n_bunch = n_bunch + 1
+  call reallocate_beam(beam, n_bunch, save = .true.)
+  call reallocate_bunch(beam%bunch(n_bunch), 1000)
+  bunch => beam%bunch(n_bunch)
+
+  cols = ''
+  p0 = coord_struct()
+
+  ! Read bunch header
   do
     read (iu, '(a)', iostat = ios) line
     if (line(1:1) /= '#') exit
     call string_trim(line(2:), line, ix)
 
     select case (line(:ix))
-    case ('n_bunch');       n_bunch      = nint(read_param(line))
-    case ('n_particle');    n_particle   = nint(read_param(line))
-    case ('bunch_charge');  bunch_charge = read_param(line)
-    case ('z_center');      z_center     = read_param(line)
-    case ('t_center');      t_center     = read_param(line)
-    case ('columns');       cols         = read_string(line)
+    case ('charge_tot');  bunch%charge_tot   = read_param(line)
+    case ('z_center');    bunch%z_center     = read_param(line)
+    case ('t_center');    bunch%t_center     = read_param(line)
+
+    case ('columns');     cols               = read_string(line)
+
+    case ('species')
+      str = read_string(line)
+      p0%species = species_id(str, positron$)
+
+    case ('location')
+      str = read_string(line)
+      call match_word(str, location_name(1:), p0%location)
+      if (p0%location <= 0) then
+        call out_io (s_error$, r_name, 'LOCATION NAME NOT RECOGNIZED: ' // str)
+        return
+      endif
+
+    case ('state')
+      str = read_string(line)
+      call match_word(str, state_name, p0%state)
+      if (p0%state <= 0) then
+        call out_io (s_error$, r_name, 'PARTICLE STATE NAME NOT RECOGNIZED: ' // str)
+        return
+      endif
+      p0%state = p0%state - 1   ! Since state_name is zero based.
+
+    case ('s_position');  p0%s           = read_param(line)
+    case ('time');        p0%t           = read_param(line)
+    case ('p0c');         p0%p0c         = read_param(line)
+    case ('charge');      p0%charge      = read_param(line)
+    case ('time_dir');    p0%time_dir    = nint(read_param(line))
+    case ('direction');   p0%direction   = nint(read_param(line))
+    case ('ix_ele');      p0%ix_ele      = nint(read_param(line))
+    case ('ix_branch');   p0%ix_branch   = nint(read_param(line))
     end select
   enddo
 
-  call reallocate_beam(beam, n_bunch, n_particle)
+  ! Read particle info
+  backspace(unit = iu)
 
-  do ip = 1, n_particle
-    p => beam%bunch(ip)%particle(ip)
-    read (iu, *, iostat = ios, err = 9000) p%vec, p%p0c
-  enddo
+  ip = 0
+  do 
+    read (iu, '(a)', iostat = ios, end = 8000) line
+    if (line == '') cycle
+    if (line(1:1) == '#') exit
+    if (ios /= 0) then
+      call out_io (s_error$, r_name, 'CANNOT READ BEAM FILE: ' // file_name)
+      return
+    endif
 
-  n_bunch = n_bunch + 1
-  !!call reallocate_bunch(
+    ip = ip + 1
+    if (ip > size(bunch%particle)) call reallocate_bunch (bunch, 2*ip, .true.)
 
+    p => bunch%particle(ip)
+    p = p0
 
-enddo
+    do
+      call string_trim(cols, cols, ix)
+      call read_particle_params(p, cols(1:ix), err); if (err) return
+      cols = cols(ix+1:)
+      if (cols == '') exit
+    enddo
+  enddo
 
+enddo
 
-9000 continue
-stop
+8000 continue
+call reallocate_bunch(bunch, ip, .true.)
+err_flag = .false.
 
 !---------------------------------------------------------------------------------------------------
 contains
 
+subroutine read_particle_params(p, col, err)
+type (coord_struct) p
+character(*) col
+logical err
+
+!
+
+select case (col)
+
+case default
+end select
+
+end subroutine read_particle_params
+
+!---------------------------------------------------------------------------------------------------
+! contains
+
 function read_param(line) result (param)
 character(*) line
 real(rp) param