Foil syntax now accepts total density.

bmad/code/attribute_bookkeeper.f90

@@ -45,9 +45,10 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
 type (converter_prob_pc_r_struct), pointer :: ppcr
 type (molecular_component_struct), allocatable :: component(:)
 type (material_struct), pointer :: material
+type (material_struct) :: materi
 
 real(rp) factor, e_factor, gc, f2, phase, E_tot, polarity, dval(num_ele_attrib$), time, beta
-real(rp) w_inv(3,3), len_old, f, dl, b_max, zmin, ky, kz
+real(rp) w_inv(3,3), len_old, f, dl, b_max, zmin, ky, kz, tot_mass
 real(rp), pointer :: val(:), tt
 real(rp) knl(0:n_pole_maxx), tilt(0:n_pole_maxx), eps6
 real(rp) kick_magnitude, bend_factor, quad_factor, radius0, step_info(7), dz_dl_max_err
@@ -472,33 +473,64 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
   call molecular_components(ele%component_name, component)
   n = size(component)
   if (.not. allocated(ele%foil%material)) allocate(ele%foil%material(n))
+  materi = ele%foil%material(1)
+
+  if (n > 1 .and. size(ele%foil%material) == 1) then
+    deallocate (ele%foil%material)
+    allocate(ele%foil%material(n))
+    ele%foil%material(1) = materi
+  endif
+
   if (n /= size(ele%foil%material)) then
     call out_io(s_error$, r_name, 'NUMBER OF COMPONENTS IN: ' // quote(ele%component_name) // ' (' // int_str(n) // &
-                    ') IS NOT THE SAME AS OTHER PARAMETERS.')
+                    ') IS NOT THE SAME AS OTHER PARAMETERS IN ' // ele%name)
     if (global_com%exit_on_error) call err_exit
     return
   endif
 
+  tot_mass = 0
   do ix = 1, n
     material => ele%foil%material(ix)
     material%species = species_id(component(ix)%atom)
+    material%number = component(ix)%number
+    tot_mass = tot_mass + mass_of(material%species) * material%number
+  enddo
+
+  if (n > 1) then
+    if (materi%density /= real_garbage$ .and. ele%foil%material(2)%density == real_garbage$) then
+      do ix = 1, n
+        material => ele%foil%material(ix)
+        material%density = materi%density * mass_of(material%species) * material%number / tot_mass
+      enddo  
+    endif
+
+    if (materi%area_density /= real_garbage$ .and. ele%foil%material(2)%area_density == real_garbage$) then
+      do ix = 1, n
+        material => ele%foil%material(ix)
+        material%area_density = materi%area_density * mass_of(material%species) * material%number / tot_mass
+      enddo  
+    endif
+  endif
+
+  do ix = 1, n
+    material => ele%foil%material(ix)
     z_material = atomic_number(material%species)
 
-    if (material%radiation_length == 0) then
+    if (material%radiation_length == real_garbage$) then
       material%radiation_length_used = x0_radiation_length(material%species)
     else
       material%radiation_length_used = material%radiation_length
     endif
 
