N = input('N=') ;
EL=180;
AL=-0.025E-6;
AT=27E-6;
ET=10;
RL=0.3;
RT=0.03;
GLT=5;
Q11=EL/(1-RL*RT);
Q22=ET/(1-RL*RT);
Q12=ET*RL/(1-RT*RL);
Q66=GLT;
Q=[Q11 Q12 0;Q12 Q22 0;0 0 Q66]
Q11L=zeros(1,N+1);
Q22L=zeros(1,N+1);
Q12L=zeros(1,N+1);
Q66L=zeros(1,N+1);
Q16L=zeros(1,N+1);
Q26L=zeros(1,N+1);
j=3;k=3;
T2=zeros(j,k,N+1);
T1=zeros(j,k,N+1);
for i=1:N+1;
c(i)=cos((i-1)*pi/180);
s(i)=sin((i-1)*pi/180);
Q11L(i)=Q11*(c(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*s(i)^4;
Q22L(i)=Q11*(s(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*c(i)^4;
Q12L(i)=Q12*(c(i)^4+s(i)^4)+(Q11+Q22-4*Q66)*(s(i)^2)*(c(i)^2);
Q66L(i)=Q66*(c(i)^4+s(i)^4)+(Q11+Q22-2*Q12-2*Q66)*s(i)^2*c(i)^2;
Q16L(i)=(Q11-Q12-2*Q66)*c(i)^3*s(i)+(-Q22+Q12+2*Q66)*c(i)*s(i)^3;
Q26L(i)=(Q11-Q12-2*Q66)*c(i)*s(i)^3+(-Q22+Q12+2*Q66)*c(i)^3*s(i);
Q21L(i)=Q12L(i);
Q61L(i)=Q16L(i);
Q62L(i)=Q26L(i);
T1(:,:,i)=[c(i)^2 s(i)^2 c(i)*s(i);s(i)^2 c(i)^2 -s(i)*c(i);-2*s(i)*c(i) 2*s(i)*c(i) c(i)^2-s(i)^2];
end
for i=1:N+1;
c(i)=cos((i-1)*pi/180);
s(i)=sin((i-1)*pi/180);
Q11L1(i)=Q11*(c(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*s(i)^4;
Q22L1(i)=Q11*(s(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*c(i)^4;
Q12L1(i)=Q12*(c(i)^4+s(i)^4)+(Q11+Q22-4*Q66)*(s(i)^2)*(c(i)^2);
Q66L1(i)=Q66*(c(i)^4+s(i)^4)+(Q11+Q22-2*Q12-2*Q66)*s(i)^2*c(i)^2;
Q16L1(i)=(Q11-Q12-2*Q66)*c(i)^3*s(i)+(-Q22+Q12+2*Q66)*c(i)*s(i)^3;
Q26L1(i)=(Q11-Q12-2*Q66)*c(i)*s(i)^3+(-Q22+Q12+2*Q66)*c(i)^3*s(i);
Q21L1(i)=Q12L1(i);
Q61L1(i)=Q16L1(i);
Q62L1(i)=Q26L1(i);
T2(:,:,i)=[c(i)^2 s(i)^2 c(i)*s(i);s(i)^2 c(i)^2 -s(i)*c(i);-2*s(i)*c(i) 2*s(i)*c(i) c(i)^2-s(i)^2];
end
j=3;k=3;
n=N+1;
A1=zeros(j,k,n);
A2=zeros(j,k,n);
for i=1:n;
A1(:,:,i)=[ Q11L(i) Q12L(i) Q16L(i);Q21L(i) Q22L(i) Q26L(i);Q61L(i) Q62L(i) Q66L(i)];
A2(:,:,i)=[ Q11L1(i) Q12L1(i) Q16L1(i);Q21L1(i) Q22L1(i) Q26L1(i);Q61L1(i) Q62L1(i) Q66L1(i)];
T=T1+T2;
A=A1+A2;
a(:,:,i)=inv(A(:,:,i));
u(i)=-a(1,2,i)/a(1,1,i);
end
j=3;k=1;ALF=zeros(j,k,n);W=zeros(j,k,1);
for i=1:n;
ALF(1,1,i)=AL*c(i)^2+AT*s(i)^2;
ALF(2,1,i)=AL*s(i)^2+AT*c(i)^2;
ALF(3,1,i)=2*s(i)*c(i)*(AL-AT);
alf(:,:,i)=-a(:,:,i)*T(:,:,i)*Q*ALF(:,:,i);
a1(i)=alf(1,1,i);
a2(i)=alf(2,1,i);
a3(i)=alf(3,1,i);
end
t=0:90;
subplot(2,2,1)
plot(t,u)
title('泊松比','FontWeight','bold')
xlabel('铺层角度')
ylabel('泊松比')
subplot(2,2,2)
plot(t,a1,'r');
title('x方向热膨胀系数','FontWeight','bold')
xlabel('铺层角度')
ylabel('膨胀系数')
subplot(2,2,3)
plot(t,a2,'g');
title('y方向热膨胀系数','FontWeight','bold')
xlabel('铺层角度')
ylabel('膨胀系数')
subplot(2,2,4)
plot(t,a3,'b');
title('剪切方向热膨胀系数','FontWeight','bold')
xlabel('铺层角度')
ylabel('膨胀系数')
