require("MASS")
major<-c("SiO2","TiO2","Al2O3","Fe2O3","FeO","MnO","MgO","CaO","Na2O","K2O","P2O5","Li2O")
LREE<-c("La","Ce","Pr","Nd","Sm","Gd")
WR<-data.matrix(0)

# Loading data file
load.file<-function(){
   if(.Platform$OS.type=="windows"){
      filename<-file.choose()
   }else{
      filename<-readline("Enter the filename: ")
   }
   x<-read.table(filename,sep="\t",dec=".",check.names=F,fill=T) 
      # For decimal commas set: dec="," 

   # Ensures that all the necessary variables are there, even if empty
   add.on<-function(param,what){
       if(any(colnames(WR)==param))return(WR)
       options(show.error.messages = F)
       try(WR<-cbind(WR,what))
       try(colnames(WR)[ncol(WR)]<-param)
       options(show.error.messages = T)
       return(WR)
   }   
   col.names<-c(major,LREE,"Zr")
   y<-matrix(nrow=nrow(x),ncol=length(col.names),dimnames=list(rownames(x),col.names))
   WR<-cbind(x,y[,setdiff(colnames(y),colnames(x))])
   WR<-add.on("Li2O",WR[,"Li"]/0.46452/1e4) # Recalculates Li (ppm) to LiO2 (wt%) if necessary
   WR<-add.on("FeO",WR[,"FeOt"])
   WR<-add.on("FeO",WR[,"FeO*"])
   WR<-add.on("A/CNK",acnk(millicat(WR)))
   WR[WR<=0]<-NA
   #WR[is.na(WR)]<-0 # All missing values are replaced by 0
   print(WR)
   WR<<-WR
   return(WR)
}

# Calculates millications
millicat<-function(what=WR){  
   MW<-c(60.0848,79.8988,101.96128,159.6922,71.88464,70.9374,40.3044,56.0794,61.97894,
    94.1954,141.94452,29.8814) # Molecular weights for 'major'
   names(MW)<-major
   fact<-c(1,1,2,2,1,1,1,1,2,2,2,2) # Number of cations per molecule
   names(fact)<-major
   z<-t(apply(what[,major],1,function(x){x/MW[major]*fact[major]*1000}))
   return(z)
}

# Calculates A/CNK 
acnk<-function(what){
    z<-what[,"Al2O3"]/(what[,"Na2O"] + what[,"K2O"]+2*what[,"CaO"])
    return(z)
}

# Normalizes a matrix to a given sum 
normalize2total<-function(what,total=100){
    z<-t(apply(what,1,function(x,y){x/sum(x,na.rm=T)*y},y=total))
    return(z)
}

# Filters out from matrix 'where' rows in which exist all columns specified in 'what'
filter.out<-function(where,what){
    i<-apply(where[,what],1,function(x){all(!is.na(x))})
    z.names<-rownames(where)[i]
    z<-as.matrix(where[i,what])
    rownames(z)<-z.names
    colnames(z)<-what
    mode(z)<-"numeric"
    return(z)
}

# Copies a matrix/vector 'what' to clipboard (under MS Windows)
r2clip<-function(what=results) {
    if(is.null(what)) stop("No data available!")
    filename<-file("clipboard",open="w")
    what<-as.matrix(what)
    if(ncol(what)==1)colnames(what)<-""
    write.matrix(cbind(rownames(what),what),filename,sep="\t")
    close(filename)
}

zr.saturation<-function(cats=millicat(WR), T=750, Zr = subset(WR, Zr>0, "Zr")){
    T<-T+273.15
    cats<-cats[rownames(Zr),]
    x<-normalize2total(cats)
    M<-(x[,"Na2O"]+x[,"K2O"]+2*x[,"CaO"])/(x[,"Al2O3"]*x[,"SiO2"])*100
    DZr1<-exp(-3.8-0.85*(M-1)+12900/T)
    Zr.sat<-497644/DZr1
    DZr<-497644/Zr
    TZr.sat.C<-12900/(log(DZr)+3.8+0.85*(M-1))-273.15
    y<-cbind(M,Zr,round(Zr.sat,1),round(TZr.sat.C,1))
    colnames(y)<-c("M","Zr","Zr.sat","TZr.sat.C")
    results<<-y
    return(y)
}



