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

n_1\sin(\theta_1)=n_2\sin(\theta_2)

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

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

1≤n(y)≤1.6

n(y)=y^{-1/2}

1≤y^{-1/2}≤1.6

0.39≤y≤1.6

Plan

1≤n(y)≤1.6

n(y)=y^{-1/2}

1≤y^{-1/2}≤1.6

0.39≤y≤1

-7≤y.pos≤-1

1≤|y.pos|≤7

\frac{7}{x}=1 \\ \frac{7}{1}=x \\ x=7

\frac{1}{7}≤\frac{|y.pos|}{7}≤\frac{7}{7}

0.14≤\frac{|y.pos|}{7}≤1

Plan

1≤n(y)≤1.6

n(y)=y^{-1/2}

1≤y^{-1/2}≤1.6

0.39≤y≤1

-7≤y.pos≤-1

1≤|y.pos|≤7

\frac{7}{x}=1 \\ \frac{7}{1}=x \\ x=7

\frac{1}{7}≤\frac{|y.pos|}{7}≤\frac{7}{7}

0.14≤\frac{|y.pos|}{7}≤1

n(y)=y^{-1/2}

y=0.39

y=1

n(y)=1

n(y)=1.6

|y.pos|

|y.pos| = 1

|y.pos| = 7

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=a_1

\frac{|y.pos|}{k}=a_2

Plan

n(y)=y^{-1/2}

y=0.39

y=1

n(y)=1

n(y)=1.6

|y.pos|

|y.pos| = 1

|y.pos| = 7

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=0.14

\frac{|y.pos|}{k}=1

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=0.14

\frac{|y.pos|}{k}=1

y=0.39

y=1

Plan

n(y)=y^{-1/2}

y=0.39

y=1

n(y)=1

n(y)=1.6

|y.pos|

|y.pos| = 1

|y.pos| = 7

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=0.14

\frac{|y.pos|}{k}=1

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=0.14

\frac{|y.pos|}{k}=1

y=0.39

y=1

Plan

\frac{|y.pos|}{k}

\frac{|y.pos|}{k}=0.14

\frac{|y.pos|}{k}=1

y=0.39

y=1

n(y)=y^{-1/2}

n(y)=1

n(y)=1.6

Plan

y=-1

y=-2

y=-3

y=-4

y=-5

y=-6

y=-7

n(y)=y^{-1/2}

Plan

y=-1

y=-2

y=-3

y=-4

y=-5

y=-6

y=-7

n(y)=y^{-1/2}

y=\frac{|-1|}{7}

k=7

y=\frac{|-2|}{7}

y=\frac{|-3|}{7}

y=\frac{|-4|}{7}

y=\frac{|-5|}{7}

y=\frac{|-6|}{7}

y=\frac{|-7|}{7}

y=0.14

y=0.28

y=0.42

y=0.57

y=0.71

y=0.85

y=1

n(y)=2.67

n(y)=1.88

n(y)=1.52

n(y)=1.32

n(y)=1.18

n(y)=1.08

n(y)=1

y<0.39

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!

y=-3

y=-4

y=-5

y=-6

y=-7

y=\frac{|-3|}{7}

y=\frac{|-4|}{7}

y=\frac{|-5|}{7}

y=\frac{|-6|}{7}

y=\frac{|-7|}{7}

y=0.42

y=0.57

y=0.71

y=0.85

y=1

n(y)=1.52

n(y)=1.32

n(y)=1.18

n(y)=1.08

n(y)=1

Refractive Indices Expectation

y=-3

y=-4

y=-5

y=-6

y=-7

y=0.42

y=0.57

y=0.71

y=0.85

y=1

n(y)=1.52

n(y)=1.32

n(y)=1.18

n(y)=1.08

n(y)=1

\theta=30, n = 1.6

\theta=31.75, n = 1.52

\theta=37.29, n = 1.32

\theta=42.66, n = 1.18

\theta=47.76, n = 1.08

\theta=53.08, n = 1

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

✅

n(y)=y^{-1/2}

The Ultimate Test: n = n(y)

n(y)=y^{-1/2}

The Ultimate Test: n = n(y)

n(y)=y^{-1/3}

The Ultimate Test: n = n(y)

n(y)=y^{-1/3}

n(y)=y^{-1/2}

✅

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

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

n(y)=y^{-1/2}

n(y)=c

n(y)=\lfloor{y}\rfloor

n(y)=y^{-1/8}

➡️

Final Checks

n(y)=y^{-1/2}

n(y)=c

n(y)=\lfloor{y}\rfloor

n(y)=y^{-1/8}

Expectation

n(y)=y^{-1/2}

Expectation

Reality

n(y)=y^{-1/2}

Comparison

n(y)=y^{-1/2}

Comparison

Different initial conditions

n(y)=y^{-1/2}

Comparison

Verdict: Correct!

Different initial conditions

➡️

Final Checks

n(y)=y^{-1/2}

n(y)=c

n(y)=\lfloor{y}\rfloor

n(y)=y^{-1/8}

Verdict: Correct!

➡️

Final Checks

n(y)=y^{-1/2}

n(y)=c

n(y)=\lfloor{y}\rfloor

n(y)=y^{-1/8}

Verdict: Correct!

Expectation

n(y)=c

Expectation

n(y)=c

Expectation

Reality

✅

➡️

Final Checks

n(y)=y^{-1/2}

n(y)=c

n(y)=\lfloor{y}\rfloor

n(y)=y^{-1/8}