Last active
February 23, 2024 04:22
-
-
Save neozhaoliang/15d0350450885b90c5efbc73ddcdbdde to your computer and use it in GitHub Desktop.
Caustics in a reflective ring
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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')}") |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
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 | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
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 | |
} | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* | |
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 | |
} | |
} |
Author
neozhaoliang
commented
Feb 7, 2024
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment