!==============================================================================
! ThermoFluidCalc — Calculator #27 : Mesh Aspect Ratio
!==============================================================================
! Physics : The aspect ratio (AR) measures how stretched / elongated a
!           mesh cell is.  Ideal AR = 1 (equilateral tri or square quad).
!
!   Triangle AR (inradius-based):
!     AR = max_edge / (2 * sqrt(3) * r)
!     where r = 2*Area / perimeter   (inradius)
!     For an equilateral triangle AR = 1.
!
!   Quad AR (opposite-edge-pair):
!     avg1 = (|e1| + |e3|) / 2    (pair of opposite edges)
!     avg2 = (|e2| + |e4|) / 2
!     AR   = max(avg1, avg2) / min(avg1, avg2)
!
!   Simple edge-ratio (both):
!     AR_edge = longest_edge / shortest_edge
!
! Quality rating:
!   1.0 - 1.5   Excellent
!   1.5 - 3.0   Good
!   3.0 - 5.0   Acceptable
!   5.0 - 10.0  Poor
!   > 10.0      Bad
!
! Reference : Gupta, §5.1.1
!
! Modes:
!   1 = Single triangle  (3 vertices 2-D)
!   2 = Single quad       (4 vertices 2-D)
!   3 = Batch             (N mixed cells)
!
! Build:
!   gfortran -O2 -o aspect_ratio aspect_ratio.f90
!==============================================================================
program aspect_ratio
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: SQRT3 = 1.7320508075688772_dp
  integer, parameter :: MAX_CELLS = 10000

  integer  :: mode, N, i, nv
  real(dp) :: vx(4), vy(4)
  real(dp) :: edges(4), perimeter, area
  real(dp) :: emax, emin, AR_edge, AR_main
  real(dp) :: inrad, circumR
  real(dp) :: avg1, avg2, d1, d2, diag_ratio
  character(len=20) :: rating

  ! Batch arrays
  real(dp) :: ar_arr(MAX_CELLS)
  real(dp) :: avg_ar, min_ar, max_ar, std_ar, 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)

    ! Edge lengths
    edges(1) = sqrt((vx(2)-vx(1))**2 + (vy(2)-vy(1))**2)
    edges(2) = sqrt((vx(3)-vx(2))**2 + (vy(3)-vy(2))**2)
    edges(3) = sqrt((vx(1)-vx(3))**2 + (vy(1)-vy(3))**2)

    perimeter = edges(1) + edges(2) + edges(3)
    area = 0.5_dp * abs((vx(2)-vx(1))*(vy(3)-vy(1)) - (vx(3)-vx(1))*(vy(2)-vy(1)))

    emax = max(edges(1), edges(2), edges(3))
    emin = min(edges(1), edges(2), edges(3))
    AR_edge = emax / max(emin, 1.0e-30_dp)

    ! Inradius
    inrad = 2.0_dp * area / max(perimeter, 1.0e-30_dp)

    ! Circumradius
    if (area > 0.0_dp) then
      circumR = (edges(1)*edges(2)*edges(3)) / (4.0_dp * area)
    else
      circumR = 0.0_dp
    end if

    ! Inradius-based AR (equilateral = 1)
    if (inrad > 0.0_dp) then
      AR_main = emax / (2.0_dp * SQRT3 * inrad)
    else
      AR_main = 1.0e6_dp
    end if

    call get_rating(AR_main, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Triangle'
    write(*,'(A,ES15.8)')  'EDGE_1=', edges(1)
    write(*,'(A,ES15.8)')  'EDGE_2=', edges(2)
    write(*,'(A,ES15.8)')  'EDGE_3=', edges(3)
    write(*,'(A,ES15.8)')  'PERIMETER=', perimeter
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,ES15.8)')  'INRADIUS=', inrad
    write(*,'(A,ES15.8)')  'CIRCUMRADIUS=', circumR
    write(*,'(A,F12.6)')   'EMAX=', emax
    write(*,'(A,F12.6)')   'EMIN=', emin
    write(*,'(A,F12.6)')   'AR_MAIN=', AR_main
    write(*,'(A,F12.6)')   'AR_EDGE=', AR_edge
    write(*,'(A,A)')       'RATING=', trim(rating)

    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)

    ! Edge lengths (sequential)
    do i = 1, 4
      edges(i) = sqrt((vx(mod(i,4)+1)-vx(i))**2 + (vy(mod(i,4)+1)-vy(i))**2)
    end do

    emax = max(edges(1), edges(2), edges(3), edges(4))
    emin = min(edges(1), edges(2), edges(3), edges(4))
    AR_edge = emax / max(emin, 1.0e-30_dp)

    ! Opposite-edge-pair AR
    avg1 = (edges(1) + edges(3)) / 2.0_dp
    avg2 = (edges(2) + edges(4)) / 2.0_dp
    AR_main = max(avg1, avg2) / max(min(avg1, avg2), 1.0e-30_dp)

    ! Diagonals
    d1 = sqrt((vx(3)-vx(1))**2 + (vy(3)-vy(1))**2)
    d2 = sqrt((vx(4)-vx(2))**2 + (vy(4)-vy(2))**2)
    diag_ratio = max(d1,d2) / max(min(d1,d2), 1.0e-30_dp)

    ! Area (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)) )

    perimeter = edges(1)+edges(2)+edges(3)+edges(4)

    call get_rating(AR_main, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Quad'
    write(*,'(A,ES15.8)')  'EDGE_1=', edges(1)
    write(*,'(A,ES15.8)')  'EDGE_2=', edges(2)
    write(*,'(A,ES15.8)')  'EDGE_3=', edges(3)
    write(*,'(A,ES15.8)')  'EDGE_4=', edges(4)
    write(*,'(A,ES15.8)')  'PERIMETER=', perimeter
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,F12.6)')   'AVG_PAIR1=', avg1
    write(*,'(A,F12.6)')   'AVG_PAIR2=', avg2
    write(*,'(A,F12.6)')   'DIAG_1=', d1
    write(*,'(A,F12.6)')   'DIAG_2=', d2
    write(*,'(A,F12.6)')   'DIAG_RATIO=', diag_ratio
    write(*,'(A,F12.6)')   'EMAX=', emax
    write(*,'(A,F12.6)')   'EMIN=', emin
    write(*,'(A,F12.6)')   'AR_MAIN=', AR_main
    write(*,'(A,F12.6)')   'AR_EDGE=', AR_edge
    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)
        edges(1) = sqrt((vx(2)-vx(1))**2+(vy(2)-vy(1))**2)
        edges(2) = sqrt((vx(3)-vx(2))**2+(vy(3)-vy(2))**2)
        edges(3) = sqrt((vx(1)-vx(3))**2+(vy(1)-vy(3))**2)
        perimeter = edges(1)+edges(2)+edges(3)
        area = 0.5_dp*abs((vx(2)-vx(1))*(vy(3)-vy(1))-(vx(3)-vx(1))*(vy(2)-vy(1)))
        emax = max(edges(1),edges(2),edges(3))
        inrad = 2.0_dp*area/max(perimeter,1.0e-30_dp)
        if (inrad > 0.0_dp) then
          AR_main = emax/(2.0_dp*SQRT3*inrad)
        else
          AR_main = 1.0e6_dp
        end if

      else if (nv == 4) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
        do nv = 1, 4
          edges(nv) = sqrt((vx(mod(nv,4)+1)-vx(nv))**2+(vy(mod(nv,4)+1)-vy(nv))**2)
        end do
        avg1 = (edges(1)+edges(3))/2.0_dp
        avg2 = (edges(2)+edges(4))/2.0_dp
        AR_main = max(avg1,avg2)/max(min(avg1,avg2),1.0e-30_dp)

      else
        read(*,*)
        AR_main = 1.0e6_dp
      end if

      ar_arr(i) = AR_main

      if (AR_main < 1.5_dp) then;      hist(1)=hist(1)+1
      else if (AR_main < 3.0_dp) then;  hist(2)=hist(2)+1
      else if (AR_main < 5.0_dp) then;  hist(3)=hist(3)+1
      else if (AR_main < 10.0_dp) then; hist(4)=hist(4)+1
      else;                              hist(5)=hist(5)+1
      end if
    end do

    ! Stats
    s1=0; s2=0; min_ar=ar_arr(1); max_ar=ar_arr(1)
    do i=1,N
      s1=s1+ar_arr(i); s2=s2+ar_arr(i)**2
      if(ar_arr(i)<min_ar) min_ar=ar_arr(i)
      if(ar_arr(i)>max_ar) max_ar=ar_arr(i)
    end do
    avg_ar=s1/real(N,dp)
    std_ar=sqrt(max(0.0_dp, s2/real(N,dp)-avg_ar**2))

    write(*,'(A,I1)')    'MODE=', mode
    write(*,'(A)')       'MODE_NAME=Batch'
    write(*,'(A,I6)')    'NCELLS=', N
    write(*,'(A,F12.6)') 'AVG_AR=', avg_ar
    write(*,'(A,F12.6)') 'MIN_AR=', min_ar
    write(*,'(A,F12.6)') 'MAX_AR=', max_ar
    write(*,'(A,F12.6)') 'STD_AR=', std_ar

    write(*,'(A)') 'HIST_START'
    write(*,'(A,I6)') '1.0-1.5,', hist(1)
    write(*,'(A,I6)') '1.5-3.0,', hist(2)
    write(*,'(A,I6)') '3.0-5.0,', hist(3)
    write(*,'(A,I6)') '5.0-10,',  hist(4)
    write(*,'(A,I6)') '>10,',     hist(5)
    write(*,'(A)') 'HIST_END'

    write(*,'(A)') 'DATA_START'
    do i=1,N
      call get_rating(ar_arr(i), rating)
      write(*,'(I6,A,F12.6,A,A)') i, ',', ar_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 get_rating(ar, r)
    real(dp), intent(in) :: ar
    character(len=20), intent(out) :: r
    if (ar < 1.5_dp) then;       r = 'Excellent'
    else if (ar < 3.0_dp) then;  r = 'Good'
    else if (ar < 5.0_dp) then;  r = 'Acceptable'
    else if (ar < 10.0_dp) then; r = 'Poor'
    else;                         r = 'Bad'
    end if
  end subroutine

end program aspect_ratio
