local_jcm_path = getenv('JCMROOT'); addpath(fullfile(local_jcm_path, 'ThirdPartySupport', 'Matlab'));
jcmwave_startup; jcmwave_set_num_threads(1); jcmwave_daemon_shutdown();
jcmwave_daemon_add_workstation('Multiplicity', 2, 'NThreads', 1); 

% geometry
keys.uol = 1.e-9;                                
keys.radius_scatterer = 400.0;
keys.contains_substrate = true;
keys.offset_sphere = -10.0; % offset between sphere and substrate in uol

% refractive indices
keys.n_1 = 1.0; % background
keys.n_2 = 3.0 + 0.1i; % sphere
keys.n_3 = 1.5; % substrate

% S&P plane waves incidence angle
keys.theta = 20;

keys.fem_degree = 3;
keys.n_steps = 0;

c0 = 299792458; mu0 = 4*pi*1e-7; eps0 = 1/(mu0*c0^2); Z0 = sqrt(mu0/eps0);
geo_cross_section = pi*(keys.radius_scatterer*keys.uol)^2; p_in = 0.5*keys.n_1/Z0*geo_cross_section; 

% parameter scan
job_ids = []; results = []; simulation_results = []; counter = 0;

wavelengths = [400:2:600]*1e-9;

for ii = 1:length(wavelengths)
  counter = counter + 1;
  keys.vacuum_wavelength = wavelengths(ii);

  job_ids(end + 1) = jcmwave_solve('project.jcmp', keys, 'workingdir', fullfile('tmp', ['s_' num2str(counter, '%05i')]));

  results(counter, 1) = keys.vacuum_wavelength;
  results(counter, 6) = keys.fem_degree;
  results(counter, 7) = keys.n_steps;
end;

[simulation_results, logs] = jcmwave_daemon_wait(job_ids);

for ii = 1:length(job_ids)
  this_result = simulation_results{ii};
  field_energy = cell2mat(this_result{2}.ElectricFieldEnergy);
  omega = 2*pi*c0/keys.vacuum_wavelength;         
  absorption = -2*omega*imag(field_energy(2, :))/p_in;
  flux_scat = cell2mat(this_result{4}.ElectromagneticFieldEnergyFlux);
  scattering = sum(real(flux_scat), 1)/p_in;
  results(ii, 2:3) = absorption;
  results(ii, 4:5) = scattering;
end

filename = ['results_' datestr(now,'yyyymmdd_HHMMSS') '.txt']; 
save ('-ascii', '-double', filename, 'results'); copyfile(filename, 'results.txt');
display_results;
