function [T_true, To]=hw3_get_sources(source);

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);


if source == 1
    To=zeros(I,J,1);
    To(21:22,19:20,1)=5;
    To(11:12,15:16,1)=0;
end
if source == 2
    % assume you know there are now sources on the boundary
    To=zeros(I,J,1);
    To(11:12,21:22,1)=0;
    To(9:10,7:8,1)=5;
    % 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.3*y_coord/J).*x_coord/(I/2)/1.5;
    v=-1.2*u_adv*cos(1.4*14*x_coord/J);
end
if source == 3
    % assume you know you have two sources, determine which one
    % is likely to be bigger based on the correlation coefficient
    % and based on the sum of the residual square
    To=zeros(I,J,1);
    To(23:24,21:22,1)=5;
    To(9:10,7:8,1)=3;
end

% how to run the model
T = adv_diff_model(To,K,N,u,v);
T100 = T(:,:,100);
T_true=T100(:);