EQUIPOTENTIALS OF A DIPOLE

www.geocities.co/serienumerica
www.oocities.org/reibaretti2004
e-mail reibaretti2004@yahoo.com

c Electric field lines and equipotentials
cEquipotentials /of an electric dipole
c  reference John R. Merrill , Using Computers in Physics, University
c  Press of America , Chapter 4
      dimension x(100), y(100),q(100) ,xe(1000),ye(1000),xv(1000),
     $ yv(1000)
      real k
      k=8.99E9
      pi=2.*asin(1.)
c postion of dipole charges
      x(1)=0.
      y(1)=0.
      x(2)=1.
      y(2)=0.
      n=2
c initial point
      xp=0.1
      yp=0.5
c    charges
      q(1)=1.e-9
      q(2)=-1.e-9
      xe(1)=xp
      ye(1)=yp
      xv(1)=xp
      yv(1)=yp
      deltas=5.*abs(yp)/500.
c dimension of xe , xv etc. must be greater than niter
      niter=int(1./deltas)
      print*,'niter=',niter
      kprint=int(float(niter)/20.)
      Kount=kprint
c niter =no. of points to be plotted
      do 70 it=1,niter
      ex=0.
      ey=0.
      v=0.
      do 20 i=1,n
      rpi=sqrt((xe(it)-x(i))**2 +(ye(it)-y(i))**2)
      if(rpi.le.1.e-1)goto 200
      rvi= sqrt((xv(it)-x(i))**2 +(yv(it)-y(i))**2)
      ex=ex+k*q(i)*(xe(it)-x(i))/rpi**3
      ey=ey+k*q(i)*(ye(it)-y(i))/rpi**3
      v=v+k*q(i)/rvi
20    continue
      er=(ex**2+ey**2)**.5
      theta= atan(ey/ex)
      if(it.eq.1)then
      print 110,V,xv(it),yv(it)
c      print 100, er,xe(it),ye(it),theta*180./pi
      endif
      if(it.eq.kount)then
      print 110,V,xv(it),yv(it)
c      print 100, er,xe(it),ye(it),theta*180./pi
      kount=kount+kprint
      print*,'  '
      endif
c must calculate either equipot or  e-lines but  not both // (dxe,dye)
c will follow a E-line (E is not constant )// (dxv,dyv)
c will follow an equipotential
      dxe=deltas*(ex/er)
      dye=deltas*(ey/er)
      dxv=dye
      dyv= -dxe
c coordinates to follow an E-line
c      xe(it+1)=xe(it)+dxe
c      ye(it+1)=ye(it)+dye
c      xv(it+1)=xv(it)+dxe
c      yv(it+1)=yv(it)+dye
c coordinates to follow an equipotential
      xe(it+1)=xe(it)+dxv
      ye(it+1)=ye(it)+dyv
      xv(it+1)=xv(it)+dxv
      yv(it+1)=yv(it)+dyv
70    continue
100   format(1x,'Er,xe,ye, theta(degrees)=',4(2x,e11.4))
110   format(1x,'V,x,y=',3(2x,e11.4))
      print*,'  '
200   stop
      end