Simulating the Action Principle
Refath Bari
What is the path of least time between these two points in a gravitational field?
?
Snell's Law
Snell's Law
Snell's Law
Snell's Law
Snell's Law
Snell's Law
Initial Light Ray
Snell's Law
Initial Light Ray
Interface between Mediums
Snell's Law
Initial Light Ray
Interface between Mediums
Refracted Light Ray
Snell's Law
Snell's Law
Snell's Law
Initial Light Ray
Interface between Mediums
Refracted Light Ray
Snell's Law
Initial Light Ray
Interface between Mediums
Refracted Light Ray
Initial Light Ray
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(pos = start,
axis = size * v1.rotate(radians(angle)),
shaftwidth = 0.05,
headwidth = 0.2,
headlength = 0.2,
color = color.yellow)
Initial Light Ray
Interface between Mediums
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(pos = start,
axis = size * v1.rotate(radians(angle)),
shaftwidth = 0.05,
headwidth = 0.2,
headlength = 0.2,
color = color.yellow)
Initial Light Ray
Interface between Mediums
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(pos = start,
axis = size * v1.rotate(radians(angle)),
shaftwidth = 0.05,
headwidth = 0.2,
headlength = 0.2,
color = color.yellow)
def createInterface(start, size):
v1 = vec(0,-1,0)
return arrow(pos = start,
axis = size * v1.rotate(pi/2),
shaftwidth = 0.05,
headwidth = 0.001,
headlength = 0.001,
opacity = 0.2)
Initial Light Ray
Interface between Mediums
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(...)
def createInterface(start, size):
v1 = vec(0,-1,0)
return arrow(...)
Initial Light Ray
Interface between Mediums
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(...)
def createInterface(start, size):
v1 = vec(0,-1,0)
return arrow(...)
Refracted Light Ray
Initial Light Ray
Interface between Mediums
def createRay(start, angle, size):
v1 = vec(0,-1,0)
return arrow(...)
def createInterface(start, size):
v1 = vec(0,-1,0)
return arrow(...)
Refracted Light Ray
def refractRay(refIndex1, refIndex2, theta_inc, vector_inc):
theta_ref = asin((refIndex1/refIndex2) * sin(radians(theta_inc)))
return arrow(pos = vector_inc.pos+vector_inc.axis,
axis = vector_inc.axis.rotate(radians(theta_ref)),
shaftwidth = 0.05,
headwidth = 0.2,
headlength = 0.2,
color = color.yellow)
Refracted Light Ray
def refractRay(refIndex1, refIndex2, theta_inc, vector_inc):
theta_ref = asin((refIndex1/refIndex2) * sin(radians(theta_inc)))
return arrow(pos = vector_inc.pos+vector_inc.axis,
axis = vector_inc.axis.rotate(radians(theta_ref)),
shaftwidth = 0.05,
headwidth = 0.2,
headlength = 0.2,
color = color.yellow)
Expectation:
Reality:
Comparison:
Comparison:
What's Wrong?
Aha! We are rotating the incident ray to an angle, not a standard angle
Plan
Plan
Plan
Plan
Plan
Plan
Plan
Plan
Plan
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Array
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
- Fill up refractive index array manually, with the following:
refIndices = [1.52, 1.32, 1.32, 1.18, 1.18, 1.08, 1.08, 1]
- Indices (except first and last) are repeated twice because of the nature of the refractRay() function
- Do you see a brachistochrone-looking cycloid shape? If so, all is well. If not, debugging time!
Refractive Indices Expectation
- Fill up refractive index array manually, with the following:
refIndices = [1.52, 1.32, 1.32, 1.18, 1.18, 1.08, 1.08, 1]
- Indices (except first and last) are repeated twice because of the nature of the refractRay() function
- Do you see a brachistochrone-looking cycloid shape? If so, all is well. If not, debugging time!
Refractive Indices Expectation
Results
Expected
Observed
Results
Questions
Can the results be refined by increasing the number of iterations?
Refine Attempt
refIndices = [1.52, 1.32, 1.32, 1.18, 1.18, 1.08, 1.08, 1]
Before
Refine Attempt
refIndices = [1.52, 1.32, 1.32, 1.18, 1.18, 1.08, 1.08, 1]
Before
refIndices = [1.52, 1.41, 1.41, 1.32, 1.32, 1.24, 1.24, 1.18, 1.18, 1.13, 1.13, 1.08, 1.08, 1.03, 1.03, 1]
After
2x Resolution
Higher Resolution Brachistochrone
Expected
Observed
Results
Success!
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
def refractiveProfile(y):
profile = y**(-1/2)
return profile
for i in range (0, len(n)-1):
newLight = refractRay(n[i], n[i+1], v[i])
v.append(newLight)
Automated Refractive Indices Method
def refractiveProfile(y):
profile = y**(-1/2)
return profile
x = np.linspace(0.4,0.99,8)
autoIndices = []
for i in x:
autoIndices.append(refractiveProfile(i))
AutoIndices: [1.5811388300841895, 1.4369748623974747, 1.3261952985323262, 1.2376302619150348, 1.1647270698695131, 1.103354568734741, 1.050762078860976, 1.005037815259212]
Next:
- Round Indices
- If bounds are [a, b], make [a+0.01, b-0.01] to be safe
- Two of each value except first and last
Indices: [1.52, 1.41, 1.41, 1.32, 1.32, 1.24, 1.24, 1.18, 1.18, 1.13, 1.13, 1.08, 1.08, 1.03, 1.03, 1]
Automated Refractive Indices Method
def refractiveProfile(y):
profile = y**(-1/2)
return profile
x = np.linspace(0.4,0.99,8)
indicesOne = []
for i in x:
indicesOne.append(refractiveProfile(i))
indicesFinal = (np.round(np.sort(np.asarray(indicesOne * 2)),2)[::-1])[1:-1]
print("indicesFinal: ", indicesFinal)
indicesFinal: [1.58 1.58 1.44 1.44 1.33 1.33 1.24 1.24 1.16 1.16 1.1 1.1 1.05 1.05 1.01 1.01]
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
✅
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
✅
✅
Modules
Recursive Refraction
Refractive Profile Function
Refractive Indices Method
✅
✅
✅
Checks
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
Checks
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
Checks
✅
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
Checks
✅
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
def refractiveProfile(y):
profile = math.ceil(y)/10
return profile
Checks
✅
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
✅
Checks
✅
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
✅
The Ultimate Test: n = n(y)
The Ultimate Test: n = n(y)
The Ultimate Test: n = n(y)
✅
The Ultimate Test: n = n(y)
Refraction for Step Function
No Refraction If n(y)=c
✅
✅
The Ultimate Test: n = n(y)
🎉
🎉
Wohoo!
Final Improvements
➡️ Automatically compute boundaries for y.pos
Final Improvements
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
Final Improvements
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
Final Improvements
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Throw error if user resolution > permitted resolution
Final Improvements
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Throw error if user resolution > permitted resolution
Final Improvements
➡️ Automatically compute boundaries for y.pos
➡️ Automatically compute boundaries for y.pos
y = Symbol('y')
profile = Eq(y**(-1/2), 1.6)
lowerBound = float(round(solve(profile)[0],2))
def refractiveProfile(y):
profile = y**(-1/2)
return profile
x = np.linspace(lowerBound+0.01,0.99,25)
➡️ Automatically compute boundaries for y.pos
y = Symbol('y')
profile = Eq(y**(-1/2), 1.6)
lowerBound = float(round(solve(profile)[0],2))
def refractiveProfile(y):
profile = y**(-1/2)
return profile
x = np.linspace(lowerBound+0.01,0.99,25)
[ 0.39013671875 , 1 ]
Success!
➡️ Automatically compute boundaries for y.pos
y = Symbol('y')
profile = Eq(y**(-1/2), 1.6)
lowerBound = float(round(solve(profile)[0],2))
def refractiveProfile(y):
profile = y**(-1/2)
return profile
x = np.linspace(lowerBound+0.01,0.99,25)
[ 0.39013671875 , 1 ]
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Throw error if user resolution > permitted resolution
Final Improvements
➡️ Automatically compute boundaries for y.pos
Final Improvements
y = Symbol('y')
profile = Eq(eval(i), 1.6)
lowerBound = float(round(solve(profile)[0],2))
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Throw error if user resolution > permitted resolution
✅
Final Improvements
➡️ Enable user input for n(y)
i = input("fn: ")
def refractiveProfile(y):
return eval(i, {'y': y, 'np': np})
indicesOne = []
for k in x:
indicesOne.append(refractiveProfile(k))
x = np.linspace(lowerBound+0.01,0.99,25)
print("[",lowerBound, ",",1,"]")
indicesFinal = (np.round(np.sort(np.asarray(indicesOne)),2)[::-1])[1:-1]
light = createRay(vec(0,0,0), 41)
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Throw error if user resolution > permitted resolution
✅
Final Improvements
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Throw error if user resolution > permitted resolution
✅
Final Improvements
✅
➡️ Animate light path
for i in range (0, len(indicesFinal)-1):
t=0
while t < 1:
rate(5)
newLight = refractRay(indicesFinal[i], indicesFinal[i+1], v[i])
v.append(newLight)
horSum = 0
verSum = 0
for i in range(0,len(v)):
horSum += (v[i].pos.x+(v[i].axis.x)/2)
verSum += (v[i].pos.y+(v[i].axis.y/2))
horAvg = horSum / len(v)
verAvg = verSum / len(v)
hi = vec(horAvg,verAvg,0)
scene.center = hi
t+=1
➡️ Animate light path
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
✅
Final Improvements
✅
➡️ Automatically assign resolution as max resolution mathematically permitted
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Automatically assign resolution as max resolution mathematically permitted
✅
Final Improvements
✅
✅
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Automatically assign resolution as max resolution mathematically permitted
maximumResolution = 0
for i in range(0, len(indicesFinal)-1):
if (indicesFinal[i+1]/indicesFinal[i]) == 1:
maximumResolution = i
for i in range (0, maximumResolution):
t=0
while t < 1:
rate(5)
newLight = refractRay(indicesFinal[i], indicesFinal[i+1], v[i])
v.append(newLight)
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Automatically assign resolution as max resolution mathematically permitted
✅
Final Improvements
✅
✅
➡️ Enable user input for n(y)
➡️ Automatically compute boundaries for y.pos
➡️ Auto-detect max resolution by calculating math domain error prior to simulation
➡️ Animate light path
➡️ Automatically assign resolution as max resolution mathematically permitted
✅
Final Improvements
✅
✅
✅
✅
Code (Thus Far)
from vpython import *
from scipy.optimize import fsolve
import math
import numpy as np
from sympy import Eq, Symbol, solve
scene = canvas()
def createRay(start, angle, size=2):
v1 = vec(0,-1,0)
return arrow(pos = start, axis = size * v1.rotate(radians(angle)),
shaftwidth = 0.2, headwidth = 0.3, headlength = 0.2, color = color.yellow)
def refractRay(refIndex1, refIndex2, vector_inc):
v1 = vec(0,-1,0)
theta_inc = diff_angle(v1,vector_inc.axis)
theta_ref = asin((refIndex1/refIndex2) * sin(theta_inc))
return arrow(pos = vector_inc.pos+vector_inc.axis, axis = v1.rotate(theta_ref),
shaftwidth = 0.2, headwidth = 0.3, headlength = 0.2, color = color.yellow)
# Accepts User Input for the Refractive Profile
i = input("n(y) = ")
y = Symbol('y')
medium = i
if i == "ceil(y)":
print("Match Found")
lowerBound = 99
upperBound = 159
print("Bounds: [",lowerBound, ",",upperBound,"]")
elif "y" not in i:
print("Constant Function Not Accepted")
i = input("n(y) = ")
y = Symbol('y')
medium = i
else:
print("Match Found")
lowerProfile = Eq(eval(medium), 1)
upperProfile = Eq(eval(medium), 1.6)
bounds = [abs(round(float(solve(lowerProfile)[0]),3)),
abs(round(float(solve(upperProfile)[0]),3))]
lowerBound = np.amin(bounds)
upperBound = np.amax(bounds)
print("Bounds: [",lowerBound, ",",upperBound,"]")
resolution = 50
def refractiveProfile(y):
if i == "ceil(y)":
return math.ceil(y)/100
elif "y" not in i:
return int(i)
else:
return eval(i, {'y': y, 'np': np})
def d_fun(x):
h = 1e-5
return (refractiveProfile(x+h)-refractiveProfile(x-h))/(2*h)
currentResolution = 50
x = np.linspace(lowerBound+0.01,upperBound-0.01, currentResolution)
change = d_fun(np.random.uniform(lowerBound,upperBound))
indices = []
for k in x:
indices.append(refractiveProfile(k))
if change < 0:
print("Function is decreasing. Order must be reversed. [::-1]")
indicesFinal = (np.round(np.sort(np.asarray(indices)),2)[::-1])[1:-1]
else:
print("Function is increasing. Preserve order.")
indicesFinal = (np.round(np.sort(np.asarray(indices)),2))[1:-1]
light = createRay(vec(0,0,0), 35)
v = [light]
# Calculates maximum resolution
maximumResolution = 0
if i == "ceil(y)":
maximumResolution = 47
else:
for i in range(0, len(indicesFinal)-2):
if (indicesFinal[i+2] == indicesFinal[i]):
maximumResolution = i
if maximumResolution == 0:
maximumResolution = currentResolution
print("The maximum resolution is ", maximumResolution)
for i in range (0, maximumResolution-3):
t=0
while t < 1:
rate(5)
sphere()
newLight = refractRay(indicesFinal[i], indicesFinal[i+1], v[i])
v.append(newLight)
# Updates center of screen to be in the middle of the rays
horSum = 0
verSum = 0
for i in range(0,len(v)):
horSum += (v[i].pos.x+(v[i].axis.x)/2)
verSum += (v[i].pos.y+(v[i].axis.y/2))
horAvg = horSum / len(v)
verAvg = verSum / len(v)
hi = vec(horAvg,verAvg,0)
scene.center = hi
t+=1
➡️
➡️
➡️
➡️
Final Checks
➡️
➡️
➡️
➡️
Final Checks
Expectation
Expectation
Reality
Comparison
Comparison
Different initial conditions
Comparison
Verdict: Correct!
Different initial conditions
➡️
➡️
➡️
➡️
Final Checks
Verdict: Correct!
➡️
➡️
➡️
➡️
Final Checks
Verdict: Correct!
Expectation
Expectation
Expectation
Reality
✅
✅
✅
✅
➡️
➡️
➡️
➡️
Final Checks
The Brachistochrone Problem
By Refath Bari
The Brachistochrone Problem
- 24