function [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(s_width, s_length, run, N_Ref, N_Cpt, Mu_beta, Mu_alpha, Bound_beta, Bound_alpha, MV, RV, var_n) 

rng(run+1306) % Random seed

%% Scene simulation 

phenLowerBound = 0.05;
phenUpperBound = 0.15;
n_pic = 15;
s_n = s_width*s_length; % Total number of areas in the scene
[xx,yy] = meshgrid((-1:2/(s_width-1):1),(-1:2/(s_length-1):1));
xxyy = cat(2,xx(:),yy(:));

g = zeros(1,s_n);
for pic = 1:n_pic
    mu = 2*(rand(1,2)-0.5);
    sig = diag([ phenLowerBound + (phenUpperBound - phenLowerBound)*abs(randn()+0.5) , phenLowerBound + (phenUpperBound - phenLowerBound)*abs(randn()+0.5) ]);
    g = g + mvnpdf(xxyy,mu,sig);
end

g = g-min(g);
g = .5*(g/max(g))+1e-5;
G_theo = [ones(s_n,1),g]; % Theoretical matrix G (see eq.(3) of [1])


%% Sensors simulation

F_theo = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
F_theo = [F_theo,[0;1]]; % Theoretical matrix F (see eq.(3) of [1])


%% Data simulation

X_theo = G_theo*F_theo; %  Theoretical matrix X (see eq.(2) of [1])

W = zeros(s_n,N_Cpt+1);

idx_Ref = randperm(s_n);
idx_Ref = idx_Ref(1:N_Ref); % Reference measurement locations
W(idx_Ref,end) = 1;

N_RV = round(N_Cpt*RV); % Nb. of sensors having a RendezVous
idx_CptRV = randperm(N_Cpt);
idx_CptRV = idx_CptRV(1:N_RV); % Selection of sensors having a RendezVous
idx_RefRV = randi(N_Ref,1,N_Cpt);
idx_RefRV = idx_Ref(idx_RefRV(1:N_RV)); % Selection of the references for each RendezVous

for i = 1 : N_RV
    W(idx_RefRV(i),idx_CptRV(i)) = 1;
end

N_data = round((1-MV)*(N_Cpt)*(100-N_Ref)); % Nb. of measurements in data matrix X
xCpt = 1 : s_n;
xCpt(idx_Ref) = []; % Reference free locations
[xx,yy] = meshgrid(xCpt,1:N_Cpt); % Possibly sensed locations

idx_data = randperm((s_n-N_Ref)*N_Cpt);
for i = 1 : N_data
    W(xx(idx_data(i)),yy(idx_data(i))) = 1; % Sensor measurement placement
end

N = var_n*randn(s_n,N_Cpt+1); % Noise simulation
N(:,end) = 0;
N = max(N,-X_theo);

X = W.*(X_theo+N); % Data matrix X

%% Calibration parameters

% % Common parameters
Omega_G = [ones(s_n,1),W(:,end)]; % Mask on known values in G (see eq.(14) of [1])
Omega_F = [zeros(2,N_Cpt),[1;1]]; % Mask on known values in F (see eq.(15) of [1])
Phi_G = [ones(s_n,1),X(:,end)]; % Known values in G (see eq.(14) of [1])
Phi_F = [zeros(2,N_Cpt),[0;1]]; % Known values in F (see eq.(15) of [1])
Ginit = abs(randn(s_n,2)+mean(Phi_G(idx_Ref,end))); % Initial matrix G : randn + mean of known ref values
Ginit = (1-Omega_G).*Ginit+Phi_G;
Finit = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
Finit = [Finit,[0;1]]; % Initial matrix F
Finit = (1-Omega_F).*Finit+Phi_F;

end