PROGRAM Electric_Field
IMPLICIT NONE

REAL, PARAMETER :: k = 9.0E+9      ! N.m2/C2
REAL, PARAMETER :: Lambda = 1.0E-9 ! C/m
REAL            :: L,Ex,Ey,x,y,xp,f1,f2

f1(xp) = k*Lambda*(x-xp)/( (x-xp)**2+y**2 )**1.5
f2(xp) = k*Lambda*( y  )/( (x-xp)**2+y**2 )**1.5


L = 2.0            ! Length of the wire in meter

DO x=-1.0,1.0,0.5  ! calculate electric filed at grid points
DO y=-1.0,1.0,0.5

  IF(y==0.0) CYCLE

  Ex = Integrate_Ex(-L/2,L/2)
  Ey = Integrate_Ey(-L/2,L/2)
  PRINT *,x,y,Ex,Ey

END DO
END DO


CONTAINS

REAL FUNCTION Integrate_Ex(A,B)
!------------------------------------------------------
! Returns integral of a function F1(x) between
! [a,b] by Simpson Method.
!------------------------------------------------------
REAL, INTENT(IN) :: A,B
REAL             :: h,Sum,x
INTEGER          :: I,N

 N   = 200     ! number of part
 h   = (B-A)/N
 Sum = F1(a)+F1(b)

 DO I=1,N-1
    x = a + i*h
    IF(MOD(I,2)==1) THEN
        Sum = Sum + 4.0*F1(x)
     ELSE
        Sum = Sum + 2.0*F1(x)
     END IF
 END DO

 Integrate_Ex = Sum*h/3.0

END FUNCTION Integrate_Ex

REAL FUNCTION Integrate_Ey(A,B)
!------------------------------------------------------
! Returns integral of a function F2(x) between
! [a,b] by Simpson Method.
!------------------------------------------------------
REAL, INTENT(IN) :: A,B
REAL             :: h,Sum,x
INTEGER          :: I,N

 N   = 200     ! number of part
 h   = (B-A)/N
 Sum = F2(a)+F2(b)

 DO I=1,N-1
    x = a + i*h
    IF(MOD(I,2)==1) THEN
        Sum = Sum + 4.0*F2(x)
     ELSE
        Sum = Sum + 2.0*F2(x)
     END IF
 END DO

 Integrate_Ey = Sum*h/3.0

END FUNCTION Integrate_Ey

END PROGRAM Electric_Field