% 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(i)=pi-theta'(i)
%        |    \ 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+1) = R(i) + 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) = / U(1)*...U(i)  bU(1)*...*U(i)*e1+R(i) \
%                        \       0                  1           /
%      =  / Q(i)  R(i+1) \
%         \   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)
%
%  here we implement this slightly differently:
%  Form B = [Q(1) Q(2) Q(3) ... Q(N)] e1
%   and R = [B(1) S_1^2B(i) S_1^3B(i) ... S_1^NB(i)];
%   compute flops to compare (cheaper by over twice!)
% 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:
%   
%   
% BAD ALGORITHM: it scales like N^2 ...
clear all
close all
view(25,0)
[a,b,c]=sphere(10);
N = 4
for t = 1:240
   %N = 10;
   omega =   pi/120;
   theta = 2*pi/3;
   %r1=0;r2=0;r3=0;
   R=zeros(3,N);
   B=zeros(3,N);
   B(1,:)=1;
   % Build Initial Configuration as a Spiral Staircase
   % count flops for construction of B, R
   %flops(0);
   phi0 = 0;
   Q = rotmat(phi0,theta);
   B(:,1) = Q(:,1); 
   l1 = 1;
   for l = 2:N  
      ct=cos(theta); st=sin(theta); 
      phi(l1) = omega*l1^2*t;
      cp=cos(phi(l1));sp=sin(phi(l1));
      DI = [ cp -sp ; sp  cp ]; % 2x2 rotation abt x (in yz plane)
      BA = [-ct  st ;-st -ct ]; % 2x2 rotation abt z (in xy plane)
      Q(:,2:3) = Q(:,2:3)*DI;  % same as Q = Q*/1 0 \
      %                                        \0 DI/ rot abt x-axis
      Q(:,1:2) = Q(:,1:2)*BA;  % same as Q = Q*/BA 0\
      %                                        \0  1/ rot abt z-axis                        
      B(:,l) = Q(:,1);
      l1=l;
   end
   R=B;
   k1=1;
   for k = 2:N
      R(:,k)=R(:,k1)+R(:,k);
      k1=k;
   end
   %a0(N)=flops;
%end
%plot(1:N,a0(1:N))
%figure
%clf
hold on
axis([-2.5,2.5,-2.5,2.5,-2.5,2.5])
for i = 1:N
   surf(.05*a+R(2,i),.05*b+R(3,i),.05*c+R(1,i));
end
plot3(R(2,:),R(3,:),R(1,:),'y-')
hold off
pause(1)
end
%for kk = 1:40
%   pause(1)
%  view(6*kk,50)
%end
% 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(i)=pi-theta'(i)
%        |    \ 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+1) = R(i) + 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) = / U(1)*...U(i)  bU(1)*...*U(i)*e1+R(i) \
%                        \       0                  1           /
%      =  / Q(i)  R(i+1) \
%         \   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)
%
%  here we implement this slightly differently:
%  Form B = [Q(1) Q(2) Q(3) ... Q(N)] e1
%   and R = [B(1) S_1^2B(i) S_1^3B(i) ... S_1^NB(i)];
%   compute flops to compare (cheaper by over twice!)
% 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:
%   
%   
% BAD ALGORITHM: it scales like N^2 ...
clear all
close all
view(25,0)
[a,b,c]=sphere(10);
N = 4
for t = 1:240
   %N = 10;
   omega =   pi/120;
   theta = 2*pi/3;
   %r1=0;r2=0;r3=0;
   R=zeros(3,N);
   B=zeros(3,N);
   B(1,:)=1;
   % Build Initial Configuration as a Spiral Staircase
   % count flops for construction of B, R
   %flops(0);
   phi0 = 0;
   Q = rotmat(phi0,theta);
   B(:,1) = Q(:,1); 
   l1 = 1;
   for l = 2:N  
      ct=cos(theta); st=sin(theta); 
      phi(l1) = omega*l1^2*t;
      cp=cos(phi(l1));sp=sin(phi(l1));
      DI = [ cp -sp ; sp  cp ]; % 2x2 rotation abt x (in yz plane)
      BA = [-ct  st ;-st -ct ]; % 2x2 rotation abt z (in xy plane)
      Q(:,2:3) = Q(:,2:3)*DI;  % same as Q = Q*/1 0 \
      %                                        \0 DI/ rot abt x-axis
      Q(:,1:2) = Q(:,1:2)*BA;  % same as Q = Q*/BA 0\
      %                                        \0  1/ rot abt z-axis                        
      B(:,l) = Q(:,1);
      l1=l;
   end
   R=B;
   k1=1;
   for k = 2:N
      R(:,k)=R(:,k1)+R(:,k);
      k1=k;
   end
   %a0(N)=flops;
