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

  1.  
  2. source('utilities.r')
  3.  
  4. tpf.elements <- read.csv("constants_elements_tpf.csv",stringsAsFactors=F)
  5.  
  6. tpf.oxides_ord <- c('SiO2_wtp', 'Al2O3_wtp', 'TiO2_wtp', 'Cr2O3_wtp', 'Fe2O3_wtp',
  7.                     'FeO_wtp', 'MnO_wtp', 'MgO_wtp', 'CaO_wtp', 'Na2O_wtp', 'K2O_wtp',
  8.                     'BaO_wtp', 'NiO_wtp', 'ZnO_wtp', 'BeO_wtp', 'P2O5_wtp', 'V2O5_wtp',
  9.                     'CoO_wtp', 'GeO2_wtp', 'ZrO2_wtp', 'Cs2O_wtp', 'PbO_wtp', 'SO3_wtp',
  10.                     'CO2_wtp', 'F_wtp', 'Cl_wtp');
  11.  
  12. tpf.oxides_elements_ord <- c('Si', 'Al', 'Ti', 'Cr', 'Fe3', 'Fe2', 'Mn', 'Mg',
  13.                              'Ca', 'Na', 'K', 'Ba', 'Ni', 'Zn', 'Be', 'P', 'V',
  14.                              'Co', 'Ge', 'Zr', 'Cs', 'Pb', 'S', 'C', 'F', 'Cl');
  15.  
  16. tpf.elements_ord <- c('Si', 'Al', 'Ti', 'Cr', 'Fe3', 'Fe2', 'Mn', 'Mg', 'Ca',
  17.                       'Na', 'K', 'Ba', 'Ni', 'Zn', 'Be', 'P', 'V', 'Co', 'Ge',
  18.                       'Zr', 'Cs', 'Pb', 'S', 'C', 'F', 'Cl', 'Al4', 'Al6', 'Fe');
  19.  
  20. tpf.elements_rename <- c('Si_pfu', 'Al_pfu', 'Ti_pfu', 'Cr_pfu', 'Fe_p3_pfu', 'Fe_p2_pfu',
  21.                          'Mn_pfu', 'Mg_pfu', 'Ca_pfu', 'Na_pfu', 'K_pfu', 'Ba_pfu', 'Ni_pfu',
  22.                          'Zn_pfu', 'Be_pfu', 'P_pfu', 'V_pfu', 'Co_pfu', 'Ge_pfu', 'Zr_pfu',
  23.                          'Cs_pfu', 'Pb_pfu', 'S_pfu', 'C_pfu', 'F_pfu', 'Cl_pfu',
  24.                          'Al_c4_pfu', 'Al_c6_pfu', 'Fe_pfu');
  25.  
  26.  
  27. tpf.rename <- function(wtp_data){
  28. renamed <- selectNames(wtp_data,tpf.oxides_ord);
  29. names(renamed) <- tpf.oxides_elements_ord;
  30. return(renamed);
  31. }
  32.  
  33. tpf.unrename <- function(mw, labels){
  34. renamed <- selectNames(mw,tpf.elements_ord);
  35. names(renamed) <- tpf.elements_rename;
  36. return(cbind(Name = labels, renamed[,colSums(renamed^2) !=0]));
  37. }
  38.  
  39. tpf.select_cations <- function(start,end){
  40. subs <- tpf.oxides_elements_ord[match(start,tpf.elements_ord):match(end,tpf.elements_ord)]
  41. return(subs);
  42. }
  43.  
  44. tpf.norm <- function(mw,sk){
  45.     ccat <- tpf.select_cations('Si','C');
  46.     an <- tpf.select_cations('F','Cl');
  47.     cat_ws <- subset(tpf.elements,element %in% ccat)$weight;
  48.     cat_cnt <- subset(tpf.elements,element %in% ccat)$cat_count;
  49.     an_ws <- c(1,1);
  50.     ws <- c(cat_ws/sk/cat_cnt,an_ws);
  51.     ret <- rowApply('/',mw,ws);
  52.     return(ret);
  53. }
  54.  
  55. tpf.normfcl <- function(mw,skfcl){
  56.     ccat <- tpf.select_cations('Si','C');
  57.     an <- tpf.select_cations('F','Cl');
  58.     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); 
  59.     an_ws <- subset(tpf.elements,element %in% an)$weight;
  60.     ws <- c(cat_ws,an_ws/skfcl);
  61.     ret <- rowApply('/',mw,ws);
  62.     return(ret);
  63. }
  64.  
  65. tpf.normbt <- function(mw,sk){
  66.     ccat1 <- c(tpf.select_cations('Si','Ca'));
  67.     ccat2 <- c(tpf.select_cations('Ba','C'));
  68.     an <- tpf.select_cations('F','Cl');
  69.     cat1_ws <- subset(tpf.elements,element %in% ccat1)$weight;
  70.     cat1_cnt <- subset(tpf.elements,element %in% ccat1)$cat_count;
  71.     catt_ws <- c(1,1); 
  72.     cat2_ws <- subset(tpf.elements,element %in% ccat2)$weight;
  73.     cat2_cnt <- subset(tpf.elements,element %in% ccat2)$cat_count;
  74.     an_ws <- c(1,1);
  75.     ws <- c(cat1_ws/sk/cat1_cnt,catt_ws,cat2_ws/sk/cat2_cnt,an_ws);
  76.     ret <- rowApply('/',mw,ws);
  77.     ret$Na <- 0.0;
  78.     ret$K <- 0.0;
  79.     return(ret);
  80. }
  81.  
  82. tpf.sumfcl <- function(mw,so,theta = 0.000001){
  83. cat <- c(tpf.select_cations('Si','Cl'));
  84. cat_ws <- subset(tpf.elements,element %in% cat)$weight;
  85. cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
  86. cat_ox[25] <- 1;
  87. cat_ox[26] <- 1;
  88.  
  89. ds <- rowApply('*',mw,cat_ox/cat_ws);
  90.  
  91. s <-rowSums(ds);
  92. sfc <- (ds$F + ds$Cl);
  93. s <- s - sfc/2;
  94.  
  95. #return (ifelse (s > theta, so/s, 0.0));
  96. return (so/s);
  97. }
  98.  
  99.  
  100. tpf.sum <- function(mw,so,theta = 0.000001){
  101. cat <- c(tpf.select_cations('Si','C'));
  102. cat_ws <- subset(tpf.elements,element %in% cat)$weight;
  103. cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
  104.  
  105. ds <- rowApply('*',mw,cat_ox/cat_ws);
  106.  
  107. s <-rowSums(ds);
  108.  
  109. #return (ifelse (s > theta, so/s, 0.0));
  110. return (so/s);
  111. }
  112.  
  113.  
  114. tpf.alal <- function(mw,siv){
  115. mw$Al4 <- ifelse(mw$Si+mw$Al>siv, ifelse(mw$Si<siv,siv-mw$Si,0),mw$Al);
  116. mw$Al6 <- mw$Al - mw$Al4;
  117. return(mw);
  118. }
  119.  
  120.  
  121. tpf.omesod_v <- function(mw,so,siv){
  122. rn <- tpf.rename(mw);
  123. clkm <- NULL;
  124. for(i in 1:nrow(rn))
  125. {
  126. clk <- tpf.omesod(rn[i,],so,siv);
  127. if (i == 1) { clkm = clk; }
  128. else { clkm = rbind(clkm,clk); }
  129. }
  130. return(tpf.unrename(clkm,mw$Name));
  131. }
  132.  
  133. tpf.omesod <- function(mw,so,siv){
  134. sk <- tpf.sum(mw,so);
  135. skfcl <- tpf.sumfcl(mw,so);
  136. mw1 <- tpf.norm(mw,sk);
  137. mw2 <- tpf.normfcl(mw1,skfcl);
  138. mw3 <- tpf.alal(mw2,siv);
  139. mw3$Fe = mw3$Fe3 + mw3$Fe2;
  140. return (mw3);
  141. }
  142.  
  143.  
  144. tpf.sumo <- function (mw){
  145.     s <- 0;
  146.     sfc <- 0;
  147.     
  148.     cat <- c(tpf.select_cations('Si','C'));
  149.     cat_ox <- subset(tpf.elements,element %in% cat)$an_count;
  150.     cat_cat <- subset(tpf.elements,element %in% cat)$cat_count;
  151.  
  152.     clist <-  c(c(cat_ox/cat_cat),c(1,1));
  153.  
  154.     s <- rowSums(rowApply('*',mw,clist));
  155.     sfc <- 0.5*(mw$F+mw$Cl);
  156.  
  157.     return (s-sfc);
  158. }
  159.  
  160.  
  161. tpf.fefe <- function(mw,ll,mm,theta = 0.000001){
  162.     #Non Vector Operation! Tests requried!
  163.     
  164.     os = tpf.sumo(mw) + (mw$F + mw$Cl)/2;
  165.  
  166.     if (os>theta) {
  167.         if( os<ll && mw$Fe > (mm-2*os)) {
  168.             mw$Fe3 <- mm-2*os;
  169.             mw$Fe2 <- mw$Fe - mw$Fe2;
  170.         }else{
  171.             mw <- rowApply('*',mw,ll/os);
  172.         }
  173.     }
  174.     return(mw);
  175. }
  176.  
  177.  
  178. tpf.recalc_minaral <- function(mn, mw){
  179. rn <- tpf.rename(mw);
  180. clkm <- NULL;
  181. for(i in 1:nrow(rn))
  182. {
  183. clk <- tpf.recalc_minaral_r(mn, rn[i,]);
  184. if (i == 1) { clkm = clk; }
  185. else { clkm = rbind(clkm,clk); }
  186. }
  187.  
  188. ds <- tpf.unrename(clkm,mw$Name)
  189. return(cbind(data.frame(Name=mw$Name),selectNames(ds, tpf.elements_rename)));
  190. }
  191.  
  192.  
  193.  
  194. tpf.recalc_minaral_r <- function(mn, mw) {
  195.  
  196. if (mn %in% c("SCP")) {  
  197.         mw <- tpf.omesod(mw,33,999);
  198.         sct <- mw$Si+mw$Al;
  199.         kp <- ifelse(sct != 0, 12/sct, 0);
  200.         mw <- rowApply('*',mw,kp);
  201.         mw <- tpf.alal(mw, 12);
  202.         mw$Fe <- mw$Fe3 + mw$Fe2;
  203.         return(mw);
  204. }
  205. if (mn %in% c("SPR")) {  
  206.          mw <- tpf.omesod(mw,10,999);
  207.          s <- rowSums(mw) - mw$F - mw$Cl;
  208.          s <- ifelse(s != 0, 7/s, 0);
  209.          mw <- rowApply('*',mw,s);
  210.          mw$Fe <- mw$Fe3 + mw$Fe2;
  211.          mw$Fe2 <- mw$Fe;
  212.          mw$Fe3 <- 0;
  213.          mw <- tpf.fefe(mw, 10,20);
  214.          mw <- tpf.alal(mw, 3);
  215.          return(mw);
  216. }
  217.     
  218. if (mn %in% c("ILM", "HEM")) {  
  219.          mw <- tpf.omesod(mw,3,999);
  220.          s <- rowSums(mw) - mw$F - mw$Cl;
  221.          s <- ifelse(s != 0, 2/s, 0);
  222.          mw <- rowApply('*',mw,s);
  223.          mw$Al4 <- mw$Al;
  224.          mw$Al6 <- 0;
  225.          mw$Fe <- mw$Fe3+mw$Fe2;
  226.          mw$Fe2 <- mw$Fe;
  227.          mw$Fe3 <- 0;
  228.          mw <- tpf.fefe(mw, 3,6);
  229.          return(mw);
  230. }
  231.  
  232. if (mn %in% c("MAG", "SPL")) {   
  233.          mw <- tpf.omesod(mw,4,999);
  234.          s <- rowSums(mw) - mw$F - mw$Cl;
  235.          s <- ifelse(s != 0, 3/s, 0);
  236.          mw <- rowApply('*',mw,s);
  237.          mw$Al4 <- mw$Al;
  238.          mw$Al6 <- 0;
  239.          mw$Fe <- mw$Fe3+mw$Fe2;
  240.          mw$Fe2 <- mw$Fe;
  241.          mw$Fe3 <- 0;
  242.          mw <- tpf.fefe(mw, 4,8);
  243.          return(mw);
  244. }
  245.     if (mn %in% c("OPX", "PGT", "CPX")) {  
  246.         return(tpf.omesod(mw,6,2));
  247.     }
  248.     if (mn %in% c("GRT")) {  
  249.         return(tpf.omesod(mw,12,999));
  250.     }
  251.     if (mn %in% c("EP", "ZO")) {        
  252.         return(tpf.omesod(mw,12.5,999));
  253.     }
  254.     if (mn %in% c("CRD")) {  
  255.         return(tpf.omesod(mw,18,6));
  256.     }
  257.     if (mn %in% c("CAM", "HBL", "CUM", "ATH")) {        
  258.         return(tpf.omesod(mw,23,8));
  259.     }
  260.     if (mn %in% c("KFS", "PL")) {  
  261.         return(tpf.omesod(mw,8,999));
  262.     }
  263.     if (mn %in% c("BT", "MS")) {  
  264.         return(tpf.omesod(mw,11,4));
  265.     }
  266.     if (mn %in% c("OL", "MTC")) {  
  267.         return(tpf.omesod(mw,4,1));
  268.     }
  269.     if (mn %in% c("RT", "QTZ")) {        
  270.         return(tpf.omesod(mw,2,999));
  271.     }
  272.     if (mn %in% c("KY","AND","SIL","TTN")) {        
  273.         return(tpf.omesod(mw,5,999));
  274.     }
  275.     if (mn %in% c("L")) {  
  276.        return(tpf.omesod(mw,50,999));
  277.     }
  278.     if (mn %in% c("SID", "DOL", "MGS", "CAL", "FE", "WUS")) {  
  279.        return(tpf.omesod(mw,1,999));
  280.     }
  281.     if (mn %in% c("WO")) {
  282.        return(tpf.omesod(mw,3,999));
  283.     }
  284.     if (mn %in% c("NE")) {  
  285.        return(tpf.omesod(mw,4,2));
  286.     }
  287.     if (mn %in% c("ST")) {  
  288.        return(tpf.omesod(mw,46,999));
  289.     }
  290.     if (mn %in% c("TLC")) {  
  291.        return(tpf.omesod(mw,11,999));
  292.     }
  293.     if (mn %in% c("CHL")) {  
  294.        return(tpf.omesod(mw,14,4));
  295.     }
  296.     if (mn %in% c("CHD")) {  
  297.        return(tpf.omesod(mw,12,999));
  298.     }
  299.     if (mn %in% c("OSU")) {  
  300.        return(tpf.omesod(mw,30,999));
  301.     }
  302.     if (mn %in% c("KRN")) {  
  303.        return(tpf.omesod(mw,43,999));
  304.     }
  305.     if (mn %in% c("STP")) {  
  306.        return(tpf.omesod(mw,23,999));
  307.     }
  308.     if (mn %in% c("HU")) {  
  309.        return(tpf.omesod(mw,13,999));
  310.     }
  311.     if (mn %in% c("CRN")) {
  312.         return(tpf.omesod(mw,3,999));
  313.     }
  314. }
  315.  
  316.