El "original"

El gráfico original con el que vamos a empezar es el siguiente.

Este no es el gráfico original en el artículo de 2015 pero sirve el propósito del reto. Los datos se pueden descargar aquí. Este grafico presenta la fracción de energía transmitida (\(\eta\)) para un compuesto helicoidadl con diferentes parámetros geométricos. Los parámetros que se variaron fueron los siguientes:

• Ángulo de rotación \(\alpha\): el ángulo formado entre capas consecutivas;

• Grosor del compuesto \(D\), normalizdo por la longitud de onda \(\lambda\); y

• Número de capas \(N\) en cada celda del compuesto.

El siguiente esquema ilustra estos parámetros.

Yo no diría que el gráfico es horrible, y, en comparación con lo que se encuentra en la literatura cientítica es incluso bueno. Pero … en tierra de ciegos, el tuerto es rey. Entonces, vamos a enumarara cada uno de los pecados del gráfico, y a corregirlos:

• Tiene dos ejes horizontales.

• La leyenda está encerrada en una recuadro que no es necesario.

• Los ejes de la derecha y la parte superior no contribuyen al gráfico.

• Las anotaciones están robando protagonismo a los datos.

• Se ve desordenado con las líneas y los marcadores.

• Es un gráfico espagueti.

El siguiente código genera este gráfico.

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rcParams

rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16
rcParams['legend.fontsize'] = 15
rcParams['mathtext.fontset'] = 'cm'

markers = ['o', '^', 'v', 's', '<', '>', 'd', '*', 'x', 'D', '+', 'H']
data = np.loadtxt("Energy_vs_D.csv", skiprows=1, delimiter=",")
labels = np.loadtxt("Energy_vs_D.csv", skiprows=0, delimiter=",",
                    usecols=range(1, 9))
labels = labels[0, :]

fig = plt.figure()
ax = plt.subplot(111)
for cont in range(8):
    plt.plot(data[:, 0], data[:, cont + 1], marker=markers[cont],
             label=r"$D/\lambda={:.3g}$".format(labels[cont]))

# First x-axis
xticks, xlabels = plt.xticks()
plt.xlabel(r"Number of layers - $N$", size=15)
plt.ylabel(r"Fraction of Energy - $\eta$", size=15)
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))

# Second x-axis
ax2 = ax.twiny()
ax2.set_xticks(xticks[2:])
ax2.set_xticklabels(180./xticks[2:])
plt.xlabel(r"Angle - $\alpha\ (^\circ)$", size=15)

plt.tight_layout()
plt.savefig("energy_vs_D_orig.svg")
plt.savefig("energy_vs_D_orig.png", dpi=300)

Correcciones

Voy a reivindicar el gráfico un pecado a la vez, veamos cómo resulta.

Tiene dos ejes horizontales

Originalmente, agregue dos ejes horizontales para mostrar el número de capas y el ángulo de rotación al mismo tiempo. La recomendación general es evitar esto, así que vamos a deshacernos del eje superior.

Leyenda en un recuadro

Es claro a qué se refiere …

Ejes a la derecha y arriba

Quítemoslos

Las anotaciones están robando protagonismo

Vamos a desplazar las anotaciones hacia el fondo cambiando el color del texto a un gris claro.

Líneas y marcadores

Conservemos únicamente las líneas.

E incrementemos su grosor.

Es un gráfico espagueti

Creo que una buena opción para este gráfico sería resaltar una línea al tiempo. Haciendo esto, terminamos con el siguiente gráfico.

El siguiente bloque de código genera esta versión.

import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rcParams

# Plots setup
gray = '#757575'
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["text.color"] = gray
plt.rcParams["xtick.color"] = gray
plt.rcParams["ytick.color"] = gray
plt.rcParams["axes.labelcolor"] = gray
plt.rcParams["axes.edgecolor"] = gray
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False
rcParams['font.family'] = 'serif'
rcParams['mathtext.fontset'] = 'cm'



