16
$\begingroup$

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.

enter image description here

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....

enter image description here

(* 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]
$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – Tomi Jan 8 at 19:23
7
$\begingroup$

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};
(* Conditional Particle Advancer *)
advance[r_, x_, v_, c_] := 
 Block[{xnew = x + dt v}, {rdf[xnew], xnew, v, c}] /; r > margin
advance[r_, x_, v_, c_] := 
 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[
    advance @@ # &, {rdf[sp], 
     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

Particle Animation

|improve this answer|||||
$\endgroup$
21
$\begingroup$

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]

enter image description here

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]

enter image description here

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

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ 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..... $\endgroup$ – Tomi Jan 8 at 11:54
  • $\begingroup$ @Tomi this is then just a change in your initial conditions correct? $\endgroup$ – CA Trevillian Jan 8 at 20:02
  • $\begingroup$ 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! $\endgroup$ – Rebel-Scum Jan 9 at 11:52
  • $\begingroup$ @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..... $\endgroup$ – Tomi Jan 9 at 12:26
3
$\begingroup$

@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:

enter image description here

With 100 photons:

enter image description here

(* 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}, (#[[3]] && #[[4]] > Pi/2) &] // 
     MinimalBy[#, Last] & // #[[1]] &;
  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]@(#[[2]] - #[[1]])) & /@ 
   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]@(#[[2]] - #[[1]])) & /@ 
   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[[1]]}], {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.

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.