Why doesn't my loop work? (2D ray tracer)(found fix)

Questions
, , ,
1

Hey guys, first time poster and rookie coder here.

I am writing a very basic 2D ray tracer, this is purely a physics simulation, nothing to do with computer graphics. I am trying to track a bunch of rays though space, intersecting with a bunch of circular lenses, for now I am only concerned with refracted rays so the only maths I need is Snell’s law. I have the ray tracer working. But I just cant figure out why my loop structure is not working.

Can someone please have a look at my code and tell me why?

To test the code I was using a place holder variable for j and manually incrementing it which works, but when I try and automate that incrementation in a Do loop it doesn’t work.

Any help would be greatly appreciated.

P.S I’m sure this code is horrible, inefficient and riddled with coding and mathematical missteps so advice and suggestions are gladly welcomed.

Thanks

ClearAll["Global`*"]
(*draw circles*)
(*input data*)
LcircRad = 5;
ScircRad = 1;
n1 = 1; (*refractiv index of free space*)
n2 = 1.65; (*ref index of lens*)
(*angle per circle theta*)
theta = 2 ArcSin[ScircRad/(LcircRad + ScircRad)];
(*lenses per circle*)
nlens = Floor[2 \[Pi]/theta];
(*corrected angle per circle to ensure even distriubution*)
theta = N[2 \[Pi]/nlens];
(*populating arrays with x and y centres of small circles
ofset with \[Pi] so that first circle is drawn on the left*)
xcircle = 
  Table[(LcircRad + ScircRad) Cos[theta i + \[Pi]], {i, 0, nlens, 1}];
ycircle = 
  Table[(LcircRad + ScircRad) Sin[theta i + \[Pi]], {i, 0, nlens, 1}];
(*calculates number of rays*)
totalD = 2 (LcircRad + 
    2 ScircRad);(*total diameter of lens structure*)
raystarty = LcircRad + 2 ScircRad;(*starting height of first ray*)
raystartx = -3;(*starting x of rays*)
rayspacing = 1;(*vertical spacing between rays*)
nrays = Floor[totalD/rayspacing];
rayspacing = totalD/nrays;

(*populate the first intersection points & grad for ray plotting*)
(*
Do[yint[i,1]=raystarty-rayspacing(i-1),{i,1,nrays+1,1}]
Do[xint[i,1]=0,{i,1,nrays+1,1}]
Do[raygrad[i,1]=0,{i,1,nrays+1,1}]
*)

xint[1, 1] = -10;
yint[1, 1] = 5.3;
raygrad[1, 1] = 0;
j = 2;

