function [netMoment, strainEnergy, rmsTraction, orientation, xx, yy txx,tyy] = traction_boussinesq_constrained(x,y,DispX,DispY,pois,young,FilePath,k,Rot,FilterFreq,Filter); 
%   Calculates the constrained tracdtions using FTTC
%    
%   Xavier Trepat 2007
%	Iva Marija Tolic-Norrelykke 03-22-01
%modified by Robert 2013

%{
pois=Settings.Poisson;
young = Settings.Young;
Rot = Settings.CalcRot;
FilePath = File.pathname;
%}
%====================================================================================
%  CONSTANTS

a1   = (1.0 + pois) * (1.0 - pois) / (pi * young);
b1   = (1.0 + pois) * pois / (pi * young);
c1   = (1.0 + pois) * pois / (pi * young);

%====================================================================================
%displacements evened before putting into traction subroutine

%====================================================================================
%  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);

xvorig = xv;
yvorig = yv;

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

%put x,y,ux,uy in matrix form

N   = sqrt(length(xv));

x  = reshape( xv,N,N);
y  = reshape( yv,N,N);
ux  = reshape(uxv,N,N);
uy  = reshape(uyv,N,N);


%====================================================================================  
            % filter displacments
            if(strcmpi(Filter,'LP') == 1)%='LP'
            %if(Settings.Filter.Type == 'LP')
                ux = xfilter(ux, FilterFreq);
                uy = xfilter(uy, FilterFreq);
            
            else
                
                ux = wiener2(ux, [3 3]);
                uy = wiener2(uy, [3 3]);
            end
            
                       
%====================================================================================
%  N, RESHAPE DISPLACEMENTS

%code below originally in routine, but commented out
%{
% make sure there is an even number of rows and columns
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);

% Prepare for zero padding.
N = 2*N;
%}

%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
%     CONSTRAINED FTTC
%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

%====================================================================================
%  N, RESHAPE DISPLACEMENTS

N = sqrt(length(xv)) * 2;

iteration = 100; 
the_same  = 1e-6;
   
ux    = ux * 2 * pi / (N * spacing);
uy    = uy * 2 * pi / (N * spacing);

[x,y]  = meshgrid([(min(xv) - spacing*N/4 : spacing : min(xv)), x(1,2:size(x,2)-1),...
         (max(xv) : spacing : max(xv) + spacing*N/4)],...
			[(min(yv) - spacing*N/4 : spacing : min(yv)), y(2:size(y,1)-1,1)',...
   		(max(yv) : spacing : max(yv) + spacing*N/4)]);
ux = [zeros(N/4,N); zeros(N/2,N/4), ux, zeros(N/2,N/4); zeros(N/4,N)];   
uy = [zeros(N/4,N); zeros(N/2,N/4), uy, zeros(N/2,N/4); zeros(N/4,N)];   
   
xv     =  reshape(x,N^2,1);
yv     =  reshape(y,N^2,1);
uxv    =  reshape(ux,N^2,1);
uyv    =  reshape(uy,N^2,1);

%====================================================================================
%  CELL BOUNDARY
%      start= (findstr(d2b(k).name,'Cell'));		% name has to contains 'Cell'
%      ed=    findstr(d2b(k).name,'disp');		% name has to cotains 'disp'
%             m=d2b(k).name(start+4:ed-1);
% where m is the number of the cell folder

