Retro programming nostalgia V: L’Equilibrio e la Titolazione Acido/Base (Parte II)

Questo secondo articolo continua il mio personale viaggio retro-computazionale delle titolazioni acido/base. Nell’articolo precedente, ho mostrato come calcolare un equilibrio acido base per acidi e basi forti. In questo articolo, vengono descritte anche le subroutines per le titolazioni di acidi e basi deboli monoprotici. Il metodo che ho usato risolve in modo esatto il calcolo dei pH e si basa su articolo pubblicato sulla rivista di chimica “Rassegna chimica” da Prof Luigi Campanella (e Dr G. Visco) nel 1985. Ricevetti dall’autore stesso una copia dell’articolo quando frequentavo il suo corso di chimica analitica presso l’Università “la Sapienza” di Roma. Ricordo che scrivere un programma per lo studio delle titolazioni non solo fu divertente e stimolante ma mi aiutò molto a capire a fondo l’argomento. Pertanto raccomando il giovane lettore di provare a convertire il programma in un linguaggio moderno a voi più familiare (per esempio il Python) per meglio comprenderne il funzionamento.

Titolazione di un acido debole con base forte

Consideriamo le concentrazioni dello ioni [H_3O^+] nelle varie fasi della titolazione. Indichiamo con C_a  la concentrazione dell’acido in una  soluzione  di volume V_A, con C_T la concentrazione del titolante e V_T il volume di titolante aggiunto.

La costante di dissociazione dell’acido è uguale a 

K_a=\frac{[H_3O^+] [B^-]}{[HB]}

da cui

K_a[HB]=[H_3O^+] [B^-]

Usando la relazione  C_a=[HB]+[B^{-}] e  sostituendo [HB]=C_a-[B^{-}] so ottiene

K_a(C_a-[B^{-}])=[H_3O^+] [B^-]

da cui possiamo ricavare 

[B^-]= \frac{K_aC_a}{K_a+[H_3O^+]}

Sostituendo la precendente relazione nella condizione di elettronegatività  

[H_3O^+]=[OH^{-}]+[B^{-}]

si ottiene

[H_3O^+]=\frac{k_W}{[H_3O^+]} + \frac{K_aC_a}{K_a+[H_3O^+]}

espandendo la precedente relazione, si ottiene la seguente equazione cubica

[H_3O^+]^3 +K_a[H_3O^+]^2-(K_aC_a+K_W)[H_3O^+]-K_WK_a=0

Risolvendo l’equazione usando, per esempio, il metodo di Newton si può trovare il pH iniziale.

Dopo l’addizione del titolante la condizione di elettroneutralità cambia in

[H_3O^+]+[Me^+]=[OH^{-}]+[B^{-}] 

con [Me^+]=C_b and  [HB]+[B^{-}]=C_a

I valori di queste concentrazioni cambiano con l’addizione del volume del titolante (v) come

