arclength.f90

Go to the documentation of this file.

subroutine arclength(fr, mfm, theta, dtc, wr, ws)
!-----------------------------------------------------------------------
! subroutine to calculate the arclength of the plasma boundary
! routine returns the theta values at equidistant arclength
! parameters :
!                fr    : fourier series of boundary (input)
!                mfm   : number of harmonics in fr  (input)
!                theta : resulting values of theta
!                dtc   : the derivative of theta to equidistant angle
!-----------------------------------------------------------------------
  use itm_types
  use mod_parameter
  use mod_mesh
  
  implicit none
  
  integer(itm_i4), intent(in) :: mfm
  real(r8), dimension(mfm + 2), intent(in) :: fr     
  real(r8), dimension(4), intent(in) :: wr, ws
  
  real(r8), dimension(np), intent(out) :: theta, dtc
  
  integer(itm_i4), parameter :: npts = 512
  
  real(r8), dimension(npts) :: xl, dxl
  real(r8) :: ti1, ti2
  real(r8) :: dl, dr, tk, rr
  real(r8) :: chil, dct
  integer(itm_i4) :: i, j, k, m
  
  xl(1) = 0._r8
!---------------------------------- calculate lenght(theta)
  do j = 1, npts - 1
    ti1 = (j - 1) / dble(npts - 1) * (ias + 1) * pi
    ti2 = j / dble(npts - 1) * (ias + 1) * pi
!---------------------------------- 4 point gaussian integration
    dl = 0._r8
    do k = 1, 4
      tk = ti1 + (ti2 - ti1) * (wr(k) + 1._r8) / 2._r8
      rr = fr(1) / 2._r8
      dr = 0._r8
      do m = 2, mfm / 2
        rr = rr + fr(2 * m - 1) * cos((m - 1) * tk) + fr(2 * m) &
         * sin((m - 1) * tk)
        dr = dr - (m - 1) * fr(2 * m - 1) * sin((m - 1) * tk) &
         + (m - 1) * fr(2 * m) * cos((m - 1) * tk)
      end do
      dl = dl + sqrt(rr * rr + dr * dr) * ws(k)
    end do
    xl(j + 1) = xl(j) + dl * (ias + 1) / (2._r8 * pi * (npts - 1))
    rr = fr(1) / 2._r8
    dr = 0._r8
    do m = 2, mfm / 2
      rr = rr + fr(2 * m - 1) * cos((m - 1) * ti2) + fr(2 * m) &
       * sin((m - 1._r8) * ti2)
      dr = dr - (m - 1) * fr(2 * m - 1) * sin((m - 1) * ti2) &
       + (m - 1) * fr(2 * m) * cos((m - 1) * ti2)
    end do
    dxl(j + 1) = sqrt(rr * rr + dr * dr)
  end do
  ti2 = 0._r8
  rr = fr(1) / 2._r8
  dr = 0._r8
  
  do m = 2, mfm / 2
    rr = rr + fr(2 * m - 1) * cos((m - 1) * ti2) + fr(2 * m) &
     * sin((m - 1) * ti2)
    dr = dr - (m - 1) * fr(2 * m- 1 ) * sin((m - 1) * ti2) &
     + (m - 1) * fr(2 * m) * cos((m - 1) * ti2)
  end do
  dxl(1) = sqrt(rr * rr + dr * dr)
  
  theta(1) = 0._r8
  theta(np) = (ias + 1) * pi
  dtc(1) = 1._r8 / (dxl(1) * (ias + 1) * pi / xl(npts))
  dtc(np) = 1._r8 / (dxl(npts) * (ias + 1) * pi / xl(npts))
  i = 1
  do j = 2, np - 1
    chil = (j - 1) / dble(np - 1) * xl(npts)
    do while ((chil > xl(i)) .and. (i < npts))
      i = i + 1
    end do
!--------------------------- interval located, use linear interpolation
    theta(j) = (dble(i - 1) + (chil - xl(i)) / (xl(i + 1) &
     - xl(i))) / dble(npts - 1) * (ias + 1) * pi
    dct = ((chil - xl(i)) / (xl(i + 1) - xl(i)) * (dxl(i + 1) &
     - dxl(i)) + dxl(i)) * (ias + 1) * pi / xl(npts)
    dtc(j) = 1._r8 / dct
  end do
  
  return
end subroutine arclength