combine equian and equinum_* subroutines

This change combines `equian` and `equinum_psi` into a new `pol_flux`
subroutine that computes ψ(R,z) and derivatives for either numerical or
analytical equilibria.
Similarly, `equian` (the fpol and dfpol outputs) and `equinum_fpol` are
combined intro `pol_curr` that computes F(ψ) and its derivative.

Callers of these subroutines do not select a specific version based on
the value of `iequil` anymore.

src/coreprofiles.f90

@@ -451,7 +451,7 @@ contains
     ! invalid as they are given assuming a vacuum (N=1)
     if (present(launch_pos)) then
       block
-        use equilibrium,          only : equinum_psi
+        use equilibrium,          only : pol_flux
         use const_and_precisions, only : cm
         real(wp_) :: R, Z, psi
 
@@ -461,7 +461,7 @@ contains
 
         ! Get the poloidal flux at the launcher
         ! Note: this returns -1 when the data is not available
-        call equinum_psi(R, z, psi)
+        call pol_flux(R, z, psi)
 
         if (psi > tail%start .and. psi < tail%end) then
           ! Fall back to the midpoint of ψ₀ and the launcher ψ

src/equilibrium.f90

@@ -52,9 +52,8 @@ module equilibrium
   private
   public read_eqdsk, read_equil_an           ! Reading data files
   public scale_equil, change_cocos           ! Transforming data
-  public equian                              ! Analytical model
-  public equinum_psi, equinum_fpol           ! Interpolated data
   public fq, bfield, tor_curr, tor_curr_psi  ! Accessing local quantities
+  public pol_flux, pol_curr                  !
   public frhotor, frhopol                    ! Converting between poloidal/toroidal flux
   public set_equil_spline, set_equil_an      ! Initialising internal state
   public unset_equil_spline                  ! Deinitialising internal state
@@ -661,7 +660,7 @@ contains
     ! Save Bt value on axis (required in flux_average and used in Jcd def)
     ! and vacuum value B0 at ref. radius data%rvac (used in Jcd_astra def)
 
-    call equinum_fpol(zero, btaxis)
+    call pol_curr(zero, btaxis)
     btaxis = btaxis/rmaxis
     btrcen = fpolas/data%rvac
     rcen   = data%rvac
@@ -803,94 +802,94 @@ contains
   end subroutine set_rho_spline
 
 
-  subroutine equinum_psi(R, z, psi, dpsidr, dpsidz,  &
+  subroutine pol_flux(R, z, psi_n, dpsidr, dpsidz,  &
-                         ddpsidrr, ddpsidzz, ddpsidrz)
+                      ddpsidrr, ddpsidzz, ddpsidrz)
-    ! Computes the normalised poloidal flux ψ (and its derivatives)
+    ! Computes the normalised poloidal flux ψ_n and its
-    ! as a function of (R, z) using the numerical equilibrium
+    ! derivatives wrt (R, z) up to the second order.
+    !
+    ! Note: all output arguments are optional.
+    use gray_params, only : iequil
 
     implicit none
 
