 function fit_trajectories
%fit_trajectories determines the parabolic trajectory of the electron beam in
%the presence of an E-field.
%
%The input is a picture of the electron beam and a calibration of the
%camera
%
%The user inputs include
%1) name of picture: file_name
%2) calibration of the camera in terms of: pixels_per_cm
%3) threshold value for discriminating the beam vs. the background.
%
%Output is
%1) appending to the parabolic coefficients in the file results.csv (open
%this with xls).
%2) saved figure: filename.FIG
%
%Input paragraphs are delimited using %************************
%
%sak@wpi.edu 1/10/12


%inputs are here
%****************************************
file_name = 'Bdown.jpg';
%to calibrate the picture, select two points xy1 and xy2. For each of these
%two points eneter their coordinates in pixels and in cm.
xy1_pixel = [527 124]; 
xy1_cm = [2 2];
xy2_pixel = [123 343];
xy2_cm = [10 -2];
%select threshold value of intensity to darken pixels which are not
%showing the electron path
threshold = 200;
%select ROI in cm which contains the electron path
ROI_cm = [2 0; 10 -2.5]; %[x1 y1; x2 y2]
%****************************************

pixel2cm = @(xy_pixel) repmat( xy1_cm, size( xy_pixel,1), 1) + ...
    (xy_pixel - repmat( xy1_pixel, size( xy_pixel, 1), 1)).*repmat( (xy2_cm - xy1_cm)./(xy2_pixel-xy1_pixel), size( xy_pixel, 1), 1) ;
cm2pixel = @(xy_cm) repmat( xy1_pixel, size( xy_cm,1), 1) + ...
    (xy_cm - repmat( xy1_cm, size( xy_cm, 1), 1)).*repmat( (xy2_pixel - xy1_pixel)./(xy2_cm-xy1_cm), size( xy_cm, 1), 1);

im = imread( file_name );
[x, y] = meshgrid( 1:size(im,2), 1:size(im,1) );
im = single( im(:,:,3) );
figure(sum(file_name)); set( clf, 'name', file_name );
imagesc( im ); axis image; colormap( 'gray' );
im(im<threshold) = 0;
imagesc( im );
colormap( 'gray' );
axis image
hold on;
ROI_pixel = cm2pixel( ROI_cm );
plot( ROI_pixel([1 1 2 2 1],1), ROI_pixel([1 2 2 1 1 ],2), '-r' );

mask = (x>min(ROI_pixel(:,1)))&(x<max(ROI_pixel(:,1)))&(y>min(ROI_pixel(:,2)))&(y<max(ROI_pixel(:,2)));
traj = [ sum(x.*im.*mask)./sum(im.*mask); sum(y.*im.*mask)./sum(im.*mask) ]';
traj( any( isnan(traj), 2 ), :) = [];
plot( traj(:,1), traj(:,2), '-r' );
traj_cm = pixel2cm( traj );
% conventional way
% coeff_cm = polyfit( traj_cm(:,1), traj_cm(:,2), 2 );
%better, because it forces the slope to zero.
coeff_cm = [traj_cm(:,1).^2 traj_cm(:,1).^0] \ traj_cm(:,2);
coeff_cm = [coeff_cm(1) 0 coeff_cm(2)];
plot( traj(:,1), cm2pixel( polyval( coeff_cm, traj_cm(:,[1 1]) ) )*[0;1], '-g' );
title( sprintf( 'y = %g x^2/cm + %g x + %g cm', coeff_cm ) );
saveas( gcf, strcat( file_name, '.fig' ) );
%%
fprintf( 'saving (appending ) to results.csv\n' );
fid = fopen( 'results.csv', 'a' );
fprintf( fid, '\n%s, %g, cm^-1, %g, cm^0, %g, cm^1\n', file_name, coeff_cm )
fclose( fid );
