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.