!==============================================================================
! ThermoFluidCalc — Calculator #28 : Orthogonal Quality (Ortho-Skew)
!==============================================================================
! Physics : Orthogonal quality measures how well the vector from the cell
!           centroid to each face centroid aligns with the face outward
!           normal.  Perfect alignment (orthogonal grid) gives OQ = 1.
!
!   For each face f:
!     OQ_f = |dot(A_f, d_f)| / (|A_f| * |d_f|) = |cos(theta)|
!
!   where
!     A_f = outward face normal vector
!     d_f = vector from cell centroid to face midpoint
!
!   OQ_cell = min(OQ_f)   over all faces of the cell
!
! Quality rating (0 = worst, 1 = best):
!   0.95 - 1.00  Excellent
!   0.70 - 0.95  Good
!   0.50 - 0.70  Acceptable
!   0.25 - 0.50  Poor
!   0.00 - 0.25  Bad
!
! Reference : Gupta, §5.1.1
!
! Modes:
!   1 = Single triangle  (3 vertices)
!   2 = Single quad       (4 vertices)
!   3 = Batch             (N mixed cells)
!
! Build:
!   gfortran -O2 -o ortho_quality ortho_quality.f90
!==============================================================================
program ortho_quality
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: MAX_CELLS = 10000

  integer  :: mode, N, i, j, nv, jp
  real(dp) :: vx(4), vy(4)
  real(dp) :: cx, cy           ! cell centroid
  real(dp) :: fmx, fmy         ! face midpoint
  real(dp) :: nx_f, ny_f       ! face outward normal
  real(dp) :: dx_f, dy_f       ! centroid-to-face vector
  real(dp) :: mag_n, mag_d, dot_val, oq_f
  real(dp) :: oq_cell, oq_min_face
  real(dp) :: face_oq(4), face_nx(4), face_ny(4), face_mx(4), face_my(4)
  real(dp) :: face_dx(4), face_dy(4)
  character(len=20) :: rating

  ! Batch
  real(dp) :: oq_arr(MAX_CELLS)
  real(dp) :: avg_oq, min_oq, max_oq, std_oq, s1, s2
  integer  :: hist(5)

  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_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, &
                    face_mx, face_my, face_dx, face_dy, oq_cell)
    call get_rating(oq_cell, rating)

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Single Triangle'
    write(*,'(A,F12.6)')  'CENTROID_X=', cx
    write(*,'(A,F12.6)')  'CENTROID_Y=', cy
    write(*,'(A,F10.6)')  'OQ_CELL=', oq_cell
    write(*,'(A,A)')      'RATING=', trim(rating)
    write(*,'(A,I1)')     'NFACES=', nv

    write(*,'(A)') 'FACES_START'
    do j = 1, nv
      write(*,'(I2,5(A,F12.6))') j, &
        ',', face_nx(j), ',', face_ny(j), &
        ',', face_mx(j), ',', face_my(j), ',', face_oq(j)
    end do
    write(*,'(A)') 'FACES_END'

    write(*,'(A)') 'VERTS_START'
    do j = 1, nv
      write(*,'(F12.6,A,F12.6)') vx(j), ',', vy(j)
    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_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, &
                    face_mx, face_my, face_dx, face_dy, oq_cell)
    call get_rating(oq_cell, rating)

    write(*,'(A,I1)')     'MODE=', mode
    write(*,'(A)')        'MODE_NAME=Single Quad'
    write(*,'(A,F12.6)')  'CENTROID_X=', cx
    write(*,'(A,F12.6)')  'CENTROID_Y=', cy
    write(*,'(A,F10.6)')  'OQ_CELL=', oq_cell
    write(*,'(A,A)')      'RATING=', trim(rating)
    write(*,'(A,I1)')     'NFACES=', nv

    write(*,'(A)') 'FACES_START'
    do j = 1, nv
      write(*,'(I2,5(A,F12.6))') j, &
        ',', face_nx(j), ',', face_ny(j), &
        ',', face_mx(j), ',', face_my(j), ',', face_oq(j)
    end do
    write(*,'(A)') 'FACES_END'

    write(*,'(A)') 'VERTS_START'
    do j = 1, nv
      write(*,'(F12.6,A,F12.6)') vx(j), ',', vy(j)
    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)
      else if (nv == 4) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
      else
        read(*,*)
        oq_arr(i) = 0.0_dp
        hist(5) = hist(5) + 1
        cycle
      end if

      call compute_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, &
                      face_mx, face_my, face_dx, face_dy, oq_cell)
      oq_arr(i) = oq_cell

      if (oq_cell >= 0.95_dp) then;      hist(1)=hist(1)+1
      else if (oq_cell >= 0.70_dp) then;  hist(2)=hist(2)+1
      else if (oq_cell >= 0.50_dp) then;  hist(3)=hist(3)+1
      else if (oq_cell >= 0.25_dp) then;  hist(4)=hist(4)+1
      else;                                hist(5)=hist(5)+1
      end if
    end do

    s1=0; s2=0; min_oq=oq_arr(1); max_oq=oq_arr(1)
    do i=1,N
      s1=s1+oq_arr(i); s2=s2+oq_arr(i)**2
      if(oq_arr(i)<min_oq) min_oq=oq_arr(i)
      if(oq_arr(i)>max_oq) max_oq=oq_arr(i)
    end do
    avg_oq=s1/real(N,dp)
    std_oq=sqrt(max(0.0_dp, s2/real(N,dp)-avg_oq**2))

    write(*,'(A,I1)')    'MODE=', mode
    write(*,'(A)')       'MODE_NAME=Batch'
    write(*,'(A,I6)')    'NCELLS=', N
    write(*,'(A,F10.6)') 'AVG_OQ=', avg_oq
    write(*,'(A,F10.6)') 'MIN_OQ=', min_oq
    write(*,'(A,F10.6)') 'MAX_OQ=', max_oq
    write(*,'(A,F10.6)') 'STD_OQ=', std_oq

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

    write(*,'(A)') 'DATA_START'
    do i=1,N
      call get_rating(oq_arr(i), rating)
      write(*,'(I6,A,F10.6,A,A)') i, ',', oq_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

  !------------------------------------------------------------------------
  ! Compute orthogonal quality for a 2-D polygon cell (nv = 3 or 4)
  !------------------------------------------------------------------------
  subroutine compute_oq(vx, vy, nv, cx, cy, foq, fnx, fny, fmx, fmy, fdx, fdy, oq_min)
    real(dp), intent(in)  :: vx(4), vy(4)
    integer,  intent(in)  :: nv
    real(dp), intent(out) :: cx, cy, foq(4), fnx(4), fny(4)
    real(dp), intent(out) :: fmx(4), fmy(4), fdx(4), fdy(4), oq_min
    integer :: j, jp
    real(dp) :: ex, ey, mn, md, dv

    ! Centroid
    cx = 0.0_dp; cy = 0.0_dp
    do j = 1, nv
      cx = cx + vx(j); cy = cy + vy(j)
    end do
    cx = cx / real(nv, dp); cy = cy / real(nv, dp)

    oq_min = 1.0_dp

    do j = 1, nv
      jp = mod(j, nv) + 1

      ! Face midpoint
      fmx(j) = 0.5_dp * (vx(j) + vx(jp))
      fmy(j) = 0.5_dp * (vy(j) + vy(jp))

      ! Edge vector
      ex = vx(jp) - vx(j)
      ey = vy(jp) - vy(j)

      ! Outward normal (rotate edge 90° clockwise, then check orientation)
      fnx(j) =  ey
      fny(j) = -ex
      ! Ensure outward: normal should point away from centroid
      ! dot(normal, midpoint - centroid) > 0  for outward
      if (fnx(j)*(fmx(j)-cx) + fny(j)*(fmy(j)-cy) < 0.0_dp) then
        fnx(j) = -fnx(j)
        fny(j) = -fny(j)
      end if

      ! d = centroid to face midpoint
      fdx(j) = fmx(j) - cx
      fdy(j) = fmy(j) - cy

      ! OQ_f = |cos(theta)|
      mn = sqrt(fnx(j)**2 + fny(j)**2)
      md = sqrt(fdx(j)**2 + fdy(j)**2)
      if (mn > 1.0e-30_dp .and. md > 1.0e-30_dp) then
        dv = abs(fnx(j)*fdx(j) + fny(j)*fdy(j)) / (mn * md)
        foq(j) = min(dv, 1.0_dp)
      else
        foq(j) = 0.0_dp
      end if

      if (foq(j) < oq_min) oq_min = foq(j)
    end do

    ! Zero unused entries
    do j = nv+1, 4
      foq(j) = 0.0_dp; fnx(j) = 0.0_dp; fny(j) = 0.0_dp
      fmx(j) = 0.0_dp; fmy(j) = 0.0_dp; fdx(j) = 0.0_dp; fdy(j) = 0.0_dp
    end do
  end subroutine

  !------------------------------------------------------------------------
  subroutine get_rating(oq, r)
    real(dp), intent(in)       :: oq
    character(len=20), intent(out) :: r
    if (oq >= 0.95_dp) then;      r = 'Excellent'
    else if (oq >= 0.70_dp) then;  r = 'Good'
    else if (oq >= 0.50_dp) then;  r = 'Acceptable'
    else if (oq >= 0.25_dp) then;  r = 'Poor'
    else;                           r = 'Bad'
    end if
  end subroutine

end program ortho_quality
