-
Notifications
You must be signed in to change notification settings - Fork 1
/
SphereSc.m
executable file
·100 lines (88 loc) · 2.67 KB
/
SphereSc.m
1
% SphereSc Stoke's vector of scattered light.% SphereSc(m,x,I,ang) returns the scattered Light for% a sphere, size x, refractive index relative to medium m% at angle ang.% Incident light has Stoke's Vector I.%% If m and/or x are vectors, return will 3D, with m and/or x% parameters spanning third dimension.% Written by and copyright% Dave Barnett% Optical Engineering Group% Dept. of Electrical and Electronic Engineering% Loughborough University% 20th November 1996% Corrected 4th September 1997% i missing from calculation of S34% 5th September 1997% calculation of E made more efficient.% Modified 14th December 1998% Generalised to receive vectors of m and x and return 3d matrices.% 30th December 1998% Changed to allow vector m and scalar x.% 8th January 1999% Large ranges of x now computable - NaN check added to scatcoefs% 25th January 1999% Remove zero particle size or refractive index from invalidity% criterion, to avoid infinite loop.% m = ratio of index of refraction of particle to surrounding medium% x = ka, where k = 2*pi/lambda (surrounding medium), a = radius of particle% l = 1;% ang = anglefunction [S,Int] = SphereSc(m, x, I, ang)if length(x)==1 x = x*ones(size(m));endif length(m)==1 m = m*ones(size(x));endnc = ceil(max(x)+4.05*(max(x)^(1/3))+2);n=(1:nc).';E = ((2*n+1)./(n.*(n+1)))*ones(1,length(x));[p,t] = ALegendr(ang,nc);W = warning;warning off[a,b] = ScatCoef(m,x,nc);% Check for invalid (NaN) results due to too many terms in% relatively small particles.invalid = find(any(isnan([a;b])));while ~isempty(invalid) a(:,invalid) = 0; b(:,invalid) = 0; nc2 = ceil(max(x(invalid))+4.05*(max(x(invalid))^(1/3))+2); [A,B] = ScatCoef(m(invalid),x(invalid),nc2); a(1:nc2,invalid) = A; b(1:nc2,invalid) = B; invalid = find(any(isnan([a;b]))); % remove invalidity of zero m or x % these _should_ return NaN! if length(x)>=max(invalid) invalid = invalid(x(invalid)~=0); else if x==0 invalid = []; end end if length(m)>=max(invalid) invalid = invalid(m(invalid)~=0); else if m==0 invalid = []; end endendwarning(W);a = a.*E;b = b.*E;S1 = a.'*p + b.'*t;S2 = a.'*t + b.'*p;S11 = ((S2.*conj(S2))+(S1.*conj(S1)))/2;S12 = ((S2.*conj(S2))-(S1.*conj(S1)))/2;S33 = ((S1.*conj(S2))+(S2.*conj(S1)))/2;S34 = i*((S1.*conj(S2))-(S2.*conj(S1)))/2;S = zeros(4,length(ang),length(x));S(1,:,:) = permute(I(1)*S11 + I(2)*S12,[3 2 1]);S(2,:,:) = permute(I(1)*S12 + I(2)*S11,[3 2 1]);S(3,:,:) = permute(I(3)*S33 + I(4)*S34,[3 2 1]);S(4,:,:) = permute(-I(3)*S34 + I(4)*S33,[3 2 1]);Int=1./x.^2.*(S1.*conj(S1));