求助！！！积分积不出

zhouyepeng123

1. clear
   
2. clc
   
3. close all
   
4. n1=1.33;n2=1.6;
   
5. s=linspace(0,1.4,700);
   
6. nn=0;
   
7. Qg=[];
   
8. for i=1:length(s)
   
9.     theta=0.1:0.1:70;
   
10.     a=theta*pi/180;
    
11.     r=asin(n1/n2*sin(a));%折射角r
    
12.     rp=(n2*cos(a)-n1*sqrt(1-(n1/n2*sin(a)).^2))./...
    
13.         (n2*cos(a)+n1*sqrt(1-(n1/n2*sin(a)).^2));%反射率p分量
    
14.     rs=(n1*cos(a)-n2*sqrt(1-(n1/n2*sin(a)).^2))./...
    
15.         (n1*cos(a)+n2*sqrt(1-(n1/n2*sin(a)).^2));%反射率s分量
    
16.     tp=2*n1*cos(a)./(n2*cos(a)+n1*sqrt(1-(n1/n2*sin(a)).^2));%透射率p分量
    
17.     ts=2*n1*cos(a)./(n1*cos(a)+n2*sqrt(1-(n1/n2*sin(a)).^2));%透射率s分量
    
18.     R=(abs(rp).^2+abs(rs).^2)/2;
    
19.     T=1-R;%给T赋值
    
20.     ss=s(i);
    
21.    
    
22.     Q1=(2*cos(a)./((1-(sin(a).^2)./ss.^2).*sin(a))).*...
    
23.         (R.*sin(2*a)-(T.^2.*(sin(2*a-2*r)+R.*sin(2*a))./(1+R.^2+2*R.*cos(2*r))));
    
24.     Qg(i)=trapz(a,Q1);
    
25. end
    
26. figure(1)
    
27. plot(a,Q1)
    
28. hold on
    
29. figure(2)
    
30. plot(s,Qg)

复制代码