Courbe de titrage/titration avec le logiciel R project

# Concentration en soude simulée (29 points ici avançant de 0,5 en 0,5.

soude = 0:28

soude = soude*0.5

# Valeurs de pH mesurées

pH = c(0.5 , 1.2 , 1.5 , 1.7 , 1.9 , 2 , 2.1 , 2.3 , 2.7 , 3.3 , 4.1 , 6 , 8 , 9.5 , 10.2 , 10.6 , 10.8 , 11 , 11.1 , 11.25 , 11.4 , 11.55 , 11.7 , 12.5 , 13.3 , 13.6 , 13.65 , 13.7 , 13.75)  

plot(soude,pH,type="o",ylim=c(0,14),xlab="titre de l'axe x",ylab="titre de l'axe y - à modifier")  

#Charger les fonctions d'identification des volumes à l'équivalence. Il suffit de copier-coller ce code dans la console R.

#################################

# Fonction de calcul du volume à l'équivalence par la règle des tangentes

veq_calc = function(x,y,trace=T) {                

    # On trace les tangentes et on calcule leurs équations            

    haute = local(1)            

    diff_temp = x-haute$x

    j = diff_temp[1]

    for (i in diff_temp) {            

         if ((i >0) & (j<0)) {        

            pos_sup = which(diff_temp==i)    

            pos_inf = which(diff_temp==j)    

        }        

        j = i        

    }    

    a = (y[pos_sup]-y[pos_inf])/(x[pos_sup]-x[pos_inf])   

    b= haute$y-haute$x*a       

    if (trace == T) {abline(b,a,lty=4)}            

    basse = local(2)            

    bb= basse$y-basse$x*a                      

    if (trace == T) {abline(bb,a,lty=4)}            

    bbb = (b+bb)/2            

    if (trace == T) {abline(bbb,a,lty=4)}            

    # On récupère la position et les valeurs x de tous les points dont x est situé entre les 2 points où passent les tangentes            

    # Fonction prévue pour inverser les tangentes si elles ont été faites dans un ordre différent que haut puis bas            

    position = which(x>min(basse$x,haute$x)&x<max(basse$x,haute$x))            

    if (length(position) < 3) {                  

        position = c(position[1]-1, position,position[length(position)]+1)               

    }            

    x_a_teste = x[position]            

    #x_a_teste = x[x>min(basse$x,haute$x)&x<max(basse$x,haute$x)]            

    ecart = c()            

    for (i in x_a_teste) {            

        temp = a*i+bbb        

        #points(i,temp,col="orange",pch=16)        

        ecart=c(ecart,temp)        

    }            

    ecart = ecart - y[x>min(basse$x,haute$x)&x<max(basse$x,haute$x)]                      

    position_sup = c()            

    position_inf = c()            

    j = ecart[1]            

    for (i in ecart) {            

        if ((i <0) & (j>0)) {        

            position_sup = which(ecart==i)    

            position_inf = which(ecart==j)    

        }        

        j = i        

    }            

    if ((length(position_inf) ==0) | (length(position_sup) ==0)) {            

        j = ecart[1]        

        for (i in ecart) {        

            if ((i >0) & (j<0)) {    

                position_sup = which(ecart==j)

                position_inf = which(ecart==i)

            }    

            j = i    

        }              

    }            

    position_sup = position_sup+(position[1])-1            

    position_inf = position_inf+(position[1])-1            

    a_graph = (y[position_sup]-y[position_inf])/(x[position_sup]-x[position_inf])            

    b_graph = y[position_sup] - x[position_sup]*a_graph            

    x_croisement = (bbb-b_graph)/(a_graph-a)            

    y_croisement = a_graph*x_croisement+b_graph            

    return(c(x_croisement,y_croisement))            

} 

# Fonction bornes : permet de retrouver la valeur y de la courbe quand on connait x. Permet aussi de l'afficher à l'écran.              

bornes = function(x=c(),y=c(),v=c(),h=c(),trace=T,col="blue") {                

    if (length(h)==0) {            

        lim = v[1]        

        liste = x  

    j = min(x)            

    k = max(x)       

    } else if (length(v)==0) {

    lim = h[1]

    liste = y 

    j = min(y)            

    k = max(y) 

    }            

    for (i in liste) {            

        if (i < lim) {        

            if (i > j) {    

                j = i

                position_sup = which(liste==i)

            }    

        }        

        if (i> lim) {        

            if (i < k) {    

                k = i

                position_inf = which(liste==i)

            }    

        }        

    } 

    a_graph = (y[position_sup]-y[position_inf]) / (x[position_sup]-x[position_inf])            

    b_graph = y[position_sup] - x[position_sup]*a_graph  

    if (length(h)==0) {            

        x_croisement = lim            

        y_croisement = a_graph*x_croisement+b_graph   

    }   else if (length(v)==0) {

        x_croisement = (lim - b_graph)/ a_graph          

        y_croisement = lim 

    }         

    if (trace==T) {            

        arrows(0,y_croisement,x_croisement,y_croisement,lty=3,col=col,angle=0)        

        arrows(x_croisement,0,x_croisement,y_croisement,lty=3,col=col,angle=0)        

    }    

    return(c(x_croisement,y_croisement))            

}             

local = function(num) {                

    cat("La tangeante  ",num,"va être tracée.\n")            

    readline("\nCliquer sur le point où afficher la tangeante :\n(si ne fonctionne pas, presser Enter)\n")            

    cat("La tangeante  ",num,"a été tracée.\n\n")            

    return(locator(n=1))            

}       

# Permet de placer manuellement une légende sur le graphique.         

txt_local = function(txt="txt",col="blue") {                

    readline("\nCliquer sur le point où afficher une légende :\n(si ne fonctionne pas, presser Enter)\n")      

    l = locator(n=1)            

    text(l$x,l$y, txt,col=col)            

} 

# Calculer Veq - il suffit de faire pareil pour veq2 

veq = veq_calc(soude,pH)

# Afficher Veq

veq = bornes(soude,pH,veq,col="blue")

# Indiquer Veq sur la légende en cliquant sur le graphique

txt_local("Veq")