Electromagnetic simulation at various timescales: https://twitter.com/bencbartlett/status/1369396941730312193

coherent_incoherent_light.jl

# Electromagnetic simulation at varying timescales - simulation code to make this animation:
# https://twitter.com/bencbartlett/status/1369396941730312193
#
# Animation idea was inspired by u/cenit997's (Twitter: @__cenit) Reddit post:
# https://www.reddit.com/r/Python/comments/hxmhai/a_simulation_of_how_an_incoherent_light_source/


#---
using LinearAlgebra, Base.Threads, BenchmarkTools, LoopVectorization, CSV, DataFrames, Printf, Plots

#---
const c = 3.0 * 10^8 * 10^9 # distance units are in nanometers
const maxR = 5000.0
const diskR = 1500.0
const numInternalSources = 30;
const numEdgeSources = 20;
const numSources = numInternalSources + numEdgeSources;
const maxIntensity = sqrt(numSources / diskR^2) # not used here, just in Mathematica code

#---

# Generate a list of points inside of a disk of radius R
function randomPointsInDisk(R, n=1)
    θ = 2π * rand(n)
    r = R * sqrt.(rand(n))
    x = r .* cos.(θ)
    y = r .* sin.(θ)
    return x, y
end

# Generate a list of points near the edge of a disk of radius R
function pointsAtEdgeOfDisk(R, n=1)
    θ = 2π * collect(range(0, length=n)) / n + rand(n) / (.2 * n)
    r = R * (ones(n) + .05 * rand(n))
    x = r .* cos.(θ)
    y = r .* sin.(θ)
    return x, y
end


#---

# Assemble all the properties of the sources
sourcesX_in, sourcesY_in = randomPointsInDisk(diskR, numInternalSources)
sourcesX_out, sourcesY_out = pointsAtEdgeOfDisk(diskR, numEdgeSources)
sourcesXY = (cat(dims=1, sourcesX_in, sourcesX_out), cat(dims=1, sourcesY_in, sourcesY_out))
sourcesPower = 1 .+ 0.0 * rand(numSources)
sourcesPhase = 2π * rand(numSources)
sourcesλ = 750 .+ 2.0 * rand(numSources)



#---

# Compute the field matrix of a source at point x₀, y₀ with given properties
# Julia threading is fast, so it's actually faster to thread in this function than
# in the getIntensityAverage function. Be sure to `set JULIA_NUM_THREADS=16` before
# running this script for fastest evaluation.
function getField(x₀, y₀, P, ϕ, λ, t, δx=50)

    xyiter = collect(enumerate(-maxR:δx:maxR))

    E = similar(x₀, length(xyiter), length(xyiter))

    for (i,x) in xyiter
        @threads for (j,y) in xyiter
            r = @. sqrt((x-x₀)^2 + (y-y₀)^2)
            E[i,j] = sum(@avx @. P / r * cos(ϕ + 2π/λ * (c*t-r)))
        end
    end

    E

end

# MC integration over nSamples number of points
function getIntensityAverage(x0, y0, P, ϕ, λ, t0, δt, nSamples, δx=50)

    size = length(-maxR:δx:maxR)
    intensities_sum = zeros(size, size)

    for i = 1:nSamples
        t = t0 + δt * rand()
        intensities_sum += getField(x0, y0, P, ϕ, λ, t, δx).^2
    end

    return intensities_sum / nSamples

end

#---
println("Threads: ", Threads.nthreads())
@btime getField(sourcesXY[1], sourcesXY[2], sourcesPower, sourcesPhase, sourcesλ, 0.0, 50)
@btime getIntensityAverage(sourcesXY[1], sourcesXY[2], sourcesPower, sourcesPhase, sourcesλ, 0.0, 10^-16, 10, 50)


#---

# Assemble list of sliding timesteps
anim_start = -16.0 .* ones(60 * 5)
anim_accel_1 = -16.0 .+ (-13.5 + 16.0) .* sin.(range(0.0, π/2, length=60 * 4)).^2
anim_pause_1 = -13.5 .* ones(60 * 8)
anim_accel_2 = -13.5 .+ (-6.0 + 13.5) .* sin.(range(0.0, π/2, length = 60 * 7)).^2
anim_pause_2 = -6.0 .* ones(60 * 1)

δtlist = 10.0 .^ cat(dims=1, anim_start, anim_accel_1, anim_pause_1, anim_accel_2, anim_pause_2)
n_frames = length(δtlist)

tlist = cumsum(δtlist)


#---

# Run everything and output data to CSV files to render in Mathematica
for (i, (t, δt)) in enumerate(zip(tlist, δtlist))
    println("Frame "*string(i)*"/"*string(n_frames), " t="*string(t), " δt="*string(δt))
    numSamples = floor(Int, min(ceil(2 * 10.0^5 * δt^.3), 500))
    intensities = getIntensityAverage(sourcesXY[1], sourcesXY[2], sourcesPower, sourcesPhase, sourcesλ, t, δt, numSamples, 50)
    fname = @sprintf("%.20f", t)
    CSV.write("data/" * fname * ".csv", DataFrame(intensities), header=false)
