#########################################################
# file: usr_temporal_coherence.lsf
#
# Description:
# This file calculates the incoherent reflection spectrum between
# 500 and 700 nm for a Lorentz function source spectrum using
# the cw result from temporal_coherence.fsp for coherence lengths
# Lc = 15 um and 300 um.
#	
# Copyright 2012, Lumerical Solutions, Inc.
#########################################################

# get cw reflection and frequency vector
f=getdata("R","f");
w = 2*pi*f;
R=-transmission("R");

tc_values = [0.5e-13,1e-12]; # set coherence times corresponding to coherence lengths 15 and 300 um

for(tc = tc_values) { # loop through coherence times
    f0 = linspace(c/700e-9,c/500e-9,1000); # set frequency vector for incoherent results
    w0 = 2*pi*f0;
    delta = 1/tc;  

    R0 = matrix(length(f0),length(tc)); # initialize result matrix
    
    # perform the incoherent weighting
    for (i=1:length(f0)) {
        # define a Lorentz function as the spectrum of the physical source
        h2 = delta / ( (w - w0(i))^2 + (pi*delta)^2 ); 

        # integral of the cw reflection weighted by the source spectrum divided by the 
        # integral of the source spectrum.
        R0(i) = integrate( h2*R,1,w) /
                integrate( h2  ,1,w) ;
    }

    # plot monochromatic and calculated incoherent reflection spectrum
    plotxy(c/f*1e9,R,c/f0*1e9,pinch(R0,2,1),"wavelength( nm)","reflection");
    legend("monochromatic","Lc = "+num2str(round(c*tc*1e6)) + " microns");

    # plot source spectrum at center value of w0
    w0 = 0.5*(min(w0)+max(w0));
    h2 = delta / ( (w - w0)^2 + (pi*delta)^2 );
    h2 = h2/max(h2); # normalize spectrum to peak value=1
    plot(c/f*1e9,h2,"wavelength (nm)","source spectrum","spectrum, Lc = "+num2str(round(c*tc*1e6)) + " microns");
}