/* Balloon design function by Rick Wagner.                   */
/* Copyright 1999, all rights reserved.                      */
/*                                                           */
/* Produced for USC Computer Science Department class cs101. */
/* No use without attribution.                               */

# include <stdio.h>        /* For scanf() and printf()       */
# include <math.h>         /* For pow() and fabs()           */

/* Function prototype: */
float sphericalDiameter(float p);

/* Compute the diameter of a spherical balloon based on a desired payload */
float sphericalDiameter(float p)
{
  float sfGuess = 0;                               /* Radius guess */
  float sfError = 0;                               /* Error term */
  float sfMaxError = 0;                            /* Maximum error allowed based on desired payload */
  float sfPi = (float) 3.14159;                    /* Single precision value for pi */
  float sfLift = 0;                                /* Lift term */
  float sfWeight = 0;                              /* Weight term */
  int i = 0;                                       /* Loop index */
  float sfTempP = 0;                               /* Temporary payload value used in the loop */

  /* Compute the first guess */
  sfGuess = (p * ((float) 3.0)) / (sfPi * ((float) 4.0) * ((float) 0.053));
  sfGuess = ((float) pow(sfGuess, 1.0 / 3.0));

  /* Compute the maximal error */
  sfMaxError = p / ((float) 10000);                /* One part in ten thousand */

  /* Loop until the error is small enough */
  for (i = 1; i <= 20; i++)                        /* If it doesn't converge in 20 iterations */
  {                                                /* it never will.                          */
    /* Compute payload based on the guess */
    sfLift = ((float) 4.0) * sfPi * sfGuess * sfGuess * sfGuess * ((float) 0.053) / ((float) 3.0);
    sfWeight = ((float) 4.0) * sfPi * sfGuess * sfGuess * ((float) 0.000333) * ((float) 44.0);
    sfTempP = sfLift - sfWeight;
    
    /* Compute the error */
    sfError = p - sfTempP;
    sfError = ((float) fabs(((double) sfError)));
    if (sfError < sfMaxError)
    {
      break;                                        /* Break out of the loop */
    }

    /* Compute the new guess */
    sfGuess = (sfGuess * p) / ((sfTempP + p) / ((float) 2.0));   /* 50% damping factor */

  } /* End of for loop */

  if (i > 20)
  {
    puts("Did not converge, returning zero.\n");    /* Maybe comment this line out */
    return 0;
  }
  else
  {
    printf("Converged in %d iterations.\n", i); /* Maybe comment this line out */
    sfGuess = sfGuess * ((float) 2.0);              /* Double the radius to get the diameter */
    return sfGuess;
  }
}