****************************************************************************************************************************
File for modeling Fe-Mg diffusion in olivine using finite differences implicit method (e.g. Crank, 197?)
This program is associated with the manuscript Costa, F., Dohmen, R., Chakraborty, S. (2008) MSA
short course and the uses the finite differences metod and equations explained in Appendix I.

Below the variables are 'cleared' to avoid memory effects from previous calculations.
    ******************************************************************************************************************************

Clear[Fo, PRES, T, fO2, ax, bx, cx, NNO, deltaNNO,DFoc, DFox,DFoa, DFob,DFocd,DFoxd,DFoad, DFobd,R1,δt,δx,datafo, Fomin, Fomax,ttgrid, datanfo,cFo,bFo,DIFFo,datafop, iniFop, datanfop,eq1,timeinsec ]

******************************************************************************************************
Temperature (T) in celsius,  Pressure (PRES) in bars, oxygen fugacity (fO2) in bars, forsterite content (Fo) in mol %.
 The value of Fo found below is the lowest found in the profile is used for stability check.

Distance in micrometers, time in seconds, (a,b,c)x = Cos of the angle between the electron microprobe traverse and the (a,b,c) axes 
of the crystal. fO2 shown with respect to the Log NNO oxygen buffer of Huebner and Sato (197?) 
*****************************************************************************************************

PRES=1.;T=1220.;fO2=10^-6.33;Fo=72.;ax=Cos[90Degree];bx=Cos[90Degree];cx=Cos[0.0Degree];NNO=((-24930.)/((273.+T)))+9.36+((0.046*(PRES-1))/((273.+T)));deltaNNO=Log[10,fO2]-NNO

**********************************************************************************************
    Diffusion coefficient data for Fe-Mg (DFo) from Dohmen and Chakraborty (2007a), with corrections of Dohmen and Chakraborty (2007b).
    Using diffusion parameters for fO2> (10^(-15)) bar, for lower fO2 see publications.
    
    DFo(a,b,c) = Fe-Mg diffusion coefficient parallel to (a,b,c) axis in ((μm)^2)(s^(-1)). 
    The only variable that changes with time is the concentration and thus we have created a new set of diffusion 
    coefficients without the dependence and this will be directly incorporated in the diffusion loop.
    DFox= diffusion coefficient parallel to the electron microprobe traverse.
    DFoxd= expression of DFo without the compositional dependence used to speed
    up the loop below.
    **********************************************************************************************

DFoc=(10^(-9.21))*((((fO2/(10^(-12))) ))^(1/6))*(10^(3*(0.9-0.01*Fo)))*Exp[((-(201000+(PRES-1)*7*(10^(-6))*(10^5)))/(8.314*(T+273.15)))]*(10^12);DFoa=DFoc/6.;DFob=DFoc/6.;

DFocd=(10^(-9.21))*((((fO2/(10^(-12))) ))^(1/6))*Exp[((-(201000+(PRES-1)*7*(10^(-6))*(10^5)))/(8.314*(T+273.15)))]*(10^12);DFoad=DFocd/6.;DFobd=DFocd/6.;DFoxd=DFoad*(ax^2)+DFobd*(bx^2)+DFocd*(cx^2); DFox=DFoxd*(10^(3*(0.9-0.01*Fo)))

**********************************************************************************************
Here we need to choose a time step (in seconds) = δt , and a distance step in μm = δx. The total lenght of the profile
will then be δx * ttgrid (= total number of grispoints), and the total time for the calculation will be δt * ttiter (= number
of time iterations, see below). 

The stability criteria for the implicit method is that ((δt)/((δx)^2))DFox < 0.5. R1= dummy variable to speed up the loop below.
************************************************************************************************

δx = 5.;δt =5.*(10^3);R1=((δt)/((δx)^2));R1*DFox

*****************************************************************************************************************************
Importing and plotting data electron microprobe data (Fomin= minimum Fo; Fomax= maximum Fo). Has to be the form of text with the first column distance and the second concentration.
datafo = natural data. Below the data has been already imported (you can use the command to import your file) so don't worry about error message.
******************************************************************************************************************************

