c*******************************
c Solve for charge distribution on a rod using "classical" MoM
c C.Furse 5-9-95
c /grad/furse/ee553.dir/mom2.f
c*******************************
      parameter(Nseg = 10 ) ! # of segments on rod
      real amatrix(Nseg,Nseg),b(Nseg)
      integer ipvt(Nseg) ! for linpack
      real L ! length of rod (meters)
      common dx

c     --Define problem parameters--
      L = 1.0 ! length of rod (meters)
      dx = L/Nseg ! size of each segment (dx)
      radius = 0.001 ! radius of wire (meters)

c     --Set up matrix--
      do 10 m=1,Nseg
      do 10 n=1,Nseg
         if (m.eq.n) then ! self-term
            amatrix(m,n) = 2.*log(dx/radius)
     &  * ww(m)
         else
           amatrix(m,n) = wg(m,n)
         endif
10    continue

c     --Set up b vector--
      pi = acos(-1.)
      eo = 8.854e-12
      V = 1. ! volts
      do 20 m=1,Nseg
         b(m) = 4.*pi*eo*V*ww(m)
20    continue

      lda = Nseg
      job = 0
c     --Linpack subroutine sgefa computes LU decomposition
c     call sgefa(matrix,lda,Nseg,ipvt,info)
c     --Linpack subroutine sgesl solves Ax=b
c     call sgesl(matrix,lda,Nseg,ipvt,b,job)

c     --print out for matrix solver--
      open(60,file='matrix.in')
      write(60,*) Nseg,Nseg
      do i=1,Nseg
      do j=1,Nseg
         write(60,*) amatrix(i,j)
      enddo
      enddo
      do i=1,Nseg
         write(60,*) b(i)
      enddo
      close(60)

c     --Print out answer--
      do n=1,Nseg
      enddo

      end
c*******************************
c Define basis function: un(x')
c*******************************
      function u(n,x)
      common dx

      xn = (float(n)-0.5)*dx

c     --Rectangular pulse--
      if ((x.ge.(xn-dx/2.)).and.(x.le.(xn+dx/2.))) then
         basis_fn = 1.
      else
         basis_fn = 0.
      endif
c
c     --triangular basis fuction--
c      (remember to change limits of an integration)
c     if ((x.ge.(xn-dx)).and.(x.le.xn)) then
c        basis_fn = (x-xn-dx)/dx
c     elseif ((x.ge.xn).and.(x.le.(xn+dx)) then
c        basis_fn = -(x-xn-dx)/dx
c     else
c        basis_fn = 0.
c     endif

	u = basis_fn

      return
      end
c*******************************
c Define g(x,x')
c*******************************
      function g(x,xp)

      if (x.ne.xp) then
      g_fn = 1./abs(x-xp)
      else  ! take care of occasional overlaps
      g_fn = 100.
      endif

	g = g_fn

      return
      end
c*******************************
c Define weighting function: wm(x)
c*******************************
      function w(m,x)
      common dx

      xm = (float(m)-0.5)*dx

c     --Rectangular pulse--
      if ((x.ge.(xm-dx/2.)).and.(x.le.(xm+dx/2.))) then
         wm = 1.
      else
         wm = 0.
      endif
 
c     --triangular basis fwmction--
c     if ((x.ge.(xm-dx)).and.(x.le.xm)) then
c        wm = (x-xm-dx)/dx
c     elseif ((x.ge.xm).and.(x.le.(xm+dx)) then
c        wm = -(x-xm-dx)/dx
c     else
c        wm = 0.
c     endif

	w = wm
      return
      end
c******************************
c Use trapezoidal integration to find gn(x)
c*******************************
      function gn(n,x)
      common dx

      xn = (float(n)-0.5)*dx

      Nn = 10 ! # of integration points
      an = xn - dx/2. ! limits of n integration
      bn = xn + dx/2.
      hn = (bn-an)/Nn ! resolution of integration

      gnx =  u(n,an)*g(x,an)/2. + u(n,bn)*g(x,bn)/2.
      do i=1,Nn-1
         gnx = gnx + u(n,an+i*hn)*g(x,an+i*hn)
      enddo
      gnx = gnx*hn

      gn = gnx
      return
      end
c*********************************
c Use trapezoidal integration to find int wm(x)gn(x)
c*********************************
      function wg(m,n)
      common dx
      
      xm = (float(m)-0.5)*dx

      am = xm - dx/2. ! limits of m integration
      bm = xm + dx/2. 
      Nm = 10 ! # of points in integration
      hm = (bm-am)/Nm

c     --Use trapezoidal integration to find int wm * gn--
c     (non-delta weighting function)
      wmgn = w(m,am)*gn(n,am)/2. + w(m,bm)*gn(n,bm)/2.
      do j=1,Nm-1
         wmgn = wmgn + w(m,am+j*hm)*gn(n,am+j*hm)
      enddo
      wmgn = wmgn* hm

c     --Using delta pulse weighting function: int wm*gn = gn(xm)
c     wmgn = gn(n,xm)

      wg = wmgn
      return
      end
c*********************************
c Use trapezoidal integration to find int wm(x)
c*********************************
      function ww(m)
      common dx
      
      xm = (float(m)-0.5)*dx

      am = xm - dx ! limits of m integration
      bm = xm + dx 
      Nm = 10 ! # of points in integration
      hm = (bm-am)/Nm

c     --Use trapezoidal integration to find int wm*f
c     (non-delta weighting function)
      wm = w(m,am)/2. + w(m,bm)/2.
      do j=1,Nm-1
         wm = wm + w(m,am+j*hm)
      enddo
      wm = wm* hm

c     --Assuming delta weight function--
c     wm = 1.

      ww = wm
      return
      end