function analyze_Cavendish
%analyze_Cavendish requires the angular deflection of the pendulum obtained
%from digitize_Cavendish_movie.M
%
%The input is the .MAT file from digitizing of the torsion pendulum's oscillations for the lead
%masses in one position and switching them to the other position with the
%accompanying ring-down.
%
%The physically important quantities are determined from the
%ring-down of the oscillations about their equillibrium points. The desired
%measurements are the neutral point of the oscillations, the damping 
%constant and the angular frequency. For this the user is required to
%identify the duration of the two ring-down regimes, and provide some good
%initial guesses for the curve-fitting.
%
%
%Input paragraphs are delimited using %************************
% It saves the fitting parameters appended into results.csv
% It also saves the figure as a .pdf and .fig
% sak@wpi.edu 1/10/12

%**********************
file_name = 'Cavendish2.mov';

%format is [starting time for ring-down, ending time for ring-down] (in
%sec)
T1 = [300 1700];
T2 = [2400 3500]; 

%fitting guesses for the two ring-down oscillations, which have format 
%[ampl, b, theta_mean, omega, phi, t]. You will have to play with these
param1 = [  0.0103 %ampl
    0.1250 %b
    0.01  %theta_mean
    0.0134 %omega
    2.5752];  %phi
param2 = [  0.08 %ampl
    0.09 %b
    -0.01  %theta_mean
    0.013 %omega
    1];  %phi
%***********************

load( strcat( file_name, '.mat' ) );
fnc = @(ampl, b, theta_mean, omega_0, phi, t) ampl*cos(sqrt(1-b^2/4)*omega_0*t+phi).*exp(-b*omega_0*t/2)+theta_mean;

figure( sum( file_name ) ); set( clf, 'name', file_name );
subplot(2, 1, 1);
plot( t, theta(:,1), '.' );
hold on;
plot( repmat( T1, 2, 1), repmat( ylim, 2, 1)' );
plot( t(t>=T1(1)&t<=T1(2)), fnc(param1(1), param1(2), param1(3), param1(4), param1(5), t(t>=T1(1)&t<=T1(2))), '-g' );
[param1, val] = LMFnlsq( @(param) theta(t>=T1(1)&t<=T1(2)) - fnc(param(1), param(2), param(3), param(4), param(5), t(t>=T1(1)&t<=T1(2)))', param1 );
plot( t(t>=T1(1)&t<=T1(2)), fnc(param1(1), param1(2), param1(3), param1(4), param1(5), t(t>=T1(1)&t<=T1(2))), '-k' );
title( sprintf( 'first region: b = %.6g, \\theta_0 = %.6g, \\omega_0 = %.6g', param1(2:4) ) );  

plot( repmat( T2, 2, 1), repmat( ylim, 2, 1)' );
plot( t(t>=T2(1)&t<=T2(2)), fnc(param2(1), param2(2), param2(3), param2(4), param2(5), t(t>=T2(1)&t<=T2(2))), '-g' );
[param2, val] = LMFnlsq( @(param) theta(t>=T2(1)&t<=T2(2)) - fnc(param(1), param(2), param(3), param(4), param(5), t(t>=T2(1)&t<=T2(2)))', param2 );
plot( t(t>=T2(1)&t<=T2(2)), fnc(param2(1), param2(2), param2(3), param2(4), param2(5), t(t>=T2(1)&t<=T2(2))), '-k' );
axis tight;

xlabel( 't (sec)' );    
ylabel( '\theta (rad)' );
title( {sprintf( 'first region: b = %.6g, \\theta_0 = %.6g, \\omega_0 = %.6g', param1(2:4) );  
    sprintf( 'second region: b = %.6g, \\theta_0 = %.6g, \\omega_0 = %.6g', param2(2:4) ) } );

subplot(2, 1, 2);
plot( t(T1(2)<t&T2(1)>t), theta(T1(2)<t&T2(1)>t) );
axis tight;
xlabel( 't (sec)' );
ylabel( '\theta (rad)' );

fprintf( 'saving (appending ) to results.csv\n' );
    fid = fopen( 'results.csv', 'a' );
    fprintf( fid, '\n%s, %g, %g, %g, %g, %g, %g, %g\n', file_name, T1, param2 );
    fclose( fid );

saveas( gcf, strcat( file_name, '.fig' ) );
saveas( gcf, strcat( file_name, '.pdf' ) );
return

%here is a much better function than fminsearch for fitting the parameters
%do not touch!!!!!!!!!!
