function [netMoment, strainEnergy, rmsTraction, orientation, tx, ty] = traction_boussinesq(x,y,ux,uy,pois,young); 
% Constrained, Unconstrained, labdata, Traction matrices

% FTTC calculates the tractions, using both the Unconstrained
%	and the Constrained FTTC.
%	The required input is a file of displacements (4 columns: x-coordinates,
%	y-coordinates, x-displacements, y-displacements, all in pixels), 
%	pixel to micron conversion factor, Young's modulus and Poisson's ratio 
%	of the gel, and the cell boundary file (it is needed for the 
%	Constrained FTTC only).
% 
%   Xavier Trepat 2007
%	Iva Marija Tolic-Norrelykke 03-22-01

%  CENTER DISPLACEMENTS, PIXELS TO MICRONS, SPACING
%{
%put even displacements into column
DispX  	= 	mat2vertcol(DispX);
DispY  	= 	mat2vertcol(DispY); 
xorig 	= 	mat2vertcol(x);
yorig 	= 	mat2vertcol(y); 

xv  = (xorig - mean(xorig));
yv  = (yorig - mean(yorig));
uxv = (DispX);
uyv = (DispY);
 %}  
%original code started off with this line
xvorig = x;
yvorig = y;

spacing = y(2) - y(1);

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%     UNCONSTRAINED FTTC
%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%====================================================================================
%  N, RESHAPE DISPLACEMENTS
%below code originally in but commneted out and continued here

if(mod(length(ux),2)),
    N = length(ux)-1;
else
    N   = length(ux);
end

x  = x(1:N,1:N);
y  = y(1:N,1:N);
ux  = ux(1:N,1:N);
uy  = uy(1:N,1:N);
%}
%{
%below code replaced above code
N   = sqrt(length(xv));

 x  = reshape( xv,N,N);
 y  = reshape( yv,N,N);
ux  = reshape(uxv,N,N);
uy  = reshape(uyv,N,N);
%}
%====================================================================================
%  kx, ky, AND  k_abs 

clear mala*;
for i = 1 : (N/2)
    malax(i,:) = 0 : ((N/2)-1);                     
    malay(i,1:(N/2)) = (N/2)-i;                   
end;

kx = [ malax  malax-(N/2);  malax  malax-(N/2) ];
ky = [ malay-(N/2)  malay-(N/2);  malay  malay ];
ky = flipud(ky);
kx(:,(N/2+1)) =  kx(:,(N/2+1));                     
ky((N/2+1),:) =  ky((N/2+1),:);

k_abs = sqrt(kx.^2 + ky.^2);

%====================================================================================
%  ALPHA

alpha = atan2(ky,kx);
if kx(1,1) == 0 & ky(1,1) == 0,
    alpha(1,1) = 1.57080;
end;

%====================================================================================
%  C AND D

Cx = ((k_abs * young) / (2 * (1 - pois^2))) .* (1 - pois + pois .* (cos(alpha)).^2); % note + sign is correct but inconsistent with Butler 2002
Cy = ((k_abs * young) / (2 * (1 - pois^2))) .* (1 - pois + pois .* (sin(alpha)).^2);
D  = ((k_abs * young) / (2 * (1 - pois^2))) .* (pois .* sin(alpha) .* cos(alpha));

D(:,(N/2+1)) = zeros(N,1);
D((N/2+1),:) = zeros(1,N);

%====================================================================================
%  CALCULATE THE TRACTIONS

Dx = fft2(ux * 2 * pi / (N * spacing));
Dy = fft2(uy * 2 * pi / (N * spacing));

Tx = Cx .* Dx  +  D  .* Dy;
Ty = D  .* Dx  +  Cy .* Dy;

tx = real(ifft2(Tx));
ty = real(ifft2(Ty));

%====================================================================================
% RMS TRACTION, THETA, TRACE, U

rmst_onestep = sqrt(mean2(tx.^2 + ty.^2));

clear i;
Dxx = imag(Tx(1,2)) * (N * spacing / (2*pi)) * spacing^2; 
Dyy = imag(Ty(2,1)) * (N * spacing / (2*pi)) * spacing^2; 
Dxy = imag(Ty(1,2)) * (N * spacing / (2*pi)) * spacing^2; 
Dyx = imag(Tx(2,1)) * (N * spacing / (2*pi)) * spacing^2; 

[a_sym, lambda, theta] = rotate1([Dxx Dxy; Dyx Dyy]);

Trace_moment = Dxx + Dyy;

Uonestep = 0.5 * sum(sum([tx ty] .* [ux uy])) * spacing^2;

interior = find(x ~= min2(x) & x ~= max2(x) & y ~= min2(y) & y ~= max2(y));
Uonestepint = 0.5 * sum(sum([tx(interior) ty(interior)].*[ux(interior) uy(interior)])) * spacing^2;

%====================================================================================
% FINAL RESULT OF UNCONSTRAINED METHOD

rmsTraction  = rmst_onestep; % RMS tractions (Pa)
orientation  = theta*180/pi; % Orientation principale tractions (º)
netMoment    = -Trace_moment * (1e12); % Net moment
strainEnergy = Uonestepint * (1e12); % Total strain Energy

disp('rms Traction=');
disp(rmsTraction);
disp('Net moment=');
disp(netMoment);
disp('strain energy=');
disp(strainEnergy);
disp('orientation=');
disp(orientation);

