###################################################################

# This file is unit04_elliptic_integral.py. It does a numerical
# integration to approximate the value of a
# complete elliptic integral of the second kind. This integral gives
# the perimeter of an ellipse with semi-major, minor axes a and b.

# George Gollin, University of Illinois, May 26, 2016

###################################################################

# import libraries

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

# value I obtain when running 100,000,000 terms
integral_100_million = 11.0517460961

# semimajor and minor axes:

axis_a = 2.0
axis_b = 1.5

# squares:

axis_asq = axis_a**2
axis_bsq = axis_b**2

# integral will go from 0 to pi/2. define an array holding theta values.

theta_min = 0.
theta_max = np.pi / 2

number_of_slices = 1

while number_of_slices <= 1024:
    
    number_of_slices = number_of_slices * 2

    array_length = number_of_slices + 1
    theta_array = np.linspace(theta_min, theta_max, array_length)

    # also get the interval width.
    dtheta = theta_array[1] - theta_array[0]

    # load the array of ordinate (y axis) values:

    f_of_theta_array = np.sqrt(axis_asq - (axis_asq - axis_bsq) * np.sin(theta_array)**2)

    # now sum the elements in the array. Note that we do NOT want to include the 
    # very last element in the array, since it is a right edge. Handle this by 
    # resizing the array.

    f_of_theta_array = np.resize(f_of_theta_array, number_of_slices)
    f_of_theta_array_sum = np.sum(f_of_theta_array)

    # now multiply by 4 and by dtheta to finish.

    integral = 4 * dtheta * f_of_theta_array_sum

    # now print what we found

    print("complete elliptic integral of the second kind for a, b = ", axis_a, axis_b)
    print("number of terms: ", number_of_slices, "  integral = ", integral, \
    " error: ", integral - integral_100_million)

# end of loop!

# now do a trapezoidal rule integration.

number_of_slices = 1

while number_of_slices <= 1024:
    
    # number_of_slices = number_of_slices * 2

    array_length = number_of_slices + 1
    theta_array = np.linspace(theta_min, theta_max, array_length)

    # also get the interval width.
    dtheta = theta_array[1] - theta_array[0]

    # load the array of ordinate (y axis) values:

    f_of_theta_array = np.sqrt(axis_asq - (axis_asq - axis_bsq) * np.sin(theta_array)**2)

    # now halve the very first and very last elements in the array:

    f_of_theta_array[0] = f_of_theta_array[0] / 2
    f_of_theta_array[number_of_slices] = f_of_theta_array[number_of_slices] / 2

    # now sum the elements in the array. 

    f_of_theta_array_sum = np.sum(f_of_theta_array)

    # now multiply by 4 and by dtheta to finish.

    integral = 4 * dtheta * f_of_theta_array_sum

    # now print what we found

    print("complete elliptic integral of the second kind for a, b = ", axis_a, axis_b)
    print("trap rule terms: ", number_of_slices, "  integral = ", integral, \
    " error: ", integral - integral_100_million)

    number_of_slices = number_of_slices * 2

# end of loop!

###################################