Copolymerization (polykin.copolymerization)

This module implements methods and classes to model and analyze binary and multicomponent copolymerization systems.

CopoDataset_Ff dataclass

Dataclass for instantaneous copolymerization data of the form F(f).

Source code in src/polykin/copolymerization/copodataset.py

@dataclass(frozen=True)
class CopoDataset_Ff():
    """Dataclass for instantaneous copolymerization data of the form F(f)."""
    name: str
    f1: FloatVector
    F1: FloatVector
    scale_f1: Union[FloatVector, float] = 1.0
    scale_F1: Union[FloatVector, float] = 1.0
    weight: float = 1.0

CopoDataset_Fx dataclass

Dataclass for drift copolymerization data of the form F1(x).

Source code in src/polykin/copolymerization/copodataset.py

@dataclass(frozen=True)
class CopoDataset_Fx():
    """Dataclass for drift copolymerization data of the form F1(x)."""
    name: str
    f10: float
    x: FloatVector
    F1: FloatVector
    # scale_x: Union[FloatVector, float] = 1.0
    scale_F1: Union[FloatVector, float] = 1.0
    weight: float = 1.0

CopoDataset_fx dataclass

Dataclass for drift copolymerization data of the form f1(x).

Source code in src/polykin/copolymerization/copodataset.py

@dataclass(frozen=True)
class CopoDataset_fx():
    """Dataclass for drift copolymerization data of the form f1(x)."""
    name: str
    f10: float
    x: FloatVector
    f1: FloatVector
    # scale_x: Union[FloatVector, float] = 1.0
    scale_f1: Union[FloatVector, float] = 1.0
    weight: float = 1.0

CopoFitResult dataclass

Dataclass for copolymerization fit results.

ATTRIBUTE	DESCRIPTION
method	Name of the fit method. TYPE: str
r1	Reactivity ratio of M1. TYPE: float
r2	Reactivity ratio of M2 TYPE: float
alpha	Significance level. TYPE: float
ci_r1	Confidence interval of r1. TYPE: float
ci_r2	Confidence interval of r2. TYPE: float
se_r1	Standard error of r1. TYPE: float
se_r2	Standard error of r2. TYPE: float
cov	Scaled variance-covariance matrix. TYPE: Float2x2Matrix
plots	Dictionary of plots. TYPE: dict[str, tuple[Figure, Axes]]
M1	Name of M1. TYPE: str
M2	Name of M2. TYPE: str

Source code in src/polykin/copolymerization/fitting.py

@dataclass
class CopoFitResult():
    r"""Dataclass for copolymerization fit results.

    Attributes
    ----------
    method : str
        Name of the fit method.
    r1 : float
        Reactivity ratio of M1.
    r2: float
        Reactivity ratio of M2
    alpha : float
        Significance level.
    ci_r1 : float
        Confidence interval of r1.
    ci_r2: float
        Confidence interval of r2.
    se_r1 : float
        Standard error of r1.
    se_r2 : float
        Standard error of r2.
    cov : Float2x2Matrix
        Scaled variance-covariance matrix.
    plots : dict[str, tuple[Figure, Axes]]
        Dictionary of plots.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    """

    method: str
    r1: float
    r2: float
    alpha: float
    ci_r1: float
    ci_r2: float
    se_r1: float
    se_r2: float
    cov: Optional[Float2x2Matrix]
    plots: dict[str, tuple[Figure, Axes]]
    M1: str = 'M1'
    M2: str = 'M2'

    def __repr__(self):
        s1 = \
            f"method:  {self.method}\n" \
            f"M1:      {self.M1}\n" \
            f"M2:      {self.M2}\n" \
            f"r1:      {self.r1:.2E}\n" \
            f"r2:      {self.r2:.2E}\n"
        if self.se_r1 is not None:
            s2 = \
                f"alpha:   {self.alpha:.2f}\n" \
                f"se_r1:   {self.se_r1:.2E}\n" \
                f"se_r2:   {self.se_r2:.2E}\n" \
                f"ci_r1:   {self.ci_r1:.2E}\n" \
                f"ci_r2:   {self.ci_r2:.2E}\n"
        else:
            s2 = ""
        if self.cov is not None:
            s3 = f"cov:     {pprint_matrix(self.cov, nspaces=9)}\n"
        else:
            s3 = "\n"
        return s1 + s2 + s3

DriftDataset

Dataset of binary monomer drift copolymerization data.

Container for binary monomer drift copolymerization data, \(f_1\) vs \(x\).

PARAMETER	DESCRIPTION
M1	Name of M1. TYPE: str
M2	Name of M2. TYPE: str
x	Vector of total molar conversion, \(x\). TYPE: FloatVectorLike
f1	Vector of monomer composition, \(f_1\). TYPE: FloatVectorLike
sigma_x	Absolute standard deviation of \(x\). TYPE: Union[float, FloatVectorLike] DEFAULT: 0.05
sigma_f1	Absolute standard deviation of \(f_1\). TYPE: Union[float, FloatVectorLike] DEFAULT: 0.05
weight	Relative weight of dataset for fitting. TYPE: float DEFAULT: 1
T	Temperature. Unit = Tunit. TYPE: float DEFAULT: 298.0
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
name	Name of dataset. TYPE: str DEFAULT: ''
source	Source of dataset. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/copodataset.py

class DriftDataset(CopoDataset):
    r"""Dataset of binary monomer drift copolymerization data.

    Container for binary monomer drift copolymerization data, $f_1$ vs $x$.

    Parameters
    ----------
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    x : FloatVectorLike
        Vector of total molar conversion, $x$.
    f1 : FloatVectorLike
        Vector of monomer composition, $f_1$.
    sigma_x: float | FloatVectorLike
        Absolute standard deviation of $x$.
    sigma_f1: float | FloatVectorLike
        Absolute standard deviation of $f_1$.
    weight: float
        Relative weight of dataset for fitting.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    name: str
        Name of dataset.
    source: str
        Source of dataset.
    """

    varmap = {'x': 'x', 'f': 'y'}

    def __init__(self,
                 M1: str,
                 M2: str,
                 x: FloatVectorLike,
                 f1: FloatVectorLike,
                 sigma_x: Union[float, FloatVectorLike] = 5e-2,
                 sigma_f1: Union[float, FloatVectorLike] = 5e-2,
                 weight: float = 1,
                 T: float = 298.,
                 Tunit: Literal['C', 'K'] = 'K',
                 name: str = '',
                 source: str = '',
                 ) -> None:
        """Construct `DriftDataset` with the given parameters."""
        super().__init__(M1, M2, x, f1, sigma_x, sigma_f1, weight, T, Tunit,
                         name, source)

__init__

__init__(
    M1: str,
    M2: str,
    x: FloatVectorLike,
    f1: FloatVectorLike,
    sigma_x: Union[float, FloatVectorLike] = 0.05,
    sigma_f1: Union[float, FloatVectorLike] = 0.05,
    weight: float = 1,
    T: float = 298.0,
    Tunit: Literal["C", "K"] = "K",
    name: str = "",
    source: str = "",
) -> None

Construct DriftDataset with the given parameters.

Source code in src/polykin/copolymerization/copodataset.py

def __init__(self,
             M1: str,
             M2: str,
             x: FloatVectorLike,
             f1: FloatVectorLike,
             sigma_x: Union[float, FloatVectorLike] = 5e-2,
             sigma_f1: Union[float, FloatVectorLike] = 5e-2,
             weight: float = 1,
             T: float = 298.,
             Tunit: Literal['C', 'K'] = 'K',
             name: str = '',
             source: str = '',
             ) -> None:
    """Construct `DriftDataset` with the given parameters."""
    super().__init__(M1, M2, x, f1, sigma_x, sigma_f1, weight, T, Tunit,
                     name, source)

ImplicitPenultimateModel

Implicit penultimate binary copolymerization model.

This model is a special case of the general (explicit) penultimate model, with a smaller number of independent parameters. As in the explicit version, the reactivity of a macroradical depends on the nature of the penultimate and terminal repeating units. A binary system is, thus, described by eight propagation reactions:

\[\begin{matrix} P^{\bullet}_{11} + M_1 \overset{k_{111}}{\rightarrow} P^{\bullet}_{11} \\ P^{\bullet}_{11} + M_2 \overset{k_{112}}{\rightarrow} P^{\bullet}_{12} \\ P^{\bullet}_{12} + M_1 \overset{k_{121}}{\rightarrow} P^{\bullet}_{21} \\ P^{\bullet}_{12} + M_2 \overset{k_{122}}{\rightarrow} P^{\bullet}_{22} \\ P^{\bullet}_{21} + M_1 \overset{k_{211}}{\rightarrow} P^{\bullet}_{11} \\ P^{\bullet}_{21} + M_2 \overset{k_{212}}{\rightarrow} P^{\bullet}_{12} \\ P^{\bullet}_{22} + M_1 \overset{k_{221}}{\rightarrow} P^{\bullet}_{21} \\ P^{\bullet}_{22} + M_2 \overset{k_{222}}{\rightarrow} P^{\bullet}_{22} \\ \end{matrix}\]

where \(k_{iii}\) are the homo-propagation rate coefficients and \(k_{ijk}\) are the cross-propagation coefficients. The six cross-propagation coefficients are specified via just four reactivity ratios, which are divided in two categories. There are two monomer reactivity ratios, which are defined as \(r_1=k_{111}/k_{112}=k_{211}/k_{212}\) and \(r_2=k_{222}/k_{221}=k_{122}/k_{121}\). Additionally, there are two radical reactivity ratios defined as \(s_1=k_{211}/k_{111}\) and \(s_2=k_{122}/k_{222}\). The latter influence the average propagation rate coefficient, but have no effect on the copolymer composition.

PARAMETER	DESCRIPTION
r1	Monomer reactivity ratio, \(r_1=k_{111}/k_{112}=k_{211}/k_{212}\). TYPE: float
r2	Monomer reactivity ratio, \(r_2=k_{222}/k_{221}=k_{122}/k_{121}\). TYPE: float
s1	Radical reactivity ratio, \(s_1=k_{211}/k_{111}\). TYPE: float
s2	Radical reactivity ratio, \(s_2=k_{122}/k_{222}\). TYPE: float
k1	Homopropagation rate coefficient of M1, \(k_1 \equiv k_{111}\). TYPE: Arrhenius | None DEFAULT: None
k2	Homopropagation rate coefficient of M2, \(k_2 \equiv k_{222}\). TYPE: Arrhenius | None DEFAULT: None
M1	Name of M1. TYPE: str DEFAULT: 'M1'
M2	Name of M2. TYPE: str DEFAULT: 'M2'
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/penultimate.py

class ImplicitPenultimateModel(TerminalModel):
    r"""Implicit penultimate binary copolymerization model.

    This model is a special case of the general (explicit) penultimate model,
    with a smaller number of independent parameters. As in the explicit
    version, the reactivity of a macroradical depends on the nature of the
    _penultimate_ and _terminal_ repeating units. A binary system is, thus,
    described by eight propagation reactions:

    \begin{matrix}
    P^{\bullet}_{11} + M_1 \overset{k_{111}}{\rightarrow} P^{\bullet}_{11} \\
    P^{\bullet}_{11} + M_2 \overset{k_{112}}{\rightarrow} P^{\bullet}_{12} \\
    P^{\bullet}_{12} + M_1 \overset{k_{121}}{\rightarrow} P^{\bullet}_{21} \\
    P^{\bullet}_{12} + M_2 \overset{k_{122}}{\rightarrow} P^{\bullet}_{22} \\
    P^{\bullet}_{21} + M_1 \overset{k_{211}}{\rightarrow} P^{\bullet}_{11} \\
    P^{\bullet}_{21} + M_2 \overset{k_{212}}{\rightarrow} P^{\bullet}_{12} \\
    P^{\bullet}_{22} + M_1 \overset{k_{221}}{\rightarrow} P^{\bullet}_{21} \\
    P^{\bullet}_{22} + M_2 \overset{k_{222}}{\rightarrow} P^{\bullet}_{22} \\
    \end{matrix}

    where $k_{iii}$ are the homo-propagation rate coefficients and $k_{ijk}$
    are the cross-propagation coefficients. The six cross-propagation
    coefficients are specified via just four reactivity ratios, which
    are divided in two categories. There are two monomer reactivity
    ratios, which are defined as $r_1=k_{111}/k_{112}=k_{211}/k_{212}$ and
    $r_2=k_{222}/k_{221}=k_{122}/k_{121}$. Additionally, there
    are two radical reactivity ratios defined as $s_1=k_{211}/k_{111}$ and
    $s_2=k_{122}/k_{222}$. The latter influence the average propagation rate
    coefficient, but have no effect on the copolymer composition.

    Parameters
    ----------
    r1 : float
        Monomer reactivity ratio, $r_1=k_{111}/k_{112}=k_{211}/k_{212}$.
    r2 : float
        Monomer reactivity ratio, $r_2=k_{222}/k_{221}=k_{122}/k_{121}$.
    s1 : float
        Radical reactivity ratio, $s_1=k_{211}/k_{111}$.
    s2 : float
        Radical reactivity ratio, $s_2=k_{122}/k_{222}$.
    k1 : Arrhenius | None
        Homopropagation rate coefficient of M1, $k_1 \equiv k_{111}$.
    k2 : Arrhenius | None
        Homopropagation rate coefficient of M2, $k_2 \equiv k_{222}$.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    name : str
        Name.
    """

    r1: float
    r2: float
    s1: float
    s2: float

    _pnames = ('r1', 'r2', 's1', 's2')

    def __init__(self,
                 r1: float,
                 r2: float,
                 s1: float,
                 s2: float,
                 k1: Optional[Arrhenius] = None,
                 k2: Optional[Arrhenius] = None,
                 M1: str = 'M1',
                 M2: str = 'M2',
                 name: str = ''
                 ) -> None:

        check_bounds(s1, 0., np.inf, 's1')
        check_bounds(s2, 0., np.inf, 's2')

        self.s1 = s1
        self.s2 = s2
        super().__init__(r1, r2, k1, k2, M1, M2, name)

    def kii(self,
            f1: Union[float, FloatArray],
            T: float,
            Tunit: Literal['C', 'K'] = 'K',
            ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
        r"""Pseudo-homopropagation rate coefficients.

        In the implicit penultimate model, the pseudohomopropagation rate
        coefficients depend on the instantaneous comonomer composition
        according to:

        \begin{aligned}
        \bar{k}_{11} &= k_{111} \frac{f_1 r_1 + f_2}{f_1 r_1 + f_2/s_1} \\
        \bar{k}_{22} &= k_{222} \frac{f_2 r_2 + f_1}{f_2 r_2 + f_1/s_2}
        \end{aligned}

        where $r_i$ are the monomer reactivity ratios, $s_i$ are the radical
        reactivity ratios, and $k_{iii}$ are the homopropagation rate
        coefficients.

        Parameters
        ----------
        f1 : float | FloatArray
            Molar fraction of M1.
        T : float
            Temperature. Unit = `Tunit`.
        Tunit : Literal['C', 'K']
            Temperature unit.

        Returns
        -------
        tuple[float | FloatArray, float | FloatArray]
            Tuple of pseudohomopropagation rate coefficients,
            ($\bar{k}_{11}$, $\bar{k}_{22}$).
        """
        if self.k1 is None or self.k2 is None:
            raise ValueError(
                "To use this feature, `k1` and `k2` cannot be `None`.")
        f2 = 1. - f1
        r1 = self.r1
        r2 = self.r2
        s1 = self.s1
        s2 = self.s2
        k1 = self.k1(T, Tunit)
        k2 = self.k2(T, Tunit)
        k11 = k1*(f1*r1 + f2)/(f1*r1 + f2/s1)
        k22 = k2*(f2*r2 + f1)/(f2*r2 + f1/s2)
        return (k11, k22)

F1

F1(
    f1: Union[float, FloatArrayLike]
) -> Union[float, FloatArray]

Calculate the instantaneous copolymer composition, \(F_1\).

The calculation is handled by inst_copolymer_binary.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
float | FloatArray	Instantaneous copolymer composition, \(F_1\).

Source code in src/polykin/copolymerization/terminal.py

def F1(self,
       f1: Union[float, FloatArrayLike]
       ) -> Union[float, FloatArray]:
    r"""Calculate the instantaneous copolymer composition, $F_1$.

    The calculation is handled by
    [`inst_copolymer_binary`](inst_copolymer_binary.md).

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    float | FloatArray
        Instantaneous copolymer composition, $F_1$.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return inst_copolymer_binary(f1, *self.ri(f1))

add_data

add_data(
    data: Union[CopoDataset, list[CopoDataset]]
) -> None

Add a copolymerization dataset for subsequent analysis.

PARAMETER	DESCRIPTION
data	Experimental dataset(s). TYPE: Union[CopoDataset, list[CopoDataset]]

Source code in src/polykin/copolymerization/terminal.py

def add_data(self,
             data: Union[CopoDataset, list[CopoDataset]]
             ) -> None:
    r"""Add a copolymerization dataset for subsequent analysis.

    Parameters
    ----------
    data : Union[CopoDataset, list[CopoDataset]]
        Experimental dataset(s).
    """
    if not isinstance(data, list):
        data = [data]

    valid_monomers = {self.M1, self.M2}
    for ds in data:
        ds_monomers = {ds.M1, ds.M2}
        if ds_monomers != valid_monomers:
            raise ValueError(
                f"Monomers defined in dataset `{ds.name}` are invalid: "
                f"{ds_monomers}!={valid_monomers}")
        if ds not in set(self.data):
            self.data.append(ds)
        else:
            print(f"Warning: duplicate dataset '{ds.name}' was skipped.")

azeotrope property

azeotrope: Optional[float]

Calculate the azeotrope composition.

An azeotrope (i.e., a point where \(F_1=f_1\)) only exists if both reactivity ratios are smaller than unity. In that case, the azeotrope composition is given by:

\[ f_{1,azeo} = \frac{1 - r_2}{(2 - r_1 - r_2)} \]

where \(r_1\) and \(r_2\) are the reactivity ratios.

RETURNS	DESCRIPTION
float | None	If an azeotrope exists, it returns its composition in terms of \(f_1\).

drift

drift(
    f10: Union[float, FloatVectorLike],
    x: Union[float, FloatVectorLike],
) -> FloatArray

Calculate drift of comonomer composition in a closed system for a given total monomer conversion.

In a closed binary system, the drift in monomer composition is given by the solution of the following differential equation:

\[ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} \]

with initial condition \(f_1(0)=f_{1,0}\), where \(f_1\) and \(F_1\) are, respectively, the instantaneous comonomer and copolymer composition of M1, and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f10	Initial molar fraction of M1, \(f_{1,0}=f_1(0)\). TYPE: float | FloatVectorLike
x	Value(s) of total monomer conversion values where the drift is to be evaluated. TYPE: float | FloatVectorLike

RETURNS	DESCRIPTION
FloatArray	Monomer fraction of M1 at a given conversion, \(f_1(x)\).

Source code in src/polykin/copolymerization/terminal.py

def drift(self,
          f10: Union[float, FloatVectorLike],
          x: Union[float, FloatVectorLike]
          ) -> FloatArray:
    r"""Calculate drift of comonomer composition in a closed system for a
    given total monomer conversion.

    In a closed binary system, the drift in monomer composition is given by
    the solution of the following differential equation:

    $$ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} $$

    with initial condition $f_1(0)=f_{1,0}$, where $f_1$ and $F_1$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    M1, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f10 : float | FloatVectorLike
        Initial molar fraction of M1, $f_{1,0}=f_1(0)$.
    x : float | FloatVectorLike
        Value(s) of total monomer conversion values where the drift is to
        be evaluated.

    Returns
    -------
    FloatArray
        Monomer fraction of M1 at a given conversion, $f_1(x)$.
    """
    f10, x = convert_FloatOrVectorLike_to_FloatVector([f10, x], False)

    def df1dx(x_, f1_):
        F1_ = inst_copolymer_binary(f1_, *self.ri(f1_))
        return (f1_ - F1_)/(1. - x_ + eps)

    sol = solve_ivp(df1dx,
                    (0., max(x)),
                    f10,
                    t_eval=x,
                    method='LSODA',
                    vectorized=True,
                    atol=1e-4,
                    rtol=1e-4)
    if sol.success:
        result = sol.y
        result = np.maximum(0., result)
        if result.shape[0] == 1:
            result = result[0]
    else:
        result = np.empty_like(x)
        result[:] = np.nan
        print(sol.message)

    return result

from_Qe classmethod

from_Qe(
    Qe1: tuple[float, float],
    Qe2: tuple[float, float],
    k1: Optional[Arrhenius] = None,
    k2: Optional[Arrhenius] = None,
    M1: str = "M1",
    M2: str = "M2",
    name: str = "",
) -> TerminalModel

Construct TerminalModel from Q-e values.

Alternative constructor that takes the \(Q\)-\(e\) values of the monomers as primary input instead of the reactivity ratios.

The conversion from Q-e to r is handled by convert_Qe_to_r.

PARAMETER	DESCRIPTION
Qe1	Q-e values of M1. TYPE: tuple[float, float]
Qe2	Q-e values of M2. TYPE: tuple[float, float]
k1	Homopropagation rate coefficient of M1, \(k_1 \equiv k_{11}\). TYPE: Arrhenius | None DEFAULT: None
k2	Homopropagation rate coefficient of M2, \(k_2 \equiv k_{22}\). TYPE: Arrhenius | None DEFAULT: None
M1	Name of M1. TYPE: str DEFAULT: 'M1'
M2	Name of M2. TYPE: str DEFAULT: 'M2'
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/terminal.py

@classmethod
def from_Qe(cls,
            Qe1: tuple[float, float],
            Qe2: tuple[float, float],
            k1: Optional[Arrhenius] = None,
            k2: Optional[Arrhenius] = None,
            M1: str = 'M1',
            M2: str = 'M2',
            name: str = ''
            ) -> TerminalModel:
    r"""Construct `TerminalModel` from Q-e values.

    Alternative constructor that takes the $Q$-$e$ values of the monomers
    as primary input instead of the reactivity ratios.

    The conversion from Q-e to r is handled by
    [`convert_Qe_to_r`](convert_Qe_to_r.md).

    Parameters
    ----------
    Qe1 : tuple[float, float]
        Q-e values of M1.
    Qe2 : tuple[float, float]
        Q-e values of M2.
    k1 : Arrhenius | None
        Homopropagation rate coefficient of M1, $k_1 \equiv k_{11}$.
    k2 : Arrhenius | None
        Homopropagation rate coefficient of M2, $k_2 \equiv k_{22}$.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    name : str
        Name.
    """
    r = convert_Qe_to_r([Qe1, Qe2])
    r1 = r[0, 1]
    r2 = r[1, 0]
    return cls(r1, r2, k1, k2, M1, M2, name)

kii

kii(
    f1: Union[float, FloatArray],
    T: float,
    Tunit: Literal["C", "K"] = "K",
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Pseudo-homopropagation rate coefficients.

In the implicit penultimate model, the pseudohomopropagation rate coefficients depend on the instantaneous comonomer composition according to:

\[\begin{aligned} \bar{k}_{11} &= k_{111} \frac{f_1 r_1 + f_2}{f_1 r_1 + f_2/s_1} \\ \bar{k}_{22} &= k_{222} \frac{f_2 r_2 + f_1}{f_2 r_2 + f_1/s_2} \end{aligned}\]

where \(r_i\) are the monomer reactivity ratios, \(s_i\) are the radical reactivity ratios, and \(k_{iii}\) are the homopropagation rate coefficients.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArray
T	Temperature. Unit = Tunit. TYPE: float
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'

RETURNS	DESCRIPTION
tuple[float | FloatArray, float | FloatArray]	Tuple of pseudohomopropagation rate coefficients, (\(\bar{k}_{11}\), \(\bar{k}_{22}\)).

Source code in src/polykin/copolymerization/penultimate.py

def kii(self,
        f1: Union[float, FloatArray],
        T: float,
        Tunit: Literal['C', 'K'] = 'K',
        ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    r"""Pseudo-homopropagation rate coefficients.

    In the implicit penultimate model, the pseudohomopropagation rate
    coefficients depend on the instantaneous comonomer composition
    according to:

    \begin{aligned}
    \bar{k}_{11} &= k_{111} \frac{f_1 r_1 + f_2}{f_1 r_1 + f_2/s_1} \\
    \bar{k}_{22} &= k_{222} \frac{f_2 r_2 + f_1}{f_2 r_2 + f_1/s_2}
    \end{aligned}

    where $r_i$ are the monomer reactivity ratios, $s_i$ are the radical
    reactivity ratios, and $k_{iii}$ are the homopropagation rate
    coefficients.

    Parameters
    ----------
    f1 : float | FloatArray
        Molar fraction of M1.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.

    Returns
    -------
    tuple[float | FloatArray, float | FloatArray]
        Tuple of pseudohomopropagation rate coefficients,
        ($\bar{k}_{11}$, $\bar{k}_{22}$).
    """
    if self.k1 is None or self.k2 is None:
        raise ValueError(
            "To use this feature, `k1` and `k2` cannot be `None`.")
    f2 = 1. - f1
    r1 = self.r1
    r2 = self.r2
    s1 = self.s1
    s2 = self.s2
    k1 = self.k1(T, Tunit)
    k2 = self.k2(T, Tunit)
    k11 = k1*(f1*r1 + f2)/(f1*r1 + f2/s1)
    k22 = k2*(f2*r2 + f1)/(f2*r2 + f1/s2)
    return (k11, k22)

kp

kp(
    f1: Union[float, FloatArrayLike],
    T: float,
    Tunit: Literal["C", "K"] = "K",
) -> Union[float, FloatArray]

Calculate the average propagation rate coefficient, \(\bar{k}_p\).

The calculation is handled by kp_average_binary.

Note

This feature requires the attributes k11 and k22 to be defined.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
T	Temperature. Unit = Tunit. TYPE: float
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'

RETURNS	DESCRIPTION
float | FloatArray	Average propagation rate coefficient. Unit = L/(mol·s)

Source code in src/polykin/copolymerization/terminal.py

def kp(self,
       f1: Union[float, FloatArrayLike],
       T: float,
       Tunit: Literal['C', 'K'] = 'K'
       ) -> Union[float, FloatArray]:
    r"""Calculate the average propagation rate coefficient, $\bar{k}_p$.

    The calculation is handled by
    [`kp_average_binary`](kp_average_binary.md).

    Note
    ----
    This feature requires the attributes `k11` and `k22` to be defined.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.

    Returns
    -------
    float | FloatArray
        Average propagation rate coefficient. Unit = L/(mol·s)
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return kp_average_binary(f1, *self.ri(f1), *self.kii(f1, T, Tunit))

