shape_from_points.f90

Go to the documentation of this file.

subroutine shape_from_points(R_bnd, Z_bnd, n_bnd, mfm, Rgeo, amin, fr)
!---------------------------------------------------------------------------
! subroutine interpolates the given n_bnd boundary points (R_bnd,Z_bnd)
! onto a regular grid in polar theta with mfm points
!---------------------------------------------------------------------------

  use itm_types
  use mod_dat, only : zgeo, verbosity
  use mod_parameter, only : twopi
  use mod_helena_io, only : iu6
  use qsort_c_module

  implicit none

  real(r8), intent(in) :: r_bnd(n_bnd), z_bnd(n_bnd)
  integer(itm_i4), intent(in) :: mfm
  integer(itm_i4), intent(inout) :: n_bnd
  real(r8), intent(out) :: fr(mfm)
  real(r8), intent(out) :: rgeo, amin

  real(r8), parameter :: error = 1e-8_r8
  real(r8), dimension(n_bnd) :: tht_tmp, fr_tmp
  real(r8), dimension(3 * n_bnd + 6) :: work
  real(r8), dimension(n_bnd + 2) :: tht_sort, fr_sort, dfr_sort
  integer(itm_i4), dimension(n_bnd + 2) :: index_order
  real(r8) :: reast, rwest
  real(r8) :: tht, values(4)
  integer(itm_i4) :: ieast(1), iwest(1), n_bnd_short
  integer(itm_i4) :: i

  if (verbosity > 2) &
   write(iu6, *) 'fshape : (R, Z) set given on ', n_bnd, ' points'

  reast = maxval(r_bnd(1 : n_bnd))
  rwest = minval(r_bnd(1 : n_bnd))
  ieast = maxloc(r_bnd(1 : n_bnd))
  iwest = minloc(r_bnd(1 : n_bnd))
  rgeo = (reast + rwest) / 2._r8
  zgeo = (z_bnd(ieast(1)) + z_bnd(iwest(1))) / 2._r8
  amin = (reast - rwest) / 2._r8

  if (verbosity > 3) &
   write(iu6, '(A, 3f12.8)') 'Rgeo, Zgeo, a : ', rgeo, zgeo, amin

!-- geometric axis defined by geometric center of boundary

  do i = 1, n_bnd
    tht_tmp(i) = atan2(z_bnd(i) - zgeo, r_bnd(i) - rgeo)
    fr_tmp(i) = sqrt((r_bnd(i) - rgeo)**2 + (z_bnd(i) - zgeo)**2)
  end do

  if (abs(tht_tmp(n_bnd) - tht_tmp(1)) < error) n_bnd = n_bnd - 1

  do i = 1, n_bnd
    index_order(i) = i
  end do

  call qsortc(tht_tmp(1 : n_bnd), index_order(1 : n_bnd))

  tht_sort(1 : n_bnd) = tht_tmp(1 : n_bnd)
  do i = 1, n_bnd
    fr_sort(i) = fr_tmp(index_order(i))
  end do

  n_bnd_short = n_bnd

  do i = 2, n_bnd

    if ((tht_sort(i) - tht_sort(1)) > twopi) then
      n_bnd_short = i - 1
      exit
    end if

  end do

  if (abs(tht_sort(n_bnd_short) - tht_sort(1)) < error) then
    tht_sort(n_bnd_short) = tht_sort(1) + twopi
    fr_sort(n_bnd_short) = fr_sort(1)
  else
    tht_sort(n_bnd_short + 1) = tht_sort(1) + twopi
    fr_sort(n_bnd_short + 1) = fr_sort(1)
    n_bnd_short = n_bnd_short + 1
  end if

  call tb15a(n_bnd_short, n_bnd, tht_sort, fr_sort, dfr_sort, work, iu6)

  do i = 1, mfm
    tht = twopi * float(i - 1) / float(mfm)
    if (tht < tht_sort(1)) tht = tht + twopi
    if (tht > tht_sort(n_bnd_short)) tht = tht - twopi

    call tg02a(0, n_bnd_short, n_bnd, tht_sort, fr_sort, dfr_sort, tht, &
     values)

    fr(i) = values(1)

  end do  

  return

end subroutine shape_from_points