/*

	Program to calculate force of mortality under reliability theory.

	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-line arguments */
#define COMMAND 0
#define LAMBDA 1	/* The constant force of failure for an element.			*/
#define N 2			/* The number of elements in a block. 						*/
#define M 3			/* The number of serial blocks in the system.				*/
#define P 4			/* The probability that an element is working at outset.	*/
#define TMAX 5		/* The maximum age for outputting the force of mortality.	*/
#define ARGC 6

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

  if (argc != ARGC)
    fprintf(stderr, "Usage: reliability lambda n m p tmax\n");
  else
  {
    const double lambda = atof(argv[LAMBDA]),
                 p      = atof(argv[P     ]),
                 tmax   = atof(argv[TMAX  ]);

    const unsigned int n = atoi(argv[N]),
                       m = atoi(argv[M]);


    if (lambda <= 0.0)
      fprintf(stderr, "Lambda must be greater than zero!\n");
    if ( (p > 1.0) || (p <= 0.0) )
      fprintf(stderr, "q cannot lie outside (0, 1]!\n");
    else if (tmax <= 0.0)
      fprintf(stderr, "tmax must be greater than zero!\n");
    else if (n <= 0)
      fprintf(stderr, "n must be greater than zero!\n");
    else if (m <= 0)
      fprintf(stderr, "m must be greater than zero!\n");
    else
    {
      double t, mu_t, c = 1.0, sum, factor;
      unsigned int i;

      printf("Age,Force of mortality\n");
      for (t=1.0; t<=tmax; t+=1.0)
      {
        sum = 0.0;
        factor = 1.0;
        for (i=1; i<=n; i++)
        {
          sum += factor * pow(1.0 - exp(-lambda * t), i-1.0) / (1.0 - pow(1.0 - exp(-lambda * t), i));
          factor *= n * p / i;
        }
        mu_t = n * p * m * c * lambda * exp(-p*n) * sum;
        printf("%.0f,%.6f\n", t, mu_t);
      }
    }
  }

  return rc;
}