Optical Trap Code

%Load data

xyscan10=load('xyscan-a-p10-1um-silica.txt');

xyscan20=load('xyscan-a-p20-1um-silica.txt');

xyscan30=load('xyscan-a-p30-xdirection-1um-silica.txt');

xyscan10b=load('xyscan-b-p10-1um-silica.txt');

xyscan20b=load('xyscan-b-p20-xdirection-1um-silica.txt ');

stokes10=load('oscillations-a-p10-1um-silica.txt');

stokes20=load('oscillations-a-p20-1um-silica.txt');

stokes30=load('oscillations-a-p30-1um-silica.txt');

stokes20b=load('oscillations-b-p20-1um-silica.txt');

stokes30b=load('oscillations-b-p30-1um-silica.txt');

stable10=load('p10-1um-silica.txt');

stable20=load('p20-1um-silica.txt');

stable30=load('p30-1um-silica.txt');

% 1= x motion of stage

% 2= y motion of stage

% 3= y motion of bead

% 4= x motion of bead

%Callibration of QPD in order to find V/um

%Based on the figures, the center scan was determined and the most suitable %range for a linear fit was determined. %

% figure(1)

% subplot(2,3,1)

% plot(xyscan10)

% title('xyscan10')

% % range10=[1117:1154]; %

% %Find the best linear fit in this region

% fit10=polyfit(range10',xyscan10(range10,4),1);

% % %Define the best fit function

% x=range10;

% y=fit10(1).*x+fit10(2)

% % %Find the distance traveled over the range used for best fit function using % %the change in stage voltage (1) multiplied by the conversion factor 2.2um/V % % distance10=2.2*(xyscan10(range10(end),1)-xyscan10(range10(1),1)); % % %Find the voltage per um % VpD10=(y(end)-y(1))/distance10 % % % subplot(2,3,2) % plot(xyscan20(:,1)) % title('xyscan20') % % a=[1072:1107]; % fit20=polyfit(a',xyscan20(a,4),1) % x=a; % y=fit20(1).*x+fit20(2); % distance20=2.2*(xyscan20(a(end),1)-xyscan20(a(1),1)); % VpD20=(y(end)-y(1))/distance20 % % % subplot(2,3,3) % plot(xyscan30) % title('xyscan30') % a=[779:822]; % fit30=polyfit(a',xyscan30(a,4),1) % x=a; % y=fit30(1).*x+fit30(2); % distance30=2.2*(xyscan30(a(end),1)-xyscan30(a(1),1)); % VpD30=(y(end)-y(1))/distance30 % % % subplot(2,3,4) % plot(xyscan10b) % title('xyscan10b') % a=[848:872]; % fit10b=polyfit(a',xyscan10b(a,4),1); % x=a; % y=fit10b(1).*x+fit10b(2); % distance10b=2.2*(xyscan10b(a(end),1)-xyscan10b(a(1),1)); % VpD10b=(y(end)-y(1))/distance10b % % subplot(2,3,5) % plot(xyscan20b) % title('xyscan20b') % a=[1334:1372]; % fit20b=polyfit(a',xyscan20b(a,4),1) % x=a; % y=fit20b(1).*x+fit20b(2); % distance20b=2.2*(xyscan20b(a(end),1)-xyscan20b(a(1),1)); % VpD20b=(y(end)-y(1))/distance20b % % %Plot a graph of all data points, including repeats. VpD10b seems to be a % %clear outlier and will not be used for the fit. % VpD=-[VpD10,VpD10b,VpD20,VpD20b,VpD30]; % power=[10,10,20,20,30]; % figure(2) % plot(power,VpD,'rx') % VpDfit=polyfit([10,20,20,30],-[VpD10,VpD20,VpD20b,VpD30],1) % hold on % plot([10,20,30],(VpDfit(1).*[10,20,30]+VpDfit(2))) % title('Voltage/distance vs. power when optical trapping silica beads') % ylabel('Voltage/um') % xlabel('Power') % legend('y=0.0071x+0.290')

%Finding trap stiffness

%Method 1: Equipartition function. This is based off the equation %Kb*T=k*<x^2>

% % % Kb=1.38e-23;%m^2 kg / s^2 K % % Assume temperature is about 300K % % Variance is converted from voltage to distance % figure(3) % subplot(2,2,1) % plot(stable10(:,4)) % title('X Variance, power=10') % % VoltageToDistance10=stable10(:,4).*(1/.1004); % variance10=var(VoltageToDistance10); % stiffM1_10=(Kb*300)/variance10 % % subplot(2,2,2) % plot(stable20(:,4)) % title('X Variance, power = 20') % VoltageToDistance20=stable20(:,4).*(1/.1717); % variance20=var(VoltageToDistance20); % stiffM1_20=(Kb*300)/variance20 % % subplot(2,2,3) % plot(stable30(:,4)) % title('X variance, power = 30') % VoltageToDistance30=stable30(:,4).*(1/.2431); % %variance30=var(VoltageToDistance30); % stiffM1_30=(Kb*300)/variance30

% stiffness10 = % 9.1441e-018 kg/s^2 % % stiffness20 = % 2.6742e-017 kg/s^2 % % stiffness30 = % 5.3636e-017 kg/s^2

%Method 2: PSD Roll Off. Based on fitting the equation: % Sxx=((Kb*temp)./((pi^3)*3*viscosity*diameter.*((f0^2)+(f.^2)))).^0.5

Kb=1.38e-23;%m^2 kg / s^2 K diameter=10^-6; %m viscosity=1; %kg / m s temp=300; %K

figure(15) XQPD10=stable10(:,4); [Pxx10,F10]=pwelch((XQPD10-mean(XQPD10)),[],[],[],22050); % subplot(2,2,1) loglog(F10,Pxx10) title('PSD at power 10') rolloffFit10=nlinfit(F10(780:(end-50)),Pxx10(780:(end-50)),@rolloff,1); stiffM2_10=2*pi*rolloffFit10*3*pi*viscosity*diameter

% XQPD20=stable20(:,4); % [Pxx20,F20]=pwelch((XQPD20-mean(XQPD20)),[],[],[],22050); % subplot(2,2,2) % loglog(F20,Pxx20) % rolloffFit20=nlinfit(F10(622:end),Pxx10(622:end),@rolloff,1); % stiffM2_20=2*pi*rolloffFit20*3*pi*viscosity*diameter % % % XQPD30=stable30(:,4); % [Pxx30,F30]=pwelch((XQPD30-mean(XQPD30)),[],[],[],22050); % subplot(2,2,3) % loglog(F30,Pxx30) % rolloffFit30=nlinfit(F10(724:end),Pxx10(724:end),@rolloff,1); % stiffM2_30=2*pi*rolloffFit30*3*pi*viscosity*diameter

% stiffM2_10 = % -2.4218e-007 % % stiffM2_20 = % -1.1570e-006 % % stiffM2_30 = % -1.2631e-006

% % % % figure subplot(2,3,1) plot(stokes10) title('stokes10a') xstage=stokes10(:,2); xstage=decimate(xstage,100); xQPD=stokes10(:,3); xQPD=decimate(xQPD,100); xQPDmod=xQPD(2:end); velocity=diff(xstage); plot(xQPDmod,velocity,'x')