From addbb92ff06b21d2188eeefdb17692beb85e4b1c Mon Sep 17 00:00:00 2001 From: Eugeniy Mikhailov Date: Tue, 22 Apr 2014 11:48:18 -0400 Subject: Hunter's initial release --- eit.xmds | 136 ++++++++++++++++++++++++++++++++++++++++++ extrema.m | 146 +++++++++++++++++++++++++++++++++++++++++++++ extrema2.m | 174 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ getPeak.m | 42 +++++++++++++ loadSimulation.m | 7 +++ loadSimulations.m | 27 +++++++++ simulate.m | 74 +++++++++++++++++++++++ 7 files changed, 606 insertions(+) create mode 100644 eit.xmds create mode 100644 extrema.m create mode 100644 extrema2.m create mode 100644 getPeak.m create mode 100644 loadSimulation.m create mode 100644 loadSimulations.m create mode 100644 simulate.m diff --git a/eit.xmds b/eit.xmds new file mode 100644 index 0000000..50419f5 --- /dev/null +++ b/eit.xmds @@ -0,0 +1,136 @@ + + + + + Hunter Rew + + Model of 3 state atom + + + + + double + + + + + + + + + + + + + + + + + + + + t + + + + + + + + r_ab r_cb + + + + + + + Paa Pbb Pcc Pab Pca Pcb + + + + + + + + + + 10 + + population + + Gamma + + + + + + + + + PabI + population + + + + diff --git a/extrema.m b/extrema.m new file mode 100644 index 0000000..2a9d265 --- /dev/null +++ b/extrema.m @@ -0,0 +1,146 @@ +function [xmax,imax,xmin,imin] = extrema(x) +%EXTREMA Gets the global extrema points from a time series. +% [XMAX,IMAX,XMIN,IMIN] = EXTREMA(X) returns the global minima and maxima +% points of the vector X ignoring NaN's, where +% XMAX - maxima points in descending order +% IMAX - indexes of the XMAX +% XMIN - minima points in descending order +% IMIN - indexes of the XMIN +% +% DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima): +% In mathematics, maxima and minima, also known as extrema, are points in +% the domain of a function at which the function takes a largest value +% (maximum) or smallest value (minimum), either within a given +% neighbourhood (local extrema) or on the function domain in its entirety +% (global extrema). +% +% Example: +% x = 2*pi*linspace(-1,1); +% y = cos(x) - 0.5 + 0.5*rand(size(x)); y(40:45) = 1.85; y(50:53)=NaN; +% [ymax,imax,ymin,imin] = extrema(y); +% plot(x,y,x(imax),ymax,'g.',x(imin),ymin,'r.') +% +% See also EXTREMA2, MAX, MIN + +% Written by +% Lic. on Physics Carlos Adrián Vargas Aguilera +% Physical Oceanography MS candidate +% UNIVERSIDAD DE GUADALAJARA +% Mexico, 2004 +% +% nubeobscura@hotmail.com + +% From : http://www.mathworks.com/matlabcentral/fileexchange +% File ID : 12275 +% Submited at: 2006-09-14 +% 2006-11-11 : English translation from spanish. +% 2006-11-17 : Accept NaN's. +% 2007-04-09 : Change name to MAXIMA, and definition added. + + +xmax = []; +imax = []; +xmin = []; +imin = []; + +% Vector input? +Nt = numel(x); +if Nt ~= length(x) + error('Entry must be a vector.') +end + +% NaN's: +inan = find(isnan(x)); +indx = 1:Nt; +if ~isempty(inan) + indx(inan) = []; + x(inan) = []; + Nt = length(x); +end + +% Difference between subsequent elements: +dx = diff(x); + +% Is an horizontal line? +if ~any(dx) + return +end + +% Flat peaks? Put the middle element: +a = find(dx~=0); % Indexes where x changes +lm = find(diff(a)~=1) + 1; % Indexes where a do not changes +d = a(lm) - a(lm-1); % Number of elements in the flat peak +a(lm) = a(lm) - floor(d/2); % Save middle elements +a(end+1) = Nt; + +% Peaks? +xa = x(a); % Serie without flat peaks +b = (diff(xa) > 0); % 1 => positive slopes (minima begin) + % 0 => negative slopes (maxima begin) +xb = diff(b); % -1 => maxima indexes (but one) + % +1 => minima indexes (but one) +imax = find(xb == -1) + 1; % maxima indexes +imin = find(xb == +1) + 1; % minima indexes +imax = a(imax); +imin = a(imin); + +nmaxi = length(imax); +nmini = length(imin); + +% Maximum or minumim on a flat peak at the ends? +if (nmaxi==0) && (nmini==0) + if x(1) > x(Nt) + xmax = x(1); + imax = indx(1); + xmin = x(Nt); + imin = indx(Nt); + elseif x(1) < x(Nt) + xmax = x(Nt); + imax = indx(Nt); + xmin = x(1); + imin = indx(1); + end + return +end + +% Maximum or minumim at the ends? +if (nmaxi==0) + imax(1:2) = [1 Nt]; +elseif (nmini==0) + imin(1:2) = [1 Nt]; +else + if imax(1) < imin(1) + imin(2:nmini+1) = imin; + imin(1) = 1; + else + imax(2:nmaxi+1) = imax; + imax(1) = 1; + end + if imax(end) > imin(end) + imin(end+1) = Nt; + else + imax(end+1) = Nt; + end +end +xmax = x(imax); +xmin = x(imin); + +% NaN's: +if ~isempty(inan) + imax = indx(imax); + imin = indx(imin); +end + +% Same size as x: +imax = reshape(imax,size(xmax)); +imin = reshape(imin,size(xmin)); + +% Descending order: +[temp,inmax] = sort(-xmax); clear temp +xmax = xmax(inmax); +imax = imax(inmax); +[xmin,inmin] = sort(xmin); +imin = imin(inmin); + + +% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com \ No newline at end of file diff --git a/extrema2.m b/extrema2.m new file mode 100644 index 0000000..53a5e76 --- /dev/null +++ b/extrema2.m @@ -0,0 +1,174 @@ +function [xymax,smax,xymin,smin] = extrema2(xy,varargin) +%EXTREMA2 Gets the extrema points from a surface. +% [XMAX,IMAX,XMIN,IMIN] = EXTREMA2(X) returns the maxima and minima +% elements of the matriz X ignoring NaN's, where +% XMAX - maxima points in descending order (the bigger first and so on) +% IMAX - linear indexes of the XMAX +% XMIN - minima points in descending order +% IMIN - linear indexes of the XMIN. +% The program uses EXTREMA. +% +% The extrema points are searched only through the column, the row and +% the diagonals crossing each matrix element, so it is not a perfect +% mathematical program and for this reason it has an optional argument. +% The user should be aware of these limitations. +% +% [XMAX,IMAX,XMIN,IMIN] = EXTREMA2(X,1) does the same but without +% searching through the diagonals (less strict and perhaps the user gets +% more output points). +% +% DEFINITION (from http://en.wikipedia.org/wiki/Maxima_and_minima): +% In mathematics, maxima and minima, also known as extrema, are points in +% the domain of a function at which the function takes a largest value +% (maximum) or smallest value (minimum), either within a given +% neighbourhood (local extrema) or on the function domain in its entirety +% (global extrema). +% +% Note: To change the linear index to (i,j) use IND2SUB. +% +% Example: +% [x,y] = meshgrid(-2:.2:2,3:-.2:-2); +% z = x.*exp(-x.^2-y.^2); z(10,7)= NaN; z(16:19,13:17) = NaN; +% surf(x,y,z), shading interp +% [zmax,imax,zmin,imin] = extrema2(z); +% hold on +% plot3(x(imax),y(imax),zmax,'bo',x(imin),y(imin),zmin,'ro') +% for i = 1:length(zmax) +% text(x(imax(i)),y(imax(i)),zmax(i),[' ' num2str(zmax(i))]) +% end +% for i = 1:length(zmin) +% text(x(imin(i)),y(imin(i)),zmin(i),[' ' num2str(zmin(i))]) +% end +% hold off +% +% See also EXTREMA, MAX, MIN + +% Written by +% Lic. on Physics Carlos Adrián Vargas Aguilera +% Physical Oceanography MS candidate +% UNIVERSIDAD DE GUADALAJARA +% Mexico, 2005 +% +% nubeobscura@hotmail.com + +% From : http://www.mathworks.com/matlabcentral/fileexchange +% File ID : 12275 +% Submited at: 2006-09-14 +% 2006-11-11 : English translation from spanish. +% 2006-11-17 : Accept NaN's. +% 2006-11-22 : Fixed bug in INDX (by JaeKyu Suhr) +% 2007-04-09 : Change name to MAXIMA2, and definition added. + +M = size(xy); +if length(M) ~= 2 + error('Entry must be a matrix.') +end +N = M(2); +M = M(1); + +% Search peaks through columns: +[smaxcol,smincol] = extremos(xy); + +% Search peaks through rows, on columns with extrema points: +im = unique([smaxcol(:,1);smincol(:,1)]); % Rows with column extrema +[smaxfil,sminfil] = extremos(xy(im,:).'); + +% Convertion from 2 to 1 index: +smaxcol = sub2ind([M,N],smaxcol(:,1),smaxcol(:,2)); +smincol = sub2ind([M,N],smincol(:,1),smincol(:,2)); +smaxfil = sub2ind([M,N],im(smaxfil(:,2)),smaxfil(:,1)); +sminfil = sub2ind([M,N],im(sminfil(:,2)),sminfil(:,1)); + +% Peaks in rows and in columns: +smax = intersect(smaxcol,smaxfil); +smin = intersect(smincol,sminfil); + +% Search peaks through diagonals? +if nargin==1 + % Check peaks on down-up diagonal: + [iext,jext] = ind2sub([M,N],unique([smax;smin])); + [sextmax,sextmin] = extremos_diag(iext,jext,xy,1); + + % Check peaks on up-down diagonal: + smax = intersect(smax,[M; (N*M-M); sextmax]); + smin = intersect(smin,[M; (N*M-M); sextmin]); + + % Peaks on up-down diagonals: + [iext,jext] = ind2sub([M,N],unique([smax;smin])); + [sextmax,sextmin] = extremos_diag(iext,jext,xy,-1); + + % Peaks on columns, rows and diagonals: + smax = intersect(smax,[1; N*M; sextmax]); + smin = intersect(smin,[1; N*M; sextmin]); +end + +% Extrema points: +xymax = xy(smax); +xymin = xy(smin); + +% Descending order: +[temp,inmax] = sort(-xymax); clear temp +xymax = xymax(inmax); +smax = smax(inmax); +[xymin,inmin] = sort(xymin); +smin = smin(inmin); + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +function [smax,smin] = extremos(matriz) +% Peaks through columns or rows. + +smax = []; +smin = []; + +for n = 1:length(matriz(1,:)) + [temp,imaxfil,temp,iminfil] = extrema(matriz(:,n)); clear temp + if ~isempty(imaxfil) % Maxima indexes + imaxcol = repmat(n,length(imaxfil),1); + smax = [smax; imaxfil imaxcol]; + end + if ~isempty(iminfil) % Minima indexes + imincol = repmat(n,length(iminfil),1); + smin = [smin; iminfil imincol]; + end +end + + +function [sextmax,sextmin] = extremos_diag(iext,jext,xy,A) +% Peaks through diagonals (down-up A=-1) + +[M,N] = size(xy); +if A==-1 + iext = M-iext+1; +end +[iini,jini] = cruce(iext,jext,1,1); +[iini,jini] = ind2sub([M,N],unique(sub2ind([M,N],iini,jini))); +[ifin,jfin] = cruce(iini,jini,M,N); +sextmax = []; +sextmin = []; +for n = 1:length(iini) + ises = iini(n):ifin(n); + jses = jini(n):jfin(n); + if A==-1 + ises = M-ises+1; + end + s = sub2ind([M,N],ises,jses); + [temp,imax,temp,imin] = extrema(xy(s)); clear temp + sextmax = [sextmax; s(imax)']; + sextmin = [sextmin; s(imin)']; +end + + +function [i,j] = cruce(i0,j0,I,J) +% Indexes where the diagonal of the element io,jo crosses the left/superior +% (I=1,J=1) or right/inferior (I=M,J=N) side of an MxN matrix. + +arriba = 2*(I*J==1)-1; + +si = (arriba*(j0-J) > arriba*(i0-I)); +i = (I - (J+i0-j0)).*si + J+i0-j0; +j = (I+j0-i0-(J)).*si + J; + + +% Carlos Adrián Vargas Aguilera. nubeobscura@hotmail.com \ No newline at end of file diff --git a/getPeak.m b/getPeak.m new file mode 100644 index 0000000..088a9ca --- /dev/null +++ b/getPeak.m @@ -0,0 +1,42 @@ +function [fwhm, absorption_min] = getPeak(detunings, absorptions) + + eit_data = [detunings, absorptions]; + + %Get local extrema values and indexes and sort them + [maxima, max_index, minima, min_index] = extrema(eit_data(:,2)); + minima = [min_index, minima]; + maxima = [max_index, maxima]; + + peaks = sortrows([maxima; minima], 1); + + % Find the set of extrema closest to being centered on 0 + center_index = find(eit_data(:,1) == 0); + differences = abs(peaks(:,1) - center_index); + [detuning_index_offset, peak_center_index] = min(differences); + peak_center = peaks(peak_center_index, :); + left_max = peaks(peak_center_index - 1, :); + right_max = peaks(peak_center_index + 1, :); + + %Define each point of interest + absorption_min = peak_center(2); + left_range = left_max(1):peak_center(1); + right_range = peak_center(1):right_max(1); + + %Calculate the width + absorption_max = min([right_max(2), left_max(2)]); % use the smallest of the maximums + absorption_mid = (absorption_max + absorption_min) / 2; + half_width_right = interp1(eit_data(right_range,2), eit_data(right_range,1), absorption_mid) - eit_data(peak_center(1),1); + half_width_left = eit_data(peak_center(1),1) - interp1(eit_data(left_range,2), eit_data(left_range,1), absorption_mid); + fwhm = half_width_right + half_width_left; + + + %Plot extrema over existing plot + % plot(eit_data(:,1), eit_data(:,2)); + % hold all + % grid on + % plot(detuning_1(left_max(1)), left_max(2), 'r*') + % plot(detuning_1(right_max(1)), right_max(2), 'r*') + % plot(detuning_1(peak_center(1)), peak_center(2), 'g*') + + +end diff --git a/loadSimulation.m b/loadSimulation.m new file mode 100644 index 0000000..e5f0be3 --- /dev/null +++ b/loadSimulation.m @@ -0,0 +1,7 @@ +function [width, minAbsorption, detunings, absorptions] = loadSimulation(dataFile) + data = load(dataFile); + detunings = data(:,1); + absorptions = data(:,2); + [width, minAbsorption] = getPeak(detunings, absorptions); + plot(detunings, absorptions); +end \ No newline at end of file diff --git a/loadSimulations.m b/loadSimulations.m new file mode 100644 index 0000000..046e336 --- /dev/null +++ b/loadSimulations.m @@ -0,0 +1,27 @@ +function [data] = loadSimulations(decay_bc, dephase_bc) + dataFiles = dir(['data/', decay_bc, '*_*_*_', dephase_bc, '*.dat']); + widths = []; + minAbsorptions = []; + drives = []; + deltas = []; + for i = 1:length(dataFiles) + dataFile = ['data/', dataFiles(i).name]; + data = load(dataFile); + detunings = data(:,1); + absorptions = data(:,2); + [width, absorption] = getPeak(detunings, absorptions); + widths = [widths; width]; + minAbsorptions = [minAbsorptions; absorption]; + sections = strfind(dataFile, '_'); + drive = dataFile(sections(1)+1:sections(2)-1); + drives = [drives; str2num(drive)]; + delta = dataFile(sections(2)+1:sections(3)-1); + deltas = [deltas; str2num(delta)]; + end + data = [drives, widths, minAbsorptions]; + data = sortrows(data, 1); + drives = data(:,1); widths = data(:,2); minAbsorptions = data(:,3); + loglog(drives.^2, widths); + hold all + loglog(drives.^2, minAbsorptions*10^5); +end \ No newline at end of file diff --git a/simulate.m b/simulate.m new file mode 100644 index 0000000..ebf882f --- /dev/null +++ b/simulate.m @@ -0,0 +1,74 @@ +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.001; +decay_ab = 6; + +widths = []; +absorptions = []; +drives = logspace(-2, -3, 10); +drives = linspace(0.003298,0.003319,20); +detunings1 = [0:2:0]; +dephase_bc = 0.05; + +for drive = drives + width_row = []; + absorption_row = []; + expected_width = (drive^2/decay_ab + dephase_bc + decay_bc); + domain_max = expected_width*10; + domain_min = -1*domain_max; + final_time = 1/expected_width*100; + for detuning1 = detunings1 + + run_eit = sprintf('time ./eit --drive_arg=%f --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, dephase_bc, 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/%f_%f_%f_%f.dat', decay_bc, drive, detuning1, dephase_bc); + csvwrite(csv_name, eit_plot); + 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']); -- cgit v1.2.3