% eval(['cd C:\C2C12\croppeddata\' cellfolder ]); pwd;
 eval(['cd ' coreDir '/croppeddata/' cellfolder ]); pwd;
% %eval(['cd /Users/chardin/Desktop/research/LAM/' datedir '/croppeddata/' cellfolder ]); pwd;
%     
%     %  filename = 'cellboundary1.mat';
%     
%  for cellbound = 1: length(min_time)
%  
%     % if str2num(savenumber) >= min_time(cellbound) & str2num(savenumber) <= max_time(cellbound)
%      filename=strcat('cellboundary',num2str(cellbound),'.mat');
%      %filename1 = strcat('Magbeadpos',num2str(cellbound),'.mat');     
%============================================================================================
eval(['cd ' coreDir]); 
load( ['croppeddata\' cellfolder '/' 'cellboundary1.mat'], '-mat');    

   xrub = (xrub - 1*mean(displ(:,1))) * pixelsize;
   yrub = (yrub - 1*mean(displ(:,2))) * pixelsize;
  % xrub = (xrub - mean(displ(:,1))) * pixelsize;
  % yrub = (yrub - mean(displ(:,2))) * pixelsize;
%end %(if loop)
%end %(cellbound)   
%    


area_cell = polyarea(xrub,yrub);

cell1 = find(xv >= min(xvorig) & xv <= max(xvorig) &...
             yv >= min(yvorig) & yv <= max(yvorig));
 cell  = inpolygon(xv(cell1),yv(cell1),xrub,yrub);
cell  = cell1(find(cell > 0));       % POINTS IN THE INTERIOR OF THE CELL

%====================================================================================
%  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);
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);

%====================================================================================
%  A AND B

if k_abs(1,1) == 0,
   k_abs(1,1) = 1;
end;

Ax = a1 * ( 2 * pi ./ k_abs) + b1 * (2 * pi ./ k_abs) .* (sin(alpha)).^2;
Ay = a1 * ( 2 * pi ./ k_abs) + b1 * (2 * pi ./ k_abs) .* (cos(alpha)).^2;
B  = c1 * (-2 * pi ./ k_abs) .* sin(alpha) .* cos(alpha);

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

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

Dx = fft2(ux);
Dy = fft2(uy);

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

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

%====================================================================================
%  ITERATIONS

step = 1;
are_the_same = 0;

while step <= iteration & are_the_same == 0,
   
%====================================================================================
%  FORM NEW (MIXED) TRACTION MATRIX

   tx2 = zeros(N,N);
   ty2 = zeros(N,N);

   tx2(cell) = tx(cell);   
   ty2(cell) = ty(cell);   
   
%====================================================================================
%  CALCULATE THE INDUCED DISPLACEMENTS

Tx = fft2(tx2);       
Ty = fft2(ty2);      

Dx = Ax .* Tx  +  B  .* Ty;  
Dy = B  .* Tx  +  Ay .* Ty;  

dx2 = real(ifft2(Dx));
dy2 = real(ifft2(Dy));

%====================================================================================
%  FORM NEW (MIXED) DISPLACEMENT MATRIX

   dx3 = dx2;                   
   dy3 = dy2;   

   dx3(cell) = ux(cell);  
   dy3(cell) = uy(cell);

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

Dx = fft2(dx3);
Dy = fft2(dy3);

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

tx3 = real(ifft2(Tx));
ty3 = real(ifft2(Ty));

%====================================================================================
%  FORM NEW (MIXED) TRACTION MATRIX
   
   tx = 0.5 * tx3 + 0.5 * tx2;
   ty = 0.5 * ty3 + 0.5 * ty2;
   
if spacing <= 4,
   tx = 0.2 * tx3 + 0.8 * tx2;
   ty = 0.2 * ty3 + 0.8 * ty2;
end;

maxt(step) = max(max(tx3.^2 + ty3.^2));

%====================================================================================
% CHECK IF THE TRACTIONS HAVE CONVERGED

if step > 3,
   if abs(maxt(step) - maxt(step-1)) <= the_same * maxt(step),
      are_the_same = 1; 
   end;
end;
step = step + 1;
end; %(iteration)
if are_the_same ~= 1, disp(' Convergence was not reached!'); end;

%====================================================================================
%  SUBTRACT mean(t(cell)), BACK TO ORIGINAL ux 
  
tx = real(tx3); 
ty = real(ty3);
tx(cell) = tx(cell) - mean(tx(cell));
ty(cell) = ty(cell) - mean(ty(cell));

   
tx4      = tx(cell);
ty4      = ty(cell);

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

% % % Save tx(cell) and ty(cell) into a file
%          tracfilename = strcat('/Users/chardin/Desktop/research/LAM/', datedir ,'/Tractions/',cellfolder,'/traction_',num2str(time),'.dat');
% %         
%          tracval = [sqrt(tx(cell).^2+ty(cell).^2)*1e12];
% 
%           save(tracfilename,'tracval','-ASCII');
        
%====================================================================================
% %  RMS TRACTION, TOTAL FORCE, THETA, A_CUTS, TRACE, U

 
 
%    tx5 = zeros(1,size(cell,1));
%    ty5 = zeros(1,size(cell,1));

   count  = 1;
        for loop = 1: size(cell,1),   
          %if or(abs(tx4(loop)*1e12) >=RMSTmin,abs(ty4(loop)*1e12) >=RMSTmin),
           tx5(count)        = tx4(loop);
           ty5(count)        = ty4(loop);
           count             = count + 1;
           data6(loop,:) =  {tx4(loop), ty4(loop)};
           
       %end;
        end;
 
   
rmst_iterative = sqrt(mean2(tx(cell).^2 + ty(cell).^2));   
rmst_filter = sqrt(mean2(tx5.^2 + ty5.^2));   

% clear tx4,ty4,tx5,ty5


tot_force   = sum(sum(sqrt(tx(cell).^2 + ty(cell).^2)))*spacing^2;

Dxx = sum(sum([tx] .* [x])) * spacing^2 * 1e-6; 
Dxy = sum(sum([ty] .* [x])) * spacing^2 * 1e-6;
Dyx = sum(sum([tx] .* [y])) * spacing^2 * 1e-6;
Dyy = sum(sum([ty] .* [y])) * spacing^2 * 1e-6;

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

if abs(lambda(1,1)) > abs(lambda(2,2)), theta = -theta0;
else theta = pi/4 - theta0;
end;

a_tractions = tan(theta);
a_cuts = -1 / tan(theta);
   
Trace_moment = Dxx + Dyy;

Uecm = 0.5 * sum(sum([tx(cell) ty(cell)] .* [ux(cell) uy(cell)])) * spacing^2 * 1e-6;


%======================================================================================
%  CALCULATION OF PRESTRESS

max_height = 6;      % height of the cell at the nucleus
min_height = 0.5;    % height of the cell at the edges of the cell

%======================================================================================
%  POINTS OF CUTS ON THE LINE THROUGH (0,0)

number = 0;
clear point_x point_y

for i = floor(-sqrt(2)*N/2) : ceil(sqrt(2)*N/2),
   
   p = [i * spacing * cos(theta), i * spacing * sin(theta)];
   axbo = max([max(xrub),max(yrub),-min(xrub),-min(yrub)]);
   if inpolygon(p(1),p(2),[-axbo axbo axbo -axbo],...
         [-axbo -axbo axbo axbo]) > 0, 
   	number = number + 1;   
   	point_x(number) = p(1);
      point_y(number) = p(2);
   end; %(if)
   
end; %(for)

%======================================================================================
%  HEIGHTS IN POINTS

height = zeros(number,1);
middle = find(point_x == 0 & point_y == 0);
height(middle) = max_height;
for i = 1 : middle - 1,
   height(i) = min_height + (i-1) * (max_height - min_height) / (middle - 1);
end;
for i = number : (-1) : middle + 1,
   height(i) = min_height + (number-i) * (max_height - min_height) / (number - middle);
end;

%======================================================================================
%  PRESTRESS

clear sumforce diameter area_cross prestress

for i = 1 : number,
   
%===================================================================================
   %  B_CUTS, SUMFORCE
   
   b_cuts      = point_y(i) - a_cuts * point_x(i);
   if a_cuts > 0, znak = '>'; else znak = '<'; end; 
   eval(['halfplane = find(y ' znak ' a_cuts * x + b_cuts);']);
   cutleft     = intersect(halfplane,cell);
   if ~isempty(cutleft) & ~isempty(setdiff(cell,cutleft)),
   vec         = [sqrt(1/(1+a_tractions^2)) sqrt(1/(1+a_tractions^2))*a_tractions];
   clear vec_line; 
   for j = 1 : length(cutleft), 
      vec_line(j,:) = vec; 
   end;

   sumforce(i) = sum(sum([tx(cutleft)  ty(cutleft)]  .* vec_line)) * spacing^2 * 1e12;
   	%===================================================================================
   %  DIAMETER
   
   if abs(a_cuts) < 1,
   	xdots    = min(xrub) : 0.1 : max(xrub);
      ydots    = a_cuts .* xdots + b_cuts;
   else
   	ydots    = min(yrub) : 0.1 : max(yrub);
      xdots    = (ydots - b_cuts) ./ a_cuts;
   end;
   indots      = inpolygon(xdots,ydots,xrub,yrub);
   xindots     = xdots(find(indots > 0));
   yindots     = ydots(find(indots > 0));
   jump        = 0;
   for j = 1 : length(xindots)-1,
      if abs((abs(a_cuts) < 1) * (xindots(j) - xindots(j+1)) + ...
            (abs(a_cuts) >= 1) * (yindots(j) - yindots(j+1))) > 0.11,
         jump = [jump, j];
      end;
   end;
   if length(jump) > 1,
      if jump(2) == 1, jump(2) = []; end;
   end;
   jump = [jump, length(xindots)+1];
   diameter(i)   = 0;
   for j = 1 : length(jump)-1,
      diameter(i) = diameter(i) + sqrt((xindots(jump(j)+1) - xindots(jump(j+1)-1))^2 +...
         (yindots(jump(j)+1) - yindots(jump(j+1)-1))^2);
   end;
   
	%===================================================================================
   %  AREA_CROSS
   
   radius = (diameter(i)^2 / 4 + height(i)^2) / (2 * height(i));
   alpha  = 2 * atan2(diameter(i)/2, radius - height(i));
   area_cross(i) = 0.5 * radius^2 * (alpha - sin(alpha));
   
   %===================================================================================
   %  PRESTRESS
   
   prestress(i) = sumforce(i) / area_cross(i);
   
   
   else
   prestress(i) = NaN;
   end; %(if ~isempty(cutleft))
end; %(for i)

   %disp(' ');
   %disp(strcat('The max prestress is:',num2str(max(prestress),'%10.5f'),' (Pa)'));
   %disp(' ');
xx =x((N/4+1):(N/4*3),(N/4+1):(N/4*3)) ; %cropping data for traction.dat
yy =y((N/4+1):(N/4*3),(N/4+1):(N/4*3)) ; 
txx =tx((N/4+1):(N/4*3),(N/4+1):(N/4*3)) ; 
tyy =ty((N/4+1):(N/4*3),(N/4+1):(N/4*3)) ; 
%  FIGURES
   figure(3);
   rect = [100, 100, 800, 560]; 
%      load('MyColormaps','mycmap')
    %  set(gcf,'Position',rect,'Renderer','zbuffer'); clf; hold on;
      set(gcf,'color',[1 1 1],'Position',rect,'Renderer','zbuffer'); clf; hold on;
   %  set(figure(3),'Colormap',mycmap);
    %set(gca,'color',[0 0 0]);
  

   %CAXIS([0 100])
   set(gca,'clim',[0 fig_max])
   surf(xx,yy,sqrt(txx.^2+tyy.^2)*1e12,'EdgeColor','none');
   whitebg('black');
   
%     view(2); colormap jet; shading interp; 
        view(2); %shading interp; 
     cbh = colorbar;
%    cbh = colorbar('YTickLabel',...
%    {'','','','',...
%     '','','',''});
   set(gca,'FontSize',18,'FontWeight','bold'); box on;
    set(cbh,'FontSize',18,'FontWeight','bold');
   m3 = max2(sqrt(tx.^2+ty.^2)*1e12);
     h = plot3(xrub,yrub,m3*ones(size(xrub)),'w-');
    set(h,'LineWidth',2.0);
    %h = quiver3(x,y,m3*ones(size(x)),tx,ty,zeros(size(x)),1);
    set(h,'Color','w','LineWidth',1.0);
      axis image; 
      axis tight; 
      axis ij;
     
      
%      axis off;
%     axis([-115  115  -115  115]);
    axis([min(xvorig) max(xvorig) min(yvorig) max(yvorig)]);
%   set (gca, 'XTickLabel',{'','','','',''},'YTickLabel',{'','','','',''});
   xlabel('x (\mum)','FontSize',18,'FontWeight','bold'); ylabel('y (\mum)','FontSize',18,'FontWeight','bold');
%   set(gca,'FontSize',18,'FontWeight','bold','ActivePositionProperty','outerposition'); box on;
%   Tractions','FontSize',18,'FontWeight','bold');
%    set(gca,'color',[0 0 0]);
   title(tractitle,'FontSize',18,'FontWeight','bold');
 %   plot3(Xs,Ys,m3*ones(size(Xs)),'r-');
%    plot3(xrub1,yrub1,m3*ones(size(xrub1)),'m+');
% % 
%    plot3(Xs1,Ys1,m3*ones(size(Xs1)),'g-');
%    plot3(xrub2,yrub2,m3*ones(size(xrub2)),'g+');
%   
%  Save the constrained traction figure to a frame 



figure(9); 
   rect = [100, 100, 800, 700]; 
   set(gcf,'Position',rect,'Renderer','zbuffer'); clf; hold on;
   set(gca,'climmode','auto')
   surf(xx,-yy,sqrt(txx.^2+tyy.^2)*1e12,'EdgeColor','none'); 
   view(2); colormap jet; %shading interp; 
   cbh = colorbar; set(cbh,'FontSize',18,'FontWeight','bold');
   set(gca,'FontSize',18,'FontWeight','bold'); box on;
    set(cbh,'FontSize',18,'FontWeight','bold');
   m3 = max2(sqrt(tx.^2+ty.^2)*1e12);
     %h = plot3(xrub,yrub,m3*ones(size(xrub)),'w-');
   set(h,'Color','w','LineWidth',1);
      m1 = max2(sqrt(ux.^2+uy.^2));
   %h = quiver3(x,y,10*m1*ones(size(x)),ux,uy,zeros(size(x)),1);
   %set(h,'Color','w','LineWidth',1);
   %axis equal; axis tight; axis ij; 
%  axis([-30 30 -30 30]);    
   %axis([min2(x) max2(x) min2(y) max2(y)])
   range = 450;
   axis([min(min(xx)) max(max(xx)) min(min(yy)) max(max(yy)) ]);
   %axis([-range range -range range])
%   xlabel('x (\mum)','FontSize',18,'FontWeight','bold'); ylabel('y (\mum)','FontSize',18,'FontWeight','bold');
    h = plot3(xrub,-yrub,m3*ones(size(xrub)),'w-');
    set(h,'LineWidth',2.0);
   set(gca,'FontSize',18,'FontWeight','bold'); box on;
   trt_title = strcat('Traction(Pa)');
   title(trt_title,'FontSize',18,'FontWeight','bold'); 
  %aviobj1 = addframe(aviobj1,figure(3));  

  
%%%%%% for FEM
nXPts = size(txx,1);
  BW = imread([coreDir '\croppeddata\' cellfolder '/' 'domain.tif']);
    % reshape the domain map to the traction grid size
    BW01 = reshape(imresize(BW,[nXPts nXPts])*1,nXPts*nXPts,1);
    eval(['mkdir ' coreDir '\calc\' cellfolder '/tractions/']);
    eval(['mkdir ' coreDir '\calc\' cellfolder '/traction/']);
    save([coreDir '\calc\' cellfolder '/tractions/' 'domain.dat'],'BW01','-ASCII');
    save([coreDir '\calc\' cellfolder '/traction/' 'domain.dat'],'BW01','-ASCII');
%%%%%%%  
%==================================================================
% calc rotation moment 
%==================================================================
BW02 =  imresize(BW,[nXPts nXPts]) ; 
txx(find(BW02==0)) =0; tyy(find(BW02==0)) =0;txx = txx*1.E12; tyy = tyy*1.E12;
txx = txx - mean2(txx); tyy = tyy - mean2(tyy); txx = txx*1.E-12; tyy = tyy*1.E-12;
[cx, cy] = ait_centroid(BW02) ;  %center of mass
nzero = find(BW02~=0)   ; Nnz = size(nzero,1) ;
%cxx = (cx-1)/(nXPts-1) * (max2(xx) - min2(xx)) + min2(xx) ;  % convert center of mass coordi (pixel) to real coordi
%cyy = (cy-1)/(nXPts-1) * (max2(yy) - min2(yy)) + min2(yy) ;
cxx = (cx-1)/(nXPts-1) * (max2(xx) - min2(xx)) + min2(xx) ;  % convert center of mass coordi (pixel) to real coordi
cyy = (cy-1)/(nXPts-1) * (max2(yy) - min2(yy)) + min2(yy) ;
xx_cell = xx(nzero) - cxx ; yy_cell = yy(nzero) - cyy ; 
txx_cell = txx(nzero) ; tyy_cell = tyy(nzero) ; 
coordi = [ xx_cell  yy_cell yy_cell*0]; 
trct = [ txx_cell  tyy_cell tyy_cell*0] ; 
momt_r = cross(coordi, trct) ; tot_momt = sum( momt_r(:,3)); 
[ trct_ang  trct_rho ] = cart2pol(coordi(:,1), coordi(:,2)) ; %trct_rho_sum = sum(ones(Nnz,1)./trct_rho);
[trctx     trcty] = pol2cart( trct_ang+pi/2 , ones(Nnz,1)./trct_rho*tot_momt/Nnz) ; 
txx(nzero) = txx(nzero) - trctx/1. ; tyy(nzero) = tyy(nzero) - trcty/1. ;  
txx = txx*1.E12; tyy = tyy*1.E12;
txx = txx - mean2(txx); tyy = tyy - mean2(tyy); txx = txx*1.E-12; tyy = tyy*1.E-12;
mean2(txx)
mean2(tyy)
%==================================================================
  
tractions = [reshape(x,size(x,1)*size(x,2),1)/pixelsize reshape(y,size(x,1)*size(x,2),1)/pixelsize reshape(tx,size(x,1)*size(x,2),1)*1 reshape(ty,size(x,1)*size(x,2),1)*1];
tractions_crop = [reshape(xx,size(xx,1)*size(xx,2),1)/pixelsize reshape(yy,size(xx,1)*size(xx,2),1)/pixelsize reshape(txx,size(xx,1)*size(xx,2),1)*1.E12 reshape(tyy,size(xx,1)*size(xx,2),1)*1.E12];
mkdir([coreDir '/calc/' cellfolder '/tractions']) 
mkdir([coreDir '/calc/' cellfolder '/traction']) 
mkdir([coreDir '/calc/' cellfolder '/FEM']) 

save([coreDir '\calc\' cellfolder '/tractions/' 'ct_traction.dat'],'tractions_crop','-ASCII'); 
if monolayer == 0 
save([coreDir '\calc\' cellfolder '/tractions/' 'traction.dat'],'tractions_crop','-ASCII'); 
save([coreDir '\calc\' cellfolder '/traction/' 'traction.dat'],'tractions_crop','-ASCII'); 
end 

meantx = mean(mean(txx)) 
meanty = mean(mean(tyy)) 
Matrix_Number = N/2 

%subtract cell mean traction in the boundary for unconstrained tx ty
%========================================================
cell1 = null(1);
cell = null(1);
[N_uc, NN_uc] = size(x_uc);
x_uc_ = reshape(x_uc,N_uc^2,1);
y_uc_ = reshape(y_uc,N_uc^2,1);
tx_uc_ = reshape(tx_uc,N_uc^2,1);
ty_uc_ = reshape(ty_uc,N_uc^2,1);
tx_uc_1 = zeros(N_uc^2,1);
ty_uc_1 = zeros(N_uc^2,1);

 cell1 = find(x_uc_ >= min(x_uc_) & x_uc_ <= max(x_uc_) &...
              y_uc_ >= min(y_uc_) & y_uc_ <= max(y_uc_));
 cell2  = inpolygon(x_uc_,y_uc_,xrub,yrub);
cell  = find(cell2 > 0);       % POINTS IN THE INTERIOR OF THE CELL
tx_uc_(cell) = tx_uc_(cell) - mean(tx_uc_(cell));
ty_uc_(cell) = ty_uc_(cell) - mean(ty_uc_(cell));
tx_uc_1(cell) = tx_uc_(cell);
ty_uc_1(cell) = ty_uc_(cell);

% x_uc = reshape(x_uc_,N_uc, N_uc);
% y_uc = reshape(y_uc_,N_uc, N_uc);
 tx_uc = reshape(tx_uc_1,N_uc, N_uc);
 ty_uc = reshape(ty_uc_1,N_uc, N_uc);
%========================================================

   uc_tractions = [reshape(x_uc,size(x_uc,1)*size(x_uc,2),1)/pixelsize reshape(y_uc,size(x_uc,1)*size(x_uc,2),1)/pixelsize reshape(tx_uc,size(x_uc,1)*size(x_uc,2),1)*1.E12 reshape(ty_uc,size(x_uc,1)*size(x_uc,2),1)*1.E12];
   %uc_tractions_crop = [reshape(xx,size(xx,1)*size(xx,2),1)/pixelsize reshape(yy,size(xx,1)*size(xx,2),1)/pixelsize reshape(txx,size(xx,1)*size(xx,2),1)*1.E12 reshape(tyy,size(xx,1)*size(xx,2),1)*1.E12];
   %mkdir([coreDir 'Tractions']);   
   save([coreDir '\calc\' cellfolder '/tractions/' 'uc_traction.dat'],'uc_tractions','-ASCII'); 
   if monolayer == 1
   save([coreDir '\calc\' cellfolder '/tractions/' 'traction.dat'],'uc_tractions','-ASCII'); 
   save([coreDir '\calc\' cellfolder '/traction/' 'traction.dat'],'uc_tractions','-ASCII'); 
   end
   %save([coreDir 'Tractions/traction.dat'],'uc_tractions','-ASCII','-double'); 
   avg_uc_tr  = mean(sqrt(uc_tractions(:,3).^2+uc_tractions(:,4).^2));





 %saveloc = strcat('/Users/chardin/Desktop/research/LAM/', datedir ,'/Tractions/');
  %cd(saveloc);
  % saveas(gcf, strcat(saveloc,'tracmap_',num2str(time),'.tif'));
% default = strcat('cd /Users/chardin/Desktop/Research/Cytoskeltetal_Aspects_of_VILI/PASM+PAH_serum/', datedir ,'/croppeddata/',cellfolder);
%aviobj1 = addframe(aviobj1,figure(3));  
% 
   figure(1);
    %h = plot3(xrub,yrub,m1*ones(size(xrub)),'w-');
    %set(h,'LineWidth',2.0);
   h = plot3(xrub,yrub,m3*ones(size(xrub)),'w-');
   set(h,'LineWidth',2.0);
 %  aviobj3 = addframe(aviobj3,figure(1));  
%            
  figure(2);
    h = plot3(xrub,yrub,m2*ones(size(xrub)),'w-'); set(h,'LineWidth',2.0);
 %  aviobj2 = addframe(aviobj2,figure(2));
      
   figure(4);
   left = 0.5;
   h = text(left+mov,top+difv,'Constrained FFTC'); set(h,'FontSize',12,'FontWeight','bold');
   text(left,top-difv,'RMS traction (Pa):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-2*difv,num2str(rmst_iterative*1e12,'%10.5f'),'FontSize',12,'FontWeight','bold');
   text(left,top-3*difv,'Orientation of principle tractions (^o):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-4*difv,num2str(theta0*180/pi,'%10.5f'),'FontSize',12,'FontWeight','bold');
   text(left,top-5*difv,'Net sumforce moment (pJ):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-6*difv,num2str(-Trace_moment*1e12,'%10.5f'),'FontSize',12,'FontWeight','bold');
   text(left,top-7*difv,'Total strain energy (pJ):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-8*difv,num2str(Uecm*1e12,'%10.5f'),'FontSize',12,'FontWeight','bold');
   text(left,top-9*difv,'Max cumulative force (nN):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-10*difv,num2str(max(sumforce)*1e-3,'%10.5f'),'FontSize',12,'FontWeight','bold');
%    text(left,top-11*difv,'Prestress (Pa):','FontSize',18,'FontWeight','bold'); 
%    text(left+mov,top-12*difv,num2str(stnanmean(prestress),'%10.5f'),'FontSize',18,'FontWeight','bold');
   text(left,top-13*difv,'Projected area of the cell (\mum^2):','FontSize',12,'FontWeight','bold'); 
   text(left+mov,top-14*difv,num2str(area_cell,'%10.5f'),'FontSize',12,'FontWeight','bold');
   %disp(num2str([rmst_iterative*1e12;theta0*180/pi;-Trace_moment*1e12;...
   %   Uecm*1e12;max(sumforce)*1e-3;stnanmean(prestress);area_cell],'%10.3f'));

   data1(kloop-2,:) =  {num2str(0/10), num2str(rmst_iterative*1e12,'%10.5f'),num2str(rmst_filter*1e12,'%10.5f'), ...
                        num2str(theta*180/pi,'%10.5f'), ...
            num2str(-Trace_moment*1e12,'%10.5f'),num2str(-Dxx*1e12,'%10.5f'),num2str(-Dyy*1e12,'%10.5f'), ...
        num2str(-lambda(1,1)*1e12,'%10.5f'),num2str(-lambda(2,2)*1e12,'%10.5f'), ...
        num2str(Uecm*1e12,'%10.5f'),num2str(max(sumforce)*1e-3,'%10.5f'), ...
            num2str(tot_force*1e9,'%10.5f'),...
            num2str(area_cell,'%10.5f')};  
%     
% netmoment = -lambda(1,1)*1e12;
% netmoment_time(kloop-2,2) = -Trace_moment*1e12 ; 
% netmoment_time(kloop-2,1) = time/10;
 -Trace_moment*1e12
   %====================================================================================
%  Save figures parameters  
   
%    m = findstr(savefilename,'disp');
%    n = findstr(savefilename,'.');
%    savedate = savefilename(1:m-1);
%    savedate(findstr(savedate,'-')) = '_';
%    savenumber = savefilename(m+4:n-1);
%   disptitle = strcat('Displacement',savenumber);
%   tractitle = strcat('Constrained FTTC: Tractions',savenumber);
   
   
%====================================================================================
%  SAVE FIGURES
%    m = findstr(savefilename,'disp');
%    n = findstr(savefilename,'.');
%    savedate = savefilename(1:m-1);
%    savedate(findstr(savedate,'-')) = '_';
%    savenumber = savefilename(m+4:n-1);
%    saveas(1,strcat(savedirectory,savedate,'_Disp',savenumber,'.fig'));
%    saveas(2,strcat(savedirectory,savedate,'_Unconstraint_Traction',savenumber,'.fig'));   
%    saveas(3,strcat(savedirectory,savedate,'_Constraint_Traction',savenumber,'.fig'));
%    saveas(4,strcat(savedirectory,savedate,'_Result',savenumber,'.fig'));

% %====================================================================================
%  Save raw tractions to file

%filename=strcat('/Users/chardin/Desktop/research/LAM/',datedir ,'/Tractions/',cellfolder,'/rawtraction',num2str(time),'.dat');
header='';
colnames2={'x';'y';'tx';'ty'};
sheetname=num2str('1');
xvalue=reshape(x,size(x,1)*size(x,1),1);
yvalue=reshape(y,size(y,1)*size(y,1),1);
txvalue=reshape(tx,size(x,1)*size(x,1),1);
tyvalue=reshape(ty,size(y,1)*size(y,1),1);
xvalue = xvalue ./pixelsize;
yvalue = yvalue ./pixelsize;
txvalue = txvalue ./ 1e-12;
tyvalue = tyvalue ./ 1e-12;
tracdata=[xvalue yvalue txvalue tyvalue];
%yvalue = yvalue ./ pixelsize;
%eval('cd /Users/chardin/Desktop/Research/Cytoskeltetal_Aspects_of_VILI/2_11_11_PASM+PAH_serum/021111/Tractions/',cellfolder);

save([coreDir '\calc\' cellfolder '/tractions/' 'rawtraction.dat'],'tractions','-ASCII'); 
%save(filename,'tracdata','-ASCII');

% %====================================================================================


%end; %(kloop)

 header   = '';       
 
     colnames = {'Filename';'RMS traction (Pa)';'RMS traction with threshold (Pa)';'Orientation of principle tractions'; ...
    'Net contractile moment (pJ)';'Mxx (pJ)';'Myy (pJ)';'Lambda(1,1) (pJ)';'Lambda(2,2) (pJ)';'Total strain energy (pJ)'; 'Max cumulative force (nN)'; 'Total force (nN)';...
    'Projected area of the cell (\mum^2)'}; 
     colnames1 = {'Filename';'RMS traction (Pa)';'Orientation of principle tractions'; ...
     'Net contractive moment (pJ)'; 'Total strain energy (pJ)'};

     xlsfilename = 'ConstrainedFTTC.xls';
     xlsfilename1 = 'UnconstrainedFTTC.xls';
     xlsfilename2 = 'tractions.xls';
     
%*     cd('C:\C2C12\traction');
      
  %figure;
  %plot(netmoment_time(:,1),netmoment_time(:,2));
  %RMSTrace = [rmst_onestep*1e12, rmst_iterative*1e12, -Trace_moment*1e12];
  %avg_uc_tr = unconstrained in cell boundary
   RMSTrace = [avg_uc_tr, rmst_iterative*1e12, -Trace_moment*1e12];
      cd([coreDir '\calc\' cellfolder '/tractions/']);
     xlswrite(xlsfilename,data1,cellfolder);     
     xlswrite(xlsfilename1,data2,cellfolder);     
     xlswrite(xlsfilename2,RMSTrace,cellfolder);
 if monolayer == 1    
 Matrix_Number = (size(x_uc,1)) ; 
 elseif monolayer == 0   
 Matrix_Number = (size(xx,1))  ;
 end 
  E_cell = 1.e9 ;
  thick_cell = 5.e-6 ;
fid = fopen([coreDir '/calc/' cellfolder '/Model.in'], 'wt');
                A = [num2str(Matrix_Number) ' ' num2str(Matrix_Number) '\n'] ; fprintf(fid, A );
                B = [num2str(pixelsize) 'e-6' ' \n'] ; fprintf(fid, B );
                C = [num2str(E_cell) 'd0  ' num2str(pois) 'd0' ' \n'] ;  fprintf(fid, C );
                D = [num2str(thick_cell) ' \n'] ;  fprintf(fid, D ); fclose(fid);            
%     xlswrite(data6,cellfolder);
      
% aviobj1 = close(aviobj1);
% aviobj2 = close(aviobj2);
% aviobj3 = close(aviobj3);
clear T* t* u* U* v* x* y*; 

close all;

%RMSTout 
%netmoment_time(:,2)

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