plot

plot(
    kind: Literal["drift", "kp", "Mayo", "triads"],
    show: Literal["auto", "all", "data", "model"] = "auto",
    M: Literal[1, 2] = 1,
    f0: Optional[Union[float, FloatVectorLike]] = None,
    T: Optional[float] = None,
    Tunit: Literal["C", "K"] = "K",
    title: Optional[str] = None,
    axes: Optional[Axes] = None,
    return_objects: bool = False,
) -> Optional[tuple[Optional[Figure], Axes]]

Generate a plot of instantaneous copolymer composition, monomer composition drift, or average propagation rate coefficient.

PARAMETER	DESCRIPTION
kind	Kind of plot to be generated. TYPE: Literal['drift', 'kp', 'Mayo', 'triads']
show	What informatation is to be plotted. TYPE: Literal['auto', 'all', 'data', 'model'] DEFAULT: 'auto'
M	Index of the monomer to be used in input argument f0 and in output results. Specifically, if M=i, then f0 stands for \(f_i(0)\) and plots will be generated in terms of \(f_i\) and \(F_i\). TYPE: Literal[1, 2] DEFAULT: 1
f0	Initial monomer composition, \(f_i(0)\), as required for a monomer composition drift plot. TYPE: float | FloatVectorLike | None DEFAULT: None
T	Temperature. Unit = Tunit. TYPE: float | None DEFAULT: None
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
title	Title of plot. If None, a default title with the monomer names will be used. TYPE: str | None DEFAULT: None
axes	Matplotlib Axes object. TYPE: Axes | None DEFAULT: None
return_objects	If True, the Figure and Axes objects are returned (for saving or further manipulations). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
tuple[Figure | None, Axes] | None	Figure and Axes objects if return_objects is True.

Source code in src/polykin/copolymerization/terminal.py

