diff --git a/projects/project_populationvector/populationvector.tex b/projects/project_populationvector/populationvector.tex
index 44c1b62..d619410 100644
--- a/projects/project_populationvector/populationvector.tex
+++ b/projects/project_populationvector/populationvector.tex
@@ -55,6 +55,8 @@
   \texttt{spikes} variables of the \texttt{population*.mat} files.
   The \texttt{angle} variable holds the angle of the presented bar.
 
+NOTE: the orientation is angle plu 90 degree!!!!!!
+
   \begin{parts}
     \part Illustrate the spiking activity of the V1 cells in response
     to different orientation angles of the bars by means of spike
@@ -97,6 +99,9 @@
   \end{parts}
 \end{questions}
 
+NOTE: change data generation such that the phase variable is indeed
+the maximum response and not the minumum!
+
 \end{document}
 
 gains and angles of the 6 neurons:
diff --git a/projects/project_populationvector/solution/populationvector.m b/projects/project_populationvector/solution/populationvector.m
index 829e07d..134b8ed 100644
--- a/projects/project_populationvector/solution/populationvector.m
+++ b/projects/project_populationvector/solution/populationvector.m
@@ -16,6 +16,7 @@ end
 close all
 cosine = @(p,xdata)0.5*p(1).*(1.0-cos(2.0*pi*(xdata/180.0-p(2))));
 files = dir('../unit*.mat');
+phases = zeros(length(files), 1);
 figure()
 for j = 1:length(files)
 	file = files(j);
@@ -31,7 +32,11 @@ for j = 1:length(files)
 	p0 = [mr, angles(maxi)/180.0-0.5];
 	%p = p0;
 	p = lsqcurvefit(cosine, p0, angles, rates');
-	phase = p(2)*180.0
+	phase = p(2)*180.0 + 90.0;
+    if phase > 180.0
+        phase = phase - 180.0;
+    end
+    phases(j) = phase;
 	subplot(2, 3, j);
 	plot(angles, rates, 'b');
 	hold on;
@@ -41,3 +46,73 @@ for j = 1:length(files)
 	ylim([0.0 50.0])
 	title(sprintf('unit %d', j))
 end
+
+
+%% read out:
+a = load('../population04.mat');
+spikes = a.spikes;
+angle = a.angle;
+unitphases = a.phases*180.0 + 90.0;
+unitphases(unitphases>180.0) = unitphases(unitphases>180.0) - 180.0;
+figure();
+subplot(1, 3, 1);
+angleestimates1 = zeros(size(spikes, 2), 1);
+angleestimates2 = zeros(size(spikes, 2), 1);
+for j = 1:size(spikes, 2)
+    rates = zeros(size(spikes, 1), 1);
+    for k = 1:size(spikes, 1)
+        r = firingrate(spikes(k, j), 0.0, 0.2);
+        rates(k) = r;
+    end
+    [x, inx] = sort(phases);
+    plot(phases(inx), rates(inx), '-o');
+    hold on;
+    angleestimates1(j) = popvecangle(phases, rates);
+    [m, i] = max(rates);
+    angleestimates2(j) = phases(i);
+end
+hold off;
+subplot(1, 3, 2);
+hist(angleestimates1);
+subplot(1, 3, 3);
+hist(angleestimates2);
+angle
+mean(angleestimates1)
+mean(angleestimates2)
+
+
+%% read out robustness:
+files = dir('../population*.mat');
+angles = zeros(length(files), 1);
+e1m = zeros(length(files), 1);
+e1s = zeros(length(files), 1);
+e2m = zeros(length(files), 1);
+e2s = zeros(length(files), 1);
+for i = 1:length(files)
+	file = files(i);
+	a = load(strcat('../', file.name));
+    spikes = a.spikes;
+    angle = a.angle;
+    angleestimates1 = zeros(size(spikes, 2), 1);
+    angleestimates2 = zeros(size(spikes, 2), 1);
+    for j = 1:size(spikes, 2)
+        rates = zeros(size(spikes, 1), 1);
+        for k = 1:size(spikes, 1)
+            r = firingrate(spikes(k, j), 0.0, 0.2);
+            rates(k) = r;
+        end
+        angleestimates1(j) = popvecangle(phases, rates);
+        [m, inx] = max(rates);
+        angleestimates2(j) = phases(inx);
+    end
+    angles(i) = angle;
+    e1m(i) = mean(angleestimates1);
+    e1s(i) = std(angleestimates1);
+    e2m(i) = mean(angleestimates2);
+    e2s(i) = std(angleestimates2);
+end
+figure();
+subplot(1, 2, 1);
+scatter(angles, e1m);
+subplot(1, 2, 2);
+scatter(angles, e2m);
diff --git a/projects/project_populationvector/solution/popvecangle.m b/projects/project_populationvector/solution/popvecangle.m
new file mode 100644
index 0000000..c06a773
--- /dev/null
+++ b/projects/project_populationvector/solution/popvecangle.m
@@ -0,0 +1,13 @@
+function angle = popvecangle(phases, rates)
+% estimate population vector
+vecs = zeros(2, length(phases));
+norm = sum(rates);
+vecs(1, :) = rates.*cos(2*pi*phases/180.0)/norm;
+vecs(2, :) = rates.*sin(2*pi*phases/180.0)/norm;
+mvec = mean(vecs, 2);
+angle = atan2(mvec(2), mvec(1));
+angle = angle/2/pi*180.0;
+if angle < 0
+    angle = angle + 180.0;
+end
+end