#include <time.h>
#include <stdlib.h>
#include <cairo.h>
#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>

static double
sample_exp (double lambda)
{
    double u = drand48 ();

    return -log (1 - u) / lambda;
}

#define FALSE 0
#define TRUE 1

#define WIDTH 12000
#define HEIGHT 25
#define N_YEARS 20
#define DEATHS_PER_YEAR 7

#define BG_COLOR        1, 1, 1
#define CLUSTER_WIDTH   20.0
#define CLUSTER_COLOR   0.2, 0.5, 0.3, 0.4
#define LINE_WIDTH      1.0
#define LINE_COLOR      0, 0, 0
#define DEATH_WIDTH     1.0
#define DEATH_COLOR     0.0, 0.0, 0.0
#define DEATH_HEIGHT    (CLUSTER_WIDTH / HEIGHT)

static double *
generate_deaths (int n_years, double deaths_per_year, int *n_deaths)
{
    double w;
    double t;
    double *result = malloc (10 * sizeof (double));
    double *first;

    t = 0.0;
    *n_deaths = 0;

    w = sample_exp (deaths_per_year);
    t += w;
    while (t < n_years)
    {
        result[(*n_deaths)++] = t;
        result = realloc (result, ((*n_deaths) + 10) * sizeof (double));

        w = sample_exp (deaths_per_year);
        t += w;
    }

    return result;
}

int
main ()
{
    cairo_surface_t *surface =
        cairo_image_surface_create (CAIRO_FORMAT_ARGB32, WIDTH, HEIGHT);
    int i;
    double t = 0.0;
    cairo_t *cr = cairo_create (surface);
    double y = HEIGHT / 2.0;
    double w;
    int n_deaths;
    int cluster;
    double *deaths;
    int in_cluster;
    double cluster_x0;

    srand48 (time (0));

    deaths = generate_deaths (N_YEARS, DEATHS_PER_YEAR, &n_deaths);
    
    cairo_set_source_rgb (cr, BG_COLOR);
    cairo_paint (cr);

    cairo_save (cr);
    cairo_move_to (cr, (deaths[0] / N_YEARS) * WIDTH, y);

    for (i = 1; i < n_deaths + 1; ++i)
    {
        double d0 = deaths[i - 1];
        double d1 = i < n_deaths? deaths[i] : 0.0;
        double x1 = (d1 / N_YEARS) * WIDTH;

        if (d1 - d0 > 1.0 / DEATHS_PER_YEAR || i == n_deaths)
        {
            double x1 = (d0 / N_YEARS) * WIDTH;
            double x2 = (d1 / N_YEARS) * WIDTH;

            /* 0 length lines don't get square caps because
             * cairo would have to invent an orientation for
             * them. So add 0.002 to make sure the line has
             * a length.
             */
            cairo_line_to (cr, x1 + 0.002, y);
            cairo_set_line_width (cr, CLUSTER_WIDTH);
            cairo_set_source_rgba (cr, CLUSTER_COLOR);
            cairo_set_line_cap (cr, CAIRO_LINE_CAP_SQUARE);
            cairo_stroke (cr);
            cairo_move_to (cr, x2 - 0.002, y);
        }
    }

    cairo_restore (cr);

    for (i = 1; i < n_deaths; ++i)
    {
        double d0 = deaths[i - 1];
        double d1 = deaths[i];
        double x0 = (d0 / N_YEARS) * WIDTH;
        double x1 = (d1 / N_YEARS) * WIDTH;

        cairo_set_line_width (cr, LINE_WIDTH);
        cairo_set_source_rgb (cr, LINE_COLOR);
        cairo_move_to (cr, x0, y + 4);
        cairo_line_to (cr, x1, y + 4);
        cairo_stroke (cr);
    }

    for (i = 0; i < n_deaths; ++i)
    {
        double d = deaths[i];
        double x = (d / N_YEARS) * WIDTH;

        cairo_set_line_width (cr, DEATH_WIDTH);
        cairo_move_to (cr, x, (0.5 - DEATH_HEIGHT / 2.0) * HEIGHT);
        cairo_line_to (cr, x, (0.5 + DEATH_HEIGHT / 2.0) * HEIGHT);
        cairo_set_source_rgb (cr, DEATH_COLOR);
        cairo_stroke (cr);
    }

    cairo_surface_write_to_png (surface, "diagram.png");

    return 0;
}