clear
close all
tic

Position=[200 200 900 800];
TitleSz=24; LblSz=20; AxSz=16;
NL=sprintf('\n');

raw=load('VarDat.dat');
Ey=raw(1);
G=raw(2); 
StudyCt=raw(3);

V0=0; 

X0=load('X0.dat');
Xdat=load('X.dat');
Ydat=load('Y.dat');
Zdat=load('Z.dat');
Vx=load('Vx.dat');
Vy=load('Vy.dat'); 
% Vz=load('Vz.dat'); 
dVdx=load('dVdX.dat');
dVdy=load('dVdY.dat');
dVdz=load('dVdZ.dat');
P=load('P.dat');
rho=load('rho.dat');
Fixed=load('Fixed.dat');
%Tvm=load('VonMises.dat');
E11=load('E11.dat');
E22=load('E22.dat');
dx=load('dX.dat');
T11=load('T11.dat');
T21=load('T21.dat');
T31=load('T31.dat');
T12=load('T12.dat');
T22=load('T22.dat');
T32=load('T32.dat');
T13=load('T13.dat');
T23=load('T23.dat');
T33=load('T13.dat');
IsSolid=load('IsSolid.dat');
Mass=load('Mass.dat');
rho0=load('rho0.dat');
rho_init=load('rho_init.dat');
TimeFct=load('Time.dat');

Cylrng=find(Fixed==8);
R1=max(X0(Cylrng,2))-min(X0(Cylrng,2))+(mean(dx)/2);

foo=size(Xdat);
ct=foo(1); stop=ct; a=ct;
N=foo(2); clear foo
if stop>ct
    stop=ct;
end
plotCt=5; Range=(1:plotCt)*floor(ct/plotCt);
Vy(:,1)=-V0; Vy(:,N)=V0;
Poisson=(Ey/(2*G))-1;
% ctZ=2*length/dx;

Ered=1/((1-(Poisson^2))/Ey);

Xdat2=Xdat;
Ydat2=Ydat;
Zdat2=Zdat;
for ii=1:N
    Xdat2(:,ii)=Xdat2(:,ii)-(X0(ii,1));
    Ydat2(:,ii)=Ydat2(:,ii)-(X0(ii,2));
    Zdat2(:,ii)=Zdat2(:,ii)-(X0(ii,3));
end

Xdat3=Xdat2;
Ydat3=Ydat2;
Zdat3=Zdat2;
for ii=1:ct
    for jj=1:N
        Xdat3(ii,jj)=Xdat3(ii,jj)-mean(Xdat2(ii,:));
        Ydat3(ii,jj)=Ydat3(ii,jj)-mean(Ydat2(ii,:));
        Zdat3(ii,jj)=Zdat3(ii,jj)-mean(Zdat2(ii,:));
    end
end

NL=(['\n']); 
%dtplot=load('Time.dat'); 

Vr=sqrt((Vx.^2)+(Vy.^2));
Xr=sqrt((Xdat.^2)+(Ydat.^2));
Xr2=Xr;
for ii=1:ct
    Xr2=Xr(ii,:)-Xr(1,:);
end

save Data

toc
