clear; closeall;
#R = 625 nm, G = 525 nm, B = 460 nm
f_red = c/625e-9;
f_green = c/525e-9;
f_blue = c/460e-9;
n = [1.5; 2.13; 1.87; 1.94; 1.75; 0.644+5.28i];
d = [0; 0.09; 0.177; 0.048; 0.029; 0]*1e-6;
n_dipoles = 1; 
z = sum(d(1:3))+d(4)/2;
orientation = matrix(n_dipoles)*0;

spectrum = 1.;
angular_res = 100;

result_unpol_red = stackdipole(n,d,f_red,z,spectrum,'rand',angular_res,-1);
result_unpol_green = stackdipole(n,d,f_green,z,spectrum,'rand',angular_res,-1);
result_unpol_blue = stackdipole(n,d,f_blue,z,spectrum,'rand',angular_res,-1);

#plot results
plotxy( result_unpol_blue.theta,result_unpol_blue.radiance,
        result_unpol_green.theta,result_unpol_green.radiance,
        result_unpol_red.theta,result_unpol_red.radiance,
        "angle (degrees)","Radiance (Watts/m2/steradian)","Radiance","linewidth=2");
legend("Blue","Green","Red");

Power_Red = 2*pi*integrate(result_unpol_red.radiance*sin(pinch(result_unpol_red.theta)*pi/180),1,pinch(result_unpol_red.theta)*pi/180);
Power_Green = 2*pi*integrate(result_unpol_green.radiance*sin(pinch(result_unpol_green.theta)*pi/180),1,pinch(result_unpol_green.theta)*pi/180);
Power_Blue = 2*pi*integrate(result_unpol_blue.radiance*sin(pinch(result_unpol_blue.theta)*pi/180),1,pinch(result_unpol_blue.theta)*pi/180);

################################################################
# Export to file for Zemax
################################################################

# Red
Ntot = 400000;
Total_radiated_power = 1; # in Watts, set to 1, scale Green and Blue appropriately
xspan = 0.25e-3;
yspan = 0.4e-3;
export_filename = "zemax_rayset_red.txt";
radiant_intensity = result_unpol_red.radiance; # note that radiance may not be correct here due to missing 1/cos(theta) factor
theta = pinch(result_unpol_red.theta); #in degrees
frequency = f_red;
#randreset(12345); # set the seed if reproducibility is required
zemax_rayset_write_v1;

# Green
?Total_radiated_power = Power_Green/Power_Red; # in Watts,
export_filename = "zemax_rayset_green.txt";
radiant_intensity = result_unpol_green.radiance; # note that radiance may not be correct here due to missing 1/cos(theta) factor
theta = pinch(result_unpol_green.theta); #in degrees
frequency = f_green;
#randreset(12345); # set the seed if reproducibility is required
zemax_rayset_write_v1;

# Blue
?Total_radiated_power = Power_Blue/Power_Red; # in Watts,
export_filename = "zemax_rayset_blue.txt";
radiant_intensity = result_unpol_blue.radiance; # note that radiance may not be correct here due to missing 1/cos(theta) factor
theta = pinch(result_unpol_blue.theta); #in degrees
frequency = f_green;
#randreset(12345); # set the seed if reproducibility is required
zemax_rayset_write_v1;