diff --git a/bmad/multiparticle/beam_file_io.f90 b/bmad/multiparticle/beam_file_io.f90 index 955170de79..139945c39f 100644 --- a/bmad/multiparticle/beam_file_io.f90 +++ b/bmad/multiparticle/beam_file_io.f90 @@ -621,26 +621,20 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin 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 + call read_switch(line(:ix), p0%location, str, err) + if (err) return 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. + call read_switch(line(:ix), p0%state, str, err) + if (err) return + + case ('species') + str = read_string(line) + call read_switch(line(:ix), p0%species, str, err) + if (err) return case ('s_position'); p0%s = read_param(line) case ('time'); p0%t = read_param(line) @@ -674,7 +668,7 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin do call string_trim(cols, cols, ix) - call read_particle_params(p, cols(1:ix), err); if (err) return + call read_particle_params(p, cols(1:ix), line, err); if (err) return cols = cols(ix+1:) if (cols == '') exit enddo @@ -689,16 +683,39 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin !--------------------------------------------------------------------------------------------------- contains -subroutine read_particle_params(p, col, err) +subroutine read_particle_params(p, col, line, err) type (coord_struct) p -character(*) col +character(*) col, line logical err ! select case (col) +case ('vec'); call read_components(col, p%vec, line, err); if (err) return +case ('spin'); call read_components(col, p%spin, line, err); if (err) return +case ('field'); call read_components(col, p%field, line, err); if (err) return +case ('s'); call read_component(col, p%s, line, err); if (err) return +case ('t'); call read_component(col, p%t, line, err); if (err) return +case ('phase'); call read_components(col, p%phase, line, err); if (err) return +case ('charge'); call read_component(col, p%charge, line, err); if (err) return +case ('dt_ref'); call read_component(col, p%dt_ref, line, err); if (err) return +case ('r'); call read_component(col, p%r, line, err); if (err) return +case ('p0c'); call read_component(col, p%p0c, line, err); if (err) return +case ('E_potential'); call read_component(col, p%E_potential, line, err); if (err) return +case ('beta'); call read_component(col, p%beta, line, err); if (err) return +case ('ix_ele'); call read_component_int(col, p%ix_ele, line, err); if (err) return +case ('ix_branch'); call read_component_int(col, p%ix_branch, line, err); if (err) return +case ('ix_user'); call read_component_int(col, p%ix_user, line, err); if (err) return +case ('direction'); call read_component_int(col, p%direction, line, err); if (err) return +case ('time_dir'); call read_component_int(col, p%time_dir, line, err); if (err) return +case ('state'); call read_switch(col, p%state, line, err); if (err) return +case ('species'); call read_switch(col, p%species, line, err); if (err) return +case ('location'); call read_switch(col, p%location, line, err); if (err) return case default + err = .true. + call out_io(s_error$, r_name, 'COLUMN NAME NOT RECOGNIZED: ' // col) + return end select end subroutine read_particle_params @@ -706,6 +723,104 @@ end subroutine read_particle_params !--------------------------------------------------------------------------------------------------- ! contains +subroutine read_switch(who, switch, line, err) + +integer switch, ix +character(*) who, line +logical err + +! + +call string_trim(line, line, ix) +if (ix == 0) err = .true. + +select case (who) +case ('location') + call match_word(line(:ix), location_name(1:), switch) + if (switch <= 0) then + call out_io (s_error$, r_name, 'LOCATION NAME NOT RECOGNIZED: ' // line) + return + endif + +case ('state') + call match_word(line(:ix), state_name, switch) + if (switch <= 0) then + call out_io (s_error$, r_name, 'PARTICLE STATE NAME NOT RECOGNIZED: ' // line) + return + endif + switch = switch - 1 ! Since state_name is zero based. + +case ('species') + switch = species_id(line(:ix), positron$) + if (switch == invalid$) err = .true. +end select + +line = line(ix+1:) + +end subroutine read_switch + +!--------------------------------------------------------------------------------------------------- +! contains + +subroutine read_components(who, components, line, err) + +real(rp) components(:) +integer ix +character(*) who, line +logical err + +! + +do ix = 1, size(components) + call read_component(who, components(ix), line, err) + if (err) return +enddo + +end subroutine read_components + +!--------------------------------------------------------------------------------------------------- +! contains + +subroutine read_component(who, component, line, err) + +real(rp) component +integer ix +character(*) who, line +logical err + +! + +call string_trim(line, line, ix) +if (ix == 0) err = .true. +if (err) return +read (line, *, iostat = ios) component +err = (ios /= 0) + +end subroutine read_component + +!--------------------------------------------------------------------------------------------------- +! contains + +subroutine read_component_int(who, component, line, err) + +integer component +integer ix +character(*) who, line +logical err + +! + +call string_trim(line, line, ix) +if (ix == 0) err = .true. +if (err) return +read (line, *, iostat = ios) component +err = (ios /= 0) + +end subroutine read_component_int + +!--------------------------------------------------------------------------------------------------- +! contains + function read_param(line) result (param) character(*) line real(rp) param diff --git a/tao/code/tao_spin_polarization_calc.f90 b/tao/code/tao_spin_polarization_calc.f90 index b92dbfca71..b65c4e1a48 100644 --- a/tao/code/tao_spin_polarization_calc.f90 +++ b/tao/code/tao_spin_polarization_calc.f90 @@ -55,6 +55,11 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k if (present(err_flag)) err_flag = .true. +if (tao_branch%spin%valid .and. tao_branch%spin_map_valid) then + if (present(err_flag)) err_flag = .false. + return +endif + g_tol = 1e-4_rp * branch%param%g1_integral / branch%param%total_length g2_tol = 1e-4_rp * branch%param%g2_integral / branch%param%total_length g3_tol = 1e-4_rp * branch%param%g3_integral / branch%param%total_length @@ -73,8 +78,6 @@ subroutine tao_spin_polarization_calc (branch, tao_branch, excite_zero, ignore_k orbit => tao_branch%orbit tspin => tao_branch%spin if (.not. allocated(tao_branch%spin_ele)) allocate (tao_branch%spin_ele(0:branch%n_ele_track)) -tao_branch%spin_ele(:)%valid = .false. -tspin%valid = .false. if (branch%param%geometry == closed$ .and. orbit(branch%n_ele_track)%state /= alive$) return !