-    ! subroutine arguments
-    real(wp_), intent(in)            :: R, z
-    real(wp_), intent(out), optional :: psi, dpsidr, dpsidz
-    real(wp_), intent(out), optional :: ddpsidrr, ddpsidzz, ddpsidrz
-
-    ! Note: here lengths are measured in meters
-    if (R <= rmxm .and. R >= rmnm .and. &
-        z <= zmxm .and. z >= zmnm) then
-
-      if (present(psi))      psi      = psi_spline%eval(R, z)
-      if (present(dpsidr))   dpsidr   = psi_spline%deriv(R, z, 1, 0)
-      if (present(dpsidz))   dpsidz   = psi_spline%deriv(R, z, 0, 1)
-      if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0)
-      if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2)
-      if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1)
-    else
-      if (present(psi))      psi      = -1
-      if (present(dpsidr))   dpsidr   = 0
-      if (present(dpsidz))   dpsidz   = 0
-      if (present(ddpsidrr)) ddpsidrr = 0
-      if (present(ddpsidzz)) ddpsidzz = 0
-      if (present(ddpsidrz)) ddpsidrz = 0
-    end if
-  end subroutine equinum_psi
-
-
-  subroutine equinum_fpol(psi, fpol, dfpol)
-    ! Computes the poloidal current function F(ψ)
-    ! (and its derivative) using the numerical equilibrium
-
-    implicit none
-
-    ! subroutine arguments
-    real(wp_), intent(in)            :: psi
-    real(wp_), intent(out)           :: fpol
-    real(wp_), intent(out), optional :: dfpol
-
-    if(psi <= 1 .and. psi >= 0) then
-      fpol = fpol_spline%eval(psi)
-      if (present(dfpol)) dfpol = fpol_spline%deriv(psi) / psia
-    else
-      fpol = fpolas
-      if (present(dfpol)) dfpol = 0
-    end if
-  end subroutine equinum_fpol
-
-
-  subroutine equian(R, z, psin, fpolv, dfpv, dpsidr, dpsidz,  &
-                    ddpsidrr, ddpsidzz, ddpsidrz)
-    ! Computes all MHD equilibrium parameters from a simple analytical model
-    implicit none
-
     ! subroutine arguments
     real(wp_), intent(in) :: R, z
-    real(wp_), intent(out), optional :: psin, fpolv, dfpv, dpsidr, dpsidz,  &
+    real(wp_), intent(out), optional :: psi_n, dpsidr, dpsidz,  &
                                         ddpsidrr, ddpsidzz, ddpsidrz
 
-    ! variables
+    ! local variables
     real(wp_) :: r_p   ! poloidal radius
     real(wp_) :: rho_p ! poloidal radius normalised
 
-    ! ρ_p is just the geometrical poloidal radius divided by a
+    if (iequil < 2) then
-    r_p = hypot(R - model%R0, z - model%z0)
+      ! Analytical model
-    rho_p = r_p / model%a
+      ! ρ_p is just the geometrical poloidal radius divided by a
+      r_p = hypot(R - model%R0, z - model%z0)
+      rho_p = r_p / model%a
 
-    ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
+      ! ψ_n = ρ_p²(R,z) = [(R-R₀)² + (z-z₀)²]/a²
-    if (present(psin)) psin = rho_p**2
+      if (present(psi_n)) psi_n = rho_p**2
 
-    ! ∂ψ/∂R, ∂ψ/∂z
+      ! ∂ψ_n/∂R, ∂ψ_n/∂z
-    if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2
+      if (present(dpsidr)) dpsidr = 2*(R - model%R0)/model%a**2
-    if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2
+      if (present(dpsidz)) dpsidz = 2*(z - model%z0)/model%a**2
 
-    ! ∂²ψ/∂R², ∂²ψ/∂z², ∂²ψ/∂R∂z
+      ! ∂²ψ_n/∂R², ∂²ψ_n/∂z², ∂²ψ_n/∂R∂z
-    if (present(ddpsidrr)) ddpsidrr = 2/model%a**2
+      if (present(ddpsidrr)) ddpsidrr = 2/model%a**2
-    if (present(ddpsidzz)) ddpsidzz = 2/model%a**2
+      if (present(ddpsidzz)) ddpsidzz = 2/model%a**2
-    if (present(ddpsidrz)) ddpsidrz = 0
+      if (present(ddpsidrz)) ddpsidrz = 0
+    else
+      ! Numerical data
+      if (R <= rmxm .and. R >= rmnm .and. &
+          z <= zmxm .and. z >= zmnm) then
+        ! Within the interpolation range
+        if (present(psi_n))    psi_n    = psi_spline%eval(R, z)
+        if (present(dpsidr))   dpsidr   = psi_spline%deriv(R, z, 1, 0)
+        if (present(dpsidz))   dpsidz   = psi_spline%deriv(R, z, 0, 1)
+        if (present(ddpsidrr)) ddpsidrr = psi_spline%deriv(R, z, 2, 0)
+        if (present(ddpsidzz)) ddpsidzz = psi_spline%deriv(R, z, 0, 2)
+        if (present(ddpsidrz)) ddpsidrz = psi_spline%deriv(R, z, 1, 1)
+      else
+        ! Outside
+        if (present(psi_n))    psi_n    = -1
+        if (present(dpsidr))   dpsidr   = 0
+        if (present(dpsidz))   dpsidz   = 0
+        if (present(ddpsidrr)) ddpsidrr = 0
+        if (present(ddpsidzz)) ddpsidzz = 0
+        if (present(ddpsidrz)) ddpsidrz = 0
+      end if
+    end if
+  end subroutine pol_flux
 
