代码主体

# -*- coding: utf-8 -*-
"""
Created on Fri Feb 24 09:30:12 2023

@author: Conway
"""

'''
@Default unit: mm

'''
#%% Diagram of the optical path
'''           z3
    \      \  :
*   \      \  :
    \      \  :
    z1     z2
    
    
    
*   light source
\   surface
:   micron structure
'''
#%% import libraries
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
matplotlib.use('Qt5Agg')
import sys
sys.path.append(r'D:\Python program\Snells law')
import snell
#%% parameters
'''
z1: the front surface of substrate
z2: the rear surface of substrate
z3: the micron structrue vertex z-axis coordinate
n3  the refractive index between surface 2&3 (unable)
n2: the refractive index between surface 1&2
n1: the refractive index in front of surface 1
thickness: thickness of substrate
'''

thickness=0.25#constant
z1=2 #constant
z2=z1+thickness #iter
z3=z2+0.08  #iter

n1=1
n2=1.5
# n3=n2
#%%
'''
theta1: incident angle of substrate front surface
theta2: emergent angle of substrate front surface
h2: corresponding to z2, denoting the y coordinate
max_step: maxinum depth of the micron structure
shift: y-shift of (z2,h2)
cvector: denoting the emergent light direction on substrate rear surface
'''
theta2=np.arange(0,np.pi/180*41.5,1e-5)
theta1=np.arcsin(n2/n1*np.sin(theta2))

h2=0

y=np.array([])
z=y

max_step=0.08
shift=max_step/2
cvector=[0,1]
#%% parameters
'''
Radius: Radius of substrate
spread_degree: spread degree of the designed collimated light
center: focus coordinates of the designed collimated light
'''
radius=2
spread_degree=6
center=[0,z2-radius/np.tan(spread_degree/180*np.pi)] # collimating light focus
#%%
while h2<=radius:
    temporary=np.array([float(h2),thickness+z1])-np.array(center)
    angle=np.arctan2(temporary[0],temporary[1])
    cvector=[float(np.sin(angle)),float(np.cos(angle))]
    
    value=np.abs(z1*np.tan(theta1)+(z2-z1)*np.tan(theta2)-(h2+shift))
    index=np.where(value==np.amin(value))
    index=index[0]
    theta2_current=theta2[index]
    
    y=np.hstack([y,h2])
    z=np.hstack([z,z2])
    if z2>thickness+z1:
        z2=thickness+z1
    
    h3=(z3-z2)*np.tan(theta2_current/100*50)+h2
    
    y=np.hstack([y,h3])
    z=np.hstack([z,z3])
    
    vector=[float(np.sin(theta2_current)),float(np.cos(theta2_current))]
    norm=snell.refract(vector, cvector, n2, n1)
    orth1,orth2=snell.orthvec(norm)
    if orth1[0]>0:
        orth=orth1
    else:
        orth=orth2
    
    alpha=np.arctan2(orth[0],orth[1])
    beta=np.pi-alpha
    
    h2_step=(z3-(thickness+z1))*np.tan(beta)
    h2_last=h2
    
    if h2_step>max_step:
        h2_step=max_step
        h2=h2_step+h3
        z2=z3-h2_step/np.tan(beta)
    else:
        z2=(thickness+z1)
        h2_step=(z3-z2)*np.tan(beta)
        h2=h2_step+h3
    shift=np.abs(h2-h2_last)/2
    
#%%
z[0]=z[1]
z=np.hstack([z,z1])
y=np.hstack([y,y[np.size(y)-1]])
curve=np.stack([z,y],axis=1)
#%%
fig1=plt.figure()
plt.axis('equal')
plt.plot(z,y)
plt.show(block=True)

snell.py代码

