% This script generates is a shortened version of the localization example
% in ../ex1_localization/, which (i) generates simulated data with a single
% diffusing fluorophore in it, and writes the results to result_trj_mv/
% (old results will be overwritten)
%
% Requirements: 
% - matlab R2015a or later, with 
%   image_toolbox
%   optimization_toolbox
%   statistics_toolbox
% - the ../psfModels/ folder must be in the matlab path, e.g.,
%       addpath(genpath(fullfile(pwd,'..','psfModels')));
% - SMeagol must be installed and on the matlab path. See 
%   http://smeagol.sourceforge.net/
%
%
% ML 2016-12-27

%% parameters & output
Nphotons=200 % photons/spot
bg=1         % background photons/pixel in each frame
tE=10e-3     % image exposure time
psfIndex=2   % choice of PSF model for the analysis, chosen from the list
%  1 : RW_514nm_NA149
%  2 : RW_639nm_NA14
%  3 : RW_680nm_NA14

%% mesord simulation to generate the molecule trajectory
outDir=fullfile(pwd,'result_trj_mv');
try
    system(['rm -rf ' outDir])
catch
end
system(['mkdir  ' outDir])

% mesord parameters
modelPath =fullfile(pwd,'..','mesoRDmodel_D100_longCells','Ecoli_1state_D100.xml');
mesordOptions=['  ' modelPath ' -K -g -x 0 -c 1 -t ' num2str(Nspots*dt+0.2) ' -i 1 -I 50 -w 0.0037 -q 0.020 um -c 5 '];
SM_runMesoRD(mesordOptions,outDir)

%% SMeagol simulation to generate the movie
% first, we need to modify the SMeagol parameters somewhat
SMopt=SM_getOptions('runinput0_SMeagol.m');

SMopt.background.photons_per_pixel=bg;

SMopt.sample.dt=dt;
SMopt.sample.tE=tE;
SMopt.baseIntensity.intensity=Nphotons/tE;
SMopt.camera.alpha=1/EMgain;
SMopt.camera.sigmaReadout=readRMS;

% PSF model
psfFile_rel={...
'../psfModels/psf_RW_514nm_NA149_ni1515.mat',...
'../psfModels/psf_RW_639nm_NA14_ni1515.mat',...
'../psfModels/psf_RW_680nm_NA14_ni1515.mat'};
psfString={...
'RW_514nm_NA149',...
'RW_639nm_NA140',...
'RW_680nm_NA140'};
PSFfile=fullfile(pwd,psfFile_rel{psfIndex});
SMopt.psf.type = 'SM_psf_numericalPSF_Ruzp';
SMopt.psf.PSFfile=PSFfile;

% output file
outDir_relative=SM_relative_path(pwd,outDir);
SMopt.trj.trajectoryFile=fullfile(outDir_relative,'trajectories.txt');
SMopt.trj.reactionFile  =fullfile(outDir_relative,'reactions.txt');

SMopt.output.resultFile=fullfile(outDir_relative,'SMeagol_result.mat');
SMopt.output.maxFrames=Nspots;

SMopt=SM_runsimulation(SMopt);
R=SM_getResults(SMopt);

disp(['Done: result written to ' outDir_relative])