mz.saturation<-function(cats=millicat(WR),H2O=3,Xmz=0.83){   
    MW.REE<-c(138.9055,140.12,140.9077,144.24,151.4,154.25) # Atomic weights for LREE
    names(MW.REE)<-LREE
    REE<-filter.out(WR,LREE)
    cats<-cats[rownames(REE),] # Get only samples for which all LREE are available
    ree<-t(t(REE)/MW.REE)
    reex<-apply(ree,1,sum,na.rm=T)/Xmz
    x<-normalize2total(cats,100)
    D<-(x[,"Na2O"]+x[,"K2O"]+2*x[,"CaO"])/x[,"Al2O3"]*1/(x[,"Al2O3"]+x[,"SiO2"])*100
    T.calc<-13318/(9.5+2.34*D+0.3879*sqrt(H2O)-log(reex))-273.15
    y<-cbind(D,round(T.calc,1))
    colnames(y)<-c("Dmz","Tmz.sat.C")
    results<<-y
    return(y)
}


ap.saturation<-function(Si=WR[,"SiO2"],ACNK=WR[,"A/CNK"],P2O5=data.matrix(WR)[,"P2O5"],T=750){  
    Si<-Si/100
    T<-T+273.15
    # Harrison and Watson (1984)
    A<-8400+(Si-0.5)*26400
    B<-3.1+12.4*(Si-0.5) 
    D.HW<-exp(A/T-B)
    P2O5.HW<-42/D.HW
    T.HW<-A/(log(42/P2O5)+B)-273.15
      
    # A general routine that solves nonlinear equation for T (deg C) by bisection method
    solve.T<-function(fun,tmin=0,tmax=NULL){
        T.calc<-NULL
        for(i in 1:length(Si)){
           if(ACNK[i]>1){
            ttold<-0
            tt<-1
            if(is.null(tmax)) tt.max<-T.HW[i]+273 else tt.max<-tmax 
                # H+W temperature is the only feasible maximum estimate
            tt.min<-tmin
            while(abs(ttold-tt)>0){
                ttold<-tt
                tt<-(tt.max-tt.min)/2+tt.min
                expr<-gsub("Si",Si[i],fun)
                expr<-gsub("ACNK",ACNK[i],expr)
                expr<-gsub("T",tt,expr)
                pp<-eval(parse(text=as.expression(expr)))
                if(pp>P2O5[i])tt.max<-tt else tt.min<-tt
            } 
                T.calc[i]<-tt
            }else{
                T.calc[i]<-NA
            }
        } 
    return(T.calc-273.15)
    } 


# Pichavant et al. (1992)  
    T.PV<-solve.T("42/exp((8400+(Si-0.5)*26400)/T-3.1-12.4*(Si-0.5))+(ACNK-1)*exp(-5900/T-3.22*Si+9.31)")

# Bea et al. (1992)  
    T.Bea<-solve.T("42*exp(((ACNK-1)*6429)/(T-273.15)-(8400+(Si-0.5)*26400)/T+3.1+12.4*(Si-0.5))")
    y<-cbind(T.HW,T.Bea,T.PV)
    y<-round(y,1)
    y<-cbind(round(ACNK,2),y)
    colnames(y)<-c("A/CNK","Tap.sat.C.H+W","Tap.sat.C.Bea","Tap.sat.C.Pich")
    rownames(y)<-names(P2O5)
    results<<-y
    return(y)
 }
options(warn=-1)
load.file()
WR[WR=="b.d."]<-NA
WR[WR<0]<-NA