# -*- coding: utf-8 -*-
"""
Created on Thu Feb  9 09:49:10 2023

@author: Conway

"""
'''
-------------------------------------------------------
Defination
-------------------------------------------------------

Define the coordinates format:
    vector=[row,column]=[y,x]
Define the positive angle direction: 
    anti-clock-wise

-------------------------------------------------------
'''
import numpy as np
# import matplotlib
# matplotlib.use('Qt5Agg')
# import matplotlib.pyplot as plt
#%%
'''
--------------------------------------------------------------------------------------------------------------
@Function orthvec:
    return the orthogonal vector corresponding to incident vector and reflected vector 
    
@input format: incident:[y,x],reflected:[y,x]

@Positive direction: 
    Agnle: Anti-Clock-wise
    Normal vector: from reflection surface to the free space
-------------------------------------------------------
'''
def orthvec (vector): #input format (including emergent&incident):[y,x]
    # Normalizing the vector
    vector=vector/np.sqrt(np.sum(np.square(vector)))
    angle=np.arctan2(vector[0],vector[1])
    
    orthognal_angle=np.array([0,0])
    orthognal_angle=np.array([angle,angle])+np.array([np.pi/2,-np.pi/2])
    
    orth_vec0=np.array([np.sin(orthognal_angle[0]),np.cos(orthognal_angle[0])])
    orth_vec1=np.array([np.sin(orthognal_angle[1]),np.cos(orthognal_angle[1])])
    return orth_vec0,orth_vec1
'''
--------------------------------------------------------------------------------------------------------------
'''
#%%
'''
--------------------------------------------------------------------------------------------------------------
@Function refract0:
    return the emergent vector corresponding to incident vector and normal vector 
    
@input format: incident:[y,x],emergent:[y,x]

@Positive direction: 
    Agnle: Anti-Clock-wise
    Normal vector: from the optically thinner medium to optically denser medium
--------------------------------------------------------------------------------------------------------------
'''
def refract0 (incident,norm,ni,ne): #input format (including emergent&incident):[y,x]
    # Normalizing the vector
    incident=incident/np.sqrt(np.sum(np.square(incident)))
    norm=norm/np.sqrt(np.sum(np.square(norm)))
    
    inc=incident*ni

    # convert to float
    ni=np.float64(ni)
    ne=np.float64(ne)
    # refraction law
    eme=norm+inc
    if np.sum(np.square(norm-inc))==0:
        norm=incident
    else:
        eme=eme/np.sqrt(np.sum(np.square(eme))) # normalize the vector
    # vector norm positive direction is from optically thinner medium to optically denser medium
    return eme
'''
--------------------------------------------------------------------------------------------------------------
'''

#%%
'''
--------------------------------------------------------------------------------------------------------------
@Function refract:
    return the normal vector corresponding to incident vector and emergent 
    
@input format: incident:[y,x],emergent:[y,x]

@Positive direction: 
    Agnle: Anti-Clock-wise
    Normal vector: from the optically thinner medium to optically denser medium
--------------------------------------------------------------------------------------------------------------
'''
def refract (incident,emergent,ni,ne): #input format (including emergent&incident):[y,x]
    # Normalizing the vector
    incident=incident/np.sqrt(np.sum(np.square(incident)))
    emergent=emergent/np.sqrt(np.sum(np.square(emergent)))
    
    inc=incident*ni
    eme=emergent*ne
    if np.dot(inc,eme)<=0:
        raise Exception('The angle betweent emergent wave vector'\
                        'and incident wave vector must be less than 90 degree')
        return
    # convert to float
    ni=np.float64(ni)
    ne=np.float64(ne)
    # refraction law
    norm=eme-inc
    if np.sum(np.square(eme-inc))==0:
        norm=incident
    else:
        norm=norm/np.sqrt(np.sum(np.square(norm))) # normalize the vector
    # vector norm positive direction is from optically thinner medium to optically denser medium
    return norm
'''
--------------------------------------------------------------------------------------------------------------
'''

#%%
'''
--------------------------------------------------------------------------------------------------------------
@Function reflect:
    return the normal vector corresponding to incident vector and reflected vector 
    
@input format: incident:[y,x],reflected:[y,x]

@Positive direction: 
    Agnle: Anti-Clock-wise
    Normal vector: from reflection surface to the free space
-------------------------------------------------------
'''
def reflect (incident,reflected): #input format (including emergent&incident):[y,x]
    # Normalizing the vector
    incident=incident/np.sqrt(np.sum(np.square(incident)))
    reflected=reflected/np.sqrt(np.sum(np.square(reflected)))

    norm=reflected-incident
    
    if norm==0:
        raise Exception('Reflected wave vector and incident wave vector'\
                        ' can\'t be collinear')
        return
    
    norm=norm/np.sqrt(np.sum(np.square(norm))) # normalize the vector
    
    return norm
'''
--------------------------------------------------------------------------------------------------------------
'''