def line_plots(data, highlight_line, title):
    plt.title(title)
    for cont, datum in enumerate(data[:, 1:].T):
        if cont == highlight_line:
            plt.plot(data[:, 0], datum, zorder=3, color="#984ea3",
                     linewidth=2)
        else:
            plt.plot(data[:, 0], datum, zorder=2, color=gray,
                     linewidth=1, alpha=0.5)


data = np.loadtxt("Energy_vs_D.csv", skiprows=1, delimiter=",")
labels = np.loadtxt("Energy_vs_D.csv", skiprows=0, delimiter=",",
                    usecols=range(1, 9))
labels = labels[0, :]

plt.close("all")
plt.figure(figsize=(8, 4))
for cont in range(8):
    ax = plt.subplot(2, 4, cont + 1)
    title = r"$D/\lambda={:.3g}$".format(labels[cont])
    line_plots(data, cont, title)
    plt.ylim(0.4, 1.0)
    if cont < 4:
        plt.xlabel("")
        ax.xaxis.set_ticks([])
        ax.spines["bottom"].set_color("none")
    else:
        plt.xlabel(r"Number of layers - $N$")
    if cont % 4 > 0:
        ax.yaxis.set_ticks([])
        ax.spines["left"].set_color("none")
    else:
        plt.ylabel(r"Fraction of Energy - $\eta$")


plt.tight_layout()
plt.savefig("energy_vs_D_7.svg")
plt.savefig("energy_vs_D_7.png", dpi=300)

Ajustes finales

Añadimos algunos detalles en Inkscape para terminar con la siguiente versión.

Lecturas adicionales

• Knaflic, Cole Nussbaumer. Storytelling with data: A data visualization guide for business professionals. John Wiley & Sons, 2015.

• Nicolás Guarín-Zapata et al. "Shear wave filtering in naturally-occurring Bouligand structures." Acta biomaterialia 23 (2015): 11-20. Preprint: https://arxiv.org/pdf/1505.04203.pdf

Primer enfoque

Repetiré el resumen de la semana pasada. Hay que tener en cuenta que Inkscape trata la rotación en sentido horario como positiva.

Luego, para crear una caja con las dimensiones deseadas, primero creamos cada rectángulo con las dimensiones correctas (en proyecciones paralelas). En el siguiente ejemplo, usamos caras con relaciones de aspecto 3:2, 2:1 y 4:3. La caja de la derecha es la cifra obtenida después de aplicar las transformaciones descritas en el esquema anterior.

Ahora podemos proceder, para recrear la figura deseada. No conozco las dimensiones exactas de la caja; supongo que la sección transversal era un cuadrado y la relación de aspecto era 5:1. Después de aplicar las transformaciones para cada rectángulo obtenemos lo siguiente (agregando algunos ajustes aquí y allá).

Segundo enfoque

Para este tipo de esquema, preferiría crear una cuadrícula axonométrica (Archivo > Propiedades de documento > Cuadrículas). Después de agregar la cuadrícula a nuestro lienzo es bastante sencillo dibujar las figuras en vista isométrica. El lienzo debe ser similar a la siguiente imagen.

Luego podemos dibujar cada cara usando la cuadrícula. Si queremos ser más precisos podemos activar Ajustar a nodos cúspides. La siguiente animación muestra la construcción paso a paso.

Y obtenemos la siguiente imagen.

Conclusión

Como mencioné, Inkscape se puede usar para dibujar figuras simples en proyección isométrica. Sin embargo, sugiero utilizar un CAD como FreeCAD para geometrías más complicadas.

Transformaciones involucradas

Como estamos trabajando en una pantalla de computador, estamos hablando de 2 dimensiones. Por lo tanto, todas las transformaciones están representadas por matrices 2×2. En el ejemplo de interés en este post necesitamos las siguientes transformaciones:

1. Escalado vertical

2. Inclinación horizontal

3. Rotación

Las siguientes son las matrices de transformación.

TL;DR