-    ! F(ψ) = B₀⋅R₀, a constant
+
-    if (present(fpolv)) fpolv = model%B0 * model%R0
+  subroutine pol_curr(psi_n, fpol, dfpol)
-    if (present(dfpv))   dfpv = 0
+    ! Computes the poloidal current function F(ψ_n)
-  end subroutine equian
+    ! and (optionally) its derivative dF/dψ given ψ_n.
+    use gray_params, only : iequil
+
+    implicit none
+
+    ! function arguments
+    real(wp_), intent(in)            :: psi_n  ! normalised poloidal flux
+    real(wp_), intent(out)           :: fpol   ! poloidal current
+    real(wp_), intent(out), optional :: dfpol  ! derivative
+
+    if (iequil < 2) then
+      ! Analytical model
+      ! F(ψ) = B₀⋅R₀, a constant
+      fpol = model%B0 * model%R0
+      if (present(dfpol)) dfpol = 0
+    else
+      ! Numerical data
+      if(psi_n <= 1 .and. psi_n >= 0) then
+        fpol = fpol_spline%eval(psi_n)
+        if (present(dfpol)) dfpol = fpol_spline%deriv(psi_n) / psia
+      else
+        fpol = fpolas
+        if (present(dfpol)) dfpol = 0
+      end if
+    end if
+  end subroutine pol_curr
 
 
   function frhotor(rho_p)
@@ -1017,83 +1016,93 @@ contains
   end function fq
 
 
-  subroutine bfield(rpsim, zpsim, bphi, br, bz)
+  subroutine bfield(R, z, B_R, B_z, B_phi)
     ! Computes the magnetic field as a function of
     ! (R, z) in cylindrical coordinates
+    !
+    ! Note: all output arguments are optional.
     use gray_params, only : iequil
 
     implicit none
 
     ! subroutine arguments
-    real(wp_), intent(in) :: rpsim, zpsim
+    real(wp_), intent(in) :: R, z
-    real(wp_), intent(out), optional :: bphi,br,bz
+    real(wp_), intent(out), optional :: B_R, B_z, B_phi
 
     ! local variables
-    real(wp_) :: psin,fpol
+    real(wp_) :: psi_n, fpol, dpsidr, dpsidz
 
-    if (iequil < 2) then
+    call pol_flux(R, z, psi_n, dpsidr, dpsidz)
-      ! Analytical model
+    call pol_curr(psi_n, fpol)
-      call equian(rpsim,zpsim,fpolv=bphi,dpsidr=bz,dpsidz=br)
-      if (present(bphi)) bphi=bphi/rpsim
-    else
-      ! Numerical data
-      call equinum_psi(rpsim,zpsim,psi=bphi,dpsidr=bz,dpsidz=br)
-      if (present(bphi)) then
-        psin=bphi
-        call equinum_fpol(psin,fpol)
-        bphi=fpol/rpsim
-      end if
-    end if
 
-    if (present(br)) br=-br/rpsim
+    ! The field in cocos=3 is given by
-    if (present(bz)) bz= bz/rpsim
+    !
+    !   B = F(ψ)∇φ + ∇ψ×∇φ.
+    !
+    ! Writing the gradient of ψ=ψ(R,z) as
+    !
+    !   ∇ψ = ∂ψ/∂R ∇R + ∂ψ/∂z ∇z,
+    !
+    ! and carrying out the cross products:
+    !
+    !   B = F(ψ)∇φ - ∂ψ/∂z ∇R/R + ∂ψ/∂R ∇z/R
+    !
+    if (present(B_R)) B_R = - 1/R * dpsidz * psia
+    if (present(B_z)) B_z = + 1/R * dpsidr * psia
+    if (present(B_phi)) B_phi = fpol / R
   end subroutine bfield
 
 
