#include <math.h>

#define PI 3.15149
#define ER .0002
#define ERRVAL .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);

double radius(ratio,p,m_air,w,t,lift,er,errval)
  double ratio;
  double p;
  double m_air;
  double w;
  double t;
  double lift;
  double er;
  double errval;
  {
  double a,b;
  double c1,c2;
  double result;

  c1 = (4.0 / 3.0) * PI * ((m_air - p * HE) / (R * t));
  c2 = 2.0 * PI * w;
  /* 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,m_air,t,w,lift;
  double a,b;

  t = TMP;
  m_air = AIR;

  printf("\n");

  for(ratio = 3.0; ratio <= 4.0; ratio += 0.5)
    {
    printf("Short to long axis ratio -> 1:%1.1lf\n",ratio);
    for(p = 1.0;p  <= 2.0; p += 0.25)
      {
      printf("   Guage Pressure = %2.2lf psi\n",14.7*(p-1.0));
      for(w = 0.010; w <= 0.100; w += .010)
        {
        printf("      Envelope areal density = %1.2lf kg/m2\n",w);
        for(n = 0.0; n < 20; n += 5)
          {
          if(n >= 10.0)
            printf("         lift:  %2.1lf       %2.1lf       %2.1lf       %2.1lf       %2.1lf  kg\n",n+1,n+2,n+3,n+4,n+5);
           else
            printf("         lift:   %2.1lf        %2.1lf        %2.1lf        %2.1lf        %2.1lf  kg\n",n+1,n+2,n+3,n+4,n+5);
          printf("             ");
          for(lift = 1.0; lift <= 5.0; lift += 1.0)
            {
            a = radius(ratio,p,m_air,w,t,n+lift,ER,ERRVAL);
            b = ratio * a;
            printf("%1.2lf,%1.2lf  ",a,b);
            }
          printf("\n");
          }
        }
      }
    }
  }