diff options
author | Hunter Rew <hbrew@email.wm.edu> | 2014-04-23 23:40:07 -0400 |
---|---|---|
committer | Hunter Rew <hbrew@email.wm.edu> | 2014-04-23 23:40:07 -0400 |
commit | 2cabf85e1496bb60a3a1ab50784aa0ba2e743660 (patch) | |
tree | 0d2db0c5af59d2b49d878daf99e8e0afbe552222 /source | |
parent | 41ab7ba66c047aa27e465484e558cd11481a8e40 (diff) | |
download | eit_filter_simulations-2cabf85e1496bb60a3a1ab50784aa0ba2e743660.tar.gz eit_filter_simulations-2cabf85e1496bb60a3a1ab50784aa0ba2e743660.zip |
removed compiled files
Diffstat (limited to 'source')
-rwxr-xr-x | source/eit1 | bin | 73741 -> 0 bytes | |||
-rw-r--r-- | source/eit1.cc | 2129 | ||||
-rw-r--r-- | source/eit1.h5 | bin | 33656 -> 0 bytes | |||
-rw-r--r-- | source/eit1.py | 17 | ||||
-rw-r--r-- | source/eit1.pyc | bin | 1217 -> 0 bytes | |||
-rw-r--r-- | source/eit1.xmds | 136 | ||||
-rw-r--r-- | source/eit1.xsil | 168 | ||||
-rwxr-xr-x | source/eit2 | bin | 73741 -> 0 bytes | |||
-rw-r--r-- | source/eit2.cc | 2129 | ||||
-rw-r--r-- | source/eit2.h5 | bin | 33656 -> 0 bytes | |||
-rw-r--r-- | source/eit2.py | 17 | ||||
-rw-r--r-- | source/eit2.pyc | bin | 1217 -> 0 bytes | |||
-rw-r--r-- | source/eit2.xmds | 136 | ||||
-rw-r--r-- | source/eit2.xsil | 168 | ||||
-rwxr-xr-x | source/eit3 | bin | 73741 -> 0 bytes | |||
-rw-r--r-- | source/eit3.cc | 2129 | ||||
-rw-r--r-- | source/eit3.h5 | bin | 33656 -> 0 bytes | |||
-rw-r--r-- | source/eit3.py | 17 | ||||
-rw-r--r-- | source/eit3.pyc | bin | 1217 -> 0 bytes | |||
-rw-r--r-- | source/eit3.xmds | 136 | ||||
-rw-r--r-- | source/eit3.xsil | 168 |
21 files changed, 0 insertions, 7350 deletions
diff --git a/source/eit1 b/source/eit1 Binary files differdeleted file mode 100755 index fcdb4a5..0000000 --- a/source/eit1 +++ /dev/null diff --git a/source/eit1.cc b/source/eit1.cc deleted file mode 100644 index f940eb5..0000000 --- a/source/eit1.cc +++ /dev/null @@ -1,2129 +0,0 @@ -// ******************************************************** -// simulation logging - -#define _SAMPLE_LOG_LEVEL (1 << 0) -#define _SEGMENT_LOG_LEVEL (1 << 1) -#define _PATH_LOG_LEVEL (1 << 2) -#define _SIMULATION_LOG_LEVEL (1 << 3) -#define _WARNING_LOG_LEVEL (1 << 4) -#define _ERROR_LOG_LEVEL (1 << 5) -#define _NO_ERROR_TERMINATE_LOG_LEVEL (1 << 6) -#define _ALL_LOG_LEVELS _SAMPLE_LOG_LEVEL|_SEGMENT_LOG_LEVEL|_PATH_LOG_LEVEL|_SIMULATION_LOG_LEVEL|_WARNING_LOG_LEVEL|_ERROR_LOG_LEVEL|_NO_ERROR_TERMINATE_LOG_LEVEL -#define _LOG_LEVELS_BEING_LOGGED (_ALL_LOG_LEVELS) - -#define real Re -#define imag Im - -#include <complex> - -#undef real -#undef imag - - -#include <stdio.h> - -#define _LOG(logLevel, ...) \ - do { \ - if (logLevel & _LOG_LEVELS_BEING_LOGGED) { \ - if (logLevel & (_ERROR_LOG_LEVEL | _WARNING_LOG_LEVEL)) \ - printf("%s:%i: ", __FILE__, __LINE__); \ - printf(__VA_ARGS__); \ - fflush(stdout); \ - if (logLevel & (_ERROR_LOG_LEVEL | _NO_ERROR_TERMINATE_LOG_LEVEL)) \ - exit(logLevel == _ERROR_LOG_LEVEL); \ - } \ - } while (0) - -// ******************************************************** -// simulation includes - -#include <xpdeint_platform.h> -#include <cmath> -#include <string> -#include <cstring> -#include <fstream> -#include <sstream> -#include <cstdlib> - -#if CFG_OSAPI == CFG_OSAPI_POSIX // These are POSIX headers (i.e. non-windows) - #include <sys/time.h> -#endif // POSIX - -#ifdef __APPLE__ - #include <Availability.h> - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - #define OS_OBJECT_USE_OBJC 0 // Don't make dispatch and xpc objects Objective-C objects. - #include <IOKit/pwr_mgt/IOPMLib.h> // To disable user idle sleep on Mountain Lion - #endif -#endif - -#include <time.h> -#include <list> -#include <vector> -#include <algorithm> - -#include <utility> -#include <map> - -#include <getopt.h> - -#if (CFG_COMPILER == CFG_COMPILER_MSVC) - #define FFTW_DLL -#endif - -#include <fftw3.h> -#include <sys/stat.h> -#include <sys/types.h> - -#define _xmds_malloc fftw_malloc -#define xmds_free fftw_free - -#define H5_USE_16_API -#include <hdf5.h> - -#if !defined(HAVE_H5LEXISTS) -htri_t H5Lexists(hid_t loc_id, const char *name, hid_t lapl_id) -{ - H5E_auto_t error_func; - void* error_client_data; - // Squelch errors generated by H5Gget_objinfo. It will report errors when it can't find an object - // but that's the purpose of calling it. - H5Eget_auto(&error_func, &error_client_data); - H5Eset_auto(NULL, NULL); - herr_t err = H5Gget_objinfo(loc_id, name, false, NULL); - H5Eset_auto(error_func, error_client_data); - if (err >= 0) - return true; - else - return false; -} -#endif - -#define H5T_NATIVE_REAL H5T_NATIVE_DOUBLE -#if defined(HAVE_HDF5_HL) - #include <hdf5_hl.h> -#endif - - -typedef long integer; -typedef double real; -typedef std::complex<real> XMDSComplexType; - -#include <xpdeint.h> - -#define complex XMDSComplexType - -const complex i(0.0, 1.0); - -using namespace std; - -#if CFG_COMPILER == CFG_COMPILER_ICC - // - // Disable ICC's warning: label was declared but never referenced - // - #pragma warning ( disable : 177 ) -#endif - -inline void *xmds_malloc(size_t size); - -// ******************************************************** -// DEFINES -// ******************************************************** - -// ******************************************************** -// Simulation defines -#define _EPSILON 1e-6 -#ifndef INFINITY -#define INFINITY HUGE_VAL -#endif - -#ifndef MAX -#define MAX(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a > _b ? _a : _b; }) -#endif - -#ifndef MIN -#define MIN(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a < _b ? _a : _b; }) -#endif - - -// ******************************************************** -// Geometry defines -#define _lattice_detuning ((int)256) -#define _min_detuning ((real)domain_min) -#define _max_detuning ((real)domain_max) -#define _ddetuning ((real)((_max_detuning - _min_detuning)/_lattice_detuning)) - -#define _lattice_kdetuning ((int)256) -#define _dkdetuning (2.0*M_PI/(_max_detuning - _min_detuning)) -#define _min_kdetuning (-(_lattice_kdetuning/2) * _dkdetuning) -#define _max_kdetuning ((_lattice_kdetuning - 1)/2 * _dkdetuning) - -// ******************************************************** -// field detuning defines -#define _detuning_ndims 1 - - -// vector Gamma defines -#define _detuning_Gamma_ncomponents 2 - -// vector population defines -#define _detuning_population_ncomponents 6 - -// ******************************************************** -// field mg0_sampling defines -#define _mg0_sampling_ndims 1 - - -// ******************************************************** -// field mg0_output defines -#define _mg0_output_ndims 2 - - -#define _mg0_output_lattice_t ((int)11) -#define _mg0_output_min_t (_mg0_output_t[0]) -#define _mg0_output_max_t (_mg0_output_t[_mg0_output_lattice_t-1]) -#define _mg0_output_dt (_mg0_output_t[_index_t+1]-_mg0_output_t[_index_t]) - -// vector raw defines -#define _mg0_output_raw_ncomponents 1 - - -// ******************************************************** -// GLOBALS -// ******************************************************** - - -// ******************************************************** -// Simulation globals - -string gsArgsAndValues = ""; - -real t; - -// ******************************************************** -// Transform Multiplexer globals -typedef pair<ptrdiff_t, ptrdiff_t> _basis_pair; -typedef void (*transform_function)(bool, real, real* const __restrict__, real* const __restrict__, ptrdiff_t, ptrdiff_t); - -// Less than operator needed by the C++ map class -struct _basis_pair_less_than -{ - bool operator()(const _basis_pair& _x, const _basis_pair& _y) const { - return (_x.first < _y.first) || ((_x.first == _y.first) && (_x.second < _y.second)); - } -}; - -struct _transform_step -{ - transform_function _func; - bool _forward; - bool _out_of_place; - ptrdiff_t _prefix_lattice; - ptrdiff_t _postfix_lattice; -}; - -// Structure to hold the basis change information -struct _basis_transform_t -{ - vector<_transform_step> _transform_steps; - real _multiplier; - - _basis_transform_t(real _multiplier_in = 1.0) : _multiplier(_multiplier_in) {} - - _basis_transform_t(const _basis_transform_t& _b) : _transform_steps(_b._transform_steps), _multiplier(_b._multiplier) {} - - void append(transform_function _func, bool _forward, bool _out_of_place, ptrdiff_t _prefix_lattice, ptrdiff_t _postfix_lattice) - { - _transform_steps.push_back((_transform_step){_func, _forward, _out_of_place, _prefix_lattice, _postfix_lattice}); - } -}; - -// Map type for holding (old_basis, new_basis) -> _basis_transform_t mappings -typedef map<_basis_pair, _basis_transform_t, _basis_pair_less_than> _basis_map; - - -real *_auxiliary_array = NULL; - -const char *_basis_identifiers[] = { -}; - -// ******************************************************** -// 'Globals' element globals - -#line 17 "./eit1.xmds" - -real drive; -real decay; -real decay_bc; -complex r_ca; -real probe = .0001; -real delta1; -real dephase_bc; - - -#line 271 "./eit1.cc" - -// ******************************************************** -// Command line argument processing globals -real drive_arg = 0.1; -real decay_arg = 10.0; -real decay_bc_arg = 0.001; -real domain_min = -1; -real domain_max = 1; -real delta1_arg = 0; -real dephase_bc_arg = 0; -real final_time = 10000; - -// ******************************************************** -// Geometry globals -real* _detuning = NULL; - -real* _kdetuning = NULL; - -// ******************************************************** -// FFTW3 globals -const real _inverse_sqrt_2pi = 1.0 / sqrt(2.0 * M_PI); -string _fftwWisdomPath; - -// ******************************************************** -// field detuning globals -// vector Gamma globals -size_t _detuning_Gamma_alloc_size = 0; -complex* _detuning_Gamma = NULL; -complex* _active_detuning_Gamma = NULL; - -// vector population globals -size_t _detuning_population_alloc_size = 0; -complex* _detuning_population = NULL; -complex* _active_detuning_population = NULL; - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) globals -complex* _segment1_akafield_detuning_population; -complex* _segment1_akbfield_detuning_population; -complex* _segment1_akcfield_detuning_population; -complex* _segment1_akdfield_detuning_population; -complex* _segment1_akefield_detuning_population; -complex* _segment1_akffield_detuning_population; -complex* _segment1_akgfield_detuning_population; -complex* _segment1_akhfield_detuning_population; -complex* _segment1_akifield_detuning_population; -complex* _segment1_akjfield_detuning_population; -complex* _segment1_initial_detuning_population; - -// ******************************************************** -// field mg0_output globals -real* _mg0_output_t = NULL; -unsigned long _mg0_output_index_t = 0; - -// vector raw globals -size_t _mg0_output_raw_alloc_size = 0; -real* _mg0_output_raw = NULL; -real* _active_mg0_output_raw = NULL; - - -// ******************************************************** -// FUNCTION PROTOTYPES -// ******************************************************** - -// ******************************************************** -// Command line argument processing function prototypes -void _print_usage(); - -// ******************************************************** -// field detuning function prototypes -void _detuning_Gamma_initialise(); - -void _detuning_population_initialise(); - -// ******************************************************** -// segment 0 (Top level sequence) function prototypes -void _segment0(); - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function prototypes -inline void _segment1_calculate_delta_a(real _step); -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent); -void _segment1(); -inline void _segment1_ip_evolve(int _exponent); -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter); -real _segment1_detuning_population_timestep_error(complex* _checkfield); -bool _segment1_detuning_population_reset(complex* _reset_to); - -void _segment1_detuning_operators_evaluate_operator0(real _step); - -// ******************************************************** -// output function prototypes -void _write_output(); - -FILE* _open_xsil_file(const char* _filename); -void _close_xsil_file(FILE*& fp); -void _write_xsil_header(FILE* fp); -void _write_xsil_footer(FILE* fp); - -// ******************************************************** -// moment group 0 function prototypes -void _mg0_sample(); -void _mg0_process(); -void _mg0_write_out(FILE* _outfile); - -// ******************************************************** -// field mg0_output function prototypes -void _mg0_output_raw_initialise(); - -// ******************************************************** -// MAIN ROUTINE -// ******************************************************** -int main(int argc, char **argv) -{ - #ifdef __APPLE__ - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - { - IOPMAssertionID _powerAssertionID = 0; - IOReturn __io_result = IOPMAssertionCreateWithDescription( - kIOPMAssertionTypePreventUserIdleSystemSleep, - CFSTR("XMDS simulation './eit1' preventing user idle sleep"), // Assertion name - NULL, // Assertion details - NULL, // Human-readable reason - NULL, // Localization bundle path - (CFTimeInterval)0, // never timeout - kIOPMAssertionTimeoutActionRelease, - &_powerAssertionID - ); - if (__io_result != kIOReturnSuccess) { - _LOG(_WARNING_LOG_LEVEL, "Failed to disable user idle sleep\n"); - } - // Note, this power assertion is automatically released when the process quits. - } - #endif - #endif - - - // *********** Parse the command line for arguments, and set ********* - // *********** the appropriate global variables ********* - - int resp; - std::map<string, string> mInputArgsAndValues; - - while (1) { - static struct option long_options[] = - { - {"help", no_argument, 0, 'h'}, - {"drive_arg", required_argument, 0, 'd'}, - {"decay_arg", required_argument, 0, 'e'}, - {"decay_bc_arg", required_argument, 0, 'c'}, - {"domain_min", required_argument, 0, 'o'}, - {"domain_max", required_argument, 0, 'm'}, - {"delta1_arg", required_argument, 0, 'l'}, - {"dephase_bc_arg", required_argument, 0, 'p'}, - {"final_time", required_argument, 0, 'f'}, - {NULL, 0, 0, 0} - }; - - int option_index = 0; - - resp = getopt_long(argc, argv, "hd:e:c:o:m:l:p:f:", long_options, &option_index); - - if (resp == -1) - break; - - switch (resp) { - case '?': - // An unknown option was passed. Show allowed options and exit. - _print_usage(); // This causes the simulation to exit - - case 'h': - _print_usage(); // This causes the simulation to exit - - case 'd': - drive_arg = strtod(optarg, NULL); - break; - - case 'e': - decay_arg = strtod(optarg, NULL); - break; - - case 'c': - decay_bc_arg = strtod(optarg, NULL); - break; - - case 'o': - domain_min = strtod(optarg, NULL); - break; - - case 'm': - domain_max = strtod(optarg, NULL); - break; - - case 'l': - delta1_arg = strtod(optarg, NULL); - break; - - case 'p': - dephase_bc_arg = strtod(optarg, NULL); - break; - - case 'f': - final_time = strtod(optarg, NULL); - break; - - default: - _LOG(_ERROR_LOG_LEVEL, "Internal error in processing arguments.\n"); - } - } - - - if (optind < argc) - _print_usage(); // This causes the simulation to exit. - - // ******** Argument post-processing code ******* - #line 38 "./eit1.xmds" - - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - - #line 496 "./eit1.cc" - // ********************************************** - - - - _mg0_output_raw_alloc_size = MAX(_mg0_output_raw_alloc_size, (_mg0_output_lattice_t * _lattice_detuning) * _mg0_output_raw_ncomponents); - _detuning_Gamma_alloc_size = MAX(_detuning_Gamma_alloc_size, (_lattice_detuning) * _detuning_Gamma_ncomponents); - _detuning_population_alloc_size = MAX(_detuning_population_alloc_size, (_lattice_detuning) * _detuning_population_ncomponents); - _detuning = (real*) xmds_malloc(sizeof(real) * (_lattice_detuning+1)); - - _kdetuning = (real*) xmds_malloc(sizeof(real) * (_lattice_kdetuning+1)); - - _detuning_Gamma = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_Gamma_alloc_size,1)); - _active_detuning_Gamma = _detuning_Gamma; - - - _detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _active_detuning_population = _detuning_population; - _mg0_output_t = (real*) xmds_malloc(sizeof(real) * (_mg0_output_lattice_t+1)); - - - _mg0_output_raw = (real*) xmds_malloc(sizeof(real) * MAX(_mg0_output_raw_alloc_size,1)); - _active_mg0_output_raw = _mg0_output_raw; - - - // Run-time validation checks - - if (domain_min >= domain_max) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The end point of the dimension 'domain_max' must be " - "greater than the start point.\n" - "Start = %e, End = %e\n", (real)domain_min, (real)domain_max); - - if (final_time <= 0.0) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The interval for segment 1 is not positive!\n" - "Interval = %e\n", final_time); - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) - _detuning[_index_detuning] = _min_detuning + _index_detuning*_ddetuning; - for (long _index_kdetuning = 0; _index_kdetuning < (_lattice_kdetuning+1)/2; _index_kdetuning++) - _kdetuning[_index_kdetuning] = _index_kdetuning*_dkdetuning; - for (long _index_kdetuning = (_lattice_kdetuning+1)/2; _index_kdetuning < _lattice_kdetuning; _index_kdetuning++) - _kdetuning[_index_kdetuning] = -(_lattice_kdetuning - _index_kdetuning) * _dkdetuning; - _active_mg0_output_raw = _mg0_output_raw; - _mg0_output_raw_initialise(); - // load wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX // Don't load wisdom on windows - { - char _hostName[256]; - gethostname(_hostName, 256); - _hostName[255] = '\0'; // just in case - - string _pathToWisdom = getenv("HOME"); - _pathToWisdom += "/.xmds/wisdom/"; - - string _wisdomFileName = _hostName; - _wisdomFileName += ".wisdom"; - _wisdomFileName += ".fftw3"; - - FILE *_fp = NULL; - - _fp = fopen(_pathToWisdom.c_str(), "r"); - if (_fp) { - fclose(_fp); - } else { - int _result = mkdir((string(getenv("HOME")) + "/.xmds").c_str(), S_IRWXU); - if (mkdir(_pathToWisdom.c_str(), S_IRWXU)) { - // We failed to create the ~/.xmds/wisdom directory - _LOG(_WARNING_LOG_LEVEL, "Warning: Cannot find enlightenment, the path to wisdom ~/.xmds/wisdom doesn't seem to exist and we couldn't create it.\n" - " I'll use the current path instead.\n"); - _pathToWisdom = ""; // present directory - } - - } - - _fftwWisdomPath = _pathToWisdom + _wisdomFileName; - - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "r")) != NULL) { - _LOG(_SIMULATION_LOG_LEVEL, "Found enlightenment... (Importing wisdom)\n"); - fftw_import_wisdom_from_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - _basis_transform_t *_basis_transform = NULL; - ptrdiff_t _auxiliary_array_size = 0; - ptrdiff_t _max_vector_size = 0; - real* _max_vector_array = NULL; - - if (_auxiliary_array_size) { - _auxiliary_array = (real*) xmds_malloc(sizeof(real) * _auxiliary_array_size); - } - - bool _allocated_temporary_array = false; - if (!_max_vector_array && _max_vector_size > 0) { - _max_vector_array = (real*) xmds_malloc(sizeof(real) * _max_vector_size); - _allocated_temporary_array = true; - } - - // Make all geometry-dependent transformations prepare plans, etc. - - if (_allocated_temporary_array) { - xmds_free(_max_vector_array); - } - - /* Code that actually does stuff goes here */ - _segment0(); - - - _write_output(); - if (_auxiliary_array) { - xmds_free(_auxiliary_array); - } - - // Save wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX - { - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "w")) != NULL) { - fftw_export_wisdom_to_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - fftw_cleanup(); - - xmds_free(_detuning_Gamma); - _active_detuning_Gamma = _detuning_Gamma = NULL; - - - xmds_free(_detuning_population); - _active_detuning_population = _detuning_population = NULL; - - xmds_free(_mg0_output_raw); - _active_mg0_output_raw = _mg0_output_raw = NULL; - - - return 0; -} - -// ******************************************************** -// FUNCTION IMPLEMENTATIONS -// ******************************************************** - -inline void *xmds_malloc(size_t size) -{ - void *retPointer = _xmds_malloc(size); - if ( !retPointer ) - _LOG(_ERROR_LOG_LEVEL, "ERROR: Couldn't allocate %zu bytes of memory!", size); - return retPointer; -} - - -// ******************************************************** -// Command line argument processing function implementations -void _print_usage() -{ - // This function does not return. - _LOG(_NO_ERROR_TERMINATE_LOG_LEVEL, "\n\nUsage: ./eit1 --drive_arg <real> --decay_arg <real> --decay_bc_arg <real> --domain_min <real> --domain_max <real> --delta1_arg <real> --dephase_bc_arg <real> --final_time <real>\n\n" - "Details:\n" - "Option\t\tType\t\tDefault value\n" - "-d, --drive_arg\treal \t\t0.1\n" - "-e, --decay_arg\treal \t\t10.0\n" - "-c, --decay_bc_arg\treal \t\t0.001\n" - "-o, --domain_min\treal \t\t-1\n" - "-m, --domain_max\treal \t\t1\n" - "-l, --delta1_arg\treal \t\t0\n" - "-p, --dephase_bc_arg\treal \t\t0\n" - "-f, --final_time\treal \t\t10000\n" - ); - // _LOG terminates the simulation. -} - -// ******************************************************** -// Default Simulation Driver function implementations -void _segment0() -{ - t = 0.0; - - _mg0_output_raw_initialise(); - _active_detuning_Gamma = _detuning_Gamma; - _detuning_Gamma_initialise(); - _active_detuning_population = _detuning_population; - _detuning_population_initialise(); - _mg0_output_index_t = 0; - _mg0_sample(); - _segment1(); - - _mg0_process(); -} - - -// ******************************************************** -// field detuning function implementations -// initialisation for vector Gamma -void _detuning_Gamma_initialise() -{ - - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 64 "./eit1.xmds" - - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - - #line 714 "./eit1.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb -} - -// initialisation for vector population -void _detuning_population_initialise() -{ - - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 75 "./eit1.xmds" - - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - - #line 756 "./eit1.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function implementations -inline void _segment1_calculate_delta_a(real _step) -{ - - - // Delta A propagation operator for field detuning - _segment1_detuning_operators_evaluate_operator0(_step); - -} - - -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent) -{ -} - - -void _segment1() -{ - real _step = final_time/(real)10; - real _old_step = _step; - real _min_step = _step; - real _max_step = _step; - long _attempted_steps = 0; - long _unsuccessful_steps = 0; - - real _tolerance = 1e-10; - - real _error, _last_norm_error = 1.0; - real _segment1_detuning_population_error; - - bool _discard = false; - bool _break_next = false; - - bool _next_sample_flag[3]; - for (long _i0 = 0; _i0 < 3; _i0++) - _next_sample_flag[_i0] = false; - - long _next_sample_counter[1]; - for (long _i0 = 0; _i0 < 1; _i0++) - _next_sample_counter[_i0] = 1; - - real _t_local = 0.0; - - real _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - - if ( (_t_local + _step)*(1.0 + _EPSILON) >= _t_break_next) { - _break_next = true; - _step = _t_break_next - _t_local; - } - - _segment1_akafield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akbfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akcfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akdfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akefield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akffield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akgfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akhfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akifield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akjfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_initial_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - complex* _akafield_detuning_population = _segment1_akafield_detuning_population; - complex* _akbfield_detuning_population = _segment1_akbfield_detuning_population; - complex* _akcfield_detuning_population = _segment1_akcfield_detuning_population; - complex* _akdfield_detuning_population = _segment1_akdfield_detuning_population; - complex* _akefield_detuning_population = _segment1_akefield_detuning_population; - complex* _akffield_detuning_population = _segment1_akffield_detuning_population; - complex* _akgfield_detuning_population = _segment1_akgfield_detuning_population; - complex* _akhfield_detuning_population = _segment1_akhfield_detuning_population; - complex* _akifield_detuning_population = _segment1_akifield_detuning_population; - complex* _akjfield_detuning_population = _segment1_akjfield_detuning_population; - complex* _initial_detuning_population = _segment1_initial_detuning_population; - - - // Runge Kutta method constants - real _a_raw[16]; - real _a[16]; - real _b[16][16]; - real _c[16]; - real _cs[16]; - real _d[16]; - - for (unsigned long _i0 = 0; _i0 < 16; _i0++) { - _a_raw[_i0] = _c[_i0] = _d[_i0] = 0.0; - for (unsigned long _i1 = 0; _i1 < 16; _i1++) - _b[_i0][_i1] = 0.0; - } - - _a_raw[1] = 0.02173913043478260869565217391304347; - _a_raw[2] = 0.09629581047800066670113001679819925; - _a_raw[3] = 0.14444371571700100005169502519729888; - _a_raw[4] = 0.52205882352941176470588235294117647; - _a_raw[5] = 0.22842443612863469578031459099794265; - _a_raw[6] = 0.54360353589933733219171338103002937; - _a_raw[7] = 0.64335664335664335664335664335664335; - _a_raw[8] = 0.48251748251748251748251748251748251; - _a_raw[9] = 0.06818181818181818181818181818181818; - _a_raw[10] = 0.25060827250608272506082725060827250; - _a_raw[11] = 0.66736715965600568968278165443304378; - _a_raw[12] = 0.85507246376811594202898550724637681; - _a_raw[13] = 0.89795918367346938775510204081632653; - _a_raw[14] = 1.0; - _a_raw[15] = 1.0; - - _a[0] = 0.0; - for (unsigned long _i0 = 1; _i0 < 16; _i0++) - _a[_i0] = _a_raw[_i0] - _a_raw[_i0 - 1]; - - _b[1][0] = 1.0/46.0; - _b[2][0] =-0.11698050118114486205818241524969622; - _b[2][1] = 0.21327631165914552875931243204789548; - _b[3][0] = 0.03611092892925025001292375629932472; - _b[3][2] = 0.10833278678775075003877126889797416; - _b[4][0] = 1.57329743908138605107331820072051125; - _b[4][2] =-5.98400943754042002888532938159655553; - _b[4][3] = 4.93277082198844574251789353381722074; - _b[5][0] = 0.05052046351120380909008334360006234; - _b[5][3] = 0.17686653884807108146683657390397612; - _b[5][4] = 0.00103743376935980522339467349390418; - _b[6][0] = 0.10543148021953768958529340893598138; - _b[6][3] =-0.16042415162569842979496486916719383; - _b[6][4] = 0.11643956912829316045688724281285250; - _b[6][5] = 0.48215663817720491194449759844838932; - _b[7][0] = 0.07148407148407148407148407148407148; - _b[7][5] = 0.32971116090443908023196389566296464; - _b[7][6] = 0.24216141096813279233990867620960722; - _b[8][0] = 0.07162368881118881118881118881118881; - _b[8][5] = 0.32859867301674234161492268975519694; - _b[8][6] = 0.11622213117906185418927311444060725; - _b[8][7] =-0.03392701048951048951048951048951048; - _b[9][0] = 0.04861540768024729180628870095388582; - _b[9][5] = 0.03998502200331629058445317782406268; - _b[9][6] = 0.10715724786209388876739304914053506; - _b[9][7] =-0.02177735985419485163815426357369818; - _b[9][8] =-0.10579849950964443770179884616296721; - _b[10][0] =-0.02540141041535143673515871979014924; - _b[10][5] = 1.0/30.0; - _b[10][6] =-0.16404854760069182073503553020238782; - _b[10][7] = 0.03410548898794737788891414566528526; - _b[10][8] = 0.15836825014108792658008718465091487; - _b[10][9] = 0.21425115805975734472868683695127609; - _b[11][0] = 0.00584833331460742801095934302256470; - _b[11][5] =-0.53954170547283522916525526480339109; - _b[11][6] = 0.20128430845560909506500331018201158; - _b[11][7] = 0.04347222773254789483240207937678906; - _b[11][8] =-0.00402998571475307250775349983910179; - _b[11][9] = 0.16541535721570612771420482097898952; - _b[11][10] = 0.79491862412512344573322086551518180; - _b[12][0] =-0.39964965968794892497157706711861448; - _b[12][5] =-3.79096577568393158554742638116249372; - _b[12][6] =-0.40349325653530103387515807815498044; - _b[12][7] =-2.82463879530435263378049668286220715; - _b[12][8] = 1.04226892772185985533374283289821416; - _b[12][9] = 1.12510956420436603974237036536924078; - _b[12][10] = 3.32746188718986816186934832571938138; - _b[12][11] = 2.77897957186355606325818219255783627; - _b[13][0] = 0.39545306350085237157098218205756922; - _b[13][5] = 5.82534730759650564865380791881446903; - _b[13][6] =-0.36527452339161313311889856846974452; - _b[13][7] = 1.18860324058346533283780076203192232; - _b[13][8] = 0.57970467638357921347110271762687972; - _b[13][9] =-0.86824862589087693262676988867897834; - _b[13][10] =-5.20227677296454721392873650976792184; - _b[13][11] =-0.79895541420753382543211121058675915; - _b[13][12] = 0.14360623206363792632792463778889008; - _b[14][0] = 8.49173149061346398013352206978380938; - _b[14][5] = 86.32213734729036800877634194386790750; - _b[14][6] = 1.02560575501091662034511526187393241; - _b[14][7] = 85.77427969817339941806831550695235092; - _b[14][8] =-13.98699305104110611795532466113248067; - _b[14][9] =-20.71537405501426352265946477613161883; - _b[14][10] =-72.16597156619946800281180102605140463; - _b[14][11] =-76.71211139107806345587696023064419687; - _b[14][12] = 4.22319427707298828839851258893735507; - _b[14][13] =-1.25649850482823521641825667745565428; - _b[15][0] =-0.42892119881959353241190195318730008; - _b[15][5] =-9.16865700950084689999297912545025359; - _b[15][6] = 1.08317616770620939241547721530003920; - _b[15][7] =-1.23501525358323653198215832293981810; - _b[15][8] =-1.21438272617593906232943856422371019; - _b[15][9] = 1.37226168507232166621351243731869914; - _b[15][10] = 9.15723239697162418155377135344394113; - _b[15][11] = 1.30616301842220047563298585480401671; - _b[15][12] =-0.25285618808937955976690569433069974; - _b[15][13] = 0.38099910799663987066763679926508552; - - _c[0] = 0.01490902081978461022483617102382552; - _c[7] =-0.20408044692054151258349120934134791; - _c[8] = 0.22901438600570447264772469337066476; - _c[9] = 0.12800558251147375669208211573729202; - _c[10] = 0.22380626846054143649770066956485937; - _c[11] = 0.39553165293700054420552389156421651; - _c[12] = 0.05416646758806981196568364538360743; - _c[13] = 0.12691439652445903685643385312168037; - _c[14] =-0.00052539244262118876455834655383035; - _c[15] = 1.0/31.0; - - _cs[0] = 0.00653047880643482012034413441159249; - _cs[7] =-2.31471038197461347517552506241529830; - _cs[8] = 0.43528227238866280799530900822377013; - _cs[9] = 0.14907947287101933118545845390618763; - _cs[10] = 0.17905535442235532311850533252768020; - _cs[11] = 2.53400872222767706921176214508820825; - _cs[12] =-0.55430437423209112896721332268159015; - _cs[13] = 0.56924788787870083224213506297615260; - _cs[14] =-0.03644749690427461198884026816573513; - _cs[15] = 1.0/31.0; - - _d[0] = 1.0-_b[15][5]/_b[14][5]; - _d[1] = _b[15][0]-_b[14][0]*_b[15][5]/_b[14][5]; - _d[2] = _b[15][5]/_b[14][5]; - _d[3] = _b[15][6]-_b[14][6]*_b[15][5]/_b[14][5]; - _d[4] = _b[15][7]-_b[14][7]*_b[15][5]/_b[14][5]; - _d[5] = _b[15][8]-_b[14][8]*_b[15][5]/_b[14][5]; - _d[6] = _b[15][9]-_b[14][9]*_b[15][5]/_b[14][5]; - _d[7] = _b[15][10]-_b[14][10]*_b[15][5]/_b[14][5]; - _d[8] = _b[15][11]-_b[14][11]*_b[15][5]/_b[14][5]; - _d[9] = _b[15][12]-_b[14][12]*_b[15][5]/_b[14][5]; - _d[10] = _b[15][13]-_b[14][13]*_b[15][5]/_b[14][5]; - - do { - - do { - - - // Step 1 - - memcpy(_akafield_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - _segment1_calculate_nonconstant_ip_fields(_step, 1); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - _active_detuning_population = _akafield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - // Step 2 - - t += _a[1] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[1][0]*_akafield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 2); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-2); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(2); - - // Step 3 - - t += _a[2] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[2][0]*_akafield_detuning_population[_i0] + _b[2][1]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 3); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-3); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(3); - - // Step 4 - - t += _a[3] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[3][0]*_akafield_detuning_population[_i0] + _b[3][1]*_akbfield_detuning_population[_i0] - + _b[3][2]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 4); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-4); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(4); - - // Step 5 - - t += _a[4] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[4][0]*_akafield_detuning_population[_i0] + _b[4][1]*_akbfield_detuning_population[_i0] - + _b[4][2]*_akcfield_detuning_population[_i0] + _b[4][3]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 5); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-5); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(5); - - // Step 6 - - t += _a[5] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[5][0]*_akafield_detuning_population[_i0] + _b[5][3]*_akdfield_detuning_population[_i0] - + _b[5][4]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akifield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 6); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-6); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(6); - - // Step 7 - - t += _a[6] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _detuning_population[_i0] + _b[6][0]*_akafield_detuning_population[_i0] + _b[6][3]*_akdfield_detuning_population[_i0] - + _b[6][4]*_akefield_detuning_population[_i0] + _b[6][5]*_akifield_detuning_population[_i0]; - } - - - _active_detuning_population = _akjfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 7); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-7); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(7); - - // Step 8 - - t += _a[7] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[7][0]*_akafield_detuning_population[_i0] + _b[7][5]*_akifield_detuning_population[_i0] - + _b[7][6]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 8); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-8); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(8); - - // Step 9 - - t += _a[8] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[8][0]*_akafield_detuning_population[_i0] + _b[8][5]*_akifield_detuning_population[_i0] - + _b[8][6]*_akjfield_detuning_population[_i0]+ _b[8][7]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 9); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-9); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(9); - - // Step 10 - - t += _a[9] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[9][0]*_akafield_detuning_population[_i0] + _b[9][5]*_akifield_detuning_population[_i0] - + _b[9][6]*_akjfield_detuning_population[_i0]+ _b[9][7]*_akbfield_detuning_population[_i0]+ _b[9][8]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 10); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-10); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(10); - - // Step 11 - - t += _a[10] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[10][0]*_akafield_detuning_population[_i0] + _b[10][5]*_akifield_detuning_population[_i0] - + _b[10][6]*_akjfield_detuning_population[_i0]+ _b[10][7]*_akbfield_detuning_population[_i0] + _b[10][8]*_akcfield_detuning_population[_i0] - + _b[10][9]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 11); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-11); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(11); - - // Step 12 - - t += _a[11] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akffield_detuning_population[_i0] = _detuning_population[_i0] + _b[11][0]*_akafield_detuning_population[_i0] + _b[11][5]*_akifield_detuning_population[_i0] - + _b[11][6]*_akjfield_detuning_population[_i0] + _b[11][7]*_akbfield_detuning_population[_i0] + _b[11][8]*_akcfield_detuning_population[_i0] - + _b[11][9]*_akdfield_detuning_population[_i0] + _b[11][10]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akffield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 12); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-12); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(12); - - // Step 13 - - t += _a[12] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akgfield_detuning_population[_i0] = _detuning_population[_i0] + _b[12][0]*_akafield_detuning_population[_i0] + _b[12][5]*_akifield_detuning_population[_i0] - + _b[12][6]*_akjfield_detuning_population[_i0]+ _b[12][7]*_akbfield_detuning_population[_i0] + _b[12][8]*_akcfield_detuning_population[_i0] - + _b[12][9]*_akdfield_detuning_population[_i0] + _b[12][10]*_akefield_detuning_population[_i0] + _b[12][11]*_akffield_detuning_population[_i0]; - } - - - _active_detuning_population = _akgfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 13); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-13); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(13); - - // Step 14 - - t += _a[13] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akhfield_detuning_population[_i0] = _detuning_population[_i0] + _b[13][0]*_akafield_detuning_population[_i0] + _b[13][5]*_akifield_detuning_population[_i0] - + _b[13][6]*_akjfield_detuning_population[_i0]+ _b[13][7]*_akbfield_detuning_population[_i0] + _b[13][8]*_akcfield_detuning_population[_i0] - + _b[13][9]*_akdfield_detuning_population[_i0] + _b[13][10]*_akefield_detuning_population[_i0] + _b[13][11]*_akffield_detuning_population[_i0] - + _b[13][12]*_akgfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akhfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 14); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-14); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(14); - - // Step 15 and 16 combined to reduce memory use - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[14][0]*_akafield_detuning_population[_i0] + _b[14][5]*_akifield_detuning_population[_i0] - + _b[14][6]*_akjfield_detuning_population[_i0]+ _b[14][7]*_akbfield_detuning_population[_i0] + _b[14][8]*_akcfield_detuning_population[_i0] - + _b[14][9]*_akdfield_detuning_population[_i0] + _b[14][10]*_akefield_detuning_population[_i0] + _b[14][11]*_akffield_detuning_population[_i0] - + _b[14][12]*_akgfield_detuning_population[_i0] + _b[14][13]*_akhfield_detuning_population[_i0]; - } - - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _d[0]*_detuning_population[_i0] - + _d[1]*_akafield_detuning_population[_i0] - + _d[2]*_akifield_detuning_population[_i0] - + _d[3]*_akjfield_detuning_population[_i0] - + _d[4]*_akbfield_detuning_population[_i0] - + _d[5]*_akcfield_detuning_population[_i0] - + _d[6]*_akdfield_detuning_population[_i0] - + _d[7]*_akefield_detuning_population[_i0] - + _d[8]*_akffield_detuning_population[_i0] - + _d[9]*_akgfield_detuning_population[_i0] - + _d[10]*_akhfield_detuning_population[_i0]; - } - - - t += _a[14] * _step; - - _active_detuning_population = _akifield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - t += _a[15] * _step; - - _active_detuning_population = _akjfield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // Take full step - - // ai = a - memcpy(_initial_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - // a = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _detuning_population[_i0] += _c[0]*_akafield_detuning_population[_i0] + _c[7]*_akbfield_detuning_population[_i0] + _c[8]*_akcfield_detuning_population[_i0] - + _c[9]*_akdfield_detuning_population[_i0] + _c[10]*_akefield_detuning_population[_i0] + _c[11]*_akffield_detuning_population[_i0] - + _c[12]*_akgfield_detuning_population[_i0] + _c[13]*_akhfield_detuning_population[_i0] + _c[14]*_akifield_detuning_population[_i0] - + _c[15]*_akjfield_detuning_population[_i0]; - } - - - // a* = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akafield_detuning_population[_i0] = _initial_detuning_population[_i0] + _cs[0]*_akafield_detuning_population[_i0] + _cs[7]*_akbfield_detuning_population[_i0] - + _cs[8]*_akcfield_detuning_population[_i0] + _cs[9]*_akdfield_detuning_population[_i0] + _cs[10]*_akefield_detuning_population[_i0] - + _cs[11]*_akffield_detuning_population[_i0] + _cs[12]*_akgfield_detuning_population[_i0] + _cs[13]*_akhfield_detuning_population[_i0] - + _cs[14]*_akifield_detuning_population[_i0] + _cs[15]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _detuning_population; - - - - _error = 0.0; - - _segment1_detuning_population_error = _segment1_detuning_population_timestep_error(_akafield_detuning_population); - if (_segment1_detuning_population_error > _error) - _error = _segment1_detuning_population_error; - - _attempted_steps++; - - if (_error < _tolerance) { - _t_local += _step; - if (_step > _max_step) - _max_step = _step; - if (!_break_next && _step < _min_step) - _min_step = _step; - _discard = false; - } else { - t -= _step; - - if (_segment1_detuning_population_reset(_initial_detuning_population) == false) { - - _LOG(_WARNING_LOG_LEVEL, "WARNING: NaN present. Integration halted at t = %e.\n" - " Non-finite number in integration vector \"population\" in segment 1.\n", t); - if (_mg0_output_index_t < _mg0_output_lattice_t) - _mg0_sample(); - - goto _SEGMENT1_END; - } - - _segment1_ip_evolve(-1); - - _discard = true; - _break_next = false; - _unsuccessful_steps++; - } - - _old_step = _step; - - // Resize step - if (_error < 0.5*_tolerance || _error > _tolerance) { - const real _safetyFactor = 0.90; - real _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-0.7/9.0)) * pow(_last_norm_error, real(0.4/9.0)); - _scalingFactor = MAX(_scalingFactor, 1.0/5.0); - _scalingFactor = MIN(_scalingFactor, 7.0); - if (_error > _tolerance && _scalingFactor > 1.0) { - // If our step failed don't try and increase our step size. That would be silly. - _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-1.0/9.0)); - } - _old_step = _step; - _last_norm_error = pow(_safetyFactor/_scalingFactor*pow(_last_norm_error, real(0.4/9.0)), real(9.0/0.7)); - _step *= _scalingFactor; - } - - } while (_discard); - - if (_break_next) { - if (_next_sample_flag[0]) { - _mg0_sample(); - _next_sample_counter[0]++; - } - if (_next_sample_flag[1]) - _next_sample_flag[2] = true; - else { - _break_next = false; - _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - } - } - - if ( (_t_local + _step)*(1.0 + _EPSILON) > _t_break_next) { - _break_next = true; - _LOG(_SAMPLE_LOG_LEVEL, "Current timestep: %e\n", _old_step); - _step = _t_break_next - _t_local; - } - } while (!_next_sample_flag[2]); - - _SEGMENT1_END:; - xmds_free(_segment1_akafield_detuning_population); - xmds_free(_segment1_akbfield_detuning_population); - xmds_free(_segment1_akcfield_detuning_population); - xmds_free(_segment1_akdfield_detuning_population); - xmds_free(_segment1_akefield_detuning_population); - xmds_free(_segment1_akffield_detuning_population); - xmds_free(_segment1_akgfield_detuning_population); - xmds_free(_segment1_akhfield_detuning_population); - xmds_free(_segment1_akifield_detuning_population); - xmds_free(_segment1_akjfield_detuning_population); - xmds_free(_segment1_initial_detuning_population); - - _LOG(_SEGMENT_LOG_LEVEL, "Segment 1: minimum timestep: %e maximum timestep: %e\n", _min_step, _max_step); - _LOG(_SEGMENT_LOG_LEVEL, " Attempted %li steps, %.2f%% steps failed.\n", _attempted_steps, (100.0*_unsuccessful_steps)/_attempted_steps); -} - - -inline void _segment1_ip_evolve(int _exponent) -{ -} -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter) -{ - // The numbers of the moment groups that need to be sampled at the next sampling point. - // An entry of N+1 means "reached end of integration interval" - long _momentGroupNumbersNeedingSamplingNext[2]; - long _numberOfMomentGroupsToBeSampledNext = 1; - - long _previous_m = 1; - long _previous_M = 1; - - real _t_break_next = (real)final_time; - _momentGroupNumbersNeedingSamplingNext[0] = 1; - - // initialise all flags to false - for (long _i0 = 0; _i0 < 2; _i0++) - _next_sample_flag[_i0] = false; - - /* Check if moment group needs sampling at the same time as another already discovered sample (or the final time). - * If so, add this moment group to the to-be-sampled list. If moment group demands sampling earlier than all - * previously noted moment groups, erase all previous ones from list and set the sample time to this earlier one. - */ - if (_next_sample_counter[0] * _previous_M == _previous_m * 10) { - _momentGroupNumbersNeedingSamplingNext[_numberOfMomentGroupsToBeSampledNext] = 0; - _numberOfMomentGroupsToBeSampledNext++; - } else if (_next_sample_counter[0] * _previous_M < _previous_m * 10) { - _t_break_next = _next_sample_counter[0] * ((real)final_time) / ((real)10); - _numberOfMomentGroupsToBeSampledNext = 1; - _momentGroupNumbersNeedingSamplingNext[0] = 0; - _previous_M = 10; - _previous_m = _next_sample_counter[0]; - } - - // _momentGroupNumbersNeedingSamplingNext now contains the complete list of moment groups that need - // to be sampled at the next sampling point. Set their flags to true. - for (long _i0 = 0; _i0 < _numberOfMomentGroupsToBeSampledNext; _i0++) - _next_sample_flag[_momentGroupNumbersNeedingSamplingNext[_i0]] = true; - - return _t_break_next; -} - -real _segment1_detuning_population_timestep_error(complex* _checkfield) -{ - real _error = 1e-24; - real _temp_error = 0.0; - real _temp_mod = 0.0; - - - // Find the peak value for each component of the field - real _cutoff[_detuning_population_ncomponents]; - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) - _cutoff[_i0] = 0.0; - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - _temp_mod = mod2(_detuning_population[_detuning_population_index_pointer + _i1]); - if (_xmds_isnonfinite(_temp_mod)) - _cutoff[_i1] = INFINITY; - else if (_cutoff[_i1] < _temp_mod) - _cutoff[_i1] = _temp_mod; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) { - if (_xmds_isnonfinite(_cutoff[_i0])) - // Return an error two times the tolerance in this case because the timestep must be reduced. - return 2.0*1e-10; - _cutoff[_i0] *= 0.001; - _cutoff[_i0] *= 0.001; - } - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (mod2(_detuning_population[_detuning_population_index_pointer + _i1]) > _cutoff[_i1]) { - _temp_error = abs(_detuning_population[_detuning_population_index_pointer + _i1] - _checkfield[_detuning_population_index_pointer + _i1]) / (0.5*abs(_detuning_population[_detuning_population_index_pointer + _i1]) + 0.5*abs(_checkfield[_detuning_population_index_pointer + _i1])); - - if (_xmds_isnonfinite(_temp_error)) { - /* For _temp_error to be NaN, both the absolute value of the higher and lower order solutions - must BOTH be zero. This therefore implies that their difference is zero, and that there is no error. */ - _temp_error = 0.0; - } - - if (_error < _temp_error) // UNVECTORISABLE - _error = _temp_error; - } - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - return _error; -} - -bool _segment1_detuning_population_reset(complex* _reset_to_detuning_population) -{ - memcpy(_detuning_population, _reset_to_detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - /* return false if there's a NaN somewhere in the vector, otherwise return true */ - bool bNoNaNsPresent = true; - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (_xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Re()) - || _xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Im())) bNoNaNsPresent = false; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - return bNoNaNsPresent; -} - -// Delta A propagation operator for field detuning -void _segment1_detuning_operators_evaluate_operator0(real _step) -{ - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - complex dPca_dt; - complex dPab_dt; - complex dPaa_dt; - complex dPbb_dt; - complex dPcc_dt; - complex dPcb_dt; - - #define dt _step - - // ************* Propagation code *************** - #line 92 "./eit1.xmds" - - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - - #line 1678 "./eit1.cc" - // ********************************************** - - #undef dt - - - _active_detuning_population[_detuning_population_index_pointer + 4] = dPca_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 3] = dPab_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 0] = dPaa_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 1] = dPbb_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 2] = dPcc_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 5] = dPcb_dt * _step; - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// output function implementations -void _write_output() -{ - _LOG(_SIMULATION_LOG_LEVEL, "Generating output for ./eit1\n"); - - - char *_xsilFilename = (char*)malloc(256); - snprintf(_xsilFilename, 256, "%s.xsil", ("./eit1" + gsArgsAndValues).c_str()); - - FILE* _outfile = _open_xsil_file(_xsilFilename); - - if (_outfile) { - _write_xsil_header(_outfile); - char _dataFilename[200]; - snprintf(_dataFilename, 200, "%s.h5", ("./eit1" + gsArgsAndValues).c_str()); - - H5Fclose(H5Fcreate(_dataFilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); - } - _mg0_write_out(_outfile); - - _write_xsil_footer(_outfile); - _close_xsil_file(_outfile); - free(_xsilFilename); - _xsilFilename = NULL; - _outfile = NULL; - -} - - -FILE* _open_xsil_file(const char* _filename) -{ - - FILE* fp = fopen(_filename, "w"); - - if (fp == NULL) - // _LOG will cause the simulation to exit - _LOG(_ERROR_LOG_LEVEL, "Unable to open output file '%s'.\n" - "Exiting.\n", _filename); - - return fp; -} - -void _close_xsil_file(FILE*& fp) -{ - if (fp) - fclose(fp); - fp = NULL; - -} - -void _write_xsil_header(FILE* fp) -{ - if (!fp) - return; - fprintf(fp, "<?xml version=\"1.0\" ?><simulation xmds-version=\"2\">\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- While not strictly necessary, the following two tags are handy. -->\n"); - fprintf(fp, " <author>Hunter Rew</author>\n"); - fprintf(fp, " <description>\n"); - fprintf(fp, " Model of 3 state atom\n"); - fprintf(fp, " </description>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This element defines some constants. It can be used for other\n"); - fprintf(fp, " features as well, but we will go into that in more detail later.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <features>\n"); - fprintf(fp, " <precision> double </precision>\n"); - fprintf(fp, " <globals>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " real drive;\n"); - fprintf(fp, " real decay;\n"); - fprintf(fp, " real decay_bc;\n"); - fprintf(fp, " complex r_ca;\n"); - fprintf(fp, " real probe = .0001;\n"); - fprintf(fp, " real delta1;\n"); - fprintf(fp, " real dephase_bc;\n"); - fprintf(fp, " \n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </globals>\n"); - fprintf(fp, " <validation kind=\"run-time\"/>\n"); - fprintf(fp, " <arguments>\n"); - fprintf(fp, " <argument default_value=\"0.1\" name=\"drive_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10.0\" name=\"decay_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0.001\" name=\"decay_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"-1\" name=\"domain_min\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"1\" name=\"domain_max\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"delta1_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"dephase_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10000\" name=\"final_time\" type=\"real\"/>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " drive = drive_arg;\n"); - fprintf(fp, " decay = decay_arg;\n"); - fprintf(fp, " decay_bc = decay_bc_arg;\n"); - fprintf(fp, " delta1 = delta1_arg;\n"); - fprintf(fp, " dephase_bc = dephase_bc_arg;\n"); - fprintf(fp, " r_ca = decay - i*delta1;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </arguments>\n"); - fprintf(fp, " </features>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This part defines all of the dimensions used in the problem,\n"); - fprintf(fp, " in this case, only the dimension of 'time' is needed.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <geometry>\n"); - fprintf(fp, " <propagation_dimension> t </propagation_dimension>\n"); - fprintf(fp, " <transverse_dimensions>\n"); - fprintf(fp, " <dimension domain=\"(domain_min, domain_max)\" lattice=\"256\" name=\"detuning\"/>\n"); - fprintf(fp, " </transverse_dimensions>\n"); - fprintf(fp, " </geometry>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- A 'vector' describes the variables that we will be evolving. -->\n"); - fprintf(fp, " <vector name=\"Gamma\" type=\"complex\">\n"); - fprintf(fp, " <components> r_ab r_cb </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " r_ab = decay + i*(delta1 + detuning);\n"); - fprintf(fp, " r_cb = decay_bc + i*detuning + dephase_bc;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, " <vector dimensions=\"detuning\" name=\"population\" type=\"complex\">\n"); - fprintf(fp, " <components>\n"); - fprintf(fp, " Paa Pbb Pcc Pab Pca Pcb \n"); - fprintf(fp, " </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " Pcc = .5;\n"); - fprintf(fp, " Pbb = .5;\n"); - fprintf(fp, " Paa = Pab = Pca = Pcb = 0;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <sequence>\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " Here we define what differential equations need to be solved\n"); - fprintf(fp, " and what algorithm we want to use.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <integrate algorithm=\"ARK89\" interval=\"final_time\" tolerance=\"1e-10\">\n"); - fprintf(fp, " <samples>10</samples>\n"); - fprintf(fp, " <operators>\n"); - fprintf(fp, " <integration_vectors>population</integration_vectors>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " <dependencies>Gamma</dependencies>\n"); - fprintf(fp, " </operators>\n"); - fprintf(fp, " </integrate>\n"); - fprintf(fp, " </sequence>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- This part defines what data will be saved in the output file -->\n"); - fprintf(fp, " <output format=\"hdf5\">\n"); - fprintf(fp, " <sampling_group basis=\"detuning\" initial_sample=\"yes\">\n"); - fprintf(fp, " <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PaaR = Paa.Re();\n"); - fprintf(fp, " PaaI = Paa.Im();\n"); - fprintf(fp, " PbbR = Pbb.Re();\n"); - fprintf(fp, " PbbI = Pbb.Im();\n"); - fprintf(fp, " PccR = Pcc.Re();\n"); - fprintf(fp, " PccI = Pcc.Im();\n"); - fprintf(fp, " PabR = Pab.Re();\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " PcaR = Pca.Re();\n"); - fprintf(fp, " PcaI = Pca.Im();\n"); - fprintf(fp, " PcbR = Pcb.Re();\n"); - fprintf(fp, " PcbI = Pcb.Im();\n"); - fprintf(fp, " ]]> -->\n"); - fprintf(fp, " <moments>PabI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </sampling_group>\n"); - fprintf(fp, " </output>\n"); - - fprintf(fp, "\n<info>\n"); - fprintf(fp, "Script compiled with XMDS2 version 2.2.0 \"Out of cheese error\" (r2941)\n"); - fprintf(fp, "See http://www.xmds.org for more information.\n"); - fprintf(fp, "\nVariables that can be specified on the command line:\n"); - - fprintf(fp, " Command line argument drive_arg = %e\n", drive_arg); - - fprintf(fp, " Command line argument decay_arg = %e\n", decay_arg); - - fprintf(fp, " Command line argument decay_bc_arg = %e\n", decay_bc_arg); - - fprintf(fp, " Command line argument domain_min = %e\n", domain_min); - - fprintf(fp, " Command line argument domain_max = %e\n", domain_max); - - fprintf(fp, " Command line argument delta1_arg = %e\n", delta1_arg); - - fprintf(fp, " Command line argument dephase_bc_arg = %e\n", dephase_bc_arg); - - fprintf(fp, " Command line argument final_time = %e\n", final_time); - fprintf(fp, "</info>\n"); - -} - -// In addition to writing the footer (if 'fp' is not NULL) -// this function closes the fp file pointer. -void _write_xsil_footer(FILE* fp) -{ - if (fp) { - fprintf(fp, "</simulation>\n"); - } -} - -// ******************************************************** -// moment group 0 function implementations -void _mg0_sample() -{ - - long _mg0_output_raw_index_pointer = 0; - #define PabI _active_mg0_output_raw[_mg0_output_raw_index_pointer + 0] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // Set index pointers explicitly for (some) vectors - _mg0_output_raw_index_pointer = ( 0 - + _mg0_output_index_t * _lattice_detuning - + _index_detuning * 1 ) * _mg0_output_raw_ncomponents; - #define _SAMPLE_COMPLEX(variable) \ - variable ## R = variable.Re(); variable ## I = variable.Im(); - - // *************** Sampling code **************** - #line 131 "./eit1.xmds" - - PabI = Pab.Im(); - - #line 1960 "./eit1.cc" - // ********************************************** - - #undef _SAMPLE_COMPLEX - // Increment index pointers for vectors in field mg0_sampling (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef PabI - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - - _mg0_output_t[0 + _mg0_output_index_t++] = t; - - _LOG(_SAMPLE_LOG_LEVEL, "Sampled field (for moment group #1) at t = %e\n", t); - -} - - -void _mg0_process() -{ - // No post processing needs to be done -} - - -void _mg0_write_out(FILE* _outfile) -{ - - if (_outfile) { - fprintf(_outfile, "\n"); - fprintf(_outfile, "<XSIL Name=\"moment_group_1\">\n"); - fprintf(_outfile, " <Param Name=\"n_independent\">2</Param>\n"); - fprintf(_outfile, " <Array Name=\"variables\" Type=\"Text\">\n"); - fprintf(_outfile, " <Dim>3</Dim>\n"); - fprintf(_outfile, " <Stream><Metalink Format=\"Text\" Delimiter=\" \\n\"/>\n"); - fprintf(_outfile, "t detuning PabI \n"); - fprintf(_outfile, " </Stream>\n"); - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, " <Array Name=\"data\" Type=\"double\">\n"); - fprintf(_outfile, " <Dim>%i</Dim>\n", _mg0_output_lattice_t); - fprintf(_outfile, " <Dim>%i</Dim>\n", _lattice_detuning); - fprintf(_outfile, " <Dim>3</Dim>\n"); - } - - - char _h5Filename[200]; - snprintf(_h5Filename, 200, "%s.h5", ("./eit1" + gsArgsAndValues).c_str()); - - /* Open the file */ - hid_t hdf5_file = H5Fopen(_h5Filename, H5F_ACC_RDWR, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_WARNING_LOG_LEVEL, "Failed to open HDF5 file '%s', will try to create it.", _h5Filename); - hdf5_file = H5Fcreate(_h5Filename, H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_ERROR_LOG_LEVEL, "Failed to create HDF5 file '%s'. Bailing.", _h5Filename); - } - } - - /* Create the group for this data */ - hid_t group; - if (!H5Lexists(hdf5_file, "/1", H5P_DEFAULT)) - group = H5Gcreate(hdf5_file, "/1", H5P_DEFAULT); - else - group = H5Gopen(hdf5_file, "/1"); - - if (_outfile) { - fprintf(_outfile, " <Stream><Metalink Format=\"HDF5\" Type=\"Remote\" Group=\"/1\"/>\n"); - fprintf(_outfile, "%s.h5\n", ("./eit1" + gsArgsAndValues).c_str()); - fprintf(_outfile, " </Stream>\n"); - } - - /* Create the coordinate data sets */ - hsize_t coordinate_length; - hid_t coordinate_dataspace; - coordinate_length = _mg0_output_lattice_t; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_t; - if (!H5Lexists(hdf5_file, "/1/t", H5P_DEFAULT)) - dataset_t = H5Dcreate(hdf5_file, "/1/t", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_t = H5Dopen(hdf5_file, "/1/t"); - H5Dwrite(dataset_t, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _mg0_output_t); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_t, "t"); - #endif - - H5Sclose(coordinate_dataspace); - coordinate_length = _lattice_detuning; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_detuning; - if (!H5Lexists(hdf5_file, "/1/detuning", H5P_DEFAULT)) - dataset_detuning = H5Dcreate(hdf5_file, "/1/detuning", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_detuning = H5Dopen(hdf5_file, "/1/detuning"); - H5Dwrite(dataset_detuning, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _detuning); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_detuning, "detuning"); - #endif - - H5Sclose(coordinate_dataspace); - - hsize_t file_dims[] = {_mg0_output_lattice_t, _lattice_detuning}; - hid_t file_dataspace = H5Screate_simple(2, file_dims, NULL); - - hid_t dataset_PabI; - if (!H5Lexists(hdf5_file, "/1/PabI", H5P_DEFAULT)) - dataset_PabI = H5Dcreate(hdf5_file, "/1/PabI", H5T_NATIVE_REAL, file_dataspace, H5P_DEFAULT); - else - dataset_PabI = H5Dopen(hdf5_file, "/1/PabI"); - #if defined(HAVE_HDF5_HL) - H5DSattach_scale(dataset_PabI, dataset_t, 0); - H5DSattach_scale(dataset_PabI, dataset_detuning, 1); - #endif - H5Dclose(dataset_t); - H5Dclose(dataset_detuning); - - - if ((_mg0_output_lattice_t * _lattice_detuning)) { - /* Create the data space */ - hsize_t file_start[2] = {0, 0}; - hsize_t mem_dims[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - hsize_t mem_start[3] = {0, 0, 0}; - hsize_t mem_stride[3] = {1, 1, 1}; - hsize_t mem_count[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - - - hid_t mem_dataspace; - mem_dims[2] = 1; - mem_dataspace = H5Screate_simple(3, mem_dims, NULL); - mem_stride[2] = 1; - - // Select hyperslabs of memory and file data spaces for data transfer operation - mem_start[2] = 0; - H5Sselect_hyperslab(mem_dataspace, H5S_SELECT_SET, mem_start, mem_stride, mem_count, NULL); - H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, file_start, mem_stride, mem_count, NULL); - - if (dataset_PabI) - H5Dwrite(dataset_PabI, H5T_NATIVE_REAL, mem_dataspace, file_dataspace, H5P_DEFAULT, _active_mg0_output_raw); - - H5Sclose(mem_dataspace); - } - - - H5Dclose(dataset_PabI); - - H5Sclose(file_dataspace); - H5Gclose(group); - H5Fclose(hdf5_file); - - - if (_outfile) { - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, "</XSIL>\n"); - } -} - -// ******************************************************** -// field mg0_output function implementations -// initialisation for vector raw -void _mg0_output_raw_initialise() -{ - - bzero(_active_mg0_output_raw, sizeof(real) * _mg0_output_raw_alloc_size); -} - diff --git a/source/eit1.h5 b/source/eit1.h5 Binary files differdeleted file mode 100644 index 17cfda9..0000000 --- a/source/eit1.h5 +++ /dev/null diff --git a/source/eit1.py b/source/eit1.py deleted file mode 100644 index 0f5c1a6..0000000 --- a/source/eit1.py +++ /dev/null @@ -1,17 +0,0 @@ -#!/usr/bin/env python -from xpdeint.XSILFile import XSILFile - -xsilFile = XSILFile("./eit1.xsil") - -def firstElementOrNone(enumerable): - for element in enumerable: - return element - return None - -t_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "t") -detuning_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "detuning") -PabI_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].dependentVariables if _["name"] == "PabI") - -# Write your plotting commands here. -# You may want to import pylab (from pylab import *) or matplotlib (from matplotlib import *) - diff --git a/source/eit1.pyc b/source/eit1.pyc Binary files differdeleted file mode 100644 index 509c401..0000000 --- a/source/eit1.pyc +++ /dev/null diff --git a/source/eit1.xmds b/source/eit1.xmds deleted file mode 100644 index 50419f5..0000000 --- a/source/eit1.xmds +++ /dev/null @@ -1,136 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument name="drive_arg" type="real" default_value="0.1" /> - <argument name="decay_arg" type="real" default_value="10.0" /> - <argument name="decay_bc_arg" type="real" default_value="0.001" /> - <argument name="domain_min" type="real" default_value="-1" /> - <argument name="domain_max" type="real" default_value="1" /> - <argument name="delta1_arg" type="real" default_value="0" /> - <argument name="dephase_bc_arg" type="real" default_value="0" /> - <argument name="final_time" type="real" default_value="10000" /> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension name="detuning" lattice="256" domain="(domain_min, domain_max)" /> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector name="population" type="complex" dimensions="detuning"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5" > - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> -</simulation> diff --git a/source/eit1.xsil b/source/eit1.xsil deleted file mode 100644 index 10b9bd1..0000000 --- a/source/eit1.xsil +++ /dev/null @@ -1,168 +0,0 @@ -<?xml version="1.0" ?><simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument default_value="0.1" name="drive_arg" type="real"/> - <argument default_value="10.0" name="decay_arg" type="real"/> - <argument default_value="0.001" name="decay_bc_arg" type="real"/> - <argument default_value="-1" name="domain_min" type="real"/> - <argument default_value="1" name="domain_max" type="real"/> - <argument default_value="0" name="delta1_arg" type="real"/> - <argument default_value="0" name="dephase_bc_arg" type="real"/> - <argument default_value="10000" name="final_time" type="real"/> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension domain="(domain_min, domain_max)" lattice="256" name="detuning"/> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector dimensions="detuning" name="population" type="complex"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5"> - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> - -<info> -Script compiled with XMDS2 version 2.2.0 "Out of cheese error" (r2941) -See http://www.xmds.org for more information. - -Variables that can be specified on the command line: - Command line argument drive_arg = 1.000000e+00 - Command line argument decay_arg = 6.000000e+00 - Command line argument decay_bc_arg = 1.000000e-05 - Command line argument domain_min = -1.666767e+00 - Command line argument domain_max = 1.666767e+00 - Command line argument delta1_arg = -1.500000e+01 - Command line argument dephase_bc_arg = 0.000000e+00 - Command line argument final_time = 5.999640e+02 -</info> - -<XSIL Name="moment_group_1"> - <Param Name="n_independent">2</Param> - <Array Name="variables" Type="Text"> - <Dim>3</Dim> - <Stream><Metalink Format="Text" Delimiter=" \n"/> -t detuning PabI - </Stream> - </Array> - <Array Name="data" Type="double"> - <Dim>11</Dim> - <Dim>256</Dim> - <Dim>3</Dim> - <Stream><Metalink Format="HDF5" Type="Remote" Group="/1"/> -./eit1.h5 - </Stream> - </Array> -</XSIL> -</simulation> diff --git a/source/eit2 b/source/eit2 Binary files differdeleted file mode 100755 index 0e07940..0000000 --- a/source/eit2 +++ /dev/null diff --git a/source/eit2.cc b/source/eit2.cc deleted file mode 100644 index 06001da..0000000 --- a/source/eit2.cc +++ /dev/null @@ -1,2129 +0,0 @@ -// ******************************************************** -// simulation logging - -#define _SAMPLE_LOG_LEVEL (1 << 0) -#define _SEGMENT_LOG_LEVEL (1 << 1) -#define _PATH_LOG_LEVEL (1 << 2) -#define _SIMULATION_LOG_LEVEL (1 << 3) -#define _WARNING_LOG_LEVEL (1 << 4) -#define _ERROR_LOG_LEVEL (1 << 5) -#define _NO_ERROR_TERMINATE_LOG_LEVEL (1 << 6) -#define _ALL_LOG_LEVELS _SAMPLE_LOG_LEVEL|_SEGMENT_LOG_LEVEL|_PATH_LOG_LEVEL|_SIMULATION_LOG_LEVEL|_WARNING_LOG_LEVEL|_ERROR_LOG_LEVEL|_NO_ERROR_TERMINATE_LOG_LEVEL -#define _LOG_LEVELS_BEING_LOGGED (_ALL_LOG_LEVELS) - -#define real Re -#define imag Im - -#include <complex> - -#undef real -#undef imag - - -#include <stdio.h> - -#define _LOG(logLevel, ...) \ - do { \ - if (logLevel & _LOG_LEVELS_BEING_LOGGED) { \ - if (logLevel & (_ERROR_LOG_LEVEL | _WARNING_LOG_LEVEL)) \ - printf("%s:%i: ", __FILE__, __LINE__); \ - printf(__VA_ARGS__); \ - fflush(stdout); \ - if (logLevel & (_ERROR_LOG_LEVEL | _NO_ERROR_TERMINATE_LOG_LEVEL)) \ - exit(logLevel == _ERROR_LOG_LEVEL); \ - } \ - } while (0) - -// ******************************************************** -// simulation includes - -#include <xpdeint_platform.h> -#include <cmath> -#include <string> -#include <cstring> -#include <fstream> -#include <sstream> -#include <cstdlib> - -#if CFG_OSAPI == CFG_OSAPI_POSIX // These are POSIX headers (i.e. non-windows) - #include <sys/time.h> -#endif // POSIX - -#ifdef __APPLE__ - #include <Availability.h> - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - #define OS_OBJECT_USE_OBJC 0 // Don't make dispatch and xpc objects Objective-C objects. - #include <IOKit/pwr_mgt/IOPMLib.h> // To disable user idle sleep on Mountain Lion - #endif -#endif - -#include <time.h> -#include <list> -#include <vector> -#include <algorithm> - -#include <utility> -#include <map> - -#include <getopt.h> - -#if (CFG_COMPILER == CFG_COMPILER_MSVC) - #define FFTW_DLL -#endif - -#include <fftw3.h> -#include <sys/stat.h> -#include <sys/types.h> - -#define _xmds_malloc fftw_malloc -#define xmds_free fftw_free - -#define H5_USE_16_API -#include <hdf5.h> - -#if !defined(HAVE_H5LEXISTS) -htri_t H5Lexists(hid_t loc_id, const char *name, hid_t lapl_id) -{ - H5E_auto_t error_func; - void* error_client_data; - // Squelch errors generated by H5Gget_objinfo. It will report errors when it can't find an object - // but that's the purpose of calling it. - H5Eget_auto(&error_func, &error_client_data); - H5Eset_auto(NULL, NULL); - herr_t err = H5Gget_objinfo(loc_id, name, false, NULL); - H5Eset_auto(error_func, error_client_data); - if (err >= 0) - return true; - else - return false; -} -#endif - -#define H5T_NATIVE_REAL H5T_NATIVE_DOUBLE -#if defined(HAVE_HDF5_HL) - #include <hdf5_hl.h> -#endif - - -typedef long integer; -typedef double real; -typedef std::complex<real> XMDSComplexType; - -#include <xpdeint.h> - -#define complex XMDSComplexType - -const complex i(0.0, 1.0); - -using namespace std; - -#if CFG_COMPILER == CFG_COMPILER_ICC - // - // Disable ICC's warning: label was declared but never referenced - // - #pragma warning ( disable : 177 ) -#endif - -inline void *xmds_malloc(size_t size); - -// ******************************************************** -// DEFINES -// ******************************************************** - -// ******************************************************** -// Simulation defines -#define _EPSILON 1e-6 -#ifndef INFINITY -#define INFINITY HUGE_VAL -#endif - -#ifndef MAX -#define MAX(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a > _b ? _a : _b; }) -#endif - -#ifndef MIN -#define MIN(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a < _b ? _a : _b; }) -#endif - - -// ******************************************************** -// Geometry defines -#define _lattice_detuning ((int)256) -#define _min_detuning ((real)domain_min) -#define _max_detuning ((real)domain_max) -#define _ddetuning ((real)((_max_detuning - _min_detuning)/_lattice_detuning)) - -#define _lattice_kdetuning ((int)256) -#define _dkdetuning (2.0*M_PI/(_max_detuning - _min_detuning)) -#define _min_kdetuning (-(_lattice_kdetuning/2) * _dkdetuning) -#define _max_kdetuning ((_lattice_kdetuning - 1)/2 * _dkdetuning) - -// ******************************************************** -// field detuning defines -#define _detuning_ndims 1 - - -// vector Gamma defines -#define _detuning_Gamma_ncomponents 2 - -// vector population defines -#define _detuning_population_ncomponents 6 - -// ******************************************************** -// field mg0_sampling defines -#define _mg0_sampling_ndims 1 - - -// ******************************************************** -// field mg0_output defines -#define _mg0_output_ndims 2 - - -#define _mg0_output_lattice_t ((int)11) -#define _mg0_output_min_t (_mg0_output_t[0]) -#define _mg0_output_max_t (_mg0_output_t[_mg0_output_lattice_t-1]) -#define _mg0_output_dt (_mg0_output_t[_index_t+1]-_mg0_output_t[_index_t]) - -// vector raw defines -#define _mg0_output_raw_ncomponents 1 - - -// ******************************************************** -// GLOBALS -// ******************************************************** - - -// ******************************************************** -// Simulation globals - -string gsArgsAndValues = ""; - -real t; - -// ******************************************************** -// Transform Multiplexer globals -typedef pair<ptrdiff_t, ptrdiff_t> _basis_pair; -typedef void (*transform_function)(bool, real, real* const __restrict__, real* const __restrict__, ptrdiff_t, ptrdiff_t); - -// Less than operator needed by the C++ map class -struct _basis_pair_less_than -{ - bool operator()(const _basis_pair& _x, const _basis_pair& _y) const { - return (_x.first < _y.first) || ((_x.first == _y.first) && (_x.second < _y.second)); - } -}; - -struct _transform_step -{ - transform_function _func; - bool _forward; - bool _out_of_place; - ptrdiff_t _prefix_lattice; - ptrdiff_t _postfix_lattice; -}; - -// Structure to hold the basis change information -struct _basis_transform_t -{ - vector<_transform_step> _transform_steps; - real _multiplier; - - _basis_transform_t(real _multiplier_in = 1.0) : _multiplier(_multiplier_in) {} - - _basis_transform_t(const _basis_transform_t& _b) : _transform_steps(_b._transform_steps), _multiplier(_b._multiplier) {} - - void append(transform_function _func, bool _forward, bool _out_of_place, ptrdiff_t _prefix_lattice, ptrdiff_t _postfix_lattice) - { - _transform_steps.push_back((_transform_step){_func, _forward, _out_of_place, _prefix_lattice, _postfix_lattice}); - } -}; - -// Map type for holding (old_basis, new_basis) -> _basis_transform_t mappings -typedef map<_basis_pair, _basis_transform_t, _basis_pair_less_than> _basis_map; - - -real *_auxiliary_array = NULL; - -const char *_basis_identifiers[] = { -}; - -// ******************************************************** -// 'Globals' element globals - -#line 17 "./eit2.xmds" - -real drive; -real decay; -real decay_bc; -complex r_ca; -real probe = .0001; -real delta1; -real dephase_bc; - - -#line 271 "./eit2.cc" - -// ******************************************************** -// Command line argument processing globals -real drive_arg = 0.1; -real decay_arg = 10.0; -real decay_bc_arg = 0.001; -real domain_min = -1; -real domain_max = 1; -real delta1_arg = 0; -real dephase_bc_arg = 0; -real final_time = 10000; - -// ******************************************************** -// Geometry globals -real* _detuning = NULL; - -real* _kdetuning = NULL; - -// ******************************************************** -// FFTW3 globals -const real _inverse_sqrt_2pi = 1.0 / sqrt(2.0 * M_PI); -string _fftwWisdomPath; - -// ******************************************************** -// field detuning globals -// vector Gamma globals -size_t _detuning_Gamma_alloc_size = 0; -complex* _detuning_Gamma = NULL; -complex* _active_detuning_Gamma = NULL; - -// vector population globals -size_t _detuning_population_alloc_size = 0; -complex* _detuning_population = NULL; -complex* _active_detuning_population = NULL; - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) globals -complex* _segment1_akafield_detuning_population; -complex* _segment1_akbfield_detuning_population; -complex* _segment1_akcfield_detuning_population; -complex* _segment1_akdfield_detuning_population; -complex* _segment1_akefield_detuning_population; -complex* _segment1_akffield_detuning_population; -complex* _segment1_akgfield_detuning_population; -complex* _segment1_akhfield_detuning_population; -complex* _segment1_akifield_detuning_population; -complex* _segment1_akjfield_detuning_population; -complex* _segment1_initial_detuning_population; - -// ******************************************************** -// field mg0_output globals -real* _mg0_output_t = NULL; -unsigned long _mg0_output_index_t = 0; - -// vector raw globals -size_t _mg0_output_raw_alloc_size = 0; -real* _mg0_output_raw = NULL; -real* _active_mg0_output_raw = NULL; - - -// ******************************************************** -// FUNCTION PROTOTYPES -// ******************************************************** - -// ******************************************************** -// Command line argument processing function prototypes -void _print_usage(); - -// ******************************************************** -// field detuning function prototypes -void _detuning_Gamma_initialise(); - -void _detuning_population_initialise(); - -// ******************************************************** -// segment 0 (Top level sequence) function prototypes -void _segment0(); - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function prototypes -inline void _segment1_calculate_delta_a(real _step); -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent); -void _segment1(); -inline void _segment1_ip_evolve(int _exponent); -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter); -real _segment1_detuning_population_timestep_error(complex* _checkfield); -bool _segment1_detuning_population_reset(complex* _reset_to); - -void _segment1_detuning_operators_evaluate_operator0(real _step); - -// ******************************************************** -// output function prototypes -void _write_output(); - -FILE* _open_xsil_file(const char* _filename); -void _close_xsil_file(FILE*& fp); -void _write_xsil_header(FILE* fp); -void _write_xsil_footer(FILE* fp); - -// ******************************************************** -// moment group 0 function prototypes -void _mg0_sample(); -void _mg0_process(); -void _mg0_write_out(FILE* _outfile); - -// ******************************************************** -// field mg0_output function prototypes -void _mg0_output_raw_initialise(); - -// ******************************************************** -// MAIN ROUTINE -// ******************************************************** -int main(int argc, char **argv) -{ - #ifdef __APPLE__ - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - { - IOPMAssertionID _powerAssertionID = 0; - IOReturn __io_result = IOPMAssertionCreateWithDescription( - kIOPMAssertionTypePreventUserIdleSystemSleep, - CFSTR("XMDS simulation './eit2' preventing user idle sleep"), // Assertion name - NULL, // Assertion details - NULL, // Human-readable reason - NULL, // Localization bundle path - (CFTimeInterval)0, // never timeout - kIOPMAssertionTimeoutActionRelease, - &_powerAssertionID - ); - if (__io_result != kIOReturnSuccess) { - _LOG(_WARNING_LOG_LEVEL, "Failed to disable user idle sleep\n"); - } - // Note, this power assertion is automatically released when the process quits. - } - #endif - #endif - - - // *********** Parse the command line for arguments, and set ********* - // *********** the appropriate global variables ********* - - int resp; - std::map<string, string> mInputArgsAndValues; - - while (1) { - static struct option long_options[] = - { - {"help", no_argument, 0, 'h'}, - {"drive_arg", required_argument, 0, 'd'}, - {"decay_arg", required_argument, 0, 'e'}, - {"decay_bc_arg", required_argument, 0, 'c'}, - {"domain_min", required_argument, 0, 'o'}, - {"domain_max", required_argument, 0, 'm'}, - {"delta1_arg", required_argument, 0, 'l'}, - {"dephase_bc_arg", required_argument, 0, 'p'}, - {"final_time", required_argument, 0, 'f'}, - {NULL, 0, 0, 0} - }; - - int option_index = 0; - - resp = getopt_long(argc, argv, "hd:e:c:o:m:l:p:f:", long_options, &option_index); - - if (resp == -1) - break; - - switch (resp) { - case '?': - // An unknown option was passed. Show allowed options and exit. - _print_usage(); // This causes the simulation to exit - - case 'h': - _print_usage(); // This causes the simulation to exit - - case 'd': - drive_arg = strtod(optarg, NULL); - break; - - case 'e': - decay_arg = strtod(optarg, NULL); - break; - - case 'c': - decay_bc_arg = strtod(optarg, NULL); - break; - - case 'o': - domain_min = strtod(optarg, NULL); - break; - - case 'm': - domain_max = strtod(optarg, NULL); - break; - - case 'l': - delta1_arg = strtod(optarg, NULL); - break; - - case 'p': - dephase_bc_arg = strtod(optarg, NULL); - break; - - case 'f': - final_time = strtod(optarg, NULL); - break; - - default: - _LOG(_ERROR_LOG_LEVEL, "Internal error in processing arguments.\n"); - } - } - - - if (optind < argc) - _print_usage(); // This causes the simulation to exit. - - // ******** Argument post-processing code ******* - #line 38 "./eit2.xmds" - - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - - #line 496 "./eit2.cc" - // ********************************************** - - - - _mg0_output_raw_alloc_size = MAX(_mg0_output_raw_alloc_size, (_mg0_output_lattice_t * _lattice_detuning) * _mg0_output_raw_ncomponents); - _detuning_Gamma_alloc_size = MAX(_detuning_Gamma_alloc_size, (_lattice_detuning) * _detuning_Gamma_ncomponents); - _detuning_population_alloc_size = MAX(_detuning_population_alloc_size, (_lattice_detuning) * _detuning_population_ncomponents); - _detuning = (real*) xmds_malloc(sizeof(real) * (_lattice_detuning+1)); - - _kdetuning = (real*) xmds_malloc(sizeof(real) * (_lattice_kdetuning+1)); - - _detuning_Gamma = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_Gamma_alloc_size,1)); - _active_detuning_Gamma = _detuning_Gamma; - - - _detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _active_detuning_population = _detuning_population; - _mg0_output_t = (real*) xmds_malloc(sizeof(real) * (_mg0_output_lattice_t+1)); - - - _mg0_output_raw = (real*) xmds_malloc(sizeof(real) * MAX(_mg0_output_raw_alloc_size,1)); - _active_mg0_output_raw = _mg0_output_raw; - - - // Run-time validation checks - - if (domain_min >= domain_max) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The end point of the dimension 'domain_max' must be " - "greater than the start point.\n" - "Start = %e, End = %e\n", (real)domain_min, (real)domain_max); - - if (final_time <= 0.0) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The interval for segment 1 is not positive!\n" - "Interval = %e\n", final_time); - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) - _detuning[_index_detuning] = _min_detuning + _index_detuning*_ddetuning; - for (long _index_kdetuning = 0; _index_kdetuning < (_lattice_kdetuning+1)/2; _index_kdetuning++) - _kdetuning[_index_kdetuning] = _index_kdetuning*_dkdetuning; - for (long _index_kdetuning = (_lattice_kdetuning+1)/2; _index_kdetuning < _lattice_kdetuning; _index_kdetuning++) - _kdetuning[_index_kdetuning] = -(_lattice_kdetuning - _index_kdetuning) * _dkdetuning; - _active_mg0_output_raw = _mg0_output_raw; - _mg0_output_raw_initialise(); - // load wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX // Don't load wisdom on windows - { - char _hostName[256]; - gethostname(_hostName, 256); - _hostName[255] = '\0'; // just in case - - string _pathToWisdom = getenv("HOME"); - _pathToWisdom += "/.xmds/wisdom/"; - - string _wisdomFileName = _hostName; - _wisdomFileName += ".wisdom"; - _wisdomFileName += ".fftw3"; - - FILE *_fp = NULL; - - _fp = fopen(_pathToWisdom.c_str(), "r"); - if (_fp) { - fclose(_fp); - } else { - int _result = mkdir((string(getenv("HOME")) + "/.xmds").c_str(), S_IRWXU); - if (mkdir(_pathToWisdom.c_str(), S_IRWXU)) { - // We failed to create the ~/.xmds/wisdom directory - _LOG(_WARNING_LOG_LEVEL, "Warning: Cannot find enlightenment, the path to wisdom ~/.xmds/wisdom doesn't seem to exist and we couldn't create it.\n" - " I'll use the current path instead.\n"); - _pathToWisdom = ""; // present directory - } - - } - - _fftwWisdomPath = _pathToWisdom + _wisdomFileName; - - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "r")) != NULL) { - _LOG(_SIMULATION_LOG_LEVEL, "Found enlightenment... (Importing wisdom)\n"); - fftw_import_wisdom_from_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - _basis_transform_t *_basis_transform = NULL; - ptrdiff_t _auxiliary_array_size = 0; - ptrdiff_t _max_vector_size = 0; - real* _max_vector_array = NULL; - - if (_auxiliary_array_size) { - _auxiliary_array = (real*) xmds_malloc(sizeof(real) * _auxiliary_array_size); - } - - bool _allocated_temporary_array = false; - if (!_max_vector_array && _max_vector_size > 0) { - _max_vector_array = (real*) xmds_malloc(sizeof(real) * _max_vector_size); - _allocated_temporary_array = true; - } - - // Make all geometry-dependent transformations prepare plans, etc. - - if (_allocated_temporary_array) { - xmds_free(_max_vector_array); - } - - /* Code that actually does stuff goes here */ - _segment0(); - - - _write_output(); - if (_auxiliary_array) { - xmds_free(_auxiliary_array); - } - - // Save wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX - { - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "w")) != NULL) { - fftw_export_wisdom_to_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - fftw_cleanup(); - - xmds_free(_detuning_Gamma); - _active_detuning_Gamma = _detuning_Gamma = NULL; - - - xmds_free(_detuning_population); - _active_detuning_population = _detuning_population = NULL; - - xmds_free(_mg0_output_raw); - _active_mg0_output_raw = _mg0_output_raw = NULL; - - - return 0; -} - -// ******************************************************** -// FUNCTION IMPLEMENTATIONS -// ******************************************************** - -inline void *xmds_malloc(size_t size) -{ - void *retPointer = _xmds_malloc(size); - if ( !retPointer ) - _LOG(_ERROR_LOG_LEVEL, "ERROR: Couldn't allocate %zu bytes of memory!", size); - return retPointer; -} - - -// ******************************************************** -// Command line argument processing function implementations -void _print_usage() -{ - // This function does not return. - _LOG(_NO_ERROR_TERMINATE_LOG_LEVEL, "\n\nUsage: ./eit2 --drive_arg <real> --decay_arg <real> --decay_bc_arg <real> --domain_min <real> --domain_max <real> --delta1_arg <real> --dephase_bc_arg <real> --final_time <real>\n\n" - "Details:\n" - "Option\t\tType\t\tDefault value\n" - "-d, --drive_arg\treal \t\t0.1\n" - "-e, --decay_arg\treal \t\t10.0\n" - "-c, --decay_bc_arg\treal \t\t0.001\n" - "-o, --domain_min\treal \t\t-1\n" - "-m, --domain_max\treal \t\t1\n" - "-l, --delta1_arg\treal \t\t0\n" - "-p, --dephase_bc_arg\treal \t\t0\n" - "-f, --final_time\treal \t\t10000\n" - ); - // _LOG terminates the simulation. -} - -// ******************************************************** -// Default Simulation Driver function implementations -void _segment0() -{ - t = 0.0; - - _mg0_output_raw_initialise(); - _active_detuning_Gamma = _detuning_Gamma; - _detuning_Gamma_initialise(); - _active_detuning_population = _detuning_population; - _detuning_population_initialise(); - _mg0_output_index_t = 0; - _mg0_sample(); - _segment1(); - - _mg0_process(); -} - - -// ******************************************************** -// field detuning function implementations -// initialisation for vector Gamma -void _detuning_Gamma_initialise() -{ - - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 64 "./eit2.xmds" - - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - - #line 714 "./eit2.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb -} - -// initialisation for vector population -void _detuning_population_initialise() -{ - - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 75 "./eit2.xmds" - - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - - #line 756 "./eit2.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function implementations -inline void _segment1_calculate_delta_a(real _step) -{ - - - // Delta A propagation operator for field detuning - _segment1_detuning_operators_evaluate_operator0(_step); - -} - - -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent) -{ -} - - -void _segment1() -{ - real _step = final_time/(real)10; - real _old_step = _step; - real _min_step = _step; - real _max_step = _step; - long _attempted_steps = 0; - long _unsuccessful_steps = 0; - - real _tolerance = 1e-10; - - real _error, _last_norm_error = 1.0; - real _segment1_detuning_population_error; - - bool _discard = false; - bool _break_next = false; - - bool _next_sample_flag[3]; - for (long _i0 = 0; _i0 < 3; _i0++) - _next_sample_flag[_i0] = false; - - long _next_sample_counter[1]; - for (long _i0 = 0; _i0 < 1; _i0++) - _next_sample_counter[_i0] = 1; - - real _t_local = 0.0; - - real _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - - if ( (_t_local + _step)*(1.0 + _EPSILON) >= _t_break_next) { - _break_next = true; - _step = _t_break_next - _t_local; - } - - _segment1_akafield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akbfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akcfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akdfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akefield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akffield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akgfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akhfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akifield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akjfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_initial_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - complex* _akafield_detuning_population = _segment1_akafield_detuning_population; - complex* _akbfield_detuning_population = _segment1_akbfield_detuning_population; - complex* _akcfield_detuning_population = _segment1_akcfield_detuning_population; - complex* _akdfield_detuning_population = _segment1_akdfield_detuning_population; - complex* _akefield_detuning_population = _segment1_akefield_detuning_population; - complex* _akffield_detuning_population = _segment1_akffield_detuning_population; - complex* _akgfield_detuning_population = _segment1_akgfield_detuning_population; - complex* _akhfield_detuning_population = _segment1_akhfield_detuning_population; - complex* _akifield_detuning_population = _segment1_akifield_detuning_population; - complex* _akjfield_detuning_population = _segment1_akjfield_detuning_population; - complex* _initial_detuning_population = _segment1_initial_detuning_population; - - - // Runge Kutta method constants - real _a_raw[16]; - real _a[16]; - real _b[16][16]; - real _c[16]; - real _cs[16]; - real _d[16]; - - for (unsigned long _i0 = 0; _i0 < 16; _i0++) { - _a_raw[_i0] = _c[_i0] = _d[_i0] = 0.0; - for (unsigned long _i1 = 0; _i1 < 16; _i1++) - _b[_i0][_i1] = 0.0; - } - - _a_raw[1] = 0.02173913043478260869565217391304347; - _a_raw[2] = 0.09629581047800066670113001679819925; - _a_raw[3] = 0.14444371571700100005169502519729888; - _a_raw[4] = 0.52205882352941176470588235294117647; - _a_raw[5] = 0.22842443612863469578031459099794265; - _a_raw[6] = 0.54360353589933733219171338103002937; - _a_raw[7] = 0.64335664335664335664335664335664335; - _a_raw[8] = 0.48251748251748251748251748251748251; - _a_raw[9] = 0.06818181818181818181818181818181818; - _a_raw[10] = 0.25060827250608272506082725060827250; - _a_raw[11] = 0.66736715965600568968278165443304378; - _a_raw[12] = 0.85507246376811594202898550724637681; - _a_raw[13] = 0.89795918367346938775510204081632653; - _a_raw[14] = 1.0; - _a_raw[15] = 1.0; - - _a[0] = 0.0; - for (unsigned long _i0 = 1; _i0 < 16; _i0++) - _a[_i0] = _a_raw[_i0] - _a_raw[_i0 - 1]; - - _b[1][0] = 1.0/46.0; - _b[2][0] =-0.11698050118114486205818241524969622; - _b[2][1] = 0.21327631165914552875931243204789548; - _b[3][0] = 0.03611092892925025001292375629932472; - _b[3][2] = 0.10833278678775075003877126889797416; - _b[4][0] = 1.57329743908138605107331820072051125; - _b[4][2] =-5.98400943754042002888532938159655553; - _b[4][3] = 4.93277082198844574251789353381722074; - _b[5][0] = 0.05052046351120380909008334360006234; - _b[5][3] = 0.17686653884807108146683657390397612; - _b[5][4] = 0.00103743376935980522339467349390418; - _b[6][0] = 0.10543148021953768958529340893598138; - _b[6][3] =-0.16042415162569842979496486916719383; - _b[6][4] = 0.11643956912829316045688724281285250; - _b[6][5] = 0.48215663817720491194449759844838932; - _b[7][0] = 0.07148407148407148407148407148407148; - _b[7][5] = 0.32971116090443908023196389566296464; - _b[7][6] = 0.24216141096813279233990867620960722; - _b[8][0] = 0.07162368881118881118881118881118881; - _b[8][5] = 0.32859867301674234161492268975519694; - _b[8][6] = 0.11622213117906185418927311444060725; - _b[8][7] =-0.03392701048951048951048951048951048; - _b[9][0] = 0.04861540768024729180628870095388582; - _b[9][5] = 0.03998502200331629058445317782406268; - _b[9][6] = 0.10715724786209388876739304914053506; - _b[9][7] =-0.02177735985419485163815426357369818; - _b[9][8] =-0.10579849950964443770179884616296721; - _b[10][0] =-0.02540141041535143673515871979014924; - _b[10][5] = 1.0/30.0; - _b[10][6] =-0.16404854760069182073503553020238782; - _b[10][7] = 0.03410548898794737788891414566528526; - _b[10][8] = 0.15836825014108792658008718465091487; - _b[10][9] = 0.21425115805975734472868683695127609; - _b[11][0] = 0.00584833331460742801095934302256470; - _b[11][5] =-0.53954170547283522916525526480339109; - _b[11][6] = 0.20128430845560909506500331018201158; - _b[11][7] = 0.04347222773254789483240207937678906; - _b[11][8] =-0.00402998571475307250775349983910179; - _b[11][9] = 0.16541535721570612771420482097898952; - _b[11][10] = 0.79491862412512344573322086551518180; - _b[12][0] =-0.39964965968794892497157706711861448; - _b[12][5] =-3.79096577568393158554742638116249372; - _b[12][6] =-0.40349325653530103387515807815498044; - _b[12][7] =-2.82463879530435263378049668286220715; - _b[12][8] = 1.04226892772185985533374283289821416; - _b[12][9] = 1.12510956420436603974237036536924078; - _b[12][10] = 3.32746188718986816186934832571938138; - _b[12][11] = 2.77897957186355606325818219255783627; - _b[13][0] = 0.39545306350085237157098218205756922; - _b[13][5] = 5.82534730759650564865380791881446903; - _b[13][6] =-0.36527452339161313311889856846974452; - _b[13][7] = 1.18860324058346533283780076203192232; - _b[13][8] = 0.57970467638357921347110271762687972; - _b[13][9] =-0.86824862589087693262676988867897834; - _b[13][10] =-5.20227677296454721392873650976792184; - _b[13][11] =-0.79895541420753382543211121058675915; - _b[13][12] = 0.14360623206363792632792463778889008; - _b[14][0] = 8.49173149061346398013352206978380938; - _b[14][5] = 86.32213734729036800877634194386790750; - _b[14][6] = 1.02560575501091662034511526187393241; - _b[14][7] = 85.77427969817339941806831550695235092; - _b[14][8] =-13.98699305104110611795532466113248067; - _b[14][9] =-20.71537405501426352265946477613161883; - _b[14][10] =-72.16597156619946800281180102605140463; - _b[14][11] =-76.71211139107806345587696023064419687; - _b[14][12] = 4.22319427707298828839851258893735507; - _b[14][13] =-1.25649850482823521641825667745565428; - _b[15][0] =-0.42892119881959353241190195318730008; - _b[15][5] =-9.16865700950084689999297912545025359; - _b[15][6] = 1.08317616770620939241547721530003920; - _b[15][7] =-1.23501525358323653198215832293981810; - _b[15][8] =-1.21438272617593906232943856422371019; - _b[15][9] = 1.37226168507232166621351243731869914; - _b[15][10] = 9.15723239697162418155377135344394113; - _b[15][11] = 1.30616301842220047563298585480401671; - _b[15][12] =-0.25285618808937955976690569433069974; - _b[15][13] = 0.38099910799663987066763679926508552; - - _c[0] = 0.01490902081978461022483617102382552; - _c[7] =-0.20408044692054151258349120934134791; - _c[8] = 0.22901438600570447264772469337066476; - _c[9] = 0.12800558251147375669208211573729202; - _c[10] = 0.22380626846054143649770066956485937; - _c[11] = 0.39553165293700054420552389156421651; - _c[12] = 0.05416646758806981196568364538360743; - _c[13] = 0.12691439652445903685643385312168037; - _c[14] =-0.00052539244262118876455834655383035; - _c[15] = 1.0/31.0; - - _cs[0] = 0.00653047880643482012034413441159249; - _cs[7] =-2.31471038197461347517552506241529830; - _cs[8] = 0.43528227238866280799530900822377013; - _cs[9] = 0.14907947287101933118545845390618763; - _cs[10] = 0.17905535442235532311850533252768020; - _cs[11] = 2.53400872222767706921176214508820825; - _cs[12] =-0.55430437423209112896721332268159015; - _cs[13] = 0.56924788787870083224213506297615260; - _cs[14] =-0.03644749690427461198884026816573513; - _cs[15] = 1.0/31.0; - - _d[0] = 1.0-_b[15][5]/_b[14][5]; - _d[1] = _b[15][0]-_b[14][0]*_b[15][5]/_b[14][5]; - _d[2] = _b[15][5]/_b[14][5]; - _d[3] = _b[15][6]-_b[14][6]*_b[15][5]/_b[14][5]; - _d[4] = _b[15][7]-_b[14][7]*_b[15][5]/_b[14][5]; - _d[5] = _b[15][8]-_b[14][8]*_b[15][5]/_b[14][5]; - _d[6] = _b[15][9]-_b[14][9]*_b[15][5]/_b[14][5]; - _d[7] = _b[15][10]-_b[14][10]*_b[15][5]/_b[14][5]; - _d[8] = _b[15][11]-_b[14][11]*_b[15][5]/_b[14][5]; - _d[9] = _b[15][12]-_b[14][12]*_b[15][5]/_b[14][5]; - _d[10] = _b[15][13]-_b[14][13]*_b[15][5]/_b[14][5]; - - do { - - do { - - - // Step 1 - - memcpy(_akafield_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - _segment1_calculate_nonconstant_ip_fields(_step, 1); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - _active_detuning_population = _akafield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - // Step 2 - - t += _a[1] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[1][0]*_akafield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 2); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-2); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(2); - - // Step 3 - - t += _a[2] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[2][0]*_akafield_detuning_population[_i0] + _b[2][1]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 3); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-3); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(3); - - // Step 4 - - t += _a[3] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[3][0]*_akafield_detuning_population[_i0] + _b[3][1]*_akbfield_detuning_population[_i0] - + _b[3][2]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 4); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-4); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(4); - - // Step 5 - - t += _a[4] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[4][0]*_akafield_detuning_population[_i0] + _b[4][1]*_akbfield_detuning_population[_i0] - + _b[4][2]*_akcfield_detuning_population[_i0] + _b[4][3]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 5); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-5); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(5); - - // Step 6 - - t += _a[5] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[5][0]*_akafield_detuning_population[_i0] + _b[5][3]*_akdfield_detuning_population[_i0] - + _b[5][4]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akifield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 6); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-6); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(6); - - // Step 7 - - t += _a[6] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _detuning_population[_i0] + _b[6][0]*_akafield_detuning_population[_i0] + _b[6][3]*_akdfield_detuning_population[_i0] - + _b[6][4]*_akefield_detuning_population[_i0] + _b[6][5]*_akifield_detuning_population[_i0]; - } - - - _active_detuning_population = _akjfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 7); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-7); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(7); - - // Step 8 - - t += _a[7] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[7][0]*_akafield_detuning_population[_i0] + _b[7][5]*_akifield_detuning_population[_i0] - + _b[7][6]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 8); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-8); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(8); - - // Step 9 - - t += _a[8] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[8][0]*_akafield_detuning_population[_i0] + _b[8][5]*_akifield_detuning_population[_i0] - + _b[8][6]*_akjfield_detuning_population[_i0]+ _b[8][7]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 9); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-9); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(9); - - // Step 10 - - t += _a[9] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[9][0]*_akafield_detuning_population[_i0] + _b[9][5]*_akifield_detuning_population[_i0] - + _b[9][6]*_akjfield_detuning_population[_i0]+ _b[9][7]*_akbfield_detuning_population[_i0]+ _b[9][8]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 10); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-10); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(10); - - // Step 11 - - t += _a[10] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[10][0]*_akafield_detuning_population[_i0] + _b[10][5]*_akifield_detuning_population[_i0] - + _b[10][6]*_akjfield_detuning_population[_i0]+ _b[10][7]*_akbfield_detuning_population[_i0] + _b[10][8]*_akcfield_detuning_population[_i0] - + _b[10][9]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 11); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-11); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(11); - - // Step 12 - - t += _a[11] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akffield_detuning_population[_i0] = _detuning_population[_i0] + _b[11][0]*_akafield_detuning_population[_i0] + _b[11][5]*_akifield_detuning_population[_i0] - + _b[11][6]*_akjfield_detuning_population[_i0] + _b[11][7]*_akbfield_detuning_population[_i0] + _b[11][8]*_akcfield_detuning_population[_i0] - + _b[11][9]*_akdfield_detuning_population[_i0] + _b[11][10]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akffield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 12); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-12); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(12); - - // Step 13 - - t += _a[12] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akgfield_detuning_population[_i0] = _detuning_population[_i0] + _b[12][0]*_akafield_detuning_population[_i0] + _b[12][5]*_akifield_detuning_population[_i0] - + _b[12][6]*_akjfield_detuning_population[_i0]+ _b[12][7]*_akbfield_detuning_population[_i0] + _b[12][8]*_akcfield_detuning_population[_i0] - + _b[12][9]*_akdfield_detuning_population[_i0] + _b[12][10]*_akefield_detuning_population[_i0] + _b[12][11]*_akffield_detuning_population[_i0]; - } - - - _active_detuning_population = _akgfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 13); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-13); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(13); - - // Step 14 - - t += _a[13] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akhfield_detuning_population[_i0] = _detuning_population[_i0] + _b[13][0]*_akafield_detuning_population[_i0] + _b[13][5]*_akifield_detuning_population[_i0] - + _b[13][6]*_akjfield_detuning_population[_i0]+ _b[13][7]*_akbfield_detuning_population[_i0] + _b[13][8]*_akcfield_detuning_population[_i0] - + _b[13][9]*_akdfield_detuning_population[_i0] + _b[13][10]*_akefield_detuning_population[_i0] + _b[13][11]*_akffield_detuning_population[_i0] - + _b[13][12]*_akgfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akhfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 14); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-14); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(14); - - // Step 15 and 16 combined to reduce memory use - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[14][0]*_akafield_detuning_population[_i0] + _b[14][5]*_akifield_detuning_population[_i0] - + _b[14][6]*_akjfield_detuning_population[_i0]+ _b[14][7]*_akbfield_detuning_population[_i0] + _b[14][8]*_akcfield_detuning_population[_i0] - + _b[14][9]*_akdfield_detuning_population[_i0] + _b[14][10]*_akefield_detuning_population[_i0] + _b[14][11]*_akffield_detuning_population[_i0] - + _b[14][12]*_akgfield_detuning_population[_i0] + _b[14][13]*_akhfield_detuning_population[_i0]; - } - - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _d[0]*_detuning_population[_i0] - + _d[1]*_akafield_detuning_population[_i0] - + _d[2]*_akifield_detuning_population[_i0] - + _d[3]*_akjfield_detuning_population[_i0] - + _d[4]*_akbfield_detuning_population[_i0] - + _d[5]*_akcfield_detuning_population[_i0] - + _d[6]*_akdfield_detuning_population[_i0] - + _d[7]*_akefield_detuning_population[_i0] - + _d[8]*_akffield_detuning_population[_i0] - + _d[9]*_akgfield_detuning_population[_i0] - + _d[10]*_akhfield_detuning_population[_i0]; - } - - - t += _a[14] * _step; - - _active_detuning_population = _akifield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - t += _a[15] * _step; - - _active_detuning_population = _akjfield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // Take full step - - // ai = a - memcpy(_initial_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - // a = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _detuning_population[_i0] += _c[0]*_akafield_detuning_population[_i0] + _c[7]*_akbfield_detuning_population[_i0] + _c[8]*_akcfield_detuning_population[_i0] - + _c[9]*_akdfield_detuning_population[_i0] + _c[10]*_akefield_detuning_population[_i0] + _c[11]*_akffield_detuning_population[_i0] - + _c[12]*_akgfield_detuning_population[_i0] + _c[13]*_akhfield_detuning_population[_i0] + _c[14]*_akifield_detuning_population[_i0] - + _c[15]*_akjfield_detuning_population[_i0]; - } - - - // a* = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akafield_detuning_population[_i0] = _initial_detuning_population[_i0] + _cs[0]*_akafield_detuning_population[_i0] + _cs[7]*_akbfield_detuning_population[_i0] - + _cs[8]*_akcfield_detuning_population[_i0] + _cs[9]*_akdfield_detuning_population[_i0] + _cs[10]*_akefield_detuning_population[_i0] - + _cs[11]*_akffield_detuning_population[_i0] + _cs[12]*_akgfield_detuning_population[_i0] + _cs[13]*_akhfield_detuning_population[_i0] - + _cs[14]*_akifield_detuning_population[_i0] + _cs[15]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _detuning_population; - - - - _error = 0.0; - - _segment1_detuning_population_error = _segment1_detuning_population_timestep_error(_akafield_detuning_population); - if (_segment1_detuning_population_error > _error) - _error = _segment1_detuning_population_error; - - _attempted_steps++; - - if (_error < _tolerance) { - _t_local += _step; - if (_step > _max_step) - _max_step = _step; - if (!_break_next && _step < _min_step) - _min_step = _step; - _discard = false; - } else { - t -= _step; - - if (_segment1_detuning_population_reset(_initial_detuning_population) == false) { - - _LOG(_WARNING_LOG_LEVEL, "WARNING: NaN present. Integration halted at t = %e.\n" - " Non-finite number in integration vector \"population\" in segment 1.\n", t); - if (_mg0_output_index_t < _mg0_output_lattice_t) - _mg0_sample(); - - goto _SEGMENT1_END; - } - - _segment1_ip_evolve(-1); - - _discard = true; - _break_next = false; - _unsuccessful_steps++; - } - - _old_step = _step; - - // Resize step - if (_error < 0.5*_tolerance || _error > _tolerance) { - const real _safetyFactor = 0.90; - real _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-0.7/9.0)) * pow(_last_norm_error, real(0.4/9.0)); - _scalingFactor = MAX(_scalingFactor, 1.0/5.0); - _scalingFactor = MIN(_scalingFactor, 7.0); - if (_error > _tolerance && _scalingFactor > 1.0) { - // If our step failed don't try and increase our step size. That would be silly. - _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-1.0/9.0)); - } - _old_step = _step; - _last_norm_error = pow(_safetyFactor/_scalingFactor*pow(_last_norm_error, real(0.4/9.0)), real(9.0/0.7)); - _step *= _scalingFactor; - } - - } while (_discard); - - if (_break_next) { - if (_next_sample_flag[0]) { - _mg0_sample(); - _next_sample_counter[0]++; - } - if (_next_sample_flag[1]) - _next_sample_flag[2] = true; - else { - _break_next = false; - _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - } - } - - if ( (_t_local + _step)*(1.0 + _EPSILON) > _t_break_next) { - _break_next = true; - _LOG(_SAMPLE_LOG_LEVEL, "Current timestep: %e\n", _old_step); - _step = _t_break_next - _t_local; - } - } while (!_next_sample_flag[2]); - - _SEGMENT1_END:; - xmds_free(_segment1_akafield_detuning_population); - xmds_free(_segment1_akbfield_detuning_population); - xmds_free(_segment1_akcfield_detuning_population); - xmds_free(_segment1_akdfield_detuning_population); - xmds_free(_segment1_akefield_detuning_population); - xmds_free(_segment1_akffield_detuning_population); - xmds_free(_segment1_akgfield_detuning_population); - xmds_free(_segment1_akhfield_detuning_population); - xmds_free(_segment1_akifield_detuning_population); - xmds_free(_segment1_akjfield_detuning_population); - xmds_free(_segment1_initial_detuning_population); - - _LOG(_SEGMENT_LOG_LEVEL, "Segment 1: minimum timestep: %e maximum timestep: %e\n", _min_step, _max_step); - _LOG(_SEGMENT_LOG_LEVEL, " Attempted %li steps, %.2f%% steps failed.\n", _attempted_steps, (100.0*_unsuccessful_steps)/_attempted_steps); -} - - -inline void _segment1_ip_evolve(int _exponent) -{ -} -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter) -{ - // The numbers of the moment groups that need to be sampled at the next sampling point. - // An entry of N+1 means "reached end of integration interval" - long _momentGroupNumbersNeedingSamplingNext[2]; - long _numberOfMomentGroupsToBeSampledNext = 1; - - long _previous_m = 1; - long _previous_M = 1; - - real _t_break_next = (real)final_time; - _momentGroupNumbersNeedingSamplingNext[0] = 1; - - // initialise all flags to false - for (long _i0 = 0; _i0 < 2; _i0++) - _next_sample_flag[_i0] = false; - - /* Check if moment group needs sampling at the same time as another already discovered sample (or the final time). - * If so, add this moment group to the to-be-sampled list. If moment group demands sampling earlier than all - * previously noted moment groups, erase all previous ones from list and set the sample time to this earlier one. - */ - if (_next_sample_counter[0] * _previous_M == _previous_m * 10) { - _momentGroupNumbersNeedingSamplingNext[_numberOfMomentGroupsToBeSampledNext] = 0; - _numberOfMomentGroupsToBeSampledNext++; - } else if (_next_sample_counter[0] * _previous_M < _previous_m * 10) { - _t_break_next = _next_sample_counter[0] * ((real)final_time) / ((real)10); - _numberOfMomentGroupsToBeSampledNext = 1; - _momentGroupNumbersNeedingSamplingNext[0] = 0; - _previous_M = 10; - _previous_m = _next_sample_counter[0]; - } - - // _momentGroupNumbersNeedingSamplingNext now contains the complete list of moment groups that need - // to be sampled at the next sampling point. Set their flags to true. - for (long _i0 = 0; _i0 < _numberOfMomentGroupsToBeSampledNext; _i0++) - _next_sample_flag[_momentGroupNumbersNeedingSamplingNext[_i0]] = true; - - return _t_break_next; -} - -real _segment1_detuning_population_timestep_error(complex* _checkfield) -{ - real _error = 1e-24; - real _temp_error = 0.0; - real _temp_mod = 0.0; - - - // Find the peak value for each component of the field - real _cutoff[_detuning_population_ncomponents]; - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) - _cutoff[_i0] = 0.0; - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - _temp_mod = mod2(_detuning_population[_detuning_population_index_pointer + _i1]); - if (_xmds_isnonfinite(_temp_mod)) - _cutoff[_i1] = INFINITY; - else if (_cutoff[_i1] < _temp_mod) - _cutoff[_i1] = _temp_mod; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) { - if (_xmds_isnonfinite(_cutoff[_i0])) - // Return an error two times the tolerance in this case because the timestep must be reduced. - return 2.0*1e-10; - _cutoff[_i0] *= 0.001; - _cutoff[_i0] *= 0.001; - } - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (mod2(_detuning_population[_detuning_population_index_pointer + _i1]) > _cutoff[_i1]) { - _temp_error = abs(_detuning_population[_detuning_population_index_pointer + _i1] - _checkfield[_detuning_population_index_pointer + _i1]) / (0.5*abs(_detuning_population[_detuning_population_index_pointer + _i1]) + 0.5*abs(_checkfield[_detuning_population_index_pointer + _i1])); - - if (_xmds_isnonfinite(_temp_error)) { - /* For _temp_error to be NaN, both the absolute value of the higher and lower order solutions - must BOTH be zero. This therefore implies that their difference is zero, and that there is no error. */ - _temp_error = 0.0; - } - - if (_error < _temp_error) // UNVECTORISABLE - _error = _temp_error; - } - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - return _error; -} - -bool _segment1_detuning_population_reset(complex* _reset_to_detuning_population) -{ - memcpy(_detuning_population, _reset_to_detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - /* return false if there's a NaN somewhere in the vector, otherwise return true */ - bool bNoNaNsPresent = true; - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (_xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Re()) - || _xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Im())) bNoNaNsPresent = false; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - return bNoNaNsPresent; -} - -// Delta A propagation operator for field detuning -void _segment1_detuning_operators_evaluate_operator0(real _step) -{ - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - complex dPca_dt; - complex dPab_dt; - complex dPaa_dt; - complex dPbb_dt; - complex dPcc_dt; - complex dPcb_dt; - - #define dt _step - - // ************* Propagation code *************** - #line 92 "./eit2.xmds" - - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - - #line 1678 "./eit2.cc" - // ********************************************** - - #undef dt - - - _active_detuning_population[_detuning_population_index_pointer + 4] = dPca_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 3] = dPab_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 0] = dPaa_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 1] = dPbb_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 2] = dPcc_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 5] = dPcb_dt * _step; - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// output function implementations -void _write_output() -{ - _LOG(_SIMULATION_LOG_LEVEL, "Generating output for ./eit2\n"); - - - char *_xsilFilename = (char*)malloc(256); - snprintf(_xsilFilename, 256, "%s.xsil", ("./eit2" + gsArgsAndValues).c_str()); - - FILE* _outfile = _open_xsil_file(_xsilFilename); - - if (_outfile) { - _write_xsil_header(_outfile); - char _dataFilename[200]; - snprintf(_dataFilename, 200, "%s.h5", ("./eit2" + gsArgsAndValues).c_str()); - - H5Fclose(H5Fcreate(_dataFilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); - } - _mg0_write_out(_outfile); - - _write_xsil_footer(_outfile); - _close_xsil_file(_outfile); - free(_xsilFilename); - _xsilFilename = NULL; - _outfile = NULL; - -} - - -FILE* _open_xsil_file(const char* _filename) -{ - - FILE* fp = fopen(_filename, "w"); - - if (fp == NULL) - // _LOG will cause the simulation to exit - _LOG(_ERROR_LOG_LEVEL, "Unable to open output file '%s'.\n" - "Exiting.\n", _filename); - - return fp; -} - -void _close_xsil_file(FILE*& fp) -{ - if (fp) - fclose(fp); - fp = NULL; - -} - -void _write_xsil_header(FILE* fp) -{ - if (!fp) - return; - fprintf(fp, "<?xml version=\"1.0\" ?><simulation xmds-version=\"2\">\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- While not strictly necessary, the following two tags are handy. -->\n"); - fprintf(fp, " <author>Hunter Rew</author>\n"); - fprintf(fp, " <description>\n"); - fprintf(fp, " Model of 3 state atom\n"); - fprintf(fp, " </description>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This element defines some constants. It can be used for other\n"); - fprintf(fp, " features as well, but we will go into that in more detail later.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <features>\n"); - fprintf(fp, " <precision> double </precision>\n"); - fprintf(fp, " <globals>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " real drive;\n"); - fprintf(fp, " real decay;\n"); - fprintf(fp, " real decay_bc;\n"); - fprintf(fp, " complex r_ca;\n"); - fprintf(fp, " real probe = .0001;\n"); - fprintf(fp, " real delta1;\n"); - fprintf(fp, " real dephase_bc;\n"); - fprintf(fp, " \n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </globals>\n"); - fprintf(fp, " <validation kind=\"run-time\"/>\n"); - fprintf(fp, " <arguments>\n"); - fprintf(fp, " <argument default_value=\"0.1\" name=\"drive_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10.0\" name=\"decay_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0.001\" name=\"decay_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"-1\" name=\"domain_min\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"1\" name=\"domain_max\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"delta1_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"dephase_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10000\" name=\"final_time\" type=\"real\"/>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " drive = drive_arg;\n"); - fprintf(fp, " decay = decay_arg;\n"); - fprintf(fp, " decay_bc = decay_bc_arg;\n"); - fprintf(fp, " delta1 = delta1_arg;\n"); - fprintf(fp, " dephase_bc = dephase_bc_arg;\n"); - fprintf(fp, " r_ca = decay - i*delta1;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </arguments>\n"); - fprintf(fp, " </features>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This part defines all of the dimensions used in the problem,\n"); - fprintf(fp, " in this case, only the dimension of 'time' is needed.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <geometry>\n"); - fprintf(fp, " <propagation_dimension> t </propagation_dimension>\n"); - fprintf(fp, " <transverse_dimensions>\n"); - fprintf(fp, " <dimension domain=\"(domain_min, domain_max)\" lattice=\"256\" name=\"detuning\"/>\n"); - fprintf(fp, " </transverse_dimensions>\n"); - fprintf(fp, " </geometry>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- A 'vector' describes the variables that we will be evolving. -->\n"); - fprintf(fp, " <vector name=\"Gamma\" type=\"complex\">\n"); - fprintf(fp, " <components> r_ab r_cb </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " r_ab = decay + i*(delta1 + detuning);\n"); - fprintf(fp, " r_cb = decay_bc + i*detuning + dephase_bc;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, " <vector dimensions=\"detuning\" name=\"population\" type=\"complex\">\n"); - fprintf(fp, " <components>\n"); - fprintf(fp, " Paa Pbb Pcc Pab Pca Pcb \n"); - fprintf(fp, " </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " Pcc = .5;\n"); - fprintf(fp, " Pbb = .5;\n"); - fprintf(fp, " Paa = Pab = Pca = Pcb = 0;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <sequence>\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " Here we define what differential equations need to be solved\n"); - fprintf(fp, " and what algorithm we want to use.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <integrate algorithm=\"ARK89\" interval=\"final_time\" tolerance=\"1e-10\">\n"); - fprintf(fp, " <samples>10</samples>\n"); - fprintf(fp, " <operators>\n"); - fprintf(fp, " <integration_vectors>population</integration_vectors>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " <dependencies>Gamma</dependencies>\n"); - fprintf(fp, " </operators>\n"); - fprintf(fp, " </integrate>\n"); - fprintf(fp, " </sequence>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- This part defines what data will be saved in the output file -->\n"); - fprintf(fp, " <output format=\"hdf5\">\n"); - fprintf(fp, " <sampling_group basis=\"detuning\" initial_sample=\"yes\">\n"); - fprintf(fp, " <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PaaR = Paa.Re();\n"); - fprintf(fp, " PaaI = Paa.Im();\n"); - fprintf(fp, " PbbR = Pbb.Re();\n"); - fprintf(fp, " PbbI = Pbb.Im();\n"); - fprintf(fp, " PccR = Pcc.Re();\n"); - fprintf(fp, " PccI = Pcc.Im();\n"); - fprintf(fp, " PabR = Pab.Re();\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " PcaR = Pca.Re();\n"); - fprintf(fp, " PcaI = Pca.Im();\n"); - fprintf(fp, " PcbR = Pcb.Re();\n"); - fprintf(fp, " PcbI = Pcb.Im();\n"); - fprintf(fp, " ]]> -->\n"); - fprintf(fp, " <moments>PabI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </sampling_group>\n"); - fprintf(fp, " </output>\n"); - - fprintf(fp, "\n<info>\n"); - fprintf(fp, "Script compiled with XMDS2 version 2.2.0 \"Out of cheese error\" (r2941)\n"); - fprintf(fp, "See http://www.xmds.org for more information.\n"); - fprintf(fp, "\nVariables that can be specified on the command line:\n"); - - fprintf(fp, " Command line argument drive_arg = %e\n", drive_arg); - - fprintf(fp, " Command line argument decay_arg = %e\n", decay_arg); - - fprintf(fp, " Command line argument decay_bc_arg = %e\n", decay_bc_arg); - - fprintf(fp, " Command line argument domain_min = %e\n", domain_min); - - fprintf(fp, " Command line argument domain_max = %e\n", domain_max); - - fprintf(fp, " Command line argument delta1_arg = %e\n", delta1_arg); - - fprintf(fp, " Command line argument dephase_bc_arg = %e\n", dephase_bc_arg); - - fprintf(fp, " Command line argument final_time = %e\n", final_time); - fprintf(fp, "</info>\n"); - -} - -// In addition to writing the footer (if 'fp' is not NULL) -// this function closes the fp file pointer. -void _write_xsil_footer(FILE* fp) -{ - if (fp) { - fprintf(fp, "</simulation>\n"); - } -} - -// ******************************************************** -// moment group 0 function implementations -void _mg0_sample() -{ - - long _mg0_output_raw_index_pointer = 0; - #define PabI _active_mg0_output_raw[_mg0_output_raw_index_pointer + 0] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // Set index pointers explicitly for (some) vectors - _mg0_output_raw_index_pointer = ( 0 - + _mg0_output_index_t * _lattice_detuning - + _index_detuning * 1 ) * _mg0_output_raw_ncomponents; - #define _SAMPLE_COMPLEX(variable) \ - variable ## R = variable.Re(); variable ## I = variable.Im(); - - // *************** Sampling code **************** - #line 131 "./eit2.xmds" - - PabI = Pab.Im(); - - #line 1960 "./eit2.cc" - // ********************************************** - - #undef _SAMPLE_COMPLEX - // Increment index pointers for vectors in field mg0_sampling (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef PabI - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - - _mg0_output_t[0 + _mg0_output_index_t++] = t; - - _LOG(_SAMPLE_LOG_LEVEL, "Sampled field (for moment group #1) at t = %e\n", t); - -} - - -void _mg0_process() -{ - // No post processing needs to be done -} - - -void _mg0_write_out(FILE* _outfile) -{ - - if (_outfile) { - fprintf(_outfile, "\n"); - fprintf(_outfile, "<XSIL Name=\"moment_group_1\">\n"); - fprintf(_outfile, " <Param Name=\"n_independent\">2</Param>\n"); - fprintf(_outfile, " <Array Name=\"variables\" Type=\"Text\">\n"); - fprintf(_outfile, " <Dim>3</Dim>\n"); - fprintf(_outfile, " <Stream><Metalink Format=\"Text\" Delimiter=\" \\n\"/>\n"); - fprintf(_outfile, "t detuning PabI \n"); - fprintf(_outfile, " </Stream>\n"); - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, " <Array Name=\"data\" Type=\"double\">\n"); - fprintf(_outfile, " <Dim>%i</Dim>\n", _mg0_output_lattice_t); - fprintf(_outfile, " <Dim>%i</Dim>\n", _lattice_detuning); - fprintf(_outfile, " <Dim>3</Dim>\n"); - } - - - char _h5Filename[200]; - snprintf(_h5Filename, 200, "%s.h5", ("./eit2" + gsArgsAndValues).c_str()); - - /* Open the file */ - hid_t hdf5_file = H5Fopen(_h5Filename, H5F_ACC_RDWR, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_WARNING_LOG_LEVEL, "Failed to open HDF5 file '%s', will try to create it.", _h5Filename); - hdf5_file = H5Fcreate(_h5Filename, H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_ERROR_LOG_LEVEL, "Failed to create HDF5 file '%s'. Bailing.", _h5Filename); - } - } - - /* Create the group for this data */ - hid_t group; - if (!H5Lexists(hdf5_file, "/1", H5P_DEFAULT)) - group = H5Gcreate(hdf5_file, "/1", H5P_DEFAULT); - else - group = H5Gopen(hdf5_file, "/1"); - - if (_outfile) { - fprintf(_outfile, " <Stream><Metalink Format=\"HDF5\" Type=\"Remote\" Group=\"/1\"/>\n"); - fprintf(_outfile, "%s.h5\n", ("./eit2" + gsArgsAndValues).c_str()); - fprintf(_outfile, " </Stream>\n"); - } - - /* Create the coordinate data sets */ - hsize_t coordinate_length; - hid_t coordinate_dataspace; - coordinate_length = _mg0_output_lattice_t; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_t; - if (!H5Lexists(hdf5_file, "/1/t", H5P_DEFAULT)) - dataset_t = H5Dcreate(hdf5_file, "/1/t", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_t = H5Dopen(hdf5_file, "/1/t"); - H5Dwrite(dataset_t, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _mg0_output_t); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_t, "t"); - #endif - - H5Sclose(coordinate_dataspace); - coordinate_length = _lattice_detuning; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_detuning; - if (!H5Lexists(hdf5_file, "/1/detuning", H5P_DEFAULT)) - dataset_detuning = H5Dcreate(hdf5_file, "/1/detuning", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_detuning = H5Dopen(hdf5_file, "/1/detuning"); - H5Dwrite(dataset_detuning, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _detuning); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_detuning, "detuning"); - #endif - - H5Sclose(coordinate_dataspace); - - hsize_t file_dims[] = {_mg0_output_lattice_t, _lattice_detuning}; - hid_t file_dataspace = H5Screate_simple(2, file_dims, NULL); - - hid_t dataset_PabI; - if (!H5Lexists(hdf5_file, "/1/PabI", H5P_DEFAULT)) - dataset_PabI = H5Dcreate(hdf5_file, "/1/PabI", H5T_NATIVE_REAL, file_dataspace, H5P_DEFAULT); - else - dataset_PabI = H5Dopen(hdf5_file, "/1/PabI"); - #if defined(HAVE_HDF5_HL) - H5DSattach_scale(dataset_PabI, dataset_t, 0); - H5DSattach_scale(dataset_PabI, dataset_detuning, 1); - #endif - H5Dclose(dataset_t); - H5Dclose(dataset_detuning); - - - if ((_mg0_output_lattice_t * _lattice_detuning)) { - /* Create the data space */ - hsize_t file_start[2] = {0, 0}; - hsize_t mem_dims[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - hsize_t mem_start[3] = {0, 0, 0}; - hsize_t mem_stride[3] = {1, 1, 1}; - hsize_t mem_count[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - - - hid_t mem_dataspace; - mem_dims[2] = 1; - mem_dataspace = H5Screate_simple(3, mem_dims, NULL); - mem_stride[2] = 1; - - // Select hyperslabs of memory and file data spaces for data transfer operation - mem_start[2] = 0; - H5Sselect_hyperslab(mem_dataspace, H5S_SELECT_SET, mem_start, mem_stride, mem_count, NULL); - H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, file_start, mem_stride, mem_count, NULL); - - if (dataset_PabI) - H5Dwrite(dataset_PabI, H5T_NATIVE_REAL, mem_dataspace, file_dataspace, H5P_DEFAULT, _active_mg0_output_raw); - - H5Sclose(mem_dataspace); - } - - - H5Dclose(dataset_PabI); - - H5Sclose(file_dataspace); - H5Gclose(group); - H5Fclose(hdf5_file); - - - if (_outfile) { - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, "</XSIL>\n"); - } -} - -// ******************************************************** -// field mg0_output function implementations -// initialisation for vector raw -void _mg0_output_raw_initialise() -{ - - bzero(_active_mg0_output_raw, sizeof(real) * _mg0_output_raw_alloc_size); -} - diff --git a/source/eit2.h5 b/source/eit2.h5 Binary files differdeleted file mode 100644 index 3c852dd..0000000 --- a/source/eit2.h5 +++ /dev/null diff --git a/source/eit2.py b/source/eit2.py deleted file mode 100644 index f4f7022..0000000 --- a/source/eit2.py +++ /dev/null @@ -1,17 +0,0 @@ -#!/usr/bin/env python -from xpdeint.XSILFile import XSILFile - -xsilFile = XSILFile("./eit2.xsil") - -def firstElementOrNone(enumerable): - for element in enumerable: - return element - return None - -t_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "t") -detuning_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "detuning") -PabI_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].dependentVariables if _["name"] == "PabI") - -# Write your plotting commands here. -# You may want to import pylab (from pylab import *) or matplotlib (from matplotlib import *) - diff --git a/source/eit2.pyc b/source/eit2.pyc Binary files differdeleted file mode 100644 index 81a774e..0000000 --- a/source/eit2.pyc +++ /dev/null diff --git a/source/eit2.xmds b/source/eit2.xmds deleted file mode 100644 index 50419f5..0000000 --- a/source/eit2.xmds +++ /dev/null @@ -1,136 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument name="drive_arg" type="real" default_value="0.1" /> - <argument name="decay_arg" type="real" default_value="10.0" /> - <argument name="decay_bc_arg" type="real" default_value="0.001" /> - <argument name="domain_min" type="real" default_value="-1" /> - <argument name="domain_max" type="real" default_value="1" /> - <argument name="delta1_arg" type="real" default_value="0" /> - <argument name="dephase_bc_arg" type="real" default_value="0" /> - <argument name="final_time" type="real" default_value="10000" /> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension name="detuning" lattice="256" domain="(domain_min, domain_max)" /> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector name="population" type="complex" dimensions="detuning"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5" > - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> -</simulation> diff --git a/source/eit2.xsil b/source/eit2.xsil deleted file mode 100644 index 28c0ff0..0000000 --- a/source/eit2.xsil +++ /dev/null @@ -1,168 +0,0 @@ -<?xml version="1.0" ?><simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument default_value="0.1" name="drive_arg" type="real"/> - <argument default_value="10.0" name="decay_arg" type="real"/> - <argument default_value="0.001" name="decay_bc_arg" type="real"/> - <argument default_value="-1" name="domain_min" type="real"/> - <argument default_value="1" name="domain_max" type="real"/> - <argument default_value="0" name="delta1_arg" type="real"/> - <argument default_value="0" name="dephase_bc_arg" type="real"/> - <argument default_value="10000" name="final_time" type="real"/> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension domain="(domain_min, domain_max)" lattice="256" name="detuning"/> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector dimensions="detuning" name="population" type="complex"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5"> - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> - -<info> -Script compiled with XMDS2 version 2.2.0 "Out of cheese error" (r2941) -See http://www.xmds.org for more information. - -Variables that can be specified on the command line: - Command line argument drive_arg = 9.545485e-01 - Command line argument decay_arg = 6.000000e+00 - Command line argument decay_bc_arg = 1.000000e-05 - Command line argument domain_min = -1.518705e+00 - Command line argument domain_max = 1.518705e+00 - Command line argument delta1_arg = -1.500000e+01 - Command line argument dephase_bc_arg = 0.000000e+00 - Command line argument final_time = 6.584559e+02 -</info> - -<XSIL Name="moment_group_1"> - <Param Name="n_independent">2</Param> - <Array Name="variables" Type="Text"> - <Dim>3</Dim> - <Stream><Metalink Format="Text" Delimiter=" \n"/> -t detuning PabI - </Stream> - </Array> - <Array Name="data" Type="double"> - <Dim>11</Dim> - <Dim>256</Dim> - <Dim>3</Dim> - <Stream><Metalink Format="HDF5" Type="Remote" Group="/1"/> -./eit2.h5 - </Stream> - </Array> -</XSIL> -</simulation> diff --git a/source/eit3 b/source/eit3 Binary files differdeleted file mode 100755 index 64fa546..0000000 --- a/source/eit3 +++ /dev/null diff --git a/source/eit3.cc b/source/eit3.cc deleted file mode 100644 index dd04f95..0000000 --- a/source/eit3.cc +++ /dev/null @@ -1,2129 +0,0 @@ -// ******************************************************** -// simulation logging - -#define _SAMPLE_LOG_LEVEL (1 << 0) -#define _SEGMENT_LOG_LEVEL (1 << 1) -#define _PATH_LOG_LEVEL (1 << 2) -#define _SIMULATION_LOG_LEVEL (1 << 3) -#define _WARNING_LOG_LEVEL (1 << 4) -#define _ERROR_LOG_LEVEL (1 << 5) -#define _NO_ERROR_TERMINATE_LOG_LEVEL (1 << 6) -#define _ALL_LOG_LEVELS _SAMPLE_LOG_LEVEL|_SEGMENT_LOG_LEVEL|_PATH_LOG_LEVEL|_SIMULATION_LOG_LEVEL|_WARNING_LOG_LEVEL|_ERROR_LOG_LEVEL|_NO_ERROR_TERMINATE_LOG_LEVEL -#define _LOG_LEVELS_BEING_LOGGED (_ALL_LOG_LEVELS) - -#define real Re -#define imag Im - -#include <complex> - -#undef real -#undef imag - - -#include <stdio.h> - -#define _LOG(logLevel, ...) \ - do { \ - if (logLevel & _LOG_LEVELS_BEING_LOGGED) { \ - if (logLevel & (_ERROR_LOG_LEVEL | _WARNING_LOG_LEVEL)) \ - printf("%s:%i: ", __FILE__, __LINE__); \ - printf(__VA_ARGS__); \ - fflush(stdout); \ - if (logLevel & (_ERROR_LOG_LEVEL | _NO_ERROR_TERMINATE_LOG_LEVEL)) \ - exit(logLevel == _ERROR_LOG_LEVEL); \ - } \ - } while (0) - -// ******************************************************** -// simulation includes - -#include <xpdeint_platform.h> -#include <cmath> -#include <string> -#include <cstring> -#include <fstream> -#include <sstream> -#include <cstdlib> - -#if CFG_OSAPI == CFG_OSAPI_POSIX // These are POSIX headers (i.e. non-windows) - #include <sys/time.h> -#endif // POSIX - -#ifdef __APPLE__ - #include <Availability.h> - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - #define OS_OBJECT_USE_OBJC 0 // Don't make dispatch and xpc objects Objective-C objects. - #include <IOKit/pwr_mgt/IOPMLib.h> // To disable user idle sleep on Mountain Lion - #endif -#endif - -#include <time.h> -#include <list> -#include <vector> -#include <algorithm> - -#include <utility> -#include <map> - -#include <getopt.h> - -#if (CFG_COMPILER == CFG_COMPILER_MSVC) - #define FFTW_DLL -#endif - -#include <fftw3.h> -#include <sys/stat.h> -#include <sys/types.h> - -#define _xmds_malloc fftw_malloc -#define xmds_free fftw_free - -#define H5_USE_16_API -#include <hdf5.h> - -#if !defined(HAVE_H5LEXISTS) -htri_t H5Lexists(hid_t loc_id, const char *name, hid_t lapl_id) -{ - H5E_auto_t error_func; - void* error_client_data; - // Squelch errors generated by H5Gget_objinfo. It will report errors when it can't find an object - // but that's the purpose of calling it. - H5Eget_auto(&error_func, &error_client_data); - H5Eset_auto(NULL, NULL); - herr_t err = H5Gget_objinfo(loc_id, name, false, NULL); - H5Eset_auto(error_func, error_client_data); - if (err >= 0) - return true; - else - return false; -} -#endif - -#define H5T_NATIVE_REAL H5T_NATIVE_DOUBLE -#if defined(HAVE_HDF5_HL) - #include <hdf5_hl.h> -#endif - - -typedef long integer; -typedef double real; -typedef std::complex<real> XMDSComplexType; - -#include <xpdeint.h> - -#define complex XMDSComplexType - -const complex i(0.0, 1.0); - -using namespace std; - -#if CFG_COMPILER == CFG_COMPILER_ICC - // - // Disable ICC's warning: label was declared but never referenced - // - #pragma warning ( disable : 177 ) -#endif - -inline void *xmds_malloc(size_t size); - -// ******************************************************** -// DEFINES -// ******************************************************** - -// ******************************************************** -// Simulation defines -#define _EPSILON 1e-6 -#ifndef INFINITY -#define INFINITY HUGE_VAL -#endif - -#ifndef MAX -#define MAX(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a > _b ? _a : _b; }) -#endif - -#ifndef MIN -#define MIN(a, b) \ - ({ typeof(a) _a = (a); \ - typeof(b) _b = (b); \ - _a < _b ? _a : _b; }) -#endif - - -// ******************************************************** -// Geometry defines -#define _lattice_detuning ((int)256) -#define _min_detuning ((real)domain_min) -#define _max_detuning ((real)domain_max) -#define _ddetuning ((real)((_max_detuning - _min_detuning)/_lattice_detuning)) - -#define _lattice_kdetuning ((int)256) -#define _dkdetuning (2.0*M_PI/(_max_detuning - _min_detuning)) -#define _min_kdetuning (-(_lattice_kdetuning/2) * _dkdetuning) -#define _max_kdetuning ((_lattice_kdetuning - 1)/2 * _dkdetuning) - -// ******************************************************** -// field detuning defines -#define _detuning_ndims 1 - - -// vector Gamma defines -#define _detuning_Gamma_ncomponents 2 - -// vector population defines -#define _detuning_population_ncomponents 6 - -// ******************************************************** -// field mg0_sampling defines -#define _mg0_sampling_ndims 1 - - -// ******************************************************** -// field mg0_output defines -#define _mg0_output_ndims 2 - - -#define _mg0_output_lattice_t ((int)11) -#define _mg0_output_min_t (_mg0_output_t[0]) -#define _mg0_output_max_t (_mg0_output_t[_mg0_output_lattice_t-1]) -#define _mg0_output_dt (_mg0_output_t[_index_t+1]-_mg0_output_t[_index_t]) - -// vector raw defines -#define _mg0_output_raw_ncomponents 1 - - -// ******************************************************** -// GLOBALS -// ******************************************************** - - -// ******************************************************** -// Simulation globals - -string gsArgsAndValues = ""; - -real t; - -// ******************************************************** -// Transform Multiplexer globals -typedef pair<ptrdiff_t, ptrdiff_t> _basis_pair; -typedef void (*transform_function)(bool, real, real* const __restrict__, real* const __restrict__, ptrdiff_t, ptrdiff_t); - -// Less than operator needed by the C++ map class -struct _basis_pair_less_than -{ - bool operator()(const _basis_pair& _x, const _basis_pair& _y) const { - return (_x.first < _y.first) || ((_x.first == _y.first) && (_x.second < _y.second)); - } -}; - -struct _transform_step -{ - transform_function _func; - bool _forward; - bool _out_of_place; - ptrdiff_t _prefix_lattice; - ptrdiff_t _postfix_lattice; -}; - -// Structure to hold the basis change information -struct _basis_transform_t -{ - vector<_transform_step> _transform_steps; - real _multiplier; - - _basis_transform_t(real _multiplier_in = 1.0) : _multiplier(_multiplier_in) {} - - _basis_transform_t(const _basis_transform_t& _b) : _transform_steps(_b._transform_steps), _multiplier(_b._multiplier) {} - - void append(transform_function _func, bool _forward, bool _out_of_place, ptrdiff_t _prefix_lattice, ptrdiff_t _postfix_lattice) - { - _transform_steps.push_back((_transform_step){_func, _forward, _out_of_place, _prefix_lattice, _postfix_lattice}); - } -}; - -// Map type for holding (old_basis, new_basis) -> _basis_transform_t mappings -typedef map<_basis_pair, _basis_transform_t, _basis_pair_less_than> _basis_map; - - -real *_auxiliary_array = NULL; - -const char *_basis_identifiers[] = { -}; - -// ******************************************************** -// 'Globals' element globals - -#line 17 "./eit3.xmds" - -real drive; -real decay; -real decay_bc; -complex r_ca; -real probe = .0001; -real delta1; -real dephase_bc; - - -#line 271 "./eit3.cc" - -// ******************************************************** -// Command line argument processing globals -real drive_arg = 0.1; -real decay_arg = 10.0; -real decay_bc_arg = 0.001; -real domain_min = -1; -real domain_max = 1; -real delta1_arg = 0; -real dephase_bc_arg = 0; -real final_time = 10000; - -// ******************************************************** -// Geometry globals -real* _detuning = NULL; - -real* _kdetuning = NULL; - -// ******************************************************** -// FFTW3 globals -const real _inverse_sqrt_2pi = 1.0 / sqrt(2.0 * M_PI); -string _fftwWisdomPath; - -// ******************************************************** -// field detuning globals -// vector Gamma globals -size_t _detuning_Gamma_alloc_size = 0; -complex* _detuning_Gamma = NULL; -complex* _active_detuning_Gamma = NULL; - -// vector population globals -size_t _detuning_population_alloc_size = 0; -complex* _detuning_population = NULL; -complex* _active_detuning_population = NULL; - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) globals -complex* _segment1_akafield_detuning_population; -complex* _segment1_akbfield_detuning_population; -complex* _segment1_akcfield_detuning_population; -complex* _segment1_akdfield_detuning_population; -complex* _segment1_akefield_detuning_population; -complex* _segment1_akffield_detuning_population; -complex* _segment1_akgfield_detuning_population; -complex* _segment1_akhfield_detuning_population; -complex* _segment1_akifield_detuning_population; -complex* _segment1_akjfield_detuning_population; -complex* _segment1_initial_detuning_population; - -// ******************************************************** -// field mg0_output globals -real* _mg0_output_t = NULL; -unsigned long _mg0_output_index_t = 0; - -// vector raw globals -size_t _mg0_output_raw_alloc_size = 0; -real* _mg0_output_raw = NULL; -real* _active_mg0_output_raw = NULL; - - -// ******************************************************** -// FUNCTION PROTOTYPES -// ******************************************************** - -// ******************************************************** -// Command line argument processing function prototypes -void _print_usage(); - -// ******************************************************** -// field detuning function prototypes -void _detuning_Gamma_initialise(); - -void _detuning_population_initialise(); - -// ******************************************************** -// segment 0 (Top level sequence) function prototypes -void _segment0(); - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function prototypes -inline void _segment1_calculate_delta_a(real _step); -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent); -void _segment1(); -inline void _segment1_ip_evolve(int _exponent); -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter); -real _segment1_detuning_population_timestep_error(complex* _checkfield); -bool _segment1_detuning_population_reset(complex* _reset_to); - -void _segment1_detuning_operators_evaluate_operator0(real _step); - -// ******************************************************** -// output function prototypes -void _write_output(); - -FILE* _open_xsil_file(const char* _filename); -void _close_xsil_file(FILE*& fp); -void _write_xsil_header(FILE* fp); -void _write_xsil_footer(FILE* fp); - -// ******************************************************** -// moment group 0 function prototypes -void _mg0_sample(); -void _mg0_process(); -void _mg0_write_out(FILE* _outfile); - -// ******************************************************** -// field mg0_output function prototypes -void _mg0_output_raw_initialise(); - -// ******************************************************** -// MAIN ROUTINE -// ******************************************************** -int main(int argc, char **argv) -{ - #ifdef __APPLE__ - #if __MAC_OS_X_VERSION_MIN_REQUIRED >= 1080 - { - IOPMAssertionID _powerAssertionID = 0; - IOReturn __io_result = IOPMAssertionCreateWithDescription( - kIOPMAssertionTypePreventUserIdleSystemSleep, - CFSTR("XMDS simulation './eit3' preventing user idle sleep"), // Assertion name - NULL, // Assertion details - NULL, // Human-readable reason - NULL, // Localization bundle path - (CFTimeInterval)0, // never timeout - kIOPMAssertionTimeoutActionRelease, - &_powerAssertionID - ); - if (__io_result != kIOReturnSuccess) { - _LOG(_WARNING_LOG_LEVEL, "Failed to disable user idle sleep\n"); - } - // Note, this power assertion is automatically released when the process quits. - } - #endif - #endif - - - // *********** Parse the command line for arguments, and set ********* - // *********** the appropriate global variables ********* - - int resp; - std::map<string, string> mInputArgsAndValues; - - while (1) { - static struct option long_options[] = - { - {"help", no_argument, 0, 'h'}, - {"drive_arg", required_argument, 0, 'd'}, - {"decay_arg", required_argument, 0, 'e'}, - {"decay_bc_arg", required_argument, 0, 'c'}, - {"domain_min", required_argument, 0, 'o'}, - {"domain_max", required_argument, 0, 'm'}, - {"delta1_arg", required_argument, 0, 'l'}, - {"dephase_bc_arg", required_argument, 0, 'p'}, - {"final_time", required_argument, 0, 'f'}, - {NULL, 0, 0, 0} - }; - - int option_index = 0; - - resp = getopt_long(argc, argv, "hd:e:c:o:m:l:p:f:", long_options, &option_index); - - if (resp == -1) - break; - - switch (resp) { - case '?': - // An unknown option was passed. Show allowed options and exit. - _print_usage(); // This causes the simulation to exit - - case 'h': - _print_usage(); // This causes the simulation to exit - - case 'd': - drive_arg = strtod(optarg, NULL); - break; - - case 'e': - decay_arg = strtod(optarg, NULL); - break; - - case 'c': - decay_bc_arg = strtod(optarg, NULL); - break; - - case 'o': - domain_min = strtod(optarg, NULL); - break; - - case 'm': - domain_max = strtod(optarg, NULL); - break; - - case 'l': - delta1_arg = strtod(optarg, NULL); - break; - - case 'p': - dephase_bc_arg = strtod(optarg, NULL); - break; - - case 'f': - final_time = strtod(optarg, NULL); - break; - - default: - _LOG(_ERROR_LOG_LEVEL, "Internal error in processing arguments.\n"); - } - } - - - if (optind < argc) - _print_usage(); // This causes the simulation to exit. - - // ******** Argument post-processing code ******* - #line 38 "./eit3.xmds" - - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - - #line 496 "./eit3.cc" - // ********************************************** - - - - _mg0_output_raw_alloc_size = MAX(_mg0_output_raw_alloc_size, (_mg0_output_lattice_t * _lattice_detuning) * _mg0_output_raw_ncomponents); - _detuning_Gamma_alloc_size = MAX(_detuning_Gamma_alloc_size, (_lattice_detuning) * _detuning_Gamma_ncomponents); - _detuning_population_alloc_size = MAX(_detuning_population_alloc_size, (_lattice_detuning) * _detuning_population_ncomponents); - _detuning = (real*) xmds_malloc(sizeof(real) * (_lattice_detuning+1)); - - _kdetuning = (real*) xmds_malloc(sizeof(real) * (_lattice_kdetuning+1)); - - _detuning_Gamma = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_Gamma_alloc_size,1)); - _active_detuning_Gamma = _detuning_Gamma; - - - _detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _active_detuning_population = _detuning_population; - _mg0_output_t = (real*) xmds_malloc(sizeof(real) * (_mg0_output_lattice_t+1)); - - - _mg0_output_raw = (real*) xmds_malloc(sizeof(real) * MAX(_mg0_output_raw_alloc_size,1)); - _active_mg0_output_raw = _mg0_output_raw; - - - // Run-time validation checks - - if (domain_min >= domain_max) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The end point of the dimension 'domain_max' must be " - "greater than the start point.\n" - "Start = %e, End = %e\n", (real)domain_min, (real)domain_max); - - if (final_time <= 0.0) - _LOG(_ERROR_LOG_LEVEL, "ERROR: The interval for segment 1 is not positive!\n" - "Interval = %e\n", final_time); - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) - _detuning[_index_detuning] = _min_detuning + _index_detuning*_ddetuning; - for (long _index_kdetuning = 0; _index_kdetuning < (_lattice_kdetuning+1)/2; _index_kdetuning++) - _kdetuning[_index_kdetuning] = _index_kdetuning*_dkdetuning; - for (long _index_kdetuning = (_lattice_kdetuning+1)/2; _index_kdetuning < _lattice_kdetuning; _index_kdetuning++) - _kdetuning[_index_kdetuning] = -(_lattice_kdetuning - _index_kdetuning) * _dkdetuning; - _active_mg0_output_raw = _mg0_output_raw; - _mg0_output_raw_initialise(); - // load wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX // Don't load wisdom on windows - { - char _hostName[256]; - gethostname(_hostName, 256); - _hostName[255] = '\0'; // just in case - - string _pathToWisdom = getenv("HOME"); - _pathToWisdom += "/.xmds/wisdom/"; - - string _wisdomFileName = _hostName; - _wisdomFileName += ".wisdom"; - _wisdomFileName += ".fftw3"; - - FILE *_fp = NULL; - - _fp = fopen(_pathToWisdom.c_str(), "r"); - if (_fp) { - fclose(_fp); - } else { - int _result = mkdir((string(getenv("HOME")) + "/.xmds").c_str(), S_IRWXU); - if (mkdir(_pathToWisdom.c_str(), S_IRWXU)) { - // We failed to create the ~/.xmds/wisdom directory - _LOG(_WARNING_LOG_LEVEL, "Warning: Cannot find enlightenment, the path to wisdom ~/.xmds/wisdom doesn't seem to exist and we couldn't create it.\n" - " I'll use the current path instead.\n"); - _pathToWisdom = ""; // present directory - } - - } - - _fftwWisdomPath = _pathToWisdom + _wisdomFileName; - - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "r")) != NULL) { - _LOG(_SIMULATION_LOG_LEVEL, "Found enlightenment... (Importing wisdom)\n"); - fftw_import_wisdom_from_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - _basis_transform_t *_basis_transform = NULL; - ptrdiff_t _auxiliary_array_size = 0; - ptrdiff_t _max_vector_size = 0; - real* _max_vector_array = NULL; - - if (_auxiliary_array_size) { - _auxiliary_array = (real*) xmds_malloc(sizeof(real) * _auxiliary_array_size); - } - - bool _allocated_temporary_array = false; - if (!_max_vector_array && _max_vector_size > 0) { - _max_vector_array = (real*) xmds_malloc(sizeof(real) * _max_vector_size); - _allocated_temporary_array = true; - } - - // Make all geometry-dependent transformations prepare plans, etc. - - if (_allocated_temporary_array) { - xmds_free(_max_vector_array); - } - - /* Code that actually does stuff goes here */ - _segment0(); - - - _write_output(); - if (_auxiliary_array) { - xmds_free(_auxiliary_array); - } - - // Save wisdom - #if CFG_OSAPI == CFG_OSAPI_POSIX - { - FILE *_wisdomFile = NULL; - if ( (_wisdomFile = fopen(_fftwWisdomPath.c_str(), "w")) != NULL) { - fftw_export_wisdom_to_file(_wisdomFile); - fclose(_wisdomFile); - } - } - #endif // POSIX - - fftw_cleanup(); - - xmds_free(_detuning_Gamma); - _active_detuning_Gamma = _detuning_Gamma = NULL; - - - xmds_free(_detuning_population); - _active_detuning_population = _detuning_population = NULL; - - xmds_free(_mg0_output_raw); - _active_mg0_output_raw = _mg0_output_raw = NULL; - - - return 0; -} - -// ******************************************************** -// FUNCTION IMPLEMENTATIONS -// ******************************************************** - -inline void *xmds_malloc(size_t size) -{ - void *retPointer = _xmds_malloc(size); - if ( !retPointer ) - _LOG(_ERROR_LOG_LEVEL, "ERROR: Couldn't allocate %zu bytes of memory!", size); - return retPointer; -} - - -// ******************************************************** -// Command line argument processing function implementations -void _print_usage() -{ - // This function does not return. - _LOG(_NO_ERROR_TERMINATE_LOG_LEVEL, "\n\nUsage: ./eit3 --drive_arg <real> --decay_arg <real> --decay_bc_arg <real> --domain_min <real> --domain_max <real> --delta1_arg <real> --dephase_bc_arg <real> --final_time <real>\n\n" - "Details:\n" - "Option\t\tType\t\tDefault value\n" - "-d, --drive_arg\treal \t\t0.1\n" - "-e, --decay_arg\treal \t\t10.0\n" - "-c, --decay_bc_arg\treal \t\t0.001\n" - "-o, --domain_min\treal \t\t-1\n" - "-m, --domain_max\treal \t\t1\n" - "-l, --delta1_arg\treal \t\t0\n" - "-p, --dephase_bc_arg\treal \t\t0\n" - "-f, --final_time\treal \t\t10000\n" - ); - // _LOG terminates the simulation. -} - -// ******************************************************** -// Default Simulation Driver function implementations -void _segment0() -{ - t = 0.0; - - _mg0_output_raw_initialise(); - _active_detuning_Gamma = _detuning_Gamma; - _detuning_Gamma_initialise(); - _active_detuning_population = _detuning_population; - _detuning_population_initialise(); - _mg0_output_index_t = 0; - _mg0_sample(); - _segment1(); - - _mg0_process(); -} - - -// ******************************************************** -// field detuning function implementations -// initialisation for vector Gamma -void _detuning_Gamma_initialise() -{ - - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 64 "./eit3.xmds" - - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - - #line 714 "./eit3.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb -} - -// initialisation for vector population -void _detuning_population_initialise() -{ - - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // The purpose of the following define is to give a (somewhat helpful) compile-time error - // if the user has attempted to use the propagation dimension variable in the initialisation - // block of a <vector> element. If they're trying to do this, what they really want is a - // <computed_vector> instead. - #define t Dont_use_propagation_dimension_t_in_vector_element_CDATA_block___Use_a_computed_vector_instead - - // ********** Initialisation code *************** - #line 75 "./eit3.xmds" - - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - - #line 756 "./eit3.cc" - // ********************************************** - #undef t - - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// segment 1 (RK89 adaptive-step integrator) function implementations -inline void _segment1_calculate_delta_a(real _step) -{ - - - // Delta A propagation operator for field detuning - _segment1_detuning_operators_evaluate_operator0(_step); - -} - - -inline void _segment1_calculate_nonconstant_ip_fields(real _step, int _exponent) -{ -} - - -void _segment1() -{ - real _step = final_time/(real)10; - real _old_step = _step; - real _min_step = _step; - real _max_step = _step; - long _attempted_steps = 0; - long _unsuccessful_steps = 0; - - real _tolerance = 1e-10; - - real _error, _last_norm_error = 1.0; - real _segment1_detuning_population_error; - - bool _discard = false; - bool _break_next = false; - - bool _next_sample_flag[3]; - for (long _i0 = 0; _i0 < 3; _i0++) - _next_sample_flag[_i0] = false; - - long _next_sample_counter[1]; - for (long _i0 = 0; _i0 < 1; _i0++) - _next_sample_counter[_i0] = 1; - - real _t_local = 0.0; - - real _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - - if ( (_t_local + _step)*(1.0 + _EPSILON) >= _t_break_next) { - _break_next = true; - _step = _t_break_next - _t_local; - } - - _segment1_akafield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akbfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akcfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akdfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akefield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akffield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akgfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akhfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akifield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_akjfield_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - _segment1_initial_detuning_population = (complex*) xmds_malloc(sizeof(complex) * MAX(_detuning_population_alloc_size,1)); - complex* _akafield_detuning_population = _segment1_akafield_detuning_population; - complex* _akbfield_detuning_population = _segment1_akbfield_detuning_population; - complex* _akcfield_detuning_population = _segment1_akcfield_detuning_population; - complex* _akdfield_detuning_population = _segment1_akdfield_detuning_population; - complex* _akefield_detuning_population = _segment1_akefield_detuning_population; - complex* _akffield_detuning_population = _segment1_akffield_detuning_population; - complex* _akgfield_detuning_population = _segment1_akgfield_detuning_population; - complex* _akhfield_detuning_population = _segment1_akhfield_detuning_population; - complex* _akifield_detuning_population = _segment1_akifield_detuning_population; - complex* _akjfield_detuning_population = _segment1_akjfield_detuning_population; - complex* _initial_detuning_population = _segment1_initial_detuning_population; - - - // Runge Kutta method constants - real _a_raw[16]; - real _a[16]; - real _b[16][16]; - real _c[16]; - real _cs[16]; - real _d[16]; - - for (unsigned long _i0 = 0; _i0 < 16; _i0++) { - _a_raw[_i0] = _c[_i0] = _d[_i0] = 0.0; - for (unsigned long _i1 = 0; _i1 < 16; _i1++) - _b[_i0][_i1] = 0.0; - } - - _a_raw[1] = 0.02173913043478260869565217391304347; - _a_raw[2] = 0.09629581047800066670113001679819925; - _a_raw[3] = 0.14444371571700100005169502519729888; - _a_raw[4] = 0.52205882352941176470588235294117647; - _a_raw[5] = 0.22842443612863469578031459099794265; - _a_raw[6] = 0.54360353589933733219171338103002937; - _a_raw[7] = 0.64335664335664335664335664335664335; - _a_raw[8] = 0.48251748251748251748251748251748251; - _a_raw[9] = 0.06818181818181818181818181818181818; - _a_raw[10] = 0.25060827250608272506082725060827250; - _a_raw[11] = 0.66736715965600568968278165443304378; - _a_raw[12] = 0.85507246376811594202898550724637681; - _a_raw[13] = 0.89795918367346938775510204081632653; - _a_raw[14] = 1.0; - _a_raw[15] = 1.0; - - _a[0] = 0.0; - for (unsigned long _i0 = 1; _i0 < 16; _i0++) - _a[_i0] = _a_raw[_i0] - _a_raw[_i0 - 1]; - - _b[1][0] = 1.0/46.0; - _b[2][0] =-0.11698050118114486205818241524969622; - _b[2][1] = 0.21327631165914552875931243204789548; - _b[3][0] = 0.03611092892925025001292375629932472; - _b[3][2] = 0.10833278678775075003877126889797416; - _b[4][0] = 1.57329743908138605107331820072051125; - _b[4][2] =-5.98400943754042002888532938159655553; - _b[4][3] = 4.93277082198844574251789353381722074; - _b[5][0] = 0.05052046351120380909008334360006234; - _b[5][3] = 0.17686653884807108146683657390397612; - _b[5][4] = 0.00103743376935980522339467349390418; - _b[6][0] = 0.10543148021953768958529340893598138; - _b[6][3] =-0.16042415162569842979496486916719383; - _b[6][4] = 0.11643956912829316045688724281285250; - _b[6][5] = 0.48215663817720491194449759844838932; - _b[7][0] = 0.07148407148407148407148407148407148; - _b[7][5] = 0.32971116090443908023196389566296464; - _b[7][6] = 0.24216141096813279233990867620960722; - _b[8][0] = 0.07162368881118881118881118881118881; - _b[8][5] = 0.32859867301674234161492268975519694; - _b[8][6] = 0.11622213117906185418927311444060725; - _b[8][7] =-0.03392701048951048951048951048951048; - _b[9][0] = 0.04861540768024729180628870095388582; - _b[9][5] = 0.03998502200331629058445317782406268; - _b[9][6] = 0.10715724786209388876739304914053506; - _b[9][7] =-0.02177735985419485163815426357369818; - _b[9][8] =-0.10579849950964443770179884616296721; - _b[10][0] =-0.02540141041535143673515871979014924; - _b[10][5] = 1.0/30.0; - _b[10][6] =-0.16404854760069182073503553020238782; - _b[10][7] = 0.03410548898794737788891414566528526; - _b[10][8] = 0.15836825014108792658008718465091487; - _b[10][9] = 0.21425115805975734472868683695127609; - _b[11][0] = 0.00584833331460742801095934302256470; - _b[11][5] =-0.53954170547283522916525526480339109; - _b[11][6] = 0.20128430845560909506500331018201158; - _b[11][7] = 0.04347222773254789483240207937678906; - _b[11][8] =-0.00402998571475307250775349983910179; - _b[11][9] = 0.16541535721570612771420482097898952; - _b[11][10] = 0.79491862412512344573322086551518180; - _b[12][0] =-0.39964965968794892497157706711861448; - _b[12][5] =-3.79096577568393158554742638116249372; - _b[12][6] =-0.40349325653530103387515807815498044; - _b[12][7] =-2.82463879530435263378049668286220715; - _b[12][8] = 1.04226892772185985533374283289821416; - _b[12][9] = 1.12510956420436603974237036536924078; - _b[12][10] = 3.32746188718986816186934832571938138; - _b[12][11] = 2.77897957186355606325818219255783627; - _b[13][0] = 0.39545306350085237157098218205756922; - _b[13][5] = 5.82534730759650564865380791881446903; - _b[13][6] =-0.36527452339161313311889856846974452; - _b[13][7] = 1.18860324058346533283780076203192232; - _b[13][8] = 0.57970467638357921347110271762687972; - _b[13][9] =-0.86824862589087693262676988867897834; - _b[13][10] =-5.20227677296454721392873650976792184; - _b[13][11] =-0.79895541420753382543211121058675915; - _b[13][12] = 0.14360623206363792632792463778889008; - _b[14][0] = 8.49173149061346398013352206978380938; - _b[14][5] = 86.32213734729036800877634194386790750; - _b[14][6] = 1.02560575501091662034511526187393241; - _b[14][7] = 85.77427969817339941806831550695235092; - _b[14][8] =-13.98699305104110611795532466113248067; - _b[14][9] =-20.71537405501426352265946477613161883; - _b[14][10] =-72.16597156619946800281180102605140463; - _b[14][11] =-76.71211139107806345587696023064419687; - _b[14][12] = 4.22319427707298828839851258893735507; - _b[14][13] =-1.25649850482823521641825667745565428; - _b[15][0] =-0.42892119881959353241190195318730008; - _b[15][5] =-9.16865700950084689999297912545025359; - _b[15][6] = 1.08317616770620939241547721530003920; - _b[15][7] =-1.23501525358323653198215832293981810; - _b[15][8] =-1.21438272617593906232943856422371019; - _b[15][9] = 1.37226168507232166621351243731869914; - _b[15][10] = 9.15723239697162418155377135344394113; - _b[15][11] = 1.30616301842220047563298585480401671; - _b[15][12] =-0.25285618808937955976690569433069974; - _b[15][13] = 0.38099910799663987066763679926508552; - - _c[0] = 0.01490902081978461022483617102382552; - _c[7] =-0.20408044692054151258349120934134791; - _c[8] = 0.22901438600570447264772469337066476; - _c[9] = 0.12800558251147375669208211573729202; - _c[10] = 0.22380626846054143649770066956485937; - _c[11] = 0.39553165293700054420552389156421651; - _c[12] = 0.05416646758806981196568364538360743; - _c[13] = 0.12691439652445903685643385312168037; - _c[14] =-0.00052539244262118876455834655383035; - _c[15] = 1.0/31.0; - - _cs[0] = 0.00653047880643482012034413441159249; - _cs[7] =-2.31471038197461347517552506241529830; - _cs[8] = 0.43528227238866280799530900822377013; - _cs[9] = 0.14907947287101933118545845390618763; - _cs[10] = 0.17905535442235532311850533252768020; - _cs[11] = 2.53400872222767706921176214508820825; - _cs[12] =-0.55430437423209112896721332268159015; - _cs[13] = 0.56924788787870083224213506297615260; - _cs[14] =-0.03644749690427461198884026816573513; - _cs[15] = 1.0/31.0; - - _d[0] = 1.0-_b[15][5]/_b[14][5]; - _d[1] = _b[15][0]-_b[14][0]*_b[15][5]/_b[14][5]; - _d[2] = _b[15][5]/_b[14][5]; - _d[3] = _b[15][6]-_b[14][6]*_b[15][5]/_b[14][5]; - _d[4] = _b[15][7]-_b[14][7]*_b[15][5]/_b[14][5]; - _d[5] = _b[15][8]-_b[14][8]*_b[15][5]/_b[14][5]; - _d[6] = _b[15][9]-_b[14][9]*_b[15][5]/_b[14][5]; - _d[7] = _b[15][10]-_b[14][10]*_b[15][5]/_b[14][5]; - _d[8] = _b[15][11]-_b[14][11]*_b[15][5]/_b[14][5]; - _d[9] = _b[15][12]-_b[14][12]*_b[15][5]/_b[14][5]; - _d[10] = _b[15][13]-_b[14][13]*_b[15][5]/_b[14][5]; - - do { - - do { - - - // Step 1 - - memcpy(_akafield_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - _segment1_calculate_nonconstant_ip_fields(_step, 1); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - _active_detuning_population = _akafield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(1); - - // Step 2 - - t += _a[1] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[1][0]*_akafield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 2); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-2); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(2); - - // Step 3 - - t += _a[2] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[2][0]*_akafield_detuning_population[_i0] + _b[2][1]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 3); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-3); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(3); - - // Step 4 - - t += _a[3] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[3][0]*_akafield_detuning_population[_i0] + _b[3][1]*_akbfield_detuning_population[_i0] - + _b[3][2]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 4); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-4); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(4); - - // Step 5 - - t += _a[4] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[4][0]*_akafield_detuning_population[_i0] + _b[4][1]*_akbfield_detuning_population[_i0] - + _b[4][2]*_akcfield_detuning_population[_i0] + _b[4][3]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 5); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-5); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(5); - - // Step 6 - - t += _a[5] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[5][0]*_akafield_detuning_population[_i0] + _b[5][3]*_akdfield_detuning_population[_i0] - + _b[5][4]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akifield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 6); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-6); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(6); - - // Step 7 - - t += _a[6] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _detuning_population[_i0] + _b[6][0]*_akafield_detuning_population[_i0] + _b[6][3]*_akdfield_detuning_population[_i0] - + _b[6][4]*_akefield_detuning_population[_i0] + _b[6][5]*_akifield_detuning_population[_i0]; - } - - - _active_detuning_population = _akjfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 7); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-7); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(7); - - // Step 8 - - t += _a[7] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akbfield_detuning_population[_i0] = _detuning_population[_i0] + _b[7][0]*_akafield_detuning_population[_i0] + _b[7][5]*_akifield_detuning_population[_i0] - + _b[7][6]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akbfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 8); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-8); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(8); - - // Step 9 - - t += _a[8] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akcfield_detuning_population[_i0] = _detuning_population[_i0] + _b[8][0]*_akafield_detuning_population[_i0] + _b[8][5]*_akifield_detuning_population[_i0] - + _b[8][6]*_akjfield_detuning_population[_i0]+ _b[8][7]*_akbfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akcfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 9); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-9); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(9); - - // Step 10 - - t += _a[9] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akdfield_detuning_population[_i0] = _detuning_population[_i0] + _b[9][0]*_akafield_detuning_population[_i0] + _b[9][5]*_akifield_detuning_population[_i0] - + _b[9][6]*_akjfield_detuning_population[_i0]+ _b[9][7]*_akbfield_detuning_population[_i0]+ _b[9][8]*_akcfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akdfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 10); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-10); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(10); - - // Step 11 - - t += _a[10] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akefield_detuning_population[_i0] = _detuning_population[_i0] + _b[10][0]*_akafield_detuning_population[_i0] + _b[10][5]*_akifield_detuning_population[_i0] - + _b[10][6]*_akjfield_detuning_population[_i0]+ _b[10][7]*_akbfield_detuning_population[_i0] + _b[10][8]*_akcfield_detuning_population[_i0] - + _b[10][9]*_akdfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akefield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 11); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-11); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(11); - - // Step 12 - - t += _a[11] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akffield_detuning_population[_i0] = _detuning_population[_i0] + _b[11][0]*_akafield_detuning_population[_i0] + _b[11][5]*_akifield_detuning_population[_i0] - + _b[11][6]*_akjfield_detuning_population[_i0] + _b[11][7]*_akbfield_detuning_population[_i0] + _b[11][8]*_akcfield_detuning_population[_i0] - + _b[11][9]*_akdfield_detuning_population[_i0] + _b[11][10]*_akefield_detuning_population[_i0]; - } - - - _active_detuning_population = _akffield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 12); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-12); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(12); - - // Step 13 - - t += _a[12] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akgfield_detuning_population[_i0] = _detuning_population[_i0] + _b[12][0]*_akafield_detuning_population[_i0] + _b[12][5]*_akifield_detuning_population[_i0] - + _b[12][6]*_akjfield_detuning_population[_i0]+ _b[12][7]*_akbfield_detuning_population[_i0] + _b[12][8]*_akcfield_detuning_population[_i0] - + _b[12][9]*_akdfield_detuning_population[_i0] + _b[12][10]*_akefield_detuning_population[_i0] + _b[12][11]*_akffield_detuning_population[_i0]; - } - - - _active_detuning_population = _akgfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 13); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-13); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(13); - - // Step 14 - - t += _a[13] * _step; - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akhfield_detuning_population[_i0] = _detuning_population[_i0] + _b[13][0]*_akafield_detuning_population[_i0] + _b[13][5]*_akifield_detuning_population[_i0] - + _b[13][6]*_akjfield_detuning_population[_i0]+ _b[13][7]*_akbfield_detuning_population[_i0] + _b[13][8]*_akcfield_detuning_population[_i0] - + _b[13][9]*_akdfield_detuning_population[_i0] + _b[13][10]*_akefield_detuning_population[_i0] + _b[13][11]*_akffield_detuning_population[_i0] - + _b[13][12]*_akgfield_detuning_population[_i0]; - } - - - _active_detuning_population = _akhfield_detuning_population; - - _segment1_calculate_nonconstant_ip_fields(_step, 14); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(-14); - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // a_i=D(a_2*dt)[y1] - _segment1_ip_evolve(14); - - // Step 15 and 16 combined to reduce memory use - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akifield_detuning_population[_i0] = _detuning_population[_i0] + _b[14][0]*_akafield_detuning_population[_i0] + _b[14][5]*_akifield_detuning_population[_i0] - + _b[14][6]*_akjfield_detuning_population[_i0]+ _b[14][7]*_akbfield_detuning_population[_i0] + _b[14][8]*_akcfield_detuning_population[_i0] - + _b[14][9]*_akdfield_detuning_population[_i0] + _b[14][10]*_akefield_detuning_population[_i0] + _b[14][11]*_akffield_detuning_population[_i0] - + _b[14][12]*_akgfield_detuning_population[_i0] + _b[14][13]*_akhfield_detuning_population[_i0]; - } - - - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akjfield_detuning_population[_i0] = _d[0]*_detuning_population[_i0] - + _d[1]*_akafield_detuning_population[_i0] - + _d[2]*_akifield_detuning_population[_i0] - + _d[3]*_akjfield_detuning_population[_i0] - + _d[4]*_akbfield_detuning_population[_i0] - + _d[5]*_akcfield_detuning_population[_i0] - + _d[6]*_akdfield_detuning_population[_i0] - + _d[7]*_akefield_detuning_population[_i0] - + _d[8]*_akffield_detuning_population[_i0] - + _d[9]*_akgfield_detuning_population[_i0] - + _d[10]*_akhfield_detuning_population[_i0]; - } - - - t += _a[14] * _step; - - _active_detuning_population = _akifield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - t += _a[15] * _step; - - _active_detuning_population = _akjfield_detuning_population; - - // a_k=G[a_k, t] - _segment1_calculate_delta_a(_step); - - // Take full step - - // ai = a - memcpy(_initial_detuning_population, _detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - // a = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _detuning_population[_i0] += _c[0]*_akafield_detuning_population[_i0] + _c[7]*_akbfield_detuning_population[_i0] + _c[8]*_akcfield_detuning_population[_i0] - + _c[9]*_akdfield_detuning_population[_i0] + _c[10]*_akefield_detuning_population[_i0] + _c[11]*_akffield_detuning_population[_i0] - + _c[12]*_akgfield_detuning_population[_i0] + _c[13]*_akhfield_detuning_population[_i0] + _c[14]*_akifield_detuning_population[_i0] - + _c[15]*_akjfield_detuning_population[_i0]; - } - - - // a* = a + etc - for (long _i0 = 0; _i0 < (_lattice_detuning) * _detuning_population_ncomponents; _i0++) { - _akafield_detuning_population[_i0] = _initial_detuning_population[_i0] + _cs[0]*_akafield_detuning_population[_i0] + _cs[7]*_akbfield_detuning_population[_i0] - + _cs[8]*_akcfield_detuning_population[_i0] + _cs[9]*_akdfield_detuning_population[_i0] + _cs[10]*_akefield_detuning_population[_i0] - + _cs[11]*_akffield_detuning_population[_i0] + _cs[12]*_akgfield_detuning_population[_i0] + _cs[13]*_akhfield_detuning_population[_i0] - + _cs[14]*_akifield_detuning_population[_i0] + _cs[15]*_akjfield_detuning_population[_i0]; - } - - - _active_detuning_population = _detuning_population; - - - - _error = 0.0; - - _segment1_detuning_population_error = _segment1_detuning_population_timestep_error(_akafield_detuning_population); - if (_segment1_detuning_population_error > _error) - _error = _segment1_detuning_population_error; - - _attempted_steps++; - - if (_error < _tolerance) { - _t_local += _step; - if (_step > _max_step) - _max_step = _step; - if (!_break_next && _step < _min_step) - _min_step = _step; - _discard = false; - } else { - t -= _step; - - if (_segment1_detuning_population_reset(_initial_detuning_population) == false) { - - _LOG(_WARNING_LOG_LEVEL, "WARNING: NaN present. Integration halted at t = %e.\n" - " Non-finite number in integration vector \"population\" in segment 1.\n", t); - if (_mg0_output_index_t < _mg0_output_lattice_t) - _mg0_sample(); - - goto _SEGMENT1_END; - } - - _segment1_ip_evolve(-1); - - _discard = true; - _break_next = false; - _unsuccessful_steps++; - } - - _old_step = _step; - - // Resize step - if (_error < 0.5*_tolerance || _error > _tolerance) { - const real _safetyFactor = 0.90; - real _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-0.7/9.0)) * pow(_last_norm_error, real(0.4/9.0)); - _scalingFactor = MAX(_scalingFactor, 1.0/5.0); - _scalingFactor = MIN(_scalingFactor, 7.0); - if (_error > _tolerance && _scalingFactor > 1.0) { - // If our step failed don't try and increase our step size. That would be silly. - _scalingFactor = _safetyFactor * pow(abs(_error/_tolerance), real(-1.0/9.0)); - } - _old_step = _step; - _last_norm_error = pow(_safetyFactor/_scalingFactor*pow(_last_norm_error, real(0.4/9.0)), real(9.0/0.7)); - _step *= _scalingFactor; - } - - } while (_discard); - - if (_break_next) { - if (_next_sample_flag[0]) { - _mg0_sample(); - _next_sample_counter[0]++; - } - if (_next_sample_flag[1]) - _next_sample_flag[2] = true; - else { - _break_next = false; - _t_break_next = _segment1_setup_sampling(_next_sample_flag, _next_sample_counter); - } - } - - if ( (_t_local + _step)*(1.0 + _EPSILON) > _t_break_next) { - _break_next = true; - _LOG(_SAMPLE_LOG_LEVEL, "Current timestep: %e\n", _old_step); - _step = _t_break_next - _t_local; - } - } while (!_next_sample_flag[2]); - - _SEGMENT1_END:; - xmds_free(_segment1_akafield_detuning_population); - xmds_free(_segment1_akbfield_detuning_population); - xmds_free(_segment1_akcfield_detuning_population); - xmds_free(_segment1_akdfield_detuning_population); - xmds_free(_segment1_akefield_detuning_population); - xmds_free(_segment1_akffield_detuning_population); - xmds_free(_segment1_akgfield_detuning_population); - xmds_free(_segment1_akhfield_detuning_population); - xmds_free(_segment1_akifield_detuning_population); - xmds_free(_segment1_akjfield_detuning_population); - xmds_free(_segment1_initial_detuning_population); - - _LOG(_SEGMENT_LOG_LEVEL, "Segment 1: minimum timestep: %e maximum timestep: %e\n", _min_step, _max_step); - _LOG(_SEGMENT_LOG_LEVEL, " Attempted %li steps, %.2f%% steps failed.\n", _attempted_steps, (100.0*_unsuccessful_steps)/_attempted_steps); -} - - -inline void _segment1_ip_evolve(int _exponent) -{ -} -real _segment1_setup_sampling(bool* _next_sample_flag, long* _next_sample_counter) -{ - // The numbers of the moment groups that need to be sampled at the next sampling point. - // An entry of N+1 means "reached end of integration interval" - long _momentGroupNumbersNeedingSamplingNext[2]; - long _numberOfMomentGroupsToBeSampledNext = 1; - - long _previous_m = 1; - long _previous_M = 1; - - real _t_break_next = (real)final_time; - _momentGroupNumbersNeedingSamplingNext[0] = 1; - - // initialise all flags to false - for (long _i0 = 0; _i0 < 2; _i0++) - _next_sample_flag[_i0] = false; - - /* Check if moment group needs sampling at the same time as another already discovered sample (or the final time). - * If so, add this moment group to the to-be-sampled list. If moment group demands sampling earlier than all - * previously noted moment groups, erase all previous ones from list and set the sample time to this earlier one. - */ - if (_next_sample_counter[0] * _previous_M == _previous_m * 10) { - _momentGroupNumbersNeedingSamplingNext[_numberOfMomentGroupsToBeSampledNext] = 0; - _numberOfMomentGroupsToBeSampledNext++; - } else if (_next_sample_counter[0] * _previous_M < _previous_m * 10) { - _t_break_next = _next_sample_counter[0] * ((real)final_time) / ((real)10); - _numberOfMomentGroupsToBeSampledNext = 1; - _momentGroupNumbersNeedingSamplingNext[0] = 0; - _previous_M = 10; - _previous_m = _next_sample_counter[0]; - } - - // _momentGroupNumbersNeedingSamplingNext now contains the complete list of moment groups that need - // to be sampled at the next sampling point. Set their flags to true. - for (long _i0 = 0; _i0 < _numberOfMomentGroupsToBeSampledNext; _i0++) - _next_sample_flag[_momentGroupNumbersNeedingSamplingNext[_i0]] = true; - - return _t_break_next; -} - -real _segment1_detuning_population_timestep_error(complex* _checkfield) -{ - real _error = 1e-24; - real _temp_error = 0.0; - real _temp_mod = 0.0; - - - // Find the peak value for each component of the field - real _cutoff[_detuning_population_ncomponents]; - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) - _cutoff[_i0] = 0.0; - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - _temp_mod = mod2(_detuning_population[_detuning_population_index_pointer + _i1]); - if (_xmds_isnonfinite(_temp_mod)) - _cutoff[_i1] = INFINITY; - else if (_cutoff[_i1] < _temp_mod) - _cutoff[_i1] = _temp_mod; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - for (long _i0 = 0; _i0 < _detuning_population_ncomponents; _i0++) { - if (_xmds_isnonfinite(_cutoff[_i0])) - // Return an error two times the tolerance in this case because the timestep must be reduced. - return 2.0*1e-10; - _cutoff[_i0] *= 0.001; - _cutoff[_i0] *= 0.001; - } - - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (mod2(_detuning_population[_detuning_population_index_pointer + _i1]) > _cutoff[_i1]) { - _temp_error = abs(_detuning_population[_detuning_population_index_pointer + _i1] - _checkfield[_detuning_population_index_pointer + _i1]) / (0.5*abs(_detuning_population[_detuning_population_index_pointer + _i1]) + 0.5*abs(_checkfield[_detuning_population_index_pointer + _i1])); - - if (_xmds_isnonfinite(_temp_error)) { - /* For _temp_error to be NaN, both the absolute value of the higher and lower order solutions - must BOTH be zero. This therefore implies that their difference is zero, and that there is no error. */ - _temp_error = 0.0; - } - - if (_error < _temp_error) // UNVECTORISABLE - _error = _temp_error; - } - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - - return _error; -} - -bool _segment1_detuning_population_reset(complex* _reset_to_detuning_population) -{ - memcpy(_detuning_population, _reset_to_detuning_population, sizeof(complex) * _detuning_population_alloc_size); - - /* return false if there's a NaN somewhere in the vector, otherwise return true */ - bool bNoNaNsPresent = true; - { - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - for (long _i0 = 0; _i0 < (_lattice_detuning); _i0++) { - for (long _i1 = 0; _i1 < _detuning_population_ncomponents; _i1++) { - if (_xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Re()) - || _xmds_isnonfinite(_detuning_population[_detuning_population_index_pointer + _i1].Im())) bNoNaNsPresent = false; - } - - _detuning_population_index_pointer += _detuning_population_ncomponents; - } - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - } - return bNoNaNsPresent; -} - -// Delta A propagation operator for field detuning -void _segment1_detuning_operators_evaluate_operator0(real _step) -{ - long _detuning_Gamma_index_pointer = 0; - #define r_ab _active_detuning_Gamma[_detuning_Gamma_index_pointer + 0] - #define r_cb _active_detuning_Gamma[_detuning_Gamma_index_pointer + 1] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - complex dPca_dt; - complex dPab_dt; - complex dPaa_dt; - complex dPbb_dt; - complex dPcc_dt; - complex dPcb_dt; - - #define dt _step - - // ************* Propagation code *************** - #line 92 "./eit3.xmds" - - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - - #line 1678 "./eit3.cc" - // ********************************************** - - #undef dt - - - _active_detuning_population[_detuning_population_index_pointer + 4] = dPca_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 3] = dPab_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 0] = dPaa_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 1] = dPbb_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 2] = dPcc_dt * _step; - _active_detuning_population[_detuning_population_index_pointer + 5] = dPcb_dt * _step; - // Increment index pointers for vectors in field detuning (or having the same dimensions) - _detuning_Gamma_index_pointer += 1 * _detuning_Gamma_ncomponents; - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef r_ab - #undef r_cb - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb -} - -// ******************************************************** -// output function implementations -void _write_output() -{ - _LOG(_SIMULATION_LOG_LEVEL, "Generating output for ./eit3\n"); - - - char *_xsilFilename = (char*)malloc(256); - snprintf(_xsilFilename, 256, "%s.xsil", ("./eit3" + gsArgsAndValues).c_str()); - - FILE* _outfile = _open_xsil_file(_xsilFilename); - - if (_outfile) { - _write_xsil_header(_outfile); - char _dataFilename[200]; - snprintf(_dataFilename, 200, "%s.h5", ("./eit3" + gsArgsAndValues).c_str()); - - H5Fclose(H5Fcreate(_dataFilename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT)); - } - _mg0_write_out(_outfile); - - _write_xsil_footer(_outfile); - _close_xsil_file(_outfile); - free(_xsilFilename); - _xsilFilename = NULL; - _outfile = NULL; - -} - - -FILE* _open_xsil_file(const char* _filename) -{ - - FILE* fp = fopen(_filename, "w"); - - if (fp == NULL) - // _LOG will cause the simulation to exit - _LOG(_ERROR_LOG_LEVEL, "Unable to open output file '%s'.\n" - "Exiting.\n", _filename); - - return fp; -} - -void _close_xsil_file(FILE*& fp) -{ - if (fp) - fclose(fp); - fp = NULL; - -} - -void _write_xsil_header(FILE* fp) -{ - if (!fp) - return; - fprintf(fp, "<?xml version=\"1.0\" ?><simulation xmds-version=\"2\">\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- While not strictly necessary, the following two tags are handy. -->\n"); - fprintf(fp, " <author>Hunter Rew</author>\n"); - fprintf(fp, " <description>\n"); - fprintf(fp, " Model of 3 state atom\n"); - fprintf(fp, " </description>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This element defines some constants. It can be used for other\n"); - fprintf(fp, " features as well, but we will go into that in more detail later.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <features>\n"); - fprintf(fp, " <precision> double </precision>\n"); - fprintf(fp, " <globals>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " real drive;\n"); - fprintf(fp, " real decay;\n"); - fprintf(fp, " real decay_bc;\n"); - fprintf(fp, " complex r_ca;\n"); - fprintf(fp, " real probe = .0001;\n"); - fprintf(fp, " real delta1;\n"); - fprintf(fp, " real dephase_bc;\n"); - fprintf(fp, " \n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </globals>\n"); - fprintf(fp, " <validation kind=\"run-time\"/>\n"); - fprintf(fp, " <arguments>\n"); - fprintf(fp, " <argument default_value=\"0.1\" name=\"drive_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10.0\" name=\"decay_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0.001\" name=\"decay_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"-1\" name=\"domain_min\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"1\" name=\"domain_max\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"delta1_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"0\" name=\"dephase_bc_arg\" type=\"real\"/>\n"); - fprintf(fp, " <argument default_value=\"10000\" name=\"final_time\" type=\"real\"/>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " drive = drive_arg;\n"); - fprintf(fp, " decay = decay_arg;\n"); - fprintf(fp, " decay_bc = decay_bc_arg;\n"); - fprintf(fp, " delta1 = delta1_arg;\n"); - fprintf(fp, " dephase_bc = dephase_bc_arg;\n"); - fprintf(fp, " r_ca = decay - i*delta1;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </arguments>\n"); - fprintf(fp, " </features>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " This part defines all of the dimensions used in the problem,\n"); - fprintf(fp, " in this case, only the dimension of 'time' is needed.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <geometry>\n"); - fprintf(fp, " <propagation_dimension> t </propagation_dimension>\n"); - fprintf(fp, " <transverse_dimensions>\n"); - fprintf(fp, " <dimension domain=\"(domain_min, domain_max)\" lattice=\"256\" name=\"detuning\"/>\n"); - fprintf(fp, " </transverse_dimensions>\n"); - fprintf(fp, " </geometry>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- A 'vector' describes the variables that we will be evolving. -->\n"); - fprintf(fp, " <vector name=\"Gamma\" type=\"complex\">\n"); - fprintf(fp, " <components> r_ab r_cb </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " r_ab = decay + i*(delta1 + detuning);\n"); - fprintf(fp, " r_cb = decay_bc + i*detuning + dephase_bc;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, " <vector dimensions=\"detuning\" name=\"population\" type=\"complex\">\n"); - fprintf(fp, " <components>\n"); - fprintf(fp, " Paa Pbb Pcc Pab Pca Pcb \n"); - fprintf(fp, " </components>\n"); - fprintf(fp, " <initialisation>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " Pcc = .5;\n"); - fprintf(fp, " Pbb = .5;\n"); - fprintf(fp, " Paa = Pab = Pca = Pcb = 0;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </initialisation>\n"); - fprintf(fp, " </vector>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <sequence>\n"); - fprintf(fp, " <!--\n"); - fprintf(fp, " Here we define what differential equations need to be solved\n"); - fprintf(fp, " and what algorithm we want to use.\n"); - fprintf(fp, " -->\n"); - fprintf(fp, " <integrate algorithm=\"ARK89\" interval=\"final_time\" tolerance=\"1e-10\">\n"); - fprintf(fp, " <samples>10</samples>\n"); - fprintf(fp, " <operators>\n"); - fprintf(fp, " <integration_vectors>population</integration_vectors>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb;\n"); - fprintf(fp, "\n"); - fprintf(fp, " dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab;\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " <dependencies>Gamma</dependencies>\n"); - fprintf(fp, " </operators>\n"); - fprintf(fp, " </integrate>\n"); - fprintf(fp, " </sequence>\n"); - fprintf(fp, "\n"); - fprintf(fp, " <!-- This part defines what data will be saved in the output file -->\n"); - fprintf(fp, " <output format=\"hdf5\">\n"); - fprintf(fp, " <sampling_group basis=\"detuning\" initial_sample=\"yes\">\n"); - fprintf(fp, " <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PaaR = Paa.Re();\n"); - fprintf(fp, " PaaI = Paa.Im();\n"); - fprintf(fp, " PbbR = Pbb.Re();\n"); - fprintf(fp, " PbbI = Pbb.Im();\n"); - fprintf(fp, " PccR = Pcc.Re();\n"); - fprintf(fp, " PccI = Pcc.Im();\n"); - fprintf(fp, " PabR = Pab.Re();\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " PcaR = Pca.Re();\n"); - fprintf(fp, " PcaI = Pca.Im();\n"); - fprintf(fp, " PcbR = Pcb.Re();\n"); - fprintf(fp, " PcbI = Pcb.Im();\n"); - fprintf(fp, " ]]> -->\n"); - fprintf(fp, " <moments>PabI</moments>\n"); - fprintf(fp, " <dependencies>population</dependencies>\n"); - fprintf(fp, " <![CDATA[\n"); - fprintf(fp, " PabI = Pab.Im();\n"); - fprintf(fp, " ]]>\n"); - fprintf(fp, " </sampling_group>\n"); - fprintf(fp, " </output>\n"); - - fprintf(fp, "\n<info>\n"); - fprintf(fp, "Script compiled with XMDS2 version 2.2.0 \"Out of cheese error\" (r2941)\n"); - fprintf(fp, "See http://www.xmds.org for more information.\n"); - fprintf(fp, "\nVariables that can be specified on the command line:\n"); - - fprintf(fp, " Command line argument drive_arg = %e\n", drive_arg); - - fprintf(fp, " Command line argument decay_arg = %e\n", decay_arg); - - fprintf(fp, " Command line argument decay_bc_arg = %e\n", decay_bc_arg); - - fprintf(fp, " Command line argument domain_min = %e\n", domain_min); - - fprintf(fp, " Command line argument domain_max = %e\n", domain_max); - - fprintf(fp, " Command line argument delta1_arg = %e\n", delta1_arg); - - fprintf(fp, " Command line argument dephase_bc_arg = %e\n", dephase_bc_arg); - - fprintf(fp, " Command line argument final_time = %e\n", final_time); - fprintf(fp, "</info>\n"); - -} - -// In addition to writing the footer (if 'fp' is not NULL) -// this function closes the fp file pointer. -void _write_xsil_footer(FILE* fp) -{ - if (fp) { - fprintf(fp, "</simulation>\n"); - } -} - -// ******************************************************** -// moment group 0 function implementations -void _mg0_sample() -{ - - long _mg0_output_raw_index_pointer = 0; - #define PabI _active_mg0_output_raw[_mg0_output_raw_index_pointer + 0] - long _detuning_population_index_pointer = 0; - #define Paa _active_detuning_population[_detuning_population_index_pointer + 0] - #define Pbb _active_detuning_population[_detuning_population_index_pointer + 1] - #define Pcc _active_detuning_population[_detuning_population_index_pointer + 2] - #define Pab _active_detuning_population[_detuning_population_index_pointer + 3] - #define Pca _active_detuning_population[_detuning_population_index_pointer + 4] - #define Pcb _active_detuning_population[_detuning_population_index_pointer + 5] - #define detuning _detuning[_index_detuning + 0] - #define ddetuning (_ddetuning * (1.0)) - - for (long _index_detuning = 0; _index_detuning < _lattice_detuning; _index_detuning++) { - // Set index pointers explicitly for (some) vectors - _mg0_output_raw_index_pointer = ( 0 - + _mg0_output_index_t * _lattice_detuning - + _index_detuning * 1 ) * _mg0_output_raw_ncomponents; - #define _SAMPLE_COMPLEX(variable) \ - variable ## R = variable.Re(); variable ## I = variable.Im(); - - // *************** Sampling code **************** - #line 131 "./eit3.xmds" - - PabI = Pab.Im(); - - #line 1960 "./eit3.cc" - // ********************************************** - - #undef _SAMPLE_COMPLEX - // Increment index pointers for vectors in field mg0_sampling (or having the same dimensions) - _detuning_population_index_pointer += 1 * _detuning_population_ncomponents; - - } - #undef detuning - #undef ddetuning - #undef PabI - #undef Paa - #undef Pbb - #undef Pcc - #undef Pab - #undef Pca - #undef Pcb - - _mg0_output_t[0 + _mg0_output_index_t++] = t; - - _LOG(_SAMPLE_LOG_LEVEL, "Sampled field (for moment group #1) at t = %e\n", t); - -} - - -void _mg0_process() -{ - // No post processing needs to be done -} - - -void _mg0_write_out(FILE* _outfile) -{ - - if (_outfile) { - fprintf(_outfile, "\n"); - fprintf(_outfile, "<XSIL Name=\"moment_group_1\">\n"); - fprintf(_outfile, " <Param Name=\"n_independent\">2</Param>\n"); - fprintf(_outfile, " <Array Name=\"variables\" Type=\"Text\">\n"); - fprintf(_outfile, " <Dim>3</Dim>\n"); - fprintf(_outfile, " <Stream><Metalink Format=\"Text\" Delimiter=\" \\n\"/>\n"); - fprintf(_outfile, "t detuning PabI \n"); - fprintf(_outfile, " </Stream>\n"); - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, " <Array Name=\"data\" Type=\"double\">\n"); - fprintf(_outfile, " <Dim>%i</Dim>\n", _mg0_output_lattice_t); - fprintf(_outfile, " <Dim>%i</Dim>\n", _lattice_detuning); - fprintf(_outfile, " <Dim>3</Dim>\n"); - } - - - char _h5Filename[200]; - snprintf(_h5Filename, 200, "%s.h5", ("./eit3" + gsArgsAndValues).c_str()); - - /* Open the file */ - hid_t hdf5_file = H5Fopen(_h5Filename, H5F_ACC_RDWR, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_WARNING_LOG_LEVEL, "Failed to open HDF5 file '%s', will try to create it.", _h5Filename); - hdf5_file = H5Fcreate(_h5Filename, H5F_ACC_EXCL, H5P_DEFAULT, H5P_DEFAULT); - if (hdf5_file < 0) { - _LOG(_ERROR_LOG_LEVEL, "Failed to create HDF5 file '%s'. Bailing.", _h5Filename); - } - } - - /* Create the group for this data */ - hid_t group; - if (!H5Lexists(hdf5_file, "/1", H5P_DEFAULT)) - group = H5Gcreate(hdf5_file, "/1", H5P_DEFAULT); - else - group = H5Gopen(hdf5_file, "/1"); - - if (_outfile) { - fprintf(_outfile, " <Stream><Metalink Format=\"HDF5\" Type=\"Remote\" Group=\"/1\"/>\n"); - fprintf(_outfile, "%s.h5\n", ("./eit3" + gsArgsAndValues).c_str()); - fprintf(_outfile, " </Stream>\n"); - } - - /* Create the coordinate data sets */ - hsize_t coordinate_length; - hid_t coordinate_dataspace; - coordinate_length = _mg0_output_lattice_t; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_t; - if (!H5Lexists(hdf5_file, "/1/t", H5P_DEFAULT)) - dataset_t = H5Dcreate(hdf5_file, "/1/t", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_t = H5Dopen(hdf5_file, "/1/t"); - H5Dwrite(dataset_t, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _mg0_output_t); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_t, "t"); - #endif - - H5Sclose(coordinate_dataspace); - coordinate_length = _lattice_detuning; - coordinate_dataspace = H5Screate_simple(1, &coordinate_length, NULL); - hid_t dataset_detuning; - if (!H5Lexists(hdf5_file, "/1/detuning", H5P_DEFAULT)) - dataset_detuning = H5Dcreate(hdf5_file, "/1/detuning", H5T_NATIVE_REAL, coordinate_dataspace, H5P_DEFAULT); - else - dataset_detuning = H5Dopen(hdf5_file, "/1/detuning"); - H5Dwrite(dataset_detuning, H5T_NATIVE_REAL, H5S_ALL, H5S_ALL, H5P_DEFAULT, _detuning); - #if defined(HAVE_HDF5_HL) - H5DSset_scale(dataset_detuning, "detuning"); - #endif - - H5Sclose(coordinate_dataspace); - - hsize_t file_dims[] = {_mg0_output_lattice_t, _lattice_detuning}; - hid_t file_dataspace = H5Screate_simple(2, file_dims, NULL); - - hid_t dataset_PabI; - if (!H5Lexists(hdf5_file, "/1/PabI", H5P_DEFAULT)) - dataset_PabI = H5Dcreate(hdf5_file, "/1/PabI", H5T_NATIVE_REAL, file_dataspace, H5P_DEFAULT); - else - dataset_PabI = H5Dopen(hdf5_file, "/1/PabI"); - #if defined(HAVE_HDF5_HL) - H5DSattach_scale(dataset_PabI, dataset_t, 0); - H5DSattach_scale(dataset_PabI, dataset_detuning, 1); - #endif - H5Dclose(dataset_t); - H5Dclose(dataset_detuning); - - - if ((_mg0_output_lattice_t * _lattice_detuning)) { - /* Create the data space */ - hsize_t file_start[2] = {0, 0}; - hsize_t mem_dims[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - hsize_t mem_start[3] = {0, 0, 0}; - hsize_t mem_stride[3] = {1, 1, 1}; - hsize_t mem_count[3] = {_mg0_output_lattice_t, _lattice_detuning, 1}; - - - hid_t mem_dataspace; - mem_dims[2] = 1; - mem_dataspace = H5Screate_simple(3, mem_dims, NULL); - mem_stride[2] = 1; - - // Select hyperslabs of memory and file data spaces for data transfer operation - mem_start[2] = 0; - H5Sselect_hyperslab(mem_dataspace, H5S_SELECT_SET, mem_start, mem_stride, mem_count, NULL); - H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET, file_start, mem_stride, mem_count, NULL); - - if (dataset_PabI) - H5Dwrite(dataset_PabI, H5T_NATIVE_REAL, mem_dataspace, file_dataspace, H5P_DEFAULT, _active_mg0_output_raw); - - H5Sclose(mem_dataspace); - } - - - H5Dclose(dataset_PabI); - - H5Sclose(file_dataspace); - H5Gclose(group); - H5Fclose(hdf5_file); - - - if (_outfile) { - fprintf(_outfile, " </Array>\n"); - fprintf(_outfile, "</XSIL>\n"); - } -} - -// ******************************************************** -// field mg0_output function implementations -// initialisation for vector raw -void _mg0_output_raw_initialise() -{ - - bzero(_active_mg0_output_raw, sizeof(real) * _mg0_output_raw_alloc_size); -} - diff --git a/source/eit3.h5 b/source/eit3.h5 Binary files differdeleted file mode 100644 index 1274852..0000000 --- a/source/eit3.h5 +++ /dev/null diff --git a/source/eit3.py b/source/eit3.py deleted file mode 100644 index ac3ff5e..0000000 --- a/source/eit3.py +++ /dev/null @@ -1,17 +0,0 @@ -#!/usr/bin/env python -from xpdeint.XSILFile import XSILFile - -xsilFile = XSILFile("./eit3.xsil") - -def firstElementOrNone(enumerable): - for element in enumerable: - return element - return None - -t_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "t") -detuning_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].independentVariables if _["name"] == "detuning") -PabI_1 = firstElementOrNone(_["array"] for _ in xsilFile.xsilObjects[0].dependentVariables if _["name"] == "PabI") - -# Write your plotting commands here. -# You may want to import pylab (from pylab import *) or matplotlib (from matplotlib import *) - diff --git a/source/eit3.pyc b/source/eit3.pyc Binary files differdeleted file mode 100644 index e03ac5b..0000000 --- a/source/eit3.pyc +++ /dev/null diff --git a/source/eit3.xmds b/source/eit3.xmds deleted file mode 100644 index 50419f5..0000000 --- a/source/eit3.xmds +++ /dev/null @@ -1,136 +0,0 @@ -<?xml version="1.0" encoding="UTF-8"?> -<simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument name="drive_arg" type="real" default_value="0.1" /> - <argument name="decay_arg" type="real" default_value="10.0" /> - <argument name="decay_bc_arg" type="real" default_value="0.001" /> - <argument name="domain_min" type="real" default_value="-1" /> - <argument name="domain_max" type="real" default_value="1" /> - <argument name="delta1_arg" type="real" default_value="0" /> - <argument name="dephase_bc_arg" type="real" default_value="0" /> - <argument name="final_time" type="real" default_value="10000" /> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension name="detuning" lattice="256" domain="(domain_min, domain_max)" /> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector name="population" type="complex" dimensions="detuning"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5" > - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> -</simulation> diff --git a/source/eit3.xsil b/source/eit3.xsil deleted file mode 100644 index 9bef9a3..0000000 --- a/source/eit3.xsil +++ /dev/null @@ -1,168 +0,0 @@ -<?xml version="1.0" ?><simulation xmds-version="2"> - - <!-- While not strictly necessary, the following two tags are handy. --> - <author>Hunter Rew</author> - <description> - Model of 3 state atom - </description> - - <!-- - This element defines some constants. It can be used for other - features as well, but we will go into that in more detail later. - --> - <features> - <precision> double </precision> - <globals> - <![CDATA[ - real drive; - real decay; - real decay_bc; - complex r_ca; - real probe = .0001; - real delta1; - real dephase_bc; - - ]]> - </globals> - <validation kind="run-time"/> - <arguments> - <argument default_value="0.1" name="drive_arg" type="real"/> - <argument default_value="10.0" name="decay_arg" type="real"/> - <argument default_value="0.001" name="decay_bc_arg" type="real"/> - <argument default_value="-1" name="domain_min" type="real"/> - <argument default_value="1" name="domain_max" type="real"/> - <argument default_value="0" name="delta1_arg" type="real"/> - <argument default_value="0" name="dephase_bc_arg" type="real"/> - <argument default_value="10000" name="final_time" type="real"/> - <![CDATA[ - drive = drive_arg; - decay = decay_arg; - decay_bc = decay_bc_arg; - delta1 = delta1_arg; - dephase_bc = dephase_bc_arg; - r_ca = decay - i*delta1; - ]]> - </arguments> - </features> - - <!-- - This part defines all of the dimensions used in the problem, - in this case, only the dimension of 'time' is needed. - --> - <geometry> - <propagation_dimension> t </propagation_dimension> - <transverse_dimensions> - <dimension domain="(domain_min, domain_max)" lattice="256" name="detuning"/> - </transverse_dimensions> - </geometry> - - <!-- A 'vector' describes the variables that we will be evolving. --> - <vector name="Gamma" type="complex"> - <components> r_ab r_cb </components> - <initialisation> - <![CDATA[ - r_ab = decay + i*(delta1 + detuning); - r_cb = decay_bc + i*detuning + dephase_bc; - ]]> - </initialisation> - </vector> - <vector dimensions="detuning" name="population" type="complex"> - <components> - Paa Pbb Pcc Pab Pca Pcb - </components> - <initialisation> - <![CDATA[ - Pcc = .5; - Pbb = .5; - Paa = Pab = Pca = Pcb = 0; - ]]> - </initialisation> - </vector> - - <sequence> - <!-- - Here we define what differential equations need to be solved - and what algorithm we want to use. - --> - <integrate algorithm="ARK89" interval="final_time" tolerance="1e-10"> - <samples>10</samples> - <operators> - <integration_vectors>population</integration_vectors> - <![CDATA[ - dPcc_dt = i*drive*(-Pca) - i*drive*Pca + decay*Paa - decay_bc*Pcc + decay_bc*Pbb; - - dPbb_dt = i*probe*Pab - i*probe*(-Pab) + decay*Paa - decay_bc*Pbb + decay_bc*Pcc; - - dPaa_dt = -i*probe*Pab + i*probe*(-Pab) - i*drive*(-Pca) + i*drive*Pca - 2*decay*Paa; - - dPab_dt = -r_ab*Pab + i*probe*(Pbb-Paa) + i*drive*Pcb; - - dPca_dt = -r_ca*Pca + i*drive*(Paa-Pcc) - i*probe*Pcb; - - dPcb_dt = -r_cb*Pcb - i*probe*Pca + i*drive*Pab; - ]]> - <dependencies>Gamma</dependencies> - </operators> - </integrate> - </sequence> - - <!-- This part defines what data will be saved in the output file --> - <output format="hdf5"> - <sampling_group basis="detuning" initial_sample="yes"> - <!-- <moments>PaaR PaaI PbbR PbbI PccR PccI PabR PabI PcaR PcaI PcbR PcbI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PaaR = Paa.Re(); - PaaI = Paa.Im(); - PbbR = Pbb.Re(); - PbbI = Pbb.Im(); - PccR = Pcc.Re(); - PccI = Pcc.Im(); - PabR = Pab.Re(); - PabI = Pab.Im(); - PcaR = Pca.Re(); - PcaI = Pca.Im(); - PcbR = Pcb.Re(); - PcbI = Pcb.Im(); - ]]> --> - <moments>PabI</moments> - <dependencies>population</dependencies> - <![CDATA[ - PabI = Pab.Im(); - ]]> - </sampling_group> - </output> - -<info> -Script compiled with XMDS2 version 2.2.0 "Out of cheese error" (r2941) -See http://www.xmds.org for more information. - -Variables that can be specified on the command line: - Command line argument drive_arg = 9.111628e-01 - Command line argument decay_arg = 6.000000e+00 - Command line argument decay_bc_arg = 1.000000e-05 - Command line argument domain_min = -1.383796e+00 - Command line argument domain_max = 1.383796e+00 - Command line argument delta1_arg = -1.500000e+01 - Command line argument dephase_bc_arg = 0.000000e+00 - Command line argument final_time = 7.226499e+02 -</info> - -<XSIL Name="moment_group_1"> - <Param Name="n_independent">2</Param> - <Array Name="variables" Type="Text"> - <Dim>3</Dim> - <Stream><Metalink Format="Text" Delimiter=" \n"/> -t detuning PabI - </Stream> - </Array> - <Array Name="data" Type="double"> - <Dim>11</Dim> - <Dim>256</Dim> - <Dim>3</Dim> - <Stream><Metalink Format="HDF5" Type="Remote" Group="/1"/> -./eit3.h5 - </Stream> - </Array> -</XSIL> -</simulation> |