triangularity.f90

Go to the documentation of this file.

subroutine triangularity(xaxis, xell, xtria)
!----------------------------------------------------------------------
! subroutine to calculate the triangularity of flux surfaces as a function
! of the minor radius (not called by Helena)
!----------------------------------------------------------------------

  use itm_types
  use mod_mesh
  use mod_helena_io, only : out_he
  use mod_interpolation
  use mod_output
  use math_functions

  implicit none

  real(r8), intent(in) :: xaxis

  real(r8), dimension(nr), intent(out) :: xell, xtria

  real(r8), dimension(4 * (np - 1)) :: rrt
  real(r8), dimension(8 * (np - 1)) :: rr
  real(r8) :: r, s
  real(r8) :: x, y, ps
  real(r8) :: xrad, xshift, xgeo
  real(r8) :: xleft, xright
  real(r8) :: xtop, ytop, dummy
  real(r8) :: ell, tri
  real(r8), parameter :: tol = 1._r8 + 1.e-8_r8
  integer(itm_i4) :: i, j, ngs
  integer(itm_i4) :: n1, n2, n3, n4
  integer(itm_i4) :: npt
  integer(itm_i4) :: ifail

  if (standard_output) then
    write(out_he, *) '**************************************************'
    write(out_he, *) '*  ellipticity and triangularity   (fourier co.) *' 
    write(out_he, *) '**************************************************'
    write(out_he, *) '* index,  s,    radius,   shift,   ellip,   tria *'
    write(out_he, *) '**************************************************'
  end if
  do i = 1, nr - 1
    do j = 1, np - 1
      n1 = (i - 1) * np + j
      n2 = n1 + 1
      n3 = n2 + np
      n4 = n1 + np
      do ngs = 1, 4
        r = -1._r8
        s = -1._r8 + 0.5_r8 * dble(ngs - 1) 
        call interpolation(0, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
         r, s, x)
        call interpolation(0, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
         r, s, y)
        call interpolation(0, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
         psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps)
        rrt(4 * (j - 1) + ngs) = sqrt((x - xaxis)**2 + y * y)
      end do
    end do
    rrt(4 * (np - 1) + 1) = sqrt((xx(1, i * np) - xaxis)**2)
    do j = 1, 4 * (np - 1) + 1
      rr(j) = rrt(4 * (np - 1) - j + 2)
    end do
    do j = 1, 4 * (np - 1) - 1
      rr(4 * (np - 1) + j + 1) = rr(4 * (np - 1) - j + 1)
    end do
    npt = 8 * (np - 1)
    call rft2(rr, npt, 1)
    xrad = rr(1) / dble(npt)
    xshift = -2._r8 * rr(3) / dble(npt) / sqrt(ps)
    xell(i) = -4._r8 * rr(5) / dble(npt) / sqrt(ps)
    xtria(i) = 8._r8 * rr(7) / dble(npt) / ps
    if (standard_output) &
     write(out_he, 11) i, sqrt(ps), xrad, xshift, xell(i), xtria(i)
  end do

!-----------------------------------------------------------------
! alternative to calculate the ellipticity and triangularity by
! finding the point dy/ds=0 on each flux surface.
!-----------------------------------------------------------------
  if (standard_output) then
    write(out_he, *) '**************************************************'
    write(out_he, *) '*  ellipticity and triangularity (geometric co.) *' 
    write(out_he, *) '**************************************************'
  end if
  do i = 1, nr - 1
    do j = 1, np - 1
      n1 = (i - 1) * np + j
      n2 = (i - 1) * np + j + 1
      if  (yy(3, n1) * yy(3, n2) <= 0._r8) then
!------------------------------------- quad. eq for s value at minimum -
        call quad_equation(yy(1, n1), yy(3, n1), yy(1, n2), yy(3, n2), &
         tol, s, ifail)
        call cubic_interp_1d(xx(1, n1), xx(3, n1), xx(1, n2), &
         xx(3, n2), s, xtop, dummy)
        call cubic_interp_1d(yy(1, n1), yy(3, n1), yy(1, n2), &
         yy(3, n2), s, ytop, dummy)
        xleft = xx(1, (i - 1) * np + 1)
        xright = xx(1, (i - 1) * np + np)
        xgeo = (xleft + xright) / 2._r8
        xrad = (xright - xleft) / 2._r8
        xshift = (xaxis - xgeo) / xrad
        ell = ytop / xrad - 1._r8
        tri = -(xtop - xgeo) / xrad**2
        ps = psi(4 * (n1 - 1) + 1)
        if (standard_output) &
         write(out_he, 11) i, sqrt(ps), xrad, xshift, ell, tri
      end if
    end do
  end do

 11 format(i3, 5e16.8)

  return
end subroutine triangularity