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

#define PI 3.15149
#define ER 0.0002
#define ERRVAL 0.0001
#define R 0.0000821
#define HE 0.004
#define AIR 0.029
#define TMP 298.0

double radius(double,double,double,double,double,double,double);

double radius(ratio,p_air,p_he,w,t_air,t_he,lift)
  double ratio;
  double p_air;
  double p_he;
  double w;
  double t_air;
  double t_he;
  double lift;
  {
  double a,b;
  double c1,c2;
  double result;
  double d_air,d_he;
  double er;

  d_air = AIR * p_air / (R * t_air);
  d_he = HE * p_he / (R * t_he);
  c1 = (4.0 / 3.0) * PI * (d_air - d_he);
  c2 = 2.0 * PI * w;

  er = ER;

  /* note: a must be > 0 to be valid. */

  a = 0.0;

  do
    {
    b = ratio * a;
    if((result = (c1 * a * a * b - c2 * a * b - lift)) < -ERRVAL)
      {
      er *= 1.5;
      a += er;
      }
     else
       if(result > ERRVAL)
         {
         er /= 2;
         a -= er;
         if(a < 0.0)
           a = -a;
         }
    } while(!(result > -ERRVAL && result < ERRVAL));

  return(a);
  }

main()
  {
  double n,ratio,p_he,p_air,t_air,t_he,w,lift;
  double a,b;


  p_air = 1.0;
  t_air = TMP;
  ratio = 3.0;
  lift  = 5.0;

  printf("\nRatio -> 1:%1.1lf, lift = %2.1lf kg, temp of air = %2.1lf C\n\n",ratio,lift,t_air - 273.0);

  for(t_he = TMP; t_he <= TMP + 30.0; t_he += 5.0)
    {
    printf("Temperature of He = %1.1lf C\n",t_he - 273.0);
    for(p_he = 1.0;p_he  <= 4.0; p_he += 0.25)
      {
      printf("   Guage Pressure of He = %2.2lf psi\n",14.7 * (p_he - 1.0));
        for(n = 0.0; n < 0.1; n += .05)
          {
          printf("      areal d:  %1.2lf       %1.2lf       %1.2lf       %1.2lf       %1.2lf  kg/m2\n",n+0.01,n+0.02,n+0.03,n+0.04,n+0.05);
          printf("             ");
          for(w = 0.01; w <= 0.05; w += 0.01)
            {
            a = radius(ratio,p_air,p_he,w+n,t_air,t_he,lift);
            b = ratio * a;
            printf("%1.2lf,%1.2lf  ",a,b);
            }
          printf("\n");
          }
      }
    }
  }