Skip to content

Instantly share code, notes, and snippets.

@neozhaoliang
Last active February 23, 2024 04:22
Show Gist options
  • Save neozhaoliang/15d0350450885b90c5efbc73ddcdbdde to your computer and use it in GitHub Desktop.
Save neozhaoliang/15d0350450885b90c5efbc73ddcdbdde to your computer and use it in GitHub Desktop.
Caustics in a reflective ring
from sympy import *
x, y, X, Y = symbols("x y X Y")
P = x**2 + y**2 - 1
dx = diff(P, x) # gradient of P
dy = diff(P, y)
curve = Matrix([x, y])
light_source = Matrix([1, 0])
l = curve - light_source # the incident ray
n = Matrix([dx, dy]) # the normal vector
r = simplify(l - 2 * l.dot(n) * n / n.dot(n)) # the reflected ray
rx, ry = r
dxrx = diff(rx, x)
dyrx = diff(rx, y)
dxry = diff(ry, x)
dyry = diff(ry, y)
denominator = ry * (dxrx * dy - dyrx * dx) - rx * (dxry * dy - dyry * dx)
nominator = dx * l[0] + dy * l[1]
F = (X - x) * denominator - nominator * rx
G = (Y - y) * denominator - nominator * ry
eqs = [eq.as_numer_denom()[0] for eq in [F, G, P]]
gb = groebner(eqs, [x, y, X, Y])
print(gb)
from sympy import *
t, X, Y = symbols("t X Y")
C = Matrix([cos(t), sin(t)]) # curve
light = Matrix([1, 0]) # light source
# C = Matrix([(2*cos(t) + cos(2*t)) / 3, (2*sin(t) + sin(2*t)) / 3])
# light = Matrix([S('-1/3', evaluate=False), 0])
l = C - light # incident ray
dx, dy = diff(C, t)
n = Matrix([dy, -dx]) # normal vector
r = simplify(l - 2 * l.dot(n) * n / n.dot(n)) # reflected ray
F = (Y - y) * r[0] - (X - x) * r[1]
dF = diff(F, t)
result = solve((F, dF), X, Y) # solve the envelope
print(f"X(t)={trigsimp(result[X], method='groebner')}")
print(f"Y(t)={trigsimp(result[Y], method='groebner')}")
/*
Code adapted from Nicolas Rougier's "ring.pov" file.
To render this file, you need to have POV-Ray 3.7 installed.
In terminal, run:
povray ring.pov +w800 +h600 +fn +am2 +a0.01
The scene contains:
1. A plane y = -2 with a checkerboard pattern.
2. A parametric surface defined by the functions:
x(u, v) = (2*cos(u) + cos(2*u)) / 3
y(u, v) = v
z(u, v) = (2*sin(u) + sin(2*u)) / 3
This is a cylinder with a cardioid cross-section. from y=-1 to y=1.
3. Two sphere_sweep objects that are placed on the top and bottom rims of the cylinder.
Note that the spheres do not have the same radius, hence the bottom rim does not
lie on the same horizontal plane. This may cause some artifacts in the rendering.
*/
#version 3.7;
#include "colors.inc"
#include "math.inc"
global_settings{
assumed_gamma 1.0
max_trace_level 25
photons {
spacing .01
autostop 0
gather 0, 200
jitter 0.4
}
radiosity {
pretrace_start 0.08
pretrace_end 0.01
count 600
error_bound .25
nearest_count 8
recursion_limit 1
gray_threshold 0
minimum_reuse 0.015
brightness 1.0
adc_bailout 0.01/2
}
}
#include "screen.inc"
#declare EyePos = <0,60,10>;
#declare EyeLook = <0,0,0>;
#declare EyeAngle = 40;
Set_Camera(EyePos, EyeLook, EyeAngle)
#macro Cx(T)
(2*cos(T) + cos(2*T)) / 3
#end
#macro Cz(T)
(2*sin(T) + sin(2*T)) / 3
#end
#macro param_surf(sc)
parametric {
function { (2*cos(u) + cos(2*u)) / 3 }
function { v }
function { (2*sin(u) + sin(2*u)) / 3 }
<0, -1>, <2*pi, 1>
scale <sc, 1.0, sc>
}
#end
#macro get_radius(T)
#local rad = vlength(<Cx(T), 0, Cz(T)>);
rad
#end
#macro tube(sc, ht)
#local num_segments = 120;
sphere_sweep {
cubic_spline
num_segments + 3
<sc * Cx(0), ht, sc * Cz(0)>, get_radius(0)
#local j = 0;
#while (j < num_segments + 1)
#local T = 2*pi*j/num_segments;
#local P = <sc * Cx(T), ht, sc * Cz(T)>;
#local Q = <Cx(T), 0, Cz(T)>;
#local rad = vlength(Q);
#local j = j + 1;
P, rad
#end
<sc * Cx(0), ht, sc * Cz(0)>, get_radius(0)
}
#end
light_source {
<10, 16, 0> // will be rotated to <-16, 10, 0>
rgb <0.75, 1, 1>
fade_distance 50 fade_power 2
area_light <10, 0, 0> <0, 0, 10> 20,20 adaptive 0 jitter circular orient
rotate z*90
photons {
reflection on
}
}
plane {
y, -2
pigment {
checker
color rgb 1, color rgb 0.8
}
finish {
reflection 0.2
diffuse 0.3
specular 0.4
}
photons {
target
collect on
reflection off
refraction off
}
}
union {
#local sc = 14;
object { param_surf(sc) }
object { param_surf(sc + 2) }
object { tube(sc + 1, -1) }
object { tube(sc + 1, 1) }
pigment {
rgb .1
}
finish {
reflection .9
specular 3
roughness 0.0025
ambient 0
diffuse 1
}
photons {
target
reflection on
collect on
}
translate -3 * x
}
/*
Code adapted from Nicolas Rougier's "ring.pov" file:
https://www.labri.fr/perso/nrougier/downloads/ring.pov
To render this file, you need to have POV-Ray 3.7 installed.
In terminal, run:
povray ring.pov +w800 +h600 +fn +am2 +a0.01
*/
#version 3.7;
#include "colors.inc"
#include "shapes.inc"
#include "glass.inc"
global_settings{
assumed_gamma 1.0
max_trace_level 25
photons {
count 8000000
spacing .01
autostop 0
gather 0, 200
jitter 0.4
}
radiosity {
pretrace_start 0.08
pretrace_end 0.01
count 600
error_bound .25
nearest_count 8
recursion_limit 1
gray_threshold 0
minimum_reuse 0.015
brightness 1.0
adc_bailout 0.01/2
}
}
#include "screen.inc"
#declare EyePos = <0,60,10>;
#declare EyeLook = <0,0,0>;
#declare EyeAngle = 40;
Set_Camera(EyePos, EyeLook, EyeAngle)
light_source {
<20, 30, 0>
rgb <0.75, 1, 1>
fade_distance 50 fade_power 2
area_light <10, 0, 0> <0, 0, 10> 20,20 adaptive 0 jitter circular orient
rotate z*90
photons {
reflection on
}
}
plane {
y,-2
pigment {
checker
color rgb 1, color rgb 0.8
}
finish {
reflection 0.2
diffuse 0.3
specular 0.4
}
photons {
target
collect on
reflection off
refraction off
}
}
union {
torus {
11, 1
translate y
}
torus {
11, 1
translate -y
}
cylinder {
-y, y, 12
open
}
cylinder {
-y, y, 10
open
}
pigment {
rgb .1
}
finish {
reflection .9
specular 3
roughness 0.0025
ambient 0
diffuse 1
}
photons {
target
reflection on
collect on
}
}
/*
Code adapted from Nicolas Rougier's "ring.pov" file.
To render this file, you need to have POV-Ray 3.7 installed.
In terminal, run:
povray ring.pov +w800 +h600 +fn +am2 +a0.01
*/
#version 3.7;
#include "colors.inc"
#include "math.inc"
global_settings{
assumed_gamma 1.0
max_trace_level 25
photons {
spacing .01
autostop 0
gather 0, 200
jitter 0.4
}
radiosity {
pretrace_start 0.08
pretrace_end 0.01
count 600
error_bound .25
nearest_count 8
recursion_limit 1
gray_threshold 0
minimum_reuse 0.015
brightness 1.0
adc_bailout 0.01/2
}
}
#include "screen.inc"
#declare EyePos = <0,60,10>;
#declare EyeLook = <0,0,0>;
#declare EyeAngle = 40;
Set_Camera(EyePos, EyeLook, EyeAngle)
#macro MyPrism(n, sc)
prism {
linear_sweep
linear_spline
-1, 1, n,
#local j=0;
#while (j<n)
<cos(2*pi*j/n), sin(2*pi*j/n)>
#local j = j + 1;
#end
open
scale <sc, 1, sc>
}
#end
#macro rim(n, sc, ht)
sphere_sweep {
linear_spline
n + 1,
#local j=0;
#while (j<n +1)
<sc*cos(2*pi*j/n), ht, sc*sin(2*pi*j/n)>, 1
#local j = j + 1;
#end
}
#end
light_source {
<10, 30, 0> // will be rotated to <-16, 10, 0>
rgb <0.75, 1, 1>
fade_distance 50 fade_power 2
area_light <10, 0, 0> <0, 0, 10> 20,20 adaptive 0 jitter circular orient
rotate z*90
photons {
reflection on
}
}
plane {
y, -2
pigment {
checker
color rgb 1, color rgb 0.8
}
finish {
reflection 0.2
diffuse 0.3
specular 0.4
}
photons {
target
collect on
reflection off
refraction off
}
}
union {
#local sc = 14;
#local n = 6;
MyPrism(n, sc + 2)
MyPrism(n, sc)
rim(n, sc + 1, 1)
rim(n, sc + 1, -1)
pigment {
rgb .1
}
finish {
reflection .9
specular 3
roughness 0.0025
ambient 0
diffuse 1
}
photons {
target
reflection on
collect on
}
}
@neozhaoliang
Copy link
Author

caustics

@neozhaoliang
Copy link
Author

caustics

@neozhaoliang
Copy link
Author

caustic_cardioid

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment