source('utilities.r')
tpf.elements <- read.csv("constants_elements_tpf.csv",stringsAsFactors=F)
tpf.oxides_ord <- c('SiO2_wtp', 'Al2O3_wtp', 'TiO2_wtp', 'Cr2O3_wtp', 'Fe2O3_wtp',
'FeO_wtp', 'MnO_wtp', 'MgO_wtp', 'CaO_wtp', 'Na2O_wtp', 'K2O_wtp',
'BaO_wtp', 'NiO_wtp', 'ZnO_wtp', 'BeO_wtp', 'P2O5_wtp', 'V2O5_wtp',
'CoO_wtp', 'GeO2_wtp', 'ZrO2_wtp', 'Cs2O_wtp', 'PbO_wtp', 'SO3_wtp',
'CO2_wtp', 'F_wtp', 'Cl_wtp');
tpf.oxides_elements_ord <- c('Si', 'Al', 'Ti', 'Cr', 'Fe3', 'Fe2', 'Mn', 'Mg',
'Ca', 'Na', 'K', 'Ba', 'Ni', 'Zn', 'Be', 'P', 'V',
'Co', 'Ge', 'Zr', 'Cs', 'Pb', 'S', 'C', 'F', 'Cl');
tpf.elements_ord <- c('Si', 'Al', 'Ti', 'Cr', 'Fe3', 'Fe2', 'Mn', 'Mg', 'Ca',
'Na', 'K', 'Ba', 'Ni', 'Zn', 'Be', 'P', 'V', 'Co', 'Ge',
'Zr', 'Cs', 'Pb', 'S', 'C', 'F', 'Cl', 'Al4', 'Al6', 'Fe');
tpf.elements_rename <- c('Si_pfu', 'Al_pfu', 'Ti_pfu', 'Cr_pfu', 'Fe_p3_pfu', 'Fe_p2_pfu',
'Mn_pfu', 'Mg_pfu', 'Ca_pfu', 'Na_pfu', 'K_pfu', 'Ba_pfu', 'Ni_pfu',
'Zn_pfu', 'Be_pfu', 'P_pfu', 'V_pfu', 'Co_pfu', 'Ge_pfu', 'Zr_pfu',
'Cs_pfu', 'Pb_pfu', 'S_pfu', 'C_pfu', 'F_pfu', 'Cl_pfu',
'Al_c4_pfu', 'Al_c6_pfu', 'Fe_pfu');
tpf.rename <- function(wtp_data){
renamed <- selectNames(wtp_data,tpf.oxides_ord);
names(renamed) <- tpf.oxides_elements_ord;
return(renamed);
}
tpf.unrename <- function(mw, labels){
renamed <- selectNames(mw,tpf.elements_ord);
names(renamed) <- tpf.elements_rename;
return(cbind(Name = labels, renamed[,colSums(renamed^2) !=0]));
}
tpf.select_cations <- function(start,end){
subs <- tpf.oxides_elements_ord[match(start,tpf.elements_ord):match(end,tpf.elements_ord)]
return(subs);
}
tpf.norm <- function(mw,sk){
ccat <- tpf.select_cations('Si','C');
an <- tpf.select_cations('F','Cl');
cat_ws <- subset(tpf.elements,element %in% ccat)$weight;
cat_cnt <- subset(tpf.elements,element %in% ccat)$cat_count;
an_ws <- c(1,1);
ws <- c(cat_ws/sk/cat_cnt,an_ws);
ret <- rowApply('/',mw,ws);
return(ret);
}
tpf.normfcl <- function(mw,skfcl){
ccat <- tpf.select_cations('Si','C');
an <- tpf.select_cations('F','Cl');
cat_ws <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1);
an_ws <- subset(tpf.elements,element %in% an)$weight;
ws <- c(cat_ws,an_ws/skfcl);
ret <- rowApply('/',mw,ws);
return(ret);
}
tpf.normbt <- function(mw,sk){
ccat1 <- c(tpf.select_cations('Si','Ca'));
ccat2 <- c(tpf.select_cations('Ba','C'));
an <- tpf.select_cations('F','Cl');
cat1_ws <- subset(tpf.elements,element %in% ccat1)$weight;
cat1_cnt <- subset(tpf.elements,element %in% ccat1)$cat_count;
catt_ws <- c(1,1);
cat2_ws <- subset(tpf.elements,element %in% ccat2)$weight;
cat2_cnt <- subset(tpf.elements,element %in% ccat2)$cat_count;
an_ws <- c(1,1);
ws <- c(cat1_ws/sk/cat1_cnt,catt_ws,cat2_ws/sk/cat2_cnt,an_ws);
ret <- rowApply('/',mw,ws);
ret$Na <- 0.0;
ret$K <- 0.0;
return(ret);
}
tpf.sumfcl <- function(mw,so,theta = 0.000001){
cat <- c(tpf.select_cations('Si','Cl'));
cat_ws <- subset(tpf.elements,element %in% cat)$weight;
cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
cat_ox[25] <- 1;
cat_ox[26] <- 1;
ds <- rowApply('*',mw,cat_ox/cat_ws);
s <-rowSums(ds);
sfc <- (ds$F + ds$Cl);
s <- s - sfc/2;
#return (ifelse (s > theta, so/s, 0.0));
return (so/s);
}
tpf.sum <- function(mw,so,theta = 0.000001){
cat <- c(tpf.select_cations('Si','C'));
cat_ws <- subset(tpf.elements,element %in% cat)$weight;
cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
ds <- rowApply('*',mw,cat_ox/cat_ws);
s <-rowSums(ds);
#return (ifelse (s > theta, so/s, 0.0));
return (so/s);
}
tpf.alal <- function(mw,siv){
mw$Al4 <- ifelse(mw$Si+mw$Al>siv, ifelse(mw$Si<siv,siv-mw$Si,0),mw$Al);
mw$Al6 <- mw$Al - mw$Al4;
return(mw);
}
tpf.omesod_v <- function(mw,so,siv){
rn <- tpf.rename(mw);
clkm <- NULL;
for(i in 1:nrow(rn))
{
clk <- tpf.omesod(rn[i,],so,siv);
if (i == 1) { clkm = clk; }
else { clkm = rbind(clkm,clk); }
}
return(tpf.unrename(clkm,mw$Name));
}
tpf.omesod <- function(mw,so,siv){
sk <- tpf.sum(mw,so);
skfcl <- tpf.sumfcl(mw,so);
mw1 <- tpf.norm(mw,sk);
mw2 <- tpf.normfcl(mw1,skfcl);
mw3 <- tpf.alal(mw2,siv);
mw3$Fe = mw3$Fe3 + mw3$Fe2;
return (mw3);
}
tpf.sumo <- function (mw){
s <- 0;
sfc <- 0;
cat <- c(tpf.select_cations('Si','C'));
cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
cat_cat <- subset(tpf.elements,element %in% cat)$cat_count;
clist <- c(c(cat_ox/cat_cat),c(1,1));
s <- rowSums(rowApply('*',mw,clist));
sfc <- 0.5*(mw$F+mw$Cl);
return (s-sfc);
}
tpf.fefe <- function(mw,ll,mm,theta = 0.000001){
#Non Vector Operation! Tests requried!
os = tpf.sumo(mw) + (mw$F + mw$Cl)/2;
if (os>theta) {
if( os<ll && mw$Fe > (mm-2*os)) {
mw$Fe3 <- mm-2*os;
mw$Fe2 <- mw$Fe - mw$Fe2;
}else{
mw <- rowApply('*',mw,ll/os);
}
}
return(mw);
}
tpf.recalc_minaral <- function(mn, mw){
rn <- tpf.rename(mw);
clkm <- NULL;
for(i in 1:nrow(rn))
{
clk <- tpf.recalc_minaral_r(mn, rn[i,]);
if (i == 1) { clkm = clk; }
else { clkm = rbind(clkm,clk); }
}
ds <- tpf.unrename(clkm,mw$Name)
return(cbind(data.frame(Name=mw$Name),selectNames(ds, tpf.elements_rename)));
}
tpf.recalc_minaral_r <- function(mn, mw) {
if (mn %in% c("SCP")) {
mw <- tpf.omesod(mw,33,999);
sct <- mw$Si+mw$Al;
kp <- ifelse(sct != 0, 12/sct, 0);
mw <- rowApply('*',mw,kp);
mw <- tpf.alal(mw, 12);
mw$Fe <- mw$Fe3 + mw$Fe2;
return(mw);
}
if (mn %in% c("SPR")) {
mw <- tpf.omesod(mw,10,999);
s <- rowSums(mw) - mw$F - mw$Cl;
s <- ifelse(s != 0, 7/s, 0);
mw <- rowApply('*',mw,s);
mw$Fe <- mw$Fe3 + mw$Fe2;
mw$Fe2 <- mw$Fe;
mw$Fe3 <- 0;
mw <- tpf.fefe(mw, 10,20);
mw <- tpf.alal(mw, 3);
return(mw);
}
if (mn %in% c("ILM", "HEM")) {
mw <- tpf.omesod(mw,3,999);
s <- rowSums(mw) - mw$F - mw$Cl;
s <- ifelse(s != 0, 2/s, 0);
mw <- rowApply('*',mw,s);
mw$Al4 <- mw$Al;
mw$Al6 <- 0;
mw$Fe <- mw$Fe3+mw$Fe2;
mw$Fe2 <- mw$Fe;
mw$Fe3 <- 0;
mw <- tpf.fefe(mw, 3,6);
return(mw);
}
if (mn %in% c("MAG", "SPL")) {
mw <- tpf.omesod(mw,4,999);
s <- rowSums(mw) - mw$F - mw$Cl;
s <- ifelse(s != 0, 3/s, 0);
mw <- rowApply('*',mw,s);
mw$Al4 <- mw$Al;
mw$Al6 <- 0;
mw$Fe <- mw$Fe3+mw$Fe2;
mw$Fe2 <- mw$Fe;
mw$Fe3 <- 0;
mw <- tpf.fefe(mw, 4,8);
return(mw);
}
if (mn %in% c("OPX", "PGT", "CPX")) {
return(tpf.omesod(mw,6,2));
}
if (mn %in% c("GRT")) {
return(tpf.omesod(mw,12,999));
}
if (mn %in% c("EP", "ZO")) {
return(tpf.omesod(mw,12.5,999));
}
if (mn %in% c("CRD")) {
return(tpf.omesod(mw,18,6));
}
if (mn %in% c("CAM", "HBL", "CUM", "ATH")) {
return(tpf.omesod(mw,23,8));
}
if (mn %in% c("KFS", "PL")) {
return(tpf.omesod(mw,8,999));
}
if (mn %in% c("BT", "MS")) {
return(tpf.omesod(mw,11,4));
}
if (mn %in% c("OL", "MTC")) {
return(tpf.omesod(mw,4,1));
}
if (mn %in% c("RT", "QTZ")) {
return(tpf.omesod(mw,2,999));
}
if (mn %in% c("KY","AND","SIL","TTN")) {
return(tpf.omesod(mw,5,999));
}
if (mn %in% c("L")) {
return(tpf.omesod(mw,50,999));
}
if (mn %in% c("SID", "DOL", "MGS", "CAL", "FE", "WUS")) {
return(tpf.omesod(mw,1,999));
}
if (mn %in% c("WO")) {
return(tpf.omesod(mw,3,999));
}
if (mn %in% c("NE")) {
return(tpf.omesod(mw,4,2));
}
if (mn %in% c("ST")) {
return(tpf.omesod(mw,46,999));
}
if (mn %in% c("TLC")) {
return(tpf.omesod(mw,11,999));
}
if (mn %in% c("CHL")) {
return(tpf.omesod(mw,14,4));
}
if (mn %in% c("CHD")) {
return(tpf.omesod(mw,12,999));
}
if (mn %in% c("OSU")) {
return(tpf.omesod(mw,30,999));
}
if (mn %in% c("KRN")) {
return(tpf.omesod(mw,43,999));
}
if (mn %in% c("STP")) {
return(tpf.omesod(mw,23,999));
}
if (mn %in% c("HU")) {
return(tpf.omesod(mw,13,999));
}
if (mn %in% c("CRN")) {
return(tpf.omesod(mw,3,999));
}
}