example_custom_input_fields_opg.m

% Script demonstrating OPG -> OPA using snlo_opg_func opg stage has a short
% crystal, the red2 wave is dumped between stages. The leftover blue wave from
% the first stage is used as pump in second stage -- though this is not a great
% idea it is easiest is simple. A better scheme is to make a good portion of the
% blue wave bypass the first stage & use that for second stage input pump.
%

%% All the inputs

pump_pulse_energy = linspace(1e-9,1000e-9,15); % make vector of equally spaced pulse energies between 1 nJ and 1000 nJ, we'll use these as energy of blue wave in opg model
pump_pulse_energy = flip(pump_pulse_energy);
% Define the full set of input parameters for opg; these correspond to the input
% boxes visible when you run opg. for most parameters, a row vector with three
% elements is defined and corresponds to the values for the red1, red2, and blue
% waves. Calling the opg function with a set of input parameters requires
% selecting a single set of inputs and sending it to the snlo_opg_func as an
% intput argument.

% lengths of crystals
L_crystal_OPG = 0.33; % millimeters
L_crystal_OPA = 0.50; % millimeters

inputs.opg_2d_wavelengths 			= [2500,1852.4,1064];     			% wavelengths (in nm)
inputs.opg_2d_ref_inds 				= [1.6273, 1.6411, 1.6353];    	% refractive indices
inputs.opg_2d_gvi 					= [1.6902, 1.6746, 1.6530];     % group velocity indices
inputs.opg_2d_gdd 					= [-380.9, -90.4, 41.6];      	% group delay dispersions (in fs^2/mm)
inputs.opg_2d_phase 				= [0,0,0];                   	% input phases (in radians)
inputs.opg_2d_input_refl 			= [0,0,0];              		% crystal input face power reflectivity coefficients (0-1)
inputs.opg_2d_output_refl 			= [0,0,0];             			% crystal output face power reflectivity coefficients (0-1)
inputs.opg_2d_crystal_losses 		= [0,0,0];          			% linear crystal losses (in 1/mm)
inputs.opg_2d_n2_red1 				= [0,0,0];                 		% red1 nonlinear refractive indices (in sq cm/W)
inputs.opg_2d_n2_red2 				= [0,0,0];                 		% red2 nonlinear refractive indices (in sq cm/W)
inputs.opg_2d_n2_blue 				= [0,0,0];                 		% blue nonlinear refractive indices (in sq cm/W)
inputs.opg_2d_beta_red1 			= [0,0,0];               		% red1 two photon absorption coefficients (in cm/W)
inputs.opg_2d_beta_red2 			= [0,0,0];               		% red2 two photon absorption coefficients (in cm/W)
inputs.opg_2d_beta_blue 			= [0,0,0];               		% blue two photon absorption coefficients (in cm/W)
inputs.opg_2d_pulseenergy 			= [1e-7,1e-7,0];          		% input pulse energies (in J)
inputs.opg_2d_pulse_durations 		= [0.1,0.1,0.0750;1,1,1].'; 	% input pulse durations (in fwhm ps); with temporal supergaussian coefficients (positive integers)
inputs.opg_2d_pulse_delays 			= [0,0];              			% temporal delays relative to blue pulse (in ps)
inputs.opg_2d_pulse_chirps 			= [0,0,0];            			% linear chirps (in THz/ps)
inputs.opg_2d_beam_diameters 		= [0.15,0.15,0.10625;0.15,0.15,0.10625].'; 	% input beam diameters (in fwhm mm), in walkoff direction; and perpendicular to walkoff direction
inputs.opg_2d_supergaussian_coeff 	= [1,1,1];     					% spatial supergaussian coefficients (integers 1-10)
inputs.opg_2d_wo_angles 			= [0,0,0];               		% spatial walkoff angles (in milliradians)
inputs.opg_2d_offset_wodir 			= [0,0];              			% spatial offsets in walkoff directions (in mm)
inputs.opg_2d_rad_curv 				= [1e7,1e7,66.983;1e7,1e7,66.983].';  % input beam radii of curvature (in mm in air) for [red1,red2,blue] and [in the walkoff direction;parallel to the walkoff direction]
inputs.opg_2d_nt 					= 350;                          % number of time points to consider
inputs.opg_2d_nxny 					= [100,100];                    % number of spatial grid points in the walkoff direction, and perpendicular to walkoff direction
inputs.opg_2d_grid_duration 		= 0.3;               			% total time to model (in ps)
inputs.opg_2d_crystal_length 		= L_crystal_OPG;              	% length of nonlinear crystal (in mm)
inputs.opg_2d_lx_ly 				= [1,1];                     	% size of spatial grid (in mm) in walkoff direction, and perpendicular to walkoff direction
inputs.opg_2d_deff 					= 30;                         	% crystal nonlinear coefficient (in pm/V)
inputs.opg_2d_deltak 				= 0;                        	% phase mismatch (in rad/mm)
inputs.opg_2d_nz 					= 100;                          	% number of (z) integration steps through the crystal

