###########################################################################
# Scriptfile: oled.lsf
#
# Description:
# 	This file is used to calculate the incoherent angular emission
#	from an istropic ensemble of dipole emitters for oled applications.
#	For each structure it calculates an x, y and z oriented dipole 
#	and then adds |E|^2 in the far field to calculate the response of 
#	isotropic dipoles.
#
#	The data is then exported to an ASAP *.inr file as a rayset so 
#       that it can be used as a source to design the packaging and 
#	other macroscopic components of an oled.
#
#	The file is organized into 3 sections. Each section can
#	be run independently
#		1) Setup and running of FDTD simulations
#		2) Analysis of FDTD simulations
#		3) Export of results to ASAP *.inr file
#
#
#	Currently, there is an assumption that the position of the dipole
#	in the x-y plane is not important, and the structure must
#	be symmetric with respect to reflection in the x=0 plane and the
#	y=0 plane. The script could be modified to break the assumptions
#	of symmetry and to model several dipole positions with respect
#	to any patterning of the layers.
# 
#
# Copyright 2007, Lumerical Solutions, Inc.
###########################################################################

# Choose which components to run (1 for yes, 0 for no)
# The first time, each component needs to be run
rerun_simulations = 1;
rerun_analysis    = 1;
rerun_asap_export = 1;

# Choose the name of the template fsp file. In this example, the file should
# either oled1 or oled2, but can be changed to run the user's own file.
basefilename = "oled1";

# Choose the size and position of the emission region when exporting to ASAP
x1_asap = -0.25e-3;
x2_asap = 0.25e-3;
y1_asap = -0.25e-3;
y2_asap = 0.25e-3;
z_asap = 0;

# Choose how many copies of the original rayset to export
N_export = 10;


############################################################
# Setup and running of FDTD simulations
############################################################
# rerun the simulations if requested
if(rerun_simulations) {
    # Choose how the simulations are run
    #	1: standard with graphical interface
    #	2: no graphical window
    #	3: parallel
    runtype = 3;

    # load the desired template file
    load(basefilename+".fsp");
 
    # x dipole simulation, setup appropriate boundary conditions
    # and orient the dipole appropriately
    switchtolayout;
    setnamed("FDTD","x min bc","Anti-Symmetric");
    setnamed("FDTD","y min bc","Symmetric");
    setnamed("Source1","theta",90);
    setnamed("Source1","phi",0);
    save(basefilename+"_x.fsp");
    run(runtype);

    # y dipole simulation, setup appropriate boundary conditions
    # and orient the dipole appropriately
    switchtolayout;
    setnamed("FDTD","x min bc","Symmetric");
    setnamed("FDTD","y min bc","Anti-Symmetric");
    setnamed("Source1","theta",90);
    setnamed("Source1","phi",90);
    save(basefilename+"_y.fsp");
    run(runtype);

    # z dipole simulation, setup appropriate boundary conditions
    # and orient the dipole appropriately
    switchtolayout;
    setnamed("FDTD","x min bc","Symmetric");
    setnamed("FDTD","y min bc","Symmetric");
    setnamed("Source1","theta",0);
    setnamed("Source1","phi",0);
    save(basefilename+"_z.fsp");
    run(runtype);
}

############################################################
# Analysis of the results of the FDTD simulations
############################################################

