Skip to content

Instantly share code, notes, and snippets.

@bencbartlett
Created March 10, 2021 20:41
Show Gist options
  • Save bencbartlett/2bec845fe7efdf9ca4c2de2fc4aea179 to your computer and use it in GitHub Desktop.
Save bencbartlett/2bec845fe7efdf9ca4c2de2fc4aea179 to your computer and use it in GitHub Desktop.
Electromagnetic simulation at various timescales: https://twitter.com/bencbartlett/status/1369396941730312193
# 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
(*
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]}]
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment