In [1]:

    

%pylab inline
import matplotlib.pyplot as plt
import numpy as np


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import scipy.ndimage as ndi


    

In [3]:

    

img = ndi.imread("Cr_quartz2_absorption.png")
print img.shape


    

    




(297, 450, 3)

In [4]:

    

#Find indices where the image is red
c = (255, 2, 2)

indices = numpy.where(numpy.all(img == c, axis=-1))
rr = zip(indices[1],297-indices[0]) #zip them together as (x,y)
rr.sort() #sort by x value
ss = np.array(zip(*rr)) #unzip into numpy array for plotting

#Find top line pixel location
c2 = (48, 48, 48)
ind2 = numpy.where(numpy.all(img == c2, axis=-1))
#print zip(ind2[0], ind2[1])


    

In [5]:

    

#Plot the line at the locations of the indices, inverting y and starting at 0
plt.plot(ss[0], ss[1],'r-')
plt.plot(ind2[1],297-ind2[0],'ko')


    

    
Out[5]:





[<matplotlib.lines.Line2D at 0x10d53c410>]






    

In [6]:

    

plt.figure(figsize=(20,10))
plt.imshow(img)
plt.plot(ss[0][82],297-ss[1][82],'bo') #Note, to plot on the image, have to invert y values
plt.plot(ss[0][0:304], 297-ss[1][0:304],'b+') 

plt.plot(ss[0][313],297-ss[1][313],'mo')
plt.plot(ss[0][305:327],297-ss[1][305:327],'m+')

plt.plot(ss[0][328],297-ss[1][328],'go')
plt.plot(ss[0][328:405],297-ss[1][328:405],'g+')

plt.plot(ss[0][415],297-ss[1][415],'co')
plt.plot(ss[0][406:423],297-ss[1][406:423],'c+')

plt.plot(ss[0][700],297-ss[1][700],'ko')
plt.plot(ss[0][424:765],297-ss[1][424:765],'k+')


    

    
Out[6]:





[<matplotlib.lines.Line2D at 0x11011b0d0>]






    

In [7]:

    

print "(%d,%d) = (143,1)"%(ss[0][0],ss[1][0])
print "(%d,%d) = (240,88)"%(ss[0][304],ss[1][304])
x1b=float(ss[0][0]); y1b=float(ss[1][0])
x2b=float(ss[0][304]); y2b=float(ss[1][304])
xab=143.; yab=0.
xbb=240.; ybb=88.
sbx=(xbb-xab)/(x2b-x1b)
sby=(ybb-yab)/(y2b-y1b)
#print sb,sby

x = (x1b-x1b)*sbx + xab
y = (y1b-y1b)*sby + yab
print x1b,y1b," = ",x,y
x = (x2b-x1b)*sbx + xab
y = (y2b-y1b)*sby + yab
print x2b,y2b," = ",x,y


    

    




(76,40) = (143,1)
(164,252) = (240,88)
76.0 40.0  =  143.0 0.0
164.0 252.0  =  240.0 88.0

In [8]:

    

print (ss[0][82],ss[1][82])
print (ss[0][82]-x1b)*sbx + xab,(ss[1][82]-y1b)*sby + yab


    

    




(78, 115)
145.204545455 31.1320754717

In [9]:

    

print "(%d,%d) = (240,88)"%(ss[0][305],ss[1][305])
print "(%d,%d) = (300,89)"%(ss[0][327],ss[1][327])
x1m=float(ss[0][305]); y1m=float(ss[1][305])
x2m=float(ss[0][327]); y2m=float(ss[1][327])
xam=241.; yam=88.
xbm=300.; ybm=89.
smx=(xbm-xam)/(x2m-x1m)
smy=(ybm-yam)/(y2m-y1m)

x = (x1m-x1m)*smx + xam
y = (y1m-y1m)*smy + yam
print x1m,y1m," = ",x,y
x = (x2m-x1m)*smx + xam
y = (y2m-y1m)*smy + yam
print x2m,y2m," = ",x,y


    

    




(165,252) = (240,88)
(186,253) = (300,89)
165.0 252.0  =  241.0 88.0
186.0 253.0  =  300.0 89.0

In [10]:

    

print (ss[0][313],ss[1][313])
print (ss[0][313]-x1m)*smx + xam,(ss[1][313]-y1m)*smy + yam


    

    




(173, 252)
263.476190476 88.0

In [11]:

    

print "(%d,%d) = (300,89)"%(ss[0][328],ss[1][328])
print "(%d,%d) = (700,90)"%(ss[0][405],ss[1][405])
x1g=float(ss[0][328]); y1g=float(ss[1][328])
x2g=float(ss[0][405]); y2g=float(ss[1][405])
xag=301.; yag=89.
xbg=700.; ybg=90.
sgx=(xbg-xag)/(x2g-x1g)
sgy=(ybg-yag)/(y2g-y1g)

