clear

matlabPath = getenv('LD_LIBRARY_PATH');
libPath = getenv('PATH');

setenv('LD_LIBRARY_PATH', libPath);
system('xmds2 eit.xmds');
setenv('LD_LIBRARY_PATH', matlabPath);

decay_bc = 0.00001;
decay_ab = 6;

widths = [];
absorptions = [];
drives = logspace(0, -2, 100);
drives = drives(1:1);
%drives = linspace(0.0033154997,0.00331549970606061,100);
detunings1 = [-15];
dephase_bc = [0];

for drive = drives
    width_row = [];
    absorption_row = [];
    for detuning1 = detunings1
        for dephasing = dephase_bc

            expected_width = (drive^2/decay_ab + dephasing + decay_bc);
            domain_max = expected_width*10;
            domain_min = -1*domain_max;
            final_time = 1/expected_width*100;
            run_eit = sprintf('time ./eit --drive_arg=%0.20f --decay_arg=%f --decay_bc_arg=%f --delta1_arg=%f --domain_min=%f --domain_max=%f --dephase_bc_arg=%f --final_time=%f', drive, decay_ab, decay_bc, detuning1, domain_min, domain_max, dephasing, final_time);
            
            setenv('LD_LIBRARY_PATH', libPath);
            system(run_eit);
            system('xsil2graphics2 -m eit.xsil');
            setenv('LD_LIBRARY_PATH', matlabPath);
            %[drive, detuning1] % Use for debugging
            run('eit.m');
            detunings2 = detuning_1';
            PabI = PabI_1(end,:)';
            [fwhm, absorption_min] = getPeak(detunings2, PabI);
            
            width_row = [width_row, fwhm];
            absorption_row = [absorption_row, absorption_min];
            eit_plot = [detunings2, PabI];
            csv_name = sprintf('data/full/%0.10f_%0.20f_%0.10f_%0.10f.dat', decay_bc, drive, detuning1, dephasing);
            csvwrite(csv_name, eit_plot);
        end
    end
    widths = [widths; width_row];
    absorptions = [absorptions; absorption_row];
end


% plot(drives.^2, widths, drives.^2, absorptions*10^5)
% legend('width', 'absorption')
% xlabel('Drive^2')

% figure
% surf(detunings1, drives, widths)
% xlabel('detuning')
% ylabel('drive')
% zlabel('width')

% figure
% surf(detunings1, drives, absorptions)
% xlabel('detuning')
% ylabel('drive')
% zlabel('absorption')

% plot(detunings2, PabI)
% xlabel('Detuning')
% ylabel('Absorption')
% title(['Drive equals ' [int2str(drive)]])


save(['workspaces/workspace_', datestr(clock,0), '.mat']);