Overview

Teaching: 0 min
Exercises: 90 min

Questions

• How can I perform a shielding calculation using Python?

Objectives

• Produce a Python script from scratch

• Execute Python scripts in a command line

• Select appropriate data structures for an application

• Write a function in Python

• Apply best practices when writing a function

• Implement for loops to solve a problem

Shielding Problem

In this lesson we’ll be doing a live coding exercise in which we calculate the uncollided flux on the other side of a heterogeneous shield composed of lead, water, and tungsten. We’ll be performing this calculation as a 1-D problem using the following model:

$$\huge \Phi = \Phi_{o} e^{-\mu x}$$

where

• $\Phi_{o}$ - is the initial flux at the front of the shield (units: $cm^{-2}s^{-1}$)
• $\mu$ - is the attenuation coefficient of the shield (units: $cm^{-1}$)
• $x$ - is the shield thickness (units: $cm$)
• $\Phi$ - is the uncollided flux at the back of the shield (units: $cm^{-2}s^{-1}$)

To begin, with a simple script which defines these variables and performs the calculation for a shield with one material, lead, with a thickness of 15 cm.

Checkpoint 0

        
"""
Checkpoint 0: Initial Problem
------------------------------

Compute flux after 5 cm of lead shielding using simple linear
attenuation formula.
"""

# need to import the "math" module to use different math functions
import math as m

# set values of initial flux (phi_0), attenuation coefficient (mu),
# and shield thickness (x)
phi_0 = 1.e15   # photon/cm**2-s
mu = 0.5182     # attenuation coefficient 1/cm
x = 5.0         # shield thickness in cm

# calcuate phi
phi = phi_0 * m.exp(-mu*x)

# print the final result
print(phi)
 
        
      

To get the version of the script shown above, use this link

The beginning of this script is a commented section of text describing what the script does. Because this code is encased in triple-quoted (""") both above and below the text, Python will ignore this section and not treat it as code. This is a nice way to put multi-line comments into your code instead of starting each new line with a # symbol.

Next Python’s built-in math module is imported into the script. Similar to the built-in functions we discussed in episode 2, Python comes with a set of standard modules which add commonly needed capabilities to your programs. In this case, we’ll use the math module to access the exponential function (math.exp). Notice, though, that when this function is called in later lines it is called as m.exp. This is because we’ve imported the math module as m. This is known as aliasing and is often used to shorten long module names for convenience.

Aliasing

Aliasing is used to import a module under a different name, or alias. This is typically used to shorten the name of a module. The module will behave the the same way under any name. For example:

import math
math.cos(0.0)

will yield the same result as

import math as m
m.cos(0.0)

Next we initialize a few variables: phi_0, mu, and x. These represent the initial flux, attenuation coefficient of lead, and the lead shield thickness respectiely. Finally, we perform the calculation for the uncollided flux according to the model above and display the result on screen.

Two Shields

Now let’s update the script so that it reports the uncollided flux for another shield composed of water. Perhaps the most straightforward way to do this is to retype our code, but update the attenuation and thickness values to represent a water shield.

Checkpoint 1

        
"""
Checkpoint 1: Two different problems
------------------------------------

Compute two separate attenuation problems for lead and plastic using
different initial fluxes, thickness, and attenuation coefficients.
"""

import math as m

# Compute result after lead shield 5 cm thick with an initial flux of
# 1e10 photons/cm**2-s
phi_0 = 1.e10   # photons/cm**2-s
mu = 0.5182     # attenuation coefficient 1/cm
x = 5.0         # shield thickness in cm
phi_lead = phi_0 * m.exp(-mu*x)

print(phi_lead)


# Compute result after water shield 15 cm thick with an initial flux
# of 1e15 photons/cm**2-s
phi_0 = 1.e15   # photons/cm**2-s
mu = 0.0493     # attenuation coefficient 1/cm
x = 15.0        # shield thickness in cm
phi_plastic = phi_0 * m.exp(-mu*x)

print(phi_plastic)
 
        
      

To get the version of the script shown above, use this link

This is fairly straightforward, but what if we want to perform this computation 10 times? 50 times? 10000 times? This small script would become large very quickly. A good way to provide the ability to repeat this calculation is to create a function for calculating the flux given an attenuation coefficient and material thickness.

Creating a function

In the code below, a function returning the flux for a shield (a $\mu$ and $x$ value) has been implemented and documented. Underneath the line defining the calc_phi function there is a set of text enclosed with triple-quotes (""") which contains information about the purpose of the function along with its expected inputs and outputs. Triple quotes were used before to create multi-line comments at the top of our script, but this comment block is special because of its location in the code. Because this comment section appears right below the function definition, it acts as documentation and is accessible from a Python interpreter using the built-in help function. Specifically, executing the command help(calc_phi) will display this information about the function on screen.