if(rerun_analysis) {
    # Ensure that CW Normalization is used during the analysis
    # This should always be used when exporting data to ASAP
    cwnorm;

    # Choose a far field filter setting to reduce erroneous ripple
    # in the far field data due to the finite size of the simulation
    # in the x-y plane
    farfieldfilter(0.25);

    # Choose a resolution for the far field calculation
    # The far field data will have size resXres
    # The number should be odd to get a data point at precisely normal (theta=0)
    res = 201;

    # Choose the frequency point for the calculation
    # In the example files, we had 50 frequency points from 400 to 700 nm
    # Point 27 corresponds to approximately 500 nm
    fpoint = 27;

    # set the name of the monitor that records the field transmitted to the substrate
    monname = "substrate_monitor";

    # load the x-dipole results
    load(basefilename+"_x.fsp");
        # calculate the far field and near fields
        E2_near = getelectric(monname); # returns |E|^2 in near field
        E2_far = farfield3d(monname,fpoint,res,res); # return |E|^2 at 1 m distance

    # load the y-dipole results
    load(basefilename+"_y.fsp");
        # calculate the far field and near fields
        E2_near = E2_near + getelectric(monname); # add to previous result
        E2_far = E2_far + farfield3d(monname,fpoint,res,res); #add to previous result

    # load the z-dipole results
    load(basefilename+"_z.fsp");
        # calculate the far field and near fields
        E2_near = E2_near + getelectric(monname); # add to previous result
        E2_far = E2_far + farfield3d(monname,fpoint,res,res); #add to previous result

    # select the appropriate frequency point for the near field data
    E2_near = pinch(E2_near,4,fpoint);

    # normalize the far field distribution to have a maximum value of 1
    E2_far = E2_far/max(E2_far);

    # get the vector of frequencies recorded
    f = getdata(monname,"f");

    # output to the user the current wavelength of the calculation based on
    # the choise of fpoint
    ?"working at a wavelength of: " + num2str(c/f(fpoint)*1e9) + " nm";

    # collect the x and y vectors for the near field data
    x = getdata(monname,"x");
    y = getdata(monname,"y");

    # image the near field results on linear and logarithmic scale
    image(x*1e6,y*1e6,E2_near,"x (microns)","y (microns)","|E|^2 vs x,y near field");
    image(x*1e6,y*1e6,E2_near,"x (microns)","y (microns)","|E|^2 vs x,y near field, logscale","logplot");

    # get the direction cosines for the far field projections
    ux = farfieldux(monname,fpoint,res,res);
    uy = farfielduy(monname,fpoint,res,res);

    # image the far field results (note the polar option in the plot options)
    image(ux,uy,E2_far,"","","|E|^2 in far field","polar");
}

############################################################
# Export the results to ASAP
############################################################

if(rerun_asap_export) {
    # set the output file name to the same as basefilename with the *.inr extension
    filename = basefilename+".inr";

    # make sure the file is deleted if it already exists
    rm(filename);

    # create a 2D mesh of direction cosines
    Ux = meshgridx(ux,uy);
    Uy = meshgridy(ux,uy);
    Uz = real(sqrt(1-Ux^2-Uy^2));

    # create a filter to remove rays at imaginary angles
    filter = Ux^2 + Uy^2 < 0.99;

    # add the asap command EMITTING RAYS
    write(filename,"emitting rays");

    # initialize a counter of the total number of exported rays
    ray_counter = 0;
    
    # initialize a counter to output progress to user
    counter = 0;
for(ii=1:N_export) {
    
    # output progress to screen
    ?"Completed " + num2str(counter/N_export*100) +"%";
    counter = counter+1;  

    # reserve a matrix for exporting the entire rayset (once)
    outresult = matrix(length(ux)*length(uy),7);

    # fill in the matrix with the correct data
    pos=0;
    for(i=1:length(ux)) {
        for(j=1:length(uy)) {
            if(filter(i,j)) {
	        ray_counter = ray_counter+1;
	        pos = pos+1;
	        xrand = rand*(x2_asap-x1_asap) + x1_asap;
	        yrand = rand*(y2_asap-y1_asap) + y1_asap;
                outresult(pos,1) = xrand;
                outresult(pos,2) = yrand; 
                outresult(pos,3) = z_asap;
                outresult(pos,4) = Ux(i,j);
                outresult(pos,5) = Uy(i,j);
                outresult(pos,6) = Uz(i,j);
		# divide the ray flux by cos(theta) to account for the
		# the uneven angular sampling
                outresult(pos,7) = E2_far(i,j)/Uz(i,j);
            }
        }
    }
 
    # cut down the result to remove zeros at the end
    outresult = outresult(1:pos,1:7);

    # write the rayset to the file
    write(filename,num2str(outresult));
    }

    # write the ASAP RETURN command
    write(filename,"return");

    # inform the user of the total number of exported rays
    ?"output: " + num2str(ray_counter) + " rays";
}