#! /usr/bin/env python
# -*- coding: utf-8 -*-
from __future__ import print_function
#import matplotlib # uncomment these two lines if the program doesn't
#matplotlib.use("TkAgg") # throw any errors but exits without having shown any plots.
import numpy as np
import matplotlib.pyplot as plt
### Task 10.1.1:
print("\nAnswer to task 10.1.5:")
print("----------------------")
print("\tThe two oppositely charged prolate ellipsoids will try to minimize their interaction energy.\n\
\tSince the energy is proportional to the negative reciprocal average distance between the\n\
\tellipsoids, the total interaction energy is minimal when the distance between their centers\n\
\tis minimal.\n\
\tDue to the fact that they are solid, and thus, cannot interpenetrate, the ellipsoids will\n\
\teventually align in a side-by-side arrangement.\n")
### Task 10.1.2:
def cuboid(n, dimensions):
return np.random.uniform(low=-1.0, high=1.0, size=(n, 3)) * dimensions
### Task 10.1.3:
def ellipsoid(n, dimensions):
# first, we set up an empty array in which we will store the generated points:
points = np.empty((n, 3))
# then, we generate random points in the bounding box until n of them lie inside the ellipsoid:
n_matches = 0
while n_matches < n:
n_matches_left = n - n_matches
# since we know that the probability for a random point in the bounding box
# to lie within the ellipsoid is pi/3 (which is approximately 0.5),
# we generate twice as many points as required:
candidates = cuboid(2 * n_matches_left, dimensions)
# now determine the ones inside the ellipsoid:
matches = candidates[np.where(np.sum(candidates ** 2 / dimensions ** 2, axis=-1) < 1.0), :].reshape((-1, 3))
new_n_matches = n_matches + matches.shape[0]
# finally, we add the matching points to the output array while making sure we don't add too many:
points[n_matches:min(new_n_matches, n), :] = matches[:min(matches.shape[0], n_matches_left), :]
n_matches = new_n_matches
return points
### Task 10.1.4:
def compute_energy(d, n, shape, dimensions):
try:
dummy = d[0]
except:
d = np.array([d])
energy = np.empty_like(d)
points1 = shape(n, dimensions)
points2 = shape(n, dimensions)
for i in range(len(d)):
r1 = points1
r2 = points2.copy()
r2[:, 0] += d[i]
energy[i] = -np.mean(1.0/np.sqrt(np.sum((r2 - r1) ** 2, axis=-1)))
return energy
### Task 10.1.5:
n_samples = 10000
n_distances = 20
d_max = 5.0
r = 0.5
dimensions = np.array([r, r, r])
distance_cubes = np.linspace(2*dimensions[0]*1.001, d_max, n_distances)
energy_cubes = compute_energy(distance_cubes, n_samples, ellipsoid, dimensions)
r = (3.0 / (4.0 * np.pi)) ** (1.0 / 3.0)
dimensions = np.array([r, r, r])
distance_spheres = np.linspace(2*dimensions[0]*1.001, d_max, n_distances)
energy_spheres = compute_energy(distance_spheres, n_samples, ellipsoid, dimensions)
r = (3.0 / (8.0 * np.pi)) ** (1.0 / 3.0)
dimensions = np.array([2 * r, r, r])
distance_ellipsoids_t2t = np.linspace(2*dimensions[0]*1.001, d_max, n_distances)
energy_ellipsoids_t2t = compute_energy(distance_ellipsoids_t2t, n_samples, ellipsoid, dimensions)
r = (3.0 / (8.0 * np.pi)) ** (1.0 / 3.0)
dimensions = np.array([r, r, 2 * r])
distance_ellipsoids_sbs = np.linspace(2*dimensions[0]*1.001, d_max, n_distances)
energy_ellipsoids_sbs = compute_energy(distance_ellipsoids_sbs, n_samples, ellipsoid, dimensions)
plt.figure()
plt.subplot(1, 1, 1)
plt.plot(distance_cubes, energy_cubes, "o-", label="cubes")
plt.plot(distance_spheres, energy_spheres, "o-", label="spheres")
plt.plot(distance_ellipsoids_t2t, energy_ellipsoids_t2t, "o-", label="tip-to-tip ellipsoids")
plt.plot(distance_ellipsoids_sbs, energy_ellipsoids_sbs, "o-", label="side-by-side ellipsoids")
plt.plot(distance_ellipsoids_sbs, -1.0/distance_ellipsoids_sbs, "o-", label="point charges")
plt.xlabel("distance")
plt.ylabel("energy")
plt.legend()
plt.show()
print("Answer to task 10.1.5:")
print("----------------------")
print("\tClearly, the tip-to-tip arrangement of the two prolate ellipsoids has a lower interaction energy\n\
\tthan all other arrangements. This can be explained by the fact that due to the -1/d relationship, the\n\
\tcharges at those two tips that face each other more than compensate the larger distance between those\n\
\tcharges that lie in the tips that are furthest away from each other.\n\
\tThis also holds true for the ellipsoids in task 10.1.1 as long as they don't touch. Once they do, the\n\
\tonly way to further minimize their interaction energy is to rotate and align parallel to each other.\n\
\tFrom a classical point of view, if the ellipsoids could interpenetrate, they would do so without any\n\
\trotation, with their interaction energy diverging to negative infinity.\n")