implicit real (a-m,o-z)
      character*8 name(100)
      dimension ncoil(20),ngroup(100),rcoil(100),zcoil(100),turns(100)
      dimension frinf(20,20),fzinf(20,20)
      open(10,file='pfforce.d',status='old')
      open(20,file='frinfmtx.o',status='old')
      open(30,file='fzinfmtx.o',status='old')
      pi=3.141592654
      mu0=4.*pi*1.e-7
      read(10,*)ncoilgroup
      write(6,*)'ncoilgroup=',ncoilgroup
      ncoiltotal=0
      do 100 n=1,ncoilgroup
      read(10,*)ncoil(n),name(n)
      ncoiltotal=ncoiltotal+ncoil(n)
100   continue
      write(6,*)'ncoiltotal=',ncoiltotal
      do 200 n=1,ncoiltotal
      read(10,*)ngroup(n),rcoil(n),zcoil(n),turns(n)
      write(6,*)n,ngroup(n),rcoil(n),zcoil(n),turns(n)
200   continue
      do 300 nsource=1,ncoiltotal
      nsourcegroup=ngroup(nsource)
      do 299 ntarget=1,ncoiltotal
      ntargetgroup=ngroup(ntarget)
      if(ntargetgroup.eq.nsourcegroup)go to 299
      rc=rcoil(nsource)
      zc=zcoil(nsource)
      ic=turns(nsource)
      rp=rcoil(ntarget)
      zp=zcoil(ntarget)
      ip=turns(ntarget)
      call drvgrn(rp,zp,rc,zc,
     & g,gr,gz,grr,grz,gzz)
      bz=-1.*ic*mu0*gr/2./pi/rp
      br=ic*mu0*gz/2./pi/rp
      fr=ip*2.*pi*rp*bz*1e6*0.2248
      fz=-1.0*ip*2.*pi*rp*br*1e6*0.2248
      frinf(ntargetgroup,nsourcegroup)=frinf(ntargetgroup,
     &nsourcegroup)+fr
      fzinf(ntargetgroup,nsourcegroup)=fzinf(ntargetgroup,
     &nsourcegroup)+fz
299   continue
300   continue
599   format(12X,12A8)
      write(20,599)(name(n),n=1,ncoilgroup)
      write(30,599)(name(n),n=1,ncoilgroup)
      do 500 n1=1,ncoilgroup
600   format(A8,12f8.1)
      write(20,600)name(n1),(frinf(n1,n2),n2=1,ncoilgroup)
      write(30,600)name(n1),(fzinf(n1,n2),n2=1,ncoilgroup)
500   continue
      end

      subroutine drvgrn( r1, z1, r2, z2, 
     1     g, gr, gz, grr, grz, gzz )
c-----------------------------------------------------------------------
c     computes:
c     complete elliptic integral of the first kind (ke)
c     complete elliptic integral of the second kind (ee)
c     greens function of the grad shafranov operator (g)
c     first derivatives of greens function wrt r1 (gr) and z1 (gz)
c     second derivative of greens function, wrt r1 (grr),
c     wrt z1 (gzz) and wrt r1 and z1 (grz).
c     uses polynomial approximations to the complete elliptic
c     integrals from c. hastings jr.
c     approximations for digital computers
c     princeton university press
c     
c     m. f. reusch - princeton university plasma physics laboratory
c     9/8/86
c     note: no overflow protection for k = 1 ( r1 = r1, z1 = z2 )
c     if Bvector = Grad phi X Grad psi / ( 2 Pi )
c     where phi is toroidal angle and psi poloidal flux then
c     fields at ( r1, z1 ) due to coil at ( r2, z2 ) with current I
c     are:
c     Bz(r1,z1) = - mu0 I (gr) / ( 2 pi r1 )
c     Br(r1,z1) =   mu0 I (gz) / ( 2 pi r1 )
c     dBr/dr = mu0 I ( (grz) / ( 2 pi r1 ) - (gz) / ( 2 pi r1 r1 ) )
c     dBz/dr =-mu0 I ( (grr) / ( 2 pi r1 ) - (gr) / ( 2 pi r1 r1 ) )
c     dBr/dz = mu0 I ( (gzz) / ( 2 pi r1 ) )
c     dBz/dz =-mu0 I ( (grz) / ( 2 pi r1 ) )
c-----------------------------------------------------------------------
      real k, k2, k3, k4, kp2, kpp2, kp4, kln
