% =========================================================================
% MATLAB CODE TO RUN DOUBLE PENDULUM SYSTEM TO
%   - DO VARIATIONAL ASSIMILATION 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 Var.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 var analysis
% ---------------------------
[time_var, ang1_var, ang2_var, roc1_var, roc2_var] = textread('Assim.anal', '%f %f %f %f %f\n');

% 6. Read-in the forecast from the background
% -------------------------------------------
[time_bck, ang1_bck, ang2_bck, roc1_bck, roc2_bck] = textread('Assim.fore', '%f %f %f %f %f\n');

% 7. Read-in the cost function values
% -----------------------------------
%[Cycle, Loop, Jo, Jb, J] = textread('Assim.J', '%u %u %f %f %f\n');

% 8. Plot the assimilation data
% -----------------------------
figure(1);

subplot(2,2,1);
plot(time_nat, ang1_nat,  time_var, ang1_var);
legend('nature run', 'variational 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_var, ang2_var);
legend('nature run', 'variational 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_var, roc1_var);
legend('nature run', 'variational 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_var, roc2_var);
legend('nature run', 'variational 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;


% 9. Plot the background forecast data
% ------------------------------------
figure(2);

subplot(2,2,1);
plot(time_nat, ang1_nat,  time_bck, ang1_bck);
legend('nature run', 'background forecast');
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_bck, ang2_bck);
legend('nature run', 'background forecast');
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_bck, roc1_bck);
legend('nature run', 'background forecast');
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_bck, roc2_bck);
legend('nature run', 'background forecast');
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;

% 10. Plot the cost function data
% -------------------------------
%figure(3);

%plot (Jb);
%hold on;
%plot (Jo);
%plot (J);