Si solo está interesado en las transformaciones necesarias en Inkscape Puedes consultar el resumen a continuación. Utiliza el tercer cuadrante como se presenta abajo.

Resumen

Tenga en cuenta que Inkscape trata la rotación en sentido horario como positiva. Opuesto a la notación común en matemáticas.

Diferencias finitas en dominios infinitos

Gracias a mi amigo, Edward Villegas, terminé pensando acerca del uso de cambio de variables en la solución de problemas de valores propios con diferencias finitas.

Mi enfoque usual

Mi enfoque usual es aproximar la segunda derivada con una diferencia centrada para el punto \(x_i\), de la siguiente manera

\begin{equation*} \frac{\mathrm{d}^2 f(x)}{\mathrm{d}x^2} \approx \frac{f(x + \Delta x) - 2 f(x_i) + f(x - \Delta x)}{\Delta x^2}\, , \end{equation*}

con \(\Delta x\) la separación entre puntos consecutivos.

Podemos solucionar este problem en Python con el siguiente bloque de código:

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigs


def regular_FD(pot, npts=101, x_max=10, nvals=6):
    """
    Find eigenvalues/eigenfunctions for Schrodinger
    equation for the given potential `pot` using
    finite differences
    """
    x = np.linspace(-x_max, x_max, npts)
    dx = x[1] - x[0]
    D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts, npts))/dx**2
    V = diags(pot(x))
    H = -0.5*D2 + V
    vals, vecs = eigs(H, k=nvals, which="SM")
    return x, np.real(vals), vecs

Configuremos los gráficos con el siguiente bloque de código.

# Jupyter notebook plotting setup & imports
%matplotlib notebook
import matplotlib.pyplot as plt

gray = '#757575'
plt.rcParams["figure.figsize"] = 6, 4
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["text.color"] = gray
fontsize = plt.rcParams["font.size"] = 12
plt.rcParams["xtick.color"] = gray
plt.rcParams["ytick.color"] = gray
plt.rcParams["axes.labelcolor"] = gray
plt.rcParams["axes.edgecolor"] = gray
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False

Consideremos el oscilador armónico cuántico, que tiene como valores propios

\begin{equation*} E_n = n + \frac{1}{2},\quad \forall n = 0, 1, 2 \cdots \end{equation*}

Usando el método de diferencias finitas obtenemos valores que están muy cerca de los analíticos.

x, vals, vecs = regular_FD(lambda x: 0.5*x**2, npts=201)
vals

Con respuesta

array([0.4996873 , 1.49843574, 2.49593063, 3.49216962, 4.48715031,
       5.4808703 ])

Los valores analíticos son los siguientes

[0.5, 1.5, 2.5, 3.5, 4.5, 5.5])

Si graficamos estos dos conjuntos, obtenemos lo siguiente.

plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.legend(["Analytic", "Finite differences"])
plt.tight_layout();

Veamos las funciones propias

plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();

Un inconveniente con este método es el muestreo redundante hacia los extremos del intervalo mientras submuestreamos el centro.

Transformando el dominio

Ahora, consideremos el caso en el que transormamos el dominio infinito a uno finito usando un cambio de variable

\begin{equation*} \xi = \xi(x) \end{equation*}

con \(\xi \in (-1, 1)\). Dos opciones para esta transformación son:

• \(\xi = \tanh x\); y

• \(\xi = \frac{2}{\pi} \arctan x\).

Haciendo este cambio de variable la ecuación, debemos resolver la siguiente ecuación

\begin{equation*} -\frac{1}{2}\left(\frac{\mathrm{d}\xi}{\mathrm{d}x}\right)^2\frac{\mathrm{d}^2}{\mathrm{d}\xi^2}\psi(\xi) - \frac{1}{2}\frac{\mathrm{d}^2\xi}{\mathrm{d}x^2}\frac{\mathrm{d}}{\mathrm{d}\xi}\psi(\xi) + V(\xi) \psi(\xi) = E\psi(\xi) \end{equation*}

El siguiente bloque de código resuelve el problema de valores propios en el dominio transformado:

