When doing optical modal modelling you will at some point have to consider how some incident mode gets scattered by some deformation or effect into another mode. There is a lot of math to take in before you start calculating stuff but there is a good review of it all here:

https://link.springer.com/article/10.1007/s41114-016-0002-8#Sec118

The computational problem is then computing a 2D overlap integral between your field and some output mode of some new basis. Some of the simple ones are doable by hand but often it’s easier to do it via numerical integration. Here is a quick code snippet for Newton-Cotes integrating using Pykat:

import pykat
from pykat.optics.gaussian_beams import HG_mode
from pykat.math import newton_weights

q = pykat.BeamParam(w0=0.5, z=0)
x = np.linspace(-4, 4, 2000)
y = np.linspace(-4, 4, 2000)
dx = x[1] - x[0]
dy = y[1] - y[0]
hg00 = HG_mode(q, n=0, m=0)
hg10 = HG_mode(q, n=1, m=0)
Y, X = np.meshgrid(x,y)

W_nc = np.outer(newton_weights(x, 5), 
                newton_weights(y, 5))
                    
dx*dy*np.sum(W_nc * hg00.Unm(x,y).conj() * X * hg10.Unm(x,y))

This is solving the overlap integral between two modes with some additional effect, X. The integrand is weighted by the Newton-Cotes weights. You can easily add in whatever distortion you need here. Currently this code is computing the overlap with some linear X term, which is required for computing how yaw rotations affect a beam.