c-----------------------------------------------------------------------
      data
     .     ak0/1.38629436112/,ak1/0.09666344259/,ak2/0.03590092383/,
     .     ak3/0.03742563713/,ak4/0.01451196212/,bk0/0.50000000000/,
     .     bk1/0.12498593597/,bk2/0.06880248576/,bk3/0.03328355346/,
     .     bk4/0.00441787012/,
     .     ae0/1.00000000000/,ae1/0.44325141463/,ae2/0.06260601220/,
     .     ae3/0.04757383546/,ae4/0.01736506451/,be0/0.00000000000/,
     .     be1/0.24998368310/,be2/0.09200180037/,be3/0.04069697526/,
     .     be4/0.00526449639/
c-----------------------------------------------------------------------
c     
      z = z1 - z2
      zsq = z * z
      r1sq = r1 * r1
      r2sq = r2 * r2
      r1r2 = r1 * r2
      s = r1 + r2
      s2 = s * s
      a2 = 4.0 * r1r2
      a4 = a2 * a2
      a = sqrt( r1r2 )
      k2 = a2 / ( s2 + zsq )
      k = sqrt(k2)
      k3 = k * k2
      k4 = k2 * k2
      kp2 = 1.0 - k2
      kp4 = kp2 * kp2
      kpp2 = 2.0 - k2
      kln = - alog(kp2)
c     
c     k(qk) complete elliptic integral of first kind
c     
      ek=ak0+kp2*(ak1+kp2*(ak2+kp2*(ak3+kp2*ak4)))
     1     +kln*(bk0+kp2*(bk1+kp2*(bk2+kp2*(bk3+kp2*bk4))))
c     
c     e(qk) complete elliptic integral of second kind
c     
      ee=ae0+kp2*(ae1+kp2*(ae2+kp2*(ae3+kp2*ae4)))
     1     +kln*(be0+kp2*(be1+kp2*(be2+kp2*(be3+kp2*be4))))
c     
c     modulus dependent part of the greens function
c     
      gt =  ( kpp2 * ek - 2.0 * ee ) / k
c     
c     first derivative of modulus dependent part of the greens function
c     
      dgdk = ( kpp2 * ee - 2.0 * kp2 * ek ) / k2 / kp2
c     
c     second derivative of modulus dependent part of the greens function
c     
      d2gdk2 = ( kp2 * ( 5.0 * kp2 - 1.0 ) * ek
     1     - ( k4 - 7.0 * k2 + 4.0 ) * ee ) / k3 / kp4
c     
c     r1 derivative of the modulus
c     
      dkdr = 0.5 * k * ( 1.0 - 0.5 * s * k2 / r2 ) / r1
c     
c     z1 derivative of the modulus
c     
      dkdz = - z * k3 / a2
c     
c     second z1 derivative of the modulus
c     
      d2kdz2 = k3 * ( 3.0 * k2 * zsq - a2 ) / a4
c     
c     second r1 derivative of the modulus
c     
      d2kdr2 = k * ( r2sq * ( 3.0 * k4 - 4.0 * k2 - 4.0 )
     1     + k2 * r1r2 * ( 6.0 * k2 - 8.0 )
     2     + 3.0 * r1sq * k4  ) / a4
c     
c     wrt r1 and z1
c     
      d2kdrz = z * k3 * ( r2 * ( 3.0 * k2 - 2.0 ) + 3.0 * r1 * k2 )
     1     / a4
c     
c     derivatives of a
c     
      dadr = 0.5 * a / r1
      d2adr2 = -0.25 * a / r1sq
c     
c     the greens function
c     
      g = - a * gt
c     
c     r1 derivative of the greens function
c     
      gr = - dadr * gt - a * dkdr * dgdk
c     
c     z1 derivative of the greens function
c     
      gz = - a * dkdz * dgdk
c     
c     second r1 derivative of the greens function
c     
      grr = - d2adr2 * gt - 2.0 * dadr * dkdr * dgdk
     1     - a * ( d2kdr2 * dgdk + dkdr * dkdr * d2gdk2 )
c     
c     second z1 derivative of the greens function
c     
      gzz = - a * ( d2kdz2 * dgdk + dkdz * dkdz * d2gdk2 )
c     
c     wrt r1 and z1
c     
      grz = - a * ( dkdr * dkdz * d2gdk2 + d2kdrz * dgdk )
     1     - dadr * dkdz * dgdk
c     
      return
      end