def mapped_FD(pot, fun, dxdxi, dxdxi2, npts=101, nvals=6, xi_tol=1e-6):
    """
    Find eigenvalues/eigenfunctions for Schrodinger
    equation for the given potential `pot` using
    finite differences over a mapped domain on (-1, 1)
    """
    xi = np.linspace(-1 + xi_tol, 1 - xi_tol, npts)
    x = fun(xi)
    dxi = xi[1] - xi[0]
    D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts, npts))/dxi**2
    D1 = 0.5*diags([-1, 1], [-1, 1], shape=(npts, npts))/dxi
    V = diags(pot(x))
    fac1 = diags(dxdxi(xi)**2)
    fac2 = diags(dxdxi2(xi))
    H = -0.5*fac1.dot(D2) - 0.5*fac2.dot(D1) + V
    vals, vecs = eigs(H, k=nvals, which="SM")
    return x, np.real(vals), vecs

Primera transformación: \(\xi = \tanh(x)\)

Consideremos la primera transformación

\begin{equation*} \xi = \tanh(x)\, . \end{equation*}

En este caso,

\begin{equation*} \frac{\mathrm{d}\xi}{\mathrm{d}x} = 1 - \tanh^2(x) = 1 - \xi^2\, , \end{equation*}

y

\begin{equation*} \frac{\mathrm{d}^2\xi}{\mathrm{d}x^2} = -2\tanh(x)[1 - \tanh^2(x)] = -2\xi[1 - \xi^2]\, . \end{equation*}

Necesitamos definir las funciones

pot = lambda x: 0.5*x**2
fun = lambda xi: np.arctanh(xi)
dxdxi = lambda xi: 1 - xi**2
dxdxi2 = lambda xi: -2*xi*(1 - xi**2)

y correr

x, vals, vecs = mapped_FD(pot, fun, dxdxi, dxdxi2, npts=201)
vals

Y obtenemos los siguientes valores propios

array([0.49989989, 1.4984226 , 2.49003572, 3.46934257, 4.46935021,
       5.59552989])

Si los comparamos con los valores analítivos obtenemos lo siguiente.

plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.legend(["Analytic", "Finite differences"])
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.tight_layout();

Y las siguientes funciones propias.

plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();

Segunda transformación: \(\xi = \frac{2}{\pi}\mathrm{atan}(x)\)

Consideremos ahora la segunda transformación

\begin{equation*} \xi = \frac{2}{\pi}\mathrm{atan}(x)\, . \end{equation*}

En este caso,

\begin{equation*} \frac{\mathrm{d}\xi}{\mathrm{d}x} = \frac{2}{\pi(1 + x^2)} = \frac{2 \cos^2\xi}{\pi} \, , \end{equation*}

y

\begin{equation*} \frac{\mathrm{d}^2\xi}{\mathrm{d}x^2} = -\frac{4x}{\pi(1 + x^2)^2} = -\frac{4 \cos^4\xi \tan \xi}{\pi}\, . \end{equation*}

Una vez más, definimos las funciones

pot = lambda x: 0.5*x**2
fun = lambda xi: np.tan(0.5*np.pi*xi)
dxdxi = lambda xi: 2/np.pi * np.cos(0.5*np.pi*xi)**2
dxdxi2 = lambda xi: -4/np.pi * np.cos(0.5*np.pi*xi)**4 * np.tan(0.5*np.pi*xi)

y ejecutamos

x, vals, vecs = mapped_FD(pot, fun, dxdxi, dxdxi2, npts=201)
vals

para obtener los siguientes valores propios

array([0.49997815, 1.49979632, 2.49930872, 3.49824697, 4.49627555,
       5.49295665])

con la siguiente gráfica

plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.legend(["Analytic", "Finite differences"])
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.tight_layout();

y las siguientes funciones propias.

plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi|^2$", fontsize=16)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();

Comparación

Para todos los ejemplos se importaron los siguientes módulos:

