montecarlo

from math import pi, sqrt
from functools import reduce
from operator import mul
import random
import pandas as pd


points = {'x': [], 'y': [], 'z': []}


class Ball:
    def __init__(self, coords, radius):
        self.coords = coords
        self.radius = radius

    def __contains__(self, point):
        return sum((a-b)**2 for a, b in zip(self.coords, point)) <= self.radius ** 2

    def aabb(self):
        low = tuple(c-self.radius for c in self.coords)
        high = tuple(c+self.radius for c in self.coords)
        return low, high

    def volume(self):
        return 4/3 * pi * self.radius ** 3


class Torus:
    # R - расстояние от оси вращения до образующей,
    # r - радиус образующей
    def __init__(self, coords, R, r):
        self.coords = coords
        self.R = R
        self.r = r

    def __contains__(self, point):
        diff = tuple(p-s for s, p in zip(self.coords, point))
        Rs = sqrt(diff[0]**2 + diff[1]**2)
        Rr = self.R - Rs
        return diff[2]**2 + Rr**2 <= self.r ** 2

    def aabb(self):
        c = self.coords
        low = (c[0]-self.r-self.R, c[1]-self.r-self.R, c[2]-self.r)
        high = (c[0]+self.r+self.R, c[1]+self.r+self.R, c[2]+self.r)
        return low, high

    def volume(self):
        return 2 * pi * pi * self.R * self.r * self.r


class Pair:
    def __init__(self, f1, f2):
        self.f1 = f1
        self.f2 = f2

    def __contains__(self, point):
        return point in self.f1 or point in self.f2

    def aabb(self):
        a1 = self.f1.aabb()
        a2 = self.f2.aabb()
        low = tuple(min(a, b) for a, b in zip(a1[0], a2[0]))
        high = tuple(max(a, b) for a, b in zip(a1[1], a2[1]))
        return low, high

    def volume(self):
        return self.f1.volume() + self.f2.volume()

def make_point(aabb, fig):
    def foo(x, i):
        low = aabb[0][i]
        high = aabb[1][i]
        return low+(high-low)*x

    r = random.random
    p = tuple(foo(r(),i) for i in range(3))

    global points
    if p in fig:
        points['x'].append(p[0])
        points['y'].append(p[1])
        points['z'].append(p[2])

    return p

def run_tests(box, fig, n):
    S = reduce(mul, (box[1][c]-box[0][c] for c in range(3)), 1.0)
    print(f'Объём коробки: {S}')
    Ng = sum(make_point(box, fig) in fig for i in range(n))
    print(f'Удачных поподаний: {Ng}, всего поподаний: {N}')
    Sf = Ng / N * S

    print(f'СКО: {sqrt((N-Ng)/(N*Ng))}')
    print(f'Теоретический объём: {F.volume()}, вычесленный объём: {Sf}')


if __name__ == "__main__":
    random.seed()

    ball_radius = 1.0
    ball_c = (0.0, 0.0, 0.0)

    torus_R = 3.0
    torus_r = 1.0
    torus_c = (0.0, 0.0, 6.0)

    B = Ball(ball_c, ball_radius)
    T = Torus(torus_c, torus_R, torus_r)
    F = Pair(B, T)

    AABB = F.aabb()
    print(f'границы проверяемой коробки{AABB}')


    N=10**4

    run_tests(AABB, F, N)


    from mpl_toolkits.mplot3d import Axes3D
    import matplotlib.pyplot as plt
    import numpy as np


    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points['x'], points['y'], points['z'], marker='o')
    ax.set_xlabel('X ')
    ax.set_ylabel('Y ')
    ax.set_zlabel('Z ')

    plt.show()