(*
Do[
*)
y = raygrad [1, j - 1] (x - xint[1, j - 1]) + 
  yint[1, j - 
    1];(*y=m(x-x1)+y1 point grad form of straight line for ray at \
prev int point*)
(*temp stores x int of line and circles and ingnores the solution \
which is the prev intersection using cases, rounded so that cases can \
remove the right thing*)
xintemp =
 Select[
  Cases[
   Round[
    Flatten[Table[
      Select[
       NSolve[(x - xcircle[[i]])^2 + (y - ycircle[[i]])^2 == 
          ScircRad^2, x, WorkingPrecision -> 3][[All, 1, 2]],
       Element[#, Reals] &],
      {i, 1, nlens} ]],
    0.00001],
   Except[xint[1, j - 1]]]
  , # > xint[1, j - 1] &];
(*solves the line for y using x intercepts*)
yintemp = Round[
   Flatten[Table[
     Solve[
       ysol == raygrad [1, j - 1] (xintemp[[i]] - xint[1, j - 1]) + 
         yint[1, j - 1], ysol][[All, 1, 2]],
     {i, Dimensions[xintemp][[1]]}]]
   , 0.00001];
(*find coord of new intersection closest to prev*)
xint[1, j] = 
  Nearest[Table[{xintemp[[i]], yintemp[[i]]}, {i, 
      Dimensions[xintemp][[1]]}], {xint[1, j - 1], 
     yint[1, j - 1]}][[1, 1]];
yint[1, j] = 
  Nearest[Table[{xintemp[[i]], yintemp[[i]]}, {i, 
      Dimensions[xintemp][[1]]}], {xint[1, j - 1], 
     yint[1, j - 1]}][[1, 2]];
(*finds the circle number of the given intersection point*)
circno = Position[
    Round[
     Table[
      Select[
       NSolve[(x - xcircle[[i]])^2 + (y - ycircle[[i]])^2 == 
          ScircRad^2, x][[All, 1, 2]],
       Element[#, Reals] &],
      {i, 1, nlens, 1} ],
     0.00001],
    xint[1, j]][[1, 1]];
(*finds tangent and norm grad at point of intersection bw ray and \
circle*)
normgrad = (ycircle[[circno]] - yint[1, j])/(
  xcircle[[circno]] - xint[1, j]);
normang = ArcTan[normgrad];
alpha1 = ArcTan[
   Abs[(raygrad [1, j - 1] - normgrad)/(
    1 + raygrad [1, j - 1]*normgrad)]];
alpha2 = ArcSin[(n1/n2) Sin[alpha1]];

If[normang > 0, raygrad[1, j] = Tan[normang - alpha2], 
  raygrad[1, j] = Tan[normang + alpha2]];

(*
,{j,20}]
*)


?xint
?yint
?raygrad

Show[Table[
  Graphics[Circle[{xcircle[[i]], ycircle[[i]]}, ScircRad]], {i, 1, 
   nlens, 1}], Graphics[Circle[{0, 0}, LcircRad]],
 Table[
  Graphics[{Red, 
    Line[{{xint[1, o], yint[1, o]}, {xint[1, o + 1], 
       yint[1, o + 1]}}]}]
  , {o, 8}],
 PlotRange -> {{-1 totalD, 1 totalD}, {-1 totalD, 1 totalD}}, 
 Axes -> True]
2

If I insert the line

Print["i=", i];

just before the line

xintemp = Select[...NSolve[(x - xcircle[[i]])^2 +...

then mine prints

i=i

and so I assume you are not correctly initializing i inside that Do loop.

3

Hi there, I am sure you are right.

I has something to do with the Do loops inside the loop I want to implement, when I un-comment out the Do loop all the loops inside it loose their colour formatting, you know how the variable you are incrementing has a different shade of blue, that goes away.

Can you please expand on your answer as I have no idea what you mean by initialisation. How does one correctly initialise?

EDIT: Just to be clear if you meant initialise as in set a value for it, i is the incrementing variable for that Do loop, as such it shouldn’t have a value

Thanks

4

Hi there,

I just turned my machine on again after lunch and everything is working. Have no idea what was wrong before, this thread can be delete.

5

It would be cool to see what this code does if this is a ray tracer. But when I run it - it did not work - error:

Coordinate {$CellContextxint[1, 8], $CellContextyint[1, 8]} should
be a pair of numbers, or a Scaled or Offset form.

Do you have a fixed version?

6

Hello, yes I have fixed it!

If you want to change the parameters n2 is the refractive index of the lens material and the initial ray height can be set by changing yint[1,1]

What I have left to do is create another loop that will add in more rays.

If you have any feedback/improvements on the code I would love to hear it as this is my first coding project since high school, and first with Mathematica.

Here is what the code will output:

working.png
working.png540×540 6.99 KB

ClearAll["Global`*"]
(*draw circles*)
(*input data*)
LcircRad = 5;
ScircRad = 1;
n1 = 1; (*refractiv index of free space*)
n2 = 1.5; (*ref index of lens*)
(*angle per circle theta*)
theta = 2 ArcSin[ScircRad/(LcircRad + ScircRad)];
(*lenses per circle*)
nlens = Floor[2 \[Pi]/theta];
(*corrected angle per circle to ensure even distriubution*)
theta = N[2 \[Pi]/nlens];
(*populating arrays with x and y centres of small circles
ofset with \[Pi] so that first circle is drawn on the left*)
xcircle = 
  Table[(LcircRad + ScircRad) Cos[theta i + \[Pi]], {i, 0, nlens, 1}];
ycircle = 
  Table[(LcircRad + ScircRad) Sin[theta i + \[Pi]], {i, 0, nlens, 1}];
(*calculates number of rays*)
totalD = 2 (LcircRad + 
    2 ScircRad);(*total diameter of lens structure*)
raystarty = LcircRad + 2 ScircRad;(*starting height of first ray*)
raystartx = -3;(*starting x of rays*)
rayspacing = 1;(*vertical spacing between rays*)
nrays = Floor[totalD/rayspacing];
rayspacing = totalD/nrays;
(*populate the first intersection points & grad for ray plotting*)
(*
Do[yint[i,1]=raystarty-rayspacing(i-1),{i,1,nrays+1,1}]
Do[xint[i,1]=0,{i,1,nrays+1,1}]
Do[raygrad[i,1]=0,{i,1,nrays+1,1}]
*)

xint[1, 1] = -totalD;
yint[1, 1] = 5.3;
raygrad[1, 1] = 0;
j = 2;
xintemp = {1};
While[xintemp != {},
 
 y = raygrad [1, j - 1] (x - xint[1, j - 1]) + 
   yint[1, j - 
     1];(*y=m(x-x1)+y1 point grad form of straight line for ray at \
prev int point*)
 (*temp stores x int of line and circles and ingnores the solution \
which is the prev intersection using cases, rounded so that cases can \
remove the right thing*)
 xintemp =
  Select[
   Cases[
    Round[
     Flatten[Table[
       Select[
        NSolve[(x - xcircle[[i]])^2 + (y - ycircle[[i]])^2 == 
           ScircRad^2, x, WorkingPrecision -> 3][[All, 1, 2]],
        Element[#, Reals] &],
       {i, 1, nlens} ]],
     0.00001],
    Except[xint[1, j - 1]]]
   , # > xint[1, j - 1] &];
 (*solves the line for y using x intercepts*)
 yintemp = Round[
   Flatten[Table[
     Solve[
       ysol == raygrad [1, j - 1] (xintemp[[i]] - xint[1, j - 1]) + 
         yint[1, j - 1], ysol][[All, 1, 2]],
     {i, Dimensions[xintemp][[1]]}]]
   , 0.00001];
 (*find coord of new intersection closest to prev*)
 xint[1, j] = 
  Nearest[Table[{xintemp[[i]], yintemp[[i]]}, {i, 
      Dimensions[xintemp][[1]]}], {xint[1, j - 1], 
     yint[1, j - 1]}][[1, 1]];
 yint[1, j] = 
  Nearest[Table[{xintemp[[i]], yintemp[[i]]}, {i, 
      Dimensions[xintemp][[1]]}], {xint[1, j - 1], 
     yint[1, j - 1]}][[1, 2]];
 (*finds the circle number of the given intersection point*)
 circno = Position[
    Round[
     Table[
      Select[
       NSolve[(x - xcircle[[i]])^2 + (y - ycircle[[i]])^2 == 
          ScircRad^2, x][[All, 1, 2]],
       Element[#, Reals] &],
      {i, 1, nlens, 1} ],
     0.00001],
    xint[1, j]][[1, 1]];
 (*finds tangent and norm grad at point of intersection bw ray and \
circle*)
 normgrad = (ycircle[[circno]] - yint[1, j])/(
  xcircle[[circno]] - xint[1, j]);
 normang = ArcTan[normgrad];
 alpha1 = 
  ArcTan[Abs[(raygrad [1, j - 1] - normgrad)/(
    1 + raygrad [1, j - 1]*normgrad)]];
 alpha2 = ArcSin[(n1/n2) Sin[alpha1]];
 
 If[normang > 0, raygrad[1, j] = Tan[normang - alpha2], 
  raygrad[1, j] = Tan[normang + alpha2]];
 
 j = j + 1;]

(*
,{j,20}]
*)
(*after while loop is exited need to find intersection with ending \
line*)
y =.;
j = j - 1;
xint[1, j] = totalD;
yint[1, j] =
      raygrad [1, j - 1] (xint[1, j] - xint[1, j - 1]) + 
       yint[1, j - 1];
matrixer[functionName_Symbol] := 
 Normal[SparseArray[ReleaseHold[DownValues[#] /. # -> List]]] &[
  functionName]
(*turns xint and yint to matrix & removes the tailing entry which is \
bad cos of bad while loop
xintmatrix=matrixer[xint];
xintmatrix=Delete[xintmatrix,{1,Dimensions[xintmatrix][[2]]}];
yintmatrix=matrixer[yint];
yintmatrix=Delete[yintmatrix,{1,Dimensions[yintmatrix][[2]]}];*)
xintmatrix = matrixer[xint];
yintmatrix = matrixer[yint];
Show[Table[
  Graphics[Circle[{xcircle[[i]], ycircle[[i]]}, ScircRad]], {i, 1, 
   nlens, 1}], Graphics[Circle[{0, 0}, LcircRad]],
 Table[
  Graphics[{Red, 
    Line[{{xintmatrix[[1, o]], 
       yintmatrix[[1, o]]}, {xintmatrix[[1, o + 1]], 
       yintmatrix[[1, o + 1]]}}]}]
  , {o, Dimensions[xintmatrix][[2]] - 1}],
 PlotRange -> {{-1 totalD, 1 totalD}, {-1 totalD, 1 totalD}}, 
 Axes -> True]
7

This is awesome! What does it model? Are these cylindrical lenses arranged circularly? Most importantly - does it model some real device?