#!/usr/bin/python # # Same as sphere-focus-1 but the formula used is for light entering # the convex face and coming to a focus within the glass. # Too lazy today to teach the earlier script to do both, so we cut-and-paste it... # import math import re import sys output_filename = "sphere-focus-2.dat" Nl_over_Ne = 1.5 steps = 100 start_val = 0.00001 end_val = 1.0 debug = 0 def dbprint (str): if debug: print str ## dbprint class Whoops: def __init__(self,str = None, throw=None): if str: print str else: print "Whoops!" if throw or not str: raise self exit() # **************************************************************** # Find sine of xi def sin_xi (x_over_r): sqrts = math.sqrt(Nl_over_Ne**2 - x_over_r**2) - math.sqrt(1.0 - (x_over_r**2)) result = (x_over_r / Nl_over_Ne) * sqrts return result # **************************************************************** # Find tangent of xi def tan_xi (x_over_r): sxi = sin_xi (x_over_r) result = sxi / math.sqrt(1 - sxi**2) return result # **************************************************************** # Process arguments file_list_name = None have_ofile_name = 0 an = 1 while an < len(sys.argv): arg = sys.argv[an] an += 1 if arg == "-debug": debug = 1 elif arg == "-steps": steps = int(sys.argv[an]) an += 1 elif arg == "-index": Nl_over_Ne = float(sys.argv[an]) an += 1 elif arg == "-out_file": output_filename = sys.argv[an] have_ofile_name = 1 an += 1 else: Whoops ("Argument %s not understood" % arg) if not have_ofile_name: output_filename = "sphere-focus-2.N_%s" % (Nl_over_Ne) # **************************************************************** # Main program ofile = open (output_filename, "w") ofile.write ("# Sphere focus, light entering convex side, index=%s, steps=%s\n" % (Nl_over_Ne, steps)) delta = (end_val - start_val) / float(steps) had_val_err = 0 for step in xrange(steps): x = start_val + (step * delta) dbprint ("Step %s, x=%s" % (step, x)) try: txi = tan_xi (x) except ValueError: dbprint ("Value error raised at x=%s" % x) had_val_err = 1 break dbprint (" ... Tan xi = %s" % txi) if txi <= 0: dbprint ("Skipping step %s, x=%s -- tangent is zero" % (step, x)) else: f = x * (1.0 / txi) + 1 - math.sqrt(1 - x**2) dbprint (" ... f = %s" % f) ofile.write ("%s %s\n" % (x, f)) if had_val_err: ofile.write ("# Error encountered at x=%s --\n" + "# but we should not see cutoff in this orientation??\n" % x) ## ofile.write ("# Hit TIR at x=%s\n" % x) ofile.close()