function A=coup_osc_oper(wsq,iBC);

% Generate matrix operator for N ACTUAL particles. 
% BCs: 
% Rigid BC: iBC=0 and length(wsq) MUST be N+1
% |--wsq(1)--O--wsq(2)--O-- ... --wsq(N)--O--wsq(N+1)--|
%            1          2                 N         <-- particle number

% Periodic BC: iBC=1 and length(wsq) MUST be N
% O--wsq(1)--O--wsq(2)--O-- ... --wsq(N)--O--wsq(1)--O
% N          1          2                 N          1  <-- particle number
% |                                                  |
% Ficticious                                        Ficticious
% Particle N is connected to particle 1 by spring wsq(1).

% "Free stress" BC: iBC=2 and length(wsq) MUST be N-1
% O--wsq(2)--O-- ... --wsq(N)--O
% 1          2                 N   <-- particle number
% This is same as rigid BC with wsq(1)=wsq(N+1)=0

if nargin==1 
  iBC=0;
end

if iBC==1 % periodic BC
  wsq=[wsq;wsq(1)]; 
end

if iBC==2 % "free stress" BC
  wsq=[0;wsq;0]; % add fictitious springs with zero spring constant
end

N=length(wsq)-1

A=sparse(N,N);
A=spdiags(-(wsq(1:end-1)+wsq(2:end)),0,N,N);
A=spdiags(wsq(2:end),-1,A);
A=spdiags(wsq,1,A);

if iBC==1 % fix for periodic BC
  A(1,end)=wsq(1);
  A(end,1)=wsq(1);
end

