Below is the code. I got error when I tried to integrate Q_total. Below is the error i got
ERROR:
Index exceeds the number of array elements (202).
Error in main (line 95)
Q_total=integral(qtot(w(a)),0,w);
CODE:
for a=1:length(w)
betaevanmax=10/d; %upper limit of the parallel wavevector for evanescent wave integration (see ACS NANO paper method section)
betaevanmin=w(a)/c; %lower limit
%Finding spectral RHF for prop and evan fields
qevan_p(a)=(Theta(w(a), T1,0)-Theta(w(a), T2, V))*1/pi^2*integral(@(beta)sevan_p(w(a),beta, d),betaevanmin,betaevanmax,'Reltol',1E-6,'AbsTol',1E-10); %default abstol = 1E-10 reltol 1E-6
qevan_s(a)=(Theta(w(a), T1,0)-Theta(w(a), T2, V))*1/pi^2*integral(@(beta)sevan_s(w(a),beta, d),betaevanmin,betaevanmax,'Reltol',1E-6,'AbsTol',1E-10);
qprop_p(a)=(Theta(w(a), T1,0)-Theta(w(a), T2, V)).*1/pi^2.*integral(@(beta)(sprop_p(w(a),beta,d)),0,w(a)/c);
qprop_s(a)=(Theta(w(a), T1,0)-Theta(w(a), T2, V)).*1/pi^2.*integral(@(beta)(sprop_s(w(a),beta,d)),0,w(a)/c);
a/length(w)% progress indicator
end
qtot=qevan_p+qevan_s+qprop_p+qprop_s; %total intensity spectum
qtot=qtot*1.5193E15/10000; %in w/cm2/eV
qprop=(qprop_p+qprop_s)*1.5193E15/10000; %prop waves only, in w/cm2/eV
qevan=(qevan_p+qevan_s)*1.5193E15/10000; %evan waves only, in w/cm2/eV
bb=blackb(w,T1)*1.5193E15/10000; %in w/cm2/eV
%Spectral efficiency calculation
Q_total=integral(qtot(w(a)),0,w);