-  function tor_curr(rpsim,zpsim) result(jphi)
+  function tor_curr(R, z) result(J_phi)
-    ! Computes the toroidal current Jφ as a function of (R, z)
+    ! Computes the toroidal current J_φ as a function of (R, z)
-    use const_and_precisions, only : ccj=>mu0inv
+    use const_and_precisions, only : mu0_
-    use gray_params,          only : iequil
 
     implicit none
 
     ! function arguments
-    real(wp_), intent(in) :: rpsim, zpsim
+    real(wp_), intent(in) :: R, z
-    real(wp_)             :: jphi
+    real(wp_)             :: J_phi
 
     ! local variables
-    real(wp_) :: bzz,dbvcdc13,dbvcdc31
+    real(wp_) :: dB_Rdz, dB_zdR              ! derivatives of B_R, B_z
-    real(wp_) :: dpsidr,ddpsidrr,ddpsidzz
+    real(wp_) :: dpsidr, ddpsidrr, ddpsidzz  ! derivatives of ψ_n
 
-    if(iequil < 2) then
+    call pol_flux(R, z, dpsidr=dpsidr, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)
-      ! Analytical model
+
-      call equian(rpsim,zpsim,dpsidr=dpsidr, &
+    ! In the usual MHD limit we have ∇×B = μ₀J. Using the
-                  ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    ! curl in cylindrical coords the toroidal current is
-    else
+    !
-      ! Numerical data
+    !   J_φ = 1/μ₀ (∇×B)_φ = 1/μ₀ [∂B_R/∂z - ∂B_z/∂R].
-      call equinum_psi(rpsim,zpsim,dpsidr=dpsidr, &
+    !
-                       ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
+    ! Finally, from B = F(ψ)∇φ + ∇ψ×∇φ we find:
-    end if
+    !
-    bzz= dpsidr/rpsim
+    !   B_R = - 1/R ∂ψ/∂z,
-    dbvcdc13=-ddpsidzz/rpsim
+    !   B_z = + 1/R ∂ψ/∂R,
-    dbvcdc31= ddpsidrr/rpsim-bzz/rpsim
+    !
-    jphi=ccj*(dbvcdc13-dbvcdc31)
+    ! from which:
+    !
+    !   ∂B_R/∂z = - 1/R ∂²ψ/∂z²
+    !   ∂B_z/∂R = + 1/R ∂²ψ/∂R² - 1/R² ∂ψ/∂R.
+    !
+    dB_Rdz = - 1/R * ddpsidzz * psia
+    dB_zdR = + 1/R * (ddpsidrr - 1/R * dpsidr) * psia
+    J_phi  = 1/mu0_ * (dB_Rdz - dB_zdR)
   end function tor_curr
 
 
-  function tor_curr_psi(psin) result(jphi)
+  function tor_curr_psi(psi_n) result(J_phi)
-    ! Computes the toroidal current Jφ as a function of ψ
+    ! Computes the toroidal current J_φ as a function of ψ
 
     implicit none
 
     ! function arguments
-    real(wp_), intent(in) :: psin
+    real(wp_), intent(in) :: psi_n
-    real(wp_)             :: jphi
+    real(wp_)             :: J_phi
 
     ! local variables
-    real(wp_) :: r1, r2
+    real(wp_) :: R1, R2
-    call psi_raxis(psin, r1, r2)
+    call psi_raxis(psi_n, R1, R2)
-    jphi = tor_curr(r2, zmaxis)
+    J_phi = tor_curr(R2, zmaxis)
   end function tor_curr_psi
 
 
@@ -1178,7 +1187,7 @@ contains
     end if
     rf=xvec(1)
     zf=xvec(2)
-    call equinum_psi(rf,zf,psinvf)
+    call pol_flux(rf, zf, psinvf)
   end subroutine points_ox
 
 
@@ -1199,12 +1208,12 @@ contains
 
     select case(iflag)
     case(1)
-      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz)
+      call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz)
       fvec(1) = dpsidr
       fvec(2) = dpsidz
     case(2)
