!==============================================================================
! ThermoFluidCalc — Calculator #26 : Mesh Skewness
!==============================================================================
! Physics : Cell skewness measures how far an element deviates from the
!           ideal (equilateral / equiangular) shape.
!
!   Angular skewness (any 2-D cell):
!     skew = max( (theta_max - theta_eq)/(180 - theta_eq) ,
!                 (theta_eq - theta_min)/theta_eq )
!     theta_eq = (n-2)*180/n   (60 for tri, 90 for quad)
!
!   Equivolume skewness (triangle):
!     skew_vol = (A_opt - A_cell) / A_opt
!     A_opt = (3*sqrt(3)/4) * R^2   (equilateral tri with same circumradius R)
!
! Quality rating:
!   0.00 - 0.25  Excellent
!   0.25 - 0.50  Good
!   0.50 - 0.75  Acceptable
!   0.75 - 0.90  Poor
!   0.90 - 0.97  Bad
!   0.97 - 1.00  Degenerate
!
! Reference : Gupta, §5.1.1
!
! Modes:
!   1 = Single triangle (2-D, 3 vertices)
!   2 = Single quad      (2-D, 4 vertices)
!   3 = Batch            (N mixed cells)
!
! Build:
!   gfortran -O2 -o skewness skewness.f90
!==============================================================================
program skewness
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: PI = 3.141592653589793238_dp
  real(dp), parameter :: DEG = 180.0_dp / PI
  integer, parameter :: MAX_CELLS = 10000

  integer  :: mode, N, i, nv
  real(dp) :: vx(4), vy(4)
  real(dp) :: angles(4), theta_min, theta_max, theta_eq
  real(dp) :: skew_ang, skew_vol, area, circumR
  real(dp) :: a_side, b_side, c_side, A_opt
  character(len=20) :: rating

  ! Batch arrays
  real(dp) :: skew_arr(MAX_CELLS)
  real(dp) :: avg_s, min_s, max_s, std_s, sum_s, sum2
  integer  :: hist(5)  ! bins: [0,.25) [.25,.5) [.5,.75) [.75,.9) [.9,1]

  read(*,*) mode

  select case (mode)

  !=========================================================================
  ! MODE 1 : Single triangle
  !=========================================================================
  case (1)
    read(*,*) vx(1),vy(1), vx(2),vy(2), vx(3),vy(3)
    nv = 3

    call compute_angles_tri(vx, vy, angles, theta_min, theta_max)
    theta_eq = 60.0_dp
    skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
                   (theta_eq - theta_min)/theta_eq)
    skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))

    ! Area
    area = 0.5_dp * abs((vx(2)-vx(1))*(vy(3)-vy(1)) - (vx(3)-vx(1))*(vy(2)-vy(1)))

    ! Side lengths
    a_side = sqrt((vx(2)-vx(1))**2 + (vy(2)-vy(1))**2)
    b_side = sqrt((vx(3)-vx(2))**2 + (vy(3)-vy(2))**2)
    c_side = sqrt((vx(1)-vx(3))**2 + (vy(1)-vy(3))**2)

    ! Circumradius R = abc / (4*Area)
    if (area > 0.0_dp) then
      circumR = (a_side * b_side * c_side) / (4.0_dp * area)
    else
      circumR = 0.0_dp
    end if

    ! Equivolume skewness
    A_opt = (3.0_dp * sqrt(3.0_dp) / 4.0_dp) * circumR**2
    if (A_opt > 0.0_dp) then
      skew_vol = (A_opt - area) / A_opt
    else
      skew_vol = 1.0_dp
    end if
    skew_vol = max(0.0_dp, min(skew_vol, 1.0_dp))

    call get_rating(skew_ang, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Triangle'
    write(*,'(A,F10.4)')   'ANGLE_1=', angles(1)
    write(*,'(A,F10.4)')   'ANGLE_2=', angles(2)
    write(*,'(A,F10.4)')   'ANGLE_3=', angles(3)
    write(*,'(A,F10.4)')   'THETA_MIN=', theta_min
    write(*,'(A,F10.4)')   'THETA_MAX=', theta_max
    write(*,'(A,F10.4)')   'THETA_EQ=', theta_eq
    write(*,'(A,F10.6)')   'SKEW_ANGULAR=', skew_ang
    write(*,'(A,F10.6)')   'SKEW_EQUIVOL=', skew_vol
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,ES15.8)')  'CIRCUMRADIUS=', circumR
    write(*,'(A,ES15.8)')  'A_OPT=', A_opt
    write(*,'(A,A)')       'RATING=', trim(rating)

    ! Vertex data for drawing
    write(*,'(A)') 'VERTS_START'
    do i = 1, 3
      write(*,'(F12.6,A,F12.6)') vx(i), ',', vy(i)
    end do
    write(*,'(A)') 'VERTS_END'

  !=========================================================================
  ! MODE 2 : Single quad
  !=========================================================================
  case (2)
    read(*,*) vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
    nv = 4

    call compute_angles_quad(vx, vy, angles, theta_min, theta_max)
    theta_eq = 90.0_dp
    skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
                   (theta_eq - theta_min)/theta_eq)
    skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))

    ! Area via shoelace
    area = 0.5_dp * abs( (vx(1)*vy(2) - vx(2)*vy(1)) + &
                          (vx(2)*vy(3) - vx(3)*vy(2)) + &
                          (vx(3)*vy(4) - vx(4)*vy(3)) + &
                          (vx(4)*vy(1) - vx(1)*vy(4)) )

    call get_rating(skew_ang, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Quad'
    write(*,'(A,F10.4)')   'ANGLE_1=', angles(1)
    write(*,'(A,F10.4)')   'ANGLE_2=', angles(2)
    write(*,'(A,F10.4)')   'ANGLE_3=', angles(3)
    write(*,'(A,F10.4)')   'ANGLE_4=', angles(4)
    write(*,'(A,F10.4)')   'THETA_MIN=', theta_min
    write(*,'(A,F10.4)')   'THETA_MAX=', theta_max
    write(*,'(A,F10.4)')   'THETA_EQ=', theta_eq
    write(*,'(A,F10.6)')   'SKEW_ANGULAR=', skew_ang
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,A)')       'RATING=', trim(rating)

    write(*,'(A)') 'VERTS_START'
    do i = 1, 4
      write(*,'(F12.6,A,F12.6)') vx(i), ',', vy(i)
    end do
    write(*,'(A)') 'VERTS_END'

  !=========================================================================
  ! MODE 3 : Batch
  !=========================================================================
  case (3)
    backspace(5)
    read(*,*) mode, N
    if (N < 1 .or. N > MAX_CELLS) then
      write(*,'(A)') 'ERROR=N must be 1-10000.'; stop
    end if

    hist = 0
    do i = 1, N
      read(*,*) nv
      backspace(5)
      if (nv == 3) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3)
        call compute_angles_tri(vx, vy, angles, theta_min, theta_max)
        theta_eq = 60.0_dp
      else if (nv == 4) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
        call compute_angles_quad(vx, vy, angles, theta_min, theta_max)
        theta_eq = 90.0_dp
      else
        read(*,*)  ! skip line
        skew_arr(i) = 1.0_dp
        cycle
      end if

      skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
                     (theta_eq - theta_min)/theta_eq)
      skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))
      skew_arr(i) = skew_ang

      ! Histogram bin
      if (skew_ang < 0.25_dp) then
        hist(1) = hist(1) + 1
      else if (skew_ang < 0.50_dp) then
        hist(2) = hist(2) + 1
      else if (skew_ang < 0.75_dp) then
        hist(3) = hist(3) + 1
      else if (skew_ang < 0.90_dp) then
        hist(4) = hist(4) + 1
      else
        hist(5) = hist(5) + 1
      end if
    end do

    ! Statistics
    sum_s = 0.0_dp; sum2 = 0.0_dp
    min_s = skew_arr(1); max_s = skew_arr(1)
    do i = 1, N
      sum_s = sum_s + skew_arr(i)
      sum2  = sum2  + skew_arr(i)**2
      if (skew_arr(i) < min_s) min_s = skew_arr(i)
      if (skew_arr(i) > max_s) max_s = skew_arr(i)
    end do
    avg_s = sum_s / real(N, dp)
    std_s = sqrt(max(0.0_dp, sum2/real(N,dp) - avg_s**2))

    write(*,'(A,I1)')    'MODE=', mode
    write(*,'(A)')       'MODE_NAME=Batch'
    write(*,'(A,I6)')    'NCELLS=', N
    write(*,'(A,F10.6)') 'AVG_SKEW=', avg_s
    write(*,'(A,F10.6)') 'MIN_SKEW=', min_s
    write(*,'(A,F10.6)') 'MAX_SKEW=', max_s
    write(*,'(A,F10.6)') 'STD_SKEW=', std_s

    write(*,'(A)') 'HIST_START'
    write(*,'(A,I6)') '0.00-0.25,', hist(1)
    write(*,'(A,I6)') '0.25-0.50,', hist(2)
    write(*,'(A,I6)') '0.50-0.75,', hist(3)
    write(*,'(A,I6)') '0.75-0.90,', hist(4)
    write(*,'(A,I6)') '0.90-1.00,', hist(5)
    write(*,'(A)') 'HIST_END'

    ! Per-cell data
    write(*,'(A)') 'DATA_START'
    do i = 1, N
      call get_rating(skew_arr(i), rating)
      write(*,'(I6,A,F10.6,A,A)') i, ',', skew_arr(i), ',', trim(rating)
    end do
    write(*,'(A)') 'DATA_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (must be 1-3).'; stop
  end select

