/*
   gravity.cpp
   CIS 150
   7/3/08
   David Klick
   
   This program calculates the effect of gravity on
   a particle within a hollow sphere that is due
   solely to the presence of the sphere and the particle
   (external gravitational fields are ignored).
   
   This problem could be more easily done by creating a
   formula for calculating the gravity and performing 
   integration on it, but this program presents a way of
   doing the integration manually one piece at a time.
   
   A few assumptions are made to simplify the calculations:
   
   1. An imaginary straight line is drawn through the center
      of the sphere which includes the particle's position.
   
   2. It can be easily shown due to symmetries that the
      gravitational forces that are not parallel to that
      line balance out and have no net effect.
      
   3. This means that the only gravitational forces that
      must be calculated are those parallel to the imaginary
      line through the center that the particle is on.
      
   4. The next step is to divide the sphere up into a
      whole bunch of slices, calculate the effect of gravity
      on the particle due to each slice, and add all those
      effects together to see what the net effect is.
      
   5. The method chosen to slice up the sphere is to choose
      equal sized angles measured from the imaginary diameter
      line. The angles are chosen by first deciding how many
      slices are desired and then dividing 180 degrees into
      that many equal divisions.
      
   6. Each slice will include a donut-shaped piece of the
      sphere perpendicular to the imaginary diameter line
      whose center also lies on that line.
      
   7. A rough estimate of the mass of the slice and its
      gravitational effect on the particle can now be 
      calculated. The estimates, and thus the overall
      net effect calculation, should become more accurate
      as the number of slices increases.
      
   8. The program does the above process many times with
      increasing numbers of slices. The answer should be
      seen to be converging on a number. With debugging
      mode turned off, the output clearly shows the result
      converges on 0, meanin that there is no net
      gravitational effect from a uniform hollow sphere on
      a particle contained within it.
*/

#include <iostream>
#include <iomanip>
#include <cmath>
using std::cin;
using std::cout;
using std::setw;
using std::setprecision;
using std::showpoint;
using std::fixed;

const bool DEBUGGING = false;   // set to true to show calculations
const double PI = 3.1415926536;

double convert_r2d(double radians); // converts radians to degrees
double convert_d2r(double degrees); // converts degrees to radians
void cont();

