diff options
| -rw-r--r-- | eit.xmds | 136 | ||||
| -rw-r--r-- | extrema.m | 146 | ||||
| -rw-r--r-- | extrema2.m | 174 | ||||
| -rw-r--r-- | getPeak.m | 42 | ||||
| -rw-r--r-- | loadSimulation.m | 7 | ||||
| -rw-r--r-- | loadSimulations.m | 27 | ||||
| -rw-r--r-- | simulate.m | 74 |
7 files changed, 606 insertions, 0 deletions
diff --git a/eit.xmds b/eit.xmds new file mode 100644 index 0000000..50419f5 --- /dev/null +++ b/eit.xmds @@ -0,0 +1,136 @@ +<?xml version="1.0" encoding="UTF-8"?> +<simulation xmds-version="2"> + + <!-- While not strictly necessary, the following two tags are handy. --> + <author>Hunter Rew</author> + <description> + Model of 3 state atom + </description> + + <!-- + This element defines some constants. It can be used for other + features as well, but we will go into that in more detail later. + --> + <features> + <precision> double </precision> + <globals> + <![CDATA[ + real drive; + real decay; + real decay_bc; + complex r_ca; + real probe = .0001; + real delta1; + real dephase_bc; + + ]]> + </globals> + <validation kind="run-time"/> + <arguments> + <argument name="drive_arg" type="real" default_value="0.1" /> + <argument name="decay_arg" type="real" default_value="10.0" /> + <argument name="decay_bc_arg" type="real" default_value="0.001" /> + <argument name="domain_min" type="real" default_value="-1" /> + <argument name="domain_max" type="real" default_value="1" /> + <argument name="delta1_arg" type="real" default_value="0" /> + <argument name="dephase_bc_arg" type="real" default_value="0" /> + <argument name="final_time" type="real" default_value="10000" /> + <![CDATA[ + drive = drive_arg; + decay = decay_arg; + decay_bc = decay_bc_arg; + delta1 = delta1_arg; + dephase_bc = dephase_bc_arg; + r_ca = decay - i*delta1; + ]]> + </arguments> + </features> + + <!-- + This part defines all of the dimensions used in the problem, + in this case, only the dimension of 'time' is needed. + --> + <geometry> + <propagation_dimension> t </propagation_dimension> + <transverse_dimensions> + <dimension name="detuning" lattice="256" domain="(domain_min, domain_max)" /> + </transverse_dimensions> + </geometry> + + <!-- A 'vector' describes the variables that we will be evolving. --> + <vector name="Gamma" type="complex"> + <components> r_ab r_cb </components> + <initialisation> + <![CDATA[ + r_ab = decay + i*(delta1 + detuning); + r_cb = decay_bc + i*detuning + dephase_bc; + ]]> + </initialisation> + </vector> + <vector name="population" type="complex" dimensions="detuning"> + <components> + Paa Pbb Pcc Pab Pca Pcb + </components> + <initialisation> + <![CDATA[ + Pcc = .5; + Pbb = .5; + Paa = Pab = Pca = Pcb = 0; + ]]> + </initialisation> + </vector> + + <sequence> + <!-- + Here we define what differential equations need to be solved + and what algorithm we want to use. + --> + <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> + <samples>10</samples> + <operators> + <integration_vectors>population</integration_vectors> + <![CDATA[ + dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; + + dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; + + dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; + + dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; + + dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; + + dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; + ]]> + <dependencies>Gamma</dependencies> + </operators> + </integrate> + </sequence> + + <!-- This part defines what data will be saved in the output file --> + <output format="hdf5" > + <sampling_group basis="detuning" initial_sample="yes"> + <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> + <dependencies>population</dependencies> + <![CDATA[ + PaaR = Paa.Re(); + PaaI = Paa.Im(); + PbbR = Pbb.Re(); + PbbI = Pbb.Im(); + PccR = Pcc.Re(); + PccI = Pcc.Im(); + PabR = Pab.Re(); + PabI = Pab.Im(); + PcaR = Pca.Re(); + PcaI = Pca.Im(); + PcbR = Pcb.Re(); + PcbI = Pcb.Im(); + ]]> --> + <moments>PabI</moments> + <dependencies>population</dependencies> + <![CDATA[ + PabI = Pab.Im(); + ]]> + </sampling_group> + </output> +</simulation> 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']); |