#include <iostream>
#include <cmath>
using namespace std;

int main() {
	double r{ 0 };
	double forceUnit{ 0 };
	const double sphereX{ 15 };
	const double sphereY{ 15 };
	const double sphereZ{ 15 };
	const double pointX{ 15 };
	const double pointY{ 15 };
	const double pointZ{ 80 };
	const double pointMass{ 10 };
	const int sphereSize{ 31 };
	const double sphereRadius{ 15 };
	double sphereMass[sphereSize][sphereSize][sphereSize];
	double sphereForceX[sphereSize][sphereSize][sphereSize];
	double sphereForceY[sphereSize][sphereSize][sphereSize];
	double sphereForceZ[sphereSize][sphereSize][sphereSize];
	double ForceXTotal{ 0 };
	double ForceYTotal{ 0 };
	double ForceZTotal{ 0 };
	double ForceTotal{ 0 };
	double total{ 0 };
	double distance{ (sphereX - pointX)*(sphereX - pointX) + (sphereY - pointY)*(sphereY - pointY) + (sphereZ - pointZ)*(sphereZ - pointZ) };
	
	for (int i = 0; i < sphereSize; i++) {
		for (int j = 0; j < sphereSize; j++) {
			for (int k = 0; k < sphereSize; k++) {
				(sphereX - (double)i)*(sphereX - (double)i) + (sphereY - (double)j)*(sphereY - (double)j) + (sphereZ - (double)k)*(sphereZ - (double)k) <= (sphereRadius*sphereRadius) ? sphereMass[i][j][k] = 1 : sphereMass[i][j][k] = 0;
				total += sphereMass[i][j][k];
			}
		}
	}
	for (int i = 0; i < sphereSize; i++) {
		for (int j = 0; j < sphereSize; j++) {
			for (int k = 0; k < sphereSize; k++) {
				r = (pointX - (double)i)*(pointX - (double)i) + (pointY - (double)j)*(pointY - (double)j) + (pointZ - (double)k)*(pointZ - (double)k);
				forceUnit = (pointMass * sphereMass[i][j][k] / r) / pow(r,.5);
				sphereForceX[i][j][k] = forceUnit * (pointX - (double)i);
				ForceXTotal += sphereForceX[i][j][k];
				sphereForceY[i][j][k] = forceUnit * (pointY - (double)j);
				ForceYTotal += sphereForceY[i][j][k];
				sphereForceZ[i][j][k] = forceUnit * (pointZ - (double)k);
				ForceZTotal += sphereForceZ[i][j][k];
			}
		}
	}
	ForceTotal = ForceXTotal * ForceXTotal + ForceYTotal * ForceYTotal + ForceZTotal * ForceZTotal;
	cout << "Total Mass of sphere " << total << endl;
	cout << "Total force in the X direction " << ForceXTotal << endl;
	cout << "Total force in the Y direction " << ForceYTotal << endl;
	cout << "Total force in the Z direction " << ForceZTotal << endl;
	cout << "Total combined force " << pow(ForceTotal, .5) << endl;
	cout << "Force if all mass was located at point inside sphere " << total * pointMass / distance << endl;
	system("pause");
}