imaging.m
Fokker-Planck imaging simulation context. Generates the Hamiltonian, the relaxation superoperator, the kinetics superoperator, the Fokker-Planck spatial dynamics generator (including diffusion and flow), gradient operators, and passes all of that to the pulse sequence, which should be supplied as a handle.
Syntax
answer=imaging(spin_system,pulse_sequence,parameters)
Parameters
pulse_sequence - pulse sequence function handle. See the
experiments directory for the list of
pulse sequences that ship with Spinach.
parameters.u - X components of the velocity vectors
for each point in the sample, m/s
parameters.v - Y components of the velocity vectors
for each point in the sample, m/s
parameters.w - Z components of the velocity vectors
for each point in the sample, m/s
parameters.diff - diffusion coefficient or 3x3 tensor, m^2/s
for situations when this parameter is the
same in every voxel
parameters.dxx - Cartesian components of the diffusion
parameters.dxy tensor for each voxel of the sample
...
parameters.dzz
parameters.dims - dimensions of the 3D box, meters
parameters.npts - number of points in each dimension
of the 3D box
parameters.deriv - {'fourier'} uses Fourier diffe-
rentiation matrices; {'period',n}
requests n-point central finite-
difference matrices with periodic
boundary conditions
Three types of phantoms must be specified. The relaxation theory phantom contains relaxation superoperators and their coefficients in each voxel, specified in the following way:
parameters.rlx_ph={Ph1,Ph2,...,PhN}
parameters.rlx_op={R1,R2,...,RN}
where PhN have the same dimension as the sample voxel grid and RN are relaxation superoperators. The initial condition phantom reflects the fact that different voxels might start off in a different spin state. It must be specified in the following way:
parameters.rho0_ph={Ph1,Ph2,...,PhN}
parameters.rho0_op={rho1,rho2,...,rhoN}
where PhN have the same dimension as the sample voxel grid and rhoN are spin states obtained from state() function. The detection state phantom reflects the fact that different voxels might be detected at different ngles and with different sensitivity. It must be specified in the following way:
parameters.coil_ph={Ph1,Ph2,...,PhN}
parameters.coil_op={rho1,rho2,...,rhoN}
where PhN have the same dimension as the sample voxel grid and rhoN are spin states obtained from state() function.
Outputs
This function returns whatever the pulse sequence returns.
Notes
- The direct product order is Z(x)Y(x)X(x)Spin, this corresponds to a column-wise vectorization of a 3D array with dimensions ordered as [X Y Z].
- Very large imaging calculations always benefit from Tesla cards (add 'gpu' to sys.enable) and may benefit from polyadic array processing (add 'polyadic' to sys.enable).
- Fourier derivatives are precise but expensive. If your accuracy requirements are not seven decimal places, consider using finite difference derivatives.
See also
crystal.m, doublerot.m, floquet.m, gridfree.m, liquid.m, meshflow.m, powder.m, singlerot.m, Kernel contexts, Built-in_experiments
Version 2.8, authors: Ilya Kuprov, Ahmed Allami, Maria Grazia Concilio