!  nldp4.f90 - stellarator symmetry added to nldp2.f90 version
 !              (at present, forces stellarator symmetry)
module mydata
  real*8, dimension(:,:), allocatable :: A    ! coefficient matrix Ax=b
  real*8, dimension(:)  , allocatable :: bvac    ! vacuum field
  real*8, dimension(:,:)  , allocatable :: h    ! total field intensity, A/m
  real*8, dimension(:,:), allocatable :: r    ! location of spheres/lumped masses
  real*8, dimension(:)  , allocatable :: xmur ! relative permeabilty
  real*8, dimension(:)  , allocatable :: vol  ! volume of spheres/lumped masses
  real*8, dimension(:,:), allocatable :: rp   ! remote points
  real*8, dimension(:,:), allocatable :: bp   ! field at remote points
  real*8, dimension(:,:), allocatable :: xmag ! magnetization of spheres/lumped masses
  real*8 :: r1, r2, r3, pi, xmu0, rij(3), fac, xmag_net, bp_net, db, btot(3), bpmax
  real*8 :: xp, yp, zp, bx, by, bz            ! point and field for coilfld
  real*8 :: bsat, hsat, xmur_h0
  real*8 :: ang, co, si, Aij, rj(3), ri(3), alp, xmagi(3)
  !
  integer :: is, iper, nper
  integer :: n, n3, i, ii, j, jj, ia, ja, ib, np, ip, iopt
