#include <stdio.h>
#include <math.h>


#define k      9.0e+9   /* N.m2/C2 */
#define Lambda 2.0e-10  /* C/m     */

float L,Ex,Ey,x,y,xp;

float f1(float xp) { return (k*Lambda*(x-xp)/pow((pow(x-xp,2)+pow(y,2)),1.5)); }
float f2(float xp) { return (k*Lambda*( y  )/pow((pow(x-xp,2)+pow(y,2)),1.5)); }
float integrate_Ex(float a,float b);
float integrate_Ey(float a,float b);

void main(void)
{

  L = 2.0;  /* Length of the wire in meter */

  /* calculate electric filed at grid points */

  for(x=-1.0;x<=1.0;x+=0.5)  
  for(y=-1.0;y<=1.0;y+=0.5)
  {
     if(y==0.0) continue;
     Ex = integrate_Ex(-L/2,L/2);
     Ey = integrate_Ey(-L/2,L/2);
     printf("%f\t%f\t%f\t%f\n",x,y,Ex,Ey);
  }

} /* main */



/* Returns integral of a function f1(x) between [a,b] by Simpson Method. */
float integrate_Ex(float a,float b)
{
  float h,sum,x,integral;
  int   i,n;

  n   = 200;     /* number of part */
  h   = (b-a)/n;
  sum = f1(a)+f1(b);

  for(i=1;i<n;i++){
    x = a + i*h;
    if( i%2 ) sum += 4.0*f1(x);
    else      sum += 2.0*f1(x);
  }
  integral = sum*h/3.0;

  return integral;
}



/* Returns integral of a function f2(x) between [a,b] by Simpson Method. */
float integrate_Ey(float a,float b)
{
  float h,sum,x,integral;
  int   i,n;

  n   = 200;     /* number of part */
  h   = (b-a)/n;
  sum = f2(a)+f2(b);

  for(i=1;i<n;i++){
    x = a + i*h;
    if( i%2 ) sum += 4.0*f2(x);
    else      sum += 2.0*f2(x);
  }
  integral = sum*h/3.0;

  return integral;
}