################################################
# filename: bsdf_simple.lsf
#
# Description: This file does a simple
# analyis of light scattering from a rough
# surface at normal incidence.
#  
# Copyright 2009 Lumerical Solutions Inc
################################################

# run the analyis if necessary
if(!haveresult("BSDF")) {
  run;
  runanalysis;
  save;
}


# collect the far field results
BSDF_up = getresult("BSDF","BSDF_up");

f = BSDF_up.f;
ux = BSDF_up.ux;
uy = BSDF_up.uy;
ES = BSDF_up.ES;
EP = BSDF_up.EP;
E2 = (BSDF_up.E2) *(1+transmission("BSDF::up"));
n = BSDF_up.n;
m = BSDF_up.m;

# remove the specular
spec_R = E2(find(n,0),find(m,0));
E2(find(n,0),find(m,0)) = 0.25*(E2(find(n,0),find(m,1)) +
                                E2(find(n,1),find(m,0)) +
                                E2(find(n,0),find(m,-1)) +
                                E2(find(n,-1),find(m,0)));

spec_R = spec_R - E2(find(n,0),find(m,0));
total_power_R = spec_R + sum(E2);

# image the far field
image(ux,uy,E2,"","","log10(order strengths)","polar,logplot");

# calculate scattering as a function of theta,
# average over many aziumuthal angles
E2 = pinch(E2);
theta = linspace(0,89,100);
phi = linspace(0,360,50);

Ux = meshgridx(ux,uy);
Uy = meshgridy(ux,uy);
cos_theta = sqrt(1-Ux^2-Uy^2);

result = 0;
for(i=1:(length(phi)-1)) {
    result = result + farfield3dintegrate(E2*cos_theta,ux,uy,5,theta,phi(i))/
                     (farfield3dintegrate(0*E2+1,ux,uy,5,theta,phi(i))+1e-20);
}
result = result/(length(phi)-1);

# compare with the theoretical result
k0 = 2*pi*f/c;
select("rough_surf");
cor_len=get("corr length x"); 
krho = k0 * sin(theta*pi/180);

BRDF_theory = exp(-(1/8)*krho^2*cor_len^2);
BRDF_sim = result/sum(result*sin(theta*pi/180))*sum(BRDF_theory*sin(theta*pi/180));

plot(theta,real(BRDF_sim),BRDF_theory,"theta (degrees)","BRDF (a.u.)");
legend("simulated","theory");