%% call 2d mix sp, run model, load output
clear problem opa_problem

powers              = cell(size(pump_pulse_energy));
opg_output_energies = zeros(length(pump_pulse_energy),3);
opa_output_energies = zeros(length(pump_pulse_energy),3);

for K = 1:length(pump_pulse_energy)
    fprintf(1,' Running %i of %i: %g nJ\n', K, length(pump_pulse_energy), pump_pulse_energy(K)*1e9);
    %% first stage
    problem(K)                          = inputs; % copy the inputs specified earlier as the starting point for this problem's inputs
    problem(K).opg_2d_pulseenergy(1:2)  = 0;                     
    problem(K).opg_2d_pulseenergy(3)    = pump_pulse_energy(K); 
    fcn_handles                         = snlo_opg_func(problem(K));
    run_handle                          = fcn_handles{1};    % the function handle of the 'Run' button callback is the first returned; the 'run' button callback function is always the first element of the cell array of returned function handles
    accept_handle                       = fcn_handles{2}; % the function handle of the 'Run' button callback is the first returned; the 'accept' button callback function is always the second element of the cell array of returned function handles, if the function has an 'accept' button
    close_handle                        = fcn_handles{end};% the function handle of the 'Close figure' callback is last returned
    accept_handle();                    % call the accept button callback - equivalent to clicking the accept button
    run_handle();                       % call the 'Run' button callback - equivalent to clicking the run button
    
    data                                = load('opg_2d_output.mat');
    powers{K}                           = data.power;                           % assign cell to contain powers, an nt*4 array; in the second dimension, the first column is time, second is red1 power, third is red2 power, and fourth is blue power
    opg_output_energies(K,1)            = trapz(powers{K}(:,1),powers{K}(:,2)); % Red1 energy - integrate power
    opg_output_energies(K,2)            = trapz(powers{K}(:,1),powers{K}(:,3)); % Red2 energy - integrate power
    opg_output_energies(K,3)            = trapz(powers{K}(:,1),powers{K}(:,4)); % Blue energy - integrate power
    fprintf(1, '  Stage 1 output energies: %g, %g, %g (total %g) nJ\n', opg_output_energies(K,:)*1e9, sum(opg_output_energies(K,:)*1e9));

    %% second stage - take the opg output, amplify in opa stage -- dump red2 wave

    clear('opa_input_fields');
    opa_input_fields.xmat = data.xgridmat;
    opa_input_fields.ymat = data.ygridmat;
    opa_input_fields.tvec = data.tgrid;
    opa_input_fields.input_field_xyt_red1   = data.field_red1;
    opa_input_fields.input_field_xyt_red2   = data.field_red2*1e-9;  % dump the red2 wave's light
    opa_input_fields.input_field_xyt_blue   = data.field_blue;
    if K == 1
        opa_problem(K)                      = problem(K);
    else
        opa_problem(K)                      = opa_problem(1);
    end
    opa_problem(K).opg_2d_pulseenergy       = opg_output_energies(K,:);
    opa_problem(K).opg_2d_pulseenergy(2)    = opa_problem(K).opg_2d_pulseenergy(2) * 1e-9;
    opa_problem(K).opg_2d_crystal_length    = L_crystal_OPA;
    opa_problem(K).custom_fields            = true;
    opa_problem(K).input_fields             = opa_input_fields;
    new_fcn_handles                         = snlo_opg_func(opa_problem(K));
    run_handle      = new_fcn_handles{1};
    accept_handle   = new_fcn_handles{2};
    accept_handle();
    run_handle();
    % keyboard;
    opa_data = load('opg_2d_output.mat');
    opa_output_energies(K,:) = [trapz(opa_data.power(:,1),opa_data.power(:,2)), ...
        trapz(opa_data.power(:,1),opa_data.power(:,3)), ...
        trapz(opa_data.power(:,1),opa_data.power(:,4))];
    fprintf(1, '  Stage 2 output energies: %g, %g, %g (total %g) nJ\n', opa_output_energies(K,:)*1e9, sum(opa_output_energies(K,:)*1e9));
end


%% plot output energies as a function of input energies
fh      = figure; % make new output figure
ax      = axes('Parent',fh); % make new axes in the new figure
ph      = plot(ax, pump_pulse_energy, opa_output_energies(:,1), 'r-', ...
            pump_pulse_energy, opa_output_energies(:,2), 'b-', ...
            pump_pulse_energy, opa_output_energies(:,3), 'g-');
xlabel(ax, 'Input blue pulse energies [J]');
ylabel(ax, 'Output pulse energies [J]');
legend(ax, 'Red1', 'Red2', 'Blue');