x = (x1g-x1g)*sgx + xag
y = (y1g-y1g)*sgy + yag
print x1g,y1g," = ",x,y
x = (x2g-x1g)*sgx + xag
y = (y2g-y1g)*sgy + yag
print x2g,y2g," = ",x,y


    

    




(187,253) = (300,89)
(257,256) = (700,90)
187.0 253.0  =  301.0 89.0
257.0 256.0  =  700.0 90.0

In [12]:

    

print (ss[0][328],ss[1][328])
print (ss[0][328]-x1g)*sgx + xag,(ss[1][328]-y1g)*sgy + yag


    

    




(187, 253)
301.0 89.0

In [13]:

    

print "(%d,%d) = (700,90)"%(ss[0][406],ss[1][406])
print "(%d,%d) = (1000,91)"%(ss[0][423],ss[1][406])
x1c=float(ss[0][406]); y1c= float(ss[1][406])
x2c=float(ss[0][423]); y2c= float(ss[1][423])
xac=701.; yac=90.
xbc=1000.; ybc=91.
scx=(xbc-xac)/(x2c-x1c)
scy=(ybc-yac)/(y2c-y1c)

x = (x1c-x1c)*scx + xac
y = (y1c-y1c)*scy + yac
print x1c,y1c," = ",x,y
x = (x2c-x1c)*scx + xac
y = (y2c-y1c)*scy + yac
print x2c,y2c," = ",x,y


    

    




(258,256) = (700,90)
(274,256) = (1000,91)
258.0 256.0  =  701.0 90.0
274.0 257.0  =  1000.0 91.0

In [14]:

    

print (ss[0][415],ss[1][415])
print (ss[0][415]-x1c)*scx + xac,(ss[1][415]-y1c)*scy + yac


    

    




(267, 256)
869.1875 90.0

In [15]:

    

print "(%d,%d) = (1000,91)"%(ss[0][424],ss[1][424])
print "(%d,%d) = (4750,0)"%(ss[0][765],ss[1][765])
x1k=float(ss[0][424]); y1k=float(ss[1][424])
x2k=float(ss[0][765]); y2k=float(ss[1][765])
xak=1001.; yak=91.
xbk=4750.; ybk=0.
skx=(xbk-xak)/(x2k-x1k)
sky=(ybk-yak)/(y2k-y1k)

x = (x1k-x1k)*skx + xak
y = (y1k-y1k)*sky + yak
print x1k,y1k," = ",x,y
x = (x2k-x1k)*skx + xak
y = (y2k-y1k)*sky + yak
print x2k,y2k," = ",x,y


    

    




(275,257) = (1000,91)
(387,39) = (4750,0)
275.0 257.0  =  1001.0 91.0
387.0 39.0  =  4750.0 0.0

In [16]:

    

print (ss[0][524],ss[1][524])
print (ss[0][524]-x1k)*skx + xak,(ss[1][524]-y1k)*sky + yak


    

    




(342, 240)
3243.70535714 83.9036697248

In [110]:

    

newx = np.zeros(780,dtype=float64)
newy = np.zeros(780,dtype=float64)

for i in range(0,305):
    x1b=float(ss[0][0]); y1b=float(ss[1][0])
    x2b=float(ss[0][304]); y2b=float(ss[1][304])
    xab=143.; yab=0.
    xbb=240.; ybb=88.
    sbx=(xbb-xab)/(x2b-x1b)
    sby=(ybb-yab)/(y2b-y1b)
    newx[i] = (ss[0][i]-x1b)*sbx + xab
    newy[i] = (ss[1][i]-y1b)*sby + yab
    
for i in range(305,328):
    x1m=float(ss[0][305]); y1m=float(ss[1][305])
    x2m=float(ss[0][327]); y2m=float(ss[1][327])
    xam=241.; yam=88.
    xbm=300.; ybm=89.
    smx=(xbm-xam)/(x2m-x1m)
    smy=(ybm-yam)/(y2m-y1m)

    newx[i] = (ss[0][i]-x1m)*smx + xam
    newy[i] = (ss[1][i]-y1m)*smy + yam
    
for i in range(328,406):
    x1g=float(ss[0][328]); y1g=float(ss[1][328])
    x2g=float(ss[0][405]); y2g=float(ss[1][405])
    xag=301.; yag=89.
    xbg=700.; ybg=90.
    sgx=(xbg-xag)/(x2g-x1g)
    sgy=(ybg-yag)/(y2g-y1g)

    newx[i] = (ss[0][i]-x1g)*sgx + xag
    newy[i] = (ss[1][i]-y1g)*sgy + yag
    
