!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
!
!_______________________________________________________________________
! Donev Aleksandar, PHY201, a Fortran 90 program
! that find the magnetic field of a circular ring (11/6/99)
!_______________________________________________________________________
!
!_______________________________________________________________________
! This module contains the functions that give the
! magnetic field elements dBx and dBz for a current-carrying ring
! placed in the xy plane and for y=0.
! Both functions are Fortran 90 versions of elemental and pure
! functions in Fortran 95 (not supported by the F compiler)
!_______________________________________________________________________
MODULE dB
      USE Numerical_Integration, ONLY: wp
      IMPLICIT NONE
      PUBLIC :: dBx, dBz
      REAL (KIND=wp), PUBLIC :: x, z
!
CONTAINS
  !___________________________________
  ! This is the x component of dB:
  !___________________________________
      FUNCTION dBx (theta) RESULT (Bx)
         REAL (KIND=wp), INTENT (IN) :: theta
         REAL (KIND=wp) :: Bx
         Bx = z * Cos (theta) / (1+x**2+z**2-2*x*Cos(theta)) ** 1.5
      END FUNCTION dBx
!
  !___________________________________
  ! This is the z component of dB:
  !___________________________________
      FUNCTION dBz (theta) RESULT (Bz)
         REAL (KIND=wp), INTENT (IN) :: theta
         REAL (KIND=wp) :: Bz
         Bz = (1-x*Cos(theta)) / (1+x**2+z**2-2*x*Cos(theta)) ** 1.5
      END FUNCTION dBz
!
END MODULE dB
!_______________________________________________________________________
!
!_______________________________________________________________________
MODULE B
      USE dB
      USE Numerical_Integration, ONLY: AdaptiveIntegral, error_flag, wp
      IMPLICIT NONE
      PUBLIC
!
      INTEGER, SAVE :: error_status
      REAL (KIND=wp), PARAMETER, PRIVATE :: pi = 3.141592653589793116_wp
!
CONTAINS
  !___________________________________
  ! This is the x component of B:
  !___________________________________
      FUNCTION Bx (x_, z_) RESULT (Bx_)
         REAL (KIND=wp), INTENT (IN) :: x_, z_
         REAL (KIND=wp) :: Bx_
         REAL (KIND=wp) :: error
    ! I must use these dirty methods because in Fortran
    ! we are restricted about the interface of the function
    ! we can pass on to ready-made routines:
         x = x_
         z = z_
         Bx_ = AdaptiveIntegral (f=dBx, a=0.0_wp, B=2.0_wp*pi, error_estimate=error)
         error_status = error_flag
      END FUNCTION Bx
!
  !___________________________________
  ! This is the z component of B:
  !___________________________________
      FUNCTION Bz (x_, z_) RESULT (Bz_)
         REAL (KIND=wp), INTENT (IN) :: x_, z_
         REAL (KIND=wp) :: Bz_
         REAL (KIND=wp) :: error
         x = x_
         z = z_
         Bz_ = AdaptiveIntegral (f=dBz, a=0.0_wp, B=2.0_wp*pi, error_estimate=error)
         error_status = error_flag
      END FUNCTION Bz
!
      FUNCTION B_ (x_, z_) RESULT (Bxz)
         USE Precision, ONLY : sp, dp
         REAL (KIND=sp), INTENT (IN) :: x_, z_
         REAL (KIND=sp), DIMENSION (2) :: Bxz
         REAL (KIND=wp) :: error
         x = real (x_, wp)
         z = real (z_, wp)
         Bxz (1) = real (AdaptiveIntegral(f=dBx, a=0.0_wp, B=2.0_wp*pi, error_estimate=error), sp)
         Bxz (2) = real (AdaptiveIntegral(f=dBz, a=0.0_wp, B=2.0_wp*pi, error_estimate=error), sp)
      END FUNCTION B_
!
END MODULE B
!_______________________________________________________________________
!
!_______________________________________________________________________
! The main program B_ring gives the value of the magnetic field
! Bx and Bz for a current-carrying ring placed in the xy plane for y=0
! by calculating a numerical integral using Newton-Cottes formulae.
!_______________________________________________________________________
PROGRAM B_ring
      USE Precision, ONLY : sp, dp
      USE B, ONLY: Bx, Bz, B_, error_status, wp
      USE SimpleGraphics
      USE FunGraphics
      IMPLICIT NONE
!
      REAL (KIND=wp) :: x, z
      REAL (KIND=sp), PARAMETER :: length = 2.5644_wp
      CHARACTER (LEN=100), DIMENSION (2) :: title
!
      WRITE (UNIT=*, FMT=*) "Enter the values of x and z:"
      READ (UNIT=*, FMT=*) x, z
!
      WRITE (UNIT=*, FMT=*) "Bx=", Bx (x, z), " with error status", error_status
      WRITE (UNIT=*, FMT=*) "Bz=", Bz (x, z), " with error status", error_status
!
      title (1) = "The magnetic field of a current-carrying ring"
      title (2) = "A. Donev, 10/19/00"
      CALL InitGraphics (file="ring.png", file_type="PNG", plot_title=title, x_label="x", y_label="z", z_label="|B(x,z)|")
      CALL FunVectorPlot2D (f_xy=B_, xy_range= (/-length, length,-length, length /), num_points= (/ 25, 25 /), plot_spec="11G", &
                            & axis= (/-length, length,-length, length, 0.1, 5.0 /), zoom=0.2)
      CALL EndGraphics ()
!
END PROGRAM B_ring