% Run automatic calibration on all fibers.

% Clear workspace.
clear;

% Process all fibers.
% Each fiber must have an associated '.mat' and '.mapping' file in data/adc-50kHz directory.
fibers = {'fiber10'}; % list of all fibers to be processed

for fiber = 1:length(fibers)
  fibername = fibers{fiber};

  % Generate expected filename of mapping file that describes which traces go with which original force data.
  mapping_filename = sprintf('data/adc-50kHz/%s.mapping', fibername);

  % Read corresponding forcedata directory name.
  forcedata_directory = textread(mapping_filename, '%s', 1);
  
  % Read mapping.
  [forcenames, adcfilenames, traceindices] = textread(mapping_filename, '%s %s %d', 'headerlines', 1);

  % Iterate over each listed trace.
  N = length(forcenames);
  a_n = zeros(N,1);
  b_n = zeros(N,1);

  % Storage for calibrated force traces.
  trace = {};
  ntraces = 0;

  % DEBUG
  avgforce_n = zeros(N,1);
  avgadc_n = zeros(N,1);

  for n = 1:N
    forcetrace_filename = sprintf('data/forcedata/%s/%sA.txt', char(forcedata_directory), char(forcenames{n}));
    adc_filename = sprintf('data/adc-50kHz/%s.mat', char(adcfilenames{n}));
    trace_index = traceindices(n);

    show_plots = true; % generate plots
%     try
      % Attempt automatic force calibration.
      [a, b, F, avgforce, avgadc] = force_calibration_automatic(forcetrace_filename, adc_filename, trace_index, show_plots);

      % Store force conversion parameters.
      a_n(n) = a;
      b_n(n) = b;
      
      % Store calibrated force trace.
      ntraces = ntraces + 1;
      trace{ntraces} = F;

      % DEBUG
      avgforce_n(n) = avgforce;
      avgadc_n(n) = avgadc;

      % Save plot.
      filename = sprintf('plots/%s-%d.eps', fibername, n);
      print('-depsc', filename);
      command = sprintf('epstopdf %s', filename);
      system(command);
%     catch exception
%       % Automatic force calibration failed for this trace.
%       disp(exception);
%       fprintf('WARNING: Processing of trace %d unsuccessful.\n', n);
%     end

    fprintf('\n');
  end

  % Compute statistics over all traces from this fiber.
  a = mean(a_n);
  da = std(a_n) / sqrt(N);
  b = mean(b_n);
  db = std(b_n) / sqrt(N);  

  fprintf('STASTICS OVER ALL TRACES FOR FIBER %d\n', fiber);
  fprintf('a = %f +- %f\n', a, da);
  fprintf('b = %f +- %f\n', b, db);

  % Make plot of statistics over all traces from this fiber.
  clf;
  subplot(2,1,1);
  plot(a_n, 'k.');
  ylabel('a / pN');
  title(sprintf('fiber %d', fiber));
  subplot(2,1,2);
  plot(b_n, 'k.');
  ylabel('b / pN');
  xlabel('trace');

  % Save plot.
  filename = sprintf('plots/%s-tracefits.eps', fibername);
  print('-depsc', filename);
  command = sprintf('epstopdf %s', filename);
  system(command);
  
  % Write data.
  filename = sprintf('calibrated/%s.mat', fibername);
  save(filename, 'trace');
end


