import sys
import math
import random
import matplotlib.pyplot as plt

# Calculate the area and mass of an irregular 2D region defined by
# multiple constraints.

def sigma(x, y):
    return x**2+4*y**2

N = 10000
seed = 42

if len(sys.argv) > 1: N = int(sys.argv[1])
if len(sys.argv) > 2: seed = int(sys.argv[2])

random.seed(seed)

sum0 = 0
sum1 = 0
sum2 = 0
xp = []
yp = []
for i in range(N):
    x = random.uniform(-1,1)
    y = random.uniform(-1,1)

    # Retain points inside the unit circle, below one line, above
    # another line, and to the left of a third.
    
    if x**2 + y**2 <= 1 \
       and y**2 < x \
       and x+y < 0.5:
        xp.append(x)
        yp.append(y)
        sum0 += 1
        s = sigma(x,y)
        sum1 += s
        sum2 += s**2

area = 4*float(sum0)/N
print('area =', area)
sbar = sum1/sum0
sbar2 = sum2/sum0
print('mass =', area*sbar , '+-', \
       area*sbar/math.sqrt(sum0) + area*math.sqrt((sbar2-sbar**2)/sum0))

plt.xlim(-1,1)
plt.ylim(-1,1)
ax = plt.gca()
ax.set_aspect('equal')
plt.plot(xp, yp, 'b.', markersize=0.5)
plt.show()