function [xy_mean,xy_var] = monte_carlo(fun,rp_true,R,N)
% MONTE_CARLO - conducts a Monte Carlo simulation on the function
% 'fun' with the true state 'rp_true' and covariance 'R'.
%   
% xy_mean - mean of the transformed state
% xy_var  - variance of the transformed state
%
  
what = 'plain';  
switch what
  case 'plain'
   for mc = 1:N
     rp_meas = rp_true + sqrt(R)*randn(size(rp_true));
     xy_meas(:,mc) = feval(fun,rp_meas);
   end
 case 'debiased'
   for mc = 1:N
     rp_meas = rp_true + sqrt(R)*randn(size(rp_true));
     xy_meas(:,mc) = debias_lerro(rp_meas,R);
   end
end

  xy_mean = mean(xy_meas,2);
  xy_var  = cov(xy_meas');