%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulate the multi-channel method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clf;
clear;

% number of channels
nc = 100;

% number of eigenvalues
ne = 30*nc;

a1 = 0.1 * nc;
a2 = 0.4 * nc;
a3 = 0.5 * nc;

b1 = 0.7 * nc;
b2 = 0.9 * nc;

p1 = 0.5 * (a2-a1) / (a3-a1);

% generate the eigenvalues
for k = 1:ne,
% p = rand(1);
p = rand(1)*0.5+mod(k,2)*0.5;
if p < p1 
 we(k) = a1 + (a2-a1)*sqrt(p/p1);
elseif p < 0.5 
 we(k) = a3 - (a3-a2)*sqrt((0.5-p)/(0.5-p1));
elseif p < 0.75 
 we(k) = b1 + (b2-b1)*sqrt((p-0.5)/0.25)/2;
else
 we(k) = b2 - (b2-b1)*sqrt((1-p)/0.25)/2;
end;
end;
we = sort(we)';
hist(we,1:nc);
hold on;

% generate the A(m) series
% uniform +++++ series:
% A = ones(nc,1);
% random phase complex series:
% phi = 2*pi*rand(nc,1);
% A = cos(phi) + i*sin(phi);
% spin-wave series
spinwave;

% calculate the response on each channel

sin2w = sin(pi*we).^2;
we2 = we.^2;
ones_e = ones(ne,1);
for m=1:nc,
 rf(:,m) = 2*m/pi*ones_e ./ (ones_e*(m^2) - we2);
end;

for m=1:nc,
 r(m) = 0;
% no parity control 
for mp = 1:nc,
% parity control
% for mp = 2-mod(m,2):2:nc,
  r(m) = r(m) + sum(rf(:,m) .* rf(:,mp) .* sin2w) * A(mp);
 end;
 r(m) = r(m)/A(m);
 r(m) = abs(r(m));
 if m==1 
  acc(m) = 0.5 * r(m);
 else
  acc(m) = acc(m-1) + 0.5*(r(m)+r(m-1));
 end;
end;

plot(r,'o');
 
ss = input('continue... ');

clf;
plot(acc,'r');
hold on;
plot([0 we' nc], [0 0.5:ne-0.5 ne]);

fprintf('Total accumulant = %f \n', ne);
fprintf('Actual accumulant = %f \n', sum(r));
fprintf('Relative error = %f per cent\n\n', abs(ne-sum(r))/ne*100 );

fprintf('Halfway accumulant = %f\n', ne/2);
fprintf('Actual accumulant = %f \n', acc(a3+1));
fprintf('Relative error = %f per cent\n\n', abs(ne-2*acc(a3))/ne*100 );