%end
%plot(1:N,a0(1:N))
%figure
%clf
hold on
axis([-2.5,2.5,-2.5,2.5,-2.5,2.5])
for i = 1:N
   surf(.05*a+R(1,i),.05*b+R(2,i),.05*c+R(3,i));
end
plot3(R(1,:),R(2,:),R(3,:),'y-')
hold off
pause(1)
end
%for kk = 1:40
%   pause(1)
%  view(6*kk,50)
%end
%function qc = quatconj(q)
% produce conjugate quaternion
% if a scalar , b row vector (1x3 array)
% then    q  = (a, b)
%         qc = (a,-b)
function qc = quatconj(q)
         qc = [q(1), -q(2:4)];function q = quaternion(vector)
%function q = quaternion(vector)
%inputs:  vector  (3X1 array)
%output    q       1X4 vector (pure quaternion)
  q = [0, vector'];
  
  function w = quatmul(u,v)
%inputs: u,v quaternions to be multiplied   (1x4 arrays)
%output:  w = u^v their quaternion product  (1x4 array)
%   if    u = [u0, U], U=[u1,u2,u3]
%         v = [v0, V], V=[v1,v2,v3]
% then
%        W = [ u0*v0-U.V, u0*U+v0*V+UXV]
% where:   * is scalar multiplication
%          . is the dot product of U and V (=u1*v1+u2*v2+u3*v3)
%          X is the cross product of U and V (intrinsic)
%   UXV= [u2*v3-u3*v2,u3*v1-u1*v3,u1*v2-u2*v1]
% humor matlab's idiocyncracies by adding 1 to all indices
   U  = u(2:4);   V = v(2:4); 
   u0 = u(1);    v0 = v(1);
   w0 = u0*v0-dot(U,V);
   W  = u0*V+v0*U + cross1(U,V);
   w  = [w0,W];
   

function p = cross1(u,v)
%  p,a,b are 1X3 arrays
p = [u(2)*v(3)-u(3)*v(2),u(3)*v(1)-u(1)*v(3),u(1)*v(2)-u(2)*v(1)]; 
% Note that cross is an intrinsic function in matlab! (help cross)%function y = rotate(x,axis,angle)
%rotate vector x about axis by angle
% x, y, axis are vextors (3X1 arrays), angle is real scalar
function y = rotate(x,axis,angle)
   q  = rotator(axis,angle);
   qc = quatconj(q);
   qx = quaternion(x);
   qx = quatmul(q ,qx);
   qx = quatmul(qx,qc);
   y  = qx(2:4)';
   

 
function q = rotator(axis,angle)
%function q = rotator(axis,angle)
%inputs:  angle, real number
%          axis,   3X1 vector
%output    q       1X4 vector (unit quaternion)
t2 = angle/2;
cs = cos(t2); si = sin(t2);
nv = norm(axis);
  if nv == 0
      printf('need a nonzero axis vector \n')
      return
  end
  axis = axis/nv;
  q = si*axis';
  q = [cs, q];
  
  