Import["E:\mikeolis\ol1t6fon.txt","Table"];

datafo={{0,82.21674401`},{9.2`,82.0097873`},{18.4`,82.05928987`},{27.6`,82.28215166`},{36.8`,82.18992248`},{46,82.12052165`},{55.2`,82.19986456`},{64.4`,82.04918427`},{73.6`,82.08912176`},{82.8`,82.14371661`},{92,82.12896324`},{101.2`,82.11251099`},{110.4`,81.99594438`},{119.6`,82.02138934`},{128.8`,81.95267698`},{138,82.23249988`},{147.2`,81.9788329`},{156.4`,81.97165385`},{165.6`,82.04087549`},{174.8`,82.17687733`},{184,82.3232813`},{193.2`,82.16433838`},{202.4`,82.02996617`},{211.6`,81.989791`},{220.8`,82.08355266`},{230,82.13975904`},{239.2`,82.15998837`},{248.4`,81.93029491`},{257.6`,82.14613416`},{266.8`,81.98421975`},{276,82.04830918`},{285.2`,81.88675596`},{294.4`,82.08563938`},{303.6`,81.74396302`},{312.8`,82.01037676`},{331.2`,82.15008188`},{340.4`,82.16704289`},{349.6`,82.0232366`},{358.8`,81.93140355`},{368,82.19950822`},{377.2`,82.03162899`},{386.4`,82.07638152`},{395.6`,81.88294186`},{404.8`,82.08566467`},{414,81.97837212`},{423.2`,82.0214454`},{432.4`,82.02409639`},{441.6`,82.00483092`},{450.8`,81.79832014`},{460,82.24783862`},{469.2`,81.94518554`},{478.4`,82.08609112`},{487.6`,82.25790848`},{496.8`,82.29727762`},{506,81.91788799`},{515.2`,81.86059936`},{524.4`,81.74890272`},{533.6`,82.09575342`},{542.8`,81.73493976`},{552,82.0626935`},{561.2`,82.06334578`},{570.4`,82.10556845`},{579.6`,81.82390851`},{588.8`,81.93358432`},{598,81.96998214`},{607.2`,81.95786625`},{616.4`,81.90835632`},{625.6`,81.76518765`},{634.8`,81.94417523`},{644,81.9788156`},{653.2`,81.94264303`},{662.4`,81.80938558`},{671.6`,81.68952859`},{680.8`,81.69649283`},{690,81.90673675`},{699.2`,81.82522017`},{708.4`,81.79232549`},{717.6`,81.60671463`},{726.8`,81.61605467`},{736,81.69352899`},{745.2`,81.78088803`},{754.4`,81.71694185`},{763.6`,81.55956679`},{772.8`,81.6035463`},{782,81.5419655`},{791.2`,81.62724601`},{800.4`,81.25`},{809.6`,81.16939365`},{818.8`,81.2841464`},{828,81.35707411`},{837.2`,81.17067992`},{864.8`,80.46758255`},{874,80.49171566`},{892.4`,80.0935399`},{901.6`,80.17265914`},{910.8`,79.72044652`},{920,79.77609419`},{938.4`,79.32611312`},{947.6`,79.10821885`},{956.8`,79.16949976`},{966,78.81899253`},{975.2`,78.67401712`},{1002.8`,77.82932251`},{1012,77.73826116`},{1021.2`,77.50518147`},{1030.4`,77.50871755`},{1039.6`,76.8151134`},{1048.8`,76.84119678`},{1058,76.49158347`},{1067.2`,75.99270178`},{1076.4`,75.93207293`},{1085.6`,75.47844578`},{1094.8`,75.21995553`},{1104,74.95055714`},{1113.2`,74.61948619`},{1131.6`,73.90651489`},{1140.8`,73.64374819`},{1150,73.09521499`},{1159.2`,72.73431448`}};ol9=ListPlot[datafo,PlotStyle ->RGBColor[0,0,1],PlotJoined->True];Fomin=70.;Fomax=85;datafop=ListPlot[datafo,PlotRange->{Fomin-3,Fomax+3}, PlotStyle ->RGBColor[0,0,1],PlotJoined->True]

*****************************************************************************************************************************
Choosing the total number of gridpoints (ttgrid) which will be the total lenght of the profile in μm  (=ttlenght) divided by δx. 
The larger the number of gridpoints the more detailed will be the result but the longer time it will take
to do the calculation. The number of gridpoints has to be an integer number so the command Round is applied.

Other lines normalize and plot the natural data to the number of gridpoints we have chosen (=datanfo).
******************************************************************************************************************************

ttlenght=1160;ttgrid=Round[ttlenght/δx]

datanfo=datafo;datanfo[[All,1]]= datafo[[All,1]]/δx ;datanfop=ListPlot[datanfo,PlotRange->{Fomin-3,Fomax+3},PlotStyle ->RGBColor[0,0,1],PlotJoined->True]

***************************************************************************************************************************
Choosing the initial distribution (in this case homogeneous and equal to Foini) of Fo (= cFo) and seeting the
diffusion coefficients for the initial distribution (= DIFFo)
****************************************************************************************************************************

Foini=82.;cFo=Table[Foini,{ttgrid}];iniFop=ListPlot[cFo, PlotRange->{Fomin-3,Fomax+3}];Show[iniFop,datanfop]

DIFFo=Table[0.,{ttgrid}];Do[DIFFo[[i]]=DFoxd*(10^(3*(0.9-0.01*cFo[[i]]))), {i,ttgrid}];ListPlot[DIFFo]

*****************************************************************************************************************************
Diffusion algorithm

The first Do loop is for the time and executes the number of iterations (ttiter) that we have chosen.
It is typically better to set 'ttiter' to a moderate number and then hit shift + return a few times rather
than having a very large 'ttiter' value and only doing shift + return once. One needs to keep track
of how many times one hits shift + return, since total number of iterations will be this number
times 'ttiter'.
 
The first equality sets the variable bFo which stores the concentration of Fo after every time step. 
The two equalities that follow are the boundary conditions; For this case the left boundary is  closed
as it would be for the case of a profile that has reached the core, and the right represents the rim
of the crystal and is constant and would be a rim in equilibrium with an infinite liquid.


The second Do loop calculates the diffusion coefficient. This needs to be done at every time step
since DFo depends on the Fo concentration which is changing with time.

The third Do loop calculates the profile and contains the finite element formulation for the
diffusion coefficient which depends on concentration (see Short course Appendix I notes).

Plot shows a comparison of natural data and diffusion model.
*****************************************************************************************************************************

ttiter=2800;

Do[bFo=cFo;cFo[[1]]=cFo[[2]];cFo[[ttgrid]]=72.7;Do[DIFFo[[i]]=DFoxd*(10^(3*(0.9-0.01*cFo[[i]]))),{i,1,ttgrid}];Do[cFo[[i]]=bFo[[i]]+R1*(((DIFFo[[i+1]]-DIFFo[[i]])*(bFo[[i+1]]-bFo[[i]]))+(DIFFo[[i]]*(bFo[[i+1]]-2*bFo[[i]]+bFo[[i-1]]))),{i,2,ttgrid-1}],{ttiter}];

eq1=ListPlot[cFo, PlotRange->{Fomin-3,Fomax+3}, PlotStyle ->RGBColor[1,0,1],PlotJoined->True];Show[eq1,datanfop]

***************************************************************************************************************************
The total time in seconds associated with each time one hits shift + return is given below (= timeinsec).
The calculated concentration is exported as dat file (=model1.dat) in your default directory. Remember
that the spacing between each point is equal to δx.
****************************************************************************************************************************

timeinsec=ttiter*δt

timeinyears=timeinsec/(3600*24*365)

Export["model1.dat", cFo]