/***************************************************************
 *     gravity1.c: 地球を微小立方体に分割して万有引力を計算する。
 *
 ***************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define G       6.67259e-11
#define RHO     5.525e+3 
#define M       5.974e+24 
#define MM       0.1e+01
#define R	6.367e+6
#define PI	0.31415926e+01 

/**************************************************************/
main(int argc, char *argv[])
{

    int    status, iter;
	double divnum;
	double  	dx, dy, dz;
	double	px, py, pz;
	double	hx, hy, hz, hz2;
	double 	distp1, distp2, r, l, rr, h, costh;
	double 	df, dfsum, dv, v, vv, f1, f2;
	double  fn, fs, dfnsum, dfssum;
	double  ns_rate;
	double 	err, err2;
	double  ncount, scount;
	double lcube, mcube;
	   
	if (argc != 3)
	{
		printf("使用法: $ gravity 高度(Km) 微小立方体１辺の分割数\n");
		printf("使用例: $ gravity 0 268\n");
		exit(-1);
	}

	divnum = (double)atof(argv[2]);
	dfsum = df = f1 = f2 = 0;
	dv = v = vv = 0;
	dfnsum = dfssum = 0;
	ncount = scount = 0;
	iter = 0;

	rr = R;
	h = (float)atof(argv[1])*1000;
	dx = dy = dz = R / divnum;
	px = py = pz = -rr;
	hz = R+h;  
	hz2 = hz*hz;
	hx = hy = 0;  
	lcube = dx/1000;
	iter=0;
	do { /* z loop */
	do { /* y loop */
	do { /* x loop */
 
		distp1 = sqrt(px*px + py*py + pz*pz);
		if (distp1 > rr)
			goto label;
		r = distp1;
		distp1 = px*px + py*py + pz*pz;

		distp2 = (hx-px)*(hx-px)+(hy-py)*(hy-py)+(hz-pz)*(hz-pz);
		l = sqrt(distp2);

		iter++;	
		dv = dx*dy*dz;
		df = (distp2+hz2-distp1)/(2*l*hz)/distp2*dx*dy*dz;
		dfsum += df;
		if (pz > 0) 
		{
			dfnsum += df;
			ncount++;
		} else {
			dfssum += df;
			scount++;
		}

		v += dv;
label:
		px += dx;
	} while( px < rr );
		px = -rr;
		py += dy;
	} while( py < rr );
		px = -rr;
		py = -rr;
		pz += dz;
	} while( pz < rr);

	vv = 4*PI*rr*rr*rr/3;
	f1 = vv*RHO*MM*G/(rr+h)/(rr+h);
	f2 = RHO*MM*G*dfsum;
	vv = 4*PI*rr*rr*rr/3;
	err =  (v/vv-1)*100;
	err2 = (f2/f1-1)*100;
	mcube = M/iter/1000;
	ns_rate = dfnsum/dfsum ;

	printf("----------------------------------------\n");
	printf("観測点の地表からの高度(Km)...: %d\n",(int)h/1000);
	printf("ニュートン方式での引力.......: %f\n",f1);
	printf("微小立方体分割方式での引力...: %f\n",f2); 
	printf("上の2方式での誤差(%).........: %f\n",-err2);
	printf("地球の分割個数...............: %d\n",iter);
	printf("微小立方体の１辺の長さ (Km)..: %f\n",lcube);
	printf("微小立方体１個の質量 (兆トン): %f\n",mcube/1e+12);
	printf("北半球が引力に寄与する割合(%): %f\n",dfnsum/dfsum*100);
	printf("地球体積の数値積分での誤差(%): %f\n",-err);
}