Code covered by the BSD License  

Highlights from
Cavendish Balance Experiment

Cavendish Balance Experiment

by

 

This function to used for the undergraduation Cavendish gravitational experiment.

cavendish_ver_2
function cavendish_ver_2

% This script reads from a three column Excel file of type 'xlsx'
% and reduces the data taken from a Cavendish Balance 
% experiment to compute G and other interesting information.

clc   % clear the screen display
%=============================================================
%=============================================================
% distance = the distance from the mirror to the screen in meters

fid=fopen('cavendish_out.txt','w+');

%=============================================================
%=============================================================
%cave_file='cavendish1';    % Excel data file name (omit the xlxs extention)
%===================================================================
%cave_file='Data_91712_2';
%cave_file='data_out';
%cave_file='sample_data';
%cave_file='test_sample_data';
cave_file='cjs_cavendish_experiment';
%===================================================================
cave_file=strcat(cave_file,'.xlsx');

disp(strcat('Reading Data File--',cave_file))

y=xlsread(cave_file);% The entire path name is not required if the 
                     % folder containing the file is in the search path
                     
%distance=input('Enter the distance between the Cavendish dumbbell mirror and the screen (meters) ');
distance=y(1,5);

fprintf(1,'The distance from the Cavendish mirror to the screen = %3.2f meters\n',distance);                   
fprintf(fid,'The distance from the Cavendish mirror to the screen = %3.2f meters\r\n',distance);
disp('===================================================');

x1=[(y(:,1))*60+y(:,2),y(:,3)];       % create a single time
                                      % column in seconds 
x2=[x1(:,1)-x1(1,1),x1(:,2)-x1(1,2)]; % set the starting values of time
                                      % and position to zero

time=x2(:,1);       % time column vector - seconds
position=x2(:,2);   % position column vector - centimeters

if mean(position) < 0   % if the positions data points are decending,
    position=-position; % reverse the order
end % if mean(position)<0

plot(time,position,'o'); grid on;
axis([0 time(end) 0 1.1*max(position)]);
title('Cavendish Balance Screen Spot Location vs. Time','fontsize',17);
xlabel('Time  (sec)','fontsize',13);   ylabel('Spot Position  (cm)','fontsize',13);
%===============================================================
peaks=find_peaks(position); % custom function "find_peaks"
[~,n_peaks]=size(peaks);
l_time=length(time); % determine the dimensions of time


fprintf(1,'The number of data points = %3i\n',l_time);
fprintf(1,'The average time between data points = %4.2f seconds\n',time(l_time)./l_time);
fprintf(1,'The total time of taking data = %4.2f minutes\n',time(end)./60);

fprintf(fid,'The number of data points = %3i\r\n',l_time);
fprintf(fid,'The average time between data points = %4.2f seconds\r\n',time(l_time)./l_time);
fprintf(fid,'The total time of taking data = %4.2f minutes\r\n',time(end)./60);
disp('===================================================');

disp(' Measured Extremum Data Points');
fprintf(1,'%5d',peaks(1,:)); fprintf(1,' -- Extremum Point Indices\n');
fprintf(1,' %8.2f',time(peaks(1,:))); fprintf(1,' -- Ext. Point Times (sec)\n');
fprintf(1,'%8.2f',peaks(2,:)); fprintf(1,' -- Ext. Point Magnitude (cm)\n');

fprintf(fid,'%5d',peaks(1,:)); fprintf(fid,' -- Extremum Point Indices\r\n');
fprintf(fid,' %8.2f',time(peaks(1,:))); fprintf(fid,' -- Ext. Point Times (sec)\r\n');
fprintf(fid,'%8.2f',peaks(2,:)); fprintf(fid,' -- Ext. Point Magnitude (cm)\r\n');
disp('===================================================');

peaks_interp=find_peaks_interp(time,position,peaks); 
% custom function "find_peaks_interp"

disp(' Interpolated Extremum Points')
fprintf(1,' %8.2f',peaks_interp(1,:)); fprintf(1,' -- Interp. Ext. Times (sec)\n');
fprintf(1,'%8.2f',peaks_interp(2,:)); fprintf(1,' -- Interp. Ext. Magnitude (cm)\n');

fprintf(fid,' %8.2f',peaks_interp(1,:)); fprintf(fid,' -- Interp. Ext. Times (sec)\r\n');
fprintf(fid,'%8.2f',peaks_interp(2,:)); fprintf(fid,' -- Interp. Ext. Magnitude (cm)\r\n');
disp('====================================================');


