*** Copyright (c) Giulio Stancari
*** This file can be copied and used, but it cannot be modified
*** without consent of the author.


*** This macro defines the function FORCE(R,Z), which is the force between two
*** charged conducting spheres (of equal charge and radius a), divided by the
*** force between two pointlike objects at the same distance R and of the same
*** charge.

*** FORCE = (force between two spheres) / (force between pointlike charges)
*** R = (distance between centers of spheres) / (radius of spheres)
*** Z = order of approximation, Z <= 100

*** The following vectors are also filled. They are useful for judging the
*** quality of the numerical approximation:

*** Q(i) = i-th image charge
*** X(i) = distance of the i-th image charge
           from center of sphere (units of radius)
*** F(i) = force calculated at the i-th iteration
*** E(i) = estimated error on F(i)

*** Example of usage:
*** PAW>exec supercoulomb.kumac
*** PAW>function/plot force(x,10.) 2 10
*** PAW>function/plot force(3.,x) 3 30

*** For a description of the algorithm, see "The Feynman Lectures on Physics",
*** vol. II, par. 6-9.


APPLICATION COMIS QUIT

  REAL FUNCTION FORCE(R,Z)

  REAL R,Z

  VECTOR Q(100),X(100),F(100),E(100)
  VECTOR CHARGE(1)
  VECTOR F0(1)
  INTEGER I,J

  IF(R.LT.2.OR.Z.LT.2.OR.Z.GT.100)THEN
    FORCE=0.
  ELSE

    Q(1) = 1.
    X(1) = 0.

    CHARGE(1) = Q(1)
    FORCE = Q(1)*Q(1)/R/R

    DO I=2,INT(Z)
      Q(I) = Q(I-1) / ( R - X(I-1) )
      X(I) = 1. / ( R - X(I-1) )
      CHARGE(1) = CHARGE(1) + Q(I)
      F0(1) = CHARGE(1)*CHARGE(1)/R/R
      DO J=1,I
        FORCE = FORCE + Q(I)*Q(J)/(R-X(I)-X(J))**2.
      ENDDO
      F(I) = FORCE/F0(1)
      E(I) = F(I) - F(I-1)
    ENDDO

    FORCE = FORCE/F0(1)

  ENDIF

  RETURN
  END

QUIT