# -*- coding: utf-8 -*-
"""
Very simple code to throw darts at a circle to find 'pi'
@author: neil 2022
"""

import numpy as np
# Very straighforward soln to this weeks lab work (some of you are making it too difficult)
# lets integrate the area of a diameter 1 circle using the trapiziod and using basic Monte Carlo
# we only need to integrate a 1/4 circle (x,y > 0) and multiply by 4, that way we stay away from sign problems

# 1-Analytic soln, for a circle of radius 1, it is just 'pi'
print('Analytic solution of area of a circle (r=1), {:.5}'.format( np.pi ))

# 2-numerical integration, trapizoid rule, 10 points, for a circle, y = sqrt(r**2-x^2)
slices = 10     # number of trapizoids
r      = 1.0    # distance from 0 to max radius
delx   = r/(slices) #width of trapizoids
area   = 0      # place to keep the area calcs
for x in np.arange(slices):
    x1 = x*delx
    x2 = x1 + delx
    area += delx*(np.sqrt(1. - x1**2) + np.sqrt(1. - x2**2))/2        
print('Numerical trapizoidal integration area,      {:.5}'.format(4*area))
    
# 3-there are several approaches to MC integration, here we throw darts (random guesses)
# at a circle inscribed in a square and calculate the number of hits vs the total throws, which will approx
# the ratio of the size of the circle to the size of the square.

# Define number of samples (darts)
N = 1000
# We only use the quadrant x,y > 0 (quadrant's area = 1)
# Sample from uniform distribution on the interval 0-1
mc = np.random.rand(N,2)
# Define the points under the curve
points_under = 0
for i in range(len(mc)):
    if mc[i,1] <= np.sqrt(1-mc[i,0]**2) :
        points_under += 1
    
# Calculate integral by taking the ratio of points in and out to the area of the quadrant of the 
# unit square
integral = points_under/ N
print("Monte Carlo integration of unit circle area: {:.5}".format(integral*4))