-      call equinum_psi(x(1),x(2),ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz, &
+      call pol_flux(x(1), x(2), ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz, &
-           ddpsidrz=ddpsidrz)
+                    ddpsidrz=ddpsidrz)
       fjac(1,1) = ddpsidrr
       fjac(1,2) = ddpsidrz
       fjac(2,1) = ddpsidrz
@@ -1272,12 +1281,12 @@ contains
 
     select case(iflag)
     case(1)
-      call equinum_psi(x(1),x(2),psinv,dpsidr)
+      call pol_flux(x(1), x(2), psinv, dpsidr)
       fvec(1) = psinv-f0(1)
       fvec(2) = dpsidr-f0(2)
     case(2)
-      call equinum_psi(x(1),x(2),dpsidr=dpsidr,dpsidz=dpsidz, &
+      call pol_flux(x(1), x(2), dpsidr=dpsidr, dpsidz=dpsidz, &
-           ddpsidrr=ddpsidrr,ddpsidrz=ddpsidrz)
+                    ddpsidrr=ddpsidrr, ddpsidrz=ddpsidrz)
       fjac(1,1) = dpsidr
       fjac(1,2) = dpsidz
       fjac(2,1) = ddpsidrr

src/gray_core.f90

@@ -1260,7 +1260,7 @@ contains
                         xg, yg, derxg, deryg, psinv_opt)
     use const_and_precisions,  only : zero, cm
     use gray_params,   only : iequil
-    use equilibrium,   only : psia, equinum_fpol, equinum_psi, equian, sgnbphi
+    use equilibrium,   only : psia, pol_flux, pol_curr, sgnbphi
     use coreprofiles,  only : density
 
     implicit none
@@ -1322,13 +1322,8 @@ contains
     zzm = 1.0e-2_wp_*zz
     rrm = 1.0e-2_wp_*rr
 
-    if(iequil==1) then
+    call pol_flux(rrm, zzm, psinv, dpsidr, dpsidz, ddpsidrr, ddpsidzz, ddpsidrz)
-      call equian(rrm,zzm,psinv,fpolv,dfpv,dpsidr,dpsidz, &
+    call pol_curr(psinv, fpolv, dfpv)
-                  ddpsidrr,ddpsidzz,ddpsidrz)
-    else
-      call equinum_psi(rrm,zzm,psinv,dpsidr,dpsidz,ddpsidrr,ddpsidzz,ddpsidrz)
-      call equinum_fpol(psinv,fpolv,dfpv)
-    end if
 
     ! copy optional output
     if (present(psinv_opt)) psinv_opt = psinv
@@ -1886,7 +1881,7 @@ contains
         rm=sqrt(xmv(1)**2+xmv(2)**2)
         csphi=xmv(1)/rm
         snphi=xmv(2)/rm
-        call bfield(rm,xmv(3),bphi,br,bz)
+        call bfield(rm,xmv(3),br,bz,bphi)
         ! bv(i)  = B_i in cartesian coordinates
         bv(1)=br*csphi-bphi*snphi
         bv(2)=br*snphi+bphi*csphi
@@ -2266,7 +2261,7 @@ bb: do
       zzk = zv(k)
       do j = 1, q_spline%ndata
         rrj = rv(j)
-        call bfield(rrj, zzk, bbphi, bbr, bbz)
+        call bfield(rrj, zzk, bbr, bbz, bbphi)
         btotal(j,k) = sqrt(bbr**2 + bbz**2 + bbphi**2)
         if(btotal(j,k) >= btmx) btmx = btotal(j,k)
         if(btotal(j,k) <= btmn) btmn = btotal(j,k)
@@ -2295,9 +2290,8 @@ bb: do
   subroutine print_maps(bres, xgcn, r0, anpl0)
     ! Prints several input quantities on a regular grid (unit 72)
 
-    use gray_params,  only : iequil
+    use equilibrium,  only : rmnm, rmxm, zmnm, zmxm, &
-    use equilibrium,  only : rmnm, rmxm, zmnm, zmxm, equian, equinum_psi, &
+                             pol_flux, pol_curr, q_spline
-                             equinum_fpol, q_spline
     use coreprofiles, only : density, temp
     use units,        only : umaps, unit_active
 