def plot(self,
         kind: Literal['drift', 'kp', 'Mayo', 'triads'],
         show: Literal['auto', 'all', 'data', 'model'] = 'auto',
         M: Literal[1, 2] = 1,
         f0: Optional[Union[float, FloatVectorLike]] = None,
         T: Optional[float] = None,
         Tunit: Literal['C', 'K'] = 'K',
         title: Optional[str] = None,
         axes: Optional[Axes] = None,
         return_objects: bool = False
         ) -> Optional[tuple[Optional[Figure], Axes]]:
    r"""Generate a plot of instantaneous copolymer composition, monomer
    composition drift, or average propagation rate coefficient.

    Parameters
    ----------
    kind : Literal['drift', 'kp', 'Mayo', 'triads']
        Kind of plot to be generated.
    show : Literal['auto', 'all', 'data', 'model']
        What informatation is to be plotted.
    M : Literal[1, 2]
        Index of the monomer to be used in input argument `f0` and in
        output results. Specifically, if `M=i`, then `f0` stands for
        $f_i(0)$ and plots will be generated in terms of $f_i$ and $F_i$.
    f0 : float | FloatVectorLike | None
        Initial monomer composition, $f_i(0)$, as required for a monomer
        composition drift plot.
    T : float | None
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    title : str | None
        Title of plot. If `None`, a default title with the monomer names
        will be used.
    axes : Axes | None
        Matplotlib Axes object.
    return_objects : bool
        If `True`, the Figure and Axes objects are returned (for saving or
        further manipulations).

    Returns
    -------
    tuple[Figure | None, Axes] | None
        Figure and Axes objects if return_objects is `True`.
    """
    check_in_set(M, {1, 2}, 'M')
    check_in_set(kind, {'Mayo', 'kp', 'drift', 'triads'}, 'kind')
    check_in_set(show, {'auto', 'all', 'data', 'model'}, 'show')

    label_model = None
    if axes is None:
        fig, ax = plt.subplots()
        if title is None:
            titles = {'Mayo': "Mayo-Lewis diagram",
                      'drift': "Monomer composition drift",
                      'kp': "Average propagation coefficient",
                      'triads': "Triad fractions"}
            title = titles[kind] + f" {self.M1}(1)-{self.M2}(2)"
        if title:
            fig.suptitle(title)
    else:
        ax = axes
        fig = None
        if self.name:
            label_model = self.name

    unit_range = (0., 1.)
    npoints = 1000
    Mname = self.M1 if M == 1 else self.M2
    ndataseries = 0

    if show == 'auto':
        if not self.data:
            show = 'model'
        else:
            show = 'all'

    if kind == 'Mayo':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(fr"$F_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        ax.plot(unit_range, unit_range, color='black', linewidth=0.5)

        if show in {'model', 'all'}:
            x = np.linspace(*unit_range, npoints, dtype=np.float64)
            y = self.F1(x)
            if M == 2:
                x[:] = 1. - x
                y[:] = 1. - y  # type: ignore
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:
            for ds in self.data:
                if not isinstance(ds, MayoDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('F', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

    elif kind == 'drift':

        ax.set_xlabel(r"Total molar monomer conversion, $x$")
        ax.set_ylabel(fr"$f_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        if show == 'model':

            if f0 is None:
                raise ValueError("`f0` is required for a `drift` plot.")
            else:
                if isinstance(f0, (int, float)):
                    f0 = [f0]
                f0 = np.array(f0, dtype=np.float64)
                check_bounds(f0, *unit_range, 'f0')

            x = np.linspace(*unit_range, 1000)
            if M == 2:
                f0[:] = 1. - f0
            y = self.drift(f0, x)
            if M == 2:
                y[:] = 1. - y
            if y.ndim == 1:
                y = y[np.newaxis, :]
            for i in range(y.shape[0]):
                ax.plot(x, y[i, :], label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, DriftDataset):
                    continue
                x = ds.getvar('x')
                y = ds.getvar('f', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.drift(ds.getvar('f', self.M1)[0], x)
                    if M == 2:
                        y[:] = 1. - y
                    ax.plot(x, y, label=label_model)

    elif kind == 'kp':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(r"$\bar{k}_p$")
        ax.set_xlim(*unit_range)

        if show == 'model':

            if T is None:
                raise ValueError("`T` is required for a `kp` plot.")

            x = np.linspace(*unit_range, npoints)
            y = self.kp(x, T, Tunit)
            if M == 2:
                x[:] = 1. - x
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, kpDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('kp')
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)
                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.kp(x, ds.T, ds.Tunit)
                    if M == 2:
                        x[:] = 1. - x
                    ax.plot(x, y, label=label_model)

    ax.grid(True)

    if axes is not None or ndataseries:
        ax.legend(bbox_to_anchor=(1.05, 1.), loc="upper left")

    if return_objects:
        return (fig, ax)

ri

ri(
    _,
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Return the evaluated reactivity ratios at the given conditions.

Source code in src/polykin/copolymerization/terminal.py

def ri(self,
       _
       ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    return (self.r1, self.r2)

sequence

sequence(
    f1: Union[float, FloatArrayLike],
    k: Optional[Union[int, IntArrayLike]] = None,
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous sequence length probability or the number-average sequence length.

For a binary system, the probability of finding \(k\) consecutive units of monomer \(i\) in a chain is:

\[ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} \]

and the corresponding number-average sequence length is:

\[ \bar{S}_i = \sum_k k S_{i,k} = \frac{1}{1 - P_{ii}} \]

where \(P_{ii}\) is the transition probability \(i \rightarrow i\), which is a function of the monomer composition and the reactivity ratios.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 177.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
k	Sequence length, i.e., number of consecutive units in a chain. If None, the number-average sequence length will be computed. TYPE: int | IntArrayLike | None DEFAULT: None

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	If k is None, the number-average sequence lengths, {'1': \(\bar{S}_1\), '2': \(\bar{S}_2\)}. Otherwise, the sequence probabilities, {'1': \(S_{1,k}\), '2': \(S_{2,k}\)}.

Source code in src/polykin/copolymerization/terminal.py

def sequence(self,
             f1: Union[float, FloatArrayLike],
             k: Optional[Union[int, IntArrayLike]] = None,
             ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous sequence length probability or the
    number-average sequence length.

    For a binary system, the probability of finding $k$ consecutive units
    of monomer $i$ in a chain is:

    $$ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} $$

    and the corresponding number-average sequence length is:

    $$ \bar{S}_i = \sum_k k S_{i,k} = \frac{1}{1 - P_{ii}} $$

    where $P_{ii}$ is the transition probability $i \rightarrow i$, which
    is a function of the monomer composition and the reactivity ratios.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 177.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    k : int | IntArrayLike | None
        Sequence length, i.e., number of consecutive units in a chain.
        If `None`, the number-average sequence length will be computed.

    Returns
    -------
    dict[str, float | FloatArray]
        If `k is None`, the number-average sequence lengths,
        {'1': $\bar{S}_1$, '2': $\bar{S}_2$}. Otherwise, the
        sequence probabilities, {'1': $S_{1,k}$, '2': $S_{2,k}$}.
    """

    P11, _, _, P22 = self.transitions(f1).values()

    if k is None:
        result = {str(i + 1): 1/(1 - P + eps)
                  for i, P in enumerate([P11, P22])}
    else:
        if isinstance(k, (list, tuple)):
            k = np.array(k, dtype=np.int32)
        result = {str(i + 1): (1. - P)*P**(k - 1)
                  for i, P in enumerate([P11, P22])}

    return result

transitions

transitions(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous transition probabilities.

For a binary system, the transition probabilities are given by:

\[\begin{aligned} P_{ii} &= \frac{r_i f_i}{r_i f_i + (1 - f_i)} \\ P_{ij} &= 1 - P_{ii} \end{aligned}\]

where \(i,j=1,2, i \neq j\).

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 178.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Transition probabilities, {'11': \(P_{11}\), '12': \(P_{12}\), ... }.

Source code in src/polykin/copolymerization/terminal.py

def transitions(self,
                f1: Union[float, FloatArrayLike]
                ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous transition probabilities.

    For a binary system, the transition probabilities are given by:

    \begin{aligned}
        P_{ii} &= \frac{r_i f_i}{r_i f_i + (1 - f_i)} \\
        P_{ij} &= 1 - P_{ii}
    \end{aligned}

    where $i,j=1,2, i \neq j$.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 178.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Transition probabilities, {'11': $P_{11}$, '12': $P_{12}$, ... }.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')

    f2 = 1. - f1
    r1 = self.r1
    r2 = self.r2
    P11 = r1*f1/(r1*f1 + f2)
    P22 = r2*f2/(r2*f2 + f1)
    P12 = 1. - P11
    P21 = 1. - P22

    result = {'11': P11,
              '12': P12,
              '21': P21,
              '22': P22}

    return result

triads

triads(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous triad fractions.

For a binary system, the triad fractions are given by:

\[\begin{aligned} A_{iii} &= F_i P_{ii}^2 \\ A_{iij} &= 2 F_i P_{ii} P_{ij} \\ A_{jij} &= F_i P_{ij}^2 \end{aligned}\]

where \(P_{ij}\) is the transition probability \(i \rightarrow j\) and \(F_i\) is the instantaneous copolymer composition.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 179.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Triad fractions, {'111': \(A_{111}\), '112': \(A_{112}\), '212': \(A_{212}\), ... }.

Source code in src/polykin/copolymerization/terminal.py

def triads(self,
           f1: Union[float, FloatArrayLike]
           ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous triad fractions.

    For a binary system, the triad fractions are given by:

    \begin{aligned}
        A_{iii} &= F_i P_{ii}^2 \\
        A_{iij} &= 2 F_i P_{ii} P_{ij} \\
        A_{jij} &= F_i P_{ij}^2
    \end{aligned}

    where $P_{ij}$ is the transition probability $i \rightarrow j$ and
    $F_i$ is the instantaneous copolymer composition.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 179.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Triad fractions,
        {'111': $A_{111}$, '112': $A_{112}$, '212': $A_{212}$, ... }.
    """

    P11, P12, P21, P22 = self.transitions(f1).values()

    F1 = P21/(P12 + P21)
    F2 = 1. - F1

    A111 = F1*P11**2
    A112 = F1*2*P11*P12
    A212 = F1*P12**2

    A222 = F2*P22**2
    A221 = F2*2*P22*P21
    A121 = F2*P21**2

    result = {'111': A111,
              '112': A112,
              '212': A212,
              '222': A222,
              '221': A221,
              '121': A121}

    return result

MayoDataset

Dataset of binary instantaneous copolymerization data.

Container for binary instantaneous copolymerization data, \(F_1\) vs \(f_1\), as usually obtained from low-conversion experiments.

PARAMETER	DESCRIPTION
M1	Name of M1. TYPE: str
M2	Name of M2. TYPE: str
f1	Vector of monomer composition, \(f_1\). TYPE: FloatVectorLike
F1	Vector of instantaneous copolymer composition, \(F_1\). TYPE: FloatVectorLike
sigma_f1	Absolute standard deviation of \(f_1\) (\(\sigma_{f_1} \equiv \sigma_{f_2}\)). TYPE: Union[float, FloatVectorLike] DEFAULT: 0.01
sigma_F1	Absolute standard deviation of \(F_1\) (\(\sigma_{F_1} \equiv \sigma_{F_2}\)). TYPE: Union[float, FloatVectorLike] DEFAULT: 0.05
weight	Relative weight of dataset for fitting. TYPE: float DEFAULT: 1
T	Temperature. Unit = Tunit. TYPE: float DEFAULT: 298.0
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
name	Name of dataset. TYPE: str DEFAULT: ''
source	Source of dataset. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/copodataset.py

class MayoDataset(CopoDataset):
    r"""Dataset of binary instantaneous copolymerization data.

    Container for binary instantaneous copolymerization data, $F_1$ vs $f_1$,
    as usually obtained from low-conversion experiments.

    Parameters
    ----------
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    f1 : FloatVectorLike
        Vector of monomer composition, $f_1$.
    F1 : FloatVectorLike
        Vector of instantaneous copolymer composition, $F_1$.
    sigma_f1: float | FloatVectorLike
        Absolute standard deviation of $f_1$
        ($\sigma_{f_1} \equiv \sigma_{f_2}$).
    sigma_F1: float | FloatVectorLike
        Absolute standard deviation of $F_1$
        ($\sigma_{F_1} \equiv \sigma_{F_2}$).
    weight: float
        Relative weight of dataset for fitting.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    name: str
        Name of dataset.
    source: str
        Source of dataset.
    """
    varmap = {'f': 'x', 'F': 'y'}

    def __init__(self,
                 M1: str,
                 M2: str,
                 f1: FloatVectorLike,
                 F1: FloatVectorLike,
                 sigma_f1: Union[float, FloatVectorLike] = 1e-2,
                 sigma_F1: Union[float, FloatVectorLike] = 5e-2,
                 weight: float = 1,
                 T: float = 298.,
                 Tunit: Literal['C', 'K'] = 'K',
                 name: str = '',
                 source: str = '',
                 ) -> None:
        """Construct `MayoDataset` with the given parameters."""
        super().__init__(M1, M2, f1, F1, sigma_f1, sigma_F1, weight, T, Tunit,
                         name, source)

__init__

__init__(
    M1: str,
    M2: str,
    f1: FloatVectorLike,
    F1: FloatVectorLike,
    sigma_f1: Union[float, FloatVectorLike] = 0.01,
    sigma_F1: Union[float, FloatVectorLike] = 0.05,
    weight: float = 1,
    T: float = 298.0,
    Tunit: Literal["C", "K"] = "K",
    name: str = "",
    source: str = "",
) -> None

Construct MayoDataset with the given parameters.

Source code in src/polykin/copolymerization/copodataset.py

def __init__(self,
             M1: str,
             M2: str,
             f1: FloatVectorLike,
             F1: FloatVectorLike,
             sigma_f1: Union[float, FloatVectorLike] = 1e-2,
             sigma_F1: Union[float, FloatVectorLike] = 5e-2,
             weight: float = 1,
             T: float = 298.,
             Tunit: Literal['C', 'K'] = 'K',
             name: str = '',
             source: str = '',
             ) -> None:
    """Construct `MayoDataset` with the given parameters."""
    super().__init__(M1, M2, f1, F1, sigma_f1, sigma_F1, weight, T, Tunit,
                     name, source)

PenultimateModel

Penultimate binary copolymerization model.

According to this model, the reactivity of a macroradical depends on the nature of the penultimate and terminal repeating units. A binary system is, thus, described by eight propagation reactions:

\[\begin{matrix} P^{\bullet}_{11} + M_1 \overset{k_{111}}{\rightarrow} P^{\bullet}_{11} \\ P^{\bullet}_{11} + M_2 \overset{k_{112}}{\rightarrow} P^{\bullet}_{12} \\ P^{\bullet}_{12} + M_1 \overset{k_{121}}{\rightarrow} P^{\bullet}_{21} \\ P^{\bullet}_{12} + M_2 \overset{k_{122}}{\rightarrow} P^{\bullet}_{22} \\ P^{\bullet}_{21} + M_1 \overset{k_{211}}{\rightarrow} P^{\bullet}_{11} \\ P^{\bullet}_{21} + M_2 \overset{k_{212}}{\rightarrow} P^{\bullet}_{12} \\ P^{\bullet}_{22} + M_1 \overset{k_{221}}{\rightarrow} P^{\bullet}_{21} \\ P^{\bullet}_{22} + M_2 \overset{k_{222}}{\rightarrow} P^{\bullet}_{22} \\ \end{matrix}\]

where \(k_{iii}\) are the homo-propagation rate coefficients and \(k_{ijk}\) are the cross-propagation coefficients. The six cross-propagation coefficients are specified via an equal number of reactivity ratios, which are divided in two categories. There are four monomer reactivity ratios, defined as \(r_{11}=k_{111}/k_{112}\), \(r_{12}=k_{122}/k_{121}\), \(r_{21}=k_{211}/k_{212}\) and \(r_{22}=k_{222}/k_{221}\). Additionally, there are two radical reactivity ratios defined as \(s_1=k_{211}/k_{111}\) and \(s_2=k_{122}/k_{222}\). The latter influence the average propagation rate coefficient, but have no effect on the copolymer composition.

PARAMETER	DESCRIPTION
r11	Monomer reactivity ratio, \(r_{11}=k_{111}/k_{112}\). TYPE: float
r12	Monomer reactivity ratio, \(r_{12}=k_{122}/k_{121}\). TYPE: float
r21	Monomer reactivity ratio, \(r_{21}=k_{211}/k_{212}\). TYPE: float
r22	Monomer reactivity ratio, \(r_{22}=k_{222}/k_{221}\). TYPE: float
s1	Radical reactivity ratio, \(s_1=k_{211}/k_{111}\). TYPE: float
s2	Radical reactivity ratio, \(s_2=k_{122}/k_{222}\). TYPE: float
k1	Homopropagation rate coefficient of M1, \(k_1 \equiv k_{111}\). TYPE: Arrhenius | None DEFAULT: None
k2	Homopropagation rate coefficient of M2, \(k_2 \equiv k_{222}\). TYPE: Arrhenius | None DEFAULT: None
M1	Name of M1. TYPE: str DEFAULT: 'M1'
M2	Name of M2. TYPE: str DEFAULT: 'M2'
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/penultimate.py

class PenultimateModel(CopoModel):
    r"""Penultimate binary copolymerization model.

    According to this model, the reactivity of a macroradical depends on the
    nature of the _penultimate_ and _terminal_ repeating units. A binary system
    is, thus, described by eight propagation reactions:

    \begin{matrix}
    P^{\bullet}_{11} + M_1 \overset{k_{111}}{\rightarrow} P^{\bullet}_{11} \\
    P^{\bullet}_{11} + M_2 \overset{k_{112}}{\rightarrow} P^{\bullet}_{12} \\
    P^{\bullet}_{12} + M_1 \overset{k_{121}}{\rightarrow} P^{\bullet}_{21} \\
    P^{\bullet}_{12} + M_2 \overset{k_{122}}{\rightarrow} P^{\bullet}_{22} \\
    P^{\bullet}_{21} + M_1 \overset{k_{211}}{\rightarrow} P^{\bullet}_{11} \\
    P^{\bullet}_{21} + M_2 \overset{k_{212}}{\rightarrow} P^{\bullet}_{12} \\
    P^{\bullet}_{22} + M_1 \overset{k_{221}}{\rightarrow} P^{\bullet}_{21} \\
    P^{\bullet}_{22} + M_2 \overset{k_{222}}{\rightarrow} P^{\bullet}_{22} \\
    \end{matrix}

    where $k_{iii}$ are the homo-propagation rate coefficients and $k_{ijk}$
    are the cross-propagation coefficients. The six cross-propagation
    coefficients are specified via an equal number of reactivity ratios, which
    are divided in two categories. There are four monomer reactivity ratios,
    defined as $r_{11}=k_{111}/k_{112}$, $r_{12}=k_{122}/k_{121}$,
    $r_{21}=k_{211}/k_{212}$ and $r_{22}=k_{222}/k_{221}$. Additionally, there
    are two radical reactivity ratios defined as $s_1=k_{211}/k_{111}$ and
    $s_2=k_{122}/k_{222}$. The latter influence the average propagation rate
    coefficient, but have no effect on the copolymer composition.

    Parameters
    ----------
    r11 : float
        Monomer reactivity ratio, $r_{11}=k_{111}/k_{112}$.
    r12 : float
        Monomer reactivity ratio, $r_{12}=k_{122}/k_{121}$.
    r21 : float
        Monomer reactivity ratio, $r_{21}=k_{211}/k_{212}$.
    r22 : float
        Monomer reactivity ratio, $r_{22}=k_{222}/k_{221}$.
    s1 : float
        Radical reactivity ratio, $s_1=k_{211}/k_{111}$.
    s2 : float
        Radical reactivity ratio, $s_2=k_{122}/k_{222}$.
    k1 : Arrhenius | None
        Homopropagation rate coefficient of M1, $k_1 \equiv k_{111}$.
    k2 : Arrhenius | None
        Homopropagation rate coefficient of M2, $k_2 \equiv k_{222}$.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    name : str
        Name.
    """

    r11: float
    r12: float
    r21: float
    r22: float
    s1: float
    s2: float

    _pnames = ('r11', 'r12', 'r21', 'r22', 's1', 's2')

    def __init__(self,
                 r11: float,
                 r12: float,
                 r21: float,
                 r22: float,
                 s1: float,
                 s2: float,
                 k1: Optional[Arrhenius] = None,
                 k2: Optional[Arrhenius] = None,
                 M1: str = 'M1',
                 M2: str = 'M2',
                 name: str = ''
                 ) -> None:

        check_bounds(r11, 0., np.inf, 'r11')
        check_bounds(r12, 0., np.inf, 'r12')
        check_bounds(r21, 0., np.inf, 'r21')
        check_bounds(r22, 0., np.inf, 'r22')
        check_bounds(s1, 0., np.inf, 's1')
        check_bounds(s2, 0., np.inf, 's2')

        self.r11 = r11
        self.r12 = r12
        self.r21 = r21
        self.r22 = r22
        self.s1 = s1
        self.s2 = s2
        super().__init__(k1, k2, M1, M2, name)

    def ri(self,
            f1: Union[float, FloatArray]
           ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
        r"""Pseudoreactivity ratios.

        In the penultimate model, the pseudoreactivity ratios depend on the
        instantaneous comonomer composition according to:

        \begin{aligned}
           \bar{r}_1 &= r_{21}\frac{f_1 r_{11} + f_2}{f_1 r_{21} + f_2} \\
           \bar{r}_2 &= r_{12}\frac{f_2 r_{22} + f_1}{f_2 r_{12} + f_1}
        \end{aligned}

        where $r_{ij}$ are the monomer reactivity ratios.

        Parameters
        ----------
        f1 : float | FloatArray
            Molar fraction of M1.

        Returns
        -------
        tuple[FloatOrArray, FloatOrArray]
            Pseudoreactivity ratios, ($\bar{r}_1$, $\bar{r}_2$).
        """
        f2 = 1. - f1
        r11 = self.r11
        r12 = self.r12
        r21 = self.r21
        r22 = self.r22
        r1 = r21*(f1*r11 + f2)/(f1*r21 + f2)
        r2 = r12*(f2*r22 + f1)/(f2*r12 + f1)
        return (r1, r2)

    def kii(self,
            f1: Union[float, FloatArray],
            T: float,
            Tunit: Literal['C', 'K'] = 'K',
            ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
        r"""Pseudohomopropagation rate coefficients.

        In the penultimate model, the pseudohomopropagation rate coefficients
        depend on the instantaneous comonomer composition according to:

        \begin{aligned}
        \bar{k}_{11} &= k_{111}\frac{f_1 r_{11} + f_2}{f_1 r_{11} + f_2/s_1} \\
        \bar{k}_{22} &= k_{222}\frac{f_2 r_{22} + f_1}{f_2 r_{22} + f_1/s_2}
        \end{aligned}

        where $r_{ij}$ are the monomer reactivity ratios, $s_i$ are the radical
        reactivity ratios, and $k_{iii}$ are the homopropagation rate
        coefficients.

        Parameters
        ----------
        f1 : float | FloatArray
            Molar fraction of M1.
        T : float
            Temperature. Unit = `Tunit`.
        Tunit : Literal['C', 'K']
            Temperature unit.

        Returns
        -------
        tuple[float | FloatArray, float | FloatArray]
            Pseudohomopropagation rate coefficients,
            ($\bar{k}_{11}$, $\bar{k}_{22}$).
        """
        if self.k1 is None or self.k2 is None:
            raise ValueError(
                "To use this feature, `k1` and `k2` cannot be `None`.")
        f2 = 1. - f1
        r11 = self.r11
        r22 = self.r22
        s1 = self.s1
        s2 = self.s2
        k1 = self.k1(T, Tunit)
        k2 = self.k2(T, Tunit)
        k11 = k1*(f1*r11 + f2)/(f1*r11 + f2/s1)
        k22 = k2*(f2*r22 + f1)/(f2*r22 + f1/s2)
        return (k11, k22)

    def transitions(self,
                    f1: Union[float, FloatArrayLike]
                    ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous transition probabilities.

        For a binary system, the transition probabilities are given by:

        \begin{aligned}
            P_{iii} &= \frac{r_{ii} f_i}{r_{ii} f_i + (1 - f_i)} \\
            P_{jii} &= \frac{r_{ji} f_i}{r_{ji} f_i + (1 - f_i)} \\
            P_{iij} &= 1 -  P_{iii} \\
            P_{jij} &= 1 -  P_{jii}
        \end{aligned}

        where $i,j=1,2, i \neq j$.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 181.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.

        Returns
        -------
        dict[str, float | FloatArray]
            Transition probabilities,
            {'111': $P_{111}$, '211': $P_{211}$, '121': $P_{121}$, ... }.
        """

        if isinstance(f1, (list, tuple)):
            f1 = np.array(f1, dtype=np.float64)
        check_bounds(f1, 0., 1., 'f1')

        f2 = 1. - f1
        r11 = self.r11
        r12 = self.r12
        r21 = self.r21
        r22 = self.r22

        P111 = r11*f1/(r11*f1 + f2)
        P211 = r21*f1/(r21*f1 + f2)
        P112 = 1. - P111
        P212 = 1. - P211

        P222 = r22*f2/(r22*f2 + f1)
        P122 = r12*f2/(r12*f2 + f1)
        P221 = 1. - P222
        P121 = 1. - P122

        result = {'111': P111,
                  '211': P211,
                  '112': P112,
                  '212': P212,
                  '222': P222,
                  '122': P122,
                  '221': P221,
                  '121': P121}

        return result

    def triads(self,
               f1: Union[float, FloatArrayLike],
               ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous triad fractions.

        For a binary system, the triad fractions are given by:

        \begin{aligned}
            F_{iii} &\propto  P_{jii} \frac{P_{iii}}{1 - P_{iii}} \\
            F_{iij} &\propto 2 P_{jii} \\
            F_{jij} &\propto 1 - P_{jii}
        \end{aligned}

        where $P_{ijk}$ is the transition probability
        $i \rightarrow j \rightarrow k$, which is a function of the monomer
        composition and the reactivity ratios.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 181.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.

        Returns
        -------
        dict[str, float | FloatArray]
            Triad fractions,
            {'111': $F_{111}$, '112': $F_{112}$, '212': $F_{212}$, ... }.
        """

        P111, P211, _, _, P222, P122, _, _ = self.transitions(f1).values()

        F111 = P211*P111/(1. - P111 + eps)
        F112 = 2*P211
        F212 = 1. - P211
        Fsum = F111 + F112 + F212
        F111 /= Fsum
        F112 /= Fsum
        F212 /= Fsum

        F222 = P122*P222/(1. - P222 + eps)
        F221 = 2*P122
        F121 = 1. - P122
        Fsum = F222 + F221 + F121
        F222 /= Fsum
        F221 /= Fsum
        F121 /= Fsum

        result = {'111': F111,
                  '112': F112,
                  '212': F212,
                  '222': F222,
                  '221': F221,
                  '121': F121}

        return result

    def sequence(self,
                 f1: Union[float, FloatArrayLike],
                 k: Optional[Union[int, IntArrayLike]] = None,
                 ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous sequence length probability or the
        number-average sequence length.

        For a binary system, the probability of finding $k$ consecutive units
        of monomer $i$ in a chain is:

        $$ S_{i,k} = \begin{cases}
            1 - P_{jii} & \text{if } k = 1 \\
            P_{jii}(1 - P_{iii})P_{iii}^{k-2}& \text{if } k \ge 2
        \end{cases} $$

        and the corresponding number-average sequence length is:

        $$ \bar{S}_i = \sum_k k S_{i,k} = 1 + \frac{P_{jii}}{1 - P_{iii}} $$

        where $P_{ijk}$ is the transition probability
        $i \rightarrow j \rightarrow k$, which is a function of the monomer
        composition and the reactivity ratios.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 180.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.
        k : int | IntArrayLike | None
            Sequence length, i.e., number of consecutive units in a chain.
            If `None`, the number-average sequence length will be computed.

        Returns
        -------
        dict[str, float | FloatArray]
            If `k is None`, the number-average sequence lengths,
            {'1': $\bar{S}_1$, '2': $\bar{S}_2$}. Otherwise, the
            sequence probabilities, {'1': $S_{1,k}$, '2': $S_{2,k}$}.
        """

        P111, P211, _, _, P222, P122, _, _ = self.transitions(f1).values()

        if k is None:
            S1 = 1. + P211/(1. - P111 + eps)
            S2 = 1. + P122/(1. - P222 + eps)
        else:
            if isinstance(k, (list, tuple)):
                k = np.array(k, dtype=np.int32)
            S1 = np.where(k == 1, 1. - P211, P211*(1. - P111)*P111**(k - 2))
            S2 = np.where(k == 1, 1. - P122, P122*(1. - P222)*P222**(k - 2))

        return {'1': S1, '2': S2}

F1

F1(
    f1: Union[float, FloatArrayLike]
) -> Union[float, FloatArray]

Calculate the instantaneous copolymer composition, \(F_1\).

The calculation is handled by inst_copolymer_binary.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
float | FloatArray	Instantaneous copolymer composition, \(F_1\).

Source code in src/polykin/copolymerization/terminal.py

def F1(self,
       f1: Union[float, FloatArrayLike]
       ) -> Union[float, FloatArray]:
    r"""Calculate the instantaneous copolymer composition, $F_1$.

    The calculation is handled by
    [`inst_copolymer_binary`](inst_copolymer_binary.md).

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    float | FloatArray
        Instantaneous copolymer composition, $F_1$.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return inst_copolymer_binary(f1, *self.ri(f1))

add_data

add_data(
    data: Union[CopoDataset, list[CopoDataset]]
) -> None

Add a copolymerization dataset for subsequent analysis.

PARAMETER	DESCRIPTION
data	Experimental dataset(s). TYPE: Union[CopoDataset, list[CopoDataset]]

Source code in src/polykin/copolymerization/terminal.py

def add_data(self,
             data: Union[CopoDataset, list[CopoDataset]]
             ) -> None:
    r"""Add a copolymerization dataset for subsequent analysis.

    Parameters
    ----------
    data : Union[CopoDataset, list[CopoDataset]]
        Experimental dataset(s).
    """
    if not isinstance(data, list):
        data = [data]

    valid_monomers = {self.M1, self.M2}
    for ds in data:
        ds_monomers = {ds.M1, ds.M2}
        if ds_monomers != valid_monomers:
            raise ValueError(
                f"Monomers defined in dataset `{ds.name}` are invalid: "
                f"{ds_monomers}!={valid_monomers}")
        if ds not in set(self.data):
            self.data.append(ds)
        else:
            print(f"Warning: duplicate dataset '{ds.name}' was skipped.")

azeotrope property

azeotrope: Optional[float]

Calculate the azeotrope composition.

RETURNS	DESCRIPTION
float | None	If an azeotrope exists, it returns its composition in terms of \(f_1\).

drift

drift(
    f10: Union[float, FloatVectorLike],
    x: Union[float, FloatVectorLike],
) -> FloatArray

Calculate drift of comonomer composition in a closed system for a given total monomer conversion.

In a closed binary system, the drift in monomer composition is given by the solution of the following differential equation:

\[ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} \]

with initial condition \(f_1(0)=f_{1,0}\), where \(f_1\) and \(F_1\) are, respectively, the instantaneous comonomer and copolymer composition of M1, and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f10	Initial molar fraction of M1, \(f_{1,0}=f_1(0)\). TYPE: float | FloatVectorLike
x	Value(s) of total monomer conversion values where the drift is to be evaluated. TYPE: float | FloatVectorLike

RETURNS	DESCRIPTION
FloatArray	Monomer fraction of M1 at a given conversion, \(f_1(x)\).

Source code in src/polykin/copolymerization/terminal.py

def drift(self,
          f10: Union[float, FloatVectorLike],
          x: Union[float, FloatVectorLike]
          ) -> FloatArray:
    r"""Calculate drift of comonomer composition in a closed system for a
    given total monomer conversion.

    In a closed binary system, the drift in monomer composition is given by
    the solution of the following differential equation:

    $$ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} $$

    with initial condition $f_1(0)=f_{1,0}$, where $f_1$ and $F_1$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    M1, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f10 : float | FloatVectorLike
        Initial molar fraction of M1, $f_{1,0}=f_1(0)$.
    x : float | FloatVectorLike
        Value(s) of total monomer conversion values where the drift is to
        be evaluated.

    Returns
    -------
    FloatArray
        Monomer fraction of M1 at a given conversion, $f_1(x)$.
    """
    f10, x = convert_FloatOrVectorLike_to_FloatVector([f10, x], False)

    def df1dx(x_, f1_):
        F1_ = inst_copolymer_binary(f1_, *self.ri(f1_))
        return (f1_ - F1_)/(1. - x_ + eps)

    sol = solve_ivp(df1dx,
                    (0., max(x)),
                    f10,
                    t_eval=x,
                    method='LSODA',
                    vectorized=True,
                    atol=1e-4,
                    rtol=1e-4)
    if sol.success:
        result = sol.y
        result = np.maximum(0., result)
        if result.shape[0] == 1:
            result = result[0]
    else:
        result = np.empty_like(x)
        result[:] = np.nan
        print(sol.message)

    return result

kii

kii(
    f1: Union[float, FloatArray],
    T: float,
    Tunit: Literal["C", "K"] = "K",
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Pseudohomopropagation rate coefficients.

In the penultimate model, the pseudohomopropagation rate coefficients depend on the instantaneous comonomer composition according to:

\[\begin{aligned} \bar{k}_{11} &= k_{111}\frac{f_1 r_{11} + f_2}{f_1 r_{11} + f_2/s_1} \\ \bar{k}_{22} &= k_{222}\frac{f_2 r_{22} + f_1}{f_2 r_{22} + f_1/s_2} \end{aligned}\]

where \(r_{ij}\) are the monomer reactivity ratios, \(s_i\) are the radical reactivity ratios, and \(k_{iii}\) are the homopropagation rate coefficients.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArray
T	Temperature. Unit = Tunit. TYPE: float
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'

RETURNS	DESCRIPTION
tuple[float | FloatArray, float | FloatArray]	Pseudohomopropagation rate coefficients, (\(\bar{k}_{11}\), \(\bar{k}_{22}\)).

Source code in src/polykin/copolymerization/penultimate.py

def kii(self,
        f1: Union[float, FloatArray],
        T: float,
        Tunit: Literal['C', 'K'] = 'K',
        ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    r"""Pseudohomopropagation rate coefficients.

    In the penultimate model, the pseudohomopropagation rate coefficients
    depend on the instantaneous comonomer composition according to:

    \begin{aligned}
    \bar{k}_{11} &= k_{111}\frac{f_1 r_{11} + f_2}{f_1 r_{11} + f_2/s_1} \\
    \bar{k}_{22} &= k_{222}\frac{f_2 r_{22} + f_1}{f_2 r_{22} + f_1/s_2}
    \end{aligned}

    where $r_{ij}$ are the monomer reactivity ratios, $s_i$ are the radical
    reactivity ratios, and $k_{iii}$ are the homopropagation rate
    coefficients.

    Parameters
    ----------
    f1 : float | FloatArray
        Molar fraction of M1.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.

    Returns
    -------
    tuple[float | FloatArray, float | FloatArray]
        Pseudohomopropagation rate coefficients,
        ($\bar{k}_{11}$, $\bar{k}_{22}$).
    """
    if self.k1 is None or self.k2 is None:
        raise ValueError(
            "To use this feature, `k1` and `k2` cannot be `None`.")
    f2 = 1. - f1
    r11 = self.r11
    r22 = self.r22
    s1 = self.s1
    s2 = self.s2
    k1 = self.k1(T, Tunit)
    k2 = self.k2(T, Tunit)
    k11 = k1*(f1*r11 + f2)/(f1*r11 + f2/s1)
    k22 = k2*(f2*r22 + f1)/(f2*r22 + f1/s2)
    return (k11, k22)

kp

kp(
    f1: Union[float, FloatArrayLike],
    T: float,
    Tunit: Literal["C", "K"] = "K",
) -> Union[float, FloatArray]

Calculate the average propagation rate coefficient, \(\bar{k}_p\).

The calculation is handled by kp_average_binary.

Note

This feature requires the attributes k11 and k22 to be defined.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
T	Temperature. Unit = Tunit. TYPE: float
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'

RETURNS	DESCRIPTION
float | FloatArray	Average propagation rate coefficient. Unit = L/(mol·s)

Source code in src/polykin/copolymerization/terminal.py

def kp(self,
       f1: Union[float, FloatArrayLike],
       T: float,
       Tunit: Literal['C', 'K'] = 'K'
       ) -> Union[float, FloatArray]:
    r"""Calculate the average propagation rate coefficient, $\bar{k}_p$.

    The calculation is handled by
    [`kp_average_binary`](kp_average_binary.md).

    Note
    ----
    This feature requires the attributes `k11` and `k22` to be defined.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.

    Returns
    -------
    float | FloatArray
        Average propagation rate coefficient. Unit = L/(mol·s)
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return kp_average_binary(f1, *self.ri(f1), *self.kii(f1, T, Tunit))

plot

plot(
    kind: Literal["drift", "kp", "Mayo", "triads"],
    show: Literal["auto", "all", "data", "model"] = "auto",
    M: Literal[1, 2] = 1,
    f0: Optional[Union[float, FloatVectorLike]] = None,
    T: Optional[float] = None,
    Tunit: Literal["C", "K"] = "K",
    title: Optional[str] = None,
    axes: Optional[Axes] = None,
    return_objects: bool = False,
) -> Optional[tuple[Optional[Figure], Axes]]

Generate a plot of instantaneous copolymer composition, monomer composition drift, or average propagation rate coefficient.

PARAMETER	DESCRIPTION
kind	Kind of plot to be generated. TYPE: Literal['drift', 'kp', 'Mayo', 'triads']
show	What informatation is to be plotted. TYPE: Literal['auto', 'all', 'data', 'model'] DEFAULT: 'auto'
M	Index of the monomer to be used in input argument f0 and in output results. Specifically, if M=i, then f0 stands for \(f_i(0)\) and plots will be generated in terms of \(f_i\) and \(F_i\). TYPE: Literal[1, 2] DEFAULT: 1
f0	Initial monomer composition, \(f_i(0)\), as required for a monomer composition drift plot. TYPE: float | FloatVectorLike | None DEFAULT: None
T	Temperature. Unit = Tunit. TYPE: float | None DEFAULT: None
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
title	Title of plot. If None, a default title with the monomer names will be used. TYPE: str | None DEFAULT: None
axes	Matplotlib Axes object. TYPE: Axes | None DEFAULT: None
return_objects	If True, the Figure and Axes objects are returned (for saving or further manipulations). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
tuple[Figure | None, Axes] | None	Figure and Axes objects if return_objects is True.

Source code in src/polykin/copolymerization/terminal.py

def plot(self,
         kind: Literal['drift', 'kp', 'Mayo', 'triads'],
         show: Literal['auto', 'all', 'data', 'model'] = 'auto',
         M: Literal[1, 2] = 1,
         f0: Optional[Union[float, FloatVectorLike]] = None,
         T: Optional[float] = None,
         Tunit: Literal['C', 'K'] = 'K',
         title: Optional[str] = None,
         axes: Optional[Axes] = None,
         return_objects: bool = False
         ) -> Optional[tuple[Optional[Figure], Axes]]:
    r"""Generate a plot of instantaneous copolymer composition, monomer
    composition drift, or average propagation rate coefficient.

    Parameters
    ----------
    kind : Literal['drift', 'kp', 'Mayo', 'triads']
        Kind of plot to be generated.
    show : Literal['auto', 'all', 'data', 'model']
        What informatation is to be plotted.
    M : Literal[1, 2]
        Index of the monomer to be used in input argument `f0` and in
        output results. Specifically, if `M=i`, then `f0` stands for
        $f_i(0)$ and plots will be generated in terms of $f_i$ and $F_i$.
    f0 : float | FloatVectorLike | None
        Initial monomer composition, $f_i(0)$, as required for a monomer
        composition drift plot.
    T : float | None
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    title : str | None
        Title of plot. If `None`, a default title with the monomer names
        will be used.
    axes : Axes | None
        Matplotlib Axes object.
    return_objects : bool
        If `True`, the Figure and Axes objects are returned (for saving or
        further manipulations).

    Returns
    -------
    tuple[Figure | None, Axes] | None
        Figure and Axes objects if return_objects is `True`.
    """
    check_in_set(M, {1, 2}, 'M')
    check_in_set(kind, {'Mayo', 'kp', 'drift', 'triads'}, 'kind')
    check_in_set(show, {'auto', 'all', 'data', 'model'}, 'show')

    label_model = None
    if axes is None:
        fig, ax = plt.subplots()
        if title is None:
            titles = {'Mayo': "Mayo-Lewis diagram",
                      'drift': "Monomer composition drift",
                      'kp': "Average propagation coefficient",
                      'triads': "Triad fractions"}
            title = titles[kind] + f" {self.M1}(1)-{self.M2}(2)"
        if title:
            fig.suptitle(title)
    else:
        ax = axes
        fig = None
        if self.name:
            label_model = self.name

    unit_range = (0., 1.)
    npoints = 1000
    Mname = self.M1 if M == 1 else self.M2
    ndataseries = 0

    if show == 'auto':
        if not self.data:
            show = 'model'
        else:
            show = 'all'

    if kind == 'Mayo':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(fr"$F_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        ax.plot(unit_range, unit_range, color='black', linewidth=0.5)

        if show in {'model', 'all'}:
            x = np.linspace(*unit_range, npoints, dtype=np.float64)
            y = self.F1(x)
            if M == 2:
                x[:] = 1. - x
                y[:] = 1. - y  # type: ignore
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:
            for ds in self.data:
                if not isinstance(ds, MayoDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('F', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

    elif kind == 'drift':

        ax.set_xlabel(r"Total molar monomer conversion, $x$")
        ax.set_ylabel(fr"$f_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        if show == 'model':

            if f0 is None:
                raise ValueError("`f0` is required for a `drift` plot.")
            else:
                if isinstance(f0, (int, float)):
                    f0 = [f0]
                f0 = np.array(f0, dtype=np.float64)
                check_bounds(f0, *unit_range, 'f0')

            x = np.linspace(*unit_range, 1000)
            if M == 2:
                f0[:] = 1. - f0
            y = self.drift(f0, x)
            if M == 2:
                y[:] = 1. - y
            if y.ndim == 1:
                y = y[np.newaxis, :]
            for i in range(y.shape[0]):
                ax.plot(x, y[i, :], label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, DriftDataset):
                    continue
                x = ds.getvar('x')
                y = ds.getvar('f', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.drift(ds.getvar('f', self.M1)[0], x)
                    if M == 2:
                        y[:] = 1. - y
                    ax.plot(x, y, label=label_model)

    elif kind == 'kp':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(r"$\bar{k}_p$")
        ax.set_xlim(*unit_range)

        if show == 'model':

            if T is None:
                raise ValueError("`T` is required for a `kp` plot.")

            x = np.linspace(*unit_range, npoints)
            y = self.kp(x, T, Tunit)
            if M == 2:
                x[:] = 1. - x
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, kpDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('kp')
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)
                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.kp(x, ds.T, ds.Tunit)
                    if M == 2:
                        x[:] = 1. - x
                    ax.plot(x, y, label=label_model)

    ax.grid(True)

    if axes is not None or ndataseries:
        ax.legend(bbox_to_anchor=(1.05, 1.), loc="upper left")

    if return_objects:
        return (fig, ax)

ri

ri(
    f1: Union[float, FloatArray]
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Pseudoreactivity ratios.

In the penultimate model, the pseudoreactivity ratios depend on the instantaneous comonomer composition according to:

\[\begin{aligned} \bar{r}_1 &= r_{21}\frac{f_1 r_{11} + f_2}{f_1 r_{21} + f_2} \\ \bar{r}_2 &= r_{12}\frac{f_2 r_{22} + f_1}{f_2 r_{12} + f_1} \end{aligned}\]

where \(r_{ij}\) are the monomer reactivity ratios.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArray

RETURNS	DESCRIPTION
tuple[FloatOrArray, FloatOrArray]	Pseudoreactivity ratios, (\(\bar{r}_1\), \(\bar{r}_2\)).

Source code in src/polykin/copolymerization/penultimate.py

def ri(self,
        f1: Union[float, FloatArray]
       ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    r"""Pseudoreactivity ratios.

    In the penultimate model, the pseudoreactivity ratios depend on the
    instantaneous comonomer composition according to:

    \begin{aligned}
       \bar{r}_1 &= r_{21}\frac{f_1 r_{11} + f_2}{f_1 r_{21} + f_2} \\
       \bar{r}_2 &= r_{12}\frac{f_2 r_{22} + f_1}{f_2 r_{12} + f_1}
    \end{aligned}

    where $r_{ij}$ are the monomer reactivity ratios.

    Parameters
    ----------
    f1 : float | FloatArray
        Molar fraction of M1.

    Returns
    -------
    tuple[FloatOrArray, FloatOrArray]
        Pseudoreactivity ratios, ($\bar{r}_1$, $\bar{r}_2$).
    """
    f2 = 1. - f1
    r11 = self.r11
    r12 = self.r12
    r21 = self.r21
    r22 = self.r22
    r1 = r21*(f1*r11 + f2)/(f1*r21 + f2)
    r2 = r12*(f2*r22 + f1)/(f2*r12 + f1)
    return (r1, r2)

sequence

sequence(
    f1: Union[float, FloatArrayLike],
    k: Optional[Union[int, IntArrayLike]] = None,
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous sequence length probability or the number-average sequence length.

For a binary system, the probability of finding \(k\) consecutive units of monomer \(i\) in a chain is:

\[ S_{i,k} = \begin{cases} 1 - P_{jii} & \text{if } k = 1 \\ P_{jii}(1 - P_{iii})P_{iii}^{k-2}& \text{if } k \ge 2 \end{cases} \]

and the corresponding number-average sequence length is:

\[ \bar{S}_i = \sum_k k S_{i,k} = 1 + \frac{P_{jii}}{1 - P_{iii}} \]

where \(P_{ijk}\) is the transition probability \(i \rightarrow j \rightarrow k\), which is a function of the monomer composition and the reactivity ratios.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 180.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
k	Sequence length, i.e., number of consecutive units in a chain. If None, the number-average sequence length will be computed. TYPE: int | IntArrayLike | None DEFAULT: None

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	If k is None, the number-average sequence lengths, {'1': \(\bar{S}_1\), '2': \(\bar{S}_2\)}. Otherwise, the sequence probabilities, {'1': \(S_{1,k}\), '2': \(S_{2,k}\)}.

Source code in src/polykin/copolymerization/penultimate.py

def sequence(self,
             f1: Union[float, FloatArrayLike],
             k: Optional[Union[int, IntArrayLike]] = None,
             ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous sequence length probability or the
    number-average sequence length.

    For a binary system, the probability of finding $k$ consecutive units
    of monomer $i$ in a chain is:

    $$ S_{i,k} = \begin{cases}
        1 - P_{jii} & \text{if } k = 1 \\
        P_{jii}(1 - P_{iii})P_{iii}^{k-2}& \text{if } k \ge 2
    \end{cases} $$

    and the corresponding number-average sequence length is:

    $$ \bar{S}_i = \sum_k k S_{i,k} = 1 + \frac{P_{jii}}{1 - P_{iii}} $$

    where $P_{ijk}$ is the transition probability
    $i \rightarrow j \rightarrow k$, which is a function of the monomer
    composition and the reactivity ratios.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 180.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    k : int | IntArrayLike | None
        Sequence length, i.e., number of consecutive units in a chain.
        If `None`, the number-average sequence length will be computed.

    Returns
    -------
    dict[str, float | FloatArray]
        If `k is None`, the number-average sequence lengths,
        {'1': $\bar{S}_1$, '2': $\bar{S}_2$}. Otherwise, the
        sequence probabilities, {'1': $S_{1,k}$, '2': $S_{2,k}$}.
    """

    P111, P211, _, _, P222, P122, _, _ = self.transitions(f1).values()

    if k is None:
        S1 = 1. + P211/(1. - P111 + eps)
        S2 = 1. + P122/(1. - P222 + eps)
    else:
        if isinstance(k, (list, tuple)):
            k = np.array(k, dtype=np.int32)
        S1 = np.where(k == 1, 1. - P211, P211*(1. - P111)*P111**(k - 2))
        S2 = np.where(k == 1, 1. - P122, P122*(1. - P222)*P222**(k - 2))

    return {'1': S1, '2': S2}

transitions

transitions(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous transition probabilities.

For a binary system, the transition probabilities are given by:

\[\begin{aligned} P_{iii} &= \frac{r_{ii} f_i}{r_{ii} f_i + (1 - f_i)} \\ P_{jii} &= \frac{r_{ji} f_i}{r_{ji} f_i + (1 - f_i)} \\ P_{iij} &= 1 - P_{iii} \\ P_{jij} &= 1 - P_{jii} \end{aligned}\]

where \(i,j=1,2, i \neq j\).

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 181.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Transition probabilities, {'111': \(P_{111}\), '211': \(P_{211}\), '121': \(P_{121}\), ... }.

Source code in src/polykin/copolymerization/penultimate.py

def transitions(self,
                f1: Union[float, FloatArrayLike]
                ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous transition probabilities.

    For a binary system, the transition probabilities are given by:

    \begin{aligned}
        P_{iii} &= \frac{r_{ii} f_i}{r_{ii} f_i + (1 - f_i)} \\
        P_{jii} &= \frac{r_{ji} f_i}{r_{ji} f_i + (1 - f_i)} \\
        P_{iij} &= 1 -  P_{iii} \\
        P_{jij} &= 1 -  P_{jii}
    \end{aligned}

    where $i,j=1,2, i \neq j$.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 181.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Transition probabilities,
        {'111': $P_{111}$, '211': $P_{211}$, '121': $P_{121}$, ... }.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')

    f2 = 1. - f1
    r11 = self.r11
    r12 = self.r12
    r21 = self.r21
    r22 = self.r22

    P111 = r11*f1/(r11*f1 + f2)
    P211 = r21*f1/(r21*f1 + f2)
    P112 = 1. - P111
    P212 = 1. - P211

    P222 = r22*f2/(r22*f2 + f1)
    P122 = r12*f2/(r12*f2 + f1)
    P221 = 1. - P222
    P121 = 1. - P122

    result = {'111': P111,
              '211': P211,
              '112': P112,
              '212': P212,
              '222': P222,
              '122': P122,
              '221': P221,
              '121': P121}

    return result

triads

triads(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous triad fractions.

For a binary system, the triad fractions are given by:

\[\begin{aligned} F_{iii} &\propto P_{jii} \frac{P_{iii}}{1 - P_{iii}} \\ F_{iij} &\propto 2 P_{jii} \\ F_{jij} &\propto 1 - P_{jii} \end{aligned}\]

where \(P_{ijk}\) is the transition probability \(i \rightarrow j \rightarrow k\), which is a function of the monomer composition and the reactivity ratios.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 181.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Triad fractions, {'111': \(F_{111}\), '112': \(F_{112}\), '212': \(F_{212}\), ... }.

Source code in src/polykin/copolymerization/penultimate.py

def triads(self,
           f1: Union[float, FloatArrayLike],
           ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous triad fractions.

    For a binary system, the triad fractions are given by:

    \begin{aligned}
        F_{iii} &\propto  P_{jii} \frac{P_{iii}}{1 - P_{iii}} \\
        F_{iij} &\propto 2 P_{jii} \\
        F_{jij} &\propto 1 - P_{jii}
    \end{aligned}

    where $P_{ijk}$ is the transition probability
    $i \rightarrow j \rightarrow k$, which is a function of the monomer
    composition and the reactivity ratios.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 181.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Triad fractions,
        {'111': $F_{111}$, '112': $F_{112}$, '212': $F_{212}$, ... }.
    """

    P111, P211, _, _, P222, P122, _, _ = self.transitions(f1).values()

    F111 = P211*P111/(1. - P111 + eps)
    F112 = 2*P211
    F212 = 1. - P211
    Fsum = F111 + F112 + F212
    F111 /= Fsum
    F112 /= Fsum
    F212 /= Fsum

    F222 = P122*P222/(1. - P222 + eps)
    F221 = 2*P122
    F121 = 1. - P122
    Fsum = F222 + F221 + F121
    F222 /= Fsum
    F221 /= Fsum
    F121 /= Fsum

    result = {'111': F111,
              '112': F112,
              '212': F212,
              '222': F222,
              '221': F221,
              '121': F121}

    return result

TerminalModel

Terminal binary copolymerization model.

According to this model, the reactivity of a macroradical depends uniquely on the nature of the terminal repeating unit. A binary system is, thus, described by four propagation reactions:

\[\begin{matrix} P^{\bullet}_1 + M_1 \overset{k_{11}}{\rightarrow} P^{\bullet}_1 \\ P^{\bullet}_1 + M_2 \overset{k_{12}}{\rightarrow} P^{\bullet}_2 \\ P^{\bullet}_2 + M_1 \overset{k_{21}}{\rightarrow} P^{\bullet}_1 \\ P^{\bullet}_2 + M_2 \overset{k_{22}}{\rightarrow} P^{\bullet}_2 \end{matrix}\]

where \(k_{ii}\) are the homopropagation rate coefficients and \(k_{ij}\) are the cross-propagation coefficients. The two cross-propagation coefficients are specified via an equal number of reactivity ratios, defined as \(r_1=k_{11}/k_{12}\) and \(r_2=k_{22}/k_{21}\).

PARAMETER	DESCRIPTION
r1	Reactivity ratio of M1, \(r_1=k_{11}/k_{12}\). TYPE: float
r2	Reactivity ratio of M2, \(r_2=k_{22}/k_{21}\). TYPE: float
k1	Homopropagation rate coefficient of M1, \(k_1 \equiv k_{11}\). TYPE: Arrhenius | None DEFAULT: None
k2	Homopropagation rate coefficient of M2, \(k_2 \equiv k_{22}\). TYPE: Arrhenius | None DEFAULT: None
M1	Name of M1. TYPE: str DEFAULT: 'M1'
M2	Name of M2. TYPE: str DEFAULT: 'M2'
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/terminal.py

class TerminalModel(CopoModel):
    r"""Terminal binary copolymerization model.

    According to this model, the reactivity of a macroradical depends
    uniquely on the nature of the _terminal_ repeating unit. A binary system
    is, thus, described by four propagation reactions:

    \begin{matrix}
    P^{\bullet}_1 + M_1 \overset{k_{11}}{\rightarrow} P^{\bullet}_1 \\
    P^{\bullet}_1 + M_2 \overset{k_{12}}{\rightarrow} P^{\bullet}_2 \\
    P^{\bullet}_2 + M_1 \overset{k_{21}}{\rightarrow} P^{\bullet}_1 \\
    P^{\bullet}_2 + M_2 \overset{k_{22}}{\rightarrow} P^{\bullet}_2
    \end{matrix}

    where $k_{ii}$ are the homopropagation rate coefficients and $k_{ij}$ are
    the cross-propagation coefficients. The two cross-propagation coefficients
    are specified via an equal number of reactivity ratios, defined as
    $r_1=k_{11}/k_{12}$ and $r_2=k_{22}/k_{21}$.

    Parameters
    ----------
    r1 : float
        Reactivity ratio of M1, $r_1=k_{11}/k_{12}$.
    r2 : float
        Reactivity ratio of M2, $r_2=k_{22}/k_{21}$.
    k1 : Arrhenius | None
        Homopropagation rate coefficient of M1, $k_1 \equiv k_{11}$.
    k2 : Arrhenius | None
        Homopropagation rate coefficient of M2, $k_2 \equiv k_{22}$.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    name : str
        Name.
    """

    r1: float
    r2: float

    _pnames = ('r1', 'r2')

    def __init__(self,
                 r1: float,
                 r2: float,
                 k1: Optional[Arrhenius] = None,
                 k2: Optional[Arrhenius] = None,
                 M1: str = 'M1',
                 M2: str = 'M2',
                 name: str = ''
                 ) -> None:

        check_bounds(r1, 0., np.inf, 'r1')
        check_bounds(r2, 0., np.inf, 'r2')

        # Perhaps this could be upgraded to exception, but I don't want to be
        # too restrictive (one does find literature data with (r1,r2)>1)
        if r1 > 1. and r2 > 1.:
            print(
                f"Warning: `r1`={r1} and `r2`={r2} are both greater than 1, "
                "which is deemed physically impossible.")

        self.r1 = r1
        self.r2 = r2
        super().__init__(k1, k2, M1, M2, name)

    @classmethod
    def from_Qe(cls,
                Qe1: tuple[float, float],
                Qe2: tuple[float, float],
                k1: Optional[Arrhenius] = None,
                k2: Optional[Arrhenius] = None,
                M1: str = 'M1',
                M2: str = 'M2',
                name: str = ''
                ) -> TerminalModel:
        r"""Construct `TerminalModel` from Q-e values.

        Alternative constructor that takes the $Q$-$e$ values of the monomers
        as primary input instead of the reactivity ratios.

        The conversion from Q-e to r is handled by
        [`convert_Qe_to_r`](convert_Qe_to_r.md).

        Parameters
        ----------
        Qe1 : tuple[float, float]
            Q-e values of M1.
        Qe2 : tuple[float, float]
            Q-e values of M2.
        k1 : Arrhenius | None
            Homopropagation rate coefficient of M1, $k_1 \equiv k_{11}$.
        k2 : Arrhenius | None
            Homopropagation rate coefficient of M2, $k_2 \equiv k_{22}$.
        M1 : str
            Name of M1.
        M2 : str
            Name of M2.
        name : str
            Name.
        """
        r = convert_Qe_to_r([Qe1, Qe2])
        r1 = r[0, 1]
        r2 = r[1, 0]
        return cls(r1, r2, k1, k2, M1, M2, name)

    @property
    def azeotrope(self) -> Optional[float]:
        r"""Calculate the azeotrope composition.

        An azeotrope (i.e., a point where $F_1=f_1$) only exists if both
        reactivity ratios are smaller than unity. In that case, the azeotrope
        composition is given by:

        $$ f_{1,azeo} = \frac{1 - r_2}{(2 - r_1 - r_2)} $$

        where $r_1$ and $r_2$ are the reactivity ratios.

        Returns
        -------
        float | None
            If an azeotrope exists, it returns its composition in terms of
            $f_1$.
        """
        r1 = self.r1
        r2 = self.r2
        if r1 < 1. and r2 < 1.:
            result = (1. - r2)/(2. - r1 - r2)
        else:
            result = None
        return result

    def ri(self,
           _
           ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
        return (self.r1, self.r2)

    def kii(self,
            _,
            T: float,
            Tunit: Literal['C', 'K'] = 'K'
            ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
        if self.k1 is None or self.k2 is None:
            raise ValueError(
                "To use this feature, `k1` and `k2` cannot be `None`.")
        return (self.k1(T, Tunit), self.k2(T, Tunit))

    def transitions(self,
                    f1: Union[float, FloatArrayLike]
                    ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous transition probabilities.

        For a binary system, the transition probabilities are given by:

        \begin{aligned}
            P_{ii} &= \frac{r_i f_i}{r_i f_i + (1 - f_i)} \\
            P_{ij} &= 1 - P_{ii}
        \end{aligned}

        where $i,j=1,2, i \neq j$.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 178.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.

        Returns
        -------
        dict[str, float | FloatArray]
            Transition probabilities, {'11': $P_{11}$, '12': $P_{12}$, ... }.
        """

        if isinstance(f1, (list, tuple)):
            f1 = np.array(f1, dtype=np.float64)
        check_bounds(f1, 0., 1., 'f1')

        f2 = 1. - f1
        r1 = self.r1
        r2 = self.r2
        P11 = r1*f1/(r1*f1 + f2)
        P22 = r2*f2/(r2*f2 + f1)
        P12 = 1. - P11
        P21 = 1. - P22

        result = {'11': P11,
                  '12': P12,
                  '21': P21,
                  '22': P22}

        return result

    def triads(self,
               f1: Union[float, FloatArrayLike]
               ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous triad fractions.

        For a binary system, the triad fractions are given by:

        \begin{aligned}
            A_{iii} &= F_i P_{ii}^2 \\
            A_{iij} &= 2 F_i P_{ii} P_{ij} \\
            A_{jij} &= F_i P_{ij}^2
        \end{aligned}

        where $P_{ij}$ is the transition probability $i \rightarrow j$ and
        $F_i$ is the instantaneous copolymer composition.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 179.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.

        Returns
        -------
        dict[str, float | FloatArray]
            Triad fractions,
            {'111': $A_{111}$, '112': $A_{112}$, '212': $A_{212}$, ... }.
        """

        P11, P12, P21, P22 = self.transitions(f1).values()

        F1 = P21/(P12 + P21)
        F2 = 1. - F1

        A111 = F1*P11**2
        A112 = F1*2*P11*P12
        A212 = F1*P12**2

        A222 = F2*P22**2
        A221 = F2*2*P22*P21
        A121 = F2*P21**2

        result = {'111': A111,
                  '112': A112,
                  '212': A212,
                  '222': A222,
                  '221': A221,
                  '121': A121}

        return result

    def sequence(self,
                 f1: Union[float, FloatArrayLike],
                 k: Optional[Union[int, IntArrayLike]] = None,
                 ) -> dict[str, Union[float, FloatArray]]:
        r"""Calculate the instantaneous sequence length probability or the
        number-average sequence length.

        For a binary system, the probability of finding $k$ consecutive units
        of monomer $i$ in a chain is:

        $$ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} $$

        and the corresponding number-average sequence length is:

        $$ \bar{S}_i = \sum_k k S_{i,k} = \frac{1}{1 - P_{ii}} $$

        where $P_{ii}$ is the transition probability $i \rightarrow i$, which
        is a function of the monomer composition and the reactivity ratios.

        **References**

        * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 177.

        Parameters
        ----------
        f1 : float | FloatArrayLike
            Molar fraction of M1.
        k : int | IntArrayLike | None
            Sequence length, i.e., number of consecutive units in a chain.
            If `None`, the number-average sequence length will be computed.

        Returns
        -------
        dict[str, float | FloatArray]
            If `k is None`, the number-average sequence lengths,
            {'1': $\bar{S}_1$, '2': $\bar{S}_2$}. Otherwise, the
            sequence probabilities, {'1': $S_{1,k}$, '2': $S_{2,k}$}.
        """

        P11, _, _, P22 = self.transitions(f1).values()

        if k is None:
            result = {str(i + 1): 1/(1 - P + eps)
                      for i, P in enumerate([P11, P22])}
        else:
            if isinstance(k, (list, tuple)):
                k = np.array(k, dtype=np.int32)
            result = {str(i + 1): (1. - P)*P**(k - 1)
                      for i, P in enumerate([P11, P22])}

        return result

F1

F1(
    f1: Union[float, FloatArrayLike]
) -> Union[float, FloatArray]

Calculate the instantaneous copolymer composition, \(F_1\).

The calculation is handled by inst_copolymer_binary.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
float | FloatArray	Instantaneous copolymer composition, \(F_1\).

Source code in src/polykin/copolymerization/terminal.py

def F1(self,
       f1: Union[float, FloatArrayLike]
       ) -> Union[float, FloatArray]:
    r"""Calculate the instantaneous copolymer composition, $F_1$.

    The calculation is handled by
    [`inst_copolymer_binary`](inst_copolymer_binary.md).

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    float | FloatArray
        Instantaneous copolymer composition, $F_1$.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return inst_copolymer_binary(f1, *self.ri(f1))

add_data

add_data(
    data: Union[CopoDataset, list[CopoDataset]]
) -> None

Add a copolymerization dataset for subsequent analysis.

PARAMETER	DESCRIPTION
data	Experimental dataset(s). TYPE: Union[CopoDataset, list[CopoDataset]]

Source code in src/polykin/copolymerization/terminal.py

def add_data(self,
             data: Union[CopoDataset, list[CopoDataset]]
             ) -> None:
    r"""Add a copolymerization dataset for subsequent analysis.

    Parameters
    ----------
    data : Union[CopoDataset, list[CopoDataset]]
        Experimental dataset(s).
    """
    if not isinstance(data, list):
        data = [data]

    valid_monomers = {self.M1, self.M2}
    for ds in data:
        ds_monomers = {ds.M1, ds.M2}
        if ds_monomers != valid_monomers:
            raise ValueError(
                f"Monomers defined in dataset `{ds.name}` are invalid: "
                f"{ds_monomers}!={valid_monomers}")
        if ds not in set(self.data):
            self.data.append(ds)
        else:
            print(f"Warning: duplicate dataset '{ds.name}' was skipped.")

azeotrope property

azeotrope: Optional[float]

Calculate the azeotrope composition.

An azeotrope (i.e., a point where \(F_1=f_1\)) only exists if both reactivity ratios are smaller than unity. In that case, the azeotrope composition is given by:

\[ f_{1,azeo} = \frac{1 - r_2}{(2 - r_1 - r_2)} \]

where \(r_1\) and \(r_2\) are the reactivity ratios.

RETURNS	DESCRIPTION
float | None	If an azeotrope exists, it returns its composition in terms of \(f_1\).

drift

drift(
    f10: Union[float, FloatVectorLike],
    x: Union[float, FloatVectorLike],
) -> FloatArray

Calculate drift of comonomer composition in a closed system for a given total monomer conversion.

In a closed binary system, the drift in monomer composition is given by the solution of the following differential equation:

\[ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} \]

with initial condition \(f_1(0)=f_{1,0}\), where \(f_1\) and \(F_1\) are, respectively, the instantaneous comonomer and copolymer composition of M1, and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f10	Initial molar fraction of M1, \(f_{1,0}=f_1(0)\). TYPE: float | FloatVectorLike
x	Value(s) of total monomer conversion values where the drift is to be evaluated. TYPE: float | FloatVectorLike

RETURNS	DESCRIPTION
FloatArray	Monomer fraction of M1 at a given conversion, \(f_1(x)\).

Source code in src/polykin/copolymerization/terminal.py

def drift(self,
          f10: Union[float, FloatVectorLike],
          x: Union[float, FloatVectorLike]
          ) -> FloatArray:
    r"""Calculate drift of comonomer composition in a closed system for a
    given total monomer conversion.

    In a closed binary system, the drift in monomer composition is given by
    the solution of the following differential equation:

    $$ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} $$

    with initial condition $f_1(0)=f_{1,0}$, where $f_1$ and $F_1$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    M1, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f10 : float | FloatVectorLike
        Initial molar fraction of M1, $f_{1,0}=f_1(0)$.
    x : float | FloatVectorLike
        Value(s) of total monomer conversion values where the drift is to
        be evaluated.

    Returns
    -------
    FloatArray
        Monomer fraction of M1 at a given conversion, $f_1(x)$.
    """
    f10, x = convert_FloatOrVectorLike_to_FloatVector([f10, x], False)

    def df1dx(x_, f1_):
        F1_ = inst_copolymer_binary(f1_, *self.ri(f1_))
        return (f1_ - F1_)/(1. - x_ + eps)

    sol = solve_ivp(df1dx,
                    (0., max(x)),
                    f10,
                    t_eval=x,
                    method='LSODA',
                    vectorized=True,
                    atol=1e-4,
                    rtol=1e-4)
    if sol.success:
        result = sol.y
        result = np.maximum(0., result)
        if result.shape[0] == 1:
            result = result[0]
    else:
        result = np.empty_like(x)
        result[:] = np.nan
        print(sol.message)

    return result

from_Qe classmethod

from_Qe(
    Qe1: tuple[float, float],
    Qe2: tuple[float, float],
    k1: Optional[Arrhenius] = None,
    k2: Optional[Arrhenius] = None,
    M1: str = "M1",
    M2: str = "M2",
    name: str = "",
) -> TerminalModel

Construct TerminalModel from Q-e values.

Alternative constructor that takes the \(Q\)-\(e\) values of the monomers as primary input instead of the reactivity ratios.

The conversion from Q-e to r is handled by convert_Qe_to_r.

PARAMETER	DESCRIPTION
Qe1	Q-e values of M1. TYPE: tuple[float, float]
Qe2	Q-e values of M2. TYPE: tuple[float, float]
k1	Homopropagation rate coefficient of M1, \(k_1 \equiv k_{11}\). TYPE: Arrhenius | None DEFAULT: None
k2	Homopropagation rate coefficient of M2, \(k_2 \equiv k_{22}\). TYPE: Arrhenius | None DEFAULT: None
M1	Name of M1. TYPE: str DEFAULT: 'M1'
M2	Name of M2. TYPE: str DEFAULT: 'M2'
name	Name. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/terminal.py

@classmethod
def from_Qe(cls,
            Qe1: tuple[float, float],
            Qe2: tuple[float, float],
            k1: Optional[Arrhenius] = None,
            k2: Optional[Arrhenius] = None,
            M1: str = 'M1',
            M2: str = 'M2',
            name: str = ''
            ) -> TerminalModel:
    r"""Construct `TerminalModel` from Q-e values.

    Alternative constructor that takes the $Q$-$e$ values of the monomers
    as primary input instead of the reactivity ratios.

    The conversion from Q-e to r is handled by
    [`convert_Qe_to_r`](convert_Qe_to_r.md).

    Parameters
    ----------
    Qe1 : tuple[float, float]
        Q-e values of M1.
    Qe2 : tuple[float, float]
        Q-e values of M2.
    k1 : Arrhenius | None
        Homopropagation rate coefficient of M1, $k_1 \equiv k_{11}$.
    k2 : Arrhenius | None
        Homopropagation rate coefficient of M2, $k_2 \equiv k_{22}$.
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    name : str
        Name.
    """
    r = convert_Qe_to_r([Qe1, Qe2])
    r1 = r[0, 1]
    r2 = r[1, 0]
    return cls(r1, r2, k1, k2, M1, M2, name)

kii

kii(
    _, T: float, Tunit: Literal["C", "K"] = "K"
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Return the evaluated homopropagation rate coefficients at the given conditions.

Source code in src/polykin/copolymerization/terminal.py

def kii(self,
        _,
        T: float,
        Tunit: Literal['C', 'K'] = 'K'
        ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    if self.k1 is None or self.k2 is None:
        raise ValueError(
            "To use this feature, `k1` and `k2` cannot be `None`.")
    return (self.k1(T, Tunit), self.k2(T, Tunit))

kp

kp(
    f1: Union[float, FloatArrayLike],
    T: float,
    Tunit: Literal["C", "K"] = "K",
) -> Union[float, FloatArray]

Calculate the average propagation rate coefficient, \(\bar{k}_p\).

The calculation is handled by kp_average_binary.

Note

This feature requires the attributes k11 and k22 to be defined.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
T	Temperature. Unit = Tunit. TYPE: float
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'

RETURNS	DESCRIPTION
float | FloatArray	Average propagation rate coefficient. Unit = L/(mol·s)

Source code in src/polykin/copolymerization/terminal.py

def kp(self,
       f1: Union[float, FloatArrayLike],
       T: float,
       Tunit: Literal['C', 'K'] = 'K'
       ) -> Union[float, FloatArray]:
    r"""Calculate the average propagation rate coefficient, $\bar{k}_p$.

    The calculation is handled by
    [`kp_average_binary`](kp_average_binary.md).

    Note
    ----
    This feature requires the attributes `k11` and `k22` to be defined.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.

    Returns
    -------
    float | FloatArray
        Average propagation rate coefficient. Unit = L/(mol·s)
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')
    return kp_average_binary(f1, *self.ri(f1), *self.kii(f1, T, Tunit))

plot

plot(
    kind: Literal["drift", "kp", "Mayo", "triads"],
    show: Literal["auto", "all", "data", "model"] = "auto",
    M: Literal[1, 2] = 1,
    f0: Optional[Union[float, FloatVectorLike]] = None,
    T: Optional[float] = None,
    Tunit: Literal["C", "K"] = "K",
    title: Optional[str] = None,
    axes: Optional[Axes] = None,
    return_objects: bool = False,
) -> Optional[tuple[Optional[Figure], Axes]]

Generate a plot of instantaneous copolymer composition, monomer composition drift, or average propagation rate coefficient.

PARAMETER	DESCRIPTION
kind	Kind of plot to be generated. TYPE: Literal['drift', 'kp', 'Mayo', 'triads']
show	What informatation is to be plotted. TYPE: Literal['auto', 'all', 'data', 'model'] DEFAULT: 'auto'
M	Index of the monomer to be used in input argument f0 and in output results. Specifically, if M=i, then f0 stands for \(f_i(0)\) and plots will be generated in terms of \(f_i\) and \(F_i\). TYPE: Literal[1, 2] DEFAULT: 1
f0	Initial monomer composition, \(f_i(0)\), as required for a monomer composition drift plot. TYPE: float | FloatVectorLike | None DEFAULT: None
T	Temperature. Unit = Tunit. TYPE: float | None DEFAULT: None
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
title	Title of plot. If None, a default title with the monomer names will be used. TYPE: str | None DEFAULT: None
axes	Matplotlib Axes object. TYPE: Axes | None DEFAULT: None
return_objects	If True, the Figure and Axes objects are returned (for saving or further manipulations). TYPE: bool DEFAULT: False

RETURNS	DESCRIPTION
tuple[Figure | None, Axes] | None	Figure and Axes objects if return_objects is True.

Source code in src/polykin/copolymerization/terminal.py

def plot(self,
         kind: Literal['drift', 'kp', 'Mayo', 'triads'],
         show: Literal['auto', 'all', 'data', 'model'] = 'auto',
         M: Literal[1, 2] = 1,
         f0: Optional[Union[float, FloatVectorLike]] = None,
         T: Optional[float] = None,
         Tunit: Literal['C', 'K'] = 'K',
         title: Optional[str] = None,
         axes: Optional[Axes] = None,
         return_objects: bool = False
         ) -> Optional[tuple[Optional[Figure], Axes]]:
    r"""Generate a plot of instantaneous copolymer composition, monomer
    composition drift, or average propagation rate coefficient.

    Parameters
    ----------
    kind : Literal['drift', 'kp', 'Mayo', 'triads']
        Kind of plot to be generated.
    show : Literal['auto', 'all', 'data', 'model']
        What informatation is to be plotted.
    M : Literal[1, 2]
        Index of the monomer to be used in input argument `f0` and in
        output results. Specifically, if `M=i`, then `f0` stands for
        $f_i(0)$ and plots will be generated in terms of $f_i$ and $F_i$.
    f0 : float | FloatVectorLike | None
        Initial monomer composition, $f_i(0)$, as required for a monomer
        composition drift plot.
    T : float | None
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    title : str | None
        Title of plot. If `None`, a default title with the monomer names
        will be used.
    axes : Axes | None
        Matplotlib Axes object.
    return_objects : bool
        If `True`, the Figure and Axes objects are returned (for saving or
        further manipulations).

    Returns
    -------
    tuple[Figure | None, Axes] | None
        Figure and Axes objects if return_objects is `True`.
    """
    check_in_set(M, {1, 2}, 'M')
    check_in_set(kind, {'Mayo', 'kp', 'drift', 'triads'}, 'kind')
    check_in_set(show, {'auto', 'all', 'data', 'model'}, 'show')

    label_model = None
    if axes is None:
        fig, ax = plt.subplots()
        if title is None:
            titles = {'Mayo': "Mayo-Lewis diagram",
                      'drift': "Monomer composition drift",
                      'kp': "Average propagation coefficient",
                      'triads': "Triad fractions"}
            title = titles[kind] + f" {self.M1}(1)-{self.M2}(2)"
        if title:
            fig.suptitle(title)
    else:
        ax = axes
        fig = None
        if self.name:
            label_model = self.name

    unit_range = (0., 1.)
    npoints = 1000
    Mname = self.M1 if M == 1 else self.M2
    ndataseries = 0

    if show == 'auto':
        if not self.data:
            show = 'model'
        else:
            show = 'all'

    if kind == 'Mayo':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(fr"$F_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        ax.plot(unit_range, unit_range, color='black', linewidth=0.5)

        if show in {'model', 'all'}:
            x = np.linspace(*unit_range, npoints, dtype=np.float64)
            y = self.F1(x)
            if M == 2:
                x[:] = 1. - x
                y[:] = 1. - y  # type: ignore
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:
            for ds in self.data:
                if not isinstance(ds, MayoDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('F', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

    elif kind == 'drift':

        ax.set_xlabel(r"Total molar monomer conversion, $x$")
        ax.set_ylabel(fr"$f_{M}$")
        ax.set_xlim(*unit_range)
        ax.set_ylim(*unit_range)

        if show == 'model':

            if f0 is None:
                raise ValueError("`f0` is required for a `drift` plot.")
            else:
                if isinstance(f0, (int, float)):
                    f0 = [f0]
                f0 = np.array(f0, dtype=np.float64)
                check_bounds(f0, *unit_range, 'f0')

            x = np.linspace(*unit_range, 1000)
            if M == 2:
                f0[:] = 1. - f0
            y = self.drift(f0, x)
            if M == 2:
                y[:] = 1. - y
            if y.ndim == 1:
                y = y[np.newaxis, :]
            for i in range(y.shape[0]):
                ax.plot(x, y[i, :], label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, DriftDataset):
                    continue
                x = ds.getvar('x')
                y = ds.getvar('f', Mname)
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)

                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.drift(ds.getvar('f', self.M1)[0], x)
                    if M == 2:
                        y[:] = 1. - y
                    ax.plot(x, y, label=label_model)

    elif kind == 'kp':

        ax.set_xlabel(fr"$f_{M}$")
        ax.set_ylabel(r"$\bar{k}_p$")
        ax.set_xlim(*unit_range)

        if show == 'model':

            if T is None:
                raise ValueError("`T` is required for a `kp` plot.")

            x = np.linspace(*unit_range, npoints)
            y = self.kp(x, T, Tunit)
            if M == 2:
                x[:] = 1. - x
            ax.plot(x, y, label=label_model)

        if show in {'data', 'all'}:

            for ds in self.data:
                if not isinstance(ds, kpDataset):
                    continue
                x = ds.getvar('f', Mname)
                y = ds.getvar('kp')
                ndataseries += 1
                ax.scatter(x, y,
                           label=ds.name if ds.name else None)
                if show == 'all':
                    x = np.linspace(*unit_range, npoints)
                    y = self.kp(x, ds.T, ds.Tunit)
                    if M == 2:
                        x[:] = 1. - x
                    ax.plot(x, y, label=label_model)

    ax.grid(True)

    if axes is not None or ndataseries:
        ax.legend(bbox_to_anchor=(1.05, 1.), loc="upper left")

    if return_objects:
        return (fig, ax)

ri

ri(
    _,
) -> tuple[
    Union[float, FloatArray], Union[float, FloatArray]
]

Return the evaluated reactivity ratios at the given conditions.

Source code in src/polykin/copolymerization/terminal.py

def ri(self,
       _
       ) -> tuple[Union[float, FloatArray], Union[float, FloatArray]]:
    return (self.r1, self.r2)

sequence

sequence(
    f1: Union[float, FloatArrayLike],
    k: Optional[Union[int, IntArrayLike]] = None,
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous sequence length probability or the number-average sequence length.

For a binary system, the probability of finding \(k\) consecutive units of monomer \(i\) in a chain is:

\[ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} \]

and the corresponding number-average sequence length is:

\[ \bar{S}_i = \sum_k k S_{i,k} = \frac{1}{1 - P_{ii}} \]

where \(P_{ii}\) is the transition probability \(i \rightarrow i\), which is a function of the monomer composition and the reactivity ratios.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 177.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
k	Sequence length, i.e., number of consecutive units in a chain. If None, the number-average sequence length will be computed. TYPE: int | IntArrayLike | None DEFAULT: None

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	If k is None, the number-average sequence lengths, {'1': \(\bar{S}_1\), '2': \(\bar{S}_2\)}. Otherwise, the sequence probabilities, {'1': \(S_{1,k}\), '2': \(S_{2,k}\)}.

Source code in src/polykin/copolymerization/terminal.py

def sequence(self,
             f1: Union[float, FloatArrayLike],
             k: Optional[Union[int, IntArrayLike]] = None,
             ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous sequence length probability or the
    number-average sequence length.

    For a binary system, the probability of finding $k$ consecutive units
    of monomer $i$ in a chain is:

    $$ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} $$

    and the corresponding number-average sequence length is:

    $$ \bar{S}_i = \sum_k k S_{i,k} = \frac{1}{1 - P_{ii}} $$

    where $P_{ii}$ is the transition probability $i \rightarrow i$, which
    is a function of the monomer composition and the reactivity ratios.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 177.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    k : int | IntArrayLike | None
        Sequence length, i.e., number of consecutive units in a chain.
        If `None`, the number-average sequence length will be computed.

    Returns
    -------
    dict[str, float | FloatArray]
        If `k is None`, the number-average sequence lengths,
        {'1': $\bar{S}_1$, '2': $\bar{S}_2$}. Otherwise, the
        sequence probabilities, {'1': $S_{1,k}$, '2': $S_{2,k}$}.
    """

    P11, _, _, P22 = self.transitions(f1).values()

    if k is None:
        result = {str(i + 1): 1/(1 - P + eps)
                  for i, P in enumerate([P11, P22])}
    else:
        if isinstance(k, (list, tuple)):
            k = np.array(k, dtype=np.int32)
        result = {str(i + 1): (1. - P)*P**(k - 1)
                  for i, P in enumerate([P11, P22])}

    return result

transitions

transitions(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous transition probabilities.

For a binary system, the transition probabilities are given by:

\[\begin{aligned} P_{ii} &= \frac{r_i f_i}{r_i f_i + (1 - f_i)} \\ P_{ij} &= 1 - P_{ii} \end{aligned}\]

where \(i,j=1,2, i \neq j\).

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 178.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Transition probabilities, {'11': \(P_{11}\), '12': \(P_{12}\), ... }.

Source code in src/polykin/copolymerization/terminal.py

def transitions(self,
                f1: Union[float, FloatArrayLike]
                ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous transition probabilities.

    For a binary system, the transition probabilities are given by:

    \begin{aligned}
        P_{ii} &= \frac{r_i f_i}{r_i f_i + (1 - f_i)} \\
        P_{ij} &= 1 - P_{ii}
    \end{aligned}

    where $i,j=1,2, i \neq j$.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 178.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Transition probabilities, {'11': $P_{11}$, '12': $P_{12}$, ... }.
    """

    if isinstance(f1, (list, tuple)):
        f1 = np.array(f1, dtype=np.float64)
    check_bounds(f1, 0., 1., 'f1')

    f2 = 1. - f1
    r1 = self.r1
    r2 = self.r2
    P11 = r1*f1/(r1*f1 + f2)
    P22 = r2*f2/(r2*f2 + f1)
    P12 = 1. - P11
    P21 = 1. - P22

    result = {'11': P11,
              '12': P12,
              '21': P21,
              '22': P22}

    return result

triads

triads(
    f1: Union[float, FloatArrayLike]
) -> dict[str, Union[float, FloatArray]]

Calculate the instantaneous triad fractions.

For a binary system, the triad fractions are given by:

\[\begin{aligned} A_{iii} &= F_i P_{ii}^2 \\ A_{iij} &= 2 F_i P_{ii} P_{ij} \\ A_{jij} &= F_i P_{ij}^2 \end{aligned}\]

where \(P_{ij}\) is the transition probability \(i \rightarrow j\) and \(F_i\) is the instantaneous copolymer composition.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 179.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike

RETURNS	DESCRIPTION
dict[str, float | FloatArray]	Triad fractions, {'111': \(A_{111}\), '112': \(A_{112}\), '212': \(A_{212}\), ... }.

Source code in src/polykin/copolymerization/terminal.py

def triads(self,
           f1: Union[float, FloatArrayLike]
           ) -> dict[str, Union[float, FloatArray]]:
    r"""Calculate the instantaneous triad fractions.

    For a binary system, the triad fractions are given by:

    \begin{aligned}
        A_{iii} &= F_i P_{ii}^2 \\
        A_{iij} &= 2 F_i P_{ii} P_{ij} \\
        A_{jij} &= F_i P_{ij}^2
    \end{aligned}

    where $P_{ij}$ is the transition probability $i \rightarrow j$ and
    $F_i$ is the instantaneous copolymer composition.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 179.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.

    Returns
    -------
    dict[str, float | FloatArray]
        Triad fractions,
        {'111': $A_{111}$, '112': $A_{112}$, '212': $A_{212}$, ... }.
    """

    P11, P12, P21, P22 = self.transitions(f1).values()

    F1 = P21/(P12 + P21)
    F2 = 1. - F1

    A111 = F1*P11**2
    A112 = F1*2*P11*P12
    A212 = F1*P12**2

    A222 = F2*P22**2
    A221 = F2*2*P22*P21
    A121 = F2*P21**2

    result = {'111': A111,
              '112': A112,
              '212': A212,
              '222': A222,
              '221': A221,
              '121': A121}

    return result

convert_Qe_to_r

convert_Qe_to_r(
    Qe_values: list[tuple[float, float]]
) -> FloatSquareMatrix

Convert Q-e values to reactivity ratios.

According to the Q-e scheme proposed by Alfrey and Price, the reactivity ratios of the terminal model can be estimated using the relationship:

\[ r_{ij} = \frac{Q_i}{Q_j}\exp{\left(-e_i(e_i -e_j)\right)} \]

where \(Q_i\) and \(e_i\) are monomer-specific constants, and \(r_{ij}=k_{ii}/k_{ij}\) is the multicomponent reactivity ratio matrix.

References

• T Alfrey, CC Price. J. Polym. Sci., 1947, 2: 101-106.

PARAMETER	DESCRIPTION
Qe_values	List (N) of Q-e values. TYPE: list[tuple[float, float]]

RETURNS	DESCRIPTION
FloatSquareMatrix(N, N)	Reactivity ratio matrix.

Examples:

Estimate the reactivity ratio matrix for styrene (1), methyl methacrylate (2), and vinyl acetate(3) using Q-e values from the literature.

>>> from polykin.copolymerization import convert_Qe_to_r
>>>
>>> Qe1 = (1.0, -0.80)    # Sty
>>> Qe2 = (0.78, 0.40)    # MMA
>>> Qe3 = (0.026, -0.88)  # VAc
>>>
>>> convert_Qe_to_r([Qe1, Qe2, Qe3])
array([[1.00000000e+00, 4.90888315e-01, 4.10035538e+01],
       [4.82651046e-01, 1.00000000e+00, 1.79788736e+01],
       [2.42325444e-02, 1.08066091e-02, 1.00000000e+00]])

Source code in src/polykin/copolymerization/multicomponent.py

def convert_Qe_to_r(Qe_values: list[tuple[float, float]]
                    ) -> FloatSquareMatrix:
    r"""Convert Q-e values to reactivity ratios.

    According to the Q-e scheme proposed by Alfrey and Price, the reactivity
    ratios of the terminal model can be estimated using the relationship:

    $$ r_{ij} = \frac{Q_i}{Q_j}\exp{\left(-e_i(e_i -e_j)\right)} $$

    where $Q_i$ and $e_i$ are monomer-specific constants, and
    $r_{ij}=k_{ii}/k_{ij}$ is the multicomponent reactivity ratio matrix.

    **References**

    *   T Alfrey, CC Price. J. Polym. Sci., 1947, 2: 101-106.

    Parameters
    ----------
    Qe_values : list[tuple[float, float]]
        List (N) of Q-e values.

    Returns
    -------
    FloatSquareMatrix (N, N)
        Reactivity ratio matrix.

    Examples
    --------
    Estimate the reactivity ratio matrix for styrene (1), methyl
    methacrylate (2), and vinyl acetate(3) using Q-e values from the
    literature.

    >>> from polykin.copolymerization import convert_Qe_to_r
    >>>
    >>> Qe1 = (1.0, -0.80)    # Sty
    >>> Qe2 = (0.78, 0.40)    # MMA
    >>> Qe3 = (0.026, -0.88)  # VAc
    >>>
    >>> convert_Qe_to_r([Qe1, Qe2, Qe3])
    array([[1.00000000e+00, 4.90888315e-01, 4.10035538e+01],
           [4.82651046e-01, 1.00000000e+00, 1.79788736e+01],
           [2.42325444e-02, 1.08066091e-02, 1.00000000e+00]])

    """

    Q = [x[0] for x in Qe_values]
    e = [x[1] for x in Qe_values]
    N = len(Q)

    r = np.empty((N, N))
    for i in range(N):
        for j in range(N):
            r[i, j] = Q[i]/Q[j]*exp(-e[i]*(e[i] - e[j]))

    return r

fit_Finemann_Ross

fit_Finemann_Ross(
    f1: FloatVectorLike, F1: FloatVectorLike
) -> tuple[float, float]

Fit binary copolymer composition data using the Finemann-Ross method.

\[ \left(\frac{x(y-1)}{y}\right) = -r_2 + r_1 \left(\frac{x^2}{y}\right) \]

where \(x = f_1/(1 - f_1)\), \(y = F_1/(1 - F_1)\), \(r_i\) are the reactivity ratios, \(f_1\) is the monomer composition, and \(F_1\) is the instantaneous copolymer composition.

Reference

• Fineman, M.; Ross, S. D. J. Polymer Sci. 1950, 5, 259.

Note

The Finemann-Ross method relies on a linearization procedure that can lead to significant statistical bias. The method is provided for its historical significance, but should no longer be used for fitting reactivity ratios.

PARAMETER	DESCRIPTION
f1	Vector of molar fraction of M1, \(f_1\). TYPE: FloatVectorLike
F1	Vector of instantaneous copolymer composition of M1, \(F_1\). TYPE: FloatVectorLike

RETURNS	DESCRIPTION
tuple[float, float]	Point estimates of the reactivity ratios \((r_1, r_2)\).

See also

• fit_copo_data: alternative (recommended) method.

Examples:

>>> from polykin.copolymerization.fitting import fit_Finemann_Ross
>>>
>>> f1 = [0.186, 0.299, 0.527, 0.600, 0.700, 0.798]
>>> F1 = [0.196, 0.279, 0.415, 0.473, 0.542, 0.634]
>>>
>>> r1, r2 = fit_Finemann_Ross(f1, F1)
>>> print(f"r1={r1:.2f}, r2={r2:.2f}")
r1=0.27, r2=0.84

Source code in src/polykin/copolymerization/fitting.py

def fit_Finemann_Ross(f1: FloatVectorLike,
                      F1: FloatVectorLike
                      ) -> tuple[float, float]:
    r"""Fit binary copolymer composition data using the Finemann-Ross method.

    $$ \left(\frac{x(y-1)}{y}\right) = -r_2 + r_1 \left(\frac{x^2}{y}\right) $$

    where $x = f_1/(1 - f_1)$, $y = F_1/(1 - F_1)$, $r_i$ are the reactivity
    ratios, $f_1$ is the monomer composition, and $F_1$ is the instantaneous
    copolymer composition.

    **Reference**

    *   Fineman, M.; Ross, S. D. J. Polymer Sci. 1950, 5, 259.

    Note
    ----
    The Finemann-Ross method relies on a linearization procedure that can lead
    to significant statistical bias. The method is provided for its historical
    significance, but should no longer be used for fitting reactivity ratios.


    Parameters
    ----------
    f1 : FloatVectorLike
        Vector of molar fraction of M1, $f_1$.
    F1 : FloatVectorLike
        Vector of instantaneous copolymer composition of M1, $F_1$.

    Returns
    -------
    tuple[float, float]
        Point estimates of the reactivity ratios $(r_1, r_2)$.

    See also
    --------
    * [`fit_copo_data`](fit_copo_data.md): alternative
      (recommended) method.  

    Examples
    --------
    >>> from polykin.copolymerization.fitting import fit_Finemann_Ross
    >>>
    >>> f1 = [0.186, 0.299, 0.527, 0.600, 0.700, 0.798]
    >>> F1 = [0.196, 0.279, 0.415, 0.473, 0.542, 0.634]
    >>>
    >>> r1, r2 = fit_Finemann_Ross(f1, F1)
    >>> print(f"r1={r1:.2f}, r2={r2:.2f}")
    r1=0.27, r2=0.84

    """

    f1 = np.asarray(f1)
    F1 = np.asarray(F1)

    # Variable transformation
    x = f1/(1. - f1)
    y = F1/(1. - F1)
    H = x**2/y
    G = x*(y - 1.)/y

    # Linear regression
    solution = linregress(H, G)
    r1 = solution.slope  # type: ignore
    r2 = - solution.intercept  # type: ignore

    return (r1, r2)

fit_copo_data

fit_copo_data(
    data_Ff: list[CopoDataset_Ff] = [],
    data_fx: list[CopoDataset_fx] = [],
    data_Fx: list[CopoDataset_Fx] = [],
    r_guess: tuple[float, float] = (1.0, 1.0),
    method: Literal["NLLS", "ODR"] = "NLLS",
    alpha: float = 0.05,
    plot_data: bool = True,
    JCR_linear: bool = True,
    JCR_exact: bool = False,
    JCR_npoints: int = 200,
    JCR_rtol: float = 0.01,
) -> CopoFitResult

Fit copolymer composition data and estimate reactivity ratios.

This function employs rigorous nonlinear algorithms to estimate the reactivity ratios from experimental copolymer composition data of type \(F(f)\), \(f(x;f_0)\), and \(F(x,f_0)\).

The fitting is performed using one of two methods: nonlinear least squares (NLLS) or orthogonal distance regression (ODR). NLLS considers only observational errors in the dependent variable, whereas ODR takes into account observational errors in both the dependent and independent variables. Although the ODR method is statistically more general, it is also more complex and can (at present) only be used for fitting \(F(f)\) data. Whenever composition drift data is provided, NLLS must be utilized.

The joint confidence region (JCR) of the reactivity ratios is generated using approximate (linear) and/or exact methods. In most cases, the linear method should be sufficiently accurate. Nonetheless, for these types of fits, the exact method is computationally inexpensive, making it perhaps a preferable choice.

Reference

• Van Herk, A.M. and Dröge, T. (1997), Nonlinear least squares fitting applied to copolymerization modeling. Macromol. Theory Simul., 6: 1263-1276.
• Boggs, Paul T., et al. "Algorithm 676: ODRPACK: software for weighted orthogonal distance regression." ACM Transactions on Mathematical Software (TOMS) 15.4 (1989): 348-364.

PARAMETER	DESCRIPTION
data_Ff	F(f) instantaneous composition datasets. TYPE: list[CopoDataset_Ff] DEFAULT: []
data_fx	f(x) composition drift datasets. TYPE: list[CopoDataset_fx] DEFAULT: []
data_Fx	F(x) composition drift datasets TYPE: list[CopoDataset_Fx] DEFAULT: []
r_guess	Initial guess for the reactivity ratios. TYPE: tuple[float, float] DEFAULT: (1.0, 1.0)
method	Optimization method. NLLS for nonlinear least squares or ODR for orthogonal distance regression. The ODR method is only available for F(f) data. TYPE: Literal['NLLS', 'ODR'] DEFAULT: 'NLLS'
alpha	Significance level. TYPE: float DEFAULT: 0.05
plot_data	If True, comparisons between experimental and fitted data will be plotted. TYPE: bool DEFAULT: True
JCR_linear	If True, the JCR will be computed using the linear method. TYPE: bool DEFAULT: True
JCR_exact	If True, the JCR will be computed using the exact method. TYPE: bool DEFAULT: False
JCR_npoints	Number of points where the JCR is evaluated. The computational effort increases linearly with npoints. TYPE: int DEFAULT: 200
JCR_rtol	Relative tolerance for the determination of the JCR. The default value (1e-2) should be adequate in most cases, as it implies a 1% accuracy in the JCR coordinates. TYPE: float DEFAULT: 0.01

RETURNS	DESCRIPTION
CopoFitResult	Dataclass with all fit results.

See also

• confidence_ellipse: linear method used to calculate the joint confidence region.
• confidence_region: exact method used to calculate the joint confidence region.
• fit_Finemann_Ross: alternative method.

Source code in src/polykin/copolymerization/fitting.py

def fit_copo_data(data_Ff: list[CopoDataset_Ff] = [],
                  data_fx: list[CopoDataset_fx] = [],
                  data_Fx: list[CopoDataset_Fx] = [],
                  r_guess: tuple[float, float] = (1.0, 1.0),
                  method: Literal['NLLS', 'ODR'] = 'NLLS',
                  alpha: float = 0.05,
                  plot_data: bool = True,
                  JCR_linear: bool = True,
                  JCR_exact: bool = False,
                  JCR_npoints: int = 200,
                  JCR_rtol: float = 1e-2
                  ) -> CopoFitResult:
    r"""Fit copolymer composition data and estimate reactivity ratios.

    This function employs rigorous nonlinear algorithms to estimate the
    reactivity ratios from experimental copolymer composition data of type
    $F(f)$, $f(x;f_0)$, and $F(x,f_0)$. 

    The fitting is performed using one of two methods: nonlinear least squares
    (NLLS) or orthogonal distance regression (ODR). NLLS considers only
    observational errors in the dependent variable, whereas ODR takes into
    account observational errors in both the dependent and independent
    variables. Although the ODR method is statistically more general, it is
    also more complex and can (at present) only be used for fitting $F(f)$
    data. Whenever composition drift data is provided, NLLS must be utilized.

    The joint confidence region (JCR) of the reactivity ratios is generated
    using approximate (linear) and/or exact methods. In most cases, the linear
    method should be sufficiently accurate. Nonetheless, for these types of
    fits, the exact method is computationally inexpensive, making it perhaps a
    preferable choice.

    **Reference**

    *   Van Herk, A.M. and Dröge, T. (1997), Nonlinear least squares fitting
        applied to copolymerization modeling. Macromol. Theory Simul.,
        6: 1263-1276.
    *   Boggs, Paul T., et al. "Algorithm 676: ODRPACK: software for weighted
        orthogonal distance regression." ACM Transactions on Mathematical
        Software (TOMS) 15.4 (1989): 348-364.

    Parameters
    ----------
    data_Ff : list[CopoDataset_Ff]
        F(f) instantaneous composition datasets.
    data_fx : list[CopoDataset_fx]
        f(x) composition drift datasets.
    data_Fx : list[CopoDataset_Fx]
        F(x) composition drift datasets
    r_guess : tuple[float, float]
        Initial guess for the reactivity ratios.
    method : Literal['NLLS', 'ODR']
        Optimization method. `NLLS` for nonlinear least squares or `ODR` for
        orthogonal distance regression. The `ODR` method is only available for
        F(f) data.
    alpha : float
        Significance level.
    plot_data : bool
        If `True`, comparisons between experimental and fitted data will be
        plotted.
    JCR_linear : bool, optional
        If `True`, the JCR will be computed using the linear method.
    JCR_exact : bool, optional
        If `True`, the JCR will be computed using the exact method.
    JCR_npoints : int
        Number of points where the JCR is evaluated. The computational effort
        increases linearly with `npoints`.
    JCR_rtol : float
        Relative tolerance for the determination of the JCR. The default value
        (1e-2) should be adequate in most cases, as it implies a 1% accuracy in
        the JCR coordinates. 

    Returns
    -------
    CopoFitResult
        Dataclass with all fit results.

    See also
    --------
    * [`confidence_ellipse`](../math/confidence_ellipse.md): linear method
      used to calculate the joint confidence region.
    * [`confidence_region`](../math/confidence_region.md): exact method
      used to calculate the joint confidence region.
    * [`fit_Finemann_Ross`](fit_Finemann_Ross.md): alternative method.  

    """

    # Calculate and check 'ndata'
    npar = 2
    ndata = 0
    for dataset in data_Ff:
        ndata += len(dataset.f1)
    for dataset in data_fx:
        ndata += len(dataset.x)
    for dataset in data_Fx:
        ndata += len(dataset.x)
    if ndata < npar:
        raise FitError("Not enough data to estimate reactivity ratios.")

    # Check method choice
    if method == 'ODR' and (len(data_fx) > 0 or len(data_Fx) > 0):
        raise FitError("ODR method not implemented for drift data.")

    # Fit data
    if method == 'NLLS':
        r_opt, cov, sse = _fit_copo_NLLS(data_Ff, data_fx, data_Fx, r_guess,
                                         ndata)
    elif method == 'ODR':
        r_opt, cov, sse = _fit_copo_ODR(data_Ff, r_guess)
    else:
        raise ValueError(f"Method {method} is not valid.")

    # Standard parameter errors and confidence intervals
    if ndata > npar and cov is not None:
        se_r = np.sqrt(np.diag(cov))
        ci_r = se_r*tdist.ppf(1. - alpha/2, ndata - npar)
    else:
        se_r = [np.nan, np.nan]
        ci_r = [np.nan, np.nan]

    # Plot data vs model
    plots = {}
    if plot_data:
        xmax = 1.
        npoints = 500
        # Plot F(f) data
        if data_Ff:
            plots['Ff'] = plt.subplots()
            ax = plots['Ff'][1]
            ax.set_xlabel(r"$f_1$")
            ax.set_ylabel(r"$F_1$")
            ax.set_xlim(0., 1.)
            ax.set_ylim(0., 1.)
            for dataset in data_Ff:
                ax.scatter(dataset.f1, dataset.F1, label=dataset.name)
            f1 = np.linspace(0., 1., npoints)
            F1_est = inst_copolymer_binary(f1, *r_opt)
            ax.plot(f1, F1_est)
            ax.legend(loc="best")

        # Plot f(x) data
        if data_fx:
            plots['fx'] = plt.subplots()
            ax = plots['fx'][1]
            ax.set_xlabel(r"$x$")
            ax.set_ylabel(r"$f_1$")
            ax.set_xlim(0., 1.)
            ax.set_ylim(0., 1.)
            x = np.linspace(0., xmax, npoints)
            for dataset in data_fx:
                ax.scatter(dataset.x, dataset.f1, label=dataset.name)
                f1_est = monomer_drift_binary(dataset.f10, x, *r_opt)
                ax.plot(x, f1_est)
            ax.legend(loc="best")

        # Plot F(x) data
        if data_Fx:
            plots['Fx'] = plt.subplots()
            ax = plots['Fx'][1]
            ax.set_xlabel(r"$x$")
            ax.set_ylabel(r"$F_1$")
            ax.set_xlim(0., 1.)
            ax.set_ylim(0., 1.)
            x = np.linspace(0., xmax, npoints)
            for dataset in data_Fx:
                f10 = dataset.f10
                ax.scatter(dataset.x, dataset.F1, label=dataset.name)
                f1_est = monomer_drift_binary(f10, x, *r_opt)
                F1_est = f1_est + (f10 - f1_est)/(x + 1e-10)
                F1_est[0] = inst_copolymer_binary(f10, *r_opt)
                ax.plot(x, F1_est)
            ax.legend(loc="best")

        # Parity plot
        plots['parity'] = plt.subplots()
        ax = plots['parity'][1]
        ax.set_xlabel("Observed value")
        ax.set_ylabel("Predicted value")
        ax.set_xlim(0., 1.)
        ax.set_ylim(0., 1.)
        ax.plot([0., 1.], [0., 1.], color='black', linewidth=0.5)
        for dataset in data_Ff:
            F1_est = inst_copolymer_binary(dataset.f1, *r_opt)
            ax.scatter(dataset.F1, F1_est, label=dataset.name)
        for dataset in data_fx:
            f1_est = monomer_drift_binary(dataset.f10, dataset.x, *r_opt)
            ax.scatter(dataset.f1, f1_est, label=dataset.name)
        for dataset in data_Fx:
            f1_est = monomer_drift_binary(dataset.f10, dataset.x, *r_opt)
            F1_est = f1_est + (dataset.f10 - f1_est)/(dataset.x + 1e-10)
            ax.scatter(dataset.F1, F1_est, label=dataset.name)
        ax.legend(loc="best")

    # Joint confidence region
    if (JCR_linear or JCR_exact) and cov is not None:

        plots['JCR'] = plt.subplots()
        ax = plots['JCR'][1]
        ax.set_xlabel(r"$r_1$")
        ax.set_ylabel(r"$r_2$")
        colors = ['tab:blue', 'tab:orange']
        idx = 0

        if JCR_linear:
            confidence_ellipse(ax,
                               center=r_opt,
                               cov=cov,
                               ndata=ndata,
                               alpha=alpha,
                               color=colors[idx], label='linear JCR')
            idx += 1

        if JCR_exact:
            jcr = confidence_region(center=r_opt,
                                    sse=sse,
                                    ndata=ndata,
                                    alpha=alpha,
                                    width=2*ci_r[0],
                                    rtol=JCR_rtol,
                                    npoints=JCR_npoints)

            # ax.scatter(r1, r2, c=colors[idx], s=5)
            ax.plot(*jcr, color=colors[idx], label='exact JCR')

        ax.legend(loc="best")

    result = CopoFitResult(method=method,
                           r1=r_opt[0], r2=r_opt[1],
                           alpha=alpha,
                           ci_r1=ci_r[0], ci_r2=ci_r[1],
                           se_r1=se_r[0], se_r2=se_r[1],
                           cov=cov,
                           plots=plots)

    return result

inst_copolymer_binary

inst_copolymer_binary(
    f1: Union[float, FloatArrayLike],
    r1: Union[float, FloatArray],
    r2: Union[float, FloatArray],
) -> Union[float, FloatArray]

Calculate the instantaneous copolymer composition using the Mayo-Lewis equation.

For a binary system, the instantaneous copolymer composition \(F_i\) is related to the comonomer composition \(f_i\) by:

\[ F_1=\frac{r_1 f_1^2 + f_1 f_2}{r_1 f_1^2 + 2 f_1 f_2 + r_2 f_2^2} \]

where \(r_i\) are the reactivity ratios. Although the equation is written using terminal model notation, it is equally applicable in the frame of the penultimate model if \(r_i \rightarrow \bar{r}_i\).

References

• Mayo & Lewis (1944)

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
r1	Reactivity ratio of M1, \(r_1\) or \(\bar{r}_1\). TYPE: float | FloatArray
r2	Reactivity ratio of M2, \(r_2\) or \(\bar{r}_2\). TYPE: float | FloatArray

RETURNS	DESCRIPTION
float | FloatArray	Instantaneous copolymer composition, \(F_1\).

See also

• inst_copolymer_ternary: specific method for terpolymer systems.
• inst_copolymer_multi: generic method for multicomponent systems.

Examples:

>>> from polykin.copolymerization import inst_copolymer_binary

An example with f1 as scalar.

>>> F1 = inst_copolymer_binary(f1=0.5, r1=0.16, r2=0.70)
>>> print(f"F1 = {F1:.3f}")
F1 = 0.406

An example with f1 as list.

>>> F1 = inst_copolymer_binary(f1=[0.2, 0.6, 0.8], r1=0.16, r2=0.70)
>>> F1
array([0.21487603, 0.45812808, 0.58259325])

Source code in src/polykin/copolymerization/binary.py

@jit
def inst_copolymer_binary(f1: Union[float, FloatArrayLike],
                          r1: Union[float, FloatArray],
                          r2: Union[float, FloatArray]
                          ) -> Union[float, FloatArray]:
    r"""Calculate the instantaneous copolymer composition using the
    [Mayo-Lewis](https://en.wikipedia.org/wiki/Mayo%E2%80%93Lewis_equation)
     equation.

    For a binary system, the instantaneous copolymer composition $F_i$ is
    related to the comonomer composition $f_i$ by:

    $$ F_1=\frac{r_1 f_1^2 + f_1 f_2}{r_1 f_1^2 + 2 f_1 f_2 + r_2 f_2^2} $$

    where $r_i$ are the reactivity ratios. Although the equation is written
    using terminal model notation, it is equally applicable in the frame of the
    penultimate model if $r_i \rightarrow \bar{r}_i$.

    **References**

    *   [Mayo & Lewis (1944)](https://doi.org/10.1021/ja01237a052)

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    r1 : float | FloatArray
        Reactivity ratio of M1, $r_1$ or $\bar{r}_1$.
    r2 : float | FloatArray
        Reactivity ratio of M2, $r_2$ or $\bar{r}_2$.

    Returns
    -------
    float | FloatArray
        Instantaneous copolymer composition, $F_1$.

    See also
    --------
    * [`inst_copolymer_ternary`](inst_copolymer_ternary.md):
      specific method for terpolymer systems.
    * [`inst_copolymer_multi`](inst_copolymer_multi.md):
      generic method for multicomponent systems.

    Examples
    --------
    >>> from polykin.copolymerization import inst_copolymer_binary

    An example with f1 as scalar.
    >>> F1 = inst_copolymer_binary(f1=0.5, r1=0.16, r2=0.70)
    >>> print(f"F1 = {F1:.3f}")
    F1 = 0.406

    An example with f1 as list.
    >>> F1 = inst_copolymer_binary(f1=[0.2, 0.6, 0.8], r1=0.16, r2=0.70)
    >>> F1
    array([0.21487603, 0.45812808, 0.58259325])

    """
    f1 = np.asarray(f1)
    f2 = 1 - f1
    return (r1*f1**2 + f1*f2)/(r1*f1**2 + 2*f1*f2 + r2*f2**2)

inst_copolymer_multi

inst_copolymer_multi(
    f: Optional[FloatVectorLike],
    r: Optional[FloatSquareMatrix],
    P: Optional[FloatSquareMatrix] = None,
) -> FloatVector

Calculate the instantaneous copolymer composition for a multicomponent system.

In a multicomponent system, the instantaneous copolymer composition \(F_i\) can be determined by solving the following set of linear algebraic equations:

\[ \begin{bmatrix} P_{11}-1 & P_{21} & \cdots & P_{N1} \\ P_{12} & P_{22}-1 & \cdots & P_{N2} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & \cdots & 1 \end{bmatrix} \begin{bmatrix} F_1 \\ F_2 \\ \vdots \\ F_N \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ \vdots \\ 1 \end{bmatrix} \]

where \(P_{ij}\) are the transition probabilitites, which can be computed from the instantaneous monomer composition and the reactivity matrix.

References

• H. K. Frensdorff, R. Pariser; Copolymerization as a Markov Chain. J. Chem. Phys. 1 November 1963; 39 (9): 2303-2309.

PARAMETER	DESCRIPTION
f	Vector of instantaneous monomer compositions, \(f_i\). TYPE: FloatVectorLike(N) | None
r	Matrix of reactivity ratios, \(r_{ij}=k_{ii}/k_{ij}\). TYPE: FloatSquareMatrix(N, N) | None
P	Matrix of transition probabilities, \(P_{ij}\). If None, it will be computed internally. When calculating other quantities (e.g., sequence lengths, tuples) that also depend on \(P\), it is more efficient to precompute \(P\) once and use it in all cases. TYPE: FloatSquareMatrix(N, N) | None DEFAULT: None

RETURNS	DESCRIPTION
FloatVector(N)	Vector of instantaneous copolymer compositions, \(F_i\).

See also

• inst_copolymer_binary: specific method for binary systems.
• inst_copolymer_ternary: specific method for terpolymer systems.
• monomer_drift_multi: monomer composition drift.
• transitions_multi: instantaneous transition probabilities.

Examples:

>>> from polykin.copolymerization import inst_copolymer_multi
>>> import numpy as np

Define the reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate the instantaneous copolymer composition at f1=0.5, f2=0.3, f3=0.2.

>>> f = [0.5, 0.3, 0.2]
>>> F = inst_copolymer_multi(f, r)
>>> F
array([0.32138111, 0.41041608, 0.26820282])

Source code in src/polykin/copolymerization/multicomponent.py

def inst_copolymer_multi(f: Optional[FloatVectorLike],
                         r: Optional[FloatSquareMatrix],
                         P: Optional[FloatSquareMatrix] = None
                         ) -> FloatVector:
    r"""Calculate the instantaneous copolymer composition for a multicomponent
    system.

    In a multicomponent system, the instantaneous copolymer composition $F_i$
    can be determined by solving the following set of linear algebraic
    equations:

    $$ \begin{bmatrix}
    P_{11}-1  & P_{21}   & \cdots  & P_{N1} \\
    P_{12}    & P_{22}-1 & \cdots  & P_{N2} \\
    \vdots    & \vdots   & \vdots  & \vdots \\
    1         & 1        & \cdots  & 1
    \end{bmatrix}
    \begin{bmatrix}
    F_1    \\
    F_2    \\
    \vdots \\
    F_N
    \end{bmatrix} =
    \begin{bmatrix}
    0      \\
    0      \\
    \vdots \\
    1
    \end{bmatrix} $$

    where $P_{ij}$ are the transition probabilitites, which can be computed
    from the instantaneous monomer composition and the reactivity matrix.

    **References**

    *   H. K. Frensdorff, R. Pariser; Copolymerization as a Markov Chain.
        J. Chem. Phys. 1 November 1963; 39 (9): 2303-2309.

    Parameters
    ----------
    f : FloatVectorLike (N) | None
        Vector of instantaneous monomer compositions, $f_i$.
    r : FloatSquareMatrix (N, N) | None
        Matrix of reactivity ratios, $r_{ij}=k_{ii}/k_{ij}$.
    P : FloatSquareMatrix (N, N) | None
        Matrix of transition probabilities, $P_{ij}$. If `None`, it will
        be computed internally. When calculating other quantities (e.g.,
        sequence lengths, tuples) that also depend on $P$, it is more efficient
        to precompute $P$ once and use it in all cases.

    Returns
    -------
    FloatVector (N)
        Vector of instantaneous copolymer compositions, $F_i$.

    See also
    --------
    * [`inst_copolymer_binary`](inst_copolymer_binary.md):
      specific method for binary systems.
    * [`inst_copolymer_ternary`](inst_copolymer_ternary.md):
      specific method for terpolymer systems.
    * [`monomer_drift_multi`](monomer_drift_multi.md):
      monomer composition drift.
    * [`transitions_multi`](transitions_multi.md):
      instantaneous transition probabilities.

    Examples
    --------
    >>> from polykin.copolymerization import inst_copolymer_multi
    >>> import numpy as np

    Define the reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate the instantaneous copolymer composition at f1=0.5, f2=0.3, f3=0.2.
    >>> f = [0.5, 0.3, 0.2]
    >>> F = inst_copolymer_multi(f, r)
    >>> F
    array([0.32138111, 0.41041608, 0.26820282])
    """

    if P is not None and (f is None and r is None):
        pass
    elif P is None and (f is not None and r is not None):
        P = transitions_multi(f, r)
    else:
        raise ValueError("Invalid combination of inputs.")

    N = P.shape[0]
    A = P.T - np.eye(N)
    A[-1, :] = 1.
    b = np.zeros(N)
    b[-1] = 1.
    F = np.linalg.solve(A, b)

    return F

inst_copolymer_ternary

inst_copolymer_ternary(
    f1: Union[float, FloatArrayLike],
    f2: Union[float, FloatArrayLike],
    r12: float,
    r21: float,
    r13: float,
    r31: float,
    r23: float,
    r32: float,
) -> tuple[
    Union[float, FloatArray],
    Union[float, FloatArray],
    Union[float, FloatArray],
]

Calculate the instantaneous copolymer composition for a ternary system.

In a ternary system, the instantaneous copolymer composition \(F_i\) is related to the monomer composition \(f_i\) by:

\[\begin{aligned} a &= \frac{f_1}{r_{21} r_{31}} + \frac{f_2}{r_{21} r_{32}} + \frac{f_3}{r_{31} r_{23}} \\ b &= f_1 + \frac{f_2}{r_{12}} + \frac{f_3}{r_{13}} \\ c &= \frac{f_1}{r_{12} r_{31}} + \frac{f_2}{r_{12} r_{32}} + \frac{f_3}{r_{13} r_{32}} \\ d &= f_2 + \frac{f_1}{r_{21}} + \frac{f_3}{r_{23}} \\ e &= \frac{f_1}{r_{13} r_{21}} + \frac{f_2}{r_{23} r_{12}} + \frac{f_3}{r_{13} r_{23}} \\ g &= f_3 + \frac{f_1}{r_{31}} + \frac{f_2}{r_{32}} \\ F_1 &= \frac{a b f_1}{a b f_1 + c d f_2 + e g f_3} \\ F_2 &= \frac{c d f_2}{a b f_1 + c d f_2 + e g f_3} \\ F_3 &= \frac{e g f_3}{a b f_1 + c d f_2 + e g f_3} \end{aligned}\]

where \(r_{ij}=k_{ii}/k_{ij}\) are the multicomponent reactivity ratios.

References

• Kazemi, N., Duever, T.A. and Penlidis, A. (2014), Demystifying the estimation of reactivity ratios for terpolymerization systems. AIChE J., 60: 1752-1766.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
f2	Molar fraction of M2. TYPE: float | FloatArrayLike
r12	Reactivity ratio. TYPE: float
r21	Reactivity ratio. TYPE: float
r13	Reactivity ratio. TYPE: float
r31	Reactivity ratio. TYPE: float
r23	Reactivity ratio. TYPE: float
r32	Reactivity ratio. TYPE: float

RETURNS	DESCRIPTION
tuple[float | FloatArray, ...]	Instantaneous terpolymer composition, \((F_1, F_2, F_3)\).

See also

• inst_copolymer_binary: specific method for binary systems.
• inst_copolymer_multi: generic method for multicomponent systems.

Examples:

>>> from polykin.copolymerization import inst_copolymer_ternary
>>> F1, F2, F3 = inst_copolymer_ternary(f1=0.5, f2=0.3, r12=0.2, r21=2.3,
...                                     r13=3.0, r31=0.9, r23=0.4, r32=1.5)
>>> print(f"F1 = {F1:.2f}; F2 = {F2:.2f}; F3 = {F3:.2f}")
F1 = 0.32; F2 = 0.41; F3 = 0.27

Source code in src/polykin/copolymerization/multicomponent.py

def inst_copolymer_ternary(f1: Union[float, FloatArrayLike],
                           f2: Union[float, FloatArrayLike],
                           r12: float,
                           r21: float,
                           r13: float,
                           r31: float,
                           r23: float,
                           r32: float,
                           ) -> tuple[Union[float, FloatArray],
                                      Union[float, FloatArray],
                                      Union[float, FloatArray]]:
    r"""Calculate the instantaneous copolymer composition for a ternary system.

    In a ternary system, the instantaneous copolymer composition $F_i$ is
    related to the monomer composition $f_i$ by:

    \begin{aligned}
        a &= \frac{f_1}{r_{21} r_{31}} + \frac{f_2}{r_{21} r_{32}} +
             \frac{f_3}{r_{31} r_{23}} \\
        b &= f_1 + \frac{f_2}{r_{12}} + \frac{f_3}{r_{13}} \\
        c &= \frac{f_1}{r_{12} r_{31}} + \frac{f_2}{r_{12} r_{32}} +
             \frac{f_3}{r_{13} r_{32}} \\
        d &= f_2 + \frac{f_1}{r_{21}} + \frac{f_3}{r_{23}} \\
        e &= \frac{f_1}{r_{13} r_{21}} + \frac{f_2}{r_{23} r_{12}} +
             \frac{f_3}{r_{13} r_{23}} \\
        g &= f_3 + \frac{f_1}{r_{31}} + \frac{f_2}{r_{32}} \\
        F_1 &= \frac{a b f_1}{a b f_1 + c d f_2 + e g f_3} \\
        F_2 &= \frac{c d f_2}{a b f_1 + c d f_2 + e g f_3} \\
        F_3 &= \frac{e g f_3}{a b f_1 + c d f_2 + e g f_3}
    \end{aligned}

    where $r_{ij}=k_{ii}/k_{ij}$ are the multicomponent reactivity ratios.

    **References**

    *   Kazemi, N., Duever, T.A. and Penlidis, A. (2014), Demystifying the
    estimation of reactivity ratios for terpolymerization systems. AIChE J.,
    60: 1752-1766.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    f2 : float | FloatArrayLike
        Molar fraction of M2.
    r12 : float
        Reactivity ratio.
    r21 : float
        Reactivity ratio.
    r13 : float
        Reactivity ratio.
    r31 : float
        Reactivity ratio.
    r23 : float
        Reactivity ratio.
    r32 : float
        Reactivity ratio.

    Returns
    -------
    tuple[float | FloatArray, ...]
        Instantaneous terpolymer composition, $(F_1, F_2, F_3)$.

    See also
    --------
    * [`inst_copolymer_binary`](inst_copolymer_binary.md):
      specific method for binary systems.
    * [`inst_copolymer_multi`](inst_copolymer_multi.md):
      generic method for multicomponent systems.

    Examples
    --------
    >>> from polykin.copolymerization import inst_copolymer_ternary
    >>> F1, F2, F3 = inst_copolymer_ternary(f1=0.5, f2=0.3, r12=0.2, r21=2.3,
    ...                                     r13=3.0, r31=0.9, r23=0.4, r32=1.5)
    >>> print(f"F1 = {F1:.2f}; F2 = {F2:.2f}; F3 = {F3:.2f}")
    F1 = 0.32; F2 = 0.41; F3 = 0.27
    """

    f1 = np.asarray(f1)
    f2 = np.asarray(f2)

    f3 = 1. - (f1 + f2)

    a = f1/(r21*r31) + f2/(r21*r32) + f3/(r31*r23)
    b = f1 + f2/r12 + f3/r13
    c = f1/(r12*r31) + f2/(r12*r32) + f3/(r13*r32)
    d = f2 + f1/r21 + f3/r23
    e = f1/(r13*r21) + f2/(r23*r12) + f3/(r13*r23)
    g = f3 + f1/r31 + f2/r32

    denominator = f1*a*b + f2*c*d + f3*e*g

    F1 = f1*a*b/denominator
    F2 = f2*c*d/denominator
    F3 = 1. - (F1 + F2)

    return (F1, F2, F3)

kpDataset

Dataset of average propagation rate coefficient data.

Container for average propagation rate coefficient as a function of monomer composition, \(k_p\) vs \(f_1\).

PARAMETER	DESCRIPTION
M1	Name of M1. TYPE: str
M2	Name of M2. TYPE: str
f1	Vector of monomer composition, \(f_1\). TYPE: FloatVectorLike
kp	Vector of average propagation rate coefficient, \(\bar{k}_p\). Unit = L/(mol·s). TYPE: FloatVectorLike
sigma_f1	Absolute standard deviation of \(f_1\). TYPE: Union[float, FloatVectorLike] DEFAULT: 0.05
sigma_kp	Absolute standard deviation of \(\bar{k}_p\). Unit = L/(mol·s). TYPE: Union[float, FloatVectorLike] DEFAULT: 100.0
weight	Relative weight of dataset for fitting. TYPE: float DEFAULT: 1
T	Temperature. Unit = Tunit. TYPE: float DEFAULT: 298.0
Tunit	Temperature unit. TYPE: Literal['C', 'K'] DEFAULT: 'K'
name	Name of dataset. TYPE: str DEFAULT: ''
source	Source of dataset. TYPE: str DEFAULT: ''

Source code in src/polykin/copolymerization/copodataset.py

class kpDataset(CopoDataset):
    r"""Dataset of average propagation rate coefficient data.

    Container for average propagation rate coefficient as a function of monomer
    composition, $k_p$ vs $f_1$.

    Parameters
    ----------
    M1 : str
        Name of M1.
    M2 : str
        Name of M2.
    f1 : FloatVectorLike
        Vector of monomer composition, $f_1$.
    kp : FloatVectorLike
        Vector of average propagation rate coefficient, $\bar{k}_p$.
        Unit = L/(mol·s).
    sigma_f1: float | FloatVectorLike
        Absolute standard deviation of $f_1$.
    sigma_kp: float | FloatVectorLike
        Absolute standard deviation of $\bar{k}_p$. Unit = L/(mol·s).
    weight: float
        Relative weight of dataset for fitting.
    T : float
        Temperature. Unit = `Tunit`.
    Tunit : Literal['C', 'K']
        Temperature unit.
    name: str
        Name of dataset.
    source: str
        Source of dataset.
    """

    varmap = {'f': 'x', 'kp': 'y'}

    def __init__(self,
                 M1: str,
                 M2: str,
                 f1: FloatVectorLike,
                 kp: FloatVectorLike,
                 sigma_f1: Union[float, FloatVectorLike] = 5e-2,
                 sigma_kp: Union[float, FloatVectorLike] = 1e2,
                 weight: float = 1,
                 T: float = 298.,
                 Tunit: Literal['C', 'K'] = 'K',
                 name: str = '',
                 source: str = '',
                 ) -> None:
        """Construct `DriftDataset` with the given parameters."""
        super().__init__(M1, M2, f1, kp, sigma_f1, sigma_kp, weight, T, Tunit,
                         name, source)

__init__

__init__(
    M1: str,
    M2: str,
    f1: FloatVectorLike,
    kp: FloatVectorLike,
    sigma_f1: Union[float, FloatVectorLike] = 0.05,
    sigma_kp: Union[float, FloatVectorLike] = 100.0,
    weight: float = 1,
    T: float = 298.0,
    Tunit: Literal["C", "K"] = "K",
    name: str = "",
    source: str = "",
) -> None

Construct DriftDataset with the given parameters.

Source code in src/polykin/copolymerization/copodataset.py

def __init__(self,
             M1: str,
             M2: str,
             f1: FloatVectorLike,
             kp: FloatVectorLike,
             sigma_f1: Union[float, FloatVectorLike] = 5e-2,
             sigma_kp: Union[float, FloatVectorLike] = 1e2,
             weight: float = 1,
             T: float = 298.,
             Tunit: Literal['C', 'K'] = 'K',
             name: str = '',
             source: str = '',
             ) -> None:
    """Construct `DriftDataset` with the given parameters."""
    super().__init__(M1, M2, f1, kp, sigma_f1, sigma_kp, weight, T, Tunit,
                     name, source)

kp_average_binary

kp_average_binary(
    f1: Union[float, FloatArrayLike],
    r1: Union[float, FloatArray],
    r2: Union[float, FloatArray],
    k11: Union[float, FloatArray],
    k22: Union[float, FloatArray],
) -> Union[float, FloatArray]

Calculate the average propagation rate coefficient in a copolymerization.

For a binary system, the instantaneous average propagation rate coefficient is related to the instantaneous comonomer composition by:

\[ \bar{k}_p = \frac{r_1 f_1^2 + r_2 f_2^2 + 2f_1 f_2} {(r_1 f_1/k_{11}) + (r_2 f_2/k_{22})} \]

where \(f_i\) is the instantaneous comonomer composition of monomer \(i\), \(r_i\) are the reactivity ratios, and \(k_{ii}\) are the homo-propagation rate coefficients. Although the equation is written using terminal model notation, it is equally applicable in the frame of the penultimate model if \(r_i \rightarrow \bar{r}_i\) and \(k_{ii} \rightarrow \bar{k}_{ii}\).

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
r1	Reactivity ratio of M1, \(r_1\) or \(\bar{r}_1\). TYPE: float | FloatArray
r2	Reactivity ratio of M2, \(r_2\) or \(\bar{r}_2\). TYPE: float | FloatArray
k11	Propagation rate coefficient of M1, \(k_{11}\) or \(\bar{k}_{11}\). Unit = L/(mol·s) TYPE: float | FloatArray
k22	Propagation rate coefficient of M2, \(k_{22}\) or \(\bar{k}_{22}\). Unit = L/(mol·s) TYPE: float | FloatArray

RETURNS	DESCRIPTION
float | FloatArray	Average propagation rate coefficient. Unit = L/(mol·s)

Examples:

>>> from polykin.copolymerization import kp_average_binary

An example with f1 as scalar.

>>> kp = kp_average_binary(f1=0.5, r1=0.16, r2=0.70, k11=100., k22=1000.)
>>> print(f"{kp:.0f} L/(mol·s)")
622 L/(mol·s)

An example with f1 as list.

>>> f1 = [0.2, 0.6, 0.8]
>>> kp = kp_average_binary(f1=f1, r1=0.16, r2=0.70, k11=100., k22=1000.)
>>> kp
array([880.        , 523.87096774, 317.18309859])

Source code in src/polykin/copolymerization/binary.py

def kp_average_binary(f1: Union[float, FloatArrayLike],
                      r1: Union[float, FloatArray],
                      r2: Union[float, FloatArray],
                      k11: Union[float, FloatArray],
                      k22: Union[float, FloatArray]
                      ) -> Union[float, FloatArray]:
    r"""Calculate the average propagation rate coefficient in a
    copolymerization.

    For a binary system, the instantaneous average propagation rate
    coefficient is related to the instantaneous comonomer composition by:

    $$ \bar{k}_p = \frac{r_1 f_1^2 + r_2 f_2^2 + 2f_1 f_2}
        {(r_1 f_1/k_{11}) + (r_2 f_2/k_{22})} $$

    where $f_i$ is the instantaneous comonomer composition of monomer $i$,
    $r_i$ are the reactivity ratios, and $k_{ii}$ are the homo-propagation
    rate coefficients. Although the equation is written using terminal
    model notation, it is equally applicable in the frame of the
    penultimate model if $r_i \rightarrow \bar{r}_i$ and
    $k_{ii} \rightarrow \bar{k}_{ii}$.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    r1 : float | FloatArray
        Reactivity ratio of M1, $r_1$ or $\bar{r}_1$.
    r2 : float | FloatArray
        Reactivity ratio of M2, $r_2$ or $\bar{r}_2$.
    k11 : float | FloatArray
        Propagation rate coefficient of M1, $k_{11}$ or $\bar{k}_{11}$.
        Unit = L/(mol·s)
    k22 : float | FloatArray
        Propagation rate coefficient of M2, $k_{22}$ or $\bar{k}_{22}$.
        Unit = L/(mol·s)

    Returns
    -------
    float | FloatArray
        Average propagation rate coefficient. Unit = L/(mol·s)

    Examples
    --------
    >>> from polykin.copolymerization import kp_average_binary

    An example with f1 as scalar.
    >>> kp = kp_average_binary(f1=0.5, r1=0.16, r2=0.70, k11=100., k22=1000.)
    >>> print(f"{kp:.0f} L/(mol·s)")
    622 L/(mol·s)

    An example with f1 as list.
    >>> f1 = [0.2, 0.6, 0.8]
    >>> kp = kp_average_binary(f1=f1, r1=0.16, r2=0.70, k11=100., k22=1000.)
    >>> kp
    array([880.        , 523.87096774, 317.18309859])
    """
    f1 = np.asarray(f1)
    f2 = 1 - f1
    return (r1*f1**2 + r2*f2**2 + 2*f1*f2)/((r1*f1/k11) + (r2*f2/k22))

monomer_drift_binary

monomer_drift_binary(
    f10: Union[float, FloatVectorLike],
    x: FloatVectorLike,
    r1: float,
    r2: float,
    atol: float = 0.0001,
    rtol: float = 0.0001,
    method: Literal["LSODA", "RK45"] = "LSODA",
) -> FloatMatrix

Compute the monomer composition drift for a binary system.

In a closed binary system, the drift in monomer composition is given by the solution of the following differential equation:

\[ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} \]

with initial condition \(f_1(0)=f_{1,0}\), where \(f_1\) and \(F_1\) are, respectively, the instantaneous comonomer and copolymer composition of M1, and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f10	Initial molar fraction of M1, \(f_{1,0}=f_1(0)\). TYPE: float | FloatVectorLike(N)
x	Total monomer conversion values where the drift is to be evaluated. TYPE: FloatVectorLike(M)
r1	Reactivity ratio of M1. TYPE: float | FloatArray
r2	Reactivity ratio of M2. TYPE: float | FloatArray
atol	Absolute tolerance of ODE solver. TYPE: float DEFAULT: 0.0001
rtol	Relative tolerance of ODE solver. TYPE: float DEFAULT: 0.0001
method	ODE solver. TYPE: Literal['LSODA', 'RK45'] DEFAULT: 'LSODA'

RETURNS	DESCRIPTION
FloatMatrix(M, N)	Monomer fraction of M1 at a given conversion, \(f_1(x)\).

See also

• monomer_drift_multi: generic method for multicomponent systems.

Examples:

>>> from polykin.copolymerization import monomer_drift_binary

An example with f10 as scalar.

>>> f1 = monomer_drift_binary(f10=0.5, x=[0.1, 0.5, 0.9], r1=0.16, r2=0.70)
>>> f1
array([0.51026241, 0.57810678, 0.87768138])

An example with f10 as list.

>>> f1 = monomer_drift_binary(f10=[0.2, 0.8], x=[0.1, 0.5, 0.9],
...                           r1=0.16, r2=0.70)
>>> f1
array([[0.19841009, 0.18898084, 0.15854395],
       [0.82315475, 0.94379024, 0.99996457]])

Source code in src/polykin/copolymerization/binary.py

def monomer_drift_binary(f10: Union[float, FloatVectorLike],
                         x: FloatVectorLike,
                         r1: float,
                         r2: float,
                         atol: float = 1e-4,
                         rtol: float = 1e-4,
                         method: Literal['LSODA', 'RK45'] = 'LSODA'
                         ) -> FloatMatrix:
    r"""Compute the monomer composition drift for a binary system.

    In a closed binary system, the drift in monomer composition is given by
    the solution of the following differential equation:

    $$ \frac{\textup{d} f_1}{\textup{d}x} = \frac{f_1 - F_1}{1 - x} $$

    with initial condition $f_1(0)=f_{1,0}$, where $f_1$ and $F_1$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    M1, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f10 : float | FloatVectorLike (N)
        Initial molar fraction of M1, $f_{1,0}=f_1(0)$.
    x : FloatVectorLike (M)
        Total monomer conversion values where the drift is to be evaluated.
    r1 : float | FloatArray
        Reactivity ratio of M1.
    r2 : float | FloatArray
        Reactivity ratio of M2.
    atol : float
        Absolute tolerance of ODE solver.
    rtol : float
        Relative tolerance of ODE solver.
    method : Literal['LSODA', 'RK45']
        ODE solver.

    Returns
    -------
    FloatMatrix (M, N)
        Monomer fraction of M1 at a given conversion, $f_1(x)$.

    See also
    --------
    * [`monomer_drift_multi`](monomer_drift_multi.md):
      generic method for multicomponent systems.

    Examples
    --------
    >>> from polykin.copolymerization import monomer_drift_binary

    An example with f10 as scalar.
    >>> f1 = monomer_drift_binary(f10=0.5, x=[0.1, 0.5, 0.9], r1=0.16, r2=0.70)
    >>> f1
    array([0.51026241, 0.57810678, 0.87768138])

    An example with f10 as list.
    >>> f1 = monomer_drift_binary(f10=[0.2, 0.8], x=[0.1, 0.5, 0.9],
    ...                           r1=0.16, r2=0.70)
    >>> f1
    array([[0.19841009, 0.18898084, 0.15854395],
           [0.82315475, 0.94379024, 0.99996457]])
    """

    if isinstance(f10, (int, float)):
        f10 = [f10]

    sol = solve_ivp(df1dx,
                    (0., max(x)),
                    f10,
                    args=(r1, r2),
                    t_eval=x,
                    method=method,
                    vectorized=True,
                    atol=atol,
                    rtol=rtol)

    if sol.success:
        result = sol.y
        result = np.maximum(0., result)
        if result.shape[0] == 1:
            result = result[0]
    else:
        raise ODESolverError(sol.message)

    return result

monomer_drift_multi

monomer_drift_multi(
    f0: FloatVectorLike,
    x: FloatVectorLike,
    r: FloatSquareMatrix,
    atol: float = 0.0001,
    rtol: float = 0.0001,
) -> FloatMatrix

Compute the monomer composition drift for a multicomponent system.

In a closed system, the drift in monomer composition is given by the solution of the following system of differential equations:

\[ \frac{\textup{d} f_i}{\textup{d}x} = \frac{f_i - F_i}{1 - x} \]

with initial condition \(f_i(0)=f_{i,0}\), where \(f_i\) and \(F_i\) are, respectively, the instantaneous comonomer and copolymer composition of monomer \(i\), and \(x\) is the total molar monomer conversion.

PARAMETER	DESCRIPTION
f0	Vector of initial instantaneous comonomer compositions. TYPE: FloatVectorLike(N)
x	Vector of total monomer conversion values where the drift is to be evaluated. TYPE: FloatVectorLike(M)
r	Matrix of reactivity ratios, \(r_{ij}=k_{ii}/k_{ij}\). TYPE: FloatSquareMatrix(N, N)
atol	Absolute tolerance of the ODE solver. TYPE: float DEFAULT: 0.0001
rtol	Relative tolerance of the ODE solver. TYPE: float DEFAULT: 0.0001

RETURNS	DESCRIPTION
FloatMatrix(M, N)	Matrix of monomer fraction of monomer \(i\) at the specified total monomer conversions, \(f_i(x_j)\).

See also

• inst_copolymer_multi: instantaneous copolymer composition.
• monomer_drift_multi: specific method for binary systems.

Examples:

>>> from polykin.copolymerization import monomer_drift_multi
>>> import numpy as np

Define reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate monomer drift.

>>> f0 = [0.5, 0.3, 0.2]
>>> x = [0.1, 0.5, 0.9, 0.99]
>>> f = monomer_drift_multi(f0, x, r)
>>> f
array([[5.19272893e-01, 2.87851432e-01, 1.92875675e-01],
       [6.38613228e-01, 2.04334321e-01, 1.57052451e-01],
       [8.31122266e-01, 5.58847454e-03, 1.63289259e-01],
       [4.98294381e-01, 1.22646553e-07, 5.01705497e-01]])

Source code in src/polykin/copolymerization/multicomponent.py

def monomer_drift_multi(f0: FloatVectorLike,
                        x: FloatVectorLike,
                        r: FloatSquareMatrix,
                        atol: float = 1e-4,
                        rtol: float = 1e-4
                        ) -> FloatMatrix:
    r"""Compute the monomer composition drift for a multicomponent system.

    In a closed system, the drift in monomer composition is given by
    the solution of the following system of differential equations:

    $$ \frac{\textup{d} f_i}{\textup{d}x} = \frac{f_i - F_i}{1 - x} $$

    with initial condition $f_i(0)=f_{i,0}$, where $f_i$ and $F_i$ are,
    respectively, the instantaneous comonomer and copolymer composition of
    monomer $i$, and $x$ is the total molar monomer conversion.

    Parameters
    ----------
    f0 : FloatVectorLike (N)
        Vector of initial instantaneous comonomer compositions.
    x : FloatVectorLike (M)
        Vector of total monomer conversion values where the drift is to be
        evaluated.
    r : FloatSquareMatrix (N, N)
        Matrix of reactivity ratios, $r_{ij}=k_{ii}/k_{ij}$.
    atol : float
        Absolute tolerance of the ODE solver.
    rtol : float
        Relative tolerance of the ODE solver.

    Returns
    -------
    FloatMatrix (M, N)
        Matrix of monomer fraction of monomer $i$ at the specified total
        monomer conversions, $f_i(x_j)$.

    See also
    --------
    * [`inst_copolymer_multi`](inst_copolymer_multi.md): 
      instantaneous copolymer composition.
    * [`monomer_drift_multi`](monomer_drift_multi.md):
      specific method for binary systems.

    Examples
    --------
    >>> from polykin.copolymerization import monomer_drift_multi
    >>> import numpy as np

    Define reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate monomer drift.
    >>> f0 = [0.5, 0.3, 0.2]
    >>> x = [0.1, 0.5, 0.9, 0.99]
    >>> f = monomer_drift_multi(f0, x, r)
    >>> f
    array([[5.19272893e-01, 2.87851432e-01, 1.92875675e-01],
           [6.38613228e-01, 2.04334321e-01, 1.57052451e-01],
           [8.31122266e-01, 5.58847454e-03, 1.63289259e-01],
           [4.98294381e-01, 1.22646553e-07, 5.01705497e-01]])

    """

    N = len(f0)
    if r.shape != (N, N):
        raise ShapeError("Shape mismatch between `f0` and `r`.")

    def dfdx(x: float, f: FloatVector) -> FloatVector:
        F = inst_copolymer_multi(f=np.append(f, 1. - f.sum()), r=r)
        return (f - F[:-1]) / (1. - x + eps)

    sol = solve_ivp(dfdx,
                    t_span=(0., max(x)),
                    y0=f0[:-1],
                    t_eval=x,
                    atol=atol,
                    rtol=rtol,
                    method='LSODA')

    if sol.success:
        f = np.empty((len(x), N))
        f[:, :-1] = np.transpose(sol.y)
        f[:, -1] = 1. - np.sum(f[:, :-1], axis=1)
    else:
        raise ODESolverError(sol.message)

    return f

radical_fractions_multi

radical_fractions_multi(
    f: FloatVectorLike, k: FloatSquareMatrix
) -> FloatVector

Calculate the radical fractions for a multicomponent system.

In a multicomponent system, the radical fractions \(p_i\) can be determined by solving the following set of linear algebraic equations:

\[ \sum_{j\ne i}^N k_{ij} p_i f_j = \sum_{j\ne i}^N k_{ji} p_j f_i \]

where \(k_{ij}\) are the cross-propagation rate coefficients and \(f_i\) are the monomer compositions. Note that the homo-propagation rate coefficients \(k_{ii}\) do not appear in the equations. For this reason, radical fractions cannot be evaluated from reactivity ratios alone.

PARAMETER	DESCRIPTION
f	Vector of instantaneous monomer compositions, \(f_i\). TYPE: FloatVectorLike(N)
k	Matrix of cross-propagation rate coefficients. The diagonal elements \(k_{ii}\) are not used. TYPE: FloatSquareMatrix(N, N)

RETURNS	DESCRIPTION
FloatVector(N)	Vector of radical fractions, \(p_i\).

See also

• radical_fractions_ternary: specific method for terpolymer systems.

Examples:

>>> from polykin.copolymerization import radical_fractions_multi
>>> import numpy as np

Define the cross-propagation coefficient matrix.

>>> k = np.zeros((3, 3))
>>> k[0, 1] = 500.
>>> k[1, 0] = 50.
>>> k[0, 2] = 30.
>>> k[2, 0] = 200.
>>> k[1, 2] = 300.
>>> k[2, 1] = 40.

Evaluate the radical fractions at f1=0.5, f2=0.3, f3=0.2.

>>> f = [0.5, 0.3, 0.2]
>>> p = radical_fractions_multi(f, k)
>>> p
array([0.25012791, 0.47956341, 0.27030868])

Source code in src/polykin/copolymerization/multicomponent.py

def radical_fractions_multi(f: FloatVectorLike,
                            k: FloatSquareMatrix,
                            ) -> FloatVector:
    r"""Calculate the radical fractions for a multicomponent system.

    In a multicomponent system, the radical fractions $p_i$ can be determined
    by solving the following set of linear algebraic equations:

    $$ \sum_{j\ne i}^N k_{ij} p_i f_j = \sum_{j\ne i}^N k_{ji} p_j f_i $$

    where $k_{ij}$ are the cross-propagation rate coefficients and $f_i$ are
    the monomer compositions. Note that the homo-propagation rate coefficients
    $k_{ii}$ do not appear in the equations. For this reason, radical fractions
    cannot be evaluated from reactivity ratios alone.

    Parameters
    ----------
    f : FloatVectorLike (N)
        Vector of instantaneous monomer compositions, $f_i$.
    k : FloatSquareMatrix (N, N)
        Matrix of cross-propagation rate coefficients. The diagonal elements
        $k_{ii}$ are not used.

    Returns
    -------
    FloatVector (N)
        Vector of radical fractions, $p_i$.

    See also
    --------
    * [`radical_fractions_ternary`](radical_fractions_ternary.md):
      specific method for terpolymer systems.

    Examples
    --------
    >>> from polykin.copolymerization import radical_fractions_multi
    >>> import numpy as np

    Define the cross-propagation coefficient matrix.
    >>> k = np.zeros((3, 3))
    >>> k[0, 1] = 500.
    >>> k[1, 0] = 50.
    >>> k[0, 2] = 30.
    >>> k[2, 0] = 200.
    >>> k[1, 2] = 300.
    >>> k[2, 1] = 40.

    Evaluate the radical fractions at f1=0.5, f2=0.3, f3=0.2.
    >>> f = [0.5, 0.3, 0.2]
    >>> p = radical_fractions_multi(f, k)
    >>> p
    array([0.25012791, 0.47956341, 0.27030868])
    """
    f = np.asarray(f)

    A = k.T * f[:, np.newaxis]
    x = A.sum(axis=0)
    np.fill_diagonal(A, A.diagonal() - x)
    A[-1, :] = 1.
    b = np.zeros(A.shape[0])
    b[-1] = 1.
    p = np.linalg.solve(A, b)

    return p

radical_fractions_ternary

radical_fractions_ternary(
    f1: Union[float, FloatArrayLike],
    f2: Union[float, FloatArrayLike],
    k12: float,
    k21: float,
    k13: float,
    k31: float,
    k23: float,
    k32: float,
) -> tuple[
    Union[float, FloatArray],
    Union[float, FloatArray],
    Union[float, FloatArray],
]

Calculate the radical fractions for a ternary system.

In a ternary system, the radical fractions \(p_i\) are related to the monomer composition \(f_i\) by:

\[\begin{aligned} a &= k_{21}k_{31}f_1^2+k_{21}k_{32}f_1f_2+k_{23}k_{31}f_1f_3 \\ b &= k_{12}k_{31}f_1f_2+k_{12}k_{32}f_2^2+k_{13}k_{32}f_2f_3 \\ c &= k_{12}k_{23}f_2f_3+k_{13}k_{21}f_1f_3+k_{13}k_{23}f_3^2 \\ p_1 &= \frac{a}{a+b+c} \\ p_2 &= \frac{b}{a+b+c} \\ p_3 &= \frac{c}{a+b+c} \end{aligned}\]

where \(k_{ij}\) are the cross-propagation rate coefficients. Note that the homo-propagation rate coefficients \(k_{ii}\) do not appear in the equations. For this reason, radical fractions cannot be evaluated from reactivity ratios alone.

References

• Hamielec, A.E., MacGregor, J.F. and Penlidis, A. (1987), Multicomponent free-radical polymerization in batch, semi- batch and continuous reactors. Makromolekulare Chemie. Macromolecular Symposia, 10-11: 521-570.

PARAMETER	DESCRIPTION
f1	Molar fraction of M1. TYPE: float | FloatArrayLike
f2	Molar fraction of M2. TYPE: float | FloatArrayLike
k12	Propagation rate coefficient. TYPE: float
k21	Propagation rate coefficient. TYPE: float
k13	Propagation rate coefficient. TYPE: float
k31	Propagation rate coefficient. TYPE: float
k23	Propagation rate coefficient. TYPE: float
k32	Propagation rate coefficient. TYPE: float

RETURNS	DESCRIPTION
tuple[float | FloatArray, ...]	Radical fractions, \((p_1, p_2, p_3)\).

See also

• radical_fractions_multi: generic method for multicomponent systems.

Examples:

>>> from polykin.copolymerization import radical_fractions_ternary
>>> p1, p2, p3 = radical_fractions_ternary(
...              f1=0.5, f2=0.3, k12=500., k21=50.,
...              k13=30., k31=200., k23=300., k32=40.)
>>> print(f"p1 = {p1:.2f}; p2 = {p2:.2f}; p3 = {p3:.2f}")
p1 = 0.25; p2 = 0.48; p3 = 0.27

Source code in src/polykin/copolymerization/multicomponent.py

def radical_fractions_ternary(f1: Union[float, FloatArrayLike],
                              f2: Union[float, FloatArrayLike],
                              k12: float,
                              k21: float,
                              k13: float,
                              k31: float,
                              k23: float,
                              k32: float,
                              ) -> tuple[Union[float, FloatArray],
                                         Union[float, FloatArray],
                                         Union[float, FloatArray]]:
    r"""Calculate the radical fractions for a ternary system.

    In a ternary system, the radical fractions $p_i$ are related to the
    monomer composition $f_i$ by:

    \begin{aligned}
    a &= k_{21}k_{31}f_1^2+k_{21}k_{32}f_1f_2+k_{23}k_{31}f_1f_3 \\
    b &= k_{12}k_{31}f_1f_2+k_{12}k_{32}f_2^2+k_{13}k_{32}f_2f_3 \\
    c &= k_{12}k_{23}f_2f_3+k_{13}k_{21}f_1f_3+k_{13}k_{23}f_3^2 \\
    p_1 &= \frac{a}{a+b+c} \\
    p_2 &= \frac{b}{a+b+c} \\
    p_3 &= \frac{c}{a+b+c}
    \end{aligned}

    where $k_{ij}$ are the cross-propagation rate coefficients. Note that the
    homo-propagation rate coefficients $k_{ii}$ do not appear in the equations.
    For this reason, radical fractions cannot be evaluated from reactivity
    ratios alone.

    **References**

    *   Hamielec, A.E., MacGregor, J.F. and Penlidis, A. (1987), Multicomponent
    free-radical polymerization in batch, semi- batch and continuous reactors.
    Makromolekulare Chemie. Macromolecular Symposia, 10-11: 521-570.

    Parameters
    ----------
    f1 : float | FloatArrayLike
        Molar fraction of M1.
    f2 : float | FloatArrayLike
        Molar fraction of M2.
    k12 : float
        Propagation rate coefficient.
    k21 : float
        Propagation rate coefficient.
    k13 : float
        Propagation rate coefficient.
    k31 : float
        Propagation rate coefficient.
    k23 : float
        Propagation rate coefficient.
    k32 : float
        Propagation rate coefficient.

    Returns
    -------
    tuple[float | FloatArray, ...]
        Radical fractions, $(p_1, p_2, p_3)$.

    See also
    --------
    * [`radical_fractions_multi`](radical_fractions_multi.md):
      generic method for multicomponent systems.

    Examples
    --------
    >>> from polykin.copolymerization import radical_fractions_ternary
    >>> p1, p2, p3 = radical_fractions_ternary(
    ...              f1=0.5, f2=0.3, k12=500., k21=50.,
    ...              k13=30., k31=200., k23=300., k32=40.)
    >>> print(f"p1 = {p1:.2f}; p2 = {p2:.2f}; p3 = {p3:.2f}")
    p1 = 0.25; p2 = 0.48; p3 = 0.27
    """

    f1 = np.asarray(f1)
    f2 = np.asarray(f2)
    f3 = 1. - (f1 + f2)

    p1 = k21*k31*f1**2 + k21*k32*f1*f2 + k23*k31*f1*f3
    p2 = k12*k31*f1*f2 + k12*k32*f2**2 + k13*k32*f2*f3
    p3 = k12*k23*f2*f3 + k13*k21*f1*f3 + k13*k23*f3**2

    denominator = p1 + p2 + p3

    p1 /= denominator
    p2 /= denominator
    p3 /= denominator

    return (p1, p2, p3)

sequence_multi

sequence_multi(
    Pself: FloatVectorLike,
    k: Optional[Union[int, IntArrayLike]] = None,
) -> FloatArray

Calculate the instantaneous sequence length probability or the number-average sequence length.

For a multicomponent system, the probability of finding \(k\) consecutive units of monomer \(i\) in a chain is:

\[ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} \]

and the corresponding number-average sequence length is:

\[ \bar{S}_i = \sum_{k=1}^{\infty} k S_{i,k} = \frac{1}{1 - P_{ii}} \]

where \(P_{ii}\) is the self-transition probability \(i \rightarrow i\), which is a function of the monomer composition and the reactivity ratios.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 177.

PARAMETER	DESCRIPTION
Pself	Vector of self-transition probabilities, \(P_{ii}\), corresponding to the diagonal of the matrix of transition probabilities. TYPE: FloatVectorLike(N)
k	Sequence length, i.e., number of consecutive units in a chain. If None, the number-average sequence length will be computed. TYPE: int | IntArrayLike(M) | None DEFAULT: None

RETURNS	DESCRIPTION
FloatArray(N, M)	If k is None, the number-average sequence lengths, \(\bar{S}_i\). Otherwise, the sequence probabilities, \(S_{i,k}\).

See also

• transitions_multi: instantaneous transition probabilities.
• tuples_multi: instantaneous tuple fractions.

Examples:

>>> from polykin.copolymerization import sequence_multi
>>> from polykin.copolymerization import transitions_multi
>>> import numpy as np

Define reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate self-transition probabilities.

>>> f = [0.5, 0.3, 0.2]
>>> Pself = transitions_multi(f, r).diagonal()
>>> Pself
array([0.24193548, 0.29487179, 0.20930233])

Evaluate number-average sequence lengths for all monomers.

>>> S = sequence_multi(Pself)
>>> S
array([1.31914894, 1.41818182, 1.26470588])

Evaluate probabilities for certain sequence lengths.

>>> S = sequence_multi(Pself, k=[1, 5])
>>> S
array([[0.75806452, 0.00259719],
       [0.70512821, 0.00533091],
       [0.79069767, 0.00151742]])

Source code in src/polykin/copolymerization/multicomponent.py

def sequence_multi(Pself: FloatVectorLike,
                   k: Optional[Union[int, IntArrayLike]] = None,
                   ) -> FloatArray:
    r"""Calculate the instantaneous sequence length probability or the
    number-average sequence length.

    For a multicomponent system, the probability of finding $k$ consecutive
    units of monomer $i$ in a chain is:

    $$ S_{i,k} = (1 - P_{ii})P_{ii}^{k-1} $$

    and the corresponding number-average sequence length is:

    $$ \bar{S}_i = \sum_{k=1}^{\infty} k S_{i,k} = \frac{1}{1 - P_{ii}} $$

    where $P_{ii}$ is the self-transition probability $i \rightarrow i$, which
    is a function of the monomer composition and the reactivity ratios.

    **References**

    * NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
    process modeling, Wiley, 1996, p. 177.

    Parameters
    ----------
    Pself : FloatVectorLike (N)
        Vector of self-transition probabilities, $P_{ii}$, corresponding to
        the diagonal of the matrix of transition probabilities.
    k : int | IntArrayLike (M) | None
        Sequence length, i.e., number of consecutive units in a chain.
        If `None`, the number-average sequence length will be computed.

    Returns
    -------
    FloatArray (N, M)
        If `k is None`, the number-average sequence lengths, $\bar{S}_i$.
        Otherwise, the sequence probabilities, $S_{i,k}$.

    See also
    --------
    * [`transitions_multi`](transitions_multi.md):
      instantaneous transition probabilities.
    * [`tuples_multi`](tuples_multi.md):
      instantaneous tuple fractions.

    Examples
    --------
    >>> from polykin.copolymerization import sequence_multi
    >>> from polykin.copolymerization import transitions_multi
    >>> import numpy as np

    Define reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate self-transition probabilities.
    >>> f = [0.5, 0.3, 0.2]
    >>> Pself = transitions_multi(f, r).diagonal()
    >>> Pself
    array([0.24193548, 0.29487179, 0.20930233])

    Evaluate number-average sequence lengths for all monomers.
    >>> S = sequence_multi(Pself)
    >>> S
    array([1.31914894, 1.41818182, 1.26470588])

    Evaluate probabilities for certain sequence lengths.
    >>> S = sequence_multi(Pself, k=[1, 5])
    >>> S
    array([[0.75806452, 0.00259719],
           [0.70512821, 0.00533091],
           [0.79069767, 0.00151742]])

    """

    Pself = np.asarray(Pself)

    if k is None:
        S = 1/(1. - Pself + eps)
    else:
        if isinstance(k, (list, tuple)):
            k = np.array(k, dtype=np.int32)
        if isinstance(k, np.ndarray):
            Pself = Pself.reshape(-1, 1)
        S = (1. - Pself)*Pself**(k - 1)

    return S

transitions_multi

transitions_multi(
    f: FloatVectorLike, r: FloatSquareMatrix
) -> FloatSquareMatrix

Calculate the instantaneous transition probabilities for a multicomponent system.

For a multicomponent system, the transition probabilities are given by:

\[ P_{ij} = \frac{r_{ij}^{-1} f_j} {\displaystyle \sum_{k=1}^{N} r_{ik}^{-1} f_k} \]

where \(f_i\) is the molar fraction of monomer \(i\) and \(r_{ij}=k_{ii}/k_{ij}\) is the multicomponent reactivity ratio matrix.

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 178.

PARAMETER	DESCRIPTION
f	Vector of instantaneous monomer compositions, \(f_i\). TYPE: FloatVectorLike(N)
r	Matrix of reactivity ratios, \(r_{ij}=k_{ii}/k_{ij}\). TYPE: FloatSquareMatrix(N, N)

RETURNS	DESCRIPTION
FloatSquareMatrix(N, N)	Matrix of transition probabilities, \(P_{ij}\).

See also

• inst_copolymer_multi: instantaneous copolymer composition.
• sequence_multi: instantaneous sequence lengths.
• tuples_multi: instantaneous tuple fractions.

Examples:

>>> from polykin.copolymerization import transitions_multi
>>> import numpy as np

Define reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate transition probabilities.

>>> f = [0.5, 0.3, 0.2]
>>> P = transitions_multi(f, r)
>>> P
array([[0.24193548, 0.72580645, 0.03225806],
       [0.21367521, 0.29487179, 0.49145299],
       [0.58139535, 0.20930233, 0.20930233]])

Source code in src/polykin/copolymerization/multicomponent.py

def transitions_multi(f: FloatVectorLike,
                      r: FloatSquareMatrix
                      ) -> FloatSquareMatrix:
    r"""Calculate the instantaneous transition probabilities for a
    multicomponent system.

    For a multicomponent system, the transition probabilities are given
    by:

    $$ P_{ij} = \frac{r_{ij}^{-1} f_j}
                     {\displaystyle \sum_{k=1}^{N} r_{ik}^{-1} f_k} $$

    where $f_i$ is the molar fraction of monomer $i$ and $r_{ij}=k_{ii}/k_{ij}$
    is the multicomponent reactivity ratio matrix.

    **References**

    *   NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 178.

    Parameters
    ----------
    f : FloatVectorLike (N)
        Vector of instantaneous monomer compositions, $f_i$.
    r : FloatSquareMatrix (N, N)
        Matrix of reactivity ratios, $r_{ij}=k_{ii}/k_{ij}$.

    Returns
    -------
    FloatSquareMatrix (N, N)
        Matrix of transition probabilities, $P_{ij}$.

    See also
    --------
    * [`inst_copolymer_multi`](inst_copolymer_multi.md):
      instantaneous copolymer composition.
    * [`sequence_multi`](sequence_multi.md):
      instantaneous sequence lengths.
    * [`tuples_multi`](tuples_multi.md):
      instantaneous tuple fractions.

    Examples
    --------
    >>> from polykin.copolymerization import transitions_multi
    >>> import numpy as np

    Define reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate transition probabilities.
    >>> f = [0.5, 0.3, 0.2]
    >>> P = transitions_multi(f, r)
    >>> P
    array([[0.24193548, 0.72580645, 0.03225806],
           [0.21367521, 0.29487179, 0.49145299],
           [0.58139535, 0.20930233, 0.20930233]])

    """
    f = np.asarray(f)

    # N = len(f)
    # P = np.empty((N, N))
    # for i in range(N):
    #     for j in range(N):
    #         P[i, j] = f[j]/r[i, j] / np.sum(f/r[i, :])
    P = (f/r) / np.sum(f/r, axis=1)[:, np.newaxis]

    return P

tuples_multi

tuples_multi(
    P: FloatSquareMatrix,
    n: int,
    F: Optional[FloatVectorLike] = None,
) -> dict[tuple[int, ...], float]

Calculate the instantaneous n-tuple fractions.

For a multicomponent system, the probability of finding a specific sequence \(ijk \cdots rs\) of repeating units is:

\[ A_{ijk \cdots rs} = F_i P_{ij} P_{jk} \cdots P_{rs} \]

where \(F_i\) is the instantaneous copolymer composition, and \(P_{ij}\) is the transition probability \(i \rightarrow j\). Since the direction of chain growth does not play a role, symmetric sequences are combined under the sequence with lower index (e.g., \(A_{112} \leftarrow A_{112} + A_{211}\)).

References

• NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization process modeling, Wiley, 1996, p. 179.

PARAMETER	DESCRIPTION
P	Matrix of transition probabilities, \(P_{ij}\). TYPE: FloatSquareMatrix(N, N)
n	Tuple length, e.g. monads (1), diads (2), triads (3), etc. TYPE: int
F	Vector of instantaneous copolymer composition, \(F_i\). If None, the value will be computed internally. When calculating tuples of various lengths, it is more efficient to precompute \(F\) and use it in all tuple cases. TYPE: FloatVectorLike(N) | None DEFAULT: None

RETURNS	DESCRIPTION
dict[tuple[int, ...], float]	Tuple of molar fractions.

See also

• sequence_multi: instantaneous sequence lengths.
• transitions_multi: instantaneous transition probabilities.

Examples:

>>> from polykin.copolymerization import tuples_multi
>>> from polykin.copolymerization import transitions_multi
>>> import numpy as np

Define reactivity ratio matrix.

>>> r = np.ones((3, 3))
>>> r[0, 1] = 0.2
>>> r[1, 0] = 2.3
>>> r[0, 2] = 3.0
>>> r[2, 0] = 0.9
>>> r[1, 2] = 0.4
>>> r[2, 1] = 1.5

Evaluate transition probabilities.

>>> f = [0.5, 0.3, 0.2]
>>> P = transitions_multi(f, r)

Evaluate triad fractions.

>>> A = tuples_multi(P, 3)
>>> A[(0, 0, 0)]
0.018811329044450834
>>> A[(1, 0, 1)]
0.06365013630778116

Source code in src/polykin/copolymerization/multicomponent.py

def tuples_multi(P: FloatSquareMatrix,
                 n: int,
                 F: Optional[FloatVectorLike] = None
                 ) -> dict[tuple[int, ...], float]:
    r"""Calculate the instantaneous n-tuple fractions.

    For a multicomponent system, the probability of finding a specific sequence
    $ijk \cdots rs$ of repeating units is:

    $$ A_{ijk \cdots rs} = F_i P_{ij} P_{jk} \cdots P_{rs} $$

    where $F_i$ is the instantaneous copolymer composition, and $P_{ij}$ is
    the transition probability $i \rightarrow j$. Since the direction of chain
    growth does not play a role, symmetric sequences are combined under the
    sequence with lower index (e.g., $A_{112} \leftarrow A_{112} + A_{211}$).

    **References**

    *   NA Dotson, R Galván, RL Laurence, and M Tirrel. Polymerization
        process modeling, Wiley, 1996, p. 179.

    Parameters
    ----------
    P : FloatSquareMatrix (N, N)
        Matrix of transition probabilities, $P_{ij}$.
    n : int
        Tuple length, e.g. monads (1), diads (2), triads (3), etc.
    F : FloatVectorLike (N) | None
        Vector of instantaneous copolymer composition, $F_i$. If `None`,
        the value will be computed internally. When calculating tuples of
        various lengths, it is more efficient to precompute $F$ and use it in
        all tuple cases.

    Returns
    -------
    dict[tuple[int, ...], float]
        Tuple of molar fractions.

    See also
    --------
    * [`sequence_multi`](sequence_multi.md):
      instantaneous sequence lengths.
    * [`transitions_multi`](transitions_multi.md):
      instantaneous transition probabilities.

    Examples
    --------
    >>> from polykin.copolymerization import tuples_multi
    >>> from polykin.copolymerization import transitions_multi
    >>> import numpy as np

    Define reactivity ratio matrix.
    >>> r = np.ones((3, 3))
    >>> r[0, 1] = 0.2
    >>> r[1, 0] = 2.3
    >>> r[0, 2] = 3.0
    >>> r[2, 0] = 0.9
    >>> r[1, 2] = 0.4
    >>> r[2, 1] = 1.5

    Evaluate transition probabilities.
    >>> f = [0.5, 0.3, 0.2]
    >>> P = transitions_multi(f, r)

    Evaluate triad fractions.
    >>> A = tuples_multi(P, 3)
    >>> A[(0, 0, 0)]
    0.018811329044450834
    >>> A[(1, 0, 1)]
    0.06365013630778116

    """

    # Compute F if not given
    if F is None:
        F = inst_copolymer_multi(f=None, r=None, P=P)

    # Generate all tuple combinations
    N = P.shape[0]
    indexes = list(itertools.product(range(N), repeat=n))

    result = {}
    for idx in indexes:
        # Compute tuple probability
        P_product = 1.
        for j in range(n-1):
            P_product *= P[idx[j], idx[j+1]]
        A = F[idx[0]]*P_product
        # Add probability to dict, but combine symmetric tuples: 12 <- 12 + 21
        reversed_idx = idx[::-1]
        if reversed_idx in result.keys():
            result[reversed_idx] += A
        else:
            result[idx] = A

    return result

Tutorial