This is wonderful simulation and it offers many possibilities. I tried to use it for a particular geometry and I think I have detected a small bug. Or it might be an opportunity for me to learn something new!

I am attaching an image of the setting I am using. This is the output of the lines generated:

j=0,I=1.0,c=1.7669502000928863,x=-75.51111260566398,y=-21.117714323462188

j=0,I=1.0,c=0.26179938779914935,x=-75.51111260566398,y=-21.117714323462188

j=1,I=0.057796105403213144,c=-1.8325957145940464,x=-15.075862213492771,y=-4.924137786507226

j=2,I=0.9422038945967869,c=0.44556125394332635,x=-15.075862213492771,y=-4.924137786507226

j=3,I=0.9422038945967869,c=-0.4455612539433263,x=58.049565975939494,y=30.0

j=4,I=0.0733637418622893,c=-2.6960313996464667,x=100.0,y=9.964797822930265

j=5,I=0.8688401527344976,c=-0.18820469334913387,x=100.0,y=9.964797822930265

j=6,I=0.0733637418622893,c=2.6960313996464667,x=16.320254742697305,y=-29.999999999999993

j=7,I=0.02264313020441028,c=2.016357580738223,x=4.221937396573228,y=-24.221937396573228

j=8,I=0.05072061165787902,c=2.1447458012139355,x=4.221937396573228,y=-24.221937396573228

j=9,I=0.00118179551253197,c=-2.016357580738223,x=-21.67403968794734,y=29.999999999999996

j=10,I=0.02146133469187831,c=2.1705020069039946,x=-21.67403968794734,y=29.999999999999996

j=11,I=0.006915452246752739,c=-2.1447458012139355,x=-30.83551657260938,y=30.0

j=12,I=0.04380515941112628,c=1.7669502000928863,x=-30.83551657260938,y=30.0

The first bug seems to be related to the ray corresponding to j=5. Since this represents a beam going from n=1.5 to n=1, the refracted beam should be larger than what what is shown. The same problem happens to beam j=8 when it refracts into beam j=12. The problem does not happen for beam j=7!

The other item I don't understand is the fact that we have two beams reflected at x=4.221937396573228,y=-24.221937396573228 (j=7 and j=8)

Thanks for your help clarifying this and for all of your amazing work.