# Doing calculations with Python

As an engineer I do a lot of calculations. These can be done with pen and paper and a calculator, in an IPython notebook or in a throwaway spreadsheet. All of these methods have shortcomings, though.

Pen and paper is hard to share and (in my case) hard for others to read. In IPython you can assign the results of calculations to a variable, but you have to perform a separate action to display them. And spreadsheets in general show you the results but not the calculations.

So I wrote a simple function in Python to help me with that.
Using this function I can print both simple assignments and relatively complex
calculations. And it shows *both* the calculation and the result.

For example, I could write:

```
do('x = 0.55')
do('y = 1 - x')
```

When this code is executed, it produces:

x = 0.55 y = 1 - x = 0.45

## Code

The function that does this is relatively simple.

```
local_dict = {}
def do(expr, comment='', loc=local_dict):
"""Process an assignment expression and print it plus the result.
Arguments:
expr: Assignment expression to print.
comment: Text to be printed after the assignment.
loc: Dictionary of locals. Defaults to local_dict.
"""
before = set(loc.keys())
exec(expr, None, loc)
after = set(loc.keys())
newkeys = after - before
n = len(newkeys)
if n == 0:
raise ValueError('not an assignment')
elif n > 1:
raise ValueError('multiple assignemt')
try:
v = float(expr.split('=')[1])
del v
print('{} {}'.format(expr, comment))
except ValueError:
k = newkeys.pop()
print('{} = {:.4g} {}'.format(expr, loc[k], comment))
```

The call to `exec` is the heart of the implementation. It executes the
statement given to it and updates the `locals` dictionary with variables that
are assigned. Any assigned variables are determined by subtracting the keys of
the local variables dictionary before from those after the `exec` invocation.

The `print` statements perform the necessary output. The rest
is basically housekeeping.

Previously, I’ve written a similar tool for adding calculations to LaTeX documents called texcalc. That is much more complicated. It needs to walk the abstract syntax tree of the expression to properly render it as LaTeX.

I’ve made this available on github.

## Example

Last week someone I cooperated with expressed doubt as to whether my lamprop software was producing correct results.

So I wanted to check the calculations that it does. Below is an implementation of the calculations of the properties of a unidirectional fiber reinforced layer.

```
import math
from calculations import do
print('Properties of the fiber')
do('E1f = 230000', 'MPa, carbon fiber')
do('ν12f = 0.27')
do('α1f = -0.41e-6', 'K⁻¹')
do('ρf = 1.76', 'g/cm³')
print('\nProperties of the resin')
do('Em = 2900', 'MPa, epoxy resin')
do('νm = 0.25')
do('αm = 40e-6', 'K⁻¹')
do('ρm = 1.15', 'g/cm³')
print('\nLayer parameters')
do('w = 550', 'g/m²')
do('vf = 0.75')
do('angle = 0', '°')
print('Engineering properties of the layer')
do('vm = (1 - vf)')
do('tf = w / (ρf * 1000)', 'mm')
do('tl = tf * (1 + vm / vf)', 'mm')
do('E1l = vf * E1f + (1 - vf) * Em', 'MPa')
do('E2l = 3*Em', 'MPa, see Tsai:1992, p. 3-13')
do('G12l = E2l / 2', 'MPa')
do('ν12l = 0.3', 'see Tsai:1992, p. 3-13')
do('ν21l = ν12l * E2l / E1l', 'see Nettles:1994, p. 4')
print('Auxiliary variables')
do('Q11 = E1l / (1 - ν12l * ν21l)', 'MPa')
do('Q12 = ν12l * E2l / (1 - ν12l * ν21l)', 'MPa')
do('Q22 = E2l / (1 - ν12l * ν21l)', 'MPa')
do('Q66 = G12l', 'MPa')
do('m = math.cos(math.radians(angle))')
do('m2 = m * m')
do('m3 = m2 * m')
do('m4 = m2 * m2')
do('n = math.sin(math.radians(angle))')
do('n2 = n * n')
do('n3 = n2 * n')
do('n4 = n2 * n2')
do('QA = Q11 - Q12 - 2 * Q66')
do('QB = Q12 - Q22 + 2 * Q66')
print('Components of the transformed stiffness matrix of the layer (Q̅)')
do('Q̅11 = Q11 * m4 + 2 * (Q12 + 2 * Q66) * n2 * m2 + Q22 * n4', 'MPa')
do('Q̅12 = (Q11 + Q22 - 4 * Q66) * n2 * m2 + Q12 * (n4 + m4)', 'MPa')
do('Q̅16 = QA * n * m3 + QB * n3 * m', 'MPa')
do('Q̅22 = Q11 * n4 + 2 * (Q12 + 2 * Q66) * n2 * m2 + Q22 * m4', 'MPa')
do('Q̅26 = QA * n3 * m + QB * n * m3', 'MPa')
do('Q̅66 = (Q11 + Q22 - 2 * Q12 - 2 * Q66) * n2 * m2 + Q66 * (n4 + m4)', 'MPa')
```