end

coherent_incoherent_light.nb

(*
Electromagnetic simulation at varying timescales - frame-rendering code to make this animation:
https://twitter.com/bencbartlett/status/1369396941730312193

Animation idea was inspired by u/cenit997's (Twitter: @__cenit) Reddit post:
https://www.reddit.com/r/Python/comments/hxmhai/a_simulation_of_how_an_incoherent_light_source/
*)

SetDirectory[NotebookDirectory[]];
Needs["AnimationTools`"]; (*my personal collection of animation \
tools. The only method used from here is ParallelProgressDo in the \
last cell. Replace with ParallelDo to run this code.*)

filenames = FileNames["data/*.csv"];

WINDOWS = False;(*directories separated by \ or /, change to False if \
running on Mac*)

clight = 3.0*10^8* 10^9 ;(*nm/s*)
maxR = 5000.0;(*nm*)

diskR = 1500;(*nm*)
numSources = 50;
maxIntensity = Sqrt[numSources/diskR^2];

(*like ScientificForm[] but keeps trailing zeros*)

ScientificFormFixedZeroes[n_, decimals_] := Module[{num, exp},
   exp = Floor[Log10[n]];
   num = n / 10^exp;
   If[StringLength[
      ToString[NumberForm[num, {decimals + 1, decimals}]]] > 
     decimals + 2, exp++; 
    num /= 10.0 ]; (*prevents 10.00\[Times]10^-5 bug*)
   
   ToString[NumberForm[num, {decimals + 1, decimals}]] <> "\[Times]" <>
     ToString[DisplayForm@SuperscriptBox["10", ToString[exp]], 
     StandardForm]
   ];

renderFrame[index_] := 
  Module[{tstring, frame, intensities, fname = filenames[[index]], 
    tCurr, tPrev, \[Delta]t, plot, gauge, dataDelete},
   dataDelete = If[WINDOWS, "data\\", "data/"];
   tstring = StringDelete[StringDelete[fname, dataDelete], ".csv"];
   intensities = Import[fname];
   (*Print[ StringDelete[StringDelete[fname,dataDelete],".csv"]];*)
  
    tCurr = 
    ToExpression[ 
     StringDelete[StringDelete[fname, dataDelete], ".csv"]];
   tPrev = 
    If[index > 1, 
     ToExpression[ 
      StringDelete[StringDelete[filenames[[index - 1]], dataDelete], 
       ".csv"]], 0.0];
   \[Delta]t = tCurr - tPrev;
   ListDensityPlot[Sqrt[intensities],
    DataRange -> {{-maxR, maxR}, {-maxR, maxR}},
    PlotRange -> {Automatic, Automatic, {0, .9*maxIntensity}},
    ClippingStyle -> Automatic,
    Frame -> False,
    ColorFunction -> "SunsetColors", ColorFunctionScaling -> True,
    Background -> Black,
    ImageSize -> 1080,
    Epilog -> { 
      Inset[
       VerticalGauge[Log10[\[Delta]t], {-16, -6},
        TicksStyle -> Directive[White, 28],
        LabelStyle -> Directive[White, 28],
        GaugeMarkers -> "Classic",
        GaugeStyle -> Directive[White, Scaled[.15]],
        GaugeLabels -> {
          
          Placed[Evaluate[
            "\[Delta]t=" <> 
             ScientificFormFixedZeroes[\[Delta]t, 
              1]], {Log10[\[Delta]t], {1.32, .5}}]
          },
        ImageSize -> {Automatic, 600}]
       , ImageScaled[{.86, .5}]],
      Inset[
       Text[Style[
         "\!\(\*SubscriptBox[\(log\), \(10\)]\) timestep (s)", White, 
         28]], ImageScaled[{.87, .77}]],
      Inset[
       Text[Style[
         Evaluate[
          "t = " <> ScientificFormFixedZeroes[tCurr, 2] <> "s"], 
         White, 34]], ImageScaled[{.07, .93}], Left],
      {Thickness[0.007], White, 
       Line[{{-maxR + 500, -maxR + 500}, {-maxR + 1500, -maxR + 
           500}}]},
      Inset[Text[Style["1000 nm", White, 28]], 
       ImageScaled[{.064, 1 - .95}], Left]
      }]
   ];
saveFrame[index_] := 
  Module[{tstring, frame, intensities, dataDelete, fname},
   fname = filenames[[index]];
   dataDelete = If[WINDOWS, "data\\", "data/"];
   tstring = StringDelete[StringDelete[fname, dataDelete], ".csv"];
   frame = renderFrame[index];
   Export["frames/" <> tstring <> ".png", frame];
   ];
renderFrame[42]

(*ProgressParallelDo[saveFrame[i],{i,Length[filenames]}]*) (*method \
from my AnimationTools package*)
ProgressDo[
 saveFrame[i], {i, Length[filenames]}]