Fixed axis limits

8 files changed, 103 insertions, 31 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..8f0ba77
--- /dev/null
+++ b/README
@@ -0,0 +1,53 @@
+Authors: 	
+	Hunter Rew (hbrew@email.wm.edu), 
+	Eugeniy Mikhailov (eemikh@wm.edu)
+
+Description: 
+	This program is a set of scripts for simulating 
+	and plotting atom-light configurations which produce electromagnetically 
+	induced transparency.
+
+Dependencies:
+	xmds2 and all of its dependencies, including the optional dependency numpy,
+	Octave
+
+Usage:
+	Set simulation parameters in config.m
+	Open Octave in the eit_filter_simulations directory
+	Run config: 
+	# config
+	Run simulate:
+	# simulate
+
+Data files are stored by default in the data directory.
+The naming convention is: 
+	(ground state decay)_(drive intensity)_(1-photon detuning)_(ground state dephasing).h5
+
+Data loading and plotting examples:
+	Plot 2-photon detuning vs absorption for a particular simulation:
+		# [detunings2, absorptions] = getSimulation('data/data_file_name.h5');
+		# figure;
+		# plotSimulation(detunings2, absorptions);
+
+	Plot width and transmission contrast vs drive intensity for given parameters:
+		# [drives, widths, min_absorptions, contrasts] = getSimulations('ground state decay', '1-photon detuning', 'ground state dephasing');
+		note: 
+			you only need to specify these values to the decimal place that the file names change.
+			if you have:
+				0.001000_(drive)_1.000000_0.005000.h5
+				0.001000_(drive)_2.000000_0.010000.h5
+				0.001000_(drive)_2.000000_0.000000.h5
+				0.000100_(drive)_1.000000_0.000000.h5
+			and you want all data with:
+				ground state decay = 0.001,
+				1-photon detuning = 2,
+				ground state dephasing = 0
+			you would call:
+				# getSimulations('0.0010', '2', '0.000')
+		# figure;
+		# plotSimulations(drives, widths, contrasts);
+		
+	Plot width and absorption vs drive intensity for given parameters:
+		Same as above except for the last line.
+		# plotSimulations(drives, widths, min_absorptions) 
+
diff --git a/config.m b/config.m
index 264b8ee..2592ce6 100644
--- a/config.m
+++ b/config.m
@@ -5,18 +5,20 @@
 
 decay_bc = [0.00001];	% Ground state decay
 decay_ab = [6];	% Excited state decay
-drive = [logspace(0,-2,100)];	% Drive/control field
-drive = drive(1:1); % Select which drives you want here
-detuning1 = [-15];	% Single photon detuning
+drive = [logspace(-2,-3,4)];	% Drive/control field
+drive = drive(2:4); % Select which drives you want here
+detuning1 = [0];	% Single photon detuning
 dephase_bc = [0];	% Ground state dephasing
 
-global n = 3;	%Number of processes to run
+global n;	%Number of processes to run
+n = 3;
 
 format long
 
-global source_dir = 'source/';
-global config_dir = [source_dir, 'config/'];
-global data_dir = 'data/'
+global source_dir; global config_dir; global data_dir;
+source_dir = 'source/';
+config_dir = [source_dir, 'config/'];
+data_dir = 'data/';
 
 % Write parameters to files
 dlmwrite([config_dir, 'decay_bc'], decay_bc, 'precision', '%0.10f');
diff --git a/control/getSimulation.m b/control/getSimulation.m
index bbc8b45..2cf09c6 100644
--- a/control/getSimulation.m
+++ b/control/getSimulation.m
@@ -1,7 +1,7 @@
 function [detunings2, absorptions] = getSimulation(dataFile)
 
 	data = load(dataFile);
-	detunings2 = data(:,1);
-	absorptions = data(:,2);
+	detunings2 = data._1.detuning';
+	absorptions = data._1.PabI(:,end);
 
 end
