diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..6ed4c25
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+2DMOT Simulation Code/Results
+AtomECS Simulation Code
+TimeSeriesAnalyzer/Time Series Data
+ULE Cavity Characterisitics/Data
+ULE Cavity Characterisitics/Figures
\ No newline at end of file
diff --git a/+Helper/ImageSelection.class b/2DMOT Simulation Code/+Helper/ImageSelection.class
similarity index 100%
rename from +Helper/ImageSelection.class
rename to 2DMOT Simulation Code/+Helper/ImageSelection.class
diff --git a/+Helper/ImageSelection.java b/2DMOT Simulation Code/+Helper/ImageSelection.java
similarity index 100%
rename from +Helper/ImageSelection.java
rename to 2DMOT Simulation Code/+Helper/ImageSelection.java
diff --git a/+Helper/PhysicsConstants.m b/2DMOT Simulation Code/+Helper/PhysicsConstants.m
similarity index 100%
rename from +Helper/PhysicsConstants.m
rename to 2DMOT Simulation Code/+Helper/PhysicsConstants.m
diff --git a/+Helper/bringFiguresWithTagInForeground.m b/2DMOT Simulation Code/+Helper/bringFiguresWithTagInForeground.m
similarity index 100%
rename from +Helper/bringFiguresWithTagInForeground.m
rename to 2DMOT Simulation Code/+Helper/bringFiguresWithTagInForeground.m
diff --git a/+Helper/calculateDistanceFromPointToLine.m b/2DMOT Simulation Code/+Helper/calculateDistanceFromPointToLine.m
similarity index 100%
rename from +Helper/calculateDistanceFromPointToLine.m
rename to 2DMOT Simulation Code/+Helper/calculateDistanceFromPointToLine.m
diff --git a/+Helper/convertstruct2cell.m b/2DMOT Simulation Code/+Helper/convertstruct2cell.m
similarity index 100%
rename from +Helper/convertstruct2cell.m
rename to 2DMOT Simulation Code/+Helper/convertstruct2cell.m
diff --git a/+Helper/findAllZeroCrossings.m b/2DMOT Simulation Code/+Helper/findAllZeroCrossings.m
similarity index 100%
rename from +Helper/findAllZeroCrossings.m
rename to 2DMOT Simulation Code/+Helper/findAllZeroCrossings.m
diff --git a/+Helper/getFigureByTag.m b/2DMOT Simulation Code/+Helper/getFigureByTag.m
similarity index 100%
rename from +Helper/getFigureByTag.m
rename to 2DMOT Simulation Code/+Helper/getFigureByTag.m
diff --git a/+Helper/ode5.m b/2DMOT Simulation Code/+Helper/ode5.m
similarity index 100%
rename from +Helper/ode5.m
rename to 2DMOT Simulation Code/+Helper/ode5.m
diff --git a/+Helper/onenoteccdata.m b/2DMOT Simulation Code/+Helper/onenoteccdata.m
similarity index 100%
rename from +Helper/onenoteccdata.m
rename to 2DMOT Simulation Code/+Helper/onenoteccdata.m
diff --git a/+Helper/parforNotifications.m b/2DMOT Simulation Code/+Helper/parforNotifications.m
similarity index 100%
rename from +Helper/parforNotifications.m
rename to 2DMOT Simulation Code/+Helper/parforNotifications.m
diff --git a/+Helper/screencapture.m b/2DMOT Simulation Code/+Helper/screencapture.m
similarity index 100%
rename from +Helper/screencapture.m
rename to 2DMOT Simulation Code/+Helper/screencapture.m
diff --git a/+Plotter/plotAngularDistributionForDifferentBeta.m b/2DMOT Simulation Code/+Plotter/plotAngularDistributionForDifferentBeta.m
similarity index 100%
rename from +Plotter/plotAngularDistributionForDifferentBeta.m
rename to 2DMOT Simulation Code/+Plotter/plotAngularDistributionForDifferentBeta.m
diff --git a/+Plotter/plotCaptureVelocityVsAngle.m b/2DMOT Simulation Code/+Plotter/plotCaptureVelocityVsAngle.m
similarity index 100%
rename from +Plotter/plotCaptureVelocityVsAngle.m
rename to 2DMOT Simulation Code/+Plotter/plotCaptureVelocityVsAngle.m
diff --git a/+Plotter/plotConfidenceIntervalRegion.m b/2DMOT Simulation Code/+Plotter/plotConfidenceIntervalRegion.m
similarity index 100%
rename from +Plotter/plotConfidenceIntervalRegion.m
rename to 2DMOT Simulation Code/+Plotter/plotConfidenceIntervalRegion.m
diff --git a/+Plotter/plotDynamicalQuantities.m b/2DMOT Simulation Code/+Plotter/plotDynamicalQuantities.m
similarity index 100%
rename from +Plotter/plotDynamicalQuantities.m
rename to 2DMOT Simulation Code/+Plotter/plotDynamicalQuantities.m
diff --git a/+Plotter/plotFreeMolecularFluxVsTemp.m b/2DMOT Simulation Code/+Plotter/plotFreeMolecularFluxVsTemp.m
similarity index 100%
rename from +Plotter/plotFreeMolecularFluxVsTemp.m
rename to 2DMOT Simulation Code/+Plotter/plotFreeMolecularFluxVsTemp.m
diff --git a/+Plotter/plotInitialVeloctiySamplingVsAngle.m b/2DMOT Simulation Code/+Plotter/plotInitialVeloctiySamplingVsAngle.m
similarity index 100%
rename from +Plotter/plotInitialVeloctiySamplingVsAngle.m
rename to 2DMOT Simulation Code/+Plotter/plotInitialVeloctiySamplingVsAngle.m
diff --git a/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m b/2DMOT Simulation Code/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m
similarity index 100%
rename from +Plotter/plotMeanFreePathAndVapourPressureVsTemp.m
rename to 2DMOT Simulation Code/+Plotter/plotMeanFreePathAndVapourPressureVsTemp.m
diff --git a/+Plotter/plotPhaseSpaceWithAccelerationField.m b/2DMOT Simulation Code/+Plotter/plotPhaseSpaceWithAccelerationField.m
similarity index 100%
rename from +Plotter/plotPhaseSpaceWithAccelerationField.m
rename to 2DMOT Simulation Code/+Plotter/plotPhaseSpaceWithAccelerationField.m
diff --git a/+Plotter/plotPositionAndVelocitySampling.m b/2DMOT Simulation Code/+Plotter/plotPositionAndVelocitySampling.m
similarity index 100%
rename from +Plotter/plotPositionAndVelocitySampling.m
rename to 2DMOT Simulation Code/+Plotter/plotPositionAndVelocitySampling.m
diff --git a/+Plotter/plotResultForOneParameterScan.m b/2DMOT Simulation Code/+Plotter/plotResultForOneParameterScan.m
similarity index 100%
rename from +Plotter/plotResultForOneParameterScan.m
rename to 2DMOT Simulation Code/+Plotter/plotResultForOneParameterScan.m
diff --git a/+Plotter/plotResultForThreeParameterScan.m b/2DMOT Simulation Code/+Plotter/plotResultForThreeParameterScan.m
similarity index 100%
rename from +Plotter/plotResultForThreeParameterScan.m
rename to 2DMOT Simulation Code/+Plotter/plotResultForThreeParameterScan.m
diff --git a/+Plotter/plotResultForTwoParameterScan.m b/2DMOT Simulation Code/+Plotter/plotResultForTwoParameterScan.m
similarity index 100%
rename from +Plotter/plotResultForTwoParameterScan.m
rename to 2DMOT Simulation Code/+Plotter/plotResultForTwoParameterScan.m
diff --git a/+Plotter/visualizeMagneticField.m b/2DMOT Simulation Code/+Plotter/visualizeMagneticField.m
similarity index 100%
rename from +Plotter/visualizeMagneticField.m
rename to 2DMOT Simulation Code/+Plotter/visualizeMagneticField.m
diff --git a/+Scripts/optimizingForSidebandEnhancement.m b/2DMOT Simulation Code/+Scripts/optimizingForSidebandEnhancement.m
similarity index 100%
rename from +Scripts/optimizingForSidebandEnhancement.m
rename to 2DMOT Simulation Code/+Scripts/optimizingForSidebandEnhancement.m
diff --git a/+Scripts/scanForSidebandEnhancement.m b/2DMOT Simulation Code/+Scripts/scanForSidebandEnhancement.m
similarity index 100%
rename from +Scripts/scanForSidebandEnhancement.m
rename to 2DMOT Simulation Code/+Scripts/scanForSidebandEnhancement.m
diff --git a/+Simulator/+Scan/doOneParameter.m b/2DMOT Simulation Code/+Simulator/+Scan/doOneParameter.m
similarity index 100%
rename from +Simulator/+Scan/doOneParameter.m
rename to 2DMOT Simulation Code/+Simulator/+Scan/doOneParameter.m
diff --git a/+Simulator/+Scan/doThreeParameters.m b/2DMOT Simulation Code/+Simulator/+Scan/doThreeParameters.m
similarity index 100%
rename from +Simulator/+Scan/doThreeParameters.m
rename to 2DMOT Simulation Code/+Simulator/+Scan/doThreeParameters.m
diff --git a/+Simulator/+Scan/doTwoParameters.m b/2DMOT Simulation Code/+Simulator/+Scan/doTwoParameters.m
similarity index 100%
rename from +Simulator/+Scan/doTwoParameters.m
rename to 2DMOT Simulation Code/+Simulator/+Scan/doTwoParameters.m
diff --git a/+Simulator/@Beams/Beams.m b/2DMOT Simulation Code/+Simulator/@Beams/Beams.m
similarity index 100%
rename from +Simulator/@Beams/Beams.m
rename to 2DMOT Simulation Code/+Simulator/@Beams/Beams.m
diff --git a/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m b/2DMOT Simulation Code/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m
similarity index 100%
rename from +Simulator/@MOTCaptureProcess/MOTCaptureProcess.m
rename to 2DMOT Simulation Code/+Simulator/@MOTCaptureProcess/MOTCaptureProcess.m
diff --git a/+Simulator/@Oven/Oven.m b/2DMOT Simulation Code/+Simulator/@Oven/Oven.m
similarity index 100%
rename from +Simulator/@Oven/Oven.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/Oven.m
diff --git a/+Simulator/@Oven/angularDistributionFunction.m b/2DMOT Simulation Code/+Simulator/@Oven/angularDistributionFunction.m
similarity index 100%
rename from +Simulator/@Oven/angularDistributionFunction.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/angularDistributionFunction.m
diff --git a/+Simulator/@Oven/calculateClausingFactor.m b/2DMOT Simulation Code/+Simulator/@Oven/calculateClausingFactor.m
similarity index 100%
rename from +Simulator/@Oven/calculateClausingFactor.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/calculateClausingFactor.m
diff --git a/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m b/2DMOT Simulation Code/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m
similarity index 100%
rename from +Simulator/@Oven/calculateFreeMolecularRegimeFlux.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/calculateFreeMolecularRegimeFlux.m
diff --git a/+Simulator/@Oven/calculateReducedClausingFactor.m b/2DMOT Simulation Code/+Simulator/@Oven/calculateReducedClausingFactor.m
similarity index 100%
rename from +Simulator/@Oven/calculateReducedClausingFactor.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/calculateReducedClausingFactor.m
diff --git a/+Simulator/@Oven/initialPositionSampling.m b/2DMOT Simulation Code/+Simulator/@Oven/initialPositionSampling.m
similarity index 100%
rename from +Simulator/@Oven/initialPositionSampling.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/initialPositionSampling.m
diff --git a/+Simulator/@Oven/initialVelocitySampling.m b/2DMOT Simulation Code/+Simulator/@Oven/initialVelocitySampling.m
similarity index 100%
rename from +Simulator/@Oven/initialVelocitySampling.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/initialVelocitySampling.m
diff --git a/+Simulator/@Oven/velocityDistributionFunction.m b/2DMOT Simulation Code/+Simulator/@Oven/velocityDistributionFunction.m
similarity index 100%
rename from +Simulator/@Oven/velocityDistributionFunction.m
rename to 2DMOT Simulation Code/+Simulator/@Oven/velocityDistributionFunction.m
diff --git a/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/TwoDimensionalMOT.m
diff --git a/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/accelerationDueToPushBeam.m
diff --git a/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/accelerationDueToSpontaneousEmissionProcess.m
diff --git a/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/bootstrapErrorEstimation.m
diff --git a/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateCaptureVelocity.m
diff --git a/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateLoadingRate.m
diff --git a/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateLocalSaturationIntensity.m
diff --git a/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/calculateTotalAcceleration.m
diff --git a/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/computeCollisionProbability.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/computeCollisionProbability.m
diff --git a/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/computeTimeSpentInInteractionRegion.m
diff --git a/+Simulator/@TwoDimensionalMOT/exitCondition.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/exitCondition.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/exitCondition.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/exitCondition.m
diff --git a/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/jackknifeErrorEstimation.m
diff --git a/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/magneticFieldForMOT.m
diff --git a/+Simulator/@TwoDimensionalMOT/runSimulation.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/runSimulation.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/runSimulation.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/runSimulation.m
diff --git a/+Simulator/@TwoDimensionalMOT/solver.m b/2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/solver.m
similarity index 100%
rename from +Simulator/@TwoDimensionalMOT/solver.m
rename to 2DMOT Simulation Code/+Simulator/@TwoDimensionalMOT/solver.m
diff --git a/test_MOTCaptureProcessSimulation.m b/2DMOT Simulation Code/test_MOTCaptureProcessSimulation.m
similarity index 100%
rename from test_MOTCaptureProcessSimulation.m
rename to 2DMOT Simulation Code/test_MOTCaptureProcessSimulation.m
diff --git a/GPE Solver/+Helper/ImageSelection.class b/GPE Solver/+Helper/ImageSelection.class
new file mode 100644
index 0000000..cfac41d
Binary files /dev/null and b/GPE Solver/+Helper/ImageSelection.class differ
diff --git a/GPE Solver/+Helper/ImageSelection.java b/GPE Solver/+Helper/ImageSelection.java
new file mode 100644
index 0000000..5be8fbb
--- /dev/null
+++ b/GPE Solver/+Helper/ImageSelection.java
@@ -0,0 +1,38 @@
+/*
+ * Based on code snippet from
+ * http://java.sun.com/developer/technicalArticles/releases/data/
+ *
+ * Copyright © 2008, 2010 Oracle and/or its affiliates. All rights reserved. Use is subject to license terms.
+ */
+
+import java.awt.image.BufferedImage;
+import java.awt.datatransfer.*;
+
+public class ImageSelection implements Transferable {
+
+ private static final DataFlavor flavors[] =
+ {DataFlavor.imageFlavor};
+
+ private BufferedImage image;
+
+ public ImageSelection(BufferedImage image) {
+ this.image = image;
+ }
+
+ // Transferable
+ public Object getTransferData(DataFlavor flavor) throws UnsupportedFlavorException {
+ if (flavor.equals(flavors[0]) == false) {
+ throw new UnsupportedFlavorException(flavor);
+ }
+ return image;
+ }
+
+ public DataFlavor[] getTransferDataFlavors() {
+ return flavors;
+ }
+
+ public boolean isDataFlavorSupported(DataFlavor
+ flavor) {
+ return flavor.equals(flavors[0]);
+ }
+}
\ No newline at end of file
diff --git a/GPE Solver/+Helper/PhysicsConstants.m b/GPE Solver/+Helper/PhysicsConstants.m
new file mode 100644
index 0000000..f62dd0a
--- /dev/null
+++ b/GPE Solver/+Helper/PhysicsConstants.m
@@ -0,0 +1,46 @@
+classdef PhysicsConstants < handle
+ properties (Constant)
+ % CODATA
+ PlanckConstant=6.62607015E-34;
+ PlanckConstantReduced=6.62607015E-34/(2*pi);
+ FineStructureConstant=7.2973525698E-3;
+ ElectronMass=9.10938291E-31;
+ GravitationalConstant=6.67384E-11;
+ ProtonMass=1.672621777E-27;
+ AtomicMassUnit=1.66053878283E-27;
+ BohrRadius=0.52917721092E-10;
+ BohrMagneton=927.400968E-26;
+ BoltzmannConstant=1.380649E-23;
+ StandardGravityAcceleration=9.80665;
+ SpeedOfLight=299792458;
+ StefanBoltzmannConstant=5.670373E-8;
+ ElectronCharge=1.602176634E-19;
+ VacuumPermeability=1.25663706212E-6;
+ DielectricConstant=8.8541878128E-12;
+ ElectronGyromagneticFactor=-2.00231930436153;
+ AvogadroConstant=6.02214076E23;
+ ZeroKelvin = 273.15;
+ GravitationalAcceleration = 9.80553;
+
+ % Dy specific constants
+ Dy164Mass = 163.929174751*1.66053878283E-27;
+ Dy164IsotopicAbundance = 0.2826;
+ BlueWavelength = 421.291e-9;
+ BlueLandegFactor = 1.22;
+ BlueLifetime = 4.94e-9;
+ BlueLinewidth = 1/4.94e-9;
+ RedWavelength = 626.086e-9;
+ RedLandegFactor = 1.29;
+ RedLifetime = 1.2e-6;
+ RedLinewidth = 1/1.2e-6;
+ PushBeamWaveLength = 626.086e-9;
+ PushBeamLifetime = 1.2e-6;
+ PushBeamLinewidth = 1/1.2e-6;
+ end
+
+ methods
+ function pc = PhysicsConstants()
+ end
+ end
+
+end
diff --git a/GPE Solver/+Helper/bringFiguresWithTagInForeground.m b/GPE Solver/+Helper/bringFiguresWithTagInForeground.m
new file mode 100644
index 0000000..c58a117
--- /dev/null
+++ b/GPE Solver/+Helper/bringFiguresWithTagInForeground.m
@@ -0,0 +1,15 @@
+function output = bringFiguresWithTagInForeground()
+
+figure_handles = findobj('type','figure');
+
+for idx = 1:length(figure_handles)
+ if ~isempty(figure_handles(idx).Tag)
+ figure(figure_handles(idx));
+ end
+end
+
+if nargout > 0
+ output = figure_handles;
+end
+
+end
\ No newline at end of file
diff --git a/GPE Solver/+Helper/calculateDistanceFromPointToLine.m b/GPE Solver/+Helper/calculateDistanceFromPointToLine.m
new file mode 100644
index 0000000..df5c8c6
--- /dev/null
+++ b/GPE Solver/+Helper/calculateDistanceFromPointToLine.m
@@ -0,0 +1,10 @@
+function ret = calculateDistanceFromPointToLine(p0 , p1, p2)
+ p01 = p0 - p1;
+ p12 = p2 - p1;
+ CrossProduct = [p01(2)*p12(3) - p01(3)*p12(2), p01(3)*p12(1) - p01(1)*p12(3), p01(1)*p12(2) - p01(2)*p12(1)];
+ ret = norm(CrossProduct) / norm(p12);
+
+ %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
+ %Area = One side of parallelogram X Base
+ %ret = norm(cross(one side, base))./ norm(base);
+end
\ No newline at end of file
diff --git a/GPE Solver/+Helper/convertstruct2cell.m b/GPE Solver/+Helper/convertstruct2cell.m
new file mode 100644
index 0000000..90fdf2c
--- /dev/null
+++ b/GPE Solver/+Helper/convertstruct2cell.m
@@ -0,0 +1,6 @@
+function CellOut = convertstruct2cell(StructIn)
+ % CellOut = Convertstruct2cell(StructIn)
+ % converts a struct into a cell-matrix where the first column contains
+ % the fieldnames and the second the contents
+ CellOut = [fieldnames(StructIn) struct2cell(StructIn)]';
+end
\ No newline at end of file
diff --git a/GPE Solver/+Helper/findAllZeroCrossings.m b/GPE Solver/+Helper/findAllZeroCrossings.m
new file mode 100644
index 0000000..4b8d9db
--- /dev/null
+++ b/GPE Solver/+Helper/findAllZeroCrossings.m
@@ -0,0 +1,18 @@
+function ret = findAllZeroCrossings(x,y)
+% Finds all Zero-crossing of the function y = f(x)
+ zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
+ zxidx = zci(y);
+ if ~isempty(zxidx)
+ for k1 = 1:numel(zxidx)
+ idxrng = max([1 zxidx(k1)-1]):min([zxidx(k1)+1 numel(y)]);
+ xrng = x(idxrng);
+ yrng = y(idxrng);
+ [yrng2, ~, jyrng] = unique(yrng); %yrng is a new array containing the unique values of yrng. jyrng contains the indices in yrng that correspond to the original vector. yrng = yrng2(jyrng)
+ xrng2 = accumarray(jyrng, xrng, [], @mean); %This function creates a new array "xrng2" by applying the function "@mean" to all elements in "xrng" that have identical indices in "jyrng". Any elements with identical X values will have identical indices in jyrng. Thus, this function creates a new array by averaging values with identical X values in the original array.
+ ret(k1) = interp1( yrng2(:), xrng2(:), 0, 'linear', 'extrap' );
+ end
+ else
+ warning('No zero crossings found!')
+ ret = nan;
+ end
+end
\ No newline at end of file
diff --git a/GPE Solver/+Helper/getFigureByTag.m b/GPE Solver/+Helper/getFigureByTag.m
new file mode 100644
index 0000000..8fc6bf9
--- /dev/null
+++ b/GPE Solver/+Helper/getFigureByTag.m
@@ -0,0 +1,191 @@
+function figure_handle = getFigureByTag(tag_name, varargin)
+ % figure_handle = getFigureByTag(tag_name, varargin)
+ %
+ % Example code:
+ % f_h = getFigureByTag('survivalMeasurement','Name','Survival')
+ %
+ % clf(f_h);
+ % a_h = gca(f_h);
+ % xlim(a_h,[10,100]);
+ % % custom position
+ % f_h.Position = [4052.3 719.67 560 420];
+
+ assert(nargin>=1 && ischar(tag_name),'You must specify ``tag_name'' as a string.');
+
+ f_h = findobj('type','figure','tag',tag_name);
+
+ if isempty(f_h)
+ f_h = figure('Tag',tag_name,varargin{:});
+
+ defaultNewFigProperties = {'Color','w','NumberTitle','off','Name',sprintf('Fig. %d',f_h.Number)};
+
+ varargin = [defaultNewFigProperties,varargin];
+ else
+ f_h = f_h(1);
+ end
+
+ if ~isempty(varargin)
+ set(f_h,varargin{:});
+ end
+
+ addCopyButton(f_h);
+
+ if nargout > 0
+ figure_handle = f_h;
+ else
+ set(groot,'CurrentFigure',f_h);
+ end
+
+end
+
+function addCopyButton(f_h)
+
+ if(strcmp(f_h.ToolBar,'none'))
+ return
+ end
+
+ tb = findall(f_h,'Type','uitoolbar');
+
+ pt = findall(tb, 'tag', 'Custom.CopyPlot' );
+ if isempty(pt)
+ pt = uipushtool(tb);
+ else
+ pt = pt(1);
+ end
+
+ cdata = zeros(16,16,3);
+
+ % Evernote Logo
+% cdata(:,:,1) =[255 NaN NaN NaN NaN 99 11 27 175 NaN NaN NaN NaN NaN NaN 255
+% NaN NaN NaN 251 93 14 0 0 0 66 70 106 210 NaN NaN NaN
+% NaN NaN NaN 42 0 43 0 0 0 0 0 0 20 185 NaN NaN
+% NaN 243 56 0 42 82 0 0 0 0 0 0 0 45 NaN NaN
+% NaN 156 44 64 113 65 0 0 0 0 0 0 0 32 NaN NaN
+% 136 9 26 28 11 0 0 0 0 0 0 0 0 10 188 NaN
+% 132 0 0 0 0 0 0 0 0 0 136 175 16 0 133 NaN
+% NaN 28 0 0 0 0 0 0 0 0 152 238 50 0 124 NaN
+% NaN 58 0 0 0 0 0 0 0 0 0 9 0 0 71 NaN
+% NaN 175 0 0 0 0 0 61 15 0 0 0 0 0 100 NaN
+% NaN NaN 143 12 0 0 0 210 195 87 17 0 0 0 126 NaN
+% NaN NaN NaN 183 118 50 150 NaN NaN 110 219 78 0 0 160 NaN
+% NaN NaN NaN NaN NaN NaN NaN 191 0 35 NaN 150 0 23 NaN NaN
+% NaN NaN NaN NaN NaN NaN NaN 124 0 172 NaN 81 0 93 NaN NaN
+% 255 NaN NaN NaN NaN NaN NaN 183 0 0 0 0 51 228 NaN 245
+% 253 254 NaN NaN NaN NaN NaN NaN 156 63 45 100 NaN NaN 255 255]/255.;
+%
+%
+% cdata(:,:,2) = [255 255 255 255 255 216 166 171 225 229 218 229 247 255 255 255
+% 255 255 255 255 201 166 159 157 167 188 189 200 243 255 255 255
+% 237 238 255 181 159 183 164 170 163 158 160 157 169 233 248 250
+% 224 235 188 140 182 195 161 168 168 168 168 169 147 186 244 240
+% 255 226 175 185 207 189 161 168 168 168 168 168 159 179 249 249
+% 227 172 172 179 172 163 169 168 168 170 163 155 160 173 231 237
+% 215 161 163 165 166 168 168 168 168 162 215 228 172 163 209 219
+% 248 178 159 168 168 168 168 168 168 159 220 249 185 158 208 222
+% 249 192 151 169 168 168 169 160 163 172 163 159 166 167 194 204
+% 246 229 155 157 168 169 159 188 174 154 162 167 166 166 202 214
+% 212 231 218 168 157 153 165 255 242 190 171 159 167 166 207 220
+% 218 203 251 243 206 181 230 210 208 207 242 196 154 168 223 232
+% 255 224 232 250 237 214 244 194 152 178 255 223 145 175 250 252
+% 255 255 244 239 222 213 240 214 149 228 254 199 136 203 244 232
+% 255 255 255 246 231 246 246 232 165 159 167 147 184 253 254 242
+% 253 254 255 255 254 255 255 255 231 183 178 199 249 255 255 255]/255.;
+%
+%
+% cdata(:,:,3) = [255 255 255 255 255 117 38 50 187 211 170 190 234 255 255 255
+% 255 254 255 255 120 51 27 20 39 97 98 122 220 255 255 255
+% 238 252 246 73 22 71 37 49 35 20 24 18 49 196 231 231
+% 232 242 86 0 78 108 29 45 45 45 45 46 0 82 214 201
+% 255 175 63 85 139 98 27 45 45 45 45 45 23 72 233 231
+% 167 51 57 72 55 32 47 45 45 50 34 14 27 57 201 218
+% 154 30 33 38 39 45 45 45 45 31 157 188 53 34 153 180
+% 234 67 24 45 45 45 45 44 45 24 169 241 83 20 146 182
+% 241 99 4 48 45 45 47 28 35 53 32 26 39 44 104 127
+% 238 192 14 20 45 47 27 97 56 10 29 44 41 40 127 158
+% 214 253 169 37 20 16 34 218 207 105 55 23 42 40 147 182
+% 218 214 241 201 138 71 177 225 181 130 224 107 12 45 175 197
+% 255 233 202 218 212 132 230 196 27 61 255 172 0 64 240 242
+% 255 255 219 197 176 160 237 143 0 195 245 110 0 123 230 230
+% 255 255 255 227 197 241 244 202 36 24 39 0 81 228 242 245
+% 253 254 255 255 254 255 255 255 191 78 71 121 221 255 255 255]/255.;
+
+ %OneNote logo
+
+ cdata(:,:,1) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 245 213 213 213 213 213 213 213 184 184 215 255
+ 255 255 255 255 241 213 213 213 213 213 213 213 184 184 208 255
+ 255 233 204 204 194 176 176 185 213 213 213 213 184 184 208 255
+ 255 154 101 101 101 101 101 103 213 213 213 206 162 162 193 255
+ 255 152 101 183 116 152 115 101 213 213 213 206 162 162 193 255
+ 255 152 101 207 189 178 122 101 213 213 213 206 162 162 193 255
+ 255 152 101 199 152 224 122 101 213 213 213 195 128 128 170 255
+ 255 152 101 166 101 183 115 101 213 213 213 195 128 128 170 255
+ 255 154 101 101 101 101 101 103 213 213 213 195 128 128 170 255
+ 255 233 204 204 194 176 176 185 213 213 213 183 95 95 148 255
+ 255 255 255 255 241 213 213 213 213 213 213 183 94 94 148 255
+ 255 255 255 255 245 213 213 213 213 213 213 183 94 94 163 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
+
+
+ cdata(:,:,2) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 219 112 110 110 110 110 110 134 84 84 158 255
+ 255 255 255 255 207 110 110 110 110 110 110 134 84 84 141 255
+ 255 222 178 178 146 81 81 88 110 110 110 134 84 84 141 255
+ 255 102 23 23 23 23 23 24 110 110 110 125 58 58 123 255
+ 255 100 23 147 46 100 44 23 110 110 110 125 58 58 123 255
+ 255 100 23 183 156 139 55 23 110 110 110 125 58 58 123 255
+ 255 100 23 170 99 208 55 23 110 110 110 119 38 38 109 255
+ 255 100 23 121 23 146 44 23 110 110 110 119 38 38 109 255
+ 255 102 23 23 23 23 23 24 110 110 110 119 38 38 109 255
+ 255 222 178 178 146 81 81 88 110 110 110 118 37 37 109 255
+ 255 255 255 255 207 110 110 110 110 110 110 118 37 37 110 255
+ 255 255 255 255 219 112 110 110 110 110 110 118 37 37 131 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
+
+
+ cdata(:,:,3) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 246 229 229 240 255
+ 255 255 255 255 255 255 255 255 255 255 255 246 229 229 238 255
+ 255 242 224 224 224 224 224 232 255 255 255 246 229 229 238 255
+ 255 194 163 163 163 163 163 164 255 255 255 244 223 223 234 255
+ 255 194 163 212 172 194 171 163 255 255 255 244 223 223 234 255
+ 255 194 163 226 216 209 176 163 255 255 255 244 223 223 234 255
+ 255 194 163 221 193 236 176 163 255 255 255 240 209 209 224 255
+ 255 194 163 202 163 212 171 163 255 255 255 240 209 209 224 255
+ 255 194 163 163 163 163 163 164 255 255 255 240 209 209 224 255
+ 255 242 224 224 224 224 224 232 255 255 255 223 161 161 192 255
+ 255 255 255 255 255 255 255 255 255 255 255 223 160 160 192 255
+ 255 255 255 255 255 255 255 255 255 255 255 223 160 160 201 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255
+ 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.;
+
+
+ pt.Tag = 'Custom.CopyPlot';
+ pt.CData = cdata;
+ pt.Separator = true;
+ pt.ClickedCallback = @copyToClipboard;
+
+end
+
+function copyToClipboard(~,~)
+ fig_h = get(get(gcbo,'Parent'),'Parent');
+ if strcmp(fig_h.WindowStyle,'docked')
+ if ismac || ispc
+ matlab.graphics.internal.copyFigureHelper(fig_h);
+ else
+ %warning('Copy function to the clipboard only works if the figure is undocked.');
+ Helper.screencapture(fig_h,[],'clipboard');
+ end
+ else
+ pos = fig_h.Position;
+ Helper.screencapture(fig_h,[],'clipboard','position',[7,7,pos(3)-2,pos(4)]);
+ end
+end
+
+
+
diff --git a/GPE Solver/+Helper/ode5.m b/GPE Solver/+Helper/ode5.m
new file mode 100644
index 0000000..3eb003f
--- /dev/null
+++ b/GPE Solver/+Helper/ode5.m
@@ -0,0 +1,92 @@
+function Y = ode5(odefun,tspan,y0,varargin)
+%ODE5 Solve differential equations with a non-adaptive method of order 5.
+% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates
+% the system of differential equations y' = f(t,y) by stepping from T0 to
+% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
+% The vector Y0 is the initial conditions at T0. Each row in the solution
+% array Y corresponds to a time specified in TSPAN.
+%
+% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
+% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
+%
+% This is a non-adaptive solver. The step sequence is determined by TSPAN
+% but the derivative function ODEFUN is evaluated multiple times per step.
+% The solver implements the Dormand-Prince method of order 5 in a general
+% framework of explicit Runge-Kutta methods.
+%
+% Example
+% tspan = 0:0.1:20;
+% y = ode5(@vdp1,tspan,[2 0]);
+% plot(tspan,y(:,1));
+% solves the system y' = vdp1(t,y) with a constant step size of 0.1,
+% and plots the first component of the solution.
+
+if ~isnumeric(tspan)
+ error('TSPAN should be a vector of integration steps.');
+end
+
+if ~isnumeric(y0)
+ error('Y0 should be a vector of initial conditions.');
+end
+
+h = diff(tspan);
+if any(sign(h(1))*h <= 0)
+ error('Entries of TSPAN are not in order.')
+end
+
+try
+ f0 = feval(odefun,tspan(1),y0,varargin{:});
+catch
+ msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
+ error(msg);
+end
+
+y0 = y0(:); % Make a column vector.
+if ~isequal(size(y0),size(f0))
+ error('Inconsistent sizes of Y0 and f(t0,y0).');
+end
+
+neq = length(y0);
+N = length(tspan);
+Y = zeros(neq,N);
+
+% Method coefficients -- Butcher's tableau
+%
+% C | A
+% --+---
+% | B
+
+C = [1/5; 3/10; 4/5; 8/9; 1];
+
+A = [ 1/5, 0, 0, 0, 0
+ 3/40, 9/40, 0, 0, 0
+ 44/45 -56/15, 32/9, 0, 0
+ 19372/6561, -25360/2187, 64448/6561, -212/729, 0
+ 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656];
+
+B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84];
+
+% More convenient storage
+A = A.';
+B = B(:);
+
+nstages = length(B);
+F = zeros(neq,nstages);
+
+Y(:,1) = y0;
+for i = 2:N
+ ti = tspan(i-1);
+ hi = h(i-1);
+ yi = Y(:,i-1);
+
+ % General explicit Runge-Kutta framework
+ F(:,1) = feval(odefun,ti,yi,varargin{:});
+ for stage = 2:nstages
+ tstage = ti + C(stage-1)*hi;
+ ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1));
+ F(:,stage) = feval(odefun,tstage,ystage,varargin{:});
+ end
+ Y(:,i) = yi + F*(hi*B);
+
+end
+Y = Y.';
diff --git a/GPE Solver/+Helper/onenoteccdata.m b/GPE Solver/+Helper/onenoteccdata.m
new file mode 100644
index 0000000..5f3c177
--- /dev/null
+++ b/GPE Solver/+Helper/onenoteccdata.m
@@ -0,0 +1,55 @@
+cmap = zeros(16,16,3);
+
+cmap(:,:,1) = [0.0000 0.0118 0.4510 0.0039 0.2078 0.1569 0.4078 0.4431 0.4510 0.1922 0.4235 0.4196 0.2235 0.4235 0.4039 0.4392
+ 0.4471 0.1647 0.4157 0.0000 0.0235 0.4353 0.0314 0.4314 0.0196 0.2392 0.0667 0.0392 0.4431 0.3804 0.2941 0.4275
+ 0.3686 0.3608 0.2000 0.2824 0.3059 0.0549 0.1804 0.1882 0.4392 0.4314 0.3255 0.0078 0.0902 0.1961 0.4353 0.1412
+ 0.2314 0.3647 0.0353 0.3804 0.1647 0.2431 0.1686 0.2745 0.2980 0.4235 0.3922 0.4157 0.2784 0.3333 0.2510 0.0588
+ 0.1020 0.0745 0.2549 0.0471 0.1216 0.4000 0.3961 0.2627 0.1098 0.1725 0.3098 0.4314 0.3529 0.3412 0.0784 0.0824
+ 0.4471 0.1490 0.1804 0.3529 0.2196 0.3137 0.3255 0.0941 0.0078 0.3294 0.3765 0.2706 0.0510 0.0157 0.4275 0.1176
+ 0.1294 0.1333 0.1725 0.3451 0.2118 0.3843 0.1255 0.1569 0.2118 0.1608 0.0353 0.2039 0.1608 0.4510 1.0000 0.8000
+ 0.9882 0.6510 0.9961 0.4549 0.4549 0.6824 0.7882 0.5686 0.5373 0.5490 0.7765 0.7137 0.8510 0.7176 0.5020 0.4902
+ 0.8941 0.9020 0.4745 0.8980 0.9098 0.4824 0.6471 0.6353 0.9922 0.9647 0.6353 0.4588 0.9647 0.9020 0.4980 0.8118
+ 0.5059 0.4941 0.9686 0.4863 0.5451 0.9725 0.8980 0.5451 0.5333 0.6824 0.4588 0.8196 0.8314 0.8980 0.8941 0.9961
+ 0.5255 0.8392 0.9804 0.5216 0.8588 0.8078 0.5176 0.7647 0.5608 0.9725 0.9059 0.4627 0.9882 0.8275 0.7725 0.8745
+ 0.8235 0.8431 0.7373 1.0000 0.5137 0.4706 0.4784 0.7412 0.8863 0.9373 0.5529 0.5804 0.4510 0.9255 0.8235 0.8667
+ 0.7569 0.8824 0.5294 0.5176 0.5373 0.9569 0.5294 0.4824 0.5098 0.5137 0.5569 0.8471 0.5098 0.9490 0.8706 0.9412
+ 0.4902 0.6000 0.6980 0.7882 0.5490 0.7216 0.6431 0.4824 0.5569 0.4667 0.6627 0.9922 0.7804 0.8039 0.6275 0.7333
+ 0.5725 0.5647 0.8549 0.7529 0.6235 0.8784 0.5922 0.7294 0.6118 0.7922 0.7843 0.6667 0.9294 0.6902 0.6784 0.9176
+ 0.6706 0.7490 0.7961 0.5882 0.8627 0.4627 0.6196 0.7059 0.6078 0.9765 0.6549 0.6863 0.5373 0.7098 0.7176 0.7765];
+
+cmap(:,:,2) = [0.0000 0.0078 0.2157 0.0000 0.0980 0.0745 0.1922 0.2157 0.2157 0.0902 0.2000 0.1961 0.1059 0.2039 0.1882 0.2078
+ 0.2078 0.0784 0.2000 0.0000 0.0118 0.2118 0.0157 0.2039 0.0078 0.1137 0.0314 0.0196 0.2118 0.1804 0.1373 0.2078
+ 0.1765 0.1725 0.0941 0.1333 0.1451 0.0275 0.0863 0.0902 0.2078 0.2078 0.1529 0.0039 0.0431 0.0941 0.2039 0.0667
+ 0.1098 0.1725 0.0157 0.1804 0.0784 0.1137 0.0824 0.1333 0.1412 0.2000 0.1882 0.2000 0.1333 0.1569 0.1176 0.0275
+ 0.0471 0.0353 0.1216 0.0196 0.0588 0.1922 0.1882 0.1255 0.0510 0.0824 0.1451 0.2039 0.1686 0.1647 0.0392 0.0392
+ 0.2157 0.0706 0.0863 0.1686 0.1020 0.1490 0.1529 0.0431 0.0039 0.1569 0.1804 0.1255 0.0235 0.0078 0.2000 0.0549
+ 0.0627 0.0627 0.0824 0.1647 0.1020 0.1843 0.0588 0.0745 0.1020 0.0784 0.0157 0.0980 0.0784 0.2157 1.0000 0.7137
+ 0.9843 0.4980 0.9961 0.2235 0.2196 0.5412 0.6980 0.3843 0.3373 0.3569 0.6824 0.5922 0.7843 0.6000 0.2902 0.2706
+ 0.8510 0.8588 0.2471 0.8549 0.8667 0.2627 0.4980 0.4784 0.9843 0.9490 0.4745 0.2235 0.9451 0.8627 0.2824 0.7333
+ 0.2941 0.2784 0.9529 0.2667 0.3490 0.9569 0.8510 0.3490 0.3333 0.5451 0.2275 0.7412 0.7608 0.8549 0.8471 0.9922
+ 0.3255 0.7686 0.9725 0.3176 0.8000 0.7255 0.3098 0.6627 0.3725 0.9647 0.8627 0.2314 0.9804 0.7529 0.6745 0.8235
+ 0.7451 0.7765 0.6235 0.9961 0.3020 0.2431 0.2510 0.6314 0.8392 0.9098 0.3608 0.4000 0.2196 0.8902 0.7490 0.8078
+ 0.6549 0.8353 0.3294 0.3137 0.3412 0.9373 0.3255 0.2588 0.2980 0.3059 0.3686 0.7843 0.3020 0.9255 0.8157 0.9176
+ 0.2745 0.4275 0.5686 0.6980 0.3569 0.6039 0.4863 0.2627 0.3647 0.2392 0.5137 0.9922 0.6863 0.7216 0.4706 0.6196
+ 0.3882 0.3765 0.7882 0.6471 0.4588 0.8275 0.4157 0.6118 0.4431 0.7059 0.6902 0.5255 0.8980 0.5569 0.5412 0.8824
+ 0.5333 0.6392 0.7098 0.4078 0.8039 0.2314 0.4549 0.5804 0.4392 0.9647 0.5059 0.5529 0.3373 0.5882 0.5961 0.6784];
+
+cmap(:,:,3) = [0.0000 0.0157 0.4980 0.0039 0.2314 0.1725 0.4627 0.5020 0.5020 0.2196 0.4745 0.4706 0.2510 0.4784 0.4510 0.4980
+ 0.4941 0.1882 0.4667 0.0000 0.0275 0.4941 0.0353 0.4902 0.0196 0.2667 0.0745 0.0471 0.4902 0.4314 0.3294 0.4784
+ 0.4196 0.4000 0.2235 0.3216 0.3412 0.0627 0.2039 0.2118 0.4863 0.4863 0.3608 0.0078 0.1020 0.2196 0.4824 0.1569
+ 0.2588 0.4118 0.0392 0.4235 0.1843 0.2745 0.1882 0.3059 0.3373 0.4784 0.4392 0.4627 0.3137 0.3765 0.2824 0.0667
+ 0.1137 0.0824 0.2863 0.0510 0.1373 0.4510 0.4471 0.2941 0.1216 0.1961 0.3490 0.4824 0.3961 0.3804 0.0902 0.0941
+ 0.4980 0.1647 0.2000 0.4000 0.2431 0.3529 0.3647 0.1059 0.0118 0.3686 0.4196 0.3020 0.0549 0.0196 0.4824 0.1294
+ 0.1451 0.1529 0.1922 0.3882 0.2392 0.4353 0.1412 0.1765 0.2353 0.1804 0.0353 0.2275 0.1843 0.5059 1.0000 0.8196
+ 0.9882 0.6863 0.9961 0.5098 0.5098 0.7137 0.8118 0.6118 0.5843 0.5922 0.8000 0.7412 0.8627 0.7451 0.5529 0.5412
+ 0.9059 0.9137 0.5255 0.9098 0.9176 0.5333 0.6824 0.6706 0.9922 0.9686 0.6706 0.5098 0.9647 0.9137 0.5490 0.8314
+ 0.5569 0.5451 0.9725 0.5373 0.5922 0.9725 0.9059 0.5882 0.5804 0.7137 0.5137 0.8353 0.8510 0.9059 0.9020 0.9961
+ 0.5725 0.8549 0.9843 0.5725 0.8745 0.8275 0.5647 0.7882 0.6039 0.9765 0.9137 0.5176 0.9882 0.8431 0.7961 0.8863
+ 0.8392 0.8588 0.7647 1.0000 0.5608 0.5216 0.5294 0.7686 0.8980 0.9412 0.6000 0.6235 0.5059 0.9333 0.8431 0.8784
+ 0.7804 0.8941 0.5765 0.5686 0.5843 0.9608 0.5765 0.5333 0.5569 0.5647 0.6039 0.8627 0.5608 0.9569 0.8863 0.9490
+ 0.5412 0.6392 0.7294 0.8078 0.5961 0.7490 0.6784 0.5373 0.6000 0.5216 0.6941 0.9922 0.8039 0.8235 0.6667 0.7608
+ 0.6157 0.6078 0.8667 0.7765 0.6588 0.8902 0.6314 0.7569 0.6510 0.8157 0.8039 0.7020 0.9373 0.7216 0.7098 0.9255
+ 0.7059 0.7725 0.8196 0.6314 0.8784 0.5137 0.6549 0.7373 0.6471 0.9804 0.6902 0.7176 0.5804 0.7412 0.7451 0.8000];
+
+%%
+[cdata, cmap] = imread('onenote.png');
\ No newline at end of file
diff --git a/GPE Solver/+Helper/parforNotifications.m b/GPE Solver/+Helper/parforNotifications.m
new file mode 100644
index 0000000..4ad3af4
--- /dev/null
+++ b/GPE Solver/+Helper/parforNotifications.m
@@ -0,0 +1,148 @@
+% Copyright (c) 2019 Andrea Alberti
+%
+% All rights reserved.
+classdef parforNotifications < handle
+ properties
+ N; % number of iterations
+ text = 'Please wait ...'; % text to show
+ width = 50;
+ showWarning = true;
+ end
+ properties (GetAccess = public, SetAccess = private)
+ n;
+ end
+ properties (Access = private)
+ inProgress = false;
+ percent;
+ DataQueue;
+ usePercent;
+ Nstr;
+ NstrL;
+ lastComment;
+ end
+ methods
+ function this = parforNotifications()
+ this.DataQueue = parallel.pool.DataQueue;
+ afterEach(this.DataQueue, @this.updateStatus);
+ end
+ % Start progress bar
+ function PB_start(this,N,varargin)
+ assert(isscalar(N) && isnumeric(N) && N == floor(N) && N>0, 'Error: ''N'' must be a scalar positive integer.');
+
+ this.N = N;
+
+ p = inputParser;
+ addParameter(p,'message','Please wait: ');
+ addParameter(p,'usePercentage',true);
+
+ parse(p,varargin{:});
+
+ this.text = p.Results.message;
+ assert(ischar(this.text), 'Error: ''Message'' must be a string.');
+
+ this.usePercent = p.Results.usePercentage;
+ assert(isscalar(this.usePercent) && islogical(this.usePercent), 'Error: ''usePercentage'' must be a logical scalar.');
+
+ this.percent = 0;
+ this.n = 0;
+ this.lastComment = '';
+ if this.usePercent
+ fprintf('%s [%s]: %3d%%\n',this.text, char(32*ones(1,this.width)),0);
+ else
+ this.Nstr = sprintf('%d',this.N);
+ this.NstrL = numel(this.Nstr);
+ fprintf('%s [%s]: %s/%s\n',this.text, char(32*ones(1,this.width)),[char(32*ones(1,this.NstrL-1)),'0'],this.Nstr);
+ end
+
+ this.inProgress = true;
+ end
+ % Iterate progress bar
+ function PB_iterate(this,str)
+ if nargin == 1
+ send(this.DataQueue,'');
+ else
+ send(this.DataQueue,str);
+ end
+ end
+ function warning(this,warn_id,msg)
+ if this.showWarning
+ msg = struct('Action','Warning','Id',warn_id,'Message',msg);
+ send(this.DataQueue,msg);
+ end
+ end
+ function PB_reprint(this)
+ p = round(100*this.n/this.N);
+
+ this.percent = p;
+
+ cursor_pos=1+round((this.width-1)*p/100);
+
+ if p < 100
+ sep_char = '|';
+ else
+ sep_char = '.';
+ end
+
+ if this.usePercent
+ fprintf('%s [%s%s%s]: %3d%%\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p);
+ else
+ nstr=sprintf('%d',this.n);
+ fprintf('%s [%s%s%s]: %s/%s\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr);
+ end
+ end
+ function updateStatus(this,data)
+
+ if ischar(data)
+
+ this.n = this.n + 1;
+
+ p = round(100*this.n/this.N);
+
+ if p >= this.percent+1 || this.n == this.N
+ this.percent = p;
+
+ cursor_pos=1+round((this.width-1)*p/100);
+
+ if p < 100
+ sep_char = '|';
+ else
+ sep_char = '.';
+ end
+
+ if ~isempty(data)
+ comment = [' (',data,')'];
+ else
+ comment = '';
+ end
+
+ if this.usePercent
+ fprintf('%s%s%s%s]: %3d%%%s\n',char(8*ones(1,58+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p,comment);
+ else
+ nstr=sprintf('%d',this.n);
+ fprintf('%s%s%s%s]: %s/%s%s\n',char(8*ones(1,55+2*numel(this.Nstr)+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr,comment)
+ end
+
+ this.lastComment = comment;
+
+
+ if p == 100
+ this.inProgress = false;
+ end
+ end
+
+ else
+ switch data.Action
+ case 'Warning'
+ warning(data.Id,[data.Message,newline]);
+ if this.inProgress
+ this.PB_reprint();
+ end
+ end
+
+ end
+
+ end
+ end
+end
+
+
diff --git a/GPE Solver/+Helper/screencapture.m b/GPE Solver/+Helper/screencapture.m
new file mode 100644
index 0000000..206312e
--- /dev/null
+++ b/GPE Solver/+Helper/screencapture.m
@@ -0,0 +1,820 @@
+function imageData = screencapture(varargin)
+% screencapture - get a screen-capture of a figure frame, component handle, or screen area rectangle
+%
+% ScreenCapture gets a screen-capture of any Matlab GUI handle (including desktop,
+% figure, axes, image or uicontrol), or a specified area rectangle located relative to
+% the specified handle. Screen area capture is possible by specifying the root (desktop)
+% handle (=0). The output can be either to an image file or to a Matlab matrix (useful
+% for displaying via imshow() or for further processing) or to the system clipboard.
+% This utility also enables adding a toolbar button for easy interactive screen-capture.
+%
+% Syntax:
+% imageData = screencapture(handle, position, target, 'PropName',PropValue, ...)
+%
+% Input Parameters:
+% handle - optional handle to be used for screen-capture origin.
+% If empty/unsupplied then current figure (gcf) will be used.
+% position - optional position array in pixels: [x,y,width,height].
+% If empty/unsupplied then the handle's position vector will be used.
+% If both handle and position are empty/unsupplied then the position
+% will be retrieved via interactive mouse-selection.
+% If handle is an image, then position is in data (not pixel) units, so the
+% captured region remains the same after figure/axes resize (like imcrop)
+% target - optional filename for storing the screen-capture, or the
+% 'clipboard'/'printer' strings.
+% If empty/unsupplied then no output to file will be done.
+% The file format will be determined from the extension (JPG/PNG/...).
+% Supported formats are those supported by the imwrite function.
+% 'PropName',PropValue -
+% optional list of property pairs (e.g., 'target','myImage.png','pos',[10,20,30,40],'handle',gca)
+% PropNames may be abbreviated and are case-insensitive.
+% PropNames may also be given in whichever order.
+% Supported PropNames are:
+% - 'handle' (default: gcf handle)
+% - 'position' (default: gcf position array)
+% - 'target' (default: '')
+% - 'toolbar' (figure handle; default: gcf)
+% this adds a screen-capture button to the figure's toolbar
+% If this parameter is specified, then no screen-capture
+% will take place and the returned imageData will be [].
+%
+% Output parameters:
+% imageData - image data in a format acceptable by the imshow function
+% If neither target nor imageData were specified, the user will be
+% asked to interactively specify the output file.
+%
+% Examples:
+% imageData = screencapture; % interactively select screen-capture rectangle
+% imageData = screencapture(hListbox); % capture image of a uicontrol
+% imageData = screencapture(0, [20,30,40,50]); % capture a small desktop region
+% imageData = screencapture(gcf,[20,30,40,50]); % capture a small figure region
+% imageData = screencapture(gca,[10,20,30,40]); % capture a small axes region
+% imshow(imageData); % display the captured image in a matlab figure
+% imwrite(imageData,'myImage.png'); % save the captured image to file
+% img = imread('cameraman.tif');
+% hImg = imshow(img);
+% screencapture(hImg,[60,35,140,80]); % capture a region of an image
+% screencapture(gcf,[],'myFigure.jpg'); % capture the entire figure into file
+% screencapture(gcf,[],'clipboard'); % capture the entire figure into clipboard
+% screencapture(gcf,[],'printer'); % print the entire figure
+% screencapture('handle',gcf,'target','myFigure.jpg'); % same as previous, save to file
+% screencapture('handle',gcf,'target','clipboard'); % same as previous, copy to clipboard
+% screencapture('handle',gcf,'target','printer'); % same as previous, send to printer
+% screencapture('toolbar',gcf); % adds a screen-capture button to gcf's toolbar
+% screencapture('toolbar',[],'target','sc.bmp'); % same with default output filename
+%
+% Technical description:
+% http://UndocumentedMatlab.com/blog/screencapture-utility/
+%
+% Bugs and suggestions:
+% Please send to Yair Altman (altmany at gmail dot com)
+%
+% See also:
+% imshow, imwrite, print
+%
+% Release history:
+% 1.17 2016-05-16: Fix annoying warning about JavaFrame property becoming obsolete someday (yes, we know...)
+% 1.16 2016-04-19: Fix for deployed application suggested by Dwight Bartholomew
+% 1.10 2014-11-25: Added the 'print' target
+% 1.9 2014-11-25: Fix for saving GIF files
+% 1.8 2014-11-16: Fixes for R2014b
+% 1.7 2014-04-28: Fixed bug when capturing interactive selection
+% 1.6 2014-04-22: Only enable image formats when saving to an unspecified file via uiputfile
+% 1.5 2013-04-18: Fixed bug in capture of non-square image; fixes for Win64
+% 1.4 2013-01-27: Fixed capture of Desktop (root); enabled rbbox anywhere on desktop (not necesarily in a Matlab figure); enabled output to clipboard (based on Jiro Doke's imclipboard utility); edge-case fixes; added Java compatibility check
+% 1.3 2012-07-23: Capture current object (uicontrol/axes/figure) if w=h=0 (e.g., by clicking a single point); extra input args sanity checks; fix for docked windows and image axes; include axes labels & ticks by default when capturing axes; use data-units position vector when capturing images; many edge-case fixes
+% 1.2 2011-01-16: another performance boost (thanks to Jan Simon); some compatibility fixes for Matlab 6.5 (untested)
+% 1.1 2009-06-03: Handle missing output format; performance boost (thanks to Urs); fix minor root-handle bug; added toolbar button option
+% 1.0 2009-06-02: First version posted on MathWorks File Exchange
+
+% License to use and modify this code is granted freely to all interested, as long as the original author is
+% referenced and attributed as such. The original author maintains the right to be solely associated with this work.
+
+% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com
+% $Revision: 1.17 $ $Date: 2016/05/16 17:59:36 $
+
+ % Ensure that java awt is enabled...
+ if ~usejava('awt')
+ error('YMA:screencapture:NeedAwt','ScreenCapture requires Java to run.');
+ end
+
+ % Ensure that our Java version supports the Robot class (requires JVM 1.3+)
+ try
+ robot = java.awt.Robot; %#ok
+ catch
+ uiwait(msgbox({['Your Matlab installation is so old that its Java engine (' version('-java') ...
+ ') does not have a java.awt.Robot class. '], ' ', ...
+ 'Without this class, taking a screen-capture is impossible.', ' ', ...
+ 'So, either install JVM 1.3 or higher, or use a newer Matlab release.'}, ...
+ 'ScreenCapture', 'warn'));
+ if nargout, imageData = []; end
+ return;
+ end
+
+ % Process optional arguments
+ paramsStruct = processArgs(varargin{:});
+
+ % If toolbar button requested, add it and exit
+ if ~isempty(paramsStruct.toolbar)
+
+ % Add the toolbar button
+ addToolbarButton(paramsStruct);
+
+ % Return the figure to its pre-undocked state (when relevant)
+ redockFigureIfRelevant(paramsStruct);
+
+ % Exit immediately (do NOT take a screen-capture)
+ if nargout, imageData = []; end
+ return;
+ end
+
+ % Convert position from handle-relative to desktop Java-based pixels
+ [paramsStruct, msgStr] = convertPos(paramsStruct);
+
+ % Capture the requested screen rectangle using java.awt.Robot
+ imgData = getScreenCaptureImageData(paramsStruct.position);
+
+ % Return the figure to its pre-undocked state (when relevant)
+ redockFigureIfRelevant(paramsStruct);
+
+ % Save image data in file or clipboard, if specified
+ if ~isempty(paramsStruct.target)
+ if strcmpi(paramsStruct.target,'clipboard')
+ if ~isempty(imgData)
+ imclipboard(imgData);
+ else
+ msgbox('No image area selected - not copying image to clipboard','ScreenCapture','warn');
+ end
+ elseif strncmpi(paramsStruct.target,'print',5) % 'print' or 'printer'
+ if ~isempty(imgData)
+ hNewFig = figure('visible','off');
+ imshow(imgData);
+ print(hNewFig);
+ delete(hNewFig);
+ else
+ msgbox('No image area selected - not printing screenshot','ScreenCapture','warn');
+ end
+ else % real filename
+ if ~isempty(imgData)
+ imwrite(imgData,paramsStruct.target);
+ else
+ msgbox(['No image area selected - not saving image file ' paramsStruct.target],'ScreenCapture','warn');
+ end
+ end
+ end
+
+ % Return image raster data to user, if requested
+ if nargout
+ imageData = imgData;
+
+ % If neither output formats was specified (neither target nor output data)
+ elseif isempty(paramsStruct.target) & ~isempty(imgData) %#ok ML6
+ % Ask the user to specify a file
+ %error('YMA:screencapture:noOutput','No output specified for ScreenCapture: specify the output filename and/or output data');
+ %format = '*.*';
+ formats = imformats;
+ for idx = 1 : numel(formats)
+ ext = sprintf('*.%s;',formats(idx).ext{:});
+ format(idx,1:2) = {ext(1:end-1), formats(idx).description}; %#ok
+ end
+ [filename,pathname] = uiputfile(format,'Save screen capture as');
+ if ~isequal(filename,0) & ~isequal(pathname,0) %#ok Matlab6 compatibility
+ try
+ filename = fullfile(pathname,filename);
+ imwrite(imgData,filename);
+ catch % possibly a GIF file that requires indexed colors
+ [imgData,map] = rgb2ind(imgData,256);
+ imwrite(imgData,map,filename);
+ end
+ else
+ % TODO - copy to clipboard
+ end
+ end
+
+ % Display msgStr, if relevant
+ if ~isempty(msgStr)
+ uiwait(msgbox(msgStr,'ScreenCapture'));
+ drawnow; pause(0.05); % time for the msgbox to disappear
+ end
+
+ return; % debug breakpoint
+
+%% Process optional arguments
+function paramsStruct = processArgs(varargin)
+
+ % Get the properties in either direct or P-V format
+ [regParams, pvPairs] = parseparams(varargin);
+
+ % Now process the optional P-V params
+ try
+ % Initialize
+ paramName = [];
+ paramsStruct = [];
+ paramsStruct.handle = [];
+ paramsStruct.position = [];
+ paramsStruct.target = '';
+ paramsStruct.toolbar = [];
+ paramsStruct.wasDocked = 0; % no false available in ML6
+ paramsStruct.wasInteractive = 0; % no false available in ML6
+
+ % Parse the regular (non-named) params in recption order
+ if ~isempty(regParams) & (isempty(regParams{1}) | ishandle(regParams{1}(1))) %#ok ML6
+ paramsStruct.handle = regParams{1};
+ regParams(1) = [];
+ end
+ if ~isempty(regParams) & isnumeric(regParams{1}) & (length(regParams{1}) == 4) %#ok ML6
+ paramsStruct.position = regParams{1};
+ regParams(1) = [];
+ end
+ if ~isempty(regParams) & ischar(regParams{1}) %#ok ML6
+ paramsStruct.target = regParams{1};
+ end
+
+ % Parse the optional param PV pairs
+ supportedArgs = {'handle','position','target','toolbar'};
+ while ~isempty(pvPairs)
+
+ % Disregard empty propNames (may be due to users mis-interpretting the syntax help)
+ while ~isempty(pvPairs) & isempty(pvPairs{1}) %#ok ML6
+ pvPairs(1) = [];
+ end
+ if isempty(pvPairs)
+ break;
+ end
+
+ % Ensure basic format is valid
+ paramName = '';
+ if ~ischar(pvPairs{1})
+ error('YMA:screencapture:invalidProperty','Invalid property passed to ScreenCapture');
+ elseif length(pvPairs) == 1
+ if isempty(paramsStruct.target)
+ paramsStruct.target = pvPairs{1};
+ break;
+ else
+ error('YMA:screencapture:noPropertyValue',['No value specified for property ''' pvPairs{1} '''']);
+ end
+ end
+
+ % Process parameter values
+ paramName = pvPairs{1};
+ if strcmpi(paramName,'filename') % backward compatibility
+ paramName = 'target';
+ end
+ paramValue = pvPairs{2};
+ pvPairs(1:2) = [];
+ idx = find(strncmpi(paramName,supportedArgs,length(paramName)));
+ if ~isempty(idx)
+ %paramsStruct.(lower(supportedArgs{idx(1)})) = paramValue; % incompatible with ML6
+ paramsStruct = setfield(paramsStruct, lower(supportedArgs{idx(1)}), paramValue); %#ok ML6
+
+ % If 'toolbar' param specified, then it cannot be left empty - use gcf
+ if strncmpi(paramName,'toolbar',length(paramName)) & isempty(paramsStruct.toolbar) %#ok ML6
+ paramsStruct.toolbar = getCurrentFig;
+ end
+
+ elseif isempty(paramsStruct.target)
+ paramsStruct.target = paramName;
+ pvPairs = {paramValue, pvPairs{:}}; %#ok (more readable this way, although a bit less efficient...)
+
+ else
+ supportedArgsStr = sprintf('''%s'',',supportedArgs{:});
+ error('YMA:screencapture:invalidProperty','%s \n%s', ...
+ 'Invalid property passed to ScreenCapture', ...
+ ['Supported property names are: ' supportedArgsStr(1:end-1)]);
+ end
+ end % loop pvPairs
+
+ catch
+ if ~isempty(paramName), paramName = [' ''' paramName '''']; end
+ error('YMA:screencapture:invalidProperty','Error setting ScreenCapture property %s:\n%s',paramName,lasterr); %#ok
+ end
+%end % processArgs
+
+%% Convert position from handle-relative to desktop Java-based pixels
+function [paramsStruct, msgStr] = convertPos(paramsStruct)
+ msgStr = '';
+ try
+ % Get the screen-size for later use
+ screenSize = get(0,'ScreenSize');
+
+ % Get the containing figure's handle
+ hParent = paramsStruct.handle;
+ if isempty(paramsStruct.handle)
+ paramsStruct.hFigure = getCurrentFig;
+ hParent = paramsStruct.hFigure;
+ else
+ paramsStruct.hFigure = ancestor(paramsStruct.handle,'figure');
+ end
+
+ % To get the acurate pixel position, the figure window must be undocked
+ try
+ if strcmpi(get(paramsStruct.hFigure,'WindowStyle'),'docked')
+ set(paramsStruct.hFigure,'WindowStyle','normal');
+ drawnow; pause(0.25);
+ paramsStruct.wasDocked = 1; % no true available in ML6
+ end
+ catch
+ % never mind - ignore...
+ end
+
+ % The figure (if specified) must be in focus
+ if ~isempty(paramsStruct.hFigure) & ishandle(paramsStruct.hFigure) %#ok ML6
+ isFigureValid = 1; % no true available in ML6
+ figure(paramsStruct.hFigure);
+ else
+ isFigureValid = 0; % no false available in ML6
+ end
+
+ % Flush all graphic events to ensure correct rendering
+ drawnow; pause(0.01);
+
+ % No handle specified
+ wasPositionGiven = 1; % no true available in ML6
+ if isempty(paramsStruct.handle)
+
+ % Set default handle, if not supplied
+ paramsStruct.handle = paramsStruct.hFigure;
+
+ % If position was not specified, get it interactively using RBBOX
+ if isempty(paramsStruct.position)
+ [paramsStruct.position, jFrameUsed, msgStr] = getInteractivePosition(paramsStruct.hFigure); %#ok jFrameUsed is unused
+ paramsStruct.wasInteractive = 1; % no true available in ML6
+ wasPositionGiven = 0; % no false available in ML6
+ end
+
+ elseif ~ishandle(paramsStruct.handle)
+ % Handle was supplied - ensure it is a valid handle
+ error('YMA:screencapture:invalidHandle','Invalid handle passed to ScreenCapture');
+
+ elseif isempty(paramsStruct.position)
+ % Handle was supplied but position was not, so use the handle's position
+ paramsStruct.position = getPixelPos(paramsStruct.handle);
+ paramsStruct.position(1:2) = 0;
+ wasPositionGiven = 0; % no false available in ML6
+
+ elseif ~isnumeric(paramsStruct.position) | (length(paramsStruct.position) ~= 4) %#ok ML6
+ % Both handle & position were supplied - ensure a valid pixel position vector
+ error('YMA:screencapture:invalidPosition','Invalid position vector passed to ScreenCapture: \nMust be a [x,y,w,h] numeric pixel array');
+ end
+
+ % Capture current object (uicontrol/axes/figure) if w=h=0 (single-click in interactive mode)
+ if paramsStruct.position(3)<=0 | paramsStruct.position(4)<=0 %#ok ML6
+ %TODO - find a way to single-click another Matlab figure (the following does not work)
+ %paramsStruct.position = getPixelPos(ancestor(hittest,'figure'));
+ paramsStruct.position = getPixelPos(paramsStruct.handle);
+ paramsStruct.position(1:2) = 0;
+ paramsStruct.wasInteractive = 0; % no false available in ML6
+ wasPositionGiven = 0; % no false available in ML6
+ end
+
+ % First get the parent handle's desktop-based Matlab pixel position
+ parentPos = [0,0,0,0];
+ dX = 0;
+ dY = 0;
+ dW = 0;
+ dH = 0;
+ if ~isFigure(hParent)
+ % Get the reguested component's pixel position
+ parentPos = getPixelPos(hParent, 1); % no true available in ML6
+
+ % Axes position inaccuracy estimation
+ deltaX = 3;
+ deltaY = -1;
+
+ % Fix for images
+ if isImage(hParent) % | (isAxes(hParent) & strcmpi(get(hParent,'YDir'),'reverse')) %#ok ML6
+
+ % Compensate for resized image axes
+ hAxes = get(hParent,'Parent');
+ if all(get(hAxes,'DataAspectRatio')==1) % sanity check: this is the normal behavior
+ % Note 18/4/2013: the following fails for non-square images
+ %actualImgSize = min(parentPos(3:4));
+ %dX = (parentPos(3) - actualImgSize) / 2;
+ %dY = (parentPos(4) - actualImgSize) / 2;
+ %parentPos(3:4) = actualImgSize;
+
+ % The following should work for all types of images
+ actualImgSize = size(get(hParent,'CData'));
+ dX = (parentPos(3) - min(parentPos(3),actualImgSize(2))) / 2;
+ dY = (parentPos(4) - min(parentPos(4),actualImgSize(1))) / 2;
+ parentPos(3:4) = actualImgSize([2,1]);
+ %parentPos(3) = max(parentPos(3),actualImgSize(2));
+ %parentPos(4) = max(parentPos(4),actualImgSize(1));
+ end
+
+ % Fix user-specified img positions (but not auto-inferred ones)
+ if wasPositionGiven
+
+ % In images, use data units rather than pixel units
+ % Reverse the YDir
+ ymax = max(get(hParent,'YData'));
+ paramsStruct.position(2) = ymax - paramsStruct.position(2) - paramsStruct.position(4);
+
+ % Note: it would be best to use hgconvertunits, but:
+ % ^^^^ (1) it fails on Matlab 6, and (2) it doesn't accept Data units
+ %paramsStruct.position = hgconvertunits(hFig, paramsStruct.position, 'Data', 'pixel', hParent); % fails!
+ xLims = get(hParent,'XData');
+ yLims = get(hParent,'YData');
+ xPixelsPerData = parentPos(3) / (diff(xLims) + 1);
+ yPixelsPerData = parentPos(4) / (diff(yLims) + 1);
+ paramsStruct.position(1) = round((paramsStruct.position(1)-xLims(1)) * xPixelsPerData);
+ paramsStruct.position(2) = round((paramsStruct.position(2)-yLims(1)) * yPixelsPerData + 2*dY);
+ paramsStruct.position(3) = round( paramsStruct.position(3) * xPixelsPerData);
+ paramsStruct.position(4) = round( paramsStruct.position(4) * yPixelsPerData);
+
+ % Axes position inaccuracy estimation
+ if strcmpi(computer('arch'),'win64')
+ deltaX = 7;
+ deltaY = -7;
+ else
+ deltaX = 3;
+ deltaY = -3;
+ end
+
+ else % axes/image position was auto-infered (entire image)
+ % Axes position inaccuracy estimation
+ if strcmpi(computer('arch'),'win64')
+ deltaX = 6;
+ deltaY = -6;
+ else
+ deltaX = 2;
+ deltaY = -2;
+ end
+ dW = -2*dX;
+ dH = -2*dY;
+ end
+ end
+
+ %hFig = ancestor(hParent,'figure');
+ hParent = paramsStruct.hFigure;
+
+ elseif paramsStruct.wasInteractive % interactive figure rectangle
+
+ % Compensate for 1px rbbox inaccuracies
+ deltaX = 2;
+ deltaY = -2;
+
+ else % non-interactive figure
+
+ % Compensate 4px figure boundaries = difference betweeen OuterPosition and Position
+ deltaX = -1;
+ deltaY = 1;
+ end
+ %disp(paramsStruct.position) % for debugging
+
+ % Now get the pixel position relative to the monitor
+ figurePos = getPixelPos(hParent);
+ desktopPos = figurePos + parentPos;
+
+ % Now convert to Java-based pixels based on screen size
+ % Note: multiple monitors are automatically handled correctly, since all
+ % ^^^^ Java positions are relative to the main monitor's top-left corner
+ javaX = desktopPos(1) + paramsStruct.position(1) + deltaX + dX;
+ javaY = screenSize(4) - desktopPos(2) - paramsStruct.position(2) - paramsStruct.position(4) + deltaY + dY;
+ width = paramsStruct.position(3) + dW;
+ height = paramsStruct.position(4) + dH;
+ paramsStruct.position = round([javaX, javaY, width, height]);
+ %paramsStruct.position
+
+ % Ensure the figure is at the front so it can be screen-captured
+ if isFigureValid
+ figure(hParent);
+ drawnow;
+ pause(0.02);
+ end
+ catch
+ % Maybe root/desktop handle (root does not have a 'Position' prop so getPixelPos croaks
+ if isequal(double(hParent),0) % =root/desktop handle; handles case of hParent=[]
+ javaX = paramsStruct.position(1) - 1;
+ javaY = screenSize(4) - paramsStruct.position(2) - paramsStruct.position(4) - 1;
+ paramsStruct.position = [javaX, javaY, paramsStruct.position(3:4)];
+ end
+ end
+%end % convertPos
+
+%% Interactively get the requested capture rectangle
+function [positionRect, jFrameUsed, msgStr] = getInteractivePosition(hFig)
+ msgStr = '';
+ try
+ % First try the invisible-figure approach, in order to
+ % enable rbbox outside any existing figure boundaries
+ f = figure('units','pixel','pos',[-100,-100,10,10],'HitTest','off');
+ drawnow; pause(0.01);
+ oldWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame');
+ jf = get(handle(f),'JavaFrame');
+ warning(oldWarn);
+ try
+ jWindow = jf.fFigureClient.getWindow;
+ catch
+ try
+ jWindow = jf.fHG1Client.getWindow;
+ catch
+ jWindow = jf.getFigurePanelContainer.getParent.getTopLevelAncestor;
+ end
+ end
+ com.sun.awt.AWTUtilities.setWindowOpacity(jWindow,0.05); %=nearly transparent (not fully so that mouse clicks are captured)
+ jWindow.setMaximized(1); % no true available in ML6
+ jFrameUsed = 1; % no true available in ML6
+ msg = {'Mouse-click and drag a bounding rectangle for screen-capture ' ...
+ ... %'or single-click any Matlab figure to capture the entire figure.' ...
+ };
+ catch
+ % Something failed, so revert to a simple rbbox on a visible figure
+ try delete(f); drawnow; catch, end %Cleanup...
+ jFrameUsed = 0; % no false available in ML6
+ msg = {'Mouse-click within any Matlab figure and then', ...
+ 'drag a bounding rectangle for screen-capture,', ...
+ 'or single-click to capture the entire figure'};
+ end
+ uiwait(msgbox(msg,'ScreenCapture'));
+
+ k = waitforbuttonpress; %#ok k is unused
+ %hFig = getCurrentFig;
+ %p1 = get(hFig,'CurrentPoint');
+ positionRect = rbbox;
+ %p2 = get(hFig,'CurrentPoint');
+
+ if jFrameUsed
+ jFrameOrigin = getPixelPos(f);
+ delete(f); drawnow;
+ try
+ figOrigin = getPixelPos(hFig);
+ catch % empty/invalid hFig handle
+ figOrigin = [0,0,0,0];
+ end
+ else
+ if isempty(hFig)
+ jFrameOrigin = getPixelPos(gcf);
+ else
+ jFrameOrigin = [0,0,0,0];
+ end
+ figOrigin = [0,0,0,0];
+ end
+ positionRect(1:2) = positionRect(1:2) + jFrameOrigin(1:2) - figOrigin(1:2);
+
+ if prod(positionRect(3:4)) > 0
+ msgStr = sprintf('%dx%d area captured',positionRect(3),positionRect(4));
+ end
+%end % getInteractivePosition
+
+%% Get current figure (even if its handle is hidden)
+function hFig = getCurrentFig
+ oldState = get(0,'showHiddenHandles');
+ set(0,'showHiddenHandles','on');
+ hFig = get(0,'CurrentFigure');
+ set(0,'showHiddenHandles',oldState);
+%end % getCurrentFig
+
+%% Get ancestor figure - used for old Matlab versions that don't have a built-in ancestor()
+function hObj = ancestor(hObj,type)
+ if ~isempty(hObj) & ishandle(hObj) %#ok for Matlab 6 compatibility
+ try
+ hObj = get(hObj,'Ancestor');
+ catch
+ % never mind...
+ end
+ try
+ %if ~isa(handle(hObj),type) % this is best but always returns 0 in Matlab 6!
+ %if ~isprop(hObj,'type') | ~strcmpi(get(hObj,'type'),type) % no isprop() in ML6!
+ try
+ objType = get(hObj,'type');
+ catch
+ objType = '';
+ end
+ if ~strcmpi(objType,type)
+ try
+ parent = get(handle(hObj),'parent');
+ catch
+ parent = hObj.getParent; % some objs have no 'Parent' prop, just this method...
+ end
+ if ~isempty(parent) % empty parent means root ancestor, so exit
+ hObj = ancestor(parent,type);
+ end
+ end
+ catch
+ % never mind...
+ end
+ end
+%end % ancestor
+
+%% Get position of an HG object in specified units
+function pos = getPos(hObj,field,units)
+ % Matlab 6 did not have hgconvertunits so use the old way...
+ oldUnits = get(hObj,'units');
+ if strcmpi(oldUnits,units) % don't modify units unless we must!
+ pos = get(hObj,field);
+ else
+ set(hObj,'units',units);
+ pos = get(hObj,field);
+ set(hObj,'units',oldUnits);
+ end
+%end % getPos
+
+%% Get pixel position of an HG object - for Matlab 6 compatibility
+function pos = getPixelPos(hObj,varargin)
+ persistent originalObj
+ try
+ stk = dbstack;
+ if ~strcmp(stk(2).name,'getPixelPos')
+ originalObj = hObj;
+ end
+
+ if isFigure(hObj) %| isAxes(hObj)
+ %try
+ pos = getPos(hObj,'OuterPosition','pixels');
+ else %catch
+ % getpixelposition is unvectorized unfortunately!
+ pos = getpixelposition(hObj,varargin{:});
+
+ % add the axes labels/ticks if relevant (plus a tiny margin to fix 2px label/title inconsistencies)
+ if isAxes(hObj) & ~isImage(originalObj) %#ok ML6
+ tightInsets = getPos(hObj,'TightInset','pixel');
+ pos = pos + tightInsets.*[-1,-1,1,1] + [-1,1,1+tightInsets(1:2)];
+ end
+ end
+ catch
+ try
+ % Matlab 6 did not have getpixelposition nor hgconvertunits so use the old way...
+ pos = getPos(hObj,'Position','pixels');
+ catch
+ % Maybe the handle does not have a 'Position' prop (e.g., text/line/plot) - use its parent
+ pos = getPixelPos(get(hObj,'parent'),varargin{:});
+ end
+ end
+
+ % Handle the case of missing/invalid/empty HG handle
+ if isempty(pos)
+ pos = [0,0,0,0];
+ end
+%end % getPixelPos
+
+%% Adds a ScreenCapture toolbar button
+function addToolbarButton(paramsStruct)
+ % Ensure we have a valid toolbar handle
+ hFig = ancestor(paramsStruct.toolbar,'figure');
+ if isempty(hFig)
+ error('YMA:screencapture:badToolbar','the ''Toolbar'' parameter must contain a valid GUI handle');
+ end
+ set(hFig,'ToolBar','figure');
+ hToolbar = findall(hFig,'type','uitoolbar');
+ if isempty(hToolbar)
+ error('YMA:screencapture:noToolbar','the ''Toolbar'' parameter must contain a figure handle possessing a valid toolbar');
+ end
+ hToolbar = hToolbar(1); % just in case there are several toolbars... - use only the first
+
+ % Prepare the camera icon
+ icon = ['3333333333333333'; ...
+ '3333333333333333'; ...
+ '3333300000333333'; ...
+ '3333065556033333'; ...
+ '3000000000000033'; ...
+ '3022222222222033'; ...
+ '3022220002222033'; ...
+ '3022203110222033'; ...
+ '3022201110222033'; ...
+ '3022204440222033'; ...
+ '3022220002222033'; ...
+ '3022222222222033'; ...
+ '3000000000000033'; ...
+ '3333333333333333'; ...
+ '3333333333333333'; ...
+ '3333333333333333'];
+ cm = [ 0 0 0; ... % black
+ 0 0.60 1; ... % light blue
+ 0.53 0.53 0.53; ... % light gray
+ NaN NaN NaN; ... % transparent
+ 0 0.73 0; ... % light green
+ 0.27 0.27 0.27; ... % gray
+ 0.13 0.13 0.13]; % dark gray
+ cdata = ind2rgb(uint8(icon-'0'),cm);
+
+ % If the button does not already exit
+ hButton = findall(hToolbar,'Tag','ScreenCaptureButton');
+ tooltip = 'Screen capture';
+ if ~isempty(paramsStruct.target)
+ tooltip = [tooltip ' to ' paramsStruct.target];
+ end
+ if isempty(hButton)
+ % Add the button with the icon to the figure's toolbar
+ hButton = uipushtool(hToolbar, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']); %#ok unused
+ else
+ % Otherwise, simply update the existing button
+ set(hButton, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']);
+ end
+%end % addToolbarButton
+
+%% Java-get the actual screen-capture image data
+function imgData = getScreenCaptureImageData(positionRect)
+ if isempty(positionRect) | all(positionRect==0) | positionRect(3)<=0 | positionRect(4)<=0 %#ok ML6
+ imgData = [];
+ else
+ % Use java.awt.Robot to take a screen-capture of the specified screen area
+ rect = java.awt.Rectangle(positionRect(1), positionRect(2), positionRect(3), positionRect(4));
+ robot = java.awt.Robot;
+ jImage = robot.createScreenCapture(rect);
+
+ % Convert the resulting Java image to a Matlab image
+ % Adapted for a much-improved performance from:
+ % http://www.mathworks.com/support/solutions/data/1-2WPAYR.html
+ h = jImage.getHeight;
+ w = jImage.getWidth;
+ %imgData = zeros([h,w,3],'uint8');
+ %pixelsData = uint8(jImage.getData.getPixels(0,0,w,h,[]));
+ %for i = 1 : h
+ % base = (i-1)*w*3+1;
+ % imgData(i,1:w,:) = deal(reshape(pixelsData(base:(base+3*w-1)),3,w)');
+ %end
+
+ % Performance further improved based on feedback from Urs Schwartz:
+ %pixelsData = reshape(typecast(jImage.getData.getDataStorage,'uint32'),w,h).';
+ %imgData(:,:,3) = bitshift(bitand(pixelsData,256^1-1),-8*0);
+ %imgData(:,:,2) = bitshift(bitand(pixelsData,256^2-1),-8*1);
+ %imgData(:,:,1) = bitshift(bitand(pixelsData,256^3-1),-8*2);
+
+ % Performance even further improved based on feedback from Jan Simon:
+ pixelsData = reshape(typecast(jImage.getData.getDataStorage, 'uint8'), 4, w, h);
+ imgData = cat(3, ...
+ transpose(reshape(pixelsData(3, :, :), w, h)), ...
+ transpose(reshape(pixelsData(2, :, :), w, h)), ...
+ transpose(reshape(pixelsData(1, :, :), w, h)));
+ end
+%end % getInteractivePosition
+
+%% Return the figure to its pre-undocked state (when relevant)
+function redockFigureIfRelevant(paramsStruct)
+ if paramsStruct.wasDocked
+ try
+ set(paramsStruct.hFigure,'WindowStyle','docked');
+ %drawnow;
+ catch
+ % never mind - ignore...
+ end
+ end
+%end % redockFigureIfRelevant
+
+%% Copy screen-capture to the system clipboard
+% Adapted from http://www.mathworks.com/matlabcentral/fileexchange/28708-imclipboard/content/imclipboard.m
+function imclipboard(imgData)
+ % Import necessary Java classes
+ import java.awt.Toolkit.*
+ import java.awt.image.BufferedImage
+ import java.awt.datatransfer.DataFlavor
+
+ % Add the necessary Java class (ImageSelection) to the Java classpath
+ if ~exist('ImageSelection', 'class')
+ % Obtain the directory of the executable (or of the M-file if not deployed)
+ %javaaddpath(fileparts(which(mfilename)), '-end');
+ if isdeployed % Stand-alone mode.
+ [status, result] = system('path'); %#ok
+ MatLabFilePath = char(regexpi(result, 'Path=(.*?);', 'tokens', 'once'));
+ else % MATLAB mode.
+ MatLabFilePath = fileparts(mfilename('fullpath'));
+ end
+ javaaddpath(MatLabFilePath, '-end');
+ end
+
+ % Get System Clipboard object (java.awt.Toolkit)
+ cb = getDefaultToolkit.getSystemClipboard; % can't use () in ML6!
+
+ % Get image size
+ ht = size(imgData, 1);
+ wd = size(imgData, 2);
+
+ % Convert to Blue-Green-Red format
+ imgData = imgData(:, :, [3 2 1]);
+
+ % Convert to 3xWxH format
+ imgData = permute(imgData, [3, 2, 1]);
+
+ % Append Alpha data (not used)
+ imgData = cat(1, imgData, 255*ones(1, wd, ht, 'uint8'));
+
+ % Create image buffer
+ imBuffer = BufferedImage(wd, ht, BufferedImage.TYPE_INT_RGB);
+ imBuffer.setRGB(0, 0, wd, ht, typecast(imgData(:), 'int32'), 0, wd);
+
+ % Create ImageSelection object
+ % % custom java class
+ imSelection = ImageSelection(imBuffer);
+
+ % Set clipboard content to the image
+ cb.setContents(imSelection, []);
+%end %imclipboard
+
+%% Is the provided handle a figure?
+function flag = isFigure(hObj)
+ flag = isa(handle(hObj),'figure') | isa(hObj,'matlab.ui.Figure');
+%end %isFigure
+
+%% Is the provided handle an axes?
+function flag = isAxes(hObj)
+ flag = isa(handle(hObj),'axes') | isa(hObj,'matlab.graphics.axis.Axes');
+%end %isFigure
+
+%% Is the provided handle an image?
+function flag = isImage(hObj)
+ flag = isa(handle(hObj),'image') | isa(hObj,'matlab.graphics.primitive.Image');
+%end %isFigure
+
+%%%%%%%%%%%%%%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%%%%%%%%%%%
+% find a way in interactive-mode to single-click another Matlab figure for screen-capture
\ No newline at end of file
diff --git a/GPE Solver/+Plotter/Analysis.m b/GPE Solver/+Plotter/Analysis.m
new file mode 100644
index 0000000..df85ab4
--- /dev/null
+++ b/GPE Solver/+Plotter/Analysis.m
@@ -0,0 +1,50 @@
+ set(0,'defaulttextInterpreter','latex')
+ set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
+ format long
+
+ runIdx = 6;
+
+ load(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
+
+ x = Transf.x*Params.l0*1e6;
+ y = Transf.y*Params.l0*1e6;
+ z = Transf.z*Params.l0*1e6;
+ %percentcomplete = linspace(0,1,Params.cut_off/200);
+
+ dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
+ %Plotting
+ subplot(2,3,1)
+ n = abs(psi).^2;
+ nxz = squeeze(trapz(n*dy,2));
+ nyz = squeeze(trapz(n*dx,1));
+ nxy = squeeze(trapz(n*dz,3));
+
+ plotxz = pcolor(x,z,nxz');
+ set(plotxz, 'EdgeColor', 'none');
+ xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
+
+ subplot(2,3,2)
+ plotyz = pcolor(y,z,nyz');
+ set(plotyz, 'EdgeColor', 'none');
+ xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
+
+ subplot(2,3,3)
+ plotxy = pcolor(x,y,nxy');
+ set(plotxy, 'EdgeColor', 'none');
+ xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
+
+ subplot(2,3,4)
+ plot(-log10(Observ.residual),'-b')
+ ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
+
+ subplot(2,3,5)
+ plot(Observ.EVec,'-b')
+ ylabel('$E$'); xlabel('steps');
+
+ subplot(2,3,6)
+ plot(Observ.mucVec,'-b')
+ ylabel('$\mu$'); xlabel('steps');
+% xlim([0,1]); ylim([0,8]);
+% xlim([0,1]); ylim([0,8]);
+
+ Ecomp = energy_components(psi,Params,Transf,VDk,V);
diff --git a/GPE Solver/+Plotter/MakeMovie.m b/GPE Solver/+Plotter/MakeMovie.m
new file mode 100644
index 0000000..8842cbd
--- /dev/null
+++ b/GPE Solver/+Plotter/MakeMovie.m
@@ -0,0 +1,77 @@
+set(0,'defaulttextInterpreter','latex')
+set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
+
+RunIdx = 1;
+
+FileDir = dir(sprintf('./Data/Run_%03i/TimeEvolution/*.mat',RunIdx));
+NumFiles = numel(FileDir);
+QuenchSettings = load(sprintf('./Data/Run_%03i/QuenchSettings',RunIdx),'Quench','Params','Transf','VDk','V');
+Transf = QuenchSettings.Transf; Params = QuenchSettings.Params;
+x = Transf.x; y = Transf.y; z = Transf.z;
+dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
+
+mkdir(sprintf('./Data/Run_%03i/Figures',RunIdx))
+outputVideo = VideoWriter(fullfile('./Data/Movie.avi'));
+outputVideo.FrameRate = 10;
+open(outputVideo)
+
+figure(1);
+x0 = 800;
+y0 = 200;
+width = 800;
+height = 600;
+set(gcf,'position',[x0,y0,width,height])
+
+EVecTemp = [];
+
+for ii = 2:(NumFiles-1)
+ load(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',RunIdx,ii),'psi','muchem','T','Observ','t_idx');
+
+ %Plotting
+ subplot(2,3,1)
+ n = abs(psi).^2;
+ nxz = squeeze(trapz(n*dy,2));
+ nyz = squeeze(trapz(n*dx,1));
+ nxy = squeeze(trapz(n*dz,3));
+
+ plotxz = pcolor(x,z,nxz'); shading interp
+ set(plotxz, 'EdgeColor', 'none');
+ xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
+
+ subplot(2,3,2)
+ plotyz = pcolor(y,z,nyz'); shading interp
+ set(plotyz, 'EdgeColor', 'none');
+ xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
+
+ subplot(2,3,3)
+ plotxy = pcolor(x,y,nxy'); shading interp
+ set(plotxy, 'EdgeColor', 'none');
+ xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
+
+ subplot(2,3,4)
+ plot(Observ.tVecPlot*1000/Params.w0,Observ.NormVec,'-b')
+ ylabel('Normalization'); xlabel('$t$ [$m$s]');
+
+ subplot(2,3,5)
+ plot(Observ.tVecPlot*1000/Params.w0,1-2*Observ.PCVec/pi,'-b')
+ ylabel('Coherence'); xlabel('$t$ [$m$s]');
+ ylim([0,1])
+
+ subplot(2,3,6)
+ plot(Observ.tVecPlot*1000/Params.w0,Observ.EVec,'-b')
+ ylabel('E'); xlabel('$t$ [$m$s]');
+
+ tVal = Observ.tVecPlot(end)*1000/Params.w0;
+ sgtitle(sprintf('$\\mu =%.3f \\hbar\\omega_0$, $T=%.1f$nK, $t=%.1f$ms',muchem,T,tVal))
+
+ drawnow
+ saveas(gcf,sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii))
+ img = imread(sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii));
+ writeVideo(outputVideo,img)
+% hold off;
+ clf
+end
+
+close(outputVideo)
+close(figure(1))
+delete(sprintf('./Data/Run_%03i/Figures/*.jpg',RunIdx)) % deleting images after movie is made
\ No newline at end of file
diff --git a/GPE Solver/+Plotter/OrderParameter_m.m b/GPE Solver/+Plotter/OrderParameter_m.m
new file mode 100644
index 0000000..2ce76dc
--- /dev/null
+++ b/GPE Solver/+Plotter/OrderParameter_m.m
@@ -0,0 +1,58 @@
+function [m_Order] = OrderParameter_m(psi,Transf,Params,VDk,V,T,muchem)
+
+ NumRealiz = 100;
+
+ Mx = numel(Transf.x);
+ My = numel(Transf.y);
+ Mz = numel(Transf.z);
+
+ r = normrnd(0,1,size(psi));
+ theta = rand(size(psi));
+ noise = r.*exp(2*pi*1i*theta);
+
+ KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
+ Gamma = 1-1i*Params.gamma_S;
+ dt = Params.dt;
+
+ avgpsi = 0;
+ avgpsi2 = 0;
+
+ for jj = 1:NumRealiz
+ %generate initial state
+ xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
+ swapx = randi(length(Transf.x),1,length(Transf.x));
+ swapy = randi(length(Transf.y),1,length(Transf.y));
+ swapz = randi(length(Transf.z),1,length(Transf.z));
+ psi_j = psi + xi * noise(swapx,swapy,swapz);
+
+ % --- % propagate forward in time 1 time step:
+ %kin
+ psi_j = fftn(psi_j);
+ psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
+ psi_j = ifftn(psi_j);
+
+ %DDI
+ frho = fftn(abs(psi_j).^2);
+ Phi = real(ifftn(frho.*VDk));
+
+ %Real-space
+ psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
+
+ %kin
+ psi_j = fftn(psi_j);
+ psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
+ psi_j = ifftn(psi_j);
+
+ %Projection
+ kcut = sqrt(2*Params.e_cut);
+ K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2) depends on how Vdk (DDI) is defined
+
+%Trap gamma
+Params.gx=(Params.wx/w0)^2;
+Params.gy=(Params.wy/w0)^2;
+Params.gz=(Params.wz/w0)^2;
+
+%Loading the rest into Params
+Params.hbar = hbar; Params.kbol = kbol; Params.mu0 = mu0; Params.muB = muB; Params.a0 = a0;
+Params.w0 = w0; Params.l0 = l0;
\ No newline at end of file
diff --git a/GPE Solver/+Simulator/running_file.m b/GPE Solver/+Simulator/running_file.m
new file mode 100644
index 0000000..f980611
--- /dev/null
+++ b/GPE Solver/+Simulator/running_file.m
@@ -0,0 +1,23 @@
+%-% Running file %-%
+clearvars
+
+njob = 6;
+
+mkdir(sprintf('./Data'))
+mkdir(sprintf('./Data/Run_%03i',njob))
+
+%Obtain simulation parameters
+[Params] = parameters();
+
+%Set up spatial grids and transforms
+[Transf] = setup_space(Params);
+
+[psi,V,VDk] = Initialize(Params,Transf);
+
+% --- Initialize
+Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
+t_idx = 1; %Start at t = 0;
+Observ.res_idx = 1;
+
+
+[psi] = ssfm_imag(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
\ No newline at end of file
diff --git a/GPE Solver/+Simulator/setup_space.m b/GPE Solver/+Simulator/setup_space.m
new file mode 100644
index 0000000..9929e31
--- /dev/null
+++ b/GPE Solver/+Simulator/setup_space.m
@@ -0,0 +1,32 @@
+function [Transf] = setup_space(Params)
+Transf.Xmax = 0.5*Params.Lx;
+Transf.Ymax = 0.5*Params.Ly;
+Transf.Zmax = 0.5*Params.Lz;
+
+Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
+
+% Fourier grids
+x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
+Kmax = pi*Params.Nx/Params.Lx;
+kx = linspace(-Kmax,Kmax,Nx+1);
+kx = kx(1:end-1); dkx = kx(2)-kx(1);
+kx = fftshift(kx);
+
+y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
+Kmax = pi*Params.Ny/Params.Ly;
+ky = linspace(-Kmax,Kmax,Ny+1);
+ky = ky(1:end-1); dky = ky(2)-ky(1);
+ky = fftshift(ky);
+
+z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
+Kmax = pi*Params.Nz/Params.Lz;
+kz = linspace(-Kmax,Kmax,Nz+1);
+kz = kz(1:end-1); dkz = kz(2)-kz(1);
+kz = fftshift(kz);
+
+[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
+[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
+Transf.x = x; Transf.y = y; Transf.z = z;
+Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
+Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
+Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
\ No newline at end of file
diff --git a/GPE Solver/+Simulator/setup_space_radial.m b/GPE Solver/+Simulator/setup_space_radial.m
new file mode 100644
index 0000000..3887b9e
--- /dev/null
+++ b/GPE Solver/+Simulator/setup_space_radial.m
@@ -0,0 +1,311 @@
+function [Transf] = setup_space_radial(Params,morder)
+Zmax = 0.5*Params.Lz;
+Rmax = Params.Lr;
+Nz = Params.Nz;
+Nr = Params.Nr;
+
+if(nargin==1)
+ morder=0; %only do Bessel J0
+end
+
+% Fourier grids
+z=linspace(-Zmax,Zmax,Nz+1);
+z=z(1:end-1);
+dz=z(2)-z(1);
+Kmax=Nz*2*pi/(4*Zmax);
+kz=linspace(-Kmax,Kmax,Nz+1);
+kz=kz(1:end-1);
+
+% Hankel grids and transform
+H = hankelmatrix(morder,Rmax,Nr);
+r=H.r(:);
+kr=H.kr(:);
+T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
+Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
+wr=H.wr;
+wk=H.wk;
+% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
+
+[Transf.R,Transf.Z]=ndgrid(r,z);
+[Transf.KR,Transf.KZ]=ndgrid(kr,kz);
+Transf.T=T;
+Transf.Tinv=Tinv;
+Transf.r=r;
+Transf.kr=kr;
+Transf.z=z;
+Transf.kz=kz;
+Transf.wr=wr;
+Transf.wk=wk;
+Transf.Rmax=Rmax;
+Transf.Zmax=Zmax;
+Transf.dz=z(2)-z(1);
+Transf.dkz=kz(2)-kz(1);
+%b1=Transf;
+
+function s_HT = hankelmatrix(order, rmax, Nr, eps_roots)
+%HANKEL_MATRIX: Generates data to use for Hankel Transforms
+%
+% s_HT = hankel_matrix(order, rmax, Nr, eps_roots)
+%
+% s_HT = Structure containing data to use for the pQDHT
+% order = Transform order
+% rmax = Radial extent of transform
+% Nr = Number of sample points
+% eps_roots = Error in estimation of roots of Bessel function (optional)
+%
+% s_HT:
+% order, rmax, Nr = As above
+% J_roots = Roots of the pth order Bessel fn.
+% J_roots_N1 = (N+1)th root
+% r = Radial co-ordinate vector
+% v = frequency co-ordinate vector
+% kr = Radial wave number co-ordinate vector
+% vmax = Limiting frequency
+% = roots_N1 / (2*pi*rmax)
+% S = rmax * 2*pi*vmax (S product)
+% T = Transform matrix
+% J = Scaling vector
+% = J_(order+1){roots}
+%
+% The algorithm used is that from:
+% "Computation of quasi-discrete Hankel transforms of the integer
+% order for propagating optical wave fields"
+% Manuel Guizar-Sicairos and Julio C. Guitierrez-Vega
+% J. Opt. Soc. Am. A 21(1) 53-58 (2004)
+%
+% The algorithm also calls the function:
+% zn = bessel_zeros(1, p, Nr+1, 1e-6),
+% where p and N are defined above, to calculate the roots of the bessel
+% function. This algorithm is taken from:
+% "An Algorithm with ALGOL 60 Program for the Computation of the
+% zeros of the Ordinary Bessel Functions and those of their
+% Derivatives".
+% N. M. Temme
+% Journal of Computational Physics, 32, 270-279 (1979)
+%
+% Example: Propagation of radial field
+%
+% % Note the use of matrix and element products / divisions
+% H = hankel_matrix(0, 1e-3, 512);
+% DR0 = 50e-6;
+% Ur0 = exp(-(H.r/DR0).^2);
+% Ukr0 = H.T * (Ur0./H.J);
+% k0 = 2*pi/800e-9;
+% kz = realsqrt((k0^2 - H.kr.^2).*(k0>H.kr));
+% z = (-5e-3:1e-5:5e-3);
+% Ukrz = (Ukr0*ones(1,length(z))).*exp(i*kz*z);
+% Urz = (H.T * Ukrz) .* (H.J * ones(1,length(z)));
+%
+% See also bessel_zeros, besselj
+
+if (~exist('eps_roots', 'var')||isemtpy(eps_roots))
+ s_HT.eps_roots = 1e-6;
+else
+ s_HT.eps_roots = eps_roots;
+end
+
+s_HT.order = order;
+s_HT.rmax = rmax;
+s_HT.Nr = Nr;
+
+% Calculate N+1 roots:
+J_roots = bessel_zeros(1, s_HT.order, s_HT.Nr+1, s_HT.eps_roots);
+s_HT.J_roots = J_roots(1:end-1);
+s_HT.J_roots_N1 = J_roots(end);
+
+% Calculate co-ordinate vectors
+s_HT.r = s_HT.J_roots * s_HT.rmax / s_HT.J_roots_N1;
+s_HT.v = s_HT.J_roots / (2*pi * s_HT.rmax);
+s_HT.kr = 2*pi * s_HT.v;
+s_HT.kmax = s_HT.J_roots_N1 / (s_HT.rmax);
+s_HT.vmax = s_HT.J_roots_N1 / (2*pi * s_HT.rmax);
+s_HT.S = s_HT.J_roots_N1;
+
+% Calculate hankel matrix and vectors
+% I use (p=order) and (p1=order+1)
+Jp = besselj(s_HT.order, (s_HT.J_roots) * (s_HT.J_roots.') / s_HT.S);
+Jp1 = abs(besselj(s_HT.order+1, s_HT.J_roots));
+s_HT.T = 2*Jp./(Jp1 * (Jp1.') * s_HT.S);
+s_HT.J = Jp1;
+s_HT.wr=2./((s_HT.kmax)^2*abs(Jp1).^2);
+s_HT.wk=2./((s_HT.rmax)^2*abs(Jp1).^2);
+
+return
+
+
+
+
+function z = bessel_zeros(d, a, n, e)
+%BESSEL_ZEROS: Finds the first n zeros of a bessel function
+%
+% z = bessel_zeros(d, a, n, e)
+%
+% z = zeros of the bessel function
+% d = Bessel function type:
+% 1: Ja
+% 2: Ya
+% 3: Ja'
+% 4: Ya'
+% a = Bessel order (a>=0)
+% n = Number of zeros to find
+% e = Relative error in root
+%
+% This function uses the routine described in:
+% "An Algorithm with ALGOL 60 Program for the Computation of the
+% zeros of the Ordinary Bessel Functions and those of their
+% Derivatives".
+% N. M. Temme
+% Journal of Computational Physics, 32, 270-279 (1979)
+
+z = zeros(n, 1);
+aa = a^2;
+mu = 4*aa;
+mu2 = mu^2;
+mu3 = mu^3;
+mu4 = mu^4;
+
+if (d<3)
+ p = 7*mu - 31;
+ p0 = mu - 1;
+ if ((1+p)==p)
+ p1 = 0;
+ q1 = 0;
+ else
+ p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p);
+ q1 = 1.6*(83*mu2 - 982*mu + 3779)/p;
+ end
+else
+ p = 7*mu2 + 82*mu - 9;
+ p0 = mu + 3;
+ if ((p+1)==1)
+ p1 = 0;
+ q1 = 0;
+ else
+ p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p);
+ q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p;
+ end
+end
+
+if (d==1)|(d==4)
+ t = .25;
+else
+ t = .75;
+end
+tt = 4*t;
+
+if (d<3)
+ pp1 = 5/48;
+ qq1 = -5/36;
+else
+ pp1 = -7/48;
+ qq1 = 35/288;
+end
+
+y = .375*pi;
+if (a>=3)
+ bb = a^(-2/3);
+else
+ bb = 1;
+end
+a1 = 3*a - 8;
+% psi = (.5*a + .25)*pi;
+
+for s=1:n
+ if ((a==0)&(s==1)&(d==3))
+ x = 0;
+ j = 0;
+ else
+ if (s>=a1)
+ b = (s + .5*a - t)*pi;
+ c = .015625/(b^2);
+ x = b - .125*(p0 - p1*c)/(b*(1 - q1*c));
+ else
+ if (s==1)
+ switch (d)
+ case (1)
+ x = -2.33811;
+ case (2)
+ x = -1.17371;
+ case (3)
+ x = -1.01879;
+ otherwise
+ x = -2.29444;
+ end
+ else
+ x = y*(4*s - tt);
+ v = x^(-2);
+ x = -x^(2/3) * (1 + v*(pp1 + qq1*v));
+ end
+ u = x*bb;
+ v = fi(2/3 * (-u)^1.5);
+ w = 1/cos(v);
+ xx = 1 - w^2;
+ c = sqrt(u/xx);
+ if (d<3)
+ x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u));
+ else
+ x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u));
+ end
+ end
+ j = 0;
+
+while ((j==0)|((j<5)&(abs(w/x)>e)))
+ xx = x^2;
+ x4 = x^4;
+ a2 = aa - xx;
+ r0 = bessr(d, a, x);
+ j = j+1;
+ if (d<3)
+ u = r0;
+ w = 6*x*(2*a + 1);
+ p = (1 - 4*a2)/w;
+ q = (4*(xx-mu) - 2 - 12*a)/w;
+ else
+ u = -xx*r0/a2;
+ v = 2*x*a2/(3*(aa+xx));
+ w = 64*a2^3;
+ q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w;
+ p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w);
+ end
+ w = u*(1 + p*r0)/(1 + q*r0);
+ x = x+w;
+ end
+ z(s) = x;
+ end
+end
+
+
+function FI = fi(y)
+ c1 = 1.570796;
+ if (~y)
+ FI = 0;
+ elseif (y>1e5)
+ FI = c1;
+ else
+ if (y<1)
+ p = (3*y)^(1/3);
+ pp = p^2;
+ p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575);
+ else
+ p = 1/(y + c1);
+ pp = p^2;
+ p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465);
+ end
+ pp = (y+p)^2;
+ r = (p - atan(p+y))/pp;
+ FI = p - (1+pp)*r*(1 + r/(p+y));
+ end
+return
+
+function Jr = bessr(d, a, x)
+ switch (d)
+ case (1)
+ Jr = besselj(a, x)./besselj(a+1, x);
+ case (2)
+ Jr = bessely(a, x)./bessely(a+1, x);
+ case (3)
+ Jr = a./x - besselj(a+1, x)./besselj(a, x);
+ otherwise
+ Jr = a./x - bessely(a+1, x)./bessely(a, x);
+ end
+return
\ No newline at end of file
diff --git a/GPE Solver/+Simulator/ssfm_imag.m b/GPE Solver/+Simulator/ssfm_imag.m
new file mode 100644
index 0000000..dc84401
--- /dev/null
+++ b/GPE Solver/+Simulator/ssfm_imag.m
@@ -0,0 +1,98 @@
+function [psi] = ssfm_imag(psi,Params,Transf,VDk,V,njob,t_idx,Observ)
+
+set(0,'defaulttextInterpreter','latex')
+set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
+
+dt=-1j*abs(Params.dt);
+
+KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
+Observ.residual = 1; Observ.res = 1;
+
+figure(1)
+muchem = chemicalpotential(psi,Params,Transf,VDk,V);
+runningplot(psi,Params,Transf,Observ)
+drawnow
+
+AdaptIdx = 0;
+
+while t_idx < Params.cut_off
+ %kin
+ psi = fftn(psi);
+ psi = psi.*exp(-0.5*1i*dt*KEop);
+ psi = ifftn(psi);
+
+ %DDI
+ frho = fftn(abs(psi).^2);
+ Phi = real(ifftn(frho.*VDk));
+
+ %Real-space
+ psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
+
+ %kin
+ psi = fftn(psi);
+ psi = psi.*exp(-0.5*1i*dt*KEop);
+ psi = ifftn(psi);
+
+ %Renorm
+ Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
+ psi = sqrt(Params.N)*psi/sqrt(Norm);
+
+ muchem = chemicalpotential(psi,Params,Transf,VDk,V);
+
+ %Plotting loop
+ if mod(t_idx,1000) == 0
+
+ %Change in Energy
+ E = energytotal(psi,Params,Transf,VDk,V);
+ E = E/Norm;
+ Observ.EVec = [Observ.EVec E];
+
+ %Chemical potential
+ Observ.mucVec = [Observ.mucVec muchem];
+
+ %Normalized residuals
+ res = norm_resid(psi,Params,Transf,VDk,V,muchem);
+ Observ.residual = [Observ.residual res];
+
+ Observ.res_idx = Observ.res_idx + 1;
+ figure(1)
+ runningplot(psi,Params,Transf,Observ)
+ drawnow
+
+ save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
+
+ %Adaptive time step -- Careful, this can quickly get out of control
+ relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
+ if relres <1e-5
+ if AdaptIdx > 4 && abs(dt) > Params.mindt
+ dt = dt / 2;
+ fprintf('Time step changed to '); disp(dt);
+ AdaptIdx = 0;
+ elseif AdaptIdx > 4 && abs(dt) < Params.mindt
+ break
+ else
+ AdaptIdx = AdaptIdx + 1;
+ end
+ else
+ AdaptIdx = 0;
+ end
+ end
+ if any(isnan(psi(:)))
+ disp('Idiot.')
+ break
+ end
+ t_idx=t_idx+1;
+end
+
+
+%Change in Energy
+E = energytotal(psi,Params,Transf,VDk,V)
+E = E/Norm;
+Observ.EVec = [Observ.EVec E];
+
+% Phase coherence
+[PhaseC] = PhaseCoherence(psi,Transf,Params);
+Observ.PCVec = [Observ.PCVec PhaseC];
+
+Observ.res_idx = Observ.res_idx + 1;
+save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
\ No newline at end of file
diff --git a/GaussianBeamABCD/BeamPropagation.py b/GaussianBeamABCD/BeamPropagation.py
new file mode 100644
index 0000000..2dfbb46
--- /dev/null
+++ b/GaussianBeamABCD/BeamPropagation.py
@@ -0,0 +1,233 @@
+"""
+@author: Adam Newton Wright, https://github.com/adamnewtonwright/GaussianBeamPropagation
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+import sympy as sym
+from sympy import oo
+
+
+# Input Ray parameter, i.e. height and angle
+def ray(y,theta):
+ '''
+ Parameters
+ ----------
+ y : float or integer or sympy symbol in meters
+ The vertical height of a ray.
+ theta : float or integer in radians
+ The angle of divergence of the ray.
+
+ Returns
+ -------
+ mat : 2x1 matrix
+ [
+ [y],
+ [teta]
+ ]
+
+ '''
+
+ mat = np.array([[y],[theta]])
+ return mat
+
+# Ray Transfer Matrix for ideal lens with focal length f
+def lens(f):
+ '''
+ Parameters
+ ----------
+ f : float or integer or sympy symbol in meters
+ Thin lens focal length in meters
+
+ Returns
+ -------
+ mat : 2x2 matrix
+ [
+ [ 1, 0],
+ [-1/f, 1]
+ ]
+
+ '''
+
+ mat = np.array([[1,0], [-1/f, 1]])
+ return mat
+
+# Ray Transfer Matrix for propagation of distance d
+def prop(d):
+ '''
+ Parameters
+ ----------
+ d : float or integer or sympy symbol
+ Distance light is propagating along the z-axis.
+
+ Returns
+ -------
+ mat: 2x2 matrix
+ [
+ [1, d],
+ [0, 1]
+ ]
+
+ '''
+ mat = np.array([[1,d], [0,1]])
+ return mat
+
+# multiplying the matrices together. mat1 is the last matrix the light interacts with
+def mult(mat1,*argv):
+ '''
+ Parameters
+ ----------
+ mat1 : 2x2 ABCD matrix
+ Last matrix light interacts with.
+ *argv : 2x2 ABCD matrices
+ From left to right, the matrices should be entered such that the leftmost matrix interacts
+ with light temporally after the rightmost matrix.
+
+ Returns
+ -------
+ Mat : 2x2 matrix
+ The ABCd matrix describing the whole optical system.
+
+ '''
+
+ Mat = mat1
+ for arg in argv:
+ Mat = np.dot(Mat, arg)
+ return Mat
+
+# Adding Gaussian beam parameters
+def Zr(wo, lam):
+ '''
+ Parameters
+ ----------
+ wo : float, integer, or symbol
+ Beam waist radius in meters.
+ lam : float, integer, or symbol
+ Wavelength of light in meters.
+
+ Returns
+ -------
+ zr : float, int, symbols
+ Rayleigh range for given beam waist and wavelength.
+
+ '''
+
+ zr = np.pi * wo**2 / lam
+ return zr
+
+def W0(zr, lam):
+ '''
+ Parameters
+ ----------
+ zr : float, integer, symbol
+ Rayleigh range in meters
+ lam : float, integer, symbol
+ Wavelength of light in meters
+
+ Returns
+ -------
+ w0 : float, integer, symbol
+ Beam waist radius in meters
+
+ '''
+
+ w0 = np.sqrt(lam * zr / np.pi)
+ return w0
+
+# Remember, there should be an i in front of zr
+# but this complicates the calculations, so we usually just let z = 0
+# and don't explicitly deal with the i, but still do the math accordingly
+#def q0_func(z,zr):
+# qz = z + zr
+# return qz
+
+def q1_func(z, w0, lam, mat):
+ '''
+ Parameters
+ ----------
+ z : float, int, symbol
+ Position of the beam waist in meters.
+ w0 : float, int, symbol
+ Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
+ lam : float, int, symbol
+ Wavelength of light in meters.
+ mat : float, int, symbol
+ The ABCD 2x2 matrix describing the optical system.
+
+ Returns
+ -------
+ z: float, int, symbol
+ Position of the beam waist after the optical system
+ zr: float, int, symbol
+ Rayleigh range of the beam after the optical system
+ '''
+
+ A = mat[0][0]
+ B = mat[0][1]
+ C = mat[1][0]
+ D = mat[1][1]
+ zr = Zr(w0, lam)
+ real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
+ imag = (zr * (A*D - B*C)) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
+ z = real
+ zr = imag
+ return z, zr
+
+def q1_inv_func(z, w0, lam, mat):
+ '''
+ Parameters
+ ----------
+ z : float, int, symbol
+ Position of the beam waist in meters.
+ w0 : float, int, symbol
+ Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
+ lam : float, int, symbol
+ Wavelength of light in meters.
+ mat : float, int, symbol
+ The ABCD 2x2 matrix describing the optical system.
+
+ Returns
+ -------
+ R : float, int, symbol
+ Radius of curvature of the wavefront in meters.
+ w : float, int, symbol
+ Radius of the beam in meters.
+
+ '''
+ A = mat[0][0]
+ B = mat[0][1]
+ C = mat[1,0]
+ D = mat[1][1]
+ zr = Zr(w0, lam)
+ real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (A**2*(z**2 + zr**2) + 2*A*B*z + B**2)
+ imag = -zr * (A*D-B*C) / (A**2 *(z**2 + zr**2) + 2*A*B*z + B**2)
+ R = 1/real
+ w = (-lam / imag / np.pi)**.5
+ return R, w
+
+def plot(func, var, rang = np.arange(0,3,.01)):
+ '''
+ Parameters
+ ----------
+ func : Sympy function of one variable
+ Sympy function defining the beam width after the last optical element.
+ var : sympy variable
+ Variable in func that will be plotted.
+ rang : numpy array
+ Array of the values along the optical axis to be plotted
+
+ Returns
+ -------
+ plot : matplotlib graph
+ Graph of the beam width of var
+
+
+ '''
+ func = sym.lambdify(var, func)
+ plt.figure()
+ plt.plot(rang, func(rang), color = 'b')
+ plt.plot(rang, -func(rang), color = 'b')
+ plt.grid()
+ plt.xlabel('Optic Axis (m)')
+ plt.ylabel('Beam size (m)')
+ plt.show()
\ No newline at end of file
diff --git a/GaussianBeamABCD/DMD_Setup.py b/GaussianBeamABCD/DMD_Setup.py
new file mode 100644
index 0000000..b752938
--- /dev/null
+++ b/GaussianBeamABCD/DMD_Setup.py
@@ -0,0 +1,70 @@
+import BeamPropagation as bp # This is the script that handles the propagation
+import sympy as sym # For Symbolic examples
+import numpy as np # Handling of lists and for plotting
+import matplotlib.pyplot as plt # Plotting
+
+
+"""A Gaussian beam can be defined by it's (radial) waist wo, it's Rayleigh range zR, and the location of its waist zO"""
+w0 = 5.2E-3 # 1mm beam waist
+lam = 532E-9 # wavelength of 355 nm (UV)
+zR = bp.Zr(w0, lam) # Rayleigh range in m
+z0 = 0 # location of waist in m
+
+"""Define first 4-f optical system using matrices"""
+d1, d2, d3, f1, f2 = sym.symbols('d1 d2 d3 f1 f2')
+M = bp.mult(bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
+
+"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
+R, w = bp.q1_inv_func(0, w0, lam, M)
+
+"""Substitute and extract the required separation between lenses of first 4-f system"""
+distance_between_dmd_first_lens = 250E-3
+first_focal_length = 250.9E-3
+second_focal_length = 50E-3
+demag = 1/5
+target_w0 = demag*w0
+w = w.subs(f1, first_focal_length).subs(f2, second_focal_length).subs(d1, distance_between_dmd_first_lens).subs(d3,0)
+eq = sym.solve(w - target_w0, d2)[0]
+distance_between_lens_of_first_4f = list(eq.atoms())[0]
+print('Distance between lenses of first 4-f system = {} mm'.format(distance_between_lens_of_first_4f * 1E3))
+
+# Sanity check
+# expansion_factor = w.subs(d2,distance_between_lens_of_first_4f).subs(d3,0)/ w0
+# print('beam is w = {:.2f} x w0'.format(expansion_factor))
+
+# """Plot beam propagation up to 3 m after the first 4-f system"""
+# M = bp.mult(bp.prop(d3),bp.lens(second_focal_length),bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
+# R, w = bp.q1_inv_func(0, w0, lam, M)
+# bp.plot(w,d3, rang = np.arange(0,0.050,.0005))
+
+
+"""Define the full optical system of two 4-f setups using matrices"""
+d1, d2, d3, d4, f1, f2, f3 = sym.symbols('d1 d2 d3 d4 f1 f2 f3')
+M = bp.mult(bp.prop(d4),bp.lens(f3), bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
+
+"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
+R, w = bp.q1_inv_func(0, w0, lam, M)
+
+# """Find the focal length of lens required after the first 4-f system to have a collimated beam, given a certain separation between the first 4-f system and this lens"""
+distance_between_4fs = 550E-3
+R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d3,distance_between_4fs).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length)
+f3_coll = sym.solve(1/R_coll,f3)[0]
+third_focal_length = list(f3_coll.atoms())[0]
+print('For a fixed separation between first 4-f and third lens of {:.3f} mm, f3 = {:.3f} mm for a collimated beam'.format(distance_between_4fs* 1E3, third_focal_length * 1E3))
+
+# # """Plot beam propagation up to 3 m after the first 4-f system"""
+# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
+# R, w = bp.q1_inv_func(0, w0, lam, M)
+# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
+
+third_focal_length = 501.8E-3
+R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length).subs(f3,third_focal_length)
+d3_coll = sym.solve(1/R_coll,d3)[1]
+distance_between_4fs = list(d3_coll.atoms())[0]
+print('For a fixed third focal length of {:.3f} mm, d3 = {:.3f} mm, for a collimated beam'.format(third_focal_length* 1E3, distance_between_4fs * 1E3))
+
+# """Plot beam propagation up to 3 m after the first 4-f system"""
+# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
+# R, w = bp.q1_inv_func(0, w0, lam, M)
+# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
+
diff --git a/GaussianBeamABCD/Example.ipynb b/GaussianBeamABCD/Example.ipynb
new file mode 100644
index 0000000..6ea246f
--- /dev/null
+++ b/GaussianBeamABCD/Example.ipynb
@@ -0,0 +1,589 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Gaussian Beam Propagation"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Import files"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import BeamPropagation as bs # This is the script that handles the propagation\n",
+ "import sympy as sym # For Symbolic examples\n",
+ "import numpy as np # Handling of lists and for plotting\n",
+ "import matplotlib.pyplot as plt # Plotting"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Let's show what BeamProp_Script has"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Help on module BeamProp_Script:\n",
+ "\n",
+ "NAME\n",
+ " BeamProp_Script - Created on Wed Feb 19 15:51:54 2020\n",
+ "\n",
+ "DESCRIPTION\n",
+ " @author: wrighta\n",
+ "\n",
+ "FUNCTIONS\n",
+ " W0(zr, lam)\n",
+ " Parameters\n",
+ " ----------\n",
+ " zr : float, integer, symbol\n",
+ " Rayleigh range in meters\n",
+ " lam : float, integer, symbol\n",
+ " Wavelength of light in meters\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " w0 : float, integer, symbol\n",
+ " Beam waist radius in meters\n",
+ " \n",
+ " Zr(wo, lam)\n",
+ " Parameters\n",
+ " ----------\n",
+ " wo : float, integer, or symbol\n",
+ " Beam waist radius in meters.\n",
+ " lam : float, integer, or symbol\n",
+ " Wavelength of light in meters.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " zr : float, int, symbols\n",
+ " Rayleigh range for given beam waist and wavelength.\n",
+ " \n",
+ " lens(f)\n",
+ " Parameters\n",
+ " ----------\n",
+ " f : float or integer or sympy symbol in meters\n",
+ " Thin lens focal length in meters\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " mat : 2x2 matrix\n",
+ " [\n",
+ " [ 1, 0],\n",
+ " [-1/f, 1]\n",
+ " ]\n",
+ " \n",
+ " mult(mat1, *argv)\n",
+ " Parameters\n",
+ " ----------\n",
+ " mat1 : 2x2 ABCD matrix\n",
+ " Last matrix light interacts with.\n",
+ " *argv : 2x2 ABCD matrices \n",
+ " From left to right, the matrices should be entered such that the leftmost matrix interacts\n",
+ " with light temporally after the rightmost matrix.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " Mat : 2x2 matrix\n",
+ " The ABCd matrix describing the whole optical system.\n",
+ " \n",
+ " plot(func, var, rang=array([0. , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,\n",
+ " 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,\n",
+ " 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,\n",
+ " 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,\n",
+ " 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,\n",
+ " 0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,\n",
+ " 0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,\n",
+ " 0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,\n",
+ " 0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,\n",
+ " 0.99, 1. , 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09,\n",
+ " 1.1 , 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2 ,\n",
+ " 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3 , 1.31,\n",
+ " 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4 , 1.41, 1.42,\n",
+ " 1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5 , 1.51, 1.52, 1.53,\n",
+ " 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6 , 1.61, 1.62, 1.63, 1.64,\n",
+ " 1.65, 1.66, 1.67, 1.68, 1.69, 1.7 , 1.71, 1.72, 1.73, 1.74, 1.75,\n",
+ " 1.76, 1.77, 1.78, 1.79, 1.8 , 1.81, 1.82, 1.83, 1.84, 1.85, 1.86,\n",
+ " 1.87, 1.88, 1.89, 1.9 , 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97,\n",
+ " 1.98, 1.99, 2. , 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08,\n",
+ " 2.09, 2.1 , 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19,\n",
+ " 2.2 , 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3 ,\n",
+ " 2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4 , 2.41,\n",
+ " 2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5 , 2.51, 2.52,\n",
+ " 2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6 , 2.61, 2.62, 2.63,\n",
+ " 2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7 , 2.71, 2.72, 2.73, 2.74,\n",
+ " 2.75, 2.76, 2.77, 2.78, 2.79, 2.8 , 2.81, 2.82, 2.83, 2.84, 2.85,\n",
+ " 2.86, 2.87, 2.88, 2.89, 2.9 , 2.91, 2.92, 2.93, 2.94, 2.95, 2.96,\n",
+ " 2.97, 2.98, 2.99]))\n",
+ " Parameters\n",
+ " ----------\n",
+ " func : Sympy function of one variable\n",
+ " Sympy function defining the beam width after the last optical element.\n",
+ " var : sympy variable\n",
+ " Variable in func that will be plotted.\n",
+ " rang : numpy array\n",
+ " Array of the values along the optical axis to be plotted\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " plot : matplotlib graph\n",
+ " Graph of the beam width of var\n",
+ " \n",
+ " prop(d)\n",
+ " Parameters\n",
+ " ----------\n",
+ " d : float or integer or sympy symbol\n",
+ " Distance light is propagating along the z-axis.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " mat: 2x2 matrix\n",
+ " [\n",
+ " [1, d],\n",
+ " [0, 1]\n",
+ " ]\n",
+ " \n",
+ " q1_func(z, w0, lam, mat)\n",
+ " Parameters\n",
+ " ----------\n",
+ " z : float, int, symbol\n",
+ " Position of the beam waist in meters.\n",
+ " w0 : float, int, symbol\n",
+ " Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).\n",
+ " lam : float, int, symbol\n",
+ " Wavelength of light in meters.\n",
+ " mat : float, int, symbol\n",
+ " The ABCD 2x2 matrix describing the optical system.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " z: float, int, symbol\n",
+ " Position of the beam waist after the optical system\n",
+ " zr: float, int, symbol\n",
+ " Rayleigh range of the beam after the optical system\n",
+ " \n",
+ " q1_inv_func(z, w0, lam, mat)\n",
+ " Parameters\n",
+ " ----------\n",
+ " z : float, int, symbol\n",
+ " Position of the beam waist in meters.\n",
+ " w0 : float, int, symbol\n",
+ " Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).\n",
+ " lam : float, int, symbol\n",
+ " Wavelength of light in meters.\n",
+ " mat : float, int, symbol\n",
+ " The ABCD 2x2 matrix describing the optical system.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " R : float, int, symbol\n",
+ " Radius of curvature of the wavefront in meters.\n",
+ " w : float, int, symbol\n",
+ " Radius of the beam in meters.\n",
+ " \n",
+ " ray(y, theta)\n",
+ " Parameters\n",
+ " ----------\n",
+ " y : float or integer or sympy symbol in meters\n",
+ " The vertical height of a ray.\n",
+ " theta : float or integer in radians\n",
+ " The angle of divergence of the ray.\n",
+ " \n",
+ " Returns\n",
+ " -------\n",
+ " mat : 2x1 matrix\n",
+ " [\n",
+ " [y],\n",
+ " [teta]\n",
+ " ]\n",
+ "\n",
+ "DATA\n",
+ " oo = oo\n",
+ "\n",
+ "FILE\n",
+ " c:\\users\\wrighta\\documents\\beamprop\\beamprop_script.py\n",
+ "\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "help(bs)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Let's first see how we define a beam and how we can visualize it propagating."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### A Gaussian beam can be defined by it's (radial) waist, $w_0$, it's Rayleigh range, $z_R = \\frac{\\pi * w_0^2}{\\lambda}$, and the location of its waist, $z_0$."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "w0 = 1E-3 # 1mm beam waist\n",
+ "lam = 355E-9 # wavelength of 355 nm (UV)\n",
+ "zR = bs.Zr(w0, lam) # Rayleigh range in m\n",
+ "z0 = 0 # location of waist in m"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### We now want to define our \"optical system\" using matrices. For this first example, we will just use a free space propagation matrix, and let the beam propagate a distance $d$ which we will define using a symbol."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "d = sym.symbols('d')\n",
+ "M = bs.prop(d)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### We now use the bs script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "R, w = bs.q1_inv_func(0, w0, lam, M)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "w = 0.001*(0.0127690021685256*d**2 + 1)**0.5\n"
+ ]
+ }
+ ],
+ "source": [
+ "print('w = {}'.format(w))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### And as simple as that, we have a function for our waist. Let's plot it and see what it looks like"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEGCAYAAAC+fkgiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5xdZX3v8c83M+QKJOQ2hknIhDBJCbSKTLkoB0ZBiNWXiBUNPSo9haYotp7WngrHHvCFpcW23nq8NSqCl4oUraSKQET3QS23gCgkISSBkAwJJCGBZMiNzPzOH2tt9mX2ntlz2XsnO9/367Vee+1nP2utZz2EfLPWetZaigjMzMxqYVS9G2BmZocPh46ZmdWMQ8fMzGrGoWNmZjXj0DEzs5pprncDDmZTp06Ntra2IS//8ssvM2HChJFr0CHMfVHI/ZHjvijUCP3x8MMPb4uIaaV+c+j0o62tjeXLlw95+UwmQ2dn58g16BDmvijk/shxXxRqhP6Q9Ey53+p6ek3SQkmrJa2VdFWJ38dI+l76+wOS2vJ+uzotXy3pgrzyGyVtkfR40bo+IelZSY+m0x9Uc9/MzKyvuoWOpCbgi8BbgQXAJZIWFFW7DNgREScAnwU+lS67AFgEnAQsBL6Urg/gprSslM9GxOvS6Y6R3B8zMxtYPY90TgPWRsRTEbEfuAW4sKjOhcDN6fxtwLmSlJbfEhH7IuJpYG26PiLiXmB7LXbAzMwGp57XdFqBjXnfu4DTy9WJiAOSXgKmpOX3Fy3bWsE2PyzpA8By4KMRsaO4gqTFwGKAlpYWMplMRTtTSnd397CWbyTui0Lujxz3RaFG7496ho5KlBU/CK5cnUqWLfZl4JNpvU8Cnwb+pM9KIpYASwA6OjpiOBf0GuGC4EhxXxRyf+S4Lwo1en/U8/RaFzAr7/tMYFO5OpKagYkkp84qWbZARDwfET0R0Qt8lfR0nJmZ1U49Q+choF3SHEmjSQYGLC2qsxS4NJ1/N/CzSB6LvRRYlI5umwO0Aw/2tzFJM/K+XgQ8Xq6umZlVR91Or6XXaD4M3AU0ATdGxApJ1wHLI2Ip8HXgW5LWkhzhLEqXXSHpVmAlcAC4MiJ6ACR9F+gEpkrqAq6NiK8D/yjpdSSn19YDf1a7vTUzO7j19MD69bBqVTLNng3vec/Ib6euN4emw5bvKCq7Jm9+L3BxmWWvB64vUX5JmfrvH1ZjzcwawL59sGZNLlxWrYKVK+HJJ2Hv3ly9RYsaMHTMzKw6du6EJ54oDJdVq+Cpp5Kjmqy2NjjxRDjvvOQzO02eXJ12OXTMzA5REbB1a3KkUhwuzz6bq3fEEdDeDr/3e/De9+aCZf58GD++tm126JiZHeR6e2HDhr7BsmoVbM+7FX7ChCRM3vzmwqOW449Pgudg4NAxMztI7N8P69eP5/vfLwyW1ath9+5cvalTYcECuPjiwnCZORNU6i7Gg4hDx8yshiKSU19PPplMq1fn5p9+Gnp6crcQHndcEibnnFMYLlOn1nEHhsmhY2ZWBS++WDpYnnyy8Khl/HiYNw9e//pkxFhv7yre9a4TmT8fjjyyfu2vFoeOmdkQ7dsH69aVDpctW3L1mppgzpwkXDo7kwv48+YlU2tr4SmxTOZ5Tj31xJrvS604dMzM+tHbC11dpYNl/frk96yWliRQ3vGOJFCy4XL88TB6dN124aDi0DEzIxkFln8KLBsua9bAnj25ekcemQTJaafB+96XC5b2dpg4sX7tP1Q4dMzssNDbC5s3J6fD1q5NPvOnHXkvOmlqgrlzkzA577zCo5YZMw7+EWIHM4eOmTWMZMhx30BZty65Ez//MS9NTcnzxebOTW6YPOGEXLjMmXPw3NfSaBw6ZnZI2bmzdKisWwcbNxZeYxk/PnfE8ta3JvPZ6bjjHCz14NAxs4NKBDz/fC5Iik+FbdtWWH/q1CREzjqrMFTmzk0u7PtU2MHFoWNmNffyy/DMM8l0993H8qMfFZ4Ge/nlXN1Ro2DWrCRELrqob7AcfXT99sMGz6FjZiPuxReTayvZYMmff+aZ4qOVeYwdmwwrnjsXzj23MFTa2jzcuJE4dMxsULJPNi4OkvzvO3cWLjNuXHLRfvZs6OjIzc+eDZs23ccf/uGZjKrne4ytZhw6ZlagpycZWlwcJNn5DRsK71uB5BRXW1synXNO8pkNlba25LpLuWsrmcw+B85hxKFjdpjp7k4eOPnss4VHKtlg2bgRDhwoXGbq1CQ8Tj4Z3va2wlCZPRsmTarDjtghyaFj1iB6epJRX9lAyZ82bcrNF5/6kpIbHmfPhjPOSO5ZyR6hzJ6dDC2eMKEuu2QNyKFjdgjYubN8iGSn554rvEcFkhsgZ8xIHiqZfSVxaysce2zyOXt28g6WMWPqs192+HHomNXRK68kYVEuSLJl3d19l500KQmO1lY46aTcfP40bVoSPGYHC4eO2Qh75ZVkdNeWLbnP7JT9vnr1qezcmZwOiyhcvrk5dyTyu78LCxfmQiRbfuyxPuVlhyaHjtkAenqSJxCXC5DisvwHR+Zrbk6OPKZNg0mT9nPWWaWPTqZOxaO5rGE5dOywE5FcI6kkQLZsSW5kLL5WAskF+ClTYPr0ZHrta3Pz06cn4ZL/fdKk3LDhTOYxOjs7a7rfZgcDh44dknp7Ydeu5M73HTuSz0qmbduSMNm/v/R6J07MhcQJJ8Ab3lA+RKZM8fUSs8Fy6FhdRORCo9zUX5i89FLfayHFjj46ObrITm1tcMopyUMgS4XI1KkexWVWbXUNHUkLgc8DTcDXIuKGot/HAN8ETgVeAN4bEevT364GLgN6gL+IiLvS8huBtwNbIuLkvHVNBr4HtAHrgfdERJmz7wbJX+p798Lu3cm0Z09ufqCpuO7Gja9FygXJSy+VPmWV78gj4ZhjcqExa1ZyYT0/SLJTfr1Jk5LA8VGI2cGnbqEjqQn4IvAWoAt4SNLSiFiZV+0yYEdEnCBpEfAp4L2SFgCLgJOAY4GfSpoXET3ATcAXSMIq31XAPRFxg6Sr0u8fq94e5vT2Jnd49/QkU/588feRnn/llfJhMVCI7Nkz8NFEKWPHJu8xGT8+eebW+PFw4MAo2tqSe0XKBUX+NHFicuHdzBpLPf+3Pg1YGxFPAUi6BbgQyA+dC4FPpPO3AV+QpLT8lojYBzwtaW26vvsi4l5JbSW2dyHQmc7fDGSoUugsX5682+OVV84Z8F/ztTR6dC4MiqdJkwpDYqCpXL1x40qPvMpkfu0L52ZW19BpBTbmfe8CTi9XJyIOSHoJmJKW31+0bOsA22uJiM3pujZLml6qkqTFwGKAlpYWMplMRTuTb+vWMVx0USs9PfsYO/YImpqCpqZg1CgYNSo7H6/ONzVVt3zcuF5Gj+4Z0dNNPT3JNZlduyqr393dPaS+bFTujxz3RaFG7496hk6pZ84Wn8wpV6eSZYckIpYASwA6OjpiqP86v/hiyGQy/td9yn1RyP2R474o1Oj9Uc9b0LqAWXnfZwKbytWR1AxMBLZXuGyx5yXNSNc1A9gy5JabmdmQ1DN0HgLaJc2RNJpkYMDSojpLgUvT+XcDP4uISMsXSRojaQ7QDjw4wPby13UpcPsI7IOZmQ1C3UInIg4AHwbuAlYBt0bECknXSXpHWu3rwJR0oMBfkYw4IyJWALeSDDq4E7gyHbmGpO8C9wHzJXVJuixd1w3AWyStIRkxVzA828zMqq+ug1Ij4g7gjqKya/Lm9wIXl1n2euD6EuWXlKn/AnDucNprZmbD48cKmplZzTh0zMysZhw6ZmZWMw4dMzOrGYeOmZnVjEPHzMxqxqFjZmY149AxM7OaceiYmVnNOHTMzKxmHDpmZlYzDh0zM6sZh46ZmdWMQ8fMzGrGoWNmZjXj0DEzs5px6JiZWc04dMzMrGYcOmZmVjMOHTMzqxmHjpmZ1YxDx8zMasahY2ZmNePQMTOzmnHomJlZzdQ1dCQtlLRa0lpJV5X4fYyk76W/PyCpLe+3q9Py1ZIuGGidkm6S9LSkR9PpddXePzMzK9Rcrw1LagK+CLwF6AIekrQ0IlbmVbsM2BERJ0haBHwKeK+kBcAi4CTgWOCnkualy/S3zv8VEbdVfefMzKykeh7pnAasjYinImI/cAtwYVGdC4Gb0/nbgHMlKS2/JSL2RcTTwNp0fZWs08zM6qRuRzpAK7Ax73sXcHq5OhFxQNJLwJS0/P6iZVvT+f7Web2ka4B7gKsiYl9xoyQtBhYDtLS0kMlkBrdXebq7u4e1fCNxXxRyf+S4Lwo1en/UM3RUoiwqrFOuvNSRW3adVwPPAaOBJcDHgOv6VI5Ykv5OR0dHdHZ2llhlZTKZDMNZvpG4Lwq5P3LcF4UavT8GPL0m6UxJX5T0W0lbJW2QdIekKyVNHMa2u4BZed9nApvK1ZHUDEwEtvezbNl1RsTmSOwDvkFyKs7MzGqo39CR9BPgcuAuYCEwA1gA/C0wFrhd0juGuO2HgHZJcySNJhkYsLSozlLg0nT+3cDPIiLS8kXp6LY5QDvwYH/rlDQj/RTwTuDxIbbbzMyGaKDTa++PiG1FZd3AI+n0aUlTh7Lh9BrNh0kCrQm4MSJWSLoOWB4RS4GvA9+StJbkCGdRuuwKSbcCK4EDwJUR0QNQap3pJr8jaRrJqblHgSuG0m4zMxu6fkOnOHAkHZ2/TERsLxFKFYuIO4A7isquyZvfC1xcZtnrgesrWWda/uahttPMzEZGRQMJJP0ZyUX3PeQuzAdwfJXaZWZmDajS0Wt/DZw0nKMaMzOzSm8OXQfsrmZDzMys8VV6pHM18F+SHgBevaEyIv6iKq0yM7OGVGno/CvwM+AxoLd6zTEzs0ZWaegciIi/qmpLzMys4VV6TefnkhZLmiFpcnaqasvMzKzhVHqk80fp59V5ZR4ybWZmg1JR6ETEnGo3xMzMGt9Az147a4Dfj5Z08sg2yczMGtVARzp/KOkfgTuBh4GtJA/6PAF4EzAb+GhVW2hmZg1joGev/aWkY0ie8HwxyVOm9wCrgH+NiF9Wv4lmZtYoBrymExE7gK+mk5mZ2ZBVOmTazMxs2Bw6ZmZWMw4dMzOrmYpCR9J4Sf9H0lfT7+2S3l7dppmZWaOp9EjnGyRPlz4z/d4F/F1VWmRmZg2r0tCZGxH/CLwCEBF7AFWtVWZm1pAqDZ39ksaRvqpa0lzy3qtjZmZWiUof+PkJkqcSzJL0HeCNwB9XqU1mZtagKn3g592SHgbOIDmt9pGI2FbVlpmZWcOpdPTaPcDpEfHjiPhRRGyTtKTKbTMzswZT6TWdOcDHJF2bV9ZRhfaYmVkDqzR0XgTOBVok/aekiVVsk5mZNahKQ0cRcSAiPgR8H/glML16zTIzs0ZUaeh8JTsTETeRjFy7e7gbl7RQ0mpJayVdVeL3MZK+l/7+gKS2vN+uTstXS7pgoHVKmpOuY026ztHDbb+ZmQ3OQG8OPTqd/XdJk7MT8DTw18PZsKQm4IvAW4EFwCWSFhRVuwzYEREnAJ8FPpUuuwBYBJwELAS+JKlpgHV+CvhsRLQDO9J1m5lZDQ10pPNv6efDwPL08+G878NxGrA2Ip6KiP3ALcCFRXUuBG5O528DzpWktPyWiNgXEU8Da9P1lVxnusyb03WQrvOdw2y/mZkN0kBvDn17+jmnCttuBTbmfe8CTi9XJyIOSHoJmJKW31+0bGs6X2qdU4AXI+JAifoFJC0GFgO0tLSQyWQGtVP5uru7h7V8I3FfFHJ/5LgvCjV6f1R0c6ikNwKPRsTLkt4HvB74XERsGMa2Sz27LSqsU6681JFbf/X7FkYsAZYAdHR0RGdnZ6lqFclkMgxn+Ubivijk/shxXxRq9P6odCDBl4Hdkl4L/A3wDPCtYW67C5iV930msKlcHUnNwERgez/LlivfBkxK11FuW2ZmVmWVhs6BiAiSaymfj4jPA0cNc9sPAe3pqLLRJAMDlhbVWQpcms6/G/hZ2o6lwKJ0dNscoB14sNw602V+nq6DdJ23D7P9ZmY2SJU+8HOXpKuB9wFnp6PEjhjOhtNrNB8G7gKagBsjYoWk64DlEbEU+DrwLUlrSY5wFqXLrpB0K7ASOABcGRE9AKXWmW7yY8Atkv4O+HW6bjMzq6FKQ+e9wB8Bl0XEc5KOA/5puBuPiDuAO4rKrsmb3wtcXGbZ64HrK1lnWv4Uyeg2MzOrk0qfMv0c8Jm87xuAb1arUWZm1pgqvaZjZmY2bA4dMzOrGYeOmZnVTKUvcXu7pF9L2i5pp6RdknZWu3FmZtZYKh299jngXcBj6T0vZmZmg1bp6bWNwOMOHDMzG45Kj3T+BrhD0v8D9mULI+Iz5RcxMzMrVGnoXA90A2MBv/zMzMyGpNLQmRwR51e1JWZm1vAqvabzU0kOHTMzG5ZKQ+dK4E5Jezxk2szMhqrSZ68N9zUGZmZmFV/TQdIxJO+tGZsti4h7q9EoMzNrTJW+rvpy4CMkb9x8FDgDuA94c/WaZmZmjabSazofAX4feCYi3gScAmytWqvMzKwhVRo6e9MXqiFpTEQ8AcyvXrPMzKwRVXpNp0vSJOCHwDJJO4BN1WuWmZk1okpHr12Uzn5C0s+BicCdVWuVmZk1pMGMXjsLaI+Ib0iaBrQCT1etZWZm1nAqfZ/OtcDHgKvToiOAb1erUWZm1pgqHUhwEfAO4GWAiNgE+IZRMzMblEpDZ3/6Lp0AkDShek0yM7NGVWno3CrpX4FJkv4U+Cnw1eo1y8zMGlGlo9f+WdJbgJ0k9+dcExHLqtoyMzNrOBWPXktDZpmkqcAL1WuSmZk1qn5Pr0k6Q1JG0g8knSLpceBx4HlJC4e6UUmTJS2TtCb9PKZMvUvTOmskXZpXfqqkxyStlfQvktTfeiV1SnpJ0qPpdM1Q225mZkM30DWdLwB/D3wX+BlweUS8Bjgb+IdhbPcq4J6IaAfuSb8XkDQZuBY4HTgNuDYvnL4MLCZ56nU7kA3A/tb7i4h4XTpdN4y2m5nZEA0UOs0RcXdE/DvwXETcD5A+e204LgRuTudvBt5Zos4FwLKI2B4RO4BlwEJJM4CjI+K+dETdN/OWr2S9ZmZWJwNd0+nNm99T9FsMY7stEbEZICI2S5peok4rsDHve1da1prOF5cPtN4zJf2G5Jlxfx0RK0o1TNJikqMoWlpayGQyg923V3V3dw9r+Ubivijk/shxXxRq9P4YKHRem76WWsC4vFdUi7yXuZUi6afAa0r89PEK26YSZdFPeX8eAWZHRLekPyB5cGl7qYoRsQRYAtDR0RGdnZ0VNrevTCbDcJZvJO6LQu6PHPdFoUbvj35DJyKahrriiDiv3G+Snpc0Iz0amQFsKVGtC+jM+z4TyKTlM4vKs0+8LrneiMiGJRFxh6QvSZoaEduGsGtmZjZEld4cOtKWAtnRaJcCt5eocxdwvqRj0gEE5wN3pafPdqUj6wR8IG/5kuuV9Jq8EW6nkey3h32bmdVYxffpjLAbSJ5ycBmwAbgYQFIHcEVEXB4R2yV9EngoXea6iNiezn8QuAkYB/wkncquF3g38EFJB0iuTS1KByGYmVkN1SV0IuIF4NwS5cuBy/O+3wjcWKbeyYNY7xdIhn+bmVkd1ev0mpmZHYYcOmZmVjMOHTMzqxmHjpmZ1YxDx8zMasahY2ZmNePQMTOzmnHomJlZzTh0zMysZhw6ZmZWMw4dMzOrGYeOmZnVjEPHzMxqxqFjZmY149AxM7OaceiYmVnNOHTMzKxmHDpmZlYzDh0zM6sZh46ZmdWMQ8fMzGrGoWNmZjXj0DEzs5px6JiZWc04dMzMrGbqEjqSJktaJmlN+nlMmXqXpnXWSLo0r/xUSY9JWivpXyQpLb9Y0gpJvZI6itZ1dVp/taQLqruHZmZWSr2OdK4C7omIduCe9HsBSZOBa4HTgdOAa/PC6cvAYqA9nRam5Y8D7wLuLVrXAmARcFJa90uSmkZ4n8zMbADNddruhUBnOn8zkAE+VlTnAmBZRGwHkLQMWCgpAxwdEfel5d8E3gn8JCJWpWWltndLROwDnpa0liTI7hvRvUpt2wb/+Z+wZk0LmzZBUxM0Nyef1ZwfNQr67rqZ2cGjXqHTEhGbASJis6TpJeq0AhvzvnelZa3pfHF5f1qB+ytZRtJikqMoWlpayGQyA6y6r1WrjuJDHzoVOHHQyw7XqFFBU1MwalQyNTcHY8b0MmZMD2PG9DJ2bP+f2bpjx/b/mb/OURUcL3d3dw+pLxuV+yPHfVGo0fujaqEj6afAa0r89PFKV1GiLPopH8q6+hZGLAGWAHR0dERnZ+cAq+7rzDPhrW+FX/3qATo6TqenBw4cgJ4e+sz399vQ5kVPj14tf+UV2LMHdu/uO+3Y0bdsKMaOhfHjS0/jxiWfu3Zt5sQTZzBpEn2mY47JzY8ff3gcrWUyGYbyZ6sRuS8KNXp/VC10IuK8cr9Jel7SjPQoZwawpUS1LnKn4ABmkpyG60rn88s3DdCcLmDWIJcZsjFjoK0N1q/fw/z51drKyIuAvXtLB9Tu3eXDq7+627Ylny+8MJlMZuBga27uG0r9hVTxNG7c4RFaZoeqep1eWwpcCtyQft5eos5dwN/nDR44H7g6IrZL2iXpDOAB4APA/61ge/8m6TPAsSSDDx4c/m40Fin5S3vcOJgyZWTXncncR2dnJ/v3w4svDm569tnc/J49/W/niCP6D6vJk2H69GSaNi03P378yO6vmZVWr9C5AbhV0mXABuBigHSY8xURcXkaLp8EHkqXuS47qAD4IHATMA74SToh6SKSAJoG/FjSoxFxQUSskHQrsBI4AFwZET212FErNHp07i/6odi3b/ChtXFj8rljR7J8KRMm9A2iUuE0fTpMnZrsh5kNXl1CJyJeAM4tUb4cuDzv+43AjWXqnVyi/D+A/yizzeuB64feajsYjBkDLS3JNBQvvwxbtsDWrcln/pQt6+qCRx5J5g8cKL2eSZP6D6b875MnJyMMzax+RzpmdTFhAsyZk0wDiYCXXiofTtnpiSfgF79Irl9FieEpo0YlR0f5QdTTM5eHHoLW1sJp3LiR32ezg4lDx6wMKXc9aN68gev39MALL5QPp2zZww9DV9ex3HZb33Uccwwce2zfMGptzZVPn05Fw9TNDkYOHbMR0tRU+fWqn//8F7z+9Z08+yyvTps2UfD9scfg+eeht7dw2eZmmDGjdCDlTxMmVGc/zYbDoWNWBxJMnJhMCxaUr3fgQBI8+WGUH1ArVsDdd8OuXX2XnThx4KOm17zGR01WWw4ds4NYc3MuKPqza1fpo6Xs9MQTsHlzcgow3+jRMGsWzJ6dTG1thfOtrckwdLOR4tAxawBHHQW/8zvJVE5PT3JdKRtEXV2wYQM880wy3XlnEkz5Ro1Kgqc4kLLfjzsueSKFWaUcOmaHiaam5FrQjBnQ0VG6zt69yX1NzzwD69fnAumZZ+Dee5OwKj5aamkpHUjZ+aOOqvKO2SHFoWNmrxo7Ftrbk6mUAweS4CkVSo88Aj/8IezfX7jM5MnlA2n27NLDzK1xOXTMrGLNzbmwOPvsvr/39iYDH4oDaf16ePJJWLYsuUE337hxZ9HeDiecAHPnFk6zZiXbtMbh/5xmNmJGjcqdwjvzzL6/R8D27YWh9MtfPse+fTNZuRJ+9KPCI6Xm5uTIqDiM5s6F44/3M/MORQ4dM6sZKXmY7JQpcOqpSdkpp6ylszN5cHxvb3L6bt26wmntWrj//uQJEflmzCgMovyjpcmT/cTxg5FDx8wOGqNGJafUZs2C4lfKZI+SigNp3brktN3NNxfWnzix9BHS3Lkwc6bvT6oXh46ZHRLyj5JOO63v77t3w9NP9w2kRx9NBji88kqu7ujRyfP3io+S5s1Lyn0dqXrctWbWEMaPh5NOSqZiPT3JUPDsqbr8ULr3XujuztVtbk5CaN48mD8/+czOt7T4lN1wOXTMrOE1NSUDEtra4Nyil6pEJA9iXbs2GWGXnVavTh4xlP8OpqOOKh1G7e2+H6lSDh0zO6xJuQe1vuENhb/19iZHSKtXF4bRf/0XfPe7hfcYzZhROpDmzPGjhPI5dMzMyhg1Kndf0vnnF/62d29yeq44kH7wg+TdSllNTcnw7uIwmjcvCarD7XSdQ8fMbAjGji1/DWn79r6n6p58Eu65B/bsydU78shcEGWn7u6j6OhIfmtEDh0zsxE2eTKccUYy5cveh1R8dPTgg3Drrdl3J53KFVckw8ZPPLHvNG1aPfZo5Dh0zMxqJP8+pPPOK/xt377kdN1ttz1OU9PJrFoFq1bBV7+aDAfPmjIleQdTcRjNmnVonKpz6JiZHQTGjEnC5OyztxXcGJsdzJANoez0/e8nr0fPmjAhebVFcRjNnXtwDWRw6JiZHcTyBzMsXFj429atsHJlYRhlMvDtb+fqNDcnQ7qLw2j+/Pq80tyhY2Z2iJo2Dc45J5ny7dqVvC02P4wefxxuv73wfUizZ5e+bjRlSvXa7NAxM2swRx0Fv//7yZRv377kJtj8MFq5Mjk62rs3V2/aNHj/++HTnx75tjl0zMwOE2PGlB7m3dOTvGYiP4xmzapOGxw6ZmaHuewNrMcfD297W3W3VZeHe0uaLGmZpDXp5zFl6l2a1lkj6dK88lMlPSZpraR/kZKBgpIulrRCUq+kjrz6bZL2SHo0nb5S/b00M7Ni9XqjxFXAPRHRDtyTfi8gaTJwLXA6cBpwbV44fRlYDLSnU3ZMx+PAu4B7S2xzXUS8Lp2uGMmdMTOzytQrdC4Esq9cuhl4Z4k6FwDLImJ7ROwAlgELJc0Ajo6I+yIigG9ml4+IVRGxuvrNNzOzoajXNZ2WiNgMEBGbJU0vUacV2Jj3vSsta03ni8sHMkfSr4GdwN9GxC9KVZK0mOQoipaWFjKZTAWrLq27u3tYyzcS90Uh90eO+6JQo/dH1UJH0k+B15T46eOVrqJEWfRT3p/NwHER8YKkU4EfSjopInb2WVHEEmAJQEdHR3QWvzN3EDKZDMNZvpG4Lwq5P3LcF4UavT+qFjoRcV653yQ9L2lGepQzA9hSoloX0Jn3fSaQSctnFpVvGqAt+4B96fzDktYB88iWQi0AAAblSURBVIDlA++JmZmNlHpd01kKZEejXQrcXqLOXcD5ko5JBxCcD9yVnpbbJemMdNTaB8os/ypJ0yQ1pfPHkww+eGpkdsXMzCpVr9C5AXiLpDXAW9LvSOqQ9DWAiNgOfBJ4KJ2uS8sAPgh8DVgLrAN+ki5/kaQu4Ezgx5LuSuufDfxW0m+A24Ar8tZlZmY1ooiBLoccviRtBZ4ZxiqmAtsGrHV4cF8Ucn/kuC8KNUJ/zI6Ikm/+cehUkaTlEdExcM3G574o5P7IcV8UavT+qNfpNTMzOww5dMzMrGYcOtW1pN4NOIi4Lwq5P3LcF4Uauj98TcfMzGrGRzpmZlYzDh0zM6sZh04VSFooaXX6vp8+r204nEiaJennklal7zr6SL3bVG+SmiT9WtKP6t2WepM0SdJtkp5I/4ycWe821ZOkv0z/P3lc0nclja13m0aaQ2eEpY/b+SLwVmABcImkBfVtVV0dAD4aEScCZwBXHub9AfARYFW9G3GQ+DxwZ0T8DvBaDuN+kdQK/AXQEREnA03Aovq2auQ5dEbeacDaiHgqIvYDt5C8P+iwFBGbI+KRdH4XyV8qlbyKoiFJmgm8jeQxToc1SUeTPKLq6wARsT8iXqxvq+quGRgnqRkYzwAPMz4UOXRGXrn3AB32JLUBpwAP1LcldfU54G+A3no35CBwPLAV+EZ6uvFrkibUu1H1EhHPAv8MbCB5HctLEXF3fVs18hw6I28o7/tpeJKOBL4P/M9S7zE6HEh6O7AlIh6ud1sOEs3A64EvR8QpwMuUeHX94SJ9mv6FwBzgWGCCpPfVt1Ujz6Ez8rqAWXnfB3zfT6OTdARJ4HwnIn5Q7/bU0RuBd0haT3La9c2Svl3fJtVVF9AVEdkj39tIQuhwdR7wdERsjYhXgB8Ab6hzm0acQ2fkPQS0S5ojaTTJhcCldW5T3aTvPPo6sCoiPlPv9tRTRFwdETMjoo3kz8XPIqLh/iVbqYh4DtgoaX5adC6wso5NqrcNwBmSxqf/35xLAw6sqNqbQw9XEXFA0odJXkLXBNwYESvq3Kx6eiPwfuAxSY+mZf87Iu6oY5vs4PHnwHfSf6A9BfyPOrenbiLiAUm3AY+QjPr8NQ34SBw/BsfMzGrGp9fMzKxmHDpmZlYzDh0zM6sZh46ZmdWMQ8fMzGrGoWM2CJJmSrpd0hpJ6yR9Ph3u298ykyR9KO/7senQ2MFst1nSNkn/UGH9OyRNGsT6Pyfp7EHUnybpzkrrm2U5dMwqlN6w9wPghxHRDswDjgSuH2DRScCroRMRmyLi3YPc/PnAauA9aTv6FRF/UOnDMyVNBs6IiHsrbUxEbAU2S3pjpcuYgUPHbDDeDOyNiG8AREQP8JfAn6R3kf9xehR0Z/o+pWvT5W4A5kp6VNI/SWqT9Di8+m6df5b0mKTfSvrzMtu+hOQ1ABtIXhGBpInpduan378r6U/T+fWSpkqaIOnHkn6TvqPlvSXW/W7g1aOWdNm/l3SfpOWSXi/prvTI7oq85X4I/Pch9aQdtvxEArPKnQQUPKwzInZK2gCckBadBpwM7AYekvRjkodYnhwRr4NXn7adtZjkAY+npE+zmFy8UUnjSB6J8mckR02XAPdFxEvp0y9ukvR54JiI+GrR4guBTRHxtnRdE0vs1xtJnnuWb2NEnCnps8BNaZ2xwArgK2md5cDflVifWVk+0jGrnCj9xPD88mUR8UJE7CE5FXfWAOs8D/hKRBwAiIjtJeq8Hfh5ROwmeXDqRenLAomIZcBjJC8OvLzEso8B50n6lKT/FhEvlagzg+QVA/myzwt8DHggInalp9T25l0r2kLyNGSzijl0zCq3AujIL0hfRDYLWJcWFYfSQM+ZKhdk+S4hCY71JEdaU4A3pdsfBZwI7AH6HCVFxJPAqSTh8Q+Srimx/j0kRzH59qWfvXnz2e/ZMyRj02XNKubQMavcPcB4SR+AV19N/mngpvQoBOAtkianp8TeCfwK2AUcVWaddwNXpG+KpPj0WhpqZwHHRURb+oTqK0mCCJJrSqvS7zemr5HIX/5YYHdEfJvkBWGlXh2witzpwcGYBzw+hOXsMObQMatQJE/HvQi4WNIa4ElgL/C/86r9EvgW8Cjw/YhYHhEvAL9KL+T/U9Fqv0YyOOC3kn4D/FHR7+8ieQVC/tHG7STv5ZlPckrtoxHxC+Be4G+Llv9d4MH0Cd8fp/Q1mB8DnQN2QF9vSpc1q5ifMm02QiT9MdARER+ud1sGS9IvgbdXOsw6XeZe4MKI2FG9llmj8ZGOmQF8FDiu0sqSpgGfceDYYPlIx8zMasZHOmZmVjMOHTMzqxmHjpmZ1YxDx8zMasahY2ZmNfP/AbRYoKP3tc6bAAAAAElFTkSuQmCC",
+ "text/plain": [
+ "