contains

  !------------------------------------------------------------------------
  subroutine compute_angles_tri(x, y, ang, amin, amax)
    real(dp), intent(in)  :: x(4), y(4)
    real(dp), intent(out) :: ang(4), amin, amax
    real(dp) :: dx1,dy1,dx2,dy2,dot,cross,mag1,mag2
    integer :: j, j1, j2

    do j = 1, 3
      j1 = mod(j, 3) + 1      ! next vertex
      j2 = mod(j + 1, 3) + 1  ! prev vertex
      dx1 = x(j1) - x(j);  dy1 = y(j1) - y(j)
      dx2 = x(j2) - x(j);  dy2 = y(j2) - y(j)
      mag1 = sqrt(dx1**2 + dy1**2)
      mag2 = sqrt(dx2**2 + dy2**2)
      if (mag1 < 1.0e-30_dp .or. mag2 < 1.0e-30_dp) then
        ang(j) = 0.0_dp; cycle
      end if
      dot = dx1*dx2 + dy1*dy2
      ang(j) = acos(max(-1.0_dp, min(1.0_dp, dot/(mag1*mag2)))) * DEG
    end do
    ang(4) = 0.0_dp
    amin = min(ang(1), ang(2), ang(3))
    amax = max(ang(1), ang(2), ang(3))
  end subroutine

  !------------------------------------------------------------------------
  subroutine compute_angles_quad(x, y, ang, amin, amax)
    real(dp), intent(in)  :: x(4), y(4)
    real(dp), intent(out) :: ang(4), amin, amax
    real(dp) :: dx1,dy1,dx2,dy2,mag1,mag2,dot
    integer :: j, jp, jm

    do j = 1, 4
      jm = mod(j + 2, 4) + 1  ! previous vertex
      jp = mod(j, 4) + 1      ! next vertex
      dx1 = x(jm) - x(j);  dy1 = y(jm) - y(j)
      dx2 = x(jp) - x(j);  dy2 = y(jp) - y(j)
      mag1 = sqrt(dx1**2 + dy1**2)
      mag2 = sqrt(dx2**2 + dy2**2)
      if (mag1 < 1.0e-30_dp .or. mag2 < 1.0e-30_dp) then
        ang(j) = 0.0_dp; cycle
      end if
      dot = dx1*dx2 + dy1*dy2
      ang(j) = acos(max(-1.0_dp, min(1.0_dp, dot/(mag1*mag2)))) * DEG
    end do
    amin = min(ang(1), ang(2), ang(3), ang(4))
    amax = max(ang(1), ang(2), ang(3), ang(4))
  end subroutine

  !------------------------------------------------------------------------
  subroutine get_rating(s, r)
    real(dp), intent(in)       :: s
    character(len=20), intent(out) :: r
    if (s < 0.25_dp) then
      r = 'Excellent'
    else if (s < 0.50_dp) then
      r = 'Good'
    else if (s < 0.75_dp) then
      r = 'Acceptable'
    else if (s < 0.90_dp) then
      r = 'Poor'
    else if (s < 0.97_dp) then
      r = 'Bad'
    else
      r = 'Degenerate'
    end if
  end subroutine

end program skewness