diff options
Diffstat (limited to 'index.py')
| -rw-r--r-- | index.py | 33 |
1 files changed, 26 insertions, 7 deletions
@@ -37,28 +37,47 @@ def get_index(p, wavelength, T): return np.sqrt((2*c+1)/(1-c)) +def approximate_index(x, a, b): + """ + Returns the index of refraction using the approximation: + + n = a + b/x**2 + """ + return a + b/x**2 + +def chi2(pars, x, data): + """ + Returns the squares of the difference between the data and the model for + the inverse of the index of refraction, i.e. + + 1/n = lambda**2/(a*lambda**2+b) + """ + a, b = pars + return np.sum((approximate_index(x,a,b)-data)**2) + if __name__ == '__main__': import matplotlib.pyplot as plt import argparse + from scipy.optimize import fmin - parser = argparse.ArgumentParser("fit a Taylor series approximation to the inverse of the index of refraction of water") + parser = argparse.ArgumentParser("fit for an approximation to the inverse of the index of refraction of water") parser.add_argument("--temperature", type=float, default=10.0) parser.add_argument("--density", type=float, default=1000.0) - parser.add_argument("--deg", type=int, default=2) - parser.add_argument("--x0", type=float, default=400.0) parser.add_argument("--low", type=float, default=300.0) parser.add_argument("--high", type=float, default=600.0) args = parser.parse_args() x = np.linspace(args.low,args.high,1000) n = get_index(args.density,x,args.temperature) - z = np.polyfit(x-args.x0,1/n, args.deg) - p = np.poly1d(z) + xmin = fmin(chi2,[1,1],args=(x,1/n)) + + print(xmin) - print(p) + a, b = xmin + y = approximate_index(x,a,b) plt.plot(x,1/n,label="data") - plt.plot(x,p(x-args.x0),label="fit") + plt.plot(x,y,label="fit") plt.xlabel("Wavelength (nm)") plt.ylabel("1/n") plt.title("Index of Refraction of Water at 273.15 K") |