period=(2.*(peaks_interp(1,n_peaks)-peaks_interp(1,1)))./(n_peaks-1); 
% measured period in seconds

dampcoeff=damp_coeff(peaks_interp); % custom function computes 
                                    % the damping coefficient of
                                    % the data

fprintf(1,'System Damping Coefficient = %5.3f\n',dampcoeff);
fprintf(fid,'System Damping Coefficient = %5.3f\r\n',dampcoeff);


[asym, asym_std]=find_asymptote(dampcoeff,peaks_interp);% custom function
% to compute the asymptotic value in centimeters

fprintf(1,'The estimated light spot final position is %5.3f +- %5.3f cm\n',asym,asym_std);
fprintf(fid,'The estimated light spot final position is %5.3f +- %5.3f cm\r\n',asym,asym_std);


period_0=(sqrt(1-dampcoeff.^2)).*period;
% the undamped period of oscillation

fprintf(1,'The measured period = %5.3f minutes\n',period/60);
fprintf(1,'The undamped period = %5.3f minutes\n',period_0/60);

fprintf(fid,'The measured period = %5.3f minutes\r\n',period/60);
fprintf(fid,'The undamped period = %5.3f minutes\r\n',period_0/60);



const=7.474e-4; % a constant depending on equipment dimensions and masses

gg=(const.*(asym./100))./((period.^2)*distance); 
gg=gg./(1-dampcoeff^2); %correction from the measured period
                        % i.e. ((1-dampcoeff^2)^1/2)*w0 arguements
                        % of sin(arg) and cos(arg) of response 
                        % waveform

% computed value of gravitational constant
                                                
fprintf(1,'The measured gravitational constant -- G = %5.3e Nm^2/kg^2\n',gg);
fprintf(fid,'The measured gravitational constant -- G = %5.3e Nm^2/kg^2\r\n',gg);

gg_acp=6.67*10^-11;
pe=((gg-gg_acp)/gg_acp)*100;
fprintf(1,'The Percent Error is %4.2f%%\n',pe);
fprintf(fid,'The Percent Error is %4.2f%%\r\n',pe);

fclose(fid); % close cavendish_out.txt
disp('The file ''cavendish_out.txt'' has been written to the current folder.');
disp('====================================================');
%================================================================
%================================================================
% The following does a least square fit.
disp('The following is a repeat of the calculation by');
disp('using the lsqcurvefit funcion to match the data');
disp('to a damped sinusoid objective function');

phase=atan(sqrt(1-dampcoeff^2)/dampcoeff)+pi;
x0=[asym period phase dampcoeff asym]; % original values take from oudtput
                                       % from above
                                       
[x_lsq,resnorm]=lsqcurvefit(@cavendish_obj,x0,time,position);                                       

disp('=======================================================');
disp('The results of a least square fit of a damped sinusiod');
disp('to the original data are shown below');
disp('=======================================================');
fprintf(1,'The RMSE is %3.2f\n',resnorm);
fprintf(1,'System Damping Coefficient = %5.3f\n',x_lsq(4));
fprintf(1,'The estimated light spot final position is %5.3f cm\n',x_lsq(5));
fprintf(1,'The measured period = %5.3f minutes\n',x_lsq(2)/60);
tt=sqrt(1-x_lsq(4)^2)*x_lsq(2);     % The undamped resonant period
fprintf(1,'The undamped period = %5.3f minutes\n',tt/60);

gg_lsq=const*(x_lsq(5)/((x_lsq(2)^2)*100*distance));

fprintf(1,'The measured value of G from the lsq fit is %5.3e Nm^2/kg^2\n',gg_lsq);
pe_lsq=((gg_lsq-gg_acp)/gg_acp)*100;
fprintf(1,'The percentage error for the lsq fit is %3.2f%%\n',pe_lsq);
disp('========================================================');
end
%=================================================================
%=================================================================
%=================================================================
%FUNCTIONS used in cavendish.m follow
% find_peaks
% find_peaks_interp
% damp_coeff
% find_asymptotes
%=================================================================
%=================================================================
%=================================================================
function [ extremum ] = find_peaks( yy )
% The following code finds the "peaks" or extremum of the 
% damped sinusoidal data.  These are needed to compute the period
% and the damping coeficient of the measured data pattern.   The 
% results are stored in the array "extremem" below.

n_length=length(yy); % length of the vector or sequence 'yy'

polarity=yy(n_length)-yy(1);% determine the direction of the data

