program noncircular_duct
  implicit none
  integer :: gtype,ftype,hcfg,i
  double precision :: W,H,Di,Do,V,Tin,Tw,Lp
  double precision :: rho,mu,kf,Pr,cp
  double precision :: Dh,Ac,Pw,ar,Re,Nu,Nu_lam,f_d,h_c,Q,Tout,Le,mdot
  double precision :: Vs,Res,Nus,hs,fs
  double precision, parameter :: pi=3.14159265358979d0
  read(*,*) gtype; read(*,*) W; read(*,*) H
  read(*,*) Di; read(*,*) Do
  read(*,*) V; read(*,*) Tin; read(*,*) Tw; read(*,*) Lp
  read(*,*) ftype; read(*,*) rho; read(*,*) mu
  read(*,*) kf; read(*,*) Pr; read(*,*) cp; read(*,*) hcfg
  if(ftype==1) then; rho=1.177d0;mu=1.85d-5;kf=0.0263d0;Pr=0.71d0;cp=1007d0
  elseif(ftype==2) then; rho=997d0;mu=8.9d-4;kf=0.613d0;Pr=6.13d0;cp=4180d0
  elseif(ftype==3) then; rho=870d0;mu=0.05d0;kf=0.14d0;Pr=500d0;cp=2000d0
  endif
  W=W/1000d0; H=H/1000d0; Di=Di/1000d0; Do=Do/1000d0
  select case(gtype)
  case(1)
    Ac=W*H; Pw=2d0*(W+H); Dh=4d0*Ac/Pw
    ar=max(W,H)/max(min(W,H),1d-6)
    if(hcfg==1) then
      if(ar<=1.5d0) then; Nu_lam=2.98d0
      elseif(ar<=2.5d0) then; Nu_lam=3.39d0
      elseif(ar<=3.5d0) then; Nu_lam=3.96d0
      elseif(ar<=6d0) then; Nu_lam=4.44d0
      elseif(ar<=10d0) then; Nu_lam=5.60d0
      else; Nu_lam=7.54d0; endif
    else
      if(ar<=1.5d0) then; Nu_lam=3.61d0
      elseif(ar<=2.5d0) then; Nu_lam=4.12d0
      elseif(ar<=3.5d0) then; Nu_lam=4.79d0
      elseif(ar<=6d0) then; Nu_lam=5.33d0
      elseif(ar<=10d0) then; Nu_lam=6.49d0
      else; Nu_lam=8.24d0; endif
    endif
  case(2)
    Ac=sqrt(3d0)/4d0*W*W; Pw=3d0*W; Dh=4d0*Ac/Pw; ar=1d0
    if(hcfg==1) then; Nu_lam=2.47d0; else; Nu_lam=3.11d0; endif
  case default
    Dh=Do-Di; Ac=pi/4d0*(Do*Do-Di*Di); Pw=pi*(Do+Di)
    ar=Di/max(Do,1d-6)
    if(hcfg==1) then
      if(ar<0.08d0) then; Nu_lam=17.46d0
      elseif(ar<0.15d0) then; Nu_lam=11.56d0
      elseif(ar<0.35d0) then; Nu_lam=7.37d0
      elseif(ar<0.7d0) then; Nu_lam=5.74d0
      else; Nu_lam=4.86d0; endif
    else
      if(ar<0.08d0) then; Nu_lam=0d0
      elseif(ar<0.15d0) then; Nu_lam=11.91d0
      elseif(ar<0.35d0) then; Nu_lam=7.54d0
      elseif(ar<0.7d0) then; Nu_lam=6.14d0
      else; Nu_lam=5.39d0; endif
    endif
  end select
  Re=rho*V*Dh/mu
  Nu=Nu_lam
  f_d=64d0/max(Re,1d0)
  if(Re>2300d0) then
    f_d=(0.790d0*log(Re)-1.64d0)**(-2)
    Nu=(f_d/8d0)*(Re-1000d0)*Pr/(1d0+12.7d0*sqrt(f_d/8d0)*(Pr**(2d0/3d0)-1d0))
  endif
  h_c=Nu*kf/Dh
  mdot=rho*V*Ac
  Q=h_c*Pw*Lp*abs(Tw-Tin)
  if(mdot*cp>0d0) then; Tout=Tin+Q/(mdot*cp); else; Tout=Tin; endif
  if(Re<2300d0) then; Le=0.05d0*Re*Dh; else; Le=10d0*Dh; endif
  write(*,'(A)') '============================================'
  write(*,'(A)') '  NON-CIRCULAR DUCT & ANNULUS CONVECTION'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  if(gtype==1) write(*,'(A)') '  Geometry                = Rectangular'
  if(gtype==2) write(*,'(A)') '  Geometry                = Equilateral Triangle'
  if(gtype==3) write(*,'(A)') '  Geometry                = Annulus'
  if(hcfg==1) write(*,'(A)') '  Boundary Condition      = Uniform Wall Temp (T)'
  if(hcfg==2) write(*,'(A)') '  Boundary Condition      = Uniform Heat Flux (H)'
  write(*,'(A,F10.4,A)') '  Velocity V              = ',V,' m/s'
  write(*,'(A,F10.2,A)') '  Inlet Temp T_in         = ',Tin,' C'
  write(*,'(A,F10.2,A)') '  Wall Temp T_wall        = ',Tw,' C'
  write(*,'(A,F10.2,A)') '  Pipe Length L            = ',Lp,' m'
  write(*,'(A)') ''
  write(*,'(A)') '--- GEOMETRY ---'
  write(*,'(A,F12.6,A)') '  Hydraulic Diameter D_h  = ',Dh,' m'
  write(*,'(A,ES12.4,A)') '  Cross Section Area A    = ',Ac,' m2'
  write(*,'(A,F12.6,A)') '  Wetted Perimeter P      = ',Pw,' m'
  write(*,'(A,F10.3)')   '  Aspect / Radius Ratio   = ',ar
  write(*,'(A)') ''
  write(*,'(A)') '--- FLUID PROPERTIES ---'
  write(*,'(A,F12.4,A)') '  Density rho             = ',rho,' kg/m3'
  write(*,'(A,ES12.4,A)') '  Viscosity mu            = ',mu,' Pa.s'
  write(*,'(A,F12.6,A)') '  Conductivity k          = ',kf,' W/mK'
  write(*,'(A,F12.4)')    '  Prandtl Pr              = ',Pr
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,ES14.4)')  '  Reynolds Number Re      = ',Re
  if(Re<2300d0) then
    write(*,'(A)') '  Flow Regime             = Laminar'
  else
    write(*,'(A)') '  Flow Regime             = Turbulent'
  endif
  write(*,'(A,F10.6)')   '  Friction Factor f       = ',f_d
  write(*,'(A,F12.2)')   '  Nusselt Number Nu       = ',Nu
  write(*,'(A,F12.4,A)') '  Convection Coeff h      = ',h_c,' W/m2K'
  write(*,'(A,F12.2,A)') '  Heat Transfer Q         = ',Q,' W'
  write(*,'(A,F10.2,A)') '  Outlet Temperature T_out= ',Tout,' C'
  write(*,'(A,F10.4,A)') '  Entry Length Le          = ',Le,' m'
  write(*,'(A,ES12.4,A)') '  Mass Flow Rate          = ',mdot,' kg/s'
  write(*,'(A)') ''
  write(*,'(A)') '--- VELOCITY SWEEP ---'
  write(*,'(A)') '  V[m/s]     Re           Nu        h[W/m2K]   Regime'
  write(*,'(A)') '  --------------------------------------------------------'
  do i=1,25
    Vs=0.1d0+(V*4d0-0.1d0)*dble(i-1)/24d0
    Res=rho*Vs*Dh/mu
    if(Res>2300d0) then
      fs=(0.790d0*log(Res)-1.64d0)**(-2)
      Nus=(fs/8d0)*(Res-1000d0)*Pr/(1d0+12.7d0*sqrt(fs/8d0)*(Pr**(2d0/3d0)-1d0))
    else
      Nus=Nu_lam
    endif
    hs=Nus*kf/Dh
    if(Res>2300d0) then
      write(*,'(2X,F8.3,2X,ES10.3,2X,F9.2,2X,F10.3,2X,A)') Vs,Res,Nus,hs,'Turbulent'
    else
      write(*,'(2X,F8.3,2X,ES10.3,2X,F9.2,2X,F10.3,2X,A)') Vs,Res,Nus,hs,'Laminar'
    endif
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Laminar: Nu from geometry tables (Incropera Table 8.1)'
  write(*,'(A)') '  Turbulent: Gnielinski Nu=(f/8)(Re-1000)Pr/[1+12.7sqrt(f/8)(Pr^2/3-1)]'
  write(*,'(A)') '  f = (0.790 ln Re - 1.64)^-2'
  write(*,'(A)') '  Ref: Incropera Ch.8 Sec.8.6, Kothandaraman Ch.9 Sec.9.6'
end program noncircular_duct
