function calculate_Jsc(layer_materials, layer_thicknesses, absorbing_layers, f, res) {
    # layer_materials: Material database names of the layer materials
    # layer_thicknesses: Thicknesses of the layers, in m
    # absorbing_layers: Index of the absorbing layer(s) which contribute to the current
    # f: vector of frequencies
    
    num_layers = length(layer_thicknesses);
    
    # add air layers above and below solar cell
    d = [0; layer_thicknesses; 0];
    
    # create index matrix
    n = matrix(length(d),length(f)); 
    n(1,1:length(f)) = 1; 
    for (ii=2:num_layers+1) {
        n(ii,1:length(f)) = getfdtdindex(layer_materials{ii-1}, f, min(f), max(f));
    } 
    n(num_layers+2,1:length(f)) = 1; 
    
    theta = 0; # angle of incidence

    field = stackfield(n,d,f,theta,res);
    E = pinch(field.Es(1, 1, :, :, 1, 2));
    x=field.x; y=field.y; z=field.z;
    
    # get index as a function of z and f
    # only add absorbing layers to index, so only those layers 
    # are included in the absorption integrals
    index = matrix(length(z), length(f));
    for (ii=absorbing_layers) {
        layer_min = sum(d(1:ii)); # z location of top of layer
        layer_max = sum(d(1:ii+1)); # z location of bottom of layer
        layer_mask = z >= layer_min and z < layer_max;
        index = index + meshgridy(z, n(ii+1, :))*meshgridx(layer_mask, f);
    }
    # Code below here is from solar_generation analysis group
    # Get the AM1.5 solar spectrum
    lam_nm = solar(0)*1e9;       # solar spectrum wavelength vector, in units of nm.
    Psolar  = solar(1)*1e-9;     # solar power spectrum, in units of Watts/m^2/nm
    f_solar = 1e9*c/lam_nm;      # solar spectrum frequency vector, in units of Hz
    Nf_solar = length(f_solar);
    
    # select the region of solar spectrum covered by the monitor data
    fl = max([min(f_solar),min(f)]);
    fh = min([max(f_solar),max(f)]);
    fi = find((f_solar>=fl)&(f_solar<=fh));
    f_solar = f_solar(fi);
    lam_nm = lam_nm(fi);
    Psolar = Psolar(fi);
    Nf_solar = length(f_solar);
    stack_sourceintensity = 0.5*c*eps0; # very similar value to sourceintensity
    nm = Psolar/stack_sourceintensity;   # units of 1/nm
    nm = meshgridy(z,nm);
    
    # Calculate number of absorbed photon per unit volume
    # assume this is equal to the number of generated electron/hole pairs
    g = 0.5*abs(E)^2*imag(eps0*index^2)/hbar;
    # interpolate to solar frequency vector
    g = interp(g, z, f, z, f_solar);
    # Calculate the generation rate by integrating 'g' over wavelength
    # The generate rate is the number of electron hole pairs per unit volume per second.  (units: 1/m^3/s)
    G_matrix = integrate2(g*nm,2,lam_nm);
    
    # Calculate the short circuit current (Jsc)
    Jsc = e*integrate(G_matrix,1,z); # A/m^2
    
    return Jsc;
}

clear;
AR_sweep = false; # set to true to perform sweep of GaAs AR layer thickness
importmaterialdb("solar_cell_materials.mdf"); # import custom materials required to match examples

# STACK silicon solar cell simulation
layer_materials = {"Si (Silicon) - Palik", "Al (Aluminium) - Palik"};
layer_thicknesses = [3e-6; 0.5e-6]; # 3 um of Si, 0.5 um of Al
absorbing_layers = [1]; # Si is the absorbing material
fmin = c/1.1e-06;
fmax = c/3e-07;
f = linspace(fmin,fmax,201);
res = 1000;

Jsc_Si = calculate_Jsc(layer_materials, layer_thicknesses, absorbing_layers, f, res);
?"Short circuit current for Si solar cell: " + num2str(Jsc_Si) + " A/m^2    ("+num2str(Jsc_Si/10)+" mA/cm^2)";

# GaAs solar cell simulation
layer_materials = {"ARC", "Al(0.8)Ga(0.2)As", "GaAs - Palik Copy 1", "Al(0.3)Ga(0.7)As", "Al (Aluminium) - Palik"};
layer_thicknesses = [0.1e-6; 0.03e-6; 1.65e-6; 0.02e-6; 0.5e-6];
absorbing_layers = [2,3,4]; # AlGaAs and GaAs layers are absorbing
fmin = c/1.3e-06;
fmax = c/3e-07;
f = linspace(fmin,fmax,201);
res = 7500; # higher resolution required for stacks with thin layers

Jsc_GaAs = calculate_Jsc(layer_materials, layer_thicknesses, absorbing_layers, f, res);
?"Short circuit current for GaAs solar cell: " + num2str(Jsc_GaAs) + " A/m^2    ("+num2str(Jsc_GaAs/10)+" mA/cm^2)";

if (AR_sweep) {
    # sweep AR layer thickness
    AR_thickness = linspace(0, 0.45e-6, 46);
    Jsc_AR_sweep = zeros(length(AR_thickness));
    
    for (ii=1:length(AR_thickness)) {
        layer_thicknesses = [AR_thickness(ii); 0.03e-6; 1.65e-6; 0.02e-6; 0.5e-6];
        Jsc_AR_sweep(ii) = calculate_Jsc(layer_materials, layer_thicknesses, absorbing_layers, f, res);
        ?"Sweep " + num2str(round(ii/length(AR_thickness)*100)) + "% complete.";
    } 
    
    plot(AR_thickness*1e9, Jsc_AR_sweep, "AR Layer Thickness (nm)", "Jsc (A/m^2)", "", "linewidth=3");
    setplot("x min", min(AR_thickness)*1e9); setplot("x max", max(AR_thickness)*1e9);
}