In [ ]:
In [1]:
from vpython import *
import numpy
scene = canvas()
#Classes
# A class is a template for an object
#Ray class is a template for a ray
#init: Set all of the variables for the beam
#define methods
class ray():
'''A light ray class'''
def __init__(self,startPoint,directionVector,color):
self.startPoint = startPoint
self.oldPoint = startPoint ### Position of light ray prior to current time step ###
self.currentPoint = startPoint
self.directionVector = directionVector.norm()
self.color = color
self.length = 10
self.world_n=1 ### Value of n for medium ray is currently in ###
self.ray = curve(pos=self.currentPoint, color=self.color)
#self.draw_ray() ### Commented out to let world() draw ray
#Draw all rays in Source
def draw_ray(self):
'''Draw a single light ray'''
self.oldPoint=self.currentPoint ### Set old point to the current position BEFORE the update ###
newPoint = self.currentPoint + self.length*self.directionVector
self.currentPoint = newPoint
self.ray.append(newPoint)
#Loop to check the position of newPoint. If position is greater than 5, call Snell's Law
# and assign a new direction vector
#Define normal as normal vector to object
#Using Snell's Law (numpy.arcsin((1/(1.5))*sin(vector.diff_angle(ray1,normal))))
def newDirection(self,newDirection):
'''Set a new direction for the ray'''
self.directionVector = newDirection
#Stop when encounters an object and then reset length to 10
def newLength(self,newLength):
'''Set a new length for the vector'''
self.length=newLength
class beam():
'''A Beam Class'''
def __init__(self, position, direction, width):
self.direction = direction
self.position = position
self.centerRay = ray(position, direction, color.green)
self.currentPoint = self.centerRay.currentPoint
self.beamDirection = direction
self.width = width
self.color = self.centerRay.color
self.Rays = []
### Begin changes to fix perpendicular rays ###
self.find_perp() #Find two rays perpendicular to centerRay
self.draw_beam()
def find_perp(self):
'''Find two vectors perpendicular to the center ray'''
if self.centerRay.directionVector.y==0 and self.centerRay.directionVector.x==0:
self.perp1=vec(1,0,0)
else:
self.perp1 = norm(vec(self.centerRay.directionVector.y,-self.centerRay.directionVector.x,0))
self.perp2 = norm(vec(self.centerRay.directionVector.cross(self.perp1)))
### End changes to fix perpendicular rays ###
def draw_beam(self):
'''Draw a beam; lots of parallel rays'''
for i in numpy.arange(0,self.width,.5):
for k in numpy.arange(0,2*pi+0.1,.5):
self.newStart = self.centerRay.startPoint + i*(self.perp1*cos(k)+self.perp2*sin(k))
beamRay = ray(self.newStart, self.beamDirection, self.color)
self.Rays.append(beamRay)
class pointSource():
def __init__(self, position):
self.position = position
self.centerRay = ray(vec(self.position), vec(1,0,0), color.green)
self.color = self.centerRay.color
self.Rays = []
self.source = sphere(pos=vec(self.centerRay.startPoint), radius=0.01, color=self.color)
self.draw_pointSource()
def draw_pointSource(self):
'''Draw a point source'''
self.source
for i in numpy.arange(0,2*pi+0.1,(pi/5)):
for k in numpy.arange(0,2*pi+0.1,(pi/5)):
x = self.centerRay.length*(cos(k)*cos(i))
y = self.centerRay.length*(cos(k)*sin(i))
z = self.centerRay.length*(sin(k))
pointSourceRay = ray(self.position, vec(x,y,z), self.centerRay.color)
self.Rays.append(pointSourceRay)
#####
# Start of World
#####
class world():
def __init__(self):
self.rays=[]
self.objects=[]
self.dL = 0.1 #Ray step size
self.n = 1 #Default index of refraction for the world
self.MAX_LENGTH = 20 ### Set to 10 ###
self.current_length = 0 ### Length of current ray being drawn ###
def add_ray(self,new_ray):
'''Add a ray to the world'''
self.rays.append(new_ray)
new_ray.length=self.dL
new_ray.world_n = self.n
def add_object(self,new_object):
'''Add a new object to the world'''
self.objects.append(new_object)
def add_light(self,light):
'''Add a point source to the world'''
for i in light.Rays:
self.add_ray(i)
def add_beam(self, beam):
'''Add a beam to the world'''
for i in beam.Rays:
self.add_ray(i)
def draw_rays(self):
#Loop for ceching Boundaries
for i in self.rays:
self.current_length=0
while self.current_length < self.MAX_LENGTH: ### Runs until we reach a max length
for j in self.objects: ### Loop through all objects ###
j.check_boundaries(i) ### Call check_boundaries method for each object ###
i.draw_ray()
self.current_length=self.current_length+self.dL
#Loop
#def draw_objets(self):
#######
#ball class
######
class ball():
def __init__(self,position,radius):
self.pos = position
self.r = radius
self.color=color.red
self.n = 1.5 #Let's make it out of glass for now
self.shape = 'sphere'
self.sphere = sphere(pos=self.pos, radius=self.r, color=self.color, opacity=0.5)
def check_boundaries(self,ray):
surface_norm = norm(ray.currentPoint-self.pos)
rel_pos = self.pos - ray.currentPoint
rel_old_pos=self.pos - ray.oldPoint
### The following code should be cleaned up and better commented by me ###
if abs(rel_pos.dot(surface_norm))< self.r and abs(rel_old_pos.dot(surface_norm))> self.r: ### Going into sphere ###
thetaIncident = pi-diff_angle(ray.directionVector,surface_norm)
thetaRefracted = asin((ray.world_n/self.n)*sin(thetaIncident))
a1 = ray.directionVector - surface_norm * (dot(ray.directionVector,surface_norm))
ray.directionVector = (-1*cos(thetaRefracted)*surface_norm) + (sin(thetaRefracted)*a1.norm())
if abs(rel_pos.dot(surface_norm))> self.r and abs(rel_old_pos.dot(surface_norm))< self.r: ### Going out of sphere ###
thetaIncident = pi-diff_angle(ray.directionVector,surface_norm)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
thetaRefracted = thetaIncident
a1 = ray.directionVector - surface_norm * (dot(ray.directionVector,surface_norm))
ray.directionVector = (-1*cos(thetaRefracted)*surface_norm) + (sin(thetaRefracted)*a1.norm())
else:
thetaRefracted = asin((self.n/ray.world_n)*sin(thetaIncident))
a1 = ray.directionVector - surface_norm * (dot(ray.directionVector,surface_norm))
ray.directionVector = (1*cos(thetaRefracted)*surface_norm) + (sin(thetaRefracted)*a1.norm())
class block():
'''an optical obstacle class'''
def __init__(self, position, axis, up, length, width, height, material, opacity):
self.position = position
self.axis = axis
self.up = up
self.material = material
self.length = length
self.width = width
self.height = height
if(self.material == "glass"):
self.n = 1.5
self.opacity = opacity
self.draw_block()
def draw_block(self):
self.shape = box(pos = self.position, axis = self.axis, up = self.up, length = self.length, width = self.width, height = self.height, opacity = self.opacity)
def check_boundaries(self, ray):
self.normalFront = self.axis.norm()
self.normalBack = -1 * self.normalFront
self.normalTop = self.up.norm()
self.normalBottom = -1 * self.normalTop
self.normalLSide = cross(self.axis.norm(), self.up.norm())
self.normalRSide = -1 * self.normalLSide
b = ray.currentPoint - self.position
b_old = ray.oldPoint - self.position
#If ALL of these 'check' statements are true, then the ray is within the block
frontCheck = bool(abs(dot(b,self.normalFront)) - ((self.length)/2) <= 0)
backCheck = bool(abs(dot(b,self.normalBack)) - ((self.length)/2) <= 0)
topCheck = bool(abs(dot(b,self.normalTop)) - ((self.height)/2) <= 0)
bottomCheck = bool(abs(dot(b,self.normalBottom)) - ((self.height)/2) <= 0)
lSideCheck = bool(abs(dot(b,self.normalLSide)) - ((self.width)/2) <= 0)
rSideCheck = bool(abs(dot(b,self.normalRSide)) - ((self.width)/2) <= 0)
#If ALL of these are true, then the ray was within the block on the last time step
frontCheck_old = bool(abs(dot(b_old,self.normalFront)) - ((self.length)/2) <= 0)
backCheck_old = bool(abs(dot(b_old,self.normalBack)) - ((self.length)/2) <= 0)
topCheck_old = bool(abs(dot(b_old,self.normalTop)) - ((self.height)/2) <= 0)
bottomCheck_old = bool(abs(dot(b_old,self.normalBottom)) - ((self.height)/2) <= 0)
lSideCheck_old = bool(abs(dot(b_old,self.normalLSide)) - ((self.width)/2) <= 0)
rSideCheck_old = bool(abs(dot(b_old,self.normalRSide)) - ((self.width)/2) <= 0)
def refractIn (face):
thetaRefracted = asin((ray.world_n/self.n)*sin(thetaIncident))
a1 = ray.directionVector - self.normalFront * (dot(ray.directionVector,self.normalFront))
ray.directionVector = (1*cos(thetaRefracted)*self.normalFront) + (sin(thetaRefracted)*a1.norm())
def reflectIn (face):
thetaReflected = thetaIncident
a1 = ray.directionVector - (self.normalFront * dot(ray.directionVector,self.normalFront))
ray.directionVector = (-1*cos(thetaReflected)*self.normalFront) + (sin(thetaReflected)*a1.norm())
def refractOut(face):
thetaRefracted = asin((self.n/ray.world_n)*sin(thetaIncident))
a1 = ray.directionVector - (face * dot(ray.directionVector,face))
ray.directionVector = (1*cos(thetaRefracted)*face) + (sin(thetaRefracted)*a1.norm())
def reflectOut(face):
thetaReflected = thetaIncident
a1 = ray.directionVector - (self.normalFront * dot(ray.directionVector,self.normalFront))
ray.directionVector = (-1*cos(thetaReflected)*self.normalFront) + (sin(thetaReflected)*a1.norm())
###Check if ray is exiting the block###
if((frontCheck_old and backCheck_old and topCheck_old and bottomCheck_old and lSideCheck_old and rSideCheck_old) and ((not frontCheck and frontCheck_old) or (not backCheck and backCheck_old) or (not topCheck and topCheck_old) or (not bottomCheck and bottomCheck_old) or (not lSideCheck and lSideCheck_old) or (not rSideCheck and rSideCheck_old))):
if(frontCheck_old and not frontCheck): #Exiting Block through Axis Direction
if(b.dot(self.normalBack) < b.dot(self.normalFront)): #Exiting through Front Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalBack)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalFront)
else:
refractOut(self.normalFront)
if(b.dot(self.normalFront) < b.dot(self.normalBack)): #Exiting through Back Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalFront)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalBack)
else:
refractOut(self.normalBack)
if(topCheck_old and not topCheck): #Exiting Block through Up Direction
if(b.dot(self.normalBottom) < b.dot(self.normalTop)): #Exiting through Top Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalBottom)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalTop)
else:
refractOut(self.normalTop)
if(b.dot(self.normalTop) < b.dot(self.normalBottom)): #Exiting through Bottom Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalTop)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalBottom)
else:
refractOut(self.normalBottom)
if(lSideCheck_old and not lSideCheck): #Exiting Block through Cross Direction
if(b.dot(self.normalRSide) < b.dot(self.normalLSide)): #Exiting through Left Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalRSide)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalLSide)
else:
refractOut(self.normalLSide)
if(b.dot(self.normalLSide) < b.dot(self.normalRSide)): #Exiting through Right Face
thetaIncident = pi-diff_angle(ray.directionVector,self.normalLSide)
if((self.n/ray.world_n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectOut(self.normalRSide)
else:
refractOut(self.normalRSide)
###Check if ray is entering the block###
if((frontCheck and backCheck and topCheck and bottomCheck and lSideCheck and rSideCheck) and (frontCheck_old or backCheck_old or topCheck_old or bottomCheck_old or lSideCheck_old or rSideCheck_old)):
if(frontCheck and not frontCheck_old): #Entering through the axis direction
if(b.dot(self.normalFront) < b.dot(self.normalBack)): #Entered through front
thetaIncident = pi-diff_angle(ray.directionVector, self.normalFront)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalFront)
else:
refractIn(self.normalFront)
if(b.dot(self.normalBack) < b.dot(self.normalFront)): #Entered through Back
thetaIncident = pi-diff_angle(ray.directionVector, self.normalBack)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalBack)
else:
refractIn(self.normalBack)
if(topCheck and not topCheck_old): #Entering through the up direction
if(b.dot(self.normalTop) < b.dot(self.normalBottom)): #Entering through Top
thetaIncident = pi-diff_angle(ray.directionVector, self.normalTop)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalTop)
else:
refractIn(self.normalTop)
if(b.dot(self.normalBottom) < b.dot(self.normalTop)): #Entering through Bottom
thetaIncident = pi-diff_angle(ray.directionVector, self.normalBottom)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalBottom)
else:
refractIn(self.normalBottom)
if(lSideCheck and not lSideCheck_old): #Entering through cross direction
if(b.dot(self.normalLSide) < b.dot(self.normalRSide)): #Entering through Left Side
thetaIncident = pi-diff_angle(ray.directionVector, self.normalLSide)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalLSide)
else:
refractIn(self.normalLSide)
if(b.dot(self.normalRSide) < b.dot(self.normalLSide)): #Entering through Right Side
thetaIncident = pi-diff_angle(ray.directionVector, self.normalRSide)
if((ray.world_n/self.n)*sin(thetaIncident) > 1): #Argument of arcsin > 1 --> Total Internal Reflection
reflectIn(self.normalRSide)
else:
refractIn(self.normalRSide)
#class lens():
# def __init__ (self, position, r1, r2)
# self.position = position #position of the center of the lens?
# self.radius1 = r1
# self.radius2 = r2
#######
# Test Code
#######
round_thing=ball(vec(5.1,0,0), 5)
#square_thing=block(vec(4,0,0), vec(1,0,0), vec(0,1,0), 1, 20, 20, "glass", opacity = 0.5)
#ps1 = pointSource(vec(0,0,0))
beam1 = beam(vec(-3,0,0), vec(1,0,0), 4)
#light1=ray(vec(-5,1,0),vec(1,0,0),color.blue)
#light2=ray(vec(-5,-1,0),vec(1,0,0), color.blue)
#light3=ray(vec(-5,0,0), vec(1,0,0), color=color.yellow)
#light4=ray(vec(-5,2,0), vec(1,0,0), color.red)
#light5=ray(vec(-5,-2,0), vec(1,0,0), color.red)
#light6=ray(vec(-5,0,-2),vec(1,0,0),color.red)
#light7=ray(vec(-5,0,-1),vec(1,0,0), color.blue)
#light8=ray(vec(-5,0,0), vec(1,0,0), color=color.yellow)
#light9=ray(vec(-5,0,1), vec(1,0,0), color.blue)
#light10=ray(vec(-5,0,2), vec(1,0,0), color.red)
my_world=world()
#my_world.add_ray(light1)
#my_world.add_ray(light2)
#my_world.add_ray(light3)
#my_world.add_ray(light4)
#my_world.add_ray(light5)
#my_world.add_ray(light6)
#my_world.add_ray(light7)
#my_world.add_ray(light8)
#my_world.add_ray(light9)
#my_world.add_ray(light10)
my_world.add_beam(beam1)
#my_world.add_light(ps1)
my_world.add_object(round_thing)
my_world.draw_rays()
In [ ]:
In [ ]: