/*

	Program to integrate survival function for evaluating e_x and abar_x

	See "Applying survival models to pensioner mortality data",	details of
	which can be found at http://www.richardsconsulting.co.uk/laws.html

	Copyright (c) 2008 Stephen Richards.  All rights reserved.

*/

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

#define COMMAND  0
#define ALPHA    1
#define BETA     2
#define EPSILON  3
#define RHO      4
#define X        5
#define INTEREST 6
#define ARGC     7

#define AGE_LIMIT 150	/*	Maximum age.	*/

int main(const int argc, const char * argv[])
{
  int rc = EXIT_FAILURE;

  if (argc != ARGC)
    fprintf(stderr, "Usage: integrate alpha beta epsilon rho x interest\n");
  else
  {
    const double alpha    = atof(argv[ALPHA  ]),
                 beta     = atof(argv[BETA   ]),
                 epsilon  = atof(argv[EPSILON]),
                 rho      = atof(argv[RHO    ]),
                 x        = atof(argv[X      ]),
                 interest = atof(argv[INTEREST]);

    double dt = 1.0, H, x_plus_t, e_x, a_x, tpx, sumdt, theta = log(1.0 + interest);

    while (dt >= 0.00001)
    {
      H     =
      e_x   =
      a_x   =
      sumdt = 0.0;
      x_plus_t = x;
      while (x_plus_t < AGE_LIMIT)
      {
        H += dt * (exp(epsilon) + exp(alpha + beta * x_plus_t)) / (1.0 + exp(alpha + rho + beta * x_plus_t));
	    tpx = exp(-H);
	    e_x += dt * tpx;
	    a_x += dt * tpx * exp(-theta * sumdt);
	    x_plus_t += dt;
	    sumdt += dt;
      }
      printf("e_x=%.2f and a_x=%.2f with dt=%f\n", e_x, a_x, dt);
      dt /= 10.0;
	}

    rc = EXIT_SUCCESS;
  }

  return rc;
}