# Using Regions, can I model a reflecting wavefront?

I was wondering if I could get Regions and Mathematica's shapes to do all the hard work for me in making a "droplet in a pond" simulation. I do not want the "waves" to interact. I've had moderate success. The black dot is the source. I was wondering if anyone has any good ideas about how to deal with a partial reflection, see the black line - can I get the circles to reflect from it.... (* set up the region *)
source = Point [{0, 0}];
reflectionsphere = Disk[{0, 0}, 10];

(* step size *)
stepsize = 1;
max = 20;

animation = ConstantArray[{}, max]

activesurface = Disk[{0, 0}, 1]

For[i = 1, i <= 20, i = i + stepsize,

If[RegionIntersection[activesurface, reflectionsphere] ===
reflectionsphere,
circlepoints = CirclePoints[{0, 0}, 10, 6];
circlelines2 = Disk[#, i - 10] & /@ circlepoints;
wavefronts = RegionUnion[circlelines2];
wavefrontsinsphere =
RegionIntersection[reflectionsphere, wavefronts];

animation[[i]] =
Graphics[{source,
RegionBoundary[reflectionsphere], {Opacity[0.2], Red,
MeshPrimitives[DiscretizeRegion[wavefrontsinsphere], 2]}}]

,

activesurface = Disk[{0, 0}, i];
animation[[i]] =
Graphics[{source,
RegionBoundary[reflectionsphere], {Opacity[0.2], Red,
activesurface}}];

];

];

animationgif = ListAnimate[animation]

• Imagine if you were shining a torch, and wanted to see where the light would go...or released a unidirectional fart and wanted to see how it would spread...the torch and fart don't have any coherency, and therefore wouldn't have any sort of superposition - so double fart doesn't mean anything....the two waves were just an example. I don't want any waves to interact with each other...
– Tomi
Jan 8, 2020 at 19:23

I adapted @Kuba's approach from this answer to generate a quick and dirty particle tracer.

(* Create and Discretize Region *)
region = RegionDifference[Disk[],
Rectangle[{-1/3, -1/3}, {1/3, -1/4}]];
R2 = RegionBoundary@DiscretizeRegion@region;
rdf = RegionDistance[R2];
rnf = RegionNearest[R2];
(* Time Increment *)
dt = 0.001;
(* Collision Margin *)
margin = 1.05 dt;
r0 = 1000;
(* Starting Point for Emission *)
sp = {0, 0};
Block[{xnew = x + dt v}, {rdf[xnew], xnew, v, c}] /; r > margin
Block[{xnew = x , vnew = v, normal = Normalize[x - rnf[x]]},
vnew = Normalize[v - 2 v.normal normal];
xnew += dt vnew;
{rdf[xnew], xnew, vnew, c + 1}] /; r <= margin


Now, we can run the simulation and create an animation at every 50 time steps.

nparticles = 1000;
ntimesteps = 5000;
tabres = Table[
NestList[
sp, {Cos[2 Pi #], Sin[2 Pi #]} &@RandomReal[], 0},
ntimesteps], {i, 1, nparticles}];
frames = Table[
RegionPlot[R2, Epilog -> (Disk[#, 0.01] & /@ tabres[[All, i, 2]]),
AspectRatio -> Automatic], {i, 1, ntimesteps, 50}];
ListAnimate@frames You could use NDSolve to do the hard work:

region = Disk[];
sol = NDSolveValue[{D[u[t, x, y], {t, 2}] -
Laplacian[u[t, x, y], {x, y}] == 0,
DirichletCondition[u[t, x, y] == 0, True],
u[0, x, y] == 2*Exp[-125 ((x)^2 + (y - 0.5)^2)],
Derivative[1, 0, 0][u][0, x, y] == 0}, u, {t, 0, 2},
Element[{x, y}, region]]


And then:

ListAnimate[
Table[Plot3D[sol[t, x, y], Element[{x, y}, region],
PlotRange -> {-0.75, 2}, AspectRatio -> Automatic, Boxed -> False,
Axes -> None, PlotPoints -> 33], {t, 0, 2, 1/25}],
SaveDefinitions -> True] To have an internal obstacle just change the region:

region = RegionDifference[Disk[],
Rectangle[{-1/3, -1/3}, {1/3, -1/4}]];
sol = NDSolveValue[{D[u[t, x, y], {t, 2}] -
Laplacian[u[t, x, y], {x, y}] == 0,
DirichletCondition[u[t, x, y] == 0, True],
u[0, x, y] == 2*Exp[-125 ((x)^2 + (y - 0.5)^2)],
Derivative[1, 0, 0][u][0, x, y] == 0}, u, {t, 0, 2},
Element[{x, y}, region]]


Visualize:

ListAnimate[
Table[
Plot3D[sol[t, x, y], Element[{x, y}, region],
PlotRange -> {-0.75, 2}, AspectRatio -> Automatic, Boxed -> False,
Axes -> None, PlotPoints -> 33], {t, 0, 2, 1/25}],
SaveDefinitions -> True] Also, you can find much more information on the wave equation by looking at the Acoustics in the Time Domain tutorial in the documentation system under PDEModels/tutorial/AcousticsTimeDomain

• Amazing response. My hope, however, was to avoid the superposition in involved in using the NDSolve and a differential equation. I was hoping to map out the wavefront without interference from other waves. This means I need to ignore the additivity of the waves (e.g. wave 1 - wave2 = n). I'm not sure this is possible with this solution. I want to map out the wavefronts when they don't interact. Physically, if helpful, imagine shining a single pulse of light in all direction with a diffuse reflecting sphere instead of a drop in a pond.....
– Tomi
Jan 8, 2020 at 11:54
• @Tomi this is then just a change in your initial conditions correct? Jan 8, 2020 at 20:02
• Awesome response, thanks! Silly question: I tried parallelizing the code and replaced the Table with ParallelTable (Math 11.3) and instead of making the code much faster, it made it much slower! Why was that? Thanks in advance for your response! Jan 9, 2020 at 11:52
• @CATrevillian, as a poor approximation you could make all u>0 and make reflections always positive, so you although you would still retain the linearity of the wave, it wouldn't ever be negative and would ignore some interference, but I really just wanted to have a single line moving through the sphere.....
– Tomi
Jan 9, 2020 at 12:26

@user21's solution is very impressive. However, it isn't quite what I was looking for. This is because of the interaction between the waves. They are acting - well - like waves. This means we have a linear addition of the waves. This was what the original question forbid ;) . We want a single wavefront to come from the centre of the sphere and watch what happens as it moves around objects. Imagine it is a single photon - and doesn't act like a water wave.

Of course, if we're talking about single photons - a raytracing solution would work. I've implemented one (inspired from here), however, again - it isn't what the original question is asking for. We want a single wavefront which spreads....

But, anyway, this is my ray-tracing attempt

With 3 photons: With 100 photons: (* Line Intersection *)

LLI[vi_List] :=
With[{x1 = vi[[1, 1]], y1 = vi[[1, 2]], x2 = vi[[2, 1]],
y2 = vi[[2, 2]], x3 = vi[[3, 1]], y3 = vi[[3, 2]], x4 = vi[[4, 1]],
y4 = vi[[4,
2]]}, {-((-(x3 - x4) (x2 y1 - x1 y2) + (x1 - x2) (x4 y3 -
x3 y4))/((x3 - x4) (y1 - y2) + (-x1 + x2) (y3 -
y4))), (x4 (y1 - y2) y3 + x1 y2 y3 - x3 y1 y4 - x1 y2 y4 +
x3 y2 y4 +
x2 y1 (-y3 + y4))/(-(x3 - x4) (y1 - y2) + (x1 - x2) (y3 - y4))}]

(* Consider how we bounce *)

bounce2[{p0_, d0_, i0_}] :=
Module[{idxL, pL, validL, distL, i, p1, d1, bValid, dist, angleL,
angle}, idxL =
Position[Pi/2 < VectorAngle[d0, #] < Pi 3/2 Pi & /@ norm, True] //
Flatten;
pL = Table[LLI[{p0, p0 + d0, ##}] & @@ edge[[j]], {j, idxL}];
validL =
Table[! Or @@ (Greater[#,
1] & /@ (EuclideanDistance[#, pL[[i]]]/
length[[idxL[[i]]]] & /@ edge[[idxL[[i]]]])), {i,
Length@idxL}];
distL = EuclideanDistance[#, p0] & /@ pL;
angleL =
Table[VectorAngle[norm[[idxL[[i]]]], pL[[i]] - p0], {i,
Length@idxL}];
{i, p1, bValid, angle, dist} =
Select[Transpose@{idxL, pL, validL, angleL,
distL}, (#[] && #[] > Pi/2) &] //
MinimalBy[#, Last] & // #[] &;
d1 = (ReflectionTransform[RotationTransform[-Pi/2]@(-norm[[i]]),
p1]@p0 - p1) // Normalize;
{p1, d1, i}]

(* Give our boundaries *)

boundary1 = CirclePoints[2, 100];

edge1 = Table[
RotateRight[boundary1, i][[;; 2]], {i, Length@boundary1}];
length1 = EuclideanDistance @@ # & /@ edge1;
norm1 = Normalize@(RotationTransform[Pi/2]@(#[] - #[])) & /@
edge1;

boundary2 = {{-1, -0.2}, {1, -0.2}, {1, 0}, {-1, 0}};

edge2 = Table[
RotateRight[boundary2, i][[;; 2]], {i, Length@boundary2}];
length2 = EuclideanDistance @@ # & /@ edge2;
norm2 = -Normalize@(RotationTransform[Pi/2]@(#[] - #[])) & /@
edge2;

boundary = Join[boundary1, boundary2];
edge = Join[edge1, edge2];
length = Join[length1, length2];
norm = Join[norm1, norm2];

photons = 3;
bounces = 100;
g = ConstantArray[{}, photons];

For[i = 1, i <= photons, i++,
p0 = {0, 0.1};
d0 = {Cos@#, Sin@#} &@RandomReal[{0, 2 Pi}];
r = NestList[bounce2, {p0, d0, 0}, bounces];
p = r[[All, 1]];
g[[i]] =
Table[Graphics[{FaceForm[LightBlue], EdgeForm[], Gray,
Line@p[[;; j]], Darker@Gray, Point@p[[;; j]], Red,
Point@p[]}], {j, 2, Length@r}];
]

surface =
Graphics[{{FaceForm[LightBlue], Polygon@boundary1}, FaceForm[White],
Polygon@boundary2}]
animate = Table[Show[surface, g[[;; , {i}]]], {i, 1, bounces}];

ListAnimate[animate]


This isn't a complete solution as I really looking for the propagation of the circles around the sphere.