#!/usr/bin/python
# -*- coding: utf8 -*-

import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from math import *

code_website = 'http://commons.wikimedia.org/wiki/User:Geek3/mplwp'
try:
    import mplwp
except ImportError, er:
    print 'ImportError:', er
    print 'You need to download mplwp.py from', code_website
    exit(1)

name = 'mplwp_Dickman_function.svg'
fig = mplwp.fig_standard(mpl)

xlim = 0, 5; fig.gca().set_xlim(xlim)
ylim = -0.1, 1.1; fig.gca().set_ylim(ylim)
mplwp.mark_axeszero(fig.gca())

n = 250 * 4 # multiple of 4
dx = 1. / n
x = np.arange(xlim[0], xlim[1] + dx/2., dx)
y = np.empty_like(x)

for i, xx in enumerate(x):
    if xx <= 1.0:
        y[i] = 1.0
    elif xx <= 2.0:
        y[i] = 1.0 - log(xx)
    else:
        # Simpson's rule to integrate x * y(x) = \int_{x-1}^x y(x) dx
        s = y[i-n]
        s += 2. * np.sum(y[i-n+2:i-1:2])
        s += 4. * np.sum(y[i-n+1:i:2])
        y[i] = s / (3. * xx / dx - 1.)

plt.plot(x, y, label='Dickman function')

plt.legend(loc='upper right').get_frame().set_alpha(0.9)
plt.savefig(name)
mplwp.postprocess(name)