/* Numerical solution to second Baseball Markov chain model */
/* with BOTH types of double plays allowed */

/* Written by Carey Williamson, January 2019 */

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

/* Transition rates for Markov chain model */
#define LAMBDA1 0.20   /* transition rate to get first out */
#define LAMBDA2 0.25   /* transition rate to get second out */
#define LAMBDA3 0.25   /* transition rate to get third out */
#define LAMBDA4 0.50   /* transition rate for TV commercial */

#define DEBUG 1   /* verbose debugging */

int main ()
  {
    float p0, p1, p2, p3;
    float sum, ratio;
    float lambdaD;
    float x, y, z;   

    printf("Enter double-play rate (lambdaD): ");
    scanf("%f", &lambdaD);

    /* start with p1 and compute others from it */
    x = 1.0/LAMBDA1;
    x += 1.0/LAMBDA2;
    x += 1.0/LAMBDA3;
    x += 1.0/LAMBDA4;
    x *= LAMBDA2;
#ifdef DEBUG
    printf("x is %f\n", x);
#endif

    y = 1.0/LAMBDA1;
    y += 1.0/LAMBDA4;
    y *= lambdaD;
#ifdef DEBUG
    printf("y is %f\n", y);
#endif

    z = 1.0/(LAMBDA1 * LAMBDA3);
    z += 1.0/(LAMBDA1 * LAMBDA4);
    z *= lambdaD;
    z *= (LAMBDA2 + lambdaD);
#ifdef DEBUG
    printf("z is %f\n", z);
#endif

    p1 = 1.0 / (x + y + z);
#ifdef DEBUG
    printf("p1 is %f\n", p1);
#endif

    ratio = (LAMBDA2 + lambdaD) / LAMBDA1;
    p0 = ratio * p1;
#ifdef DEBUG
    printf("p0 is %f times p1, so p0 = %f\n", ratio, p0);
#endif

    p2 = (lambdaD * p0 + LAMBDA2 * p1) / LAMBDA3;
#ifdef DEBUG
    printf("p2 is %f\n", p2);
#endif

    ratio = (LAMBDA1 + lambdaD) / LAMBDA4;
    p3 = ratio * p0;
#ifdef DEBUG
    printf("p3 is %f times p0, so p3 = %f\n", ratio, p3);
#endif

    /* Renormalize the probabilities now */
    sum = p0 + p1 + p2 + p3;
#ifdef DEBUG
    printf("sum is %f\n", sum);
#endif
    p0 /= sum;
    p1 /= sum;
    p2 /= sum;
    p3 /= sum;
    
    printf("\n");
    printf("Analytical Markov Chain Model of Baseball with Both Double Plays\n");
    printf("\n");
    printf("  Prob(nobody out): p0 = %f\n", p0);
    printf("  Prob(one out):    p1 = %f\n", p1);
    printf("  Prob(two out):    p2 = %f\n", p2);
    printf("  Prob(three out):  p3 = %f\n", p3);
    printf("\n");
    printf("Sanity check: p0 + p1 + p2 + p3 = %f\n",
	   p0 + p1 + p2 + p3);
    printf("\n");
  }