Underneath this is the implementation of the function which computes the flux value.

Checkpoint 2

        
"""
Checkpoint 2: Modularize a repeated problem
-------------------------------------------

Create a function to calculate the result for any initial flux, shield
thickness, and attenuation coefficient.
"""

import math as m


def calc_phi(phi_0, mu, x):
    """Calculates final flux after shield of thickness x with
    attenuation coefficient mu given the starting flux phi_0.

    Input:
    ------
        phi_0: float, initial flux value (photons/cm**2-s)
        mu: float, attenuation coefficient (1/cm)
        x: float, shield thickness (cm)

    Returns:
    --------
        phi: float, final flux value (photons/cm**2-s)
    """

    phi = phi_0 * m.exp(-mu*x)
    return(phi)


# Lead shield
phi_0 = 1.e10   # photons/cm**2-s
mu = 0.5182     # attenuation coefficient 1/cm
x = 5.0         # shield thickness in cm
phi_lead = calc_phi(phi_0, mu, x)
print(phi_lead)


# Water shield
phi_0 = 1.e15   # photons/cm**2-s
mu = 0.0493     # attenuation coefficient 1/cm
x = 15.0        # shield thickness in cm
phi_water = calc_phi(phi_0, mu, x)
print(phi_water)
 
        
      

To get the version of the script shown above, use this link

Correlating Shield Properties

For the purposes of today’s exercise, the two pieces of information required to define a shield are its attenuation coefficient, $\mu$, and its thickness, $x$. Rather than separate these values in to different variable, a tuple can be used to make sure these two properties will always be associated with each other. Our function for this calculation now needs to be updated as well.

Checkpoint 3

        
"""
Checkpoint 3: Passing a set of properties
-----------------------------------------

Pass a tuple representing the shielding properties (x and mu) to the
function.
"""

import math as m

def calc_phi(phi_0, mat_props):
    """Calculates final flux after shield of thickness x with
    attenuation coefficient mu given the starting flux phi_0.

    Input:
    ------
        phi_0: float, initial flux value (photons/cm**2-s)
        mat_props: tuple, a tuple of size 2 in the form (mu, x) where
            mu and x are:
                mu: float, attenuation coefficient (1/cm)
                x: float, shield thickness (cm)

    Returns:
    --------
        phi: float, final flux value (photons/cm**2-s)
    """

    # unpack properties
    mu = mat_props[0]
    x = mat_props[1]

    # calcluate phi
    phi = phi_0 * m.exp(-mu*x)

    return(phi)


# Lead shield
lead_props = (0.5182, 5.0)
phi_lead = calc_phi(1.e10, lead_props)
print(phi_lead)


# Water shield
water_props = (0.0493, 15.0)
phi_water = calc_phi(1.e15, water_props)
print(phi_water)
 
        
      

To get the version of the script shown above, use this link

As seen in the updated function documentation, the function calc_phi now takes in only two arguments - the initial flux phi_0 and an argument named mat_props. Now when this argument is passed into the function we can unpack the attenuation coefficient and shield thickness from mat_props variable and perform our calculation. Note that if these lines were not being added in a function, they would be added in the main body of the script for each shield instead.

A Two-Material Shield

Now that we have a function in place which returns the uncollided flux for a single shield. We can explore that calculation for a multi-layer shield in the code below. Because the material information is passed in as a list, it makes the extension of this function for attenuation through multiple materials fairly natural. The calculation is performed for the first material and the resulting flux value is passed into a calculation for the second material. Finally the flux value is returned from the function and printed to screen in at the end of the script. Note that this multi-material calculation now requires only one call to the calc_phi function.

Checkpoint 4

        
"""
Checkpoint 4: Shield made of two materials
------------------------------------------

Pass a set of material properties representing two materials in a shield
each with their own thickness and attenuation coefficient.
"""