end module mydata
!
!****************************************************************************
!
program multidipole
  !
  ! routine to calculate magnetization of distributed spheres of constant permeable material
  ! in a fixed vacuum background field.
  !
  ! Solves the following system of equations for spheres:
  !  Mi = 3/xmu0*(xmur-1)/(xmur+2)*{Bv + sum(Mj*3[(ud.ur)ur-ud]/rij**3*Volj*1.e-7)}
  !
  ! which is based on M = 3/xmu0*(xmur-1)/(xmur+2)*B from Jackson's "Classical Electrodynamics", eq 5.115
  !  for "soft" magnetic material (no hysteresis, ie M=0 at B=0) in an external field
  !
  !  The above is recast as Ax-b=0.
  !  where b = Bv  (Vacuum field ordered Bx(1),By(1),Bz(1),Bx(2),By(2),Bz(2),...
  !        x = M   (Magnetizations to solve for and ordered Mx(1),My(1),Mz(1),Mx(2),My(2),Mz(2),...
  !        A = 1/(3/xmu0*(xmur-1)/(xmur+2)) for i=j
  !          = -3[(ud.ur)ur-ud]/rij**3*Volj*1.e-7 for i<>j
  !            (field at dipole i do to unit dipole at j)
  !
  ! where the background field at sphere i is the Vacuum Field plus the sum of the field
  !  contributions from all other spheres (dipoles)
 !
 ! Optionally assumes all spheres independent (no field contribution from other spheres)
 !  and solves Mi = 3/xmu0*(xmur-1)/(xmur+2)*Bv
  ! iopt = 0  - linear (no coupled field, linear b-h curve ie constant mur)
  ! iopt = 1  - do full nonlinear (coupled field, nonlinear b-h curve)
  ! iopt = 2  - no field coupling, but nonlinear b-h curve with Bsat
  ! iopt = 3  - field coupling, but linear b-h curve (ie constant mur)
 !
  !
!  use nag_gen_lin_sys
  use nag_nlin_sys
  use mydata
  implicit none
  real*8, dimension(:)  , allocatable :: x    ! vacuum field / magnetization
  logical ltotal ! flag to calculate total field (can take a bit of time for all plasma points)
  real*8 :: scale

  external fun

  character*80 :: coilfile
  !
  pi = asin(1.0d0)*2.0d0
  nper = 3
  alp = pi*2.d0/nper
  xmu0 = 4.d-7*pi

  read(*,*) coilfile
    write(*,*) coilfile, "/ coilfile"
  read(*,*) n, xmur_h0, bsat, ltotal, scale, iopt
    write(*,*) n, xmur_h0, bsat , ltotal,  scale, "/ n, xmur_h0, bsat, ltotal, scale"
  hsat = 3.d0*bsat/xmu0/(xmur_h0-1.d0) ! assuming exponential form for saturuation
  allocate(vol(n),r(3,n),xmur(n))

  n3 = 3*n
  allocate(A(n3,n3), x(n3), bvac(n3), h(3,n), xmag(3,n))

  do i = 1,n
    read(*,*) r(1:3,i), vol(i)   !  , xmur(i)
    write(*,"(1p5d15.6)") r(1:3,i), vol(i) ! , xmur(i)
  enddo
!  vol = vol*1.111 ! TMP FIX
      r = r*scale
      vol = vol*scale**3

!  call filla

!  do ja=1,n3
!  write(*,"(1p12e10.2)") A(1:n3,ja)
!  enddo

  open(unit=7, file=coilfile, status="old")
  call readfila(7)

  call fillb

!  write(*,"(1p3e15.6)") bvac

  xmag=reshape(bvac,(/3,n/))

  write(*,*) "vacuum field at magnetized material"
  write(*,"(a5,7a15)")"i","x","y","z","Bvacx","Bvacy","Bvacz","Bvacnet"
  do i=1,n
  rij(:) = xmag(:,i)
  xmag_net = sqrt(dot_product(rij,rij))
  write(*,"(i5,2x,1p7e15.6,0p,f10.5)") i,r(1:3,i),xmag(1:3,i), xmag_net
  enddo

  !
  if(iopt.gt.0) then
  ! solve for magnetization m (not dipole m = m*vol)

!  call nag_gen_lin_sol(A,bvac)
     x = 0.
     call nag_nlin_sys_sol(fun,x, user_jac=.false.)

     xmag=reshape(x,(/3,n/))
  else
    do i=1,n
      do ii=1,3
        ia = ii+3*(i-1)
        xmag(ii,i) = bvac(ia)/xmu0*3.0d0*(xmur_h0-1.0d0)/(xmur_h0+2.0d0)
      enddo
    enddo
  endif

  write(*,*) " magnetization vectors"
  write(*,"(a5,8a15)")"i","x","y","z","Mx","My","Mz","Mnet", "Bi=Mnet*mu0"
  do i=1,n
  rij(:) = xmag(:,i)
  xmag_net = sqrt(dot_product(rij,rij))
  write(*,"(i5,2x,1p7e15.6,0p,f10.5)") i,r(1:3,i),xmag(1:3,i), xmag_net, xmag_net*xmu0
  enddo

  ! now calculate field at remote points
  !
  read(*,*) np ! number of remote points
  if(np.gt.0)then
  allocate(rp(3,np),bp(3,np))
  read(*,*) rp

  call remfld

  bpmax = 0.d0
  write(*,*) " field from magnetized material at remote points"
  write(*,"(a5,8a15)")"i","x","y","z","Bx","By","Bz","Bnet, T", "Bnet, G"
  do ip=1,np
  rij(:) = bp(:,ip)
  bp_net = sqrt(dot_product(rij,rij))
  write(*,"(i5,2x,1p7e15.6,0p,f10.5)") ip,rp(1:3,ip),bp(1:3,ip), bp_net, bp_net*1.d4
  bpmax = max (bpmax, bp_net)
  enddo
  write(*,*) "Max field(G) from magnetized material at remote points = ", bpmax*1.d4


  if(ltotal)then
  write(*,*) " total field at remote points"
  write(*,"(a5,7a15)")"i","x","y","z","Bx","By","Bz","Bnet"
  do ip=1,np
  call coilfld(rp(1,ip),rp(2,ip),rp(3,ip),btot(1),btot(2),btot(3))
  btot(:) = btot(:) +  bp(:,ip)
  bp_net = sqrt(dot_product(btot,btot))
  write(*,"(i5,2x,1p7e15.6)") ip,rp(1:3,ip), btot(1:3), bp_net
  enddo
  endif
  endif

end program multidipole
!
!****************************************************************************
!
module Vfilaments
  integer, parameter :: mf=1000000
  real*8, dimension(mf) :: xf1,yf1,zf1,c12
  integer :: nf
end module Vfilaments
!
!****************************************************************************
!
subroutine readfila(iunit)
  use Vfilaments
  implicit real*8(a-h,o-z)
  read(iunit,*)
  read(iunit,*)
  read(iunit,*)
  nf=0
 5    nf=nf+1
  read(iunit,*,end=99,err=99) xf1(nf),yf1(nf),zf1(nf),c12(nf)
  goto 5
 99   nf=nf-2 ! count filaments, not points
  if(nf.gt.mf)then
    write(102,*)"number of filaments in filamentdump", nf,"exceeds dimension", mf
  stop
  endif
  write(6,*) nf," background coil filaments read in "
  return
end subroutine readfila
!
!****************************************************************************
!
subroutine coilfld (xp,yp,zp, bx,by,bz)
  use Vfilaments
  implicit real*8(a-h,o-z)

  dimension dx(mf),dy(mf),dz(mf),fa(mf)
  dimension vx(mf),vy(mf),vz(mf)

  do i=1,nf
     dx(i) = (xf1(i+1) - xf1(i))*c12(i)
     dy(i) = (yf1(i+1) - yf1(i))*c12(i)
     dz(i) = (zf1(i+1) - zf1(i))*c12(i)

     vx(i) = yf1(i)*dz(i) - zf1(i)*dy(i)
     vy(i) = zf1(i)*dx(i) - xf1(i)*dz(i)
     vz(i) = xf1(i)*dy(i) - yf1(i)*dx(i)

     xp1 = xp - xf1(i)
     yp1 = yp - yf1(i)
     zp1 = zp - zf1(i)

     r1 = sqrt(xp1*xp1 + yp1*yp1 + zp1*zp1)

     xp2 = xp - xf1(i+1)
     yp2 = yp - yf1(i+1)
     zp2 = zp - zf1(i+1)

     r2  = sqrt(xp2*xp2 + yp2*yp2 + zp2*zp2)

     r12    = r2*r1
     fa(i)  = (r2 + r1)/(r12*(r12 + xp2*xp1 + yp2*yp1 + zp2*zp1))

  enddo

  ax     = dot_product(fa(1:nf),dx(1:nf))
  ay     = dot_product(fa(1:nf),dy(1:nf))
  az     = dot_product(fa(1:nf),dz(1:nf))

  bx     = dot_product(fa(1:nf),vx(1:nf)) - yp*az+zp*ay
  by     = dot_product(fa(1:nf),vy(1:nf)) - zp*ax+xp*az
  bz     = dot_product(fa(1:nf),vz(1:nf)) - xp*ay+yp*ax

  bx = bx*1.d-7
  by = by*1.d-7
  bz = bz*1.d-7

  return
end subroutine coilfld
!
!****************************************************************************
!
subroutine filla
  use mydata
  implicit none
  real*8 hnet, xmu_vs_h
  !
  ! calculate murx from H and B-H Curve
  !
  do i = 1,n
      if(iopt.eq.1 .or. iopt.eq.2)then ! use non-linear B-H curve
        hnet = sqrt( dot_product(h(:,i),h(:,i)))
        xmur(i) = xmu_vs_h(hnet)
      else
        xmur(i)=xmur_h0
      endif
  enddo
  !
  ! set up A matrix (field at ia do to dipole at ja AND ALL STELLARATOR IMAGES OF JA)
  !
  A(:,:) = 0.0d0
  do i = 1,n
    do j = 1,n
      if(i.eq.j)then
        do ii = 1,3
        ia = ii+3*(i-1)
        A(ia,ia) = xmu0/3.0d0*(xmur(i)+2.0d0)/(xmur(i)-1.0d0)
        enddo
      else
        if(iopt.eq.1 .or. iopt.eq.3)then ! do field coupling from dipoles
          do is=1,-1,-2 ! stellarator symmetry
          do iper=0,nper-1
           ang= alp*iper
           co = cos(ang)
           si = sin(ang)
           rj(1) = r(1,j)*co - r(2,j)*si
           rj(2) = r(1,j)*si + r(2,j)*co
           rj(2) = rj(2)*is
           rj(3) = r(3,j)*is
!            rij(:) = r(:,i)-r(:,j) ! calculating field at i due to material at j
 !
            rij(:) = r(:,i)-rj ! calculating field at i due to material at j
            r2 = dot_product(rij,rij)
            if(r2.eq.0.0d0)stop "coincident points"
            r1 = sqrt(r2)
            r3 = r2*r1
            rij = rij/r1
            fac = -vol(j)/r3*1.d-7
            do ii = 1,3
              ia = ii+3*(i-1)
              do jj = 1,3
                ja = jj+3*(j-1)

                Aij = 3.0d0*rij(ii)*rij(jj)
                if(ii.eq.jj)then
                  Aij = Aij-1.0d0
                endif
                A(ia,ja) = A(ia,ja)+Aij*fac

              enddo
            enddo
          enddo ! iper
          enddo ! is
        endif
      endif
    enddo
  enddo

  return
end subroutine filla
!
!****************************************************************************
!
subroutine fillb
  use mydata
  implicit none

  ! set up bvac vector - vacuum field from coils
  !
  bvac(:) = 0.0d0
  do i=1,n
    xp = r(1,i)
    yp = r(2,i)
    zp = r(3,i)
    call coilfld (xp,yp,zp, bx,by,bz)
    ii = 3*(i-1)
    bvac(ii+1) = bx
    bvac(ii+2) = by
    bvac(ii+3) = bz
  enddo

  return
end subroutine fillb
!
!****************************************************************************
!
subroutine remfld
  use mydata
  implicit none
  bp=0.d0
  !
  do is=1,-1,-2
  do iper = 0,nper-1
    ang= alp*iper
    co = cos(ang)
    si = sin(ang)
  do i=1,n
  ! image dipole geometry and magnetization vector
    ri(1) = r(1,i)*co - r(2,i)*si
    ri(2) = r(1,i)*si + r(2,i)*co
    ri(2) = ri(2)*is
    ri(3) = r(3,i)*is
    xmagi(1) = xmag(1,i)*co - xmag(2,i)*si
    xmagi(2) = xmag(1,i)*si + xmag(2,i)*co
    xmagi(1) = xmagi(1)*is
    xmagi(3) = xmag(3,i) ! *is
 !
    do ip=1,np
      rij(:)=rp(:,ip)-ri  ! field at remote point rp due to dipole at r
      r2 = dot_product(rij,rij)
      if(r2.eq.0.0d0)stop "remote point coincident with dipole"
      r1 = sqrt(r2)
      r3 = r1*r2
      rij = rij/r1
      fac = vol(i)/r3*1.d-7
      do ii = 1,3  ! dipole component
        do jj = 1,3 ! field component
          db = 3.0d0*rij(ii)*rij(jj)
          if(ii.eq.jj)then
            db = db-1.0d0
          endif
          db = db*fac*xmagi(ii)
          bp(jj,ip) =  bp(jj,ip) + db
        enddo
      enddo
    enddo
  enddo
  enddo ! iper
  enddo ! is
  return
end subroutine remfld
!
!****************************************************************************
!
subroutine fun(x,finish,f_vec,f_jac)
  use mydata
  real*8, intent(in) :: x(:)
  logical, intent(inout) :: finish
  real*8, intent(out) :: f_vec(:)
!  read*8, intent(out), optional :: f_jac(:,:)

  h = reshape( bvac/xmu0 - x/3,(/3,n/))
  call filla
  f_vec = matmul(A,x) - bvac

  return
end subroutine fun
!
!****************************************************************************
!
real*8 function xmu_vs_h(hh)
  use mydata
  real*8 hh
!  xmu_vs_h = 1.3d0
  if(hh.eq.0.d0)then
    xmu_vs_h = xmur_h0
  else
    xmu_vs_h = 1.0d0 + bsat/xmu0/hh*(1.d0-exp(-3*hh/hsat))
  endif
return
end function xmu_vs_h