/*
**	mom.c
**	caluculate moment of a charge beside a current
**
*/
#include <stdio.h>
#include <math.h>
/*
the postion of charge: (1, 0)
*/
double xc=  1.0;
double yc=  0.0;
double zc=  0.0;

struct vector{
double	x;
double	y;
double	z;
	};
typedef struct vector Vec;

Vec E(x,y,z)
double x,y,z;
	{
	Vec e;
	double r;
	x-= xc;
	/*
	y-= yc;
	z-= zc;
	*/

	r= x*x+y*y+z*z;
	r*= sqrt(r);
	e.x= x/r;
	e.y= y/r;
	e.z= z/r;
	return e;
	}

Vec B(x,y,z)
double x,y,z;
	{
	Vec b;
	double r= x*x+y*y+z*z;
	b.x= y/r;
	b.y= -x/r;
	b.z= 0;
	return b;
	}
double
PointVz(x,y,z)
double x,y,z;
	{
	Vec e;
	Vec b;
	e= E(x,y,z);
	b= B(x,y,z);
	return e.x*b.y-e.y*b.x;	/* S= E x B */
	}

main()
	{
	double x,y,z;
	double t0,t1,t2,t3,td;
	double S,S1,S2,S3;

	t0= 10000;
	t1= 100;
	t2= 10;
	t3= 2;
	td= 0.01;

	S= S1= S2= S3= 0.0;

	for(z= 0; z<t0 ; z+=100.)
	for(y= 0; y<t0 ; y+=100.)
	for(x= -t0; x<t0 ; x+=100.)
	if(x<=-t1||x>=t1||y>=t1||z>=t1)
		S1+= PointVz(x,y,z);

	S1*= 100*100*100;
	printf("1st step done S1= %f\n",S1);

	for(z= 0; z<t1 ; z+=1.)
	for(y= 0; y<t1 ; y+=1.)
	for(x= -t1; x<t1 ; x+=1.)
	if(x<=-t2||x>=t2||y>=t2||z>=t2)
		S2+= PointVz(x,y,z);

	S2*= 1*1*1;
	printf("2nd step done S2= %f\n",S2);

	for(z= 0; z<t2 ; z+=.1)
	for(y= 0; y<t2 ; y+=.1)
	for(x= -t2; x<t2 ; x+=.1)
	if(x<=-t3||x>=t3||y>=t3||z>=t3)
		S3+= PointVz(x,y,z);
	S3*= .1*.1*.1;
	printf("3rd step done S3= %f\n",S3);

	for(z= td/2.; z<t3 ; z+= td)
	for(y= td/2.; y<t3 ; y+= td)
	for(x= -t3-td/2.; x<t3 ; x+= td)
		S+= PointVz(x,y,z);
	S*= td*td*td;
	printf("4th step done S= %f\n",S);

	S+= S1+S2+S3;
	printf("sum of S is %f\n",S);
	S*= 4;
	printf("S is %f for all space\n",S);
	}