Interestingly, the symbolic evaluation mentioned in the previous post now “just works”. Not sure what I did wrong. Here’s the updated snippet.

focal_lengths = np.linspace(length_of_recycling_cavity - 10.0, length_of_recycling_cavity, 11)
trace = virgo.propagate_beam(from_node=virgo.mIM.p2.o, to_node=virgo.mPR.p1, symbolic=True)
phases = trace.total_acc_gouy.eval(subs={virgo.mPR2.f: focal_lengths})

for fl, phase in zip(focal_lengths, phases):
    print(f'f={fl}m: {phase:.2f}deg')

Also, instead of the additional bp detectors, one can get to all the interesting bits directly from the trace output, e.g.

w @ PR = {trace.w(virgo.mPR.p1.i)*1000.0:.3f}mm
w0 = {trace.w0(virgo.mPR.p1.i)*1000.0:.3f}mm
wavefront radius @ PR = {trace.q(virgo.mPR.p1.i).Rc:.3f}m