python - punto - formula de haversine ejemplos

¿Puntos que siguen una distribución de densidad 1/R en una cuadrícula XY?

Question

Tengo una cuadrícula XY con algunos puntos de la grilla que tienen ciertos valores asignados a ellos, en este caso, cada valor significa una cierta masa, por lo que básicamente las masas de puntos en una cuadrícula. Ahora quiero obtener un conjunto de puntos que siguen una distribución de densidad de 1 / R, donde R es la distancia desde el centro, entonces R = sqrt (x ^ 2 + y ^ 2). Por distribución de densidad, quiero decir que el número de puntos debe caerse como 1 / R. ¿Cómo voy a codificar esto?

Mi código está abajo:

import numpy as np
x = np.linspace(-50,50,100)
y = np.linspace(-50,50,100)
X, Y = np.meshgrid(x,y)
zeta_a = (25,25)
zeta_b = (-10,5) 
M_a = 150
M_b = 150 

Los zeta_a y zeta_b corresponden a masas de 2 puntos que tienen masas de 150 unidades. También necesito realizar cálculos de seguimiento usando estos puntos, así que también me gustaría saber cómo usar un formato más general en lugar de usar ''a'', ''b'', ''c'' para masas de n puntos.

Gracias por tu ayuda.

Answer 1

Suponiendo que entendí su pregunta (si no se aceptan comentarios):

La forma de crear cualquier distribución dada es mediante la interpolación sobre el inverso de la distribución CDF. Esta es mi función para hacerlo:

import numpy as np
import matplotlib.pyplot as plt

def randdist(PDF, x, n):
    """Create a distribution following PDF(x). PDF and x
    must be of the same length. n is the number of samples."""
    fp = np.random.rand(n,)
    CDF = np.cumsum(PDF)
    return np.interp(fp, CDF, x)

Ahora, en su caso, vamos a trabajar en coordenadas polares con R distribuido como 1 / ry Theta distribuidos uniformemente:

num = 1000   # The number of points
r = np.linspace(-50, 50, 100)
PDF = np.abs(1/r)
PDF = PDF/np.sum(PDF)    # PDF should be normalized
R = randdist(PDF, r, num)
Theta = 2*np.pi*np.random.rand(num,)

Ahora vamos a crear los vectores de puntos x y y

x = [R[k]*np.cos(Theta[k]) for k in range(num)]
y = [R[k]*np.sin(Theta[k]) for k in range(num)]

Para trazar

plot(x,y,''.'')

Tenga en cuenta que en mi respuesta hay un límite de rigidez en r = 50. Hay formas de superar esto, pero por ahora lo dejo como está.

Ahora parece que también desea insertar los puntos dentro de una cuadrícula 2D, como un histograma. Puedes hacer eso usando

z, _, _ = np.histogram2d(x, y, [100, 100])