clear;
clf;

mmax = 1000;
mint = mmax*3;
d_phase = 2*pi/mint;

ddsin = sin(d_phase*(1:mint));
ddcos = cos(d_phase*(1:mint));

% ++--++--++-- series
twist = 0;
AR(1) = 1;
AR(2) = 1;
for m = 3:mmax,
 AR(m) = -AR(m-2);
end;

% imaginary part in alternating parity
scale = 0.1;
parity = 1;
for m=1:mmax,
 AI(m) = parity*scale*AR(m);
 parity = -parity;
end;

% spin wave rotation 
twist = 8;
damax = 2*pi/twist;
phi(1) = 0;
for m = 1:mmax,
 BR(m) = AR(m)*cos(phi(m))-AI(m)*sin(phi(m));
 BI(m) = AR(m)*sin(phi(m))+AR(m)*sin(phi(m));
 phi(m+1) = phi(m)+damax*rand(1);
end;

order = (1:mmax).^0;

for ip = 1:mint,
 ipp = mod((1:mmax)*ip-1,mint)+1;
 dsin = ddsin(ipp).*order;
 fr(ip) = dsin*BR';
 fi(ip) = dsin*BI';
end;

plot(1:mint,fr);
title (['real part of force: ' num2str(mmax) ' channels, 2*pi/' num2str(twist) ' twist']);

%var(fr)

PRINT;
 
