// SCSSCalculation.sce // written by Don R. Baker, McGill and Sincrotron Trieste, April 2010 // Agrees with old C program 22 April 2010, 2 July 2010 // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // This program is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License // along with this program. If not, see . // To Use: // // To create a datafile in a spreadsheet: // Each analysis is a row, no headers (titles) // The order of the columns is: // T(K), P (GPA), log fO2 ,log fS2, SiO2, TiO2, Al2O3, FeO* ,FeO, Fe2O3, MnO, MgO, CaO, Na2O, K2O, P2O5, H2O (all wt%), S(ppm) , title // After the spreadsheet is creasted save it as a comma delimited file // Then run this program // Start Scilab // Open editor // Load the file into editor: File -> Open // Load file into Scilab: Execute -> Load into Scilab // Define molecular weights MWSIO2= 60.09; MWTIO2 =79.9; MWAL2O3 =101.94; MWFEO =71.85; MWFE2O3 =159.70; MWMNO =70.94; MWMGO =40.32; MWCAO =56.08; MWNA2O =61.982; MWK2O =94.2; MWP2O5 =141.95; MWS =32.066; MWH2O =18.015; //Define constants for calculating SCSS // --------------------------------------------- // Liu et al. model parameters (GCA 2007) for SCSS // ln S = A/T + B + C*P/T + D*ln(MFM) + E* ln X(FeO) + F*((MFM)*X-H2O) + G*ln(X-H2O) ConstTerm_Liu2007=11.35251; //A OneoverT_Liu2007=-0.44546; //B PoverTTerm_Liu2007=-0.03190; //C lnMFMTerm_Liu2007=0.71006; //D lnXFeO_termLiu2007=0.36192; //E MFMxH2Oterm_Liu2007 = -1.98063; //F lnXH2Oterm_Liu2007 = 0.21867; //G // ---------------------------------------------- // Define constants for calculating SCAS // ln S = A/T + B + C*P/T + D*ln(MFM) + E* ln X(CaO) + F*((MFM)*X-H2O) + G*ln(X-H2O) // SCAS(i) = ConstTerm_SCAS + OneoverTTerm_SCAS*10000.0/T(i) + PoverTTerm_SCAS*(P(i)*10000./T(i)) + lnMFMTerm_SCAS*log(MFM(i)) + MFMxH2OTerm_SCAS*((MFM(i))*(MH2O(i)/MolesSumF3(i))) + lnXH2OTerm_SCAS * log(MH2O(i)/MolesSumF3(i))+ lnXCaOTerm_SCAS*log(MCaO(i)/MolesSumF3(i)); ConstTerm_SCAS = 22.53502; OneoverTTerm_SCAS = -1.90738; PoverTTerm_SCAS = 0.40558; lnMFMTerm_SCAS = 0.82637; MFMxH2OTerm_SCAS = -0.79932; lnXH2OTerm_SCAS = 0.81241; lnXCaOTerm_SCAS = -0.21087; // ---------------------------------------------- //********************************************************************** function [ferric_ferrousi]=Calc_FerricFerrous(j,T,P,fO2,MFAL2O3,MFFeOStar,MFCaO,MFNa2O,MFK2O) // Following Kress and Carmichael 1991; Contrib Mineral Petrol 108:82-92. printf("Calculating Ferric iron\n"); a= 0.196; b=1.1492e4; c=-6.675; dAl=-2.243; dFeO=-1.828; dCaO=3.201; dNa2O=5.854; dK2O=6.215; e=-3.36; f=-7.01e-7; g=-1.54e-10; h=3.85e-17; To=1673.; //printf("%d, %f, %f, %e, %f, %f, %f, %f, %f \n",j, T, P, fO2,MFAL2O3, MFFeOStar, MFCaO,MFNa2O, MFK2O ); ferric_ferrousi = a*log(fO2) + b/T + c + dAl*MFAL2O3+dFeO*MFFeOStar+dCaO*MFCaO+dNa2O*MFNa2O+dK2O*MFK2O + e*(1.0-To/T - log(T/To)) +f*(P*(10.0^9.0))/T + g*((T-To)*P*(10.0^9.0))/T +h* ((P*(10.0^9.0))*(P*(10.0^9.0)))/T ferric_ferrousi = exp(ferric_ferrousi); printf ("Fe2O3/FeO ratio %f\n",ferric_ferrousi); endfunction; //********************************************************** function [MFMi] = Calc_MFM(MNa2O,MK2O,MCaO,MMgO,MFeO,MFe2O3,MSiO2,MAL2O3,MTiO2,MMnO,MP2O5,MH2O) MNa2O=2.0*MNa2O; MK2O= 2.0*MK2O; MFe2O3 = 2.0*MFe2O3; MAL2O3 = 2.0*MAL2O3; MP2O5 = 2.0*MP2O5; //calculate MFM dry, no H2O used CationSum = MNa2O+MK2O+MCaO+MMgO+MFeO+MFe2O3+MSiO2+MAL2O3+MTiO2+MMnO+MP2O5; MSiO2 = MSiO2/CationSum; MTiO2 = MTiO2/CationSum; MAL2O3 = MAL2O3/CationSum; MFe2O3 = MFe2O3/CationSum; MFeO = MFeO/CationSum; MMnO = MMnO/CationSum; MMgO = MMgO/CationSum; MCaO = MCaO/CationSum; MNa2O = MNa2O/CationSum; MK2O = MK2O/CationSum; MH2O = MH2O/CationSum; MFMi = (MNa2O+MK2O+2.0*(MCaO+MMgO+MFeO))/(MSiO2*(MAL2O3+MFe2O3)); endfunction; // *************************************************************** // *************************************************************** */ function [fS2calc]=Calc_fS2(T,P,fO2,XFeO) DeltaV = 0.4517; DeltafO2 = (log(fO2)/log(10.)) - (-24014./T + 8.555 +0.092*(P*10000.-1.)/T); //O'Neill 1987 FMQ //printf ("log FMQ fO2 %f\n",-24014./T + 8.555); printf("XFeO %f \n", XFeO); fS2calc = 6.7 - 12800./T - 2.* log(XFeO)/log(10.) + DeltafO2; printf("log fO2 %f, log fS2 at 1 bar %f\n", log(fO2)/log(10.), fS2calc); fS2calc = log((10.^fS2calc)) + (2.*DeltaV*(P*10000.-1))/(2.303*8.314*T); endfunction; // *************************************************************** */ // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! // Main Program begins here! //********************************************************************* // Read in data from a comma-delimited file // Input data file format: // T(K), P (GPA), log fO2 ,log fS2, SiO2, TiO2, Al2O3, FeO* ,FeO, Fe2O3, MnO, MgO, CaO, Na2O, K2O, P2O5, H2O (all wt%), S(ppm) , title datafilename =input("What is the complete data file name?","string"); //datafilename='1050_300MPa_FMQ.txt'; // temporary file name for testing ONLY u=mopen(datafilename,'r'); [n,TK,GPa,logfO2,logfS2,SiO2,TiO2,Al2O3,FeOstar,Fe2O3,FeO,MnO,MgO,CaO,Na2O,K2O,P2O5,H2O,S,Title]=mfscanf(-1,u,'%f,%f, %f,%f, %f,%f, %f,%f,%f,%f,%f,%f, %f,%f,%f, %f,%f,%f,%s'); mclose(u); // ********************************************************************* // All data has now been read into Scilab TotalAnalyses=length(TK); for (i=1:TotalAnalyses) T(i)=TK(i); P(i)=GPa(i); fO2(i)= (10.0^logfO2(i)); if logfS2(i) <> 0.0 fS2(i) = log((10.0^logfS2(i))); else fS2(i) = logfS2(i); end; WSiO2(i)=SiO2(i); WTiO2(i)=TiO2(i); WAL2O3(i)=Al2O3(i); WFeOStar(i)=FeOstar(i); WFeO(i) = FeO(i); WFe2O3(i)=Fe2O3(i); WMnO(i)=MnO(i); WMgO(i)=MgO(i); WCaO(i)=CaO(i); WNa2O(i)=Na2O(i); WK2O(i)=K2O(i); WP2O5(i)=P2O5(i); WS(i)=S(i)/1e4; //convert ppm to weight percent WH2O(i)=H2O(i); // calculate mols in 100 g of rock (i.e., wt% input) MSiO2(i)=WSiO2(i)/MWSIO2; MTiO2(i)=WTiO2(i)/MWTIO2; MAL2O3(i)=WAL2O3(i)/MWAL2O3; MFeOStar(i)=WFeOStar(i)/MWFEO; MFeO(i) = WFeO(i)/MWFEO; MFe2O3(i) = WFe2O3(i)/MWFE2O3; MMnO(i)=WMnO(i)/MWMNO; MMgO(i)=WMgO(i)/MWMGO; MCaO(i)=WCaO(i)/MWCAO; MNa2O(i)=WNa2O(i)/MWNA2O; MK2O(i)=WK2O(i)/MWK2O; MP2O5(i)=WP2O5(i)/MWP2O5; MS(i)=WS(i)/MWS; MH2O(i)=WH2O(i)/MWH2O; if (MH2O(i) == 0.0) MH2O(i) = 0.001; end; printf(" sample %d\n, Temp (K) %f, P(GPa) %f, fO2 %e\n,wt percent\n SiO2 %f\n TiO2 %f\n Al2O3 %f\n FeO* %f\n MnO %f\n MgO %f\n CaO %f\n Na2O %f\n K2O %f\n P2O5 %f\n S %f\n",i, T(i),P(i),fO2(i),WSiO2(i), WTiO2(i),WAL2O3(i),WFeOStar(i),WMnO(i),WMgO(i),WCaO(i),WNa2O(i),WK2O(i),WP2O5(i), WS(i)); MolesSum(i) = MSiO2(i)+MTiO2(i)+MAL2O3(i)+MFeOStar(i)+MMnO(i)+MMgO(i)+MCaO(i)+MNa2O(i)+MK2O(i)+MP2O5(i)+MH2O(i)+MS(i); printf(" sample %d\n, Temp (K) %f, P(GPa) %f, fO2 %e\n,moles \n SiO2 %f\n TiO2 %f\n Al2O3 %f\n FeO* %f\n MnO %f\n MgO %f\n CaO %f\n Na2O %f\n K2O %f\n P2O5 %f\n S %f\n",i, T(i),P(i),fO2(i),MSiO2(i), MTiO2(i),MAL2O3(i),MFeOStar(i),MMnO(i),MMgO(i),MCaO(i),MNa2O(i),MK2O(i),MP2O5(i), MS(i)); if (WFeOStar(i) <> 0.0) [ferric_ferrousi]=Calc_FerricFerrous(i, T(i),P(i),fO2(i),MAL2O3(i)/MolesSum(i),MFeOStar(i)/MolesSum(i),MCaO(i)/MolesSum(i),MNa2O(i)/MolesSum(i),MK2O(i)/MolesSum(i)); ferric_ferrous(i)=ferric_ferrousi; XFeO = (MFeOStar(i)/MolesSum(i))/(1.0 + 2.0*(ferric_ferrousi)); XFe2O3 = XFeO*ferric_ferrousi; MFeO(i) = MolesSum(i)*XFeO; MFe2O3(i)=MolesSum(i)*XFe2O3; WFeO(i) = MFeO(i)*MWFEO; WFe2O3(i) = MFe2O3(i)*MWFE2O3; //printf("Moles FeO* %f, Moles FeO %f, Moles Fe2O3 % lf\n",MFeOStar(j),MFeO(j),MFe2O3(j));MFe2O3(j)); end; MolesSumF3(i) = MSiO2(i)+MTiO2(i)+MAL2O3(i)+MFe2O3(i)+MFeO(i)+MMnO(i)+MMgO(i)+MCaO(i)+MNa2O(i)+MK2O(i)+MP2O5(i)+MH2O(i)+MS(i); [MFMi] = Calc_MFM(MNa2O(i),MK2O(i),MCaO(i),MMgO(i),MFeO(i),MFe2O3(i),MSiO2(i),MAL2O3(i), MTiO2(i), MMnO(i),MP2O5(i),MH2O(i)); MFM(i)=MFMi; printf("MFM(%d) = %f\n",i,MFM(i)); if (fS2(i) == 0.0) [fS2calc] = Calc_fS2(T(i),P(i),fO2(i),MFeO(i)/MolesSumF3(i)); fS2(i) = fS2calc; // fS2(i) is really ln(fS2(i)) end; // If fO2 <= NNO+1.5 calculate SCSS, if fO2 >= NNO + 1.5 calculate SCAS // Calculate NNO fO2 following Huebner and Sato (1971?) log_10NNO = 9.26-24930./T(i)+0.046*(10000.*P(i)-1.)/T(i); log_10NNO1_5= log_10NNO+1.5; MFeOStar = MFeO(i)+2*MFe2O3(i); [ferric_ferrousi]=Calc_FerricFerrous(i, T(i),P(i),10.0^log_10NNO1_5,MAL2O3(i)/MolesSum(i),MFeOStar/MolesSum(i),MCaO(i)/MolesSum(i),MNa2O(i)/MolesSum(i),MK2O(i)/MolesSum(i)); SCSS_Liu2007(i) = 0.0; SCAS(i) = 0.0; if (ferric_ferrous(i)<= ferric_ferrousi) or (fO2(i)<=10.0^(log_10NNO1_5)) // --------------------------------------------------------------------- // Now calculate the ln(SCSS) in ppm following Liu et al (CGA 2007)--------------- SCSS_Liu2007(i) = OneoverT_Liu2007*10000.0/T(i) + ConstTerm_Liu2007 + lnMFMTerm_Liu2007*log(MFM(i))+PoverTTerm_Liu2007*(P(i)*10000./T(i))+ lnXFeO_termLiu2007*log(MFeO(i)/MolesSumF3(i))+MFMxH2Oterm_Liu2007*((MFM(i))*(MH2O(i)/MolesSumF3(i))) + lnXH2Oterm_Liu2007 * log (MH2O(i)/MolesSumF3(i)); // printf("The SCSS calculated following Liu et al.(2007) is %f\n\n\n\n",SCSS_Liu2007(i)); printf("ln S(ppm) meas %f S(ppm) meas %f ln SCSS(ppm) Liu calc %f S(ppm) Liu calc %f %f %s \n", log(1e4*WS(i)), 1e4*WS(i), SCSS_Liu2007(i), exp(SCSS_Liu2007(i) ), T(i), Title(i)); // End calculation of SCSS following Liu et al (GCA 2007)-------------- // --------------------------------------------------------------------- end; if (ferric_ferrous(i)>= ferric_ferrousi) or (fO2(i)>=(10.0^log_10NNO1_5)) // --------------------------------------------------------------------- // Now calculate the ln(SCAS) in ppm following Moretti and Baker (2011) SCAS(i) = ConstTerm_SCAS + OneoverTTerm_SCAS*10000.0/T(i) + PoverTTerm_SCAS*10000.0*P(i)/T(i) + lnMFMTerm_SCAS*log(MFM(i)) + MFMxH2OTerm_SCAS*(MFM(i)*(MH2O(i)/MolesSumF3(i))) + lnXH2OTerm_SCAS*log(MH2O(i)/MolesSumF3(i))+ lnXCaOTerm_SCAS*log(MCaO(i)/MolesSumF3(i)); // printf("The SCAS calculated following Moretti & Baker.(2011) is %f\n\n\n\n",SCAS(i)); printf("ln S(ppm) meas %f S(ppm) meas %f ln SCAS(ppm) %f SCAS(ppm) %f %f %s \n\n\n\n", log(1e4*WS(i)), 1e4*WS(i), SCAS(i), exp(SCAS(i)), T(i), Title(i)); // End calculation of SCAS // --------------------------------------------------------------------- end; end; // end looping through all analyses //********************************************************************* // Print out data from a tab-delimited file zzz=length(datafilename)-4; name=strsplit(datafilename,[zzz]); datafilename=name(1)+'_SCSS_SCAS_Calculations.txt'; u=file('open',datafilename,'unknown'); fprintf(u, "#Calculation of the SCSS and SCAS in Scliab:\n"); fprintf(u,"#ln S = A/T + B + C*P/T + D*ln(MFM) + E* ln X(FeO) + F*((MFM)*X-H2O) + G*ln(X-H2O)\n"); fprintf(u,"#CAUTION: A value of 0 or 1 in a column means that the value was not calculated\n"); fprintf(u,"#ln S(ppm) meas S(ppm) meas ln SCSS(ppm) Liu calc SCSS(ppm) Liu calc ln SCAS(ppm) SCAS(ppm) H2O P(GPa) T(K) Title\n"); for (i=1:TotalAnalyses) fprintf(u,"%f %f %f %f %f %f %f %f %f %s\n", log(1e4*WS(i)), 1e4*WS(i), SCSS_Liu2007(i), exp(SCSS_Liu2007(i)),SCAS(i), exp(SCAS(i)),WH2O(i), P(i), T(i), Title(i) ); end; file('close',u); // *********************************************************************