% script Local_Moves; 
%                    perform localized moves in a polymer chain.
% >> Perform moves of the chain that involve only 
%     a subset of atoms. Augment to avoid steric clashes.  <<
% PROBLEM 1: given chain, deform torsion angles so first and last atoms
%            approach to distance 1
% PROBLEM 2: further deform torsion angles to allow joining of first and
%            last atoms with correct bond angle
% PROBLEM 3: introduce arbitrary moves of closed chain, centered about any
%            given torsion angle. Find maximally distorting infinitesimal 
%            deformation (leading singular vector).
%
% For simplicity, consider a polymer where successive 
% bond-angles are fixed, equal to >> theta <<, while torsion
% angles phi_j are entirely free. 
%                                               7
%        ____________i+2                       /
%        |            ^                       /
%        |            |                      /\theta=pi-theta'
%        |    \ b(i)  |____n(i+1) --------->----------> x-axis
%        i-----\----->i+1                   \    b(i;0)
%       / ^   phi(i)   \                     \
%      /   \            \                     \
%     n(i)  \            \                     V z-axis
%           i-1------------         
%                                 q(i-1) * p(i)* b(i;0) = b(i)
%      q(i)   = q(i-1)   * p(i;0)
%      p(i;0) = d(i-1;0) * a(i;0)
%      d(i;0) =    Ux(phi(i))   = /     1           0         0     \
%                                 |     0        cos(phi) -sin(phi) |
%                                 \     0        sin(phi)  cos(phi) /
%      a(i;0) =    Uz(theta-pi) = / -cos(theta)  sin(theta)   0     \
%                                 | -sin(theta) -cos(theta)   0     |
%                                 \     0           0         1     /
%  Then, 
%     if cos(theta)=ct, sin(theta)=st, cos(phi)=cp(i), sin(phi)=sp(i),   
%        W(i) = /  -ct       st        0   -b*ct     \
%               |  -st*cp   -ct*cp   -sp   -b*st*cp  |
%               |  -st*sp   -ct*sp    cp   -b*st*sp  |
%               \    0        0        0       1     /
% This latter matrix works as follows:
% Since 
%        R(i) = R(i-1) + b(i)
%  and
%        b(i) = b*U(1)*U(2)...*U(i)*e1 
%  then  
%        R(i) = W(1)*...*W(i)(1:3,4)
%  since 
%        W(1)*...*W(i) = / p(1)*...*p(i)  p(1)*...*p(i)*e1+R(i) \
%                        \       0                  1          /
%      =  / Q(1-i)  R(i) \
%         \   0      1   /
%
%      n(i)   = q(i-1)*d(i-1;0)*e3
%      b(i)   = q(i)*e1 = q(i-1)*d(i-1;0)*a(i;0)*e1
%      b(0)   = e1 (unused arm on Atom 1 afixed to the x-axis;
%                   other arm at a(1)*e1, d(0) = I )
%      phi(1), phi(N) may be defined but have no effect                        
%                                  
%                                 |sin(phi(i))| = ||n(i)Xn(i+1)||
%    R(i+1) = R(i) + b(i)
%    n(i+1) = b(i) X b(i+1)
%
% CLOSURE CONDITIONS
%
%  I.               R(1)=R(N+1)
% II.               b(N)=e1
%      
%    
% To place the chain given the dihedral angles phi(1:n)
% proceed from atoms 1 to N:
%   
%   
% Linear algorithm, scales like 139*N; can we cut down on the constant?
clear all
[a,b,c]=sphere(10);
N = 10;
omega =   pi/6;
theta = 2*pi/3;
%r1=0;r2=0;r3=0;
R=eye(4);
figure
% Build Initial Configuration as a Spiral Staircase
flops(0);
for i = 1:N
   phi(i) = omega;
   Q = affmat(phi(i),theta,1);
   R=R*Q;
   X(i) = R(1,4);Y(i) = R(2,4); Z(i) = R(3,4);
end
a1=flops
hold on
for i = 1:N
   surf(.1*a+X(i),.1*b+Y(i),.1*c+Z(i));
end
plot3(X,Y,Z,'y-')
hold off
for kk = 1:40
   pause(1)
   view(6*kk,30)
end
