# Scriptfile: usr_far_field_filter.ldf
#
# Description:
#    This file is used with usr_far_field_filter.fsp to calculate the
#    angular distribution of radiation from an isotropic
#    ensemble of incoherent dipole emitters. The dipoles
#    are contained in a halfspace of index n = 2 and emit
#    into air.
#    
#    The far field filter is used to smooth the field 
#    truncation effects

# choose to rerun the simulation or not
rerun_simulations = 1; # set to 0 to not rerun simulations

# turn on smoothing of the far field projection
farfieldfilter(0);

load("usr_far_field_filter"); # load the initial file
if(rerun_simulations) {

    # simulation 1 has a dipole angle of 0
    switchtolayout;
    select("Source1");
    set("angle",0);
    save("usr_far_field_filter_0");
    run(1);

    # simulation 2 has a dipole angle of 90
    switchtolayout;
    select("Source1");
    set("angle",90);
    save("usr_far_field_filter_90");
    run(1);
}


# set the index of the 2 media
n1 = 2;
n2 = 1;

# collect data from both polarizations in the far field
monname = "above_power";
load("usr_far_field_filter_0");
E2 = farfield2d(monname,1,1000);
T = transmission(monname);
load("usr_far_field_filter_90");
E2 = E2 + farfield2d(monname,1,1000);
theta_projection = pi/180*farfieldangle(monname,1,1000);
theta2 = linspace(-pi/2,pi/2,1000);
T = T+transmission(monname);
T = T/2;

# interpolate the data onto a linear mesh
E2 = interp(E2,theta_projection,theta2);
# normalize the data to the near field transmission
E2 = T * E2 / integrate(E2,1,theta2);

# calculate the theoretical curve with Fresnel coefficients
# S polarization
theta1 = n2*sin(theta2)/n1;
R_s = ( (n1*cos(theta1) - n2*cos(theta2))/(n1*cos(theta1) + n2*cos(theta2)) )^2;
T_s = 1 - R_s;
dtheta1_dtheta2 = n2*cos(theta2)/ (n1*cos(theta1));
E2_s = T_s * dtheta1_dtheta2;
# normalize the data (only integrate below critical angle)
E2_s = (max(theta1)-min(theta1))/(2*pi)*E2_s / integrate(E2_s,1,theta2); 

# P polarization
theta1 = asin(n2*sin(theta2)/n1);
R_p = ( (n1*cos(theta2) - n2*cos(theta1))/(n1*cos(theta2) + n2*cos(theta1)) )^2;
T_p = 1 - R_p;
dtheta1_dtheta2 = n2*cos(theta2)/ (n1*cos(theta1));
E2_p = T_p * dtheta1_dtheta2;
# normalize the data  (only integrate below critical angle)
E2_p = (max(theta1)-min(theta1))/(2*pi)*E2_p / integrate(E2_p,1,theta2);

# normalize the data to give the same transmission as the simulation
E2_p_norm = E2_p / integrate(E2_p,1,theta2) * integrate(E2,1,theta2); 
E2_s_norm = E2_s / integrate(E2_s,1,theta2) * integrate(E2,1,theta2);

# plot the results using the unnormalized theoretical curve and the normalized curve
plot(theta2*180/pi,E2,E2_p_norm,"angle (degrees)","|E|^2");
legend("simulated","theory");

# reload the template simulation
load("usr_far_field_filter");