C_a=\frac{C_a' V}{V+v}

C_b=\frac{C_b' V}{V+v}

Con V, C_a'  e C_b' il volume iniziale della soluzione, la  concentrazione iniziale dell’acido e quella della base titolante.

Quindi usando  condizione di elettroneutralità come nella  precedente derivazione 

si giunge alla nuova equazione

[H_3O^+]^3 +(K_a+C_b)[H_3O^+]^2+(K_aC_b-K_aC_a+K_W)[H_3O^+]-K_WK_a=0

Nella figura che segue viene mostrato un esempio di curva di titolazione dui un acido debole monoprotico con una base forte.

Figure 1: Esempio di curva di titolazione di un acido debole con una base forte.

La titolazione di una base debole con un acido forte

Procedendo come nel caso precedente, si ottiene dapprima un’ equazione da cui calcolare la concentrazione degli ioni idrossilici simile alla precedente

[OH^-]^3 +K_b[O^-]^2-(KbCb+K_W)[O^-]-K_WK_b=0

Con l’aggiunta del volume del titolante acido la precedente equazione cambia in 

 [OH^-]^3 +(K_b+C_a)[OH^-]^2+(K_bC_a-K_bC_b+K_W)[OH^-]-K_WK_b=0

Nella figura che segue viene mostrato un esempio di curva di titolazione dui un base debole con un acido forte.

Figure 2: Esempio di curva di titolazione di un base debole con un acido forte.

La titolazione di una base debole con un acido debole

Per completezza aggiungiamo le equazioni per la titolazione di un acido debole con una base debole. Questo tipo di titolazioni non viene quasi mai effettuata poichè è molto difficile identificare il punto finale.

La derivazione segue lo stessa procedura basata sulla condizione di elettroneutralità della soluzione. Di seguito vengono riportate solo le equazioni finali di quarto ordine  usate nel programma in linguaggio BASIC.

La concentrazione protonica iniziale dell’acido debole è data dalla seguente equazione come nel caso precedente

[H_3O^+]^3 +K_a[H_3O^+]^2-(KaCa+K_W)[H_3O^+]-K_WK_a=0

Con l’aggiunta del volume della base debole come titolante  la precedente equazione cambia in 

\begin{aligned}K_b[H_3O^+]^4 +(K_aK_b+K_bC_b+K_w)[H_3O^+]^3 \\ -(K_aK_bC_b-K_aC_aK_b +K_aK_W-K_bK_W)[H_3O^+]^2 \\-(K_aK_bK_W+K_aC_aK_W+K_W^2)[H_3O^+]-K_W^2K_a  =0 \end{aligned}

Nella figura che segue viene mostrato un esempio di curva di titolazione dui un acido debole con una base debole.

Figure 3: Exempio di curva di titolazione di un acido debole con una base debole

Il calcolo delle radici dei polinomi delle concentrazioni protoniche con il di Newton-Raphson

I valori delle concentrazione di [H_3O^+] e [OH^-] sono ottenuti usando il metodo delle tangenti o di Newton-Raphson. Il metodo consiste nel sostituire alla curva y=f([H_3O^3]) la tangente alla curva stessa, partendo da un punto di concentrazione vicino a quello cercato. Per esempio per primo punto ([H_3O^+]_0) si pu\’o considerare il valore corrispondente a pH=0.1

Si assume quindi che il valore approssimato della radice sia il valore dell’ascissa 

[H_3O^+]_1= [H_3O^+]_0 -\frac{f([H_3O^+]_0)}{ f'([H_3O^+]_0)}

in cui la retta tangente  (di inclinazione f'([H_3O^+]_0)) interseca l’asse delle [H_3O^+] nell’ intevallo corrispondente 

al range di pH 1-14.  Procedendo in iterativamente  si usa il nuovo punto a [H_3O^+]_1 per ottenere una nuova approssimazione usando la formula generale

[H_3O^+]_{n+1}= [H_3O^+]_n -\frac{f([H_3O^+]_n)}{ f'([H_3O^+]_n)}

che aumenta l’accuratezza del valore della radice. La procedure viene interrotta quando il livello di accuratezza viene raggiunto.

Il programma in linguaggio BASIC

Le curve mostrate in questo articolo sono state ottenute con il programma in BASIC riportato in appendice. Questo è una versione aggiornata del programma riportato nel precedente articolo.Il programma, originariamente sviluppato in MSX Basic, è stato adattato a MS Qbasic e può essere compilato usando il superbo compilatore open source QB64 (https://www.qb64.org/portal).

Questa versione del programma permette anche la simulazione delle titolazioni con base debole o acido (mono) protico deboli. Il programma usa valori prestabiliti per i parametri di titolazione che possono essere comunque modificati. La titolazione può essere effettuata passo passo aggiungendo aliquote di acido e basi o in un passo inserendo un numero negativo.  

In questa appendice è riportato il listato programma in BASIC. Il programma viene distribuito per scopi didattici e non viene garantito che sia esente da bug di programmazione. Per cui se trovate errori o avete suggerimenti, mandatemi pure un commento. Se poi lo trovate interessante potreste anche provare a svilupparlo ulteriormente.

Schermo iniziale del programma per la selezione della titolazione
Schermo per il cambiamento dei parametri

APPENDICE

In questa appendice è riportato il listato programma in BASIC. Per il momento il programma calcolo solo le curve di titolazioni Il programma viene distribuito per scopi didattici e non viene garantito che sia esente da bug di programmazione. Per cui se trovate e(o,;-))rrori o avete suggerimenti, mandate un commento. Se poi lo trovate interessante potreste anche provare a svilupparlo ulteriormente.

'================================================================
' PROGRAMMA  : TITOLAZIONI_AB.BAS
' DESCRIZIONE: QUESTO PROGRAMMA CALCOLA LE CURVE DI TITOLAZIONI
'              ACIDO BASE
' AUTORE     : D. ROCCATANO
' ANNO DI PRIMA IDEAZIONE (MSX BASIC)   : 1987
' ANNO DI ADATTAMENTO E AGGIORNAMENTO (QB64 BASIC): 2019
'================================================================

Dim Shared XI, YI, YF, SY, SX, XF, CI, CF, SE
Dim Shared As Double G, KW

'line colors

CF = 2 'colore delle linee del grafico
CI = 12 'Curva di titolazione

'Screen parameters
G = 1 / 448
XI = 150
YI = 485
YF = 38
SY = 448
SX = 500
XF = XI + SX
eps = 1.0E-6

Randomize Timer
Screen _NewImage(1000, 600, 256)
Color 2, 15
Out &H3C8, 0
Out &H3C9, 63
Out &H3C9, 63
Out &H3C9, 63

' PARAMETRI DI DEFAULT
VC = 50
' Concentrazione del titolante
CT = 0.02
' Concentrazione delle soluzione
CC = 0.001
' Aliquota di volume aggiunto
DV = 0.01
KW = 1.0E-14
'
GENC = 0
CMAX = 0.5
CMIN = 0.01
VMAX = 100
VMIN = 10

'    ====  PROGRAMMA  PRINCIPALE =====
'
Do
    Do:
        Cls
        Color 8: Locate 5, 5: Print "CURVE DI TITOLAZIONE ACIDO-BASE"
        Color 2: Locate 10, 10: Print "[1] ACIDO FORTE+BASE FORTE"
        Color 3: Locate 12, 10: Print "[2] ACIDO DEBOLE+BASE FORTE"
        Color 4: Locate 14, 10: Print "[3] BASE FORTE+ACIDO FORTE"
        Color 5: Locate 16, 10: Print "[4] BASE DEBOLE+ACIDO FORTE"
        Color 6: Locate 18, 10: Print "[5] ACIDO DEBOLE+BASE DEBOLE"
        Color 7: Locate 20, 10: Print "[6] CAMBIA PARAMETRI"
        Color 9: Locate 22, 10: Print "[7] ESCI "
        Color 8: Locate 25, 12: Input "Scegli una opzione [1-7]:"; SE
    Loop Until SE > 0 And SE < 8

    If SE <> 6 Then
        Cls
        If GENC = 0 Then
            Locate 20, 5: Input "CONCENTRAZIONE DEL CAMPIONE : "; CC
        Else
            CC = CMIN + Rnd * CMAX
        End If
        If SE = 2 Or SE = 5 Then
            Locate 21, 5: Input "COSTANTE ACIDO DEBOLE : "; KA
        End If
        If SE = 4 Or SE = 5 Then
            Locate 22, 5: Input "COSTANTE BASE DEBOLE : "; KB
        End If

        If SE > 6 Then End
        'Disegna la cornice e gli assi del grafico
        FRAME
    End If

    ' Seleziona una subroutine
    Select Case SE
        Case 1
            AFBF VC, CC, CT, DV
        Case 2
            ADBF VC, CC, CT, KA, DV, eps
        Case 3
            AFBF VC, CC, CT, DV
        Case 4
            BDAF VC, CC, CT, KB, DV, eps
        Case 5
            ADBD VC, CC, CT, KA, KB, DV, eps
        Case 6
            CAMBIA_PARAM VC, CT, DV, GENC
    End Select

    Locate 36, 1: Print "UN'ALTRA TITOLAZIONE? (S/N): ";
    Do: K$ = UCase$(InKey$)
    Loop Until K$ = "S" Or K$ = "N"
Loop Until K$ = "N"

End

' ==== SUBROUTINES ====
Sub CAMBIA_PARAM (VC, CT, DV, GENC)
    ' INSERIMENTO DEI DATI
    ' INIZIALI
    Cls
    Color 8
    Locate 5, 5: Print "INSERISCI DATI INIZIALI USANDO LA STESSA UNITA' DI MISURA"
    Locate 10, 5: Input "VOLUME CAMPIONE   : "; VC
    Locate 11, 5: Input "CONC. TITOLANTE   : "; CT
    Locate 12, 5: Input "ALIQUOTA AGGIUNTA : "; DV

    ' DISEGNA LO SCHERMO INIZIALE
    Do:
        Locate 13, 5: Print "GENERA CAUSALMENTE LA CONCENTRAZIONE DEL CAMPIONE? (S/N)"
        K$ = UCase$(InKey$)
    Loop Until K$ = "S" Or K$ = "N"

    If K$ = "N" Then
        GENC = 0
    End If
    If K$ = "S" Then
        GENC = 1
    End If
End Sub

Sub FRAME
    '   DRAW THE GRAPH FRAME
    Shared XI, XF, YI, YF, SX, SY, CF, SE
    Cls
    '   DRAW GRAPH FRAME
    Line (XI, YF)-(XF, YI), CF, B
    '   DRAW x-ticks
    For I = XI To XF Step 5
        Line (I, YF)-(I, YF + 5), CF
        Line (I, YI)-(I, YI - 5), CF
    Next
    For I = XI To XF Step 50
        Line (I, YF)-(I, YF + 10), CF
        Line (I, YI)-(I, YI - 10), CF
    Next

    '   DRAW y-tics
    For I = YI To YF Step -4
        Line (XI, I)-(XI + 5, I), CF
        Line (XF, I)-(XF - 5, I), CF
    Next
    For I = YI To YF Step -32
        Line (XI, I)-(XI + 10, I), CF
        Line (XF, I)-(XF - 10, I), CF
    Next
    Color 8
    Locate 15, 12: Print "pH"
    For I = 1 To 14 Step 1
        Locate 3 + (I - 1) * 2, 17: Print Using "##"; 15 - I
        Locate 3 + (I - 1) * 2, 83: Print Using "##"; 15 - I
    Next
End Sub

Sub SCALAVOL (RV)
    '
    ' STAMPA I VALORI DELLA SCALA DEI VOLUMI AGGIUNTI
    '
    inc = RV / 10
    tmp = 0
    Locate 32, 18: Print Using "###.#"; tmp;
    tmp = inc
    For i = 1 To 10 Step 1
        Locate 32, 22 + (i - 1) * 6.3: Print Using "###.#"; tmp;
        tmp = tmp + inc
    Next
    Locate 33, 70: Print "VOL. AGGIUNTO"

End Sub

Sub AFBF (VC, CC, CT, DV)

    Shared XI, SE, YF, YI
    Shared KW

    Color 8
    If SE = 1 Then
        Locate 1, 2: Print "TITOLAZIONE ACIDO FORTE - BASE FORTE"
    Else
        Locate 1, 2: Print "TITOLAZIONE BASE FORTE - ACIDO FORTE"
    End If
    Locate 5, 87: Print Using "VOLUME SOLUZIONE: ##.##"; VC
    Locate 6, 87: Print Using "CONC. TITOLANTE: ##.##"; CT
    Locate 7, 87: Print Using "DEFAULT ALIQU. DI TITOLANTE: ##.##"; DV

    VT = DV ' Volume iniziale di titolante
    AUTO = 0
    '
    ' Calcolo del volume necessario a completare la titolazione
    '
    RV = (1.1 + 0.9 * Rnd) * (VC * CC) / CT
    ' Calcola il fattore di scala per il volume

    GF = 500 / RV

    SCALAVOL RV

    PP = 0
    MAXV = RV

    ' Esegue la titoalazione
    Do:
        H = (VC * CC - VT * CT) / (VC + VT)
        f = -H

        If H = 0 Then Ph = 7
        If H < 0 Then
            If f <= KW / f Then
                Ph = 14 + Log10#(f + KW / (f + 1E-7))
            Else
                Ph = 14 + Log10#(f + KW / f)
            End If
        End If

        If H > 0 Then
            If H < KW / H Then
                Ph = -Log10#(H + KW / (H + 1E-7))
            Else
                Ph = -Log10#(H + KW / H)
            End If
        End If

        ' SE LA TITOLAZIONE E' DI UNA BASE FORTE CON UN ACIDO FORTE
        If (SE = 3) Then Ph = 14 - Ph
        ' DISEGNA LA CURVA
        If (PP = 0) Then
            PSet (XI + VT * GF, YI - SY * Ph / 14), CI
        Else
            Line -(XI + VT * GF, YI - SY * Ph / 14), CI
        End If
        PP = PP + 1

        'aggiunti aliquota
        NALI = 1
        If (AUTO = 0) Then
            Locate 34, 5: Print "NUMERO (n) DI ALIQUOTE (dv) DA AGGIUNGERE. IL NUMERO n E' UN FATTORE DI dv."
            Locate 35, 5: Print "                                                      "
            Locate 35, 5: Input "UN VALORE n<0 TERMINA LA TITOLAZIONE:"; inval
        End If
        If inval > 0 Then NALI = inval
        If inval < 0 Then AUTO = 1: NALI = 1
        VT = VT + DV * NALI
        Locate 30, 87: Print Using "VOLUME AGGIUNTO: ##.##"; VT;
    Loop Until VT > MAXV

    Locate 31, 87: Print Using "VOLUME AGGIUNTO ALL'EQUIVALENZA:##.##"; VC * CC / CT

    Exit Sub

End Sub

Sub ADBF (VC, CC, CT, KA, DV, eps)

    '
    '======================================================================
    '  DESCRIPTION: Curve di titolazione acido debole base debole
    '======================================================================
    ' Author:        Danilo Roccatano <danilo.roccatano@gmail.com>
    ' Version:       0.1
    ' Created:       2021-03-03 18:16
    '======================================================================
    '

    Dim As Double H, H1, f, fp

    Shared XI, SE, YF, YI
    Shared KW

    Color 8

    Locate 1, 2: Print "TITOLAZIONE ACIDO DEBOLE - BASE FORTE"
    Locate 5, 87: Print Using "VOLUME SOLUZIONE: ##.##"; VC
    Locate 6, 87: Print Using "CONC. TITOLANTE: ##.##"; CT
    Locate 7, 87: Print Using "COSTANTE ACIDO DEBOLE: ##.##^^^^"; KA
    Locate 8, 87: Print Using "ALIQU. DI TITOLANTE: ##.##"; DV
    Locate 9, 87: Print Using "CONC. SOLUZIONE: ##.##^^^^"; CC

    '
    ' Set up
    '
    H = 0

    VT = 0 ' Volume iniziale di titolante
    AUTO = 0

    '
    ' Calcolo del volume necessario a completare la titolazione
    '

    RV = (1.1 + 0.9 * Rnd) * (VC * CC) / CT

    ' Calcola il fattore di scala per il volume

    GF = 500 / RV

    SCALAVOL RV

    PP = 0
    MAXV = RV

    '
    ' First titolation point
    '

    pH = 1.0
    H1 = 10. ^ (-pH)

    '
    ' Use Newton-Raphson to find the first point (zero titrant)
    '
    ncycles = 0

    Do:
        H = H1
        f = H ^ 3 + KA * H ^ 2 - (KA * CC + KW) * H - KA * KW
        fp = 3 * H ^ 2 + 2 * KA * H - (KA * CC + KW)
        H1 = H - (f / fp)
        ncycles = ncycles + 1
    Loop Until (Abs(H1 - H) <= eps Or ncycles > 50)

    pH = -Log10#(H1)

    PSet (XI + VT * GF, YI - SY * pH / 14), CI
    '
    ' Adding titrant
    '

    Do:
        CA = CC * VC / (VC + VT)
        CB = CT * VT / (VC + VT)

        T2 = KA + CB
        T1 = KA * CB - KA * CA - KW
        T0 = KA * KW
        '
        ' Newton/Raphson method
        '
        ncycles = 0

        Do:
            H = H1
            f = H ^ 3 + T2 * H ^ 2 + T1 * H - T0
            fp = 3 * H ^ 2 + 2 * T2 * H + T1
            H1 = H - (f / fp)
            ncycles = ncycles + 1
        Loop Until (Abs(H1 - H) <= eps Or ncycle > 100)

        pH = -Log10#(H1)

        ' DISEGNA LA CURVA

        Line -(XI + VT * GF, YI - SY * pH / 14), CI

        PP = PP + 1

        'Aggiunge aliquote di titolante

        NALI = 1
        If (AUTO = 0) Then
            Locate 34, 5: Print "NUMERO (n) DI ALIQUOTE (dv) DA AGGIUNGERE. IL NUMERO n E' UN FATTORE DI dv."
            Locate 35, 5: Print "                                                      "
            Locate 35, 5: Input "UN VALORE n<0 TERMINA LA TITOLAZIONE:"; inval
        End If
        If inval > 0 Then NALI = inval
        If inval < 0 Then AUTO = 1: NALI = 1
        VT = VT + DV * NALI
        Locate 30, 87: Print Using "VOLUME AGGIUNTO: ##.##"; VT;
    Loop Until VT > MAXV

    Locate 31, 87: Print Using "VOLUME AGGIUNTO ALL'EQUIVALENZA:##.##"; VC * CC / CT


End Sub

Sub BDAF (VC, CC, CT, KB, DV, eps)

    '======================================================================
    '  DESCRIPTION: Curve di titolazione base debole acido forte
    '======================================================================
    ' Author:        Danilo Roccatano <danilo.roccatano@gmail.com>
    ' Version:       0.1
    ' Created:       2021-03-03 18:16
    '======================================================================

    Dim H, H1, f, fp As Double
    Shared XI, SE, YF, YI
    Shared KW

    Color 8

    Locate 1, 2: Print "TITOLAZIONE ACIDO DEBOLE - BASE FORTE"
    Locate 5, 87: Print Using "VOLUME SOLUZIONE: ##.##"; VC
    Locate 6, 87: Print Using "CONC. TITOLANTE: ##.##"; CT
    Locate 7, 87: Print Using "COSTANTE BASE DEBOLE: ##.##^^^^"; KB
    Locate 8, 87: Print Using "ALIQU. DI TITOLANTE: ##.##"; DV
    Locate 9, 87: Print Using "CONC. SOLUZIONE: ##.##^^^^"; CC


    '
    ' Set up
    '

    H = 0
    VT = DV ' Volume iniziale di titolante
    AUTO = 0

    '
    ' Calcolo del volume necessario a completare la titolazione
    '

    RV = (1.1 + 0.9 * Rnd) * (VC * CC) / CT

    ' Calcola il fattore di scala per il volume

    GF = 500 / RV

    SCALAVOL RV

    PP = 0
    MAXV = RV

    '
    ' First titolation point
    '

    Ph = 1.
    OH1 = 10 ^ (-Ph)

    '
    ' Use Newton-Raphson to find the first point (zero titrant)
    '

    ncycles = 0
    Do:
        OH = OH1
        f = OH ^ 3 + KB * OH ^ 2 - (KB * CC + KW) * OH - KB * KW
        fp = 3 * OH ^ 2 + 2 * KB * OH - (KB * CC + KW)
        OH1 = OH - (f / fp)
        ncycles = ncycles + 1
    Loop Until (Abs(OH1 - OH) <= eps Or ncycles > 50)

    Ph = 14 + Log10#(OH1)

    PSet (XI + VT * GF, YI - SY * Ph / 14), CI

    '
    ' Adding titrant
    '

    Do:
        CB = CC * VC / (VC + VT)
        CA = CT * VT / (VC + VT)

        T2 = KB + CA
        T1 = KB * CA - KB * CB - KW
        T0 = KB * KW

        '
        ' Newton-Raphson method
        '

        Do:
            OH = OH1
            f = OH ^ 3 + T2 * OH ^ 2 + T1 * OH - T0
            fp = 3 * OH ^ 2 + 2 * T2 * OH + T1
            OH1 = OH - (f / fp)
            ncycles = ncycles + 1
        Loop Until (Abs(OH1 - OH) <= eps Or ncycles > 50)

        ' DISEGNA LA CURVA

        Ph = 14 + Log10#(OH1)

        Line -(XI + VT * GF, YI - SY * Ph / 14), CI

        PP = PP + 1

        'Aggiunge aliquote di titolante

        NALI = 1
        If (AUTO = 0) Then
            Locate 34, 5: Print "NUMERO (n) DI ALIQUOTE (dv) DA AGGIUNGERE. IL NUMERO n E' UN FATTORE DI dv."
            Locate 35, 5: Print "                                                      "
            Locate 35, 5: Input "UN VALORE n<0 TERMINA LA TITOLAZIONE:"; inval
        End If
        If inval > 0 Then NALI = inval
        If inval < 0 Then AUTO = 1: NALI = 1
        VT = VT + DV * NALI
        Locate 30, 87: Print Using "VOLUME AGGIUNTO: ##.##"; VT;
    Loop Until VT > MAXV

    Locate 31, 87: Print Using "VOLUME AGGIUNTO ALL'EQUIVALENZA:##.##"; VC * CC / CT

End Sub

Sub ADBD (VC, CC, CT, KA, KB, DV, eps)
    '======================================================================
    '  DESCRIPTION: Curve di titolazione acido debole base debole
    '======================================================================
    ' Author:        Danilo Roccatano <danilo.roccatano@gmail.com>
    ' Version:       0.1
    ' Created:       2021-03-03 18:16
    '======================================================================
    '

    Dim H, H1, f, fp As Double
    Shared XI, SE, YF, YI
    Shared KW

    Color 8

    Locate 1, 2: Print "TITOLAZIONE ACIDO DEBOLE - BASE FORTE"
    Locate 5, 87: Print Using "VOLUME SOLUZIONE: ##.##"; VC
    Locate 6, 87: Print Using "CONC. TITOLANTE: ##.##"; CT
    Locate 7, 87: Print Using "COSTANTE ACIDO DEBOLE: ##.##^^^^"; KA
    Locate 8, 87: Print Using "COSTANTE BASE DEBOLE: ##.##^^^^"; KB
    Locate 9, 87: Print Using "ALIQU. DI TITOLANTE: ##.##"; DV

    '
    ' Set up
    '
    H = 0
    VT = 0 ' Volume iniziale di titolante
    AUTO = 0

    '
    ' Calcolo del volume necessario a completare la titolazione
    '

    RV = (1.1 + 0.9 * Rnd) * (VC * CC) / CT

    ' Calcola il fattore di scala per il volume

    GF = 500 / RV

    SCALAVOL RV

    PP = 0
    MAXV = RV

    '
    ' First titolation point
    '

    Ph = 1.0
    H1 = 10 ^ (-Ph)
    '
    ' Use Newton-Raphson to find the first point (zero titrant)
    '  ncycles=0

    Do:
        H = H1
        f = H ^ 3 + KA * H ^ 2 - (KA * CC + KW) * H - KA * KW
        fp = 3 * H ^ 2 + 2 * KA * H - (KA * CC + KW)
        H1 = H - (f / fp)
        ncycles = ncycles + 1
    Loop Until (Abs(H1 - H) <= eps Or ncycles > 50)

    Ph = -Log10#(H1)

    PSet (XI + VT * GF, YI - SY * Ph / 14), CI

    '
    ' Adding titrant
    '

    Do:
        CA = CC * VC / (VC + VT)
        CB = CT * VT / (VC + VT)
        T4 = KB
        T3 = KA * KB + KB * CB + KW
        T2 = KA * KB * CB - KA * CA * KB + KA * KW - KB * KW
        T1 = KA * KB * KW + KA * CA * KW + KW * KW
        T0 = KA * KW * KW

        '
        ' USA DI NUOVO IL METODO DI NEWTON-RAPHSON
        '
        Do:
            H = H1
            f = T4 * H ^ 4 + T3 * H ^ 3 + T2 * H ^ 2 - T1 * H - T0
            fp = 4 * T4 * H ^ 3 + 3 * T3 * H ^ 2 + 2 * T2 * H - T1
            H1 = H - (f / fp)
        Loop Until (Abs(H1 - H) <= eps Or ncycle > 100)

        Ph = -Log10#(H1)

        '
        ' DISEGNA LA CURVA
        '

        Line -(XI + VT * GF, YI - SY * Ph / 14), CI

        PP = PP + 1

        '
        'Aggiunge aliquote di titolante
        '

        NALI = 1
        If (AUTO = 0) Then
            Locate 34, 5: Print "NUMERO (n) DI ALIQUOTE (dv) DA AGGIUNGERE. IL NUMERO n E' UN FATTORE DI dv."
            Locate 35, 5: Print "                                                      "
            Locate 35, 5: Input "UN VALORE n<0 TERMINA LA TITOLAZIONE:"; inval
        End If
        If inval > 0 Then NALI = inval
        If inval < 0 Then AUTO = 1: NALI = 1
        VT = VT + DV * NALI
        Locate 30, 87: Print Using "VOLUME AGGIUNTO: ##.##"; VT;


    Loop Until VT > MAXV

    Locate 31, 87: Print Using "VOLUME AGGIUNTO ALL'EQUIVALENZA:##.##"; VC * CC / CT

End Sub

Function Log10# (value As Double) Static
    Log10# = Log(value) / Log(10.#)
End Function

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.