Calculation of the deflection angle in a two-body system

#ifdef _WIN32
    #define _USE_MATH_DEFINES
#endif

#include <stdio.h>
#include <math.h>

const size_t MAX_DIVISIONS = 21;
const size_t MAX_ITERATIONS = 21;
const double PRECISION_DELTA = 0.00000001;
const double DELTA = 16.0;

double potential(double r)
{
    return 1.0/r;
}

double function(double r, double b, double E)
{  
    return 1.0 / (pow(r, 2) * sqrt(1 - (pow(b, 2) / pow(r, 2)) - (potential(r) / E)));
}

double finite_integral(double up, double to, double b, double E)
{
    double cur_step,
           step,
           prev_sum,
           partial_sum,
           f_x;
    size_t i;

    prev_sum = 0.0;

    for(i = 0; i != MAX_DIVISIONS; ++i)
    {
        step = (to-up)/pow(2, i);
        partial_sum = 0.0;

        for(cur_step = up + step; cur_step <= to; cur_step += step) {
            f_x = step * function(cur_step - (step/2.0), b, E);

            // when f_x is nan
            // if f_x is nan, f_x != f_x will be true
            if(f_x == f_x) {  
                partial_sum = partial_sum + f_x;
            }
        }

        if(prev_sum != 0.0 && fabs(prev_sum - partial_sum) < PRECISION_DELTA) break;
        prev_sum = partial_sum;
    }

    return partial_sum;
}

double integral_to_infinite(double a, double b, double E)
{
    double to,
           prev_res,
           partial_res;
    size_t i;

    prev_res = 0.0;

    for(i = 1; i != pow(2, MAX_ITERATIONS); ++i)
    {
        to = a + DELTA;
        partial_res = finite_integral(a, to, b, E);

        if(prev_res != 0.0 && partial_res != 0.0 && partial_res < PRECISION_DELTA)
        {
            prev_res += partial_res;
            break;
        }

        prev_res += partial_res;
        a = to;
    }

    return (2.0 * b) * prev_res;
}

double to_degrees(double radians) {
    return 180.0 - radians * (180.0 / M_PI);
}

from cffi import FFI
from matplotlib import pyplot as plt
from matplotlib import patches as mpatches
from ctypes import c_double
from math import asin, sqrt, degrees


def analytic_evaluation(var_b, var_e):
    """Analytic function of theta integral."""
    return 2.0 * asin(1 / (sqrt(1 + 4 * float(var_b) ** 2 * float(var_e) ** 2)))


def range_f(start, end, step):
    """Create a list of double values."""
    end = c_double(end)
    ffi = FFI()
    start = c_double(start)
    while start.value <= end.value:
        yield ffi.cast("double", start.value)
        start.value = start.value + c_double(step).value


def main():
    """Program main."""
    ffi = FFI()
    ffi.cdef("""
double integral_to_infinite(double a, double b, double E);
double to_degrees(double radians);

""")
    lib = ffi.verify("""
#include "tetaQuad.h"
""", include_dirs=["."], extra_compile_args=["-O3"])

    var_e = ffi.cast("double", 0.1)
    var_a = ffi.cast("double", 0.0)

    points_x = []
    points_y = []
    points_y_a = []

    print("Calculus may take some time...")

    for var_b in range_f(0.0, 100, 0.5):
        points_x.append(var_b)
        rad = lib.integral_to_infinite(var_a, var_b, var_e)
        theta = lib.to_degrees(abs(rad))
        analytic = degrees(analytic_evaluation(var_b, var_e))

        print("b = {}\ttheta = {}\tanalytic = {}".format(
            float(var_b), theta, analytic))

        points_y.append(theta)
        points_y_a.append(analytic)

    plt.ylabel('teta')
    plt.xlabel('b')
    plt.plot(points_x, points_y, 'ro')
    plt.plot(points_x, points_y_a, 'b--')
    red_line = mpatches.Patch(color='red', label='theta value')
    blue_line = mpatches.Patch(color='blue', label='analytic value')
    plt.legend(handles=[red_line, blue_line])
    plt.grid(True)
    plt.show()


if __name__ == '__main__':
    main()