This produced the following result:

Properties of the fiber E1f = 230000 MPa, carbon fiber ν12f = 0.27 α1f = -0.41e-6 K⁻¹ ρf = 1.76 g/cm³ Properties of the resin Em = 2900 MPa, epoxy resin νm = 0.25 αm = 40e-6 K⁻¹ ρm = 1.15 g/cm³ Layer parameters w = 550 g/m² vf = 0.75 angle = 0 ° Engineering properties of the layer vm = (1 - vf) = 0.25 tf = w / (ρf * 1000) = 0.3125 mm tl = tf * (1 + vm / vf) = 0.4167 mm E1l = vf * E1f + (1 - vf) * Em = 1.732e+05 MPa E2l = 3*Em = 8700 MPa, see Tsai:1992, p. 3-13 G12l = E2l / 2 = 4350 MPa ν12l = 0.3 see Tsai:1992, p. 3-13 ν21l = ν12l * E2l / E1l = 0.01507 see Nettles:1994, p. 4 Auxiliary variables Q11 = E1l / (1 - ν12l * ν21l) = 1.74e+05 MPa Q12 = ν12l * E2l / (1 - ν12l * ν21l) = 2622 MPa Q22 = E2l / (1 - ν12l * ν21l) = 8740 MPa Q66 = G12l = 4350 MPa m = math.cos(math.radians(angle)) = 1 m2 = m * m = 1 m3 = m2 * m = 1 m4 = m2 * m2 = 1 n = math.sin(math.radians(angle)) = 0 n2 = n * n = 0 n3 = n2 * n = 0 n4 = n2 * n2 = 0 QA = Q11 - Q12 - 2 * Q66 = 1.627e+05 QB = Q12 - Q22 + 2 * Q66 = 2582 Components of the transformed stiffness matrix of the layer (Q̅) Q̅11 = Q11 * m4 + 2 * (Q12 + 2 * Q66) * n2 * m2 + Q22 * n4 = 1.74e+05 MPa Q̅12 = (Q11 + Q22 - 4 * Q66) * n2 * m2 + Q12 * (n4 + m4) = 2622 MPa Q̅16 = QA * n * m3 + QB * n3 * m = 0 MPa Q̅22 = Q11 * n4 + 2 * (Q12 + 2 * Q66) * n2 * m2 + Q22 * m4 = 8740 MPa Q̅26 = QA * n3 * m + QB * n * m3 = 0 MPa Q̅66 = (Q11 + Q22 - 2 * Q12 - 2 * Q66) * n2 * m2 + Q66 * (n4 + m4) = 4350 MPa

The `comment` parameter of the `do` function can be used to add units to
the results. Also note how I use unicode characters as variable names, thanks
to using Python 3. Unfortunately, using unicode sub- and superscripts are not
allowed in identifiers.