====== Базовые алгоритмы программы TPF ======

<Code linenums lang-r>

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));
    }
}


</Code>