Creates a plot of ray height vs. spherochromatism for a equiconvex lens.

spherochromatism.py

#!/usr/bin/python -u
# -*- coding: utf-8 -*-

"""
Creates a plot of ray height vs. spherochromatism for a equiconvex lens.
The lens parameters have been hardcoded below. You can tweak them as required.

Conventions:
Angle made with the positive axis in couter clockwise fashion is positive.
All units, unless mentioned are in SI

License: BSD
"""

__author__ = 'पुष्पक दगड़े (Pushpak Dagade) <guanidene@gmail.com>'
__date__ = 'Thu Sep 25 07:30:18 2012'
__version__ = '1.0.0'

from numpy import arcsin, arange, sin, sqrt
import matplotlib.pyplot as plt
DEBUG = False

# Lens parameters
lambda_C = 656.3
lambda_F = 486.1
C = 1e-1
t = 2/C - 2*sqrt((1/C)**2 - 1)      # I want height of marginal ray = 1metre


def n_of_lambda(lambda_):
    """
    Returns 'n' for a given lambda' in nanometers.
    Uses Cauchy's formula and returns data for hard crown glass K5.
    The parameters 'A' and 'B' for the same are defined below.
    """
    # For hard crown glass K5
    A = 1.5220
    B = 0.00459  # micron**2
    return A + 1e6*B/(lambda_**2)


def Q_U_method(Q, U, n1, n2, C):
    """
    Returns the Q' and U' of the ray after refraction from surface of
    curvature C. Refractive indices of incident and refracted ray are n1 to n2.
    """

    sin_u = sin(U)
    sin_i = sin_u + Q*C
    sin_i_prime = n1/n2*sin_i
    sin_u_prime = sin(U - arcsin(sin_i) + arcsin(sin_i_prime))

    Q_prime = (sin_i_prime + sin_u_prime)/C
    U_prime = arcsin(sin_u_prime)
    if DEBUG:
        print '\nn1,n2:', n1, n2
        print 'Q ,U :', Q, U
        print "Q',U':", Q_prime, U_prime

    return Q_prime, U_prime


def transfer_method(Q1_prime, U1_prime, t):
    """
    Transfers Q', U' obtained from one surface to Q2, U2 for other surface
    depending upon the distance between them.

    If U1_prime is negative, these would hold:
        Q2 < Q1_prime
        Q2 > 0, Q1 > 0

    XXX. This has not been verified thoroughly
    """

    Q2 = Q1_prime + t*sin(U1_prime)
    return Q2, U1_prime         # U2 = U1_prime


def spherical_aberration(n1, n2, C, h, t):
    """
    Return the lateral distance between the focal points of marginal
    and paraxial parallel incident rays.

    This function requires the refractive indices of medium and lens for the
    wavelength of incident ray, curvature of equiconvex lens and its thickness.
    """

    U = 0
    h_lower = 1e-10      # For paraxial ray
    L_primes = []

    for Q in [h_lower, h]:
        if DEBUG:
            print '\n','-'*50
        Q1_prime, U1_prime = Q_U_method(Q, U, n1, n2, C)
        Q2, U2 = transfer_method(Q1_prime, U1_prime, t)
        Q2_prime, U2_prime = Q_U_method(Q2, U2, n2, n1, -C)
        L_primes.append(-Q2_prime/sin(U2_prime))

    return L_primes[1] - L_primes[0]


def sphero_chromatism(nF1, nF2, nC1, nC2, C, h, t):
    """
    Spherochromatism: Difference between chromatic aberrations of the F and C
    light wavelengths.

    This function requires the medium and lens refractive indices for F and C
    light, curvature of equiconvex lens and its thickness.
    """
    return spherical_aberration(nF1, nF2, C, h, t) - \
           spherical_aberration(nC1, nC2, C, h, t)


if __name__ == '__main__':
    print 'Lens Parameters:'
    print 'n for F light:', n_of_lambda(lambda_F)
    print 'n for C light:', n_of_lambda(lambda_C)
    print 'Curvature (m):', C
    print 'Thickness (m):', t

    # Create the plot.
    H = arange(1e-10, 1 + 1e-10, 1e-2)
    sph_ch = [sphero_chromatism(1, n_of_lambda(lambda_F),
                                1, n_of_lambda(lambda_C),
                                C, h, t)/1e-3 for h in H]       # units in mm
    plt.plot(sph_ch, H)
    plt.grid()
    plt.xlabel('Spherochromatism (mm)')
    plt.ylabel('Height of ray from principal axis (m)')
    plt.show()