import math as m

def calc_phi(phi_0, mat_props):
    """Calculates final flux after two shields of different thicknesses
    and attenuation coefficients given the starting flux phi_0.

    Input:
    ------
        phi_0: float, initial flux value (photons/cm**2-s)
        mat_props: list of tuples, each tuple represents a set of
            properties for each material in the problem in the form
            [(mu0, x0), (mu1, x1)] where:
                mu0: float, attenuation coefficient for first shield
                    (1/cm)
                x0: float, shield thickness for first shield (cm)
                mu1: float, attenuation coefficient for second shield
                    (1/cm)
                x1: float, shield thickness for second sheild (cm)

    Returns:
    --------
        phi_2: float, final flux value after second shield
            (photons/cm**2-s)
    """

    # unpack each set of material properties
    mat0 = mat_props[0]     # first set of properties (mu0, x0)
    mat1 = mat_props[1]     # second set of properties (mu1, x1)

    # first material properties and get flux after first shield
    mu0 = mat0[0]
    x0 = mat0[1]
    phi_1 = phi_0 * m.exp(-mu0*x0)

    # get second material properties and get flux after second shield
    mu1 = mat1[0]
    x1 = mat1[1]
    phi_2 = phi_1 * m.exp(-mu1*x1)

    # return final flux value
    return(phi_2)


# Set initial flux and material properties
phi_0 = 1.e15
lead_props = (0.5182, 5.0)
water_props = (0.0493, 15.0)
mat_props = [lead_props, water_props]

phi_final = calc_phi(phi_0, mat_props)

print(phi_final)

 
        
      

To get the version of the script shown above, use this link

A Multi-Material Shield

Our previous work made possible to calculate the flux after passing through two materials, but what about 5 materials? 10 materials?

The design of the calc_phi function makes this possible with a relatively small change to the code. Rather than assuming a specific number of materials are being passed into the function, we can use a for loop to calculate the flux after each material and return the final flux value at the end of the function.

Checkpoint 5

        
"""
Checkpoint 5: Shield made of arbitrary number of materials
----------------------------------------------------------

Pass a set of material properties representing any number of materials
in a shield each with their own thickness and attenuation coefficient.
"""

import math as m

def calc_phi(phi, mat_props):
    """Calculates final flux after two shields of different thicknesses
    and attenuation coefficients given the starting flux phi_0.

    Input:
    ------
        phi: float, initial flux value (photons/cm**2-s)
        mat_props: list of tuples, each tuple represents a set of
            properties for each material in the problem in the form
            [(mu0, x0), (mu1, x1), (mu2, x2) ... (muN, xN)] for any N
            number of shields. Each shielding material (one tuple)
            should be defined by:
                mu: float, attenuation coefficient (1/cm)
                x: float, shield thickness (cm)

    Returns:
    --------
        phi: float, final flux value after all shields (photons/cm**2-s)
    """

    # iterate over each shield in the property list
    for mat in mat_props:

        # unpack current shield properties mu and x
        mu = mat[0]
        x = mat[1]

        # update phi_0 with new flux value after attenuation through
        # current shield
        phi = phi * m.exp(-mu*x)

    # return final flux value
    return(phi)


# Set properties for various shielding types:
lead_props = (0.5182, 5.0)
water_props = (0.0493, 15.0)
tungsten_props = (0.0437, 7.0)
aluminum_props = (0.1166, 10.0)

# Set initial flux
phi_0 = 1.e15

# Calculate using three shields (lead, water, and tungsten)
mat_props = [lead_props, water_props, tungsten_props]
phi_final = calc_phi(phi_0, mat_props)
print(phi_final)

# Calculate using two shields (water and tungsten)
mat_props = [water_props, tungsten_props]
phi_final = calc_phi(phi_0, mat_props)
print(phi_final)

# Calculate using one shield (aluminum)
mat_props = [aluminum_props]
phi_final = calc_phi(phi_0, mat_props)
print(phi_final)
 
        
      

To get the version of the script shown above, use this link

As demonstrated at the bottom of our code, this function can now be called to calculate the uncollided flux at the end of an arbitrary set of shields with varying thicknesses. It is also easy to set up different materials for the calculation and use them interchangeably.

Key Points

• Use Python to write a simple shielding calculation script.

• Learn some best practices when writing Python.