MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_rbsl.t
Go to the documentation of this file.
1 !> To get a RBSL magnetic flux rope in 3D (Titov 2018 ApJL 852, L21)
2 module mod_rbsl
3  implicit none
4 
5 contains
6 
7 {^ifthreed
8  subroutine rbsl(ixI^L,ixO^L,np,a,F_flx,positive_helicity,x,x_axis,Atotal,Bfr)
10  use mod_geometry
11 
12  integer, intent(in) :: ixI^L,ixO^L
13  ! resolution of flux rope axis
14  integer, intent(in) :: np
15  double precision, intent(in) :: x(ixI^S,1:ndim)
16  ! coordinates of flux rope axis (integral path)
17  double precision, intent(in) :: x_axis(np,1:ndim)
18  ! cross-sectional radius of flux rope
19  double precision, intent(in) :: a
20  ! net magnetic flux along flux rope axis
21  double precision, intent(in) :: F_flx
22  ! is positive helicity
23  logical, intent(in) :: positive_helicity
24  ! vector potential of flux rope
25  double precision, intent(out) :: Atotal(ixI^S,1:ndim)
26  ! magnetic field of flux rope
27  double precision, optional, intent(out) :: Bfr(ixI^S,1:ndir)
28 
29  ! net current along flux rope axis for azimuthal magnetic field
30  double precision :: I_cur
31  ! vector potential for azimuthal magnetic field
32  double precision :: AIx(ixI^S,1:ndim)
33  ! vector potential for axial magnetic field
34  double precision :: AFx(ixI^S,1:ndim)
35  integer :: ix^D, ixp, idirmin
36  double precision :: r_mag, KIr, KFr, dl, re_pi, sqrt1r, f52r,fsqrt6
37  double precision :: Rpl(1:ndim), r_vec(1:ndim), Rcr(1:ndim)
38 
39  if(positive_helicity) then
40  i_cur = 5.d0*sqrt(2.d0)*f_flx/(3.d0*4.d0*dpi*a)
41  else
42  i_cur =-5.d0*sqrt(2.d0)*f_flx/(3.d0*4.d0*dpi*a)
43  end if
44 
45  re_pi=1.d0/dpi
46  aix = 0.d0
47  afx = 0.d0
48  fsqrt6=1.d0/sqrt(6.d0)
49  {do ix^db=ixomin^db,ixomax^db\}
50  do ixp=1,np
51  ! position vector from source point to field point
52  r_vec(:) = (x(ix^d,:) - x_axis(ixp,:))/a
53  r_mag = sqrt(sum(r_vec(:)**2))
54  ! calculate tangential vector of the axis which is a circle
55  if (ixp == 1) then
56  rpl(:) = 0.5d0*(x_axis(ixp+1,:)-x_axis(np,:))
57  else if (ixp==np) then
58  rpl(:) = 0.5d0*(x_axis(1,:)-x_axis(ixp-1,:))
59  else
60  rpl(:) = 0.5d0*(x_axis(ixp+1,:)-x_axis(ixp-1,:))
61  end if
62  ! Rpl X r_vec
63  rcr(1) = rpl(2)*r_vec(3) - rpl(3)*r_vec(2)
64  rcr(2) = rpl(3)*r_vec(1) - rpl(1)*r_vec(3)
65  rcr(3) = rpl(1)*r_vec(2) - rpl(2)*r_vec(1)
66  if (r_mag < 1.d0) then
67  sqrt1r=sqrt(1.d0-r_mag**2)
68  f52r=5.d0-2.d0*r_mag**2
69  kir = 2.d0*re_pi*(asin(r_mag)/r_mag + f52r*third*sqrt1r)
70  kfr = 2.d0*re_pi/r_mag**2*(asin(r_mag)/r_mag-sqrt1r) + &
71  2.d0*re_pi*sqrt1r + f52r*0.5d0*fsqrt6*(1.d0 - &
72  2.d0*re_pi*asin((1.d0+2.d0*r_mag**2)/f52r))
73  else
74  kir = 1.d0/r_mag
75  kfr = kir**3
76  endif
77  aix(ix^d,:) = aix(ix^d,:) + kir*rpl(:)
78  afx(ix^d,:) = afx(ix^d,:) + kfr*rcr(:)
79  end do
80  aix(ix^d,:) = aix(ix^d,:)*i_cur/a
81  afx(ix^d,:) = afx(ix^d,:)*f_flx*0.25d0*re_pi/a**2
82  {end do\}
83  atotal=aix+afx
84  if(present(bfr)) call curlvector(atotal,ixi^l,ixo^l,bfr,idirmin,1,ndir)
85 
86  end subroutine rbsl
87 }
88 
89 end module mod_rbsl
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
This module contains definitions of global parameters and variables and some generic functions/subrou...
To get a RBSL magnetic flux rope in 3D (Titov 2018 ApJL 852, L21)
Definition: mod_rbsl.t:2
subroutine rbsl(ixIL, ixOL, np, a, F_flx, positive_helicity, x, x_axis, Atotal, Bfr)
Definition: mod_rbsl.t:9