
## Usage:
##   [force, time] = static_analysis(data, cal_start, cal_end, cal_grams, zero_start, zero_end, thrust_start, thrust_end, frequency)
##
## Performs analysis on static motor thrust curve data.  This function
## does automatic calibration, burn start, burnout and ejection charge
## detection.
##
## The data should contain a portion where a known calibration weight is
## weighed and a portion where no excess weight is present.
##
## Arguments:
##   data - the raw data as an array
##   cal_start - start index when the known calibration weight is on
##   cal_end   - end index when the known calibration weight is on
##   cal_grams - weight of the calibration weight in grams
##   zero_start - start index where no excess weight is on
##   zero_end   - end index where no excess weight is on
##   thrust_start - approximate start index for thrust curve measurement
##   thrust_end   - approximate end index for thrust curve measurement
##   frequency - sampling frequency of the data
##
## Typical usage:
##   data = load("SF-A8.csv")(:, 4)  # Load raw data into a vector
##   plot(data)   # Plot data to detect approximate indices
##   [f,t] = static_analysis(data, 18000, 25000, 886.2, 45000, 51000, 51000, 57000, 1000);
##

function [force, time] = static_analysis(data, cal_start, cal_end, cal_grams, zero_start, zero_end, thrust_start, thrust_end, frequency)

  % Gravity acceleration is 9.82 m/s^2 in Finland
  G = 9.82;

  # Find calibration values
  cal_average = mean(data(cal_start:cal_end));
  zero_average = mean(data(zero_start:zero_end));
  zero_max = max(data(zero_start:zero_end));

  thrust_data = data(thrust_start:thrust_end);
  t0_limit = zero_max + (zero_max - zero_average)/2;

  # Find thrust curve start and end positions
  indices = find(data(thrust_start:thrust_end) > t0_limit) + thrust_start;
  t0 = indices(1);
  t1 = indices(length(indices));

  start_index = t0 - frequency/2;
  end_index = t1 + frequency/2;

  # Extract force and time data
  force = data(start_index:end_index);
  time = (1:length(force)) / frequency - 0.5;

  # Check if an ejection charge was present
  cropped = force(1:(length(force) - frequency));
  indices = find(cropped > t0_limit);
  burnout_index = indices(length(indices));

  if (burnout_index < length(cropped)-frequency/4)
    # An ejection charge was present
    burnout_time = burnout_index / frequency - 0.5;

    cropped = force((length(force) - frequency):length(force));
    indices = find(cropped > t0_limit);
    ejection_index = indices(1) + (length(force) - frequency);
    ejection_time = ejection_index / frequency - 0.5;
    ejection_charge = 1;

  else

    burnout_index = length(force) - frequency/2;
    burnout_time = max(time) - 0.5;
    ejection_time = burnout_time;
    ejection_charge = 0;

  endif


  # Scale force
  force = force - zero_average;
  force = force / (cal_average - zero_average) * cal_grams / 1000 * G;

  # Calculate specifics
  total_impulse = sum(force((frequency/2):burnout_index)) / frequency;
  average_thrust = total_impulse / burnout_time;

  # Plot
  ymin = floor(min(force));
  ymax = ceil(max(force));
  vert = [-1000000, 1000000];
  plot(time, force, 'b', [0 0], vert, 'r', [burnout_time burnout_time], vert, 'r', [ejection_time ejection_time], vert, 'r', [min(time) max(time)], [0 0], 'r', [0 burnout_time], [average_thrust average_thrust], 'r');
  axis([min(time) max(time) ymin ymax]);
  xlabel("Time / s");
  ylabel("Thrust / N");

  
  printf("Burnout time:         %.2f s\n", burnout_time);
  if (ejection_charge)
    printf("Ejection charge time: %.2f s\n", ejection_time);
    printf("Delay time:           %.2f s\n", ejection_time - burnout_time);
  endif
  printf("Total impulse:        %.2f Ns\n", total_impulse);
  printf("Average thrust:       %.2f N\n", total_impulse/burnout_time);

endfunction