Hunter's initial release

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']);