Python Matrices Multiplication with a range of variable

Question

I am an international graduate student who is studying optics recently. I got an HW problem that asks me to solve an optics problem using Python.

Basically, I have to multiply three matrices and plot wavelength (X-Axis) Vs. one component of the final matrix.

So,

I have a certain range of variable (500~600) which will be my X-axis.

There are three matrices I12 L2 I23

All of them are 2x2 matrices.

in L2 Matrice, the variable lam (the x-axis variable) is used.

I think that is why I get an error message. Basically, this(L2 matrix) should be an array of matrices but I don't know how to do this while I12 and I23 matrices are only one without variables that have a range.

[So the final result should be like this photo.][1]

Wish this is still understandable. Please understand that English is not my first language.

I have attached my code below.

import numpy as np
import cmath

lmin, lmax, lstep = 500, 600, 0.1
lam = np.arange(lmin, lmax + lstep, lstep)

d = 5e-6

n1 = 1
n2 = 2.6
n3 = 1.45

r12, t12 = (n1 - n2) / (n1 + n2), (2 * n1) / (n1 + n2)
r23, t23 = (n2 - n3) / (n2 + n3), (2 * n2) / (n2 + n3)

psi = 2 * np.pi * n2 * d / (lam*1e-9)

# Transfer Matrix at the 1st Border (without 1/t12)
I12 = [[1, r12], [r12, 1]]
print(f'I12 Matrix:\n{I12}')

# Transfer Matrix in the Slab
L2 = [[cmath.exp(psi * -1j), 0], [0, cmath.exp(psi * 1j)]] << I got problem here
print(f'L2 Matrix:\n{L2}')

# Transfer Matrix at the 2nd Border (without 1/t23)
I23 = [[1, r23], [r23, 1]]
print(f'I23 Matrix:\n{I23}')

# Total Transfer Matrix: I12 X L2 X I23
coeff = (1 / t12) * (1 / t23)
M1 = np.matmul(L2, I23)

error code

I12 Matrix:
[[1, -0.4444444444444445], [-0.4444444444444445, 1]]
Traceback (most recent call last):
  File "D:\OneDrive - Northeastern University\Desktop\SynologyDrive\Northeastern\courses\EECE7284_Optical Properties of Matter\HW5_3.py", line 24, in <module>
    L2 = [[cmath.exp(psi * -1j), 0], [0, cmath.exp(psi * 1j)]]
TypeError: only length-1 arrays can be converted to Python scalars


  [1]: https://i.sstatic.net/Vp7Cn.png

Answer 1

I made

L2 = [[np.exp(psi * -1j), 0], [0, np.exp(psi * 1j)]] 

Like this and it worked

Answer 2

The problem is psi is an array itself, so you're trying to create a "3d" L2 matrix. Of course, numpy does not know what are you trying to achieve and raises an error. There are (at least) two ways to solve your problem:

Calculate the final matrix symbolically and substitute numerical values in the end. This can be achieved using sympy:

import sympy as sm
from sympy import symbols

d, lam, n1, n2, n3 = symbols('d, la, n1, n2, n3', real=True, positive=True)
r12, t12 = (n1 - n2) / (n1 + n2), (2 * n1) / (n1 + n2)
r23, t23 = (n2 - n3) / (n2 + n3), (2 * n2) / (n2 + n3)
psi = 2 * sm.pi * n2 * d / (lam)
I12 = sm.Matrix([[1, r12], [r12, 1]])
L2 = sm.Matrix([[sm.exp(psi * -sm.I), 0], [0, sm.exp(psi * sm.I)]])
I23 = sm.Matrix([[1, r23], [r23, 1]])
M1 = L2*I23
print(M1)

gives

Matrix([[exp(-2*I*pi*d*n2/la), (n2 - n3)*exp(-2*I*pi*d*n2/la)/(n2 + n3)], 
        [(n2 - n3)*exp(2*I*pi*d*n2/la)/(n2 + n3), exp(2*I*pi*d*n2/la)]])

now just use the formula you need.

Another approach is to iterate over all lam values:

import numpy as np

lmin, lmax, lstep = 500, 600, 0.1
lam = np.arange(lmin, lmax + lstep, lstep)

d = 5e-6

n1 = 1
n2 = 2.6
n3 = 1.45

r12, t12 = (n1 - n2) / (n1 + n2), (2 * n1) / (n1 + n2)
r23, t23 = (n2 - n3) / (n2 + n3), (2 * n2) / (n2 + n3)

result = np.zeros(len(lam),dtype=complex)
for ii in range(len(lam)):
    wavelength = lam[ii]
    psi = 2 * np.pi * n2 * d / (wavelength*1e-9)
    
    # Transfer Matrix at the 1st Border (without 1/t12)
    I12 = np.array([[1, r12], [r12, 1]])
    # print(f'I12 Matrix:\n{I12}')
    
    # Transfer Matrix in the Slab
    L2 = np.array([[np.exp(psi * -1j), 0], [0, np.exp(psi * 1j)]])
    # print(f'L2 Matrix:\n{L2}')
    
    # Transfer Matrix at the 2nd Border (without 1/t23)
    I23 = np.array([[1, r23], [r23, 1]])
    # print(f'I23 Matrix:\n{I23}')
    
    # Total Transfer Matrix: I12 X L2 X I23
    # coeff = (1 / t12) * (1 / t23)
    M1 = np.matmul(L2, I23)

    # get the value you need from M1
    result[ii] = M1[0][1]

the result variable now contain calculated values over the range of lam. There are more efficient and/or elegant ways to iterate over without for using numpy, but I think there is no point to use them in this case.