#########################################################
# file: run_monte_carlo_script.lsf
#
# Description: This script file
#   is used to generate and run the Monte Carlo analysis
#   and plot the results at the end
#
# Copyright 2016, Lumerical Solutions, Inc.
#########################################################

closeall;
clear;

# add a new Monte Carlo analysis and set basic properties
addsweep(2);
MC_name = "MC_script";
?setsweep("Monte Carlo analysis");
setsweep("Monte Carlo analysis", "name", MC_name);
setsweep(MC_name, "number of trials", 50);
setsweep(MC_name, "enable seed", 1);
setsweep(MC_name, "seed", 1);
setsweep(MC_name, "Variation", "Both");

# define the parameter cpl
cpl2 = struct;
cpl2.Name = "cpl_2";
cpl2.Parameter = "::Root Element::WC2::coupling coefficient 1";
cpl2.Value = getnamed("WC2", "coupling coefficient 1");
dis = struct;
dis.type = "gaussian";
dis.variation = 0.02;
cpl2.Distribution = dis;

cpl3 = struct;
cpl3.Name = "cpl_3";
cpl3.Parameter = "::Root Element::WC3::coupling coefficient 1";
cpl3.Value = getnamed("WC3", "coupling coefficient 1");
dis = struct;
dis.type = "uniform";
dis.variation = 0.04;
cpl3.Distribution = dis;

addsweepparameter(MC_name, cpl2);
addsweepparameter(MC_name, cpl3);

# define the model ng
wgd_model = struct;
wgd_model.Name = "wgd_model";
wgd_model.Model = "WGD::group index 1";
wgd_model.Value = getnamed("SW1", "group index 1");
global = struct;
global.type = "uniform";
global.variation = 1.02;
wgd_model.Global = global;
local = struct;
local.type = "gaussian";
local.variation = 0.1;
wgd_model.Local = local;

addsweepparameter(MC_name, wgd_model);

# define the correlation 
wgd_corr = struct;
wgd_corr.Name = "wgd_corr";
wgd_corr.Parameters = "SW1_group_index_1,SW2_group_index_1,SW3_group_index_1";
wgd_corr.Value = 0.95;

addsweepparameter(MC_name, wgd_corr);

# define results
fsr = struct;
fsr.Name = "fsr";
fsr.Result = "::Root Element::Optical Network Analyzer::input 2/mode 1/peak/free spectral range";
fsr.Estimation = true;
fsr.Min = 1e11;
fsr.Max = 2e11;

bd = struct;
bd.Name = "bd";
bd.Result = "::Root Element::Optical Network Analyzer::input 1/mode 1/peak/bandwidth";
bd.Estimation = false;

gain = struct;
gain.Name = "gain";
gain.Result = "::Root Element::Optical Network Analyzer::input 1/mode 1/peak/gain";
gain.Estimation = false;

# add Monte Carlo results
addsweepresult(MC_name, fsr);
addsweepresult(MC_name, bd);
addsweepresult(MC_name, gain);

# run the Monte Carlo analysis
runsweep(MC_name);

# get & view Monte Carlo results
fsr = getsweepresult(MC_name, "analysis/results/histogram/fsr");
fsrCount = fsr.count;
fsr = fsr.fsr;

pdf = getsweepresult(MC_name, "analysis/results/pdf/fsr");
pdfCount = pdf.count;
pdfFsr = pdf.fsr;


histc(fsr/1e9, fsrCount, "free spectral range (GHz)", "count", "histogram");
legend("");
plot(pdfFsr/1e9, pdfCount, "fsr (GHz)", "count", "probability density function");
legend("");