% =========================================================================
% MATLAB CODE TO RUN DOUBLE PENDULUM SYSTEM TO
%   - DO INSERTION RUN
%   - PLOT THE OUTPUT OF THE ASSIMILATION
% -------------------------------------------------------------------------
% ROSS BANNISTER, NCEO, MAY/JUNE 2010
% =========================================================================

% 1. Copy the configuration file to that which is read-in by assim code
% ---------------------------------------------------------------------
!cp Insert.conf Assim.conf

% 2. Run the Fortran code
% -----------------------
!Assim.out

% 3. Read-in the nature run forecast
% ----------------------------------
[time_nat, ang1_nat, ang2_nat, roc1_nat, roc2_nat] = textread('Assim.Nature', '%f %f %f %f %f\n');

% 4. Read-in the observations and their errors
% --------------------------------------------
dd = dir('Assim.obs1');
if (dd.bytes > 0);
  [obstime1,obsvalue1, obserr1] = textread('Assim.obs1', '%f %f %f\n');
  temp = size(obstime1);
  numobs1 = temp(1);
else;
  numobs1 = 0;
end;
dd = dir('Assim.obs2');
if (dd.bytes > 0);
  [obstime2,obsvalue2, obserr2] = textread('Assim.obs2', '%f %f %f\n');
  temp = size(obstime2);
  numobs2 = temp(1);
else;
  numobs2 = 0;
end;
dd = dir('Assim.obs3');
if (dd.bytes > 0);
  [obstime3,obsvalue3, obserr3] = textread('Assim.obs3', '%f %f %f\n');
  temp = size(obstime3);
  numobs3 = temp(1);
else;
  numobs3 = 0;
end;
dd = dir('Assim.obs4');
if (dd.bytes > 0);
  [obstime4,obsvalue4, obserr4] = textread('Assim.obs4', '%f %f %f\n');
  temp = size(obstime4);
  numobs4 = temp(1);
else;
  numobs4 = 0;
end;

% 5. Read-in the insertion run analysis
% -------------------------------------
[time_ins, ang1_ins, ang2_ins, roc1_ins, roc2_ins] = textread('Assim.ins', '%f %f %f %f %f\n');

% 6. Plot the data
% ----------------
subplot(2,2,1);
plot(time_nat, ang1_nat, time_ins, ang1_ins);
legend('nature run', 'insertion analysis');
xlabel('time');
ylabel('angle');
title('Angle 1 of pendulum');
if (numobs1>0);
  hold on;
  errorbar(obstime1, obsvalue1, obserr1, 'd');
  hold off;
end;

subplot(2,2,2);
plot(time_nat, ang2_nat, time_ins, ang2_ins);
legend('nature run', 'insertion analysis');
xlabel('time');
ylabel('angle');
title('Angle 2 of pendulum');
if (numobs2>0);
  hold on;
  errorbar(obstime2, obsvalue2, obserr2, 'd');
  hold off;
end;

subplot(2,2,3);
plot(time_nat, roc1_nat, time_ins, roc1_ins);
legend('nature run', 'insertion analysis');
xlabel('time');
ylabel('d(angle)/dt');
title('Rate of change of angle 1 of pendulum');
if (numobs3>0);
  hold on;
  errorbar(obstime3, obsvalue3, obserr3, 'd');
  hold off;
end;


subplot(2,2,4);
plot(time_nat, roc2_nat, time_ins, roc2_ins);
legend('nature run', 'insertion analysis');
xlabel('time');
ylabel('d(angle)/dt');
title('Rate of change of angle 2 of pendulum');
if (numobs3>0);
  hold on;
  errorbar(obstime4, obsvalue4, obserr4, 'd');
  hold off;
end;
