Calculating the nonlinear reflection coefficient of a crystalline silicon slab

The purpose of this script is to calculate the nonlinear reflection coefficient of a crystalline silicon slab.

It takes some input files (columns of data separated by whitespace), converts that data to NumPy matrices, then operates on that data via some formulas that are coded into different functions. Finally, it spits out the final results to a file, as columns separated by whitespace.

I inherited an old FORTRAN program that serves the same purpose. Even though it's completely illegible it is fast as a bat out of hell! It executes in ~0.020 seconds, while my script clocks in at around ~1.25 seconds.

Please offer feedback for improving the general quality overall. I keep this code on GitHub so any and all feedback will be very helpful for me and others.

```
"""
nrc.py is a python program designed to calculate the Nonlinear reflection
coefficient for silicon surfaces. It works in conjunction with the matrix
elements calculated using ABINIT, and open source ab initio software,
and TINIBA, our in-house optical calculation software.

The work codified in this software can be found in Phys.Rev.B66, 195329(2002).
"""

from math import sin, cos, radians
from scipy import constants, interpolate
from numpy import loadtxt, savetxt, column_stack, absolute, \
sqrt, linspace, ones, complex128

########### user input ###########
OUT = "data/nrc/"
CHI1 = "data/res/chi1"
ZZZ = "data/res/zzz"
ZXX = "data/res/zxx"
XXZ = "data/res/xxz"
XXX = "data/res/xxx"
# Angles
THETA_RAD = radians(65)
PHI_RAD = radians(30)
# Misc
ELEC_DENS = 1e-28 # electronic density and scaling factor (1e-7 * 1e-21)
ENERGIES = linspace(0.01, 12.00, 1200)

########### functions ###########
def nonlinear_reflection():
""" calls the different math functions and returns matrix,
which is written to file """
onee = linspace(0.01, 12.00, 1200)
twoe = 2 * onee
polarization = [["p", "p"], ["p", "s"], ["s", "p"], ["s", "s"]]
for state in polarization:
nrc = rif_consta

Avoid computing the same thing more than once. For example:

• Define global constants such as SIN_THETA_RAD_SQUARED = sin(THETA_RAD) ** 2

• Load the data files just once into global constants, or into local variables in nonlinear_reflection and pass them to reflection_components as parameters.

• Assign epsilon(energy) to a local variable in reflection_components and use that in the lengthy expression.

StackExchange Code Review Q#30109, answer score: 2