function [] = frames(shapes)  
  
%This function constructs a set of DNA base reference points and frames
%corresponding to the Tsukuba convention 
%from a given shape vector of the dimension Nx1, N=12*nbp-6, 
%where nbp is the number of base pairs in the DNA sequence.
%The parameters in the vector 'shapes' have to be ordered by basepairs 
%along the sequence, and for every basepair the order is:
%Intra rotations (Buckle, Propeller, Opening), Intra translations 
%(Shear, Stretch, Stagger), Inter rotations (Tilt, Roll,Twist) 
%and Inter translations (Shift, Slide, Rise). E.g. the entry shapes((i-1)*12+1) 
%corresponds to Buckle of the i-th basepair for i=1...nbp, and the entry 
%shapes((i-1)*12+7) contains the Tilt value for the i-th basepair
%step, i=1..nbp-1. 
%
%The output of this function can be chosen to be saved in one of the two text formats, 
%specified by the variable 'format'. 
%For the first choice, corresponding to 'format=1', 
%the reference points and frames for each base are printed as a 4x3 matrix,
%with first row being the coordinates of the reference point, and the each
%row below being a reference frame vector d_i, i=1...3.
%The bases are ordered from top to bottom, first all the bases of a reference strand,   
%then all the bases of a complementary strand. Example: file base_frames1.bfr
%
%For the second choice, corresponding to 'format=2', 
%the reference points and frames for each base are printed as a line of 14 integers,
%the first being the index of a strand (1-reference strand, 2-complementary strand),
%second being the index of a base pair, the next nine being the coordinates
%of reference frame vectors d_i, i=1...3, multiplied by 1000, and the last
%three being the coordinates of the reference point multiplied by 1000. 
%The bases are ordered from top to bottom, first all the bases of a reference strand,   
%then all the bases of a complementary strand. Example: file base_frames2.fra
%
%Setting 'format=0' results in no text output. However the results will still be saved
%in the file absolute_coord.mat 
%
%References: O. Gonzalez, D. Petkeviciute and J.H. Maddocks, A sequence-dependent
%rigid-base model of DNA, J. Chem. Phys 138 (5), 2013. 

format = 1; 
nbp = (size(shapes,1)+6)/12;

% absolute coordinates of the first base pair 
basepair(1).G = eye(3);
basepair(1).q =[0, 0, 0]'; 

% relative coordinates of the oligomer
[eta w u v] = vector2shapes(shapes);

for i=1:nbp
    
   % complimentary strand 
   basepair(i).Dc = basepair(i).G*(sqrtm(rodrigues(eta(i,:))))'; 
   basepair(i).rc = basepair(i).q - 0.5*basepair(i).G*w(i,:)';
   
   % main strand
   basepair(i).D = basepair(i).Dc*rodrigues(eta(i,:));
   basepair(i).r = basepair(i).rc + basepair(i).G*w(i,:)';
   
   if i<nbp
       
        % base pair 
        basepair(i+1).G = basepair(i).G*rodrigues(u(i,:));
        H = basepair(i).G*sqrtm(rodrigues(u(i,:)));
        basepair(i+1).q = basepair(i).q + H*v(i,:)';
   end  
end

%printing base coordinates to a text file
if (format==1)
    print_frames1(basepair,nbp);
elseif (format==2) 
    print_frames2(basepair,nbp);
end    

%saving 
save absolute_coord basepair; 
    
end

%-------------------------------------------------------------

function [Q] = rodrigues(k)

  k = k/5;
  X = [0 -k(3) k(2); k(3) 0 -k(1); -k(2) k(1) 0];
  Q = (eye(3)+0.5*X)*inv(eye(3)-0.5*X);
  
end

%-------------------------------------------------------------

function [eta,w,u,v] = vector2shapes(y)

    nn = size(y);
    n = (nn(1)+6)/12;
    eta = zeros(n,3);
    w = zeros(n,3);
    u = zeros(n-1,3);
    v = zeros(n-1,3);
    
    for i = 1:n
        j = (i-1)*12;
        eta(i,:) = y(j+1:j+3)'; 
        w(i,:) = y(j+4:j+6)';
        if (i<n)
            u(i,:) = y(j+7:j+9);
            v(i,:) = y(j+10:j+12);
        end;    
    end;    

end

%-------------------------------------------------------------

function y = shapes2vector(eta,w,u,v)
    
    nn = size(eta);
    for i = 1:nn(1)
        j = (i-1)*12;
        y(j+1:j+3) = eta(i,:);
        y(j+4:j+6) = w(i,:);
        if (i<nn(1))
            y(j+7:j+9) = u(i,:);
            y(j+10:j+12) = v(i,:);
        end;    
    end;    
    y = y';
    
end

%-------------------------------------------------------------

function [] = print_frames1(basepair,nbp)

    s = fopen('base_frames1.bfr','w');
    fprintf(s,'%d \n' ,1);
       
    for j=1:nbp
       %main strand
       fprintf(s,' %d   %d \n' ,1,j);
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).r(1), basepair(j).r(2), basepair(j).r(3));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).D(1,1), basepair(j).D(2,1), basepair(j).D(3,1));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).D(1,2), basepair(j).D(2,2), basepair(j).D(3,2));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).D(1,3), basepair(j).D(2,3), basepair(j).D(3,3));
    end
    
    for j=1:nbp   
       %complementary strand
       fprintf(s,' %d   %d \n' ,2,j);
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).rc(1), basepair(j).rc(2), basepair(j).rc(3));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,basepair(j).Dc(1,1), basepair(j).Dc(2,1), basepair(j).Dc(3,1));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,-basepair(j).Dc(1,2), -basepair(j).Dc(2,2), -basepair(j).Dc(3,2));
       fprintf(s,'  %10.6f  %10.6f  %10.6f \n' ,-basepair(j).Dc(1,3), -basepair(j).Dc(2,3), -basepair(j).Dc(3,3));     
    end
    
    fclose(s);
    
end

%-------------------------------------------------------------

function [] = print_frames2(basepair,nbp)

    s = fopen('base_frames2.fra','w');
         
    for j=1:nbp
       %main strand
       fprintf(s,' %d   %2d ' ,1,j);
       fprintf(s,' %6.0f %6.0f %6.0f' ,basepair(j).D(1,1)*1000, basepair(j).D(2,1)*1000, basepair(j).D(3,1)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f' ,basepair(j).D(1,2)*1000, basepair(j).D(2,2)*1000, basepair(j).D(3,2)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f' ,basepair(j).D(1,3)*1000, basepair(j).D(2,3)*1000, basepair(j).D(3,3)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f \n' ,basepair(j).r(1)*1000, basepair(j).r(2)*1000, basepair(j).r(3)*1000);
    end
    
    for j=1:nbp      
       %complementary strand
       fprintf(s,' %d   %2d ' ,2,j);
       fprintf(s,' %6.0f %6.0f %6.0f' ,basepair(j).Dc(1,1)*1000, basepair(j).Dc(2,1)*1000, basepair(j).Dc(3,1)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f' ,-basepair(j).Dc(1,2)*1000, -basepair(j).Dc(2,2)*1000, -basepair(j).Dc(3,2)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f' ,-basepair(j).Dc(1,3)*1000, -basepair(j).Dc(2,3)*1000, -basepair(j).Dc(3,3)*1000);
       fprintf(s,' %6.0f %6.0f %6.0f \n' ,basepair(j).rc(1)*1000, basepair(j).rc(2)*1000, basepair(j).rc(3)*1000);      
    end
 
    fclose(s);
    
end
