#include <stdio.h>
#include <math.h>
#define PI 4. * atan(1.0)

main()
{
   int mu, n, terms;
   float kappa, epsilon, gamma, alpha, biot, lambda[2000];
   float phin, psin, r1dstar, rcstar, sum, sumold;
   float dn, la, le, lg, l2, k1, k2, a, b, dval;
   float beta, b2, b4, b6, b8, t1, t2, h;

   biot = h * b * 0.001 / k1;
   gamma = t1 / b;
   alpha = (t1+t2) / b;
   kappa = k1/k2;
   epsilon = a/b;

   mu -= 2.;
   n=0;
   sum = sumold = 0.;
   for (;;){
     n++;
     /* calculate roots using Stokes approximation */
     beta = PI * (4. * n + 1.);
     b2 = beta * beta;
     b4 = b2 * b2;
     b6 = b2 * b4;
     b8 = b2 * b6;
     lambda[n] = (beta/4.) * (1. - 6./b2 + 6./b4 - 4716./(5.*b6) + 3902418./(70.*b8));
     le = lambda[n] * epsilon;
     lg = lambda[n] * gamma;
     la = lambda[n] * alpha;
     l2 = lambda[n] * lambda[n];
     psin = (lambda[n] + biot * tanh(la)) / (lambda[n] * tanh(la) + biot);
     phin = ((kappa - 1.)/kappa) * (cosh(lg)) * (cosh(lg) - psin * sinh(lg));
     switch(mu)
        {
          case -1: dn = -2. * epsilon * sin(le) / (l2 * j0(lambda[n]) * j0(lambda[n]));
                   break;
          case 0 : dn = -2. * epsilon * j1(le) / (l2 * j0(lambda[n]) * j0(lambda[n]));
                   break;
          case  1: dn = (-2. * epsilon * sin(le) / (l2 * j0(lambda[n]) * j0(lambda[n])))
                        * (1./(le * le) - 1./(le * tan(le)));
                   break;
        }
 
     sum += dn * (phin * tanh(lg) - psin) /(1. - phin) * j1(le)/le;
     if (n == terms)
        break;
     sumold = sum;
     }
   r1dstar =  (epsilon * gamma + kappa * epsilon * (alpha - gamma)
              + epsilon / biot) / atan(1.0) ;
   rcstar = 8. * (mu + 1.)/(PI * epsilon) * sum;
   printf("r1dstar = %f  rcstar = %f\n",r1dstar,rcstar);
   dval = 4. * k1 * a;

   
   printf("\n=================================================\n");
   printf("  DIMENSIONLESS VALUES\n");
   printf("          1D resistance: %f\n",r1dstar);
   printf("constriction resistance: %f\n",rcstar);
   printf("       total resistance: %f\n",rcstar+r1dstar);
   printf("  DIMENSIONAL VALUES \n");
   printf("          1D resistance: %f\n",r1dstar/dval);
   printf("constriction resistance: %f\n",rcstar/dval);
   printf("       total resistance: %f\n",(rcstar+r1dstar)/dval);
   printf("=================================================");

}