detat=4.0e-7;speed=1800;pianyi=0.005;lth=1024;
%C25的強度，波速一般是3800M/S——4170M/S 
%C30的強度，波速一般是4200M/S以上
%load alldata.mat;
number=1024;
[sigmatr,sigmatr2]=dataget(number,detat,pianyi,speed);
 for k=1:number
   n=1;
   m=1;  A0121(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0122(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0123(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A0124(k,1)=sigmatr2(k,m);                                                        
   n=3;
   m=m+1;A0319(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0320(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A0321(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0322(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0323(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0324(k,1)=sigmatr2(k,m);                                                     
   m=m+1;A0325(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0326(k,1)=sigmatr2(k,m);                                                       
   n=5;
   m=m+1;A0517(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0518(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0519(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0520(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0521(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0522(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0523(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0524(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0525(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0526(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0527(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0528(k,1)=sigmatr2(k,m);                                                     
   n=7;
   m=m+1;A0719(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0720(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0721(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0722(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0723(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0724(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0725(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0726(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0727(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0728(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0729(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0730(k,1)=sigmatr2(k,m);                                                        
   n=9;
   m=m+1;A0921(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0922(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0923(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A0924(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0925(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0926(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0927(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0928(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0929(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0930(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A0931(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A0932(k,1)=sigmatr2(k,m);                                                      
   n=11; 
   m=m+1;A1123(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1124(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1125(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1126(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1127(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1128(k,1)=sigmatr2(k,m);                                                     
   m=m+1;A1129(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1130(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1131(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1132(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1133(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A1134(k,1)=sigmatr2(k,m);                                                       
   n=13;
   m=m+1;A1325(k,1)=sigmatr2(k,m);                                                     
   m=m+1;A1326(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1327(k,1)=sigmatr2(k,m);                                                     
   m=m+1;A1328(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1329(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1330(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1331(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1332(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1333(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1334(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1335(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1336(k,1)=sigmatr2(k,m);                                                     
   
   n=15; 
   m=m+1;A1527(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1528(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1529(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1530(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1531(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A1532(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1533(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1534(k,1)=sigmatr2(k,m);                                                       
   
   n=17;
   m=m+1;A1729(k,1)=sigmatr2(k,m);                                                       
   m=m+1;A1730(k,1)=sigmatr2(k,m);                                                        
   m=m+1;A1731(k,1)=sigmatr2(k,m);                                                      
   m=m+1;A1732(k,1)=sigmatr2(k,m); 
 end                                                      


k=1:lth;
t=k*detat;

x=[20*pianyi 16*pianyi 12*pianyi]';
v(k,1)=A0121; 
v(k,2)=A0319; 
v(k,3)=A0517;
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,1)=v(k,1)+v(k,2)+v(k,3);
xmatr(k,1)=11*pianyi;
tmatr(k,1)=k*detat;

x=[21*pianyi 17*pianyi 13*pianyi]';%%%
v(k,1)=A0122; v(k,2)=A0320; v(k,3)=A0518;%%%
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,2)=v(k,1)+v(k,2)+v(k,3);
xmatr(k,2)=11.5*pianyi;
tmatr(k,2)=k*detat;



x=[22*pianyi 18*pianyi 14*pianyi]';
v(k,1)=A0123; v(k,2)=A0321; v(k,3)=A0519;
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,3)=v(k,1)+v(k,2)+v(k,3);
xmatr(k,3)=12*pianyi;
tmatr(k,3)=k*detat;


x=[23*pianyi 19*pianyi 15*pianyi]';
v(k,1)=A0124; v(k,2)=A0322; v(k,3)=A0520;
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,4)=v(k,1)+v(k,2)+v(k,3);
xmatr(k,4)=12.5*pianyi;
tmatr(k,4)=k*detat;


x=[20*pianyi 16*pianyi 12*pianyi]';%%%%
v(k,1)=A0323; v(k,2)=A0521; v(k,3)=A0719;  %%%%
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,5)=v(k,1)+v(k,2)+v(k,3);%%%%
xmatr(k,5)=13*pianyi;          %%%%
tmatr(k,5)=k*detat;%%%%


x=[21*pianyi 17*pianyi 13*pianyi]';%%%%
v(k,1)=A0324; v(k,2)=A0522; v(k,3)=A0720;  %%%%
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,6)=v(k,1)+v(k,2)+v(k,3);%%%%
xmatr(k,6)=13.5*pianyi;          %%%%
tmatr(k,6)=k*detat;%%%%



x=[22*pianyi 18*pianyi 14*pianyi]';%%%%
v(k,1)=A0325; v(k,2)=A0523; v(k,3)=A0721;  %%%%
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);
end 
for k=(nnormal+1):1024
    v(k,2)=0;
end
tmin=x(3)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(3)/2)^2)/(speed*detat));
   v(nnormal,3)=v(k,3);
end 
for k=(nnormal+1):1024
    v(k,3)=0;
end
k=1:1024;
vmatr(k,7)=v(k,1)+v(k,2)+v(k,3);%%%%
xmatr(k,7)=14*pianyi;          %%%%
tmatr(k,7)=k*detat;%%%%


x=[23*pianyi 19*pianyi 15*pianyi]';%%%%
v(k,1)=A0326; v(k,2)=A0524; v(k,3)=A0722;  %%%%
tmin=x(1)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(1)/2)^2)/(speed*detat));
   v(nnormal,1)=v(k,1);
end 
for k=(nnormal+1):1024
    v(k,1)=0;
end
tmin=x(2)/speed; 
ntmin=round(tmin/detat)+1;
for k=ntmin:1024
   nnormal=round(sqrt((k*detat*speed)^2-(x(2)/2)^2)/(speed*detat));
   v(nnormal,2)=v(k,2);