from __future__ import division, print_function
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap
from colorspacious import cspace_convert

Imagen de prueba

Peter Kovesi propuso una manera de comparar mapas de colores cíclicos en un artículo en 2015. Consta de una rampa en espiral con ondulaciones. En coordenadas polares es

\begin{equation*} c = (A \rho^p \sin(n \theta) + \theta)\, \mathrm{mod}\, 2\pi \end{equation*}

con \(A\) la amplitud de la oscilación, \(\rho\) el radio normalizado en [0, 1], \(p\) un número positivo, y \(n\) el número de ciclos.

Y la siguiente función crea la rejilla en Python.

def circle_sine_ramp(r_max=1, r_min=0.3, amp=np.pi/5, cycles=50,
                     power=2, nr=50, ntheta=1025):
r, t = np.mgrid[r_min:r_max:nr*1j, 0:2*np.pi:ntheta*1j]
r_norm = (r - r_min)/(r_max - r_min)
vals = amp * r_norm**power * np.sin(cycles*t) + t
vals = np.mod(vals, 2*np.pi)
return t, r, vals

El resultado es el siguiente

Prueba para daltonismo

t, r, vals = circle_sine_ramp(cycles=0)
cmaps = ["hsv",
         "cmocean_phase",
         "hue_L60",
         "erdc_iceFire",
         "nic_Edge",
         "colorwheel",
         "cyclic_mrybm",
         "cyclic_mygbm"]
severity = [0, 50, 50, 50]
for colormap in cmaps:
    data = np.loadtxt(colormap + ".txt")
    fig = plt.figure()
    for cont in range(4):
        cvd_space = {"name": "sRGB1+CVD",
             "cvd_type": cvd_type[cont],
             "severity": severity[cont]}
        data2 = cspace_convert(data, cvd_space, "sRGB1")
        data2 = np.clip(data2, 0, 1)
        cmap = LinearSegmentedColormap.from_list('my_colormap', data2)
        ax = plt.subplot(2, 2, 1 + cont, projection="polar")
        ax.pcolormesh(t, r, vals, cmap=cmap)
        ax.set_xticks([])
        ax.set_yticks([])
    plt.suptitle(colormap)
    plt.tight_layout()
    plt.savefig(colormap + ".png", dpi=300)

Mapas de colores cíclicos generados aleatoriamente

¿Qué pasa si generamos mapas de colores cíclicos aleatoriamente? ¿Cómo lucirían?

A continuación están las piezas de código y mapas de colores resultantes.

from __future__ import division, print_function
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.colors import LinearSegmentedColormap


# Next line to silence pyflakes. This import is needed.
Axes3D

plt.close("all")
fig = plt.figure()
fig2 = plt.figure()
nx = 4
ny = 4
azimuths = np.arange(0, 361, 1)
zeniths = np.arange(30, 70, 1)
values = azimuths * np.ones((30, 361))
for cont in range(nx * ny):
    np.random.seed(seed=cont)
    rad = np.random.uniform(0.1, 0.5)
    center = np.random.uniform(rad, 1 - rad, size=(3, 1))
    mat = np.random.rand(3, 3)
    rot_mat, _ = np.linalg.qr(mat)
    t = np.linspace(0, 2*np.pi, 256)
    x = rad*np.cos(t)
    y = rad*np.sin(t)
    z = 0.0*np.cos(t)
    X = np.vstack((x, y, z))
    X = rot_mat.dot(X) + center

    ax = fig.add_subplot(ny, nx, 1 + cont, projection='polar')
    cmap = LinearSegmentedColormap.from_list('my_colormap', X.T)
    ax.pcolormesh(azimuths*np.pi/180.0, zeniths, values, cmap=cmap)
    ax.set_xticks([])
    ax.set_yticks([])

    ax2 = fig2.add_subplot(ny, nx, 1 + cont, projection='3d')
    ax2.plot(X[0, :], X[1, :], X[2, :])
    ax2.set_xlim(0, 1)
    ax2.set_ylim(0, 1)
    ax2.set_zlim(0, 1)
    ax2.view_init(30, -60)
    ax2.set_xticks([0, 0.5, 1.0])
    ax2.set_yticks([0, 0.5, 1.0])
    ax2.set_zticks([0, 0.5, 1.0])
    ax2.set_xticklabels([])
    ax2.set_yticklabels([])
    ax2.set_zticklabels([])