peak_mag=zeros(1,10); peak_index=zeros(1,10);% initialize memory
n=1; ii=1;       %counters for the folowing while loop
if polarity  > 0 % the recorded data are increasing in value
    while ii<n_length      
        [mag,index]=max(yy(ii:n_length));
        peak_mag(n)=mag; 
        
        if n==1
            peak_index(1)=index;
        else % if n>1
            peak_index(n)=peak_index(n-1)+index;
        end % if n=1 
        
        ii=peak_index(n)+1; n=n+1;
        if ii<n_length
            [mag,index]=min(yy(ii:n_length));
            peak_mag(n)=mag;
            peak_index(n)=peak_index(n-1)+index;
            ii=peak_index(n)+1; n=n+1; 
        end % if ii<n_length
    end % while ii<n_length
else % if polarity <0 % the recorded data are decreasing in value
    while ii<n_length
       [mag,index]=min(yy(ii:n_length));
        peak_mag(n)=mag;
        
        if n==1
            peak_index(1)=index;
        else % if n>1
            peak_index(n)=peak_index(n-1)+index;
        end % if n=1
        
        ii=peak_index(n)+1; n=n+1;
        if ii<n_length
            [mag,index]=max(yy(ii:n_length));
            peak_mag(n)=mag;
            peak_index(n)=peak_index(n-1)+index;
            ii=peak_index(n)+1; n=n+1;
        end % if ii<n_length
    end % while ii<n_length
end % if polarity 

mm=find(peak_index,1,'last')-1;
peakindex=peak_index(1:mm);
peakmag=peak_mag(1:mm);
% the magnitude and position index results are stored in the
% array "extremum."
extremum=[peakindex;peakmag];
% ================================================================
end % function extremum
% ================================================================
%=================================================================
function [ peaks ] = find_peaks_interp(x,y,extremum)
% This function uses 2*pp+1 points arount the maximum or minimum 
% point in the sequence to do quadratic interperlation to find the 
% max or min location or magnitude of the underlying 
% sampled function.S

[m,n]=size(extremum);
peaks=zeros(m,n);
pp=1;   % number of points away from the maximum to use to 
        % form interperlating polynomial, 2*pp+1 points used
for i=1:n
    mm=extremum(1,i);
    x_array=x(mm-pp:mm+pp);
    y_array=y(mm-pp:mm+pp);
    x_array=x_array./100; %scaling to improve numerical precision
    p=polyfit(x_array,y_array,2);
    peaks(1,i)=(-p(2))./(2.*p(1)); % time location
    peaks(1,i)=100*peaks(1,i); % to reverse scaling above
    peaks(2,i)=(-p(2).*p(2))./(4.*p(1))+p(3); % position
end % for i=1:n

end % function find_peaks_interp
% =======================================================
%=================================================================
function [ dampcoeff ] = damp_coeff( peaks_interp )
% This function computes the damping coefficient of the
% damped sinusoidal measured data based on the
% first two extremum points.

n=peaks_interp(2,1)./peaks_interp(2,2);

x=(-1+sqrt(1+4.*n.*(n-1)))./(2.*n);

damp_sqr=log(x).*log(x);
damp_sqr=damp_sqr./(pi.*pi+log(x).*log(x));

dampcoeff=sqrt(damp_sqr);

end % function damp_coeff
% =============================================================
function [ asymptote_ave, asymptote_std ]=find_asymptote( dampcoeff,peaks_interp )
% This functions computes the assymtotic value of the damped sinusoid 
% data.

[~,n]=size(peaks_interp);   % number of columns in peaks_interp
aa=1:n;
dc=dampcoeff;   % shorten the variable namae to simplify the code
dcs=sqrt(1-dc.*dc);

peak=1+(-1).^(aa+1).*(exp(-aa.*pi.*dc./dcs)); % this computes the normalized
% valued of each maximum and minumum of the damped sinusoid

asymptote=peaks_interp(2,:)./peak;

asymptote_ave=sum(asymptote)./length(asymptote);

asymptote_std=std(asymptote);

end % function find_asymptotes
%=================================================================
function f_ydata = cavendish_obj( x,xdata )
%This is the objective function for curve
%fitting the Cavendish Balance data to a
%damped sinusoid

f_ydata=x(1).*(sin((2.*pi/x(2)).*xdata+x(3))).*exp((-x(4).*2*pi/x(2)).*xdata)+x(5);

end % function f = canvendish



Contact us