% three3.m

% Adapted from Luke's code
% NEUR 1680, Spring 2008

clear

w = -pi/2:0.005:pi/2;                    % wind directions
n_trials = 1000;                        % number of trials per wind didrection
n_w = length(w);                        % number of wind directions
errors = zeros(1,length(w));            % MSE for each direction
pref = [pi/4, 3*pi/4, 5*pi/4, 7*pi/4];  % preferred direction for each neuron

pref_repeat = repmat(pref,n_trials,1);  % will be useful

positive_part = @(s) s.*(s>0);          % defining an anonymous function

for n = 1:n_w
    
    wind = w(n);
    wind_repeat = repmat(wind,n_trials,4);
    responses = 50*cos(wind_repeat - pref_repeat);
    responses = positive_part(responses);
    responses = responses + 5*randn(n_trials,4);
    responses = positive_part(responses);
    population_vector_x = responses.*cos(pref_repeat);
    population_vector_y = responses.*sin(pref_repeat);
    x = sum(population_vector_x,2);
    y = sum(population_vector_y,2);
    [theta,rho] = cart2pol(x,y);
    recovered = theta;
    error = sqrt(mean((mean(wind_repeat,2)-recovered).^2));
    errors(n) = error;
    
end

figure(1)
scatter(rad2deg(w),rad2deg(errors),15,'r','s','filled');
grid
title('Population-Vector Simulation of Cercal System','FontSize',15)
xlim([-95 95]);
xlabel('Wind Direction (Degrees)','FontSize',15)
ylabel('MSE (Degrees)','FontSize',15)