for i in range(406,424):
    x1c=float(ss[0][406]); y1c= float(ss[1][406])
    x2c=float(ss[0][423]); y2c= float(ss[1][423])
    xac=701.; yac=90.
    xbc=1000.; ybc=91.
    scx=(xbc-xac)/(x2c-x1c)
    scy=(ybc-yac)/(y2c-y1c)

    newx[i] = (ss[0][i]-x1c)*scx + xac
    newy[i] = (ss[1][i]-y1c)*scy + yac
    
for i in range(424,766):
    x1k=float(ss[0][424]); y1k=float(ss[1][424])
    x2k=float(ss[0][765]); y2k=float(ss[1][765])
    xak=1001.; yak=91.
    xbk=4750.; ybk=0.
    skx=(xbk-xak)/(x2k-x1k)
    sky=(ybk-yak)/(y2k-y1k)

    newx[i] = (ss[0][i]-x1k)*skx + xak
    newy[i] = (ss[1][i]-y1k)*sky + yak
    
for i in range(766,780):
    newx[i] = newx[i-1]+142.
    newy[i] = 0.


    

In [112]:

    

plt.plot(newx,newy,'r-');
plt.xlim(10,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[112]:





<matplotlib.text.Text at 0x113326710>






    

In [113]:

    

plt.plot(newx,newy,'r-');
plt.xlim(2500,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[113]:





<matplotlib.text.Text at 0x117473210>






    

In [23]:

    

#Import the stats package
from scipy import stats
from scipy.optimize import curve_fit


    

In [76]:

    

print "Range of indices to fit over:"
print (np.where(newx > 4247.))


    

    




Range of indices to fit over:
(array([743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755,
       756, 757, 758, 759, 760, 761, 762, 763, 764, 765]),)

In [87]:

    

fitfun = lambda x, a, b, c, d: a + b*x + c*x**2 + d*x**3
funx = np.linspace(0,5000,5000)


    

In [88]:

    

p0 = np.array([0.,1.,1.,1.])
popt, pcov = curve_fit(fitfun, newx[561:689], newy[561:689], p0)


    

In [94]:

    

plt.plot(newx,newy,'r-');
plt.plot(funx[3390:3780],fitfun(funx[3390:3780],popt[0],popt[1],popt[2],popt[3]),'b-')
plt.xlim(2500,5000)
plt.ylim(-1,100)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[94]:





<matplotlib.text.Text at 0x1156fab90>






    

In [66]:

    

p0b = np.array([0.,1.,1.,1.])
poptb, pcovb = curve_fit(fitfun, newx[689:745], newy[689:745], p0b)


    

In [114]:

    

plt.plot(newx,newy,'r-');
plt.plot(funx[3390:3780],fitfun(funx[3390:3780],popt[0],popt[1],popt[2],popt[3]),'b-')
plt.plot(funx[3810:4200],fitfun(funx[3810:4200],poptb[0],poptb[1],poptb[2],poptb[3]),'g-')
plt.xlim(3500,5000)
plt.ylim(-1,30)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[114]:





<matplotlib.text.Text at 0x11784d450>






    

In [118]:

    

p0c = np.array([0.,1.,1.,1.])
poptc, pcovc = curve_fit(fitfun, newx[760:], newy[760:], p0c)


    

In [119]:

    

plt.plot(newx,newy,'r-');
plt.plot(funx[4200:],fitfun(funx[4200:],poptc[0],poptc[1],poptc[2],poptc[3]),'b-')
plt.xlim(4000,5000)
plt.ylim(-1,20)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[119]:





<matplotlib.text.Text at 0x117d81790>






    

In [86]:

    

funx = np.linspace(0,5000,5000)
plt.plot(newx,newy,'r-');
plt.plot(funx,fitfun(funx,poptc[0],poptc[1],poptc[2],poptc[3]),'b-')
plt.xlim(4000,5000)
plt.ylim(-1,10)
plt.xlabel("Wavelength (nm)",fontsize=20)
plt.ylabel("Transmission (%)",fontsize=20)


    

    
Out[86]:





<matplotlib.text.Text at 0x11452d4d0>






    

In [85]:

    

np.where(fitfun(funx,poptc[0],poptc[1],poptc[2],poptc[3])<0)


    

    
Out[85]:





(array([4361, 4362, 4363, 4364, 4365, 4366, 4367, 4368, 4369, 4370, 4371,
        4372, 4373, 4374, 4375, 4376, 4377, 4378, 4379, 4380, 4381, 4382,
        4383, 4384, 4385, 4386, 4387, 4388, 4389, 4390, 4391, 4392, 4393,
        4394, 4395, 4396, 4397, 4398, 4399, 4400, 4401, 4402, 4403, 4404,
        4405, 4406, 4407, 4408, 4409, 4410, 4411, 4412, 4413, 4414, 4415,
        4416, 4417, 4418, 4419, 4420, 4421, 4422, 4423, 4424, 4425, 4426,
        4427, 4428, 4429, 4430, 4431, 4432, 4433, 4434, 4435, 4436, 4437,
        4438, 4439, 4440, 4441, 4442, 4443, 4444, 4445, 4446, 4447, 4448,
        4449, 4450, 4451, 4452, 4453, 4454, 4455, 4456, 4457, 4458, 4459,
        4460, 4461, 4462, 4463, 4464, 4465, 4466, 4467, 4468, 4469, 4470,
        4471, 4472, 4473, 4474, 4475, 4476, 4477, 4478, 4479, 4480, 4481,
        4482, 4483, 4484, 4485, 4486, 4487, 4488, 4489, 4490, 4491, 4492,
        4493, 4494, 4495, 4496, 4497, 4498, 4499, 4500, 4501, 4502, 4503,
        4504, 4505, 4506, 4507, 4508, 4509, 4510, 4511, 4512, 4513, 4514,
        4515, 4516, 4517, 4518, 4519, 4520, 4521, 4522, 4523, 4726, 4727,
        4728, 4729, 4730, 4731, 4732, 4733, 4734, 4735, 4736, 4737, 4738,
        4739, 4740, 4741, 4742, 4743, 4744, 4745, 4746, 4747, 4748, 4749,
        4750, 4751, 4752, 4753, 4754, 4755, 4756, 4757, 4758, 4759, 4760,
        4761, 4762, 4763, 4764, 4765, 4766, 4767, 4768, 4769, 4770, 4771,
        4772, 4773, 4774, 4775, 4776, 4777, 4778, 4779, 4780, 4781, 4782,
        4783, 4784, 4785, 4786, 4787, 4788, 4789, 4790, 4791, 4792, 4793,
        4794, 4795, 4796, 4797, 4798, 4799, 4800, 4801, 4802, 4803, 4804,
        4805, 4806, 4807, 4808, 4809, 4810, 4811, 4812, 4813, 4814, 4815,
        4816, 4817, 4818, 4819, 4820, 4821, 4822, 4823, 4824, 4825, 4826,
        4827, 4828, 4829, 4830, 4831, 4832, 4833, 4834, 4835, 4836, 4837,
        4838, 4839, 4840, 4841, 4842, 4843, 4844, 4845, 4846, 4847, 4848,
        4849, 4850, 4851, 4852, 4853, 4854, 4855, 4856, 4857, 4858, 4859,
        4860, 4861, 4862, 4863, 4864, 4865, 4866, 4867, 4868, 4869, 4870,
        4871, 4872, 4873, 4874, 4875, 4876, 4877, 4878, 4879, 4880, 4881,
        4882, 4883, 4884, 4885, 4886, 4887, 4888, 4889, 4890, 4891, 4892,
        4893, 4894, 4895, 4896, 4897, 4898, 4899, 4900, 4901, 4902, 4903,
        4904, 4905, 4906, 4907, 4908, 4909, 4910, 4911, 4912, 4913, 4914,
        4915, 4916, 4917, 4918, 4919, 4920, 4921, 4922, 4923, 4924, 4925,
        4926, 4927, 4928, 4929, 4930, 4931, 4932, 4933, 4934, 4935, 4936,
        4937, 4938, 4939, 4940, 4941, 4942, 4943, 4944, 4945, 4946, 4947,
        4948, 4949, 4950, 4951, 4952, 4953, 4954, 4955, 4956, 4957, 4958,
        4959, 4960, 4961, 4962, 4963, 4964, 4965, 4966, 4967, 4968, 4969,
        4970, 4971, 4972, 4973, 4974, 4975, 4976, 4977, 4978, 4979, 4980,
        4981, 4982, 4983, 4984, 4985, 4986, 4987, 4988, 4989, 4990, 4991,
        4992, 4993, 4994, 4995, 4996, 4997, 4998, 4999]),)

In [280]:

    

x,y = np.loadtxt("DefaultDataset.csv",delimiter=',',dtype=float64,unpack=True)


    

In [281]:

    

plt.plot(x,y,'k-')


    

    
Out[281]:





[<matplotlib.lines.Line2D at 0x12d48ec10>]






    

In [ ]: