N=100;       % number of time to run the model
K=0.04;      % diffusion coefficient
I=30; J=30;  % domain size

% Velocity field is given below (DO NOT CHANGE)
u=zeros(I,J,1);
v=zeros(I,J,1);
[y_coord,x_coord]=meshgrid(1:I,1:J);
u_adv=1.6;
u=u_adv*cos(3*y_coord/J).*x_coord/(I/2)/1.5;
v=-1.2*u_adv*cos(2*14*x_coord/J);

% how to run the model
% assume you have two sources a t=0
To=zeros(I,J,1);
To(21:22,11:12,1)=5;
To(11:12,15:16,1)=5;
T = adv_diff_model(To,K,N,u,v);

% how to convert between gridded and array
T100 = T(:,:,100);
T_array = T100(:);  % how to convert from gridded to array
T100 = reshape(T_array, [30 30]); % how to convert from array to gridded

% how to make a plot of the model output
subplot(1,2,1); pcolor(x_coord, y_coord, To); shading flat; colorbar;
title 'DISPERSION PATTERN at t=0'
subplot(1,2,2); pcolor(x_coord, y_coord, T100); shading flat; colorbar;
title 'DISPERSION PATTERN at t=100'


% load true dispersion pattern at t=100
T_true = load('pattern1.dat');
T_true_gridded = reshape(T_true, [30 30]);
pcolor(x_coord, y_coord, T_true_gridded); shading flat; colorbar;
title 'TRUE DISPERSION PATTERN at t=100'

% below are the logical steps of the code you should produce

% generate the measurment error (the gaussian noise) with the choosen signal to 
% noise and add this to T_true to produce the synthetic observations, there is
% a MATLAB function randn.m
% e.g. T_obs = T_true + err;

% build the matrix of the adv_diff dynamics

% set your weights

% do the inversion

% re run adv-diff model with your estimate of the initial conditions

% check the posterior statistics

% do figures

return
% ---------------------------------------------------------------------------
% HINTS below
%
% HINT 1) how to add Gaussian noise with a given signal to noise ratio ?
%
% compute the STD of the true dispersion pattern, use the matlab function
% randn.m to generate Gaussian noise with zero mean and variance = 1, 
% renormalize the noise according to the signal to noise ration and the 
% STD of the true dispersion. Then add the noise to the true dispersion 
% data to generate your synthtic observations.
%
%
% HINT 2) how to build the matrix E, to solve the model y = Ex + n ?
%
% You know that sources occupy a square of 4 gridpoints with uniform value, 
% therefore your model parameters x will only be the inital values for each 
% 4 gridpoint area (which constitute the source or sources)
% This means that your matrix E will have dimensions of E(I*J, I*J/4)
% Therefore to generate E you need to run the model I*J/4 times, 
% each time initializing the model with a 4 gridpoint square impulse of 1.
% e.g. patch 1 
%      i=1; j=1;
%      To( i:i+1, j:j+1) =1; 
%      T = adv_diff_model((To,K,N,u,v);
% e.g. patch 2 
%      i=3; j=1;
%      To( i:i+1, j:j+1) =1; 
%      T = adv_diff_model((To,K,N,u,v);
%      ...... 
% then the array of the solution at t=100 for pathc 1 is going to be the first column of
% E and so on..
% Once you invert the system and find your estimate for the model parameters x, then
% you will need to remap the model parameters them to the initial condition in the same way 
% e.g. patch1
%      i=1; j=1;
%      To( i:i+1, j:j+1) =xhat(1); 
%      T = adv_diff_model((To,K,N,u,v);
% e.g. patch 2 
%      i=3; j=1;
%      To( i:i+1, j:j+1) =xhat(2); 
%      T = adv_diff_model((To,K,N,u,v);
%      ...... 
% ALSO NOTE: you need to compute E only one time, once you have it you can 
% use it for the entire exerc. becuase E is time independent.
% ---------------------------------------------------------------------------



% If you want to check your solutions below is the MATLAB code for HINT 1
% and HINT 2
% ---------------------------------------------------------------------------


% add Guassian noise with zero mean and a given signal to noise ratio 
% to the true dispersion data to generate some synthetic observations
signal_noise = 5;  
rms_T = std(T_true);
rms_err = rms_T / signal_noise;
err = randn( length(T_true), 1);
err = err*rms_err/std(err);
T_obs  = T_true + err;


%==========================================================
%	how to build E,  for the model y = E*x + n
% NOTE: becuase you know that sources occupy a square of 
% 4 gridpoints, we will build E by perturbing the model 
% initial conditions with impulses over 4 grid points
% square only.
%==========================================================
[I,J]=size(x_coord);
%load E
if ~exist('E')
   n=0;
   E = zeros(I*J, I*J/4);
   for i=1:2:I-1
      for j=1:2:J-1
         n=n+1;
         To=zeros(I,J,1);
         % perturb an inital patch of 4 gridpoints
         To(i:i+1,j:j+1)=1; 
         % run model
         T=adv_diff_model(To,K,N,u,v);
         % get time = 100
         T100=T(:,:,100);
         % save the column of E associated with the patch
         E(:,n)=T100(:);
      end
   end
end

% once you find your estimates of the model parameters x
% you can remap these on the model intial conditions with
% the reverse loop
n=0; To_hat=zeros(I,J,1);
for i=1:2:I-1
   for j=1:2:J-1
      n=n+1;
      To_hat(i:i+1,j:j+1)=xhat(n);
   end
end