\ No newline at end of file
diff --git a/control/getSimulations.m b/control/getSimulations.m
index c57d690..c072a1f 100644
--- a/control/getSimulations.m
+++ b/control/getSimulations.m
@@ -1,6 +1,7 @@
-function [widths, minAbsorptions, drives, contrast] = getSimulations(directory, decay_bc, detuning, dephase_bc)
+function [drives, widths, minAbsorptions, contrasts] = getSimulations(decay_bc, detuning, dephase_bc)
 
-	dataFiles = dir([directory, decay_bc, '*_*_', detuning, '*_', dephase_bc, '*.dat']);
+	global data_dir
+	dataFiles = dir([data_dir, decay_bc, '*_*_', detuning, '*_', dephase_bc, '*.h5']);
 	widths = [];
 	minAbsorptions = [];
 	contrasts = [];
@@ -10,7 +11,7 @@ function [widths, minAbsorptions, drives, contrast] = getSimulations(directory,
 	for i = 1:length(dataFiles)
 
 		dataFile = dataFiles(i).name;
-		[detunings2, absorptions] = getSimulation([directory, dataFile]);
+		[detunings2, absorptions] = getSimulation([data_dir, dataFile]);
 		[width, absorption, contrast] = getPeak(detunings2, absorptions);
 		widths = [widths; width];
 		minAbsorptions = [minAbsorptions; absorption];
@@ -23,6 +24,7 @@ function [widths, minAbsorptions, drives, contrast] = getSimulations(directory,
 
 	end
 
+	% Sort by drive
 	data = [drives, widths, minAbsorptions, contrasts];
 	data = sortrows(data, 1);
 	drives = data(:,1); widths = data(:,2); minAbsorptions = data(:,3); contrasts = data(:,4);
diff --git a/simulate.m b/simulate.m
index aeef180..35e6ca0 100644
--- a/simulate.m
+++ b/simulate.m
@@ -5,7 +5,7 @@ function simulate()
 
 	target = [source_dir, 'simulate.py'];
 	tic();
-	runCommand(sprintf('python %s %s %d', target, source_dir, n));
+	runCommand(sprintf('python2 %s %s %d', target, source_dir, n));
 	toc();
 
 end
diff --git a/source/simulate.py b/source/simulate.py
index d0ff380..2e35c38 100644
--- a/source/simulate.py
+++ b/source/simulate.py
@@ -29,7 +29,7 @@ def splitRange(x, n):
 # sort every other array in x in reverse order
 def mixRange(x, n):
 	i = len(x) % 2
-	for j in range(i,len(x),2):
+	for j in range((1-i),len(x),2):
 		x[j].sort(reverse=True)
 	newX = []
 	for i in range(0,len(x[0])):
@@ -64,7 +64,7 @@ def compile(target):
 
 def run(target, options):
 	command = target + ' ' + options
-	#print command + '\n'
+	print command + '\n'
 	bash(command)
 
 def xsil2graphics2(target):
@@ -73,13 +73,17 @@ def xsil2graphics2(target):
 	bash(command)
 
 def load(target):
-	data = __import__(target)
+	data = __import__(target + '.py')
 	data = [data.detuning_1.tolist(), data.PabI_1[-1,:].tolist()]
 	return np.transpose(data)
 
 def savecsv(target, data):
 	np.savetxt(target, data, delimiter=',')
 
+def saveh5(target, destination):
+	command = 'cp %s.h5 %s' % (target, destination)
+	run(command, '')
+
 def simulate((index, decays_ab, decays_bc, drives, detunings1, dephasings_bc)):
 	print 'Simulation batch ' + str(index) + ' started.' + '\n'
 	source = 'eit' + str(index)
@@ -91,20 +95,21 @@ def simulate((index, decays_ab, decays_bc, drives, detunings1, dephasings_bc)):
 			for drive in drives:
 				for detuning1 in detunings1:
 					for dephase_bc in dephasings_bc:
+						
 						expected_width = (math.pow(drive,2)/decay_ab + dephase_bc + decay_bc)
 						domain_max = expected_width*10
 						domain_min = -1*domain_max
 						final_time = 1/expected_width*time_multiplier
+						
 						options = '--drive_arg=%0.20f --decay_arg=%f --decay_bc_arg=%f --delta1_arg=%f --dephase_bc_arg=%f --domain_min=%f --domain_max=%f --final_time=%f'
 						values = (drive, decay_ab, decay_bc, detuning1, dephase_bc, domain_min, domain_max, final_time)
 						options = options % values
 						run(target, options)
-						xsil2graphics2(target + '.xsil')
-						data = load(source)
-						dataFile = '%s%0.10f_%0.10f_%0.10f_%0.10f.dat' % (data_dir, decay_bc, drive, detuning1, dephase_bc)
-						savecsv(dataFile, data)
-
+						
+						dataFile = '%s%0.10f_%0.10f_%0.10f_%0.10f.h5' % (data_dir, decay_bc, drive, detuning1, dephase_bc)
+						saveh5(target, dataFile)
 	
+
 	print 'Simulation batch ' + str(index) + ' finished.' + '\n'
 
 def main():
@@ -125,20 +130,20 @@ def main():
 	config_dir = source_dir + 'config/'
 
 	global time_multiplier
-	time_multiplier = 100
+	time_multiplier = 10 # This should be about 100 to be certain the system has settled for the smallest drives.
 	
 	n = int(sys.argv[2]) # Number of processes to run
 
 	decays_ab = readArgs('decay_ab')
 	decays_bc = readArgs('decay_bc')
 	drives = readArgs('drive')
-	drives = splitRange(drives, n)
-	drives = mixRange(drives, n)
-	drives = splitRange(drives, n)
+	drives = splitRange(drives, n)	# Split drives into n arrays
+	drives = mixRange(drives, n)	# Shuffle the n arrays back into a single array. This is so that each process will receive a mix of large and small drive values.
+	drives = splitRange(drives, n)	# Split drives back into n arrays
 	detunings1 = readArgs('detuning1')
 	dephasings_bc = readArgs('dephase_bc')
 
-	prepSource(n)
+	prepSource(n)	# Make copies of the xmds script for each process
 	options = []
 	for i in range(0, len(drives)):
 		drive = drives[i]
diff --git a/view/plotSimulation.m b/view/plotSimulation.m
index d494287..a741333 100644
--- a/view/plotSimulation.m
+++ b/view/plotSimulation.m
@@ -1,5 +1,8 @@
 function h = plotSimulation(detunings2, absorptions)
 
 	h = plot(detunings2, absorptions);
+	xlabel('2-photon Detuning')
+	ylabel('Absorption')
+	title('Absorption Vs 2-photon Detuning')
 
 end
\ No newline at end of file
diff --git a/view/plotSimulations.m b/view/plotSimulations.m
index 12273af..6197c75 100644
--- a/view/plotSimulations.m
+++ b/view/plotSimulations.m
@@ -1,11 +1,18 @@
-function plotSimulations(widths, absorptions, drives)
+function plotSimulations(drives, y1, y2)
 
-	loglog(drives, widths)
+	loglog(drives, y1)
 	hold all
-	loglog(drives, absorptions)
+	loglog(drives, y2)
 	xlabel('Drive Intensity (log scale)')
 	ylabel('(log scale)')
-	title('Width and Absorption Vs Drive Intensity')
-	legend('Width', 'Absorption')
+	title(sprintf('%s and %s Vs Drive Intensity', inputname(2), inputname(3)))
+	legend(inputname(2), inputname(3))
+
+	% When plotting contrast, the y-axis upper limit is the maximum value (1)
+	% Since contrast approaches one quickly and becomes a practically straight line, it's hard to see
+	% Here we increase the y-axis upper limit to prevent the contrast line from overlapping with the axis
+	limits = ylim;
+	ylim([limits(1), limits(2)*10]);
+	hold off
 
 end
\ No newline at end of file