int main() {
   // The following group of variables may be set to
   // other values, but that should not change the
   // basic outcome of this program, which is looking
   // for a net gravitational force of either 0, a pull
   // toward the closer end of the sphere (>0), or a pull
   // toward the farther end of the sphere (<0).
   double sphereRadius = 1.0;          // arbitrary
   double sphereThickness = 1.0;       // arbitrary
   double massPerUnit = 1.0;           // arbitrary
   double massOfParticle = 1.0;        // arbitrary
   double gravitationalConstant = 1.0; // arbitrary

   int steps;             // the number of slices to make
   double angleIncrement; // the angle size for each slice
   double angle;          // current deflection from diameter
    
   double sliceRadius;    // radius of perpendicular slice (donut)
   double sliceWidth;     // width of slice (used for mass calc)
   double sliceMass;      // calculated mass of entire slice
   double gravitationalForce;  // total g force on particle due to slice
   double xComponent;     // g force component along axis of interest
   double particlePosition = 0.5; // range: 0 to sphereRadius
   double xLength;        // distance along diameter from particle
                          // to center of slice
   double particleAngle;  // angle between diameter line and point on
                          // sphere surface where slice intersects,
                          // on the diameter line where the particle
                          // is located
   double distanceToParticle; // distance from shell to particle for
                              //    a particular slice
   double sum;            // accumulated horizontal force due to gravity

   for (steps=10; steps<=100000; steps*=2) { 
      // clear accumulated horizontal g force
      sum = 0.0;
      
      // calculate angle from center of sphere that defines
      //    each slice
      angleIncrement = PI / steps;

      // set initial angle (actually one previous)
      angle = 0 - angleIncrement/2.0;
      // calculate width of each slice (for mass calculation)
      sliceWidth = 2.0 * sphereRadius * sin(angleIncrement/2.0);
      
      if (DEBUGGING) { // display report headings
         cout << " angle  s radius    s mass      base"
              << "   p angle    p dist   g force    x comp\n";
      }
      
      // Iterate over a semicircle starting with 0 degrees
      //    and going counter-clockwise to 180 degrees
      //    (C++ wants angles in radians)
      for (int i=0; i<steps; i++) {
         // angle goes through middle of slice
         angle += angleIncrement;
         
         // radius of slice (perpendicular to diameter line)
         sliceRadius = sin(angle) * sphereRadius;
         sliceMass = 2.0 * PI * sliceRadius * sliceWidth
            * sphereThickness * massPerUnit;
            
         // calculate where evrything is in relation to particle
         xLength = sphereRadius * cos(angle) - particlePosition;
         particleAngle = atan(sliceRadius/xLength);
         if (particleAngle < 0.0) particleAngle += PI;
         distanceToParticle = pow(xLength * xLength + sliceRadius * sliceRadius, 0.5);
         
         // calculate total gravitational effect on particle
         gravitationalForce = gravitationalConstant
            * sliceMass * massOfParticle
            / pow(distanceToParticle, 2.0);
   
         // calculate the g force parallel to the diameter line
         xComponent = gravitationalForce * cos(particleAngle);
         sum += xComponent;  // accumulate total effective force
         
         if (DEBUGGING) {  // display calcs if debugging
            cout << fixed << showpoint << setprecision(1)
                 << setw(6) << convert_r2d(angle) << setprecision(5)
                 << setw(10) << sliceRadius
                 << setw(10) << sliceMass
                 << setw(10) << xLength
                 << setw(10) << convert_r2d(particleAngle)
                 << setw(10) << distanceToParticle
                 << setw(10) << gravitationalForce
                 << setw(10) << xComponent << '\n';
         }
      }
      
      // display results for this specific number of steps
      cout << "Steps: " << setw(6) << steps << fixed
           << ", Gravitational force = " << setprecision(9)
           << sum << '\n';
   }
   
   cont();
   return 0;
}

/*
   double convert_r2d(double radians);
   
   Arguments:
      radians: the radians value to be converted to degrees
   
   Returns: A double which is the number of degrees the input
      angle is equivalent to.
   
   This is a convenience function to make converting from radians
   to degrees easier.
*/
double convert_r2d(double radians) {
   return radians * 360.0 / (2.0 * PI);
}

/*
   double convert_d2r(double degrees);
   
   Arguments:
      degrees: the degrees value to be converted to radians
   
   Returns: A double which is the number of radians the input
      angle is equivalent to.
   
   This is a convenience function to make converting from degrees
   to radians easier.
*/
double convert_d2r(double degrees) {
   return 2.0 * PI * degrees / 360.0;
}

/*
   void cont();
   
   Arguments: none
   
   Returns: nothing
   
   Clears input buffer and pauses waiting for user 
   to press Enter.
*/
void cont() {
   if (cin.rdbuf()->sungetc() != -1 && cin.get() != '\n') cin.ignore(80,'\n');
   cout << "Press enter to continue...";
   cin.get();
}

/*
   Sample output:

   Steps:     10, Gravitational force = 0.100425978
   Steps:     20, Gravitational force = 0.023424269
   Steps:     40, Gravitational force = 0.005769611
   Steps:     80, Gravitational force = 0.001437195
   Steps:    160, Gravitational force = 0.000358976
   Steps:    320, Gravitational force = 0.000089724
   Steps:    640, Gravitational force = 0.000022430
   Steps:   1280, Gravitational force = 0.000005607
   Steps:   2560, Gravitational force = 0.000001402
   Steps:   5120, Gravitational force = 0.000000350
   Steps:  10240, Gravitational force = 0.000000088
   Steps:  20480, Gravitational force = 0.000000022
   Steps:  40960, Gravitational force = 0.000000005
   Steps:  81920, Gravitational force = 0.000000001
*/