% =========================================================================
% MATLAB CODE TO RUN DOUBLE PENDULUM SYSTEM TO
%   - DO NATURE RUN
%   - MAKE OBSERVATIONS FOR LATER ASSIMILATION
%   - PLOT THE ABOVE INFORMATION
% -------------------------------------------------------------------------
% ROSS BANNISTER, NCEO, MAY/JUNE 2010
% =========================================================================

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

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

% 3. Read-in the nature run forecast
% ----------------------------------
[time, ang1, ang2, roc1, roc2] = textread('Assim.Nature', '%f %f %f %f %f\n');

% 4. Read-in the simulated observations and their errors
% ------------------------------------------------------
[obsvalue, obstype, obserr, obstime] = textread('Observations', '%f %u %f %f\n');

% 5. Sort out the observations into type
% --------------------------------------
nobs(1:4) = 0;
obsvalues(1,1:4) = 0.0;
obserrs(1,1:4) = 0.0;
obstimes(1,1:4) = 0.0;
nob = 1;
type = obstype(nob);
while (type ~= 0);
  nobs(type) = nobs(type) + 1;
  thisnumber = nobs(type);
  obsvalues(thisnumber,type) = obsvalue(nob);
  obserrs(thisnumber,type) = obserr(nob);
  obstimes(thisnumber,type) = obstime(nob);
  nob = nob + 1;
  type = obstype(nob);
end;
fprintf('There are %u observation of angle 1\n', nobs(1));
fprintf('There are %u observation of angle 2\n', nobs(2));
fprintf('There are %u observation of rate of change of angle 1\n', nobs(3));
fprintf('There are %u observation of rate of change of angle 2\n', nobs(4));

% 6. Plot the data
% ----------------
subplot(2,2,1);
plot(time, ang1);
xlabel('time');
ylabel('angle');
title('Angle 1 of pendulum');
if (nobs(1)>0);
  hold on;
  errorbar(obstimes(1:nobs(1),1), obsvalues(1:nobs(1),1), obserrs(1:nobs(1),1), 'd');
  hold off;
end;

subplot(2,2,2);
plot(time, ang2);
xlabel('time');
ylabel('angle');
title('Angle 2 of pendulum');
if (nobs(2)>0);
  hold on;
  errorbar(obstimes(1:nobs(2),2), obsvalues(1:nobs(2),2), obserrs(1:nobs(2),2), 'd');
  hold off;
end;

subplot(2,2,3);
plot(time, roc1);
xlabel('time');
ylabel('d(angle)/dt');
title('Rate of change of angle 1 of pendulum');
if (nobs(3)>0);
  hold on;
  errorbar(obstimes(1:nobs(3),3), obsvalues(1:nobs(3),3), obserrs(1:nobs(3),3), 'd');
  hold off;
end;

subplot(2,2,4);
plot(time, roc2);
xlabel('time');
ylabel('d(angle)/dt');
title('Rate of change of angle 2 of pendulum');
if (nobs(4)>0);
  hold on;
  errorbar(obstimes(1:nobs(4),4), obsvalues(1:nobs(4),4), obserrs(1:nobs(4),4), 'd');
  hold off;
end;
