mod_dipole.t

Go to the documentation of this file.

!> sets up a magnetic dipole in a 3D cartesian box
!  parameters for dipole moment: 
!      strength mu_dipole (factor 4pi included)
!      polar angle in degrees theta_dipole [0 to 90]
!      azimuthal angle in degrees phi_dipole [0 to 360]
!  center of dipole at x1_dip, x2_dip, x3_dip
module mod_dipole
  implicit none
  double precision :: mu_dipole,theta_dipole,phi_dipole,x1_dip,x2_dip,x3_dip
  double precision :: mx, my, mz
 
contains
 
{^ifthreed
  subroutine dipole(ixI^L,ixO^L,x,B)
    use mod_global_parameters
    use mod_comm_lib, only: mpistop
 
    integer, intent(in) :: ixI^L, ixO^L
    double precision, intent(in) :: x(ixI^S,1:ndim)
    double precision, intent(inout) :: B(ixI^S,1:ndir)
 
    double precision :: xsincos(ixI^S),xsinsin(ixI^S),xcos(ixI^S)
    double precision :: rr0(ixI^S),mdotr(ixI^S)
 
    b=0.0d0
    mx=mu_dipole*dsin(theta_dipole*dpi/180.0d0)*dcos(phi_dipole*dpi/180.d0)
    my=mu_dipole*dsin(theta_dipole*dpi/180.0d0)*dsin(phi_dipole*dpi/180.0d0)
    mz=mu_dipole*dcos(theta_dipole*dpi/180.0d0)
    rr0(ixo^s)=dsqrt((x(ixo^s,1)-x1_dip)**2 &
                    +(x(ixo^s,2)-x2_dip)**2 &
                    +(x(ixo^s,3)-x3_dip)**2)
    if(any(rr0(ixo^s)<smalldouble)) call mpistop("dipole center at cell center")
    !where(rr0(ixO^S)<smalldouble) 
    !   rr0(ixO^S)=bigdouble
    !endwhere
    xsincos(ixo^s)=(x(ixo^s,1)-x1_dip)/rr0(ixo^s)
    xsinsin(ixo^s)=(x(ixo^s,2)-x2_dip)/rr0(ixo^s)
    xcos(ixo^s)   =(x(ixo^s,3)-x3_dip)/rr0(ixo^s)
    mdotr(ixo^s)=mx*xsincos(ixo^s)+my*xsinsin(ixo^s)+mz*xcos(ixo^s)
    b(ixo^s,1)=(3.0d0*xsincos(ixo^s)*mdotr(ixo^s)-mx)/rr0(ixo^s)**3
    b(ixo^s,2)=(3.0d0*xsinsin(ixo^s)*mdotr(ixo^s)-my)/rr0(ixo^s)**3
    b(ixo^s,3)=(3.0d0*xcos(ixo^s)   *mdotr(ixo^s)-mz)/rr0(ixo^s)**3
 
  end subroutine dipole
}
end module mod_dipole