Solution to Worksheet 8 of phy201

!_______________________________________________________________________
! 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 Precision
  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 Precision
  use dB
  use Numerical_Integration, only: AdaptiveIntegral,error_flag
  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)
    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
  use B, only: Bx, Bz, B_, error_status
  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="XWIN", & 
                    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