fig.savefig("random_cmaps.png", dpi=300, transparent=True)
fig2.savefig("random_cmaps_traj.svg", transparent=True)

Una idea interesante sería tomar estos mapas de colores y optimizar algunos parámetros perceptuales como luminosidad para obtener mapas de colores útiles.

Conclusiones

Probablemente, yo usaría phase, colorwheel, o mrybm en el futuro.

Referencias

1. Peter Kovesi. Good Colour Maps: How to Design Them. arXiv:1509.03700 [cs.GR] 2015

El problema

Entremos en materia. El problema fue originalmente parte de la Competencia 53 de Putnam en 1992. Su planteamiento es el siguiente:

Se eligen cuatro puntos al azar en la superficie de un esfera. ¿Cuál es la probabilidad de que el centro de la la esfera se encuentra dentro del tetraedro cuyos vértices están en los cuatro puntos? (Se entiende que cada punto es elegido de forma dependiente en relación con una distribución uniforme en la esfera.)

Como se muestra en el video mencionado, la probabilidad es \(1/8\). Vamos a escribir un algoritmo para obtener este resultado, aproximadamente, al menos.

El enfoque propuesto.

El enfoque que vamos a utilizar es bastante sencillo. Vamos a obtener una muestra de conjuntos aleatorios (independientes), con cuatro puntos cada uno, y vamos a verificar cuántos satisfacen la condición de estar adentro del tetraedro con los puntos como vértices.

Para que este enfoque funcione, necesitamos dos cosas:

1. Una forma de generar números aleatorios distribuidos uniformemente. Esto ya lo tenemos en numpy.random.uniform, por lo que no necesitamos hacer mucho al respecto.

2. Una forma de verificar si un punto está dentro de un tetraedro.

El algoritmo

A continuación se muestra una implementación de Python del enfoque discutido anteriormente

from __future__ import division, print_function
from numpy import (random, pi, cos, sin, sign, hstack,
                   column_stack, logspace)
from numpy.linalg import det
import matplotlib.pyplot as plt


def in_tet(x, y, z, pt):
    """
    Determina si un punto pt está al interior de
    un tetraedro con vértices x, y, z
    """
    mat0 = column_stack((x, y, z, [1, 1, 1, 1]))
    det0 = det(mat0)
    for cont in range(4):
        mat = mat0.copy()
        mat[cont] = hstack((pt, 1))
        if sign(det(mat)*det0) < 0:
            inside = False
            break
        else:
            inside = True
    return inside


#%% Cálculo
prob = []
random.seed(seed=2)
N_min = 1
N_max = 5
N_vals = logspace(N_min, N_max, 100, dtype=int)
for N in N_vals:
    inside_cont = 0
    for cont_pts in range(N):
        phi = random.uniform(low=0.0, high=2*pi, size=4)
        theta = random.uniform(low=0.0, high=pi, size=4)
        x = sin(theta)*cos(phi)
        y = sin(theta)*sin(phi)
        z = cos(theta)
        if in_tet(x, y, z, [0, 0, 0]):
            inside_cont += 1

    prob.append(inside_cont/N)


#%% Graficación
plt.figure(figsize=(4, 3))
plt.hlines(0.125, 10**N_min, 10**N_max, color="#3f3f3f")
plt.semilogx(N_vals, prob, "o", alpha=0.5)
plt.xlabel("Number of trials")
plt.ylabel("Computed probability")
plt.tight_layout()
plt.show()

Como se esperaba, cuando el número de muestras es suficientemente grande, la probabilidad estimada es cercana al valor teórico: 0,125. Esto se puede ver en la siguiente figura.