-    if (material%density /= 0 .and. material%area_density /= 0) then
+    if (material%density /= real_garbage$ .and. material%area_density /= real_garbage$) then
       call out_io(s_error$, r_name, 'SETTING BOTH DENSITY AND AREA_DENSITY IS NOT PERMITTED FOR: ' // ele%name)
       return
     endif
 
-    if (material%density == 0) then
-      if (material%area_density /= 0) then
+    if (material%density == real_garbage$) then
+      if (material%area_density /= real_garbage$) then
         if (ele%value(thickness$) == 0) then
-          material%density_used = 0
+          material%density_used = real_garbage$
         else
           material%density_used = material%area_density / ele%value(thickness$)
         endif
@@ -509,13 +541,13 @@ subroutine attribute_bookkeeper (ele, force_bookkeeping)
       material%density_used = material%density
     endif
 
-    if (material%area_density == 0) then
+    if (material%area_density == real_garbage$) then
       material%area_density_used = material%density_used * ele%value(thickness$)
     else
       material%area_density_used = material%area_density
     endif
 
-    if (material%density == 0 .and. material%area_density == 0 .and. n > 1) then
+    if (material%density == real_garbage$ .and. material%area_density == real_garbage$ .and. n > 1) then
       call out_io(s_warn$, r_name, 'Neither foil DENSITY(s) nor AREA_DENSITY(s) set for compound material for: ' // ele%name, &
                                    'This will produce HIGHLY inaccurate results!')
     endif

bmad/doc/elements.tex

@@ -2290,10 +2290,16 @@ \section{Foil}
 \begin{example}
   f1: foil, thickness = 1e-2, material_type = "B4C", 
           density = (2e3, 1e3), radiation_length = (5.49, 4.26)
+  f1: foil, thickness = 1e-2, material_type = "B4C", 
+          density = 9e3, radiation_length = (5.49, 4.26)  ! Same as above
 \end{example}
 Here, since \vn{material_type} is set to \vn{B4C}, there are two components: Boron and Carbon. If
-\vn{density}, \vn{area_density}, or \vn{radiation_length} are present, they must have the same
-number of values as material components. Values are set in order so, in the above example, the
+\vn{radiation_length} is used, it must have the same number of values as the number of material
+components. If \vn{density}, or \vn{area_density} is used they either must have the same number of
+values as the number of material components or they can have a single value.  If there is a single
+value, this value represents the sum of the individual component values.
+
+Values are set in order so, in the above example, the
 Carbon component has a density of 1e3 and a radiation length of 4.26.\footnote
   {
 From a computational standpoint it does not matter which parameter is associated with which component.

bmad/modules/bmad_struct.f90

@@ -1280,9 +1280,10 @@ module bmad_struct
 
 type material_struct
   integer :: species = not_set$
-  real(rp) :: density = 0, density_used = 0
-  real(rp) :: area_density = 0, area_density_used = 0
-  real(rp) :: radiation_length = 0, radiation_length_used = 0
+  integer :: number = int_garbage$                      ! Relative number
+  real(rp) :: density = real_garbage$, density_used = real_garbage$
+  real(rp) :: area_density = real_garbage$, area_density_used = real_garbage$
+  real(rp) :: radiation_length = real_garbage$, radiation_length_used = real_garbage$
 end type
 
 type foil_struct

bmad/output/type_ele.f90

@@ -430,10 +430,11 @@ subroutine type_ele (ele, type_zero_attrib, type_mat6, type_taylor, twiss_out, t
                           attribute_units('DENSITY'), attribute_units('AREA_DENSITY'), attribute_units('RADIATION_LENGTH')
       endif
 
-      nl=nl+1; write(li(nl), '(3(a, es14.6))') '  Density      =', matter%density, &
-                  '  Area_Density      =', matter%area_density,      '  Radiation_Length      =', matter%radiation_length
-      nl=nl+1; write(li(nl), '(3(a, es14.6))') '  Density_Used =', matter%density_used, &
-                  '  Area_Density_Used =', matter%area_density_used, '  Radiation_Length_Used =', matter%radiation_length_used
+      nl=nl+1; write(li(nl), '(3(a, a14))') '  Density      =', this_real(matter%density, 'es14.6', '  Not_Set'), &
+             '  Area_Density      =', this_real(matter%area_density, 'es14.6', '  Not_Set'), &
+             '  Radiation_Length      =', this_real(matter%radiation_length, 'es14.6', '  Not_Set')
+      nl=nl+1; write(li(nl), '((a, a14), 2(a, es14.6))') '  Density_Used =', this_real(matter%density_used, 'es14.6', '  Not_Used'), &
+             '  Area_Density_Used =', matter%area_density_used, '  Radiation_Length_Used =', matter%radiation_length_used
     enddo
 endif
 
@@ -1811,4 +1812,24 @@ subroutine print_this_stack(nl, li, stack, print_it, str_index)
 
 end subroutine print_this_stack
 
+!--------------------------------------------------------------------------
+! contains
+
+function this_real(value, fmt, garbage_str) result (out_str)
+
+real(rp) value
+character(*) fmt, garbage_str
+character(30) out_str
+
+!
+
+if (value == real_garbage$) then
+  out_str = garbage_str
+  return
+endif
+
+write (out_str, '(' // fmt // ')') value
+
+end function this_real
+
 end subroutine type_ele

bmad/output/write_bmad_lattice_file.f90

@@ -526,7 +526,7 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
 
     if (associated(ele%foil)) then
       if (size(ele%foil%material) > 1) then
-        if (any(ele%foil%material%density /= 0)) then
+        if (any(ele%foil%material%density /= real_garbage$)) then
           line = trim(line) // ', density = ('
           do n = 1, size(ele%foil%material)
             if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%density)
@@ -536,7 +536,7 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
           line = trim(line) // ')'
         endif
 
-        if (any(ele%foil%material%area_density /= 0)) then
+        if (any(ele%foil%material%area_density /= real_garbage$)) then
           line = trim(line) // ', area_density = ('
           do n = 1, size(ele%foil%material)
             if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%area_density)
@@ -546,7 +546,7 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
           line = trim(line) // ')'
         endif
 
-        if (any(ele%foil%material%radiation_length /= 0)) then
+        if (any(ele%foil%material%radiation_length /= real_garbage$)) then
           line = trim(line) // ', radiation_length = ('
           do n = 1, size(ele%foil%material)
             if (n == 1) then; line = trim(line) // re_str(ele%foil%material(n)%radiation_length)
@@ -558,9 +558,9 @@ subroutine write_bmad_lattice_file (bmad_file, lat, err, output_form, orbit0)
 
       else
         material => ele%foil%material(1)
-        if (material%density /= 0)          line = trim(line) // ', density = ' // re_str(material%density)
-        if (material%area_density /= 0)     line = trim(line) // ', area_density = ' // re_str(material%area_density)
-        if (material%radiation_length /= 0) line = trim(line) // ', radiation_length = ' // re_str(material%radiation_length)
+        if (material%density /= real_garbage$)          line = trim(line) // ', density = ' // re_str(material%density)
+        if (material%area_density /= real_garbage$)     line = trim(line) // ', area_density = ' // re_str(material%area_density)
+        if (material%radiation_length /= real_garbage$) line = trim(line) // ', radiation_length = ' // re_str(material%radiation_length)
       endif
     endif
 

bmad/parsing/parser_set_attribute.f90

@@ -58,6 +58,7 @@ subroutine parser_set_attribute (how, ele, delim, delim_found, err_flag, pele, c
 type (photon_element_struct), pointer :: ph
 type (photon_reflect_table_struct), allocatable :: rt_save(:)
 type (photon_reflect_table_struct), pointer :: rt
+type (material_struct) :: material
 
 real(rp) kx, ky, kz, tol, value, coef, r_vec(10), r0(2), vec(1000)
 real(rp), allocatable :: table(:,:), arr(:)
@@ -1783,18 +1784,26 @@ subroutine parser_set_attribute (how, ele, delim, delim_found, err_flag, pele, c
   if (.not. ok) return
 
   if (allocated(ele%foil%material)) then
-    if (size(ele%foil%material) /= n) then
+    if (size(ele%foil%material) == 1 .and. n > 1) then
+      material = ele%foil%material(1)
+      deallocate (ele%foil%material)
+      allocate(ele%foil%material(n))
+      ele%foil%material(1) = material
+    endif
+
+    if (size(ele%foil%material) /= n .and. (attrib_word == 'RADIATION_LENGTH' .or. (n > 1 .and. size(ele%foil%material) > 1))) then
       call parser_error('MATERIAL_TYPE, DENSITY, AREA_DENSITY, AND RADIATION_LENGTH MUST ALL BE THE SAME SIZE VECTORS FOR ELE: ' // ele%name)
       return
     endif
+
   else
     allocate(ele%foil%material(n))
   endif
 
   select case (attrib_word)
-  case ('DENSITY');           ele%foil%material(:)%density = arr
-  case ('AREA_DENSITY');      ele%foil%material(:)%area_density = arr
-  case ('RADIATION_LENGTH');  ele%foil%material(:)%radiation_length = arr
+  case ('DENSITY');           ele%foil%material(1:n)%density = arr(1:n)
+  case ('AREA_DENSITY');      ele%foil%material(1:n)%area_density = arr(1:n)
+  case ('RADIATION_LENGTH');  ele%foil%material(1:n)%radiation_length = arr(1:n)
   end select
 
   if (delim == ')') then