@@ -2327,12 +2321,8 @@ bb: do
       anpl = anpl0 * r0/rj
       do k = 1, q_spline%ndata
         zk = z(k)
-        if (iequil < 2) then
+        call pol_flux(rj, zk, psin, dpsidr=bz, dpsidz=br)
-          call equian(rj, zk, psin=psin, fpolv=bphi, dpsidr=bz, dpsidz=br)
+        call pol_curr(psin, fpol=bphi)
-        else
-          call equinum_psi(rj, zk, psi=psin, dpsidr=bz, dpsidz=br)
-          call equinum_fpol(psin, fpol=bphi)
-        end if
         br   = -br/rj
         bphi =  bphi/rj
         bz   =  bz/rj

src/magsurf_data.f90

@@ -101,9 +101,8 @@ contains
 
   subroutine flux_average
     use const_and_precisions, only : wp_,zero,one,pi,ccj=>mu0inv
-    use gray_params, only : iequil
     use equilibrium, only : btrcen,btaxis,rmaxis,zmaxis,phitedge,zbsup,zbinf, &
-                            equian,equinum_psi,bfield,frhotor,fq,tor_curr,psia
+                            bfield,frhotor,fq,tor_curr,psia,pol_flux
     use dierckx,     only : regrid,coeff_parder
     implicit none
 
@@ -170,11 +169,7 @@ contains
     ratio_cdbtor=1.0_wp_
     ratio_pltor=1.0_wp_
     fc=1.0_wp_
-    if(iequil < 2) then
+    call pol_flux(rmaxis, zmaxis, ddpsidrr=ddpsidrr, ddpsidzz=ddpsidzz)
-      call equian(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
-    else
-      call equinum_psi(rmaxis,zmaxis,ddpsidrr=ddpsidrr,ddpsidzz=ddpsidzz)
-    end if       
     qq=abs(btaxis)/sqrt(ddpsidrr*ddpsidzz)/psia
     ajphiav=-ccj*(ddpsidrr+ddpsidzz)*psia/rmaxis
     
@@ -230,7 +225,7 @@ contains
       bmmn=1.0e+30_wp_
 
       ajphi0 = tor_curr(rctemp(1),zctemp(1))
-      call bfield(rctemp(1),zctemp(1),bphi,br=brr,bz=bzz)
+      call bfield(rctemp(1), zctemp(1), brr, bzz, bphi)
       fpolv=bphi*rctemp(1)
       btot0=sqrt(bphi**2+brr**2+bzz**2)
       bpoloid0=sqrt(brr**2+bzz**2)
@@ -256,7 +251,7 @@ contains
 ! compute line integrals on the contour psi=psinjp=height^2
         rpsim=rctemp(inc1)
         zpsim=zctemp(inc1)
-        call bfield(rpsim,zpsim,br=brr,bz=bzz)
+        call bfield(rpsim, zpsim, brr, bzz)
         ajphi = tor_curr(rpsim,zpsim)
         bphi=fpolv/rpsim
         btot=sqrt(bphi**2+brr**2+bzz**2)

src/multipass.f90

@@ -37,7 +37,7 @@ contains
     rm=sqrt(xmv(1)**2+xmv(2)**2)
     csphi=xmv(1)/rm
     snphi=xmv(2)/rm
-    call bfield(rm,xmv(3),bphi,br,bz)
+    call bfield(rm,xmv(3),br,bz,bphi)
     bv(1)=br*csphi-bphi*snphi
     bv(2)=br*snphi+bphi*csphi
     bv(3)=bz
@@ -82,7 +82,7 @@ contains
     rm=sqrt(xmv(1)**2+xmv(2)**2)
     csphi=xmv(1)/rm
     snphi=xmv(2)/rm
-    call bfield(rm,xmv(3),bphi,br,bz)
+    call bfield(rm,xmv(3),br,bz,bphi)
     bv(1)=br*csphi-bphi*snphi
     bv(2)=br*snphi+bphi*csphi
     bv(3)=bz