Interested to model diffusion in garnet using a simple MATLAB script? The full multi-component case in 3D is available here. Please don’t judge: it’s the very first code I had written many years ago. But maybe it’s still useful for you?
diffuseo9.m is a MATLAB script designed to simulate multi-component diffusion (Mn, Fe, Mg, Ca) in garnet at a constant pressure, temperature, and fO2 assuming that the rim composition of garnet is fixed during the diffusion (open system). The kinetic data are those of Chakraborty and Ganguly (1992) with the option to modify fO2; by default, equilibration with graphite is assumed. Input is a compositional profile (core to rim), with the distance data (in cm) contained in vector Xnodecm, and the corresponding compositional data in the vectors Xmn, Xfe, and Xmg. Output from Theria_G (e.g., garnet_gen001.txt) can be used as input (and is included as example). However, other input can also be used if the script is adjusted.
Run the script and enter the required information into the Command Window (temperature, pressure, fO2, duration of diffusion). Also, when requested by the script, specify in the Command Window the spacing between isochrones to be plotted, as well as the size of the time step for the simulation. Ideally, the size of the time step should not exceed the number given. However, it is advisable to start the modelling with larger time steps to obtain first estimates quickly, and to decrease these time steps to the size given for final results.
This MATLAB script creates the files spss_profiles.mat, alm_profiles.mat, grs_profiles.mat, and py_profiles.mat which can be edited with text editors and spreadsheet software. These files contain the calculated compositional profiles after diffusion, with the first column containing the duration of diffusion (in years) and the first row containing the radius (in cm). The files will be replaced during each run of this script.
This model is part of Theria_G (Gaidies et al., 2008), so please refer to it if this script turns out to be useful for you. Happy modelling!
Update (June 2022):
Some people have expressed an interest to use this multi-component diffusion script in their teaching, so I updated it a little (diffuse_oMC.m). I also added a script to model the 1-component case (diffuse_o1C.m). Both scripts can be found here: diffuse_o.zip. The 1C model allows to quickly account for the influence of pre-exponential constant and activation energy, which may be useful when testing proposed Y,REE diffusion data. As always: feedback and questions are welcome.