program regenerator
  implicit none
  integer :: i,rtype
  double precision :: Ch,Cc,Thi,Tci,NTU,Cr_star,rpm,m_mat,cp_mat
  double precision :: Cmin,Cmax,Cr,eps,Q,Tho,Tco,eps_cf
  double precision :: leak,NTUs,epss,Qs
  read(*,*) rtype; read(*,*) Ch; read(*,*) Cc
  read(*,*) Thi; read(*,*) Tci; read(*,*) NTU; read(*,*) Cr_star
  read(*,*) rpm; read(*,*) m_mat; read(*,*) cp_mat
  if(Ch<Cc) then; Cmin=Ch; Cmax=Cc; else; Cmin=Cc; Cmax=Ch; endif
  Cr=Cmin/Cmax
  if(abs(Cr-1d0)<1d-6) then
    eps_cf=NTU/(1d0+NTU)
  else
    eps_cf=(1d0-exp(-NTU*(1d0-Cr)))/(1d0-Cr*exp(-NTU*(1d0-Cr)))
  endif
  if(rtype==1.and.Cr_star>0d0) then
    eps=eps_cf*(1d0-1d0/(9d0*Cr_star**1.93d0))
    if(eps<0d0) eps=0d0
  else
    eps=eps_cf
  endif
  Q=eps*Cmin*(Thi-Tci)
  Tho=Thi-Q/Ch
  Tco=Tci+Q/Cc
  if(rtype==1) then; leak=0.03d0; else; leak=0d0; endif
  write(*,'(A)') '============================================'
  write(*,'(A)') '  REGENERATIVE HEAT EXCHANGER'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  if(rtype==1) write(*,'(A)') '  Type                    = Rotary (wheel)'
  if(rtype==2) write(*,'(A)') '  Type                    = Fixed matrix'
  write(*,'(A,F10.2,A)') '  C_hot                   = ',Ch,' W/K'
  write(*,'(A,F10.2,A)') '  C_cold                  = ',Cc,' W/K'
  write(*,'(A,F10.2,A)') '  T_hot_in                = ',Thi,' C'
  write(*,'(A,F10.2,A)') '  T_cold_in               = ',Tci,' C'
  write(*,'(A,F10.4)')    '  NTU                     = ',NTU
  write(*,'(A,F10.4)')    '  Cr* (matrix cap ratio)  = ',Cr_star
  if(rtype==1) write(*,'(A,F10.2,A)') '  Rotation speed          = ',rpm,' rpm'
  write(*,'(A,F10.2,A)') '  Matrix mass             = ',m_mat,' kg'
  write(*,'(A,F10.1,A)') '  Matrix cp               = ',cp_mat,' J/kgK'
  write(*,'(A)') ''
  write(*,'(A)') '--- RESULTS ---'
  write(*,'(A,F10.4)')    '  C_min                   = ',Cmin
  write(*,'(A,F10.4)')    '  C_max                   = ',Cmax
  write(*,'(A,F10.4)')    '  C_r = Cmin/Cmax         = ',Cr
  write(*,'(A,F10.4)')    '  eps_counterflow          = ',eps_cf
  write(*,'(A,F10.4)')    '  eps (with Cr* corr)     = ',eps
  write(*,'(A,F12.2,A)') '  Heat transfer Q         = ',Q,' W'
  write(*,'(A,F10.2,A)') '  T_hot_out               = ',Tho,' C'
  write(*,'(A,F10.2,A)') '  T_cold_out              = ',Tco,' C'
  if(rtype==1) write(*,'(A,F8.2,A)') '  Carryover leakage est   = ',leak*100d0,' %'
  write(*,'(A)') ''
  write(*,'(A)') '--- NTU SWEEP ---'
  write(*,'(A)') '  NTU     eps_cf    eps       Q[W]         Tho[C]   Tco[C]'
  write(*,'(A)') '  ---------------------------------------------------------------'
  do i=1,25
    NTUs=0.5d0+9.5d0*dble(i-1)/24d0
    if(abs(Cr-1d0)<1d-6) then
      epss=NTUs/(1d0+NTUs)
    else
      epss=(1d0-exp(-NTUs*(1d0-Cr)))/(1d0-Cr*exp(-NTUs*(1d0-Cr)))
    endif
    if(rtype==1.and.Cr_star>0d0) epss=epss*(1d0-1d0/(9d0*Cr_star**1.93d0))
    if(epss<0d0) epss=0d0
    Qs=epss*Cmin*(Thi-Tci)
    write(*,'(2X,F6.2,2X,F8.4,2X,F8.4,2X,F12.2,2X,F8.2,2X,F8.2)') NTUs,NTUs/(1d0+NTUs),epss,Qs,Thi-Qs/Ch,Tci+Qs/Cc
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Counterflow: eps=(1-exp(-NTU(1-Cr)))/(1-Cr*exp(-NTU(1-Cr)))'
  write(*,'(A)') '  Rotary Cr* correction: eps=eps_cf*(1-1/(9*Cr*^1.93))'
  write(*,'(A)') '  Valid for Cr* > 2, NTU up to ~10'
  write(*,'(A)') '  Ref: Kays & London, Compact Heat Exchangers (1984)'
  write(*,'(A)') '       Shah & Sekulic, Fundamentals of HX Design (2003)'
end program regenerator