你这个程序运行太慢,达15秒之多,我改了一下,你试试
c=3e8;
n1=1.33;
n2=1.6;tic
sat=0:0.0003:pi/2;
r=asin(sin(sat)*n1/n2);
A=tan(sat+r);
B=sat+r;
Rp=(sin(sat-r)./sin(sat+r)).^2;
Rs=(tan(sat-r)./tan(sat+r)).^2;
fp=1/2;fs=1/2;
R=Rp.*fp+Rs.*fs;
T=1-R;
Qs=1+R.*cos(2.*sat)-(T.^2.*(cos(2.*sat-2.*r)+R.*cos(2.*sat)))./(1+R.^2+2.*R.*cos(2.*r));
Qg=-(R.*sin(2.*sat)-(T.^2.*(sin(2.*sat-2.*r)+R.*sin(2.*sat)))./(1+R.^2+2.*R.*cos(2.*r)));
Qmag=sqrt(Qs.^2+Qg.^2);
plot(sat*180/pi,Qs,'r');
title('单根光线的梯度力和散射力值大小的图形关系')
grid on
hold on
plot(sat*180/pi,Qg,'b');hold on
plot(sat*180/pi,Qmag,'g');legend('Qs','Qg','Qmag',2)
hold off
toc
r=asin(sin(sat)*n1/n2);
A=tan(sat+r)
当sat=pi/2的时候显然tansat+r会超过pi/2,导致A会小于0