program porous_media_cond
  implicit none
  integer :: i, j
  double precision :: ks,kf,eps,Th,Tc,Lth,A
  double precision :: dT,kpar,kser,kHSu,kHSl,kEMT,kgeo,keff
  double precision :: Q,Rth,es,f1,f2,df,dk
  read(*,*) ks; read(*,*) kf; read(*,*) eps
  read(*,*) Th; read(*,*) Tc; read(*,*) Lth; read(*,*) A
  dT=abs(Th-Tc)
  kpar=eps*kf+(1d0-eps)*ks
  kser=1d0/(eps/kf+(1d0-eps)/ks)
  if(abs(kf-ks)>1d-12) then
    kHSu=ks+eps/(1d0/(kf-ks)+(1d0-eps)/(3d0*ks))
    kHSl=kf+(1d0-eps)/(1d0/(ks-kf)+eps/(3d0*kf))
  else; kHSu=ks; kHSl=ks; endif
  kgeo=ks**(1d0-eps)*kf**eps
  kEMT=(kpar+kser)/2d0
  do i=1,200
    f1=(1d0-eps)*(ks-kEMT)/(ks+2d0*kEMT)+eps*(kf-kEMT)/(kf+2d0*kEMT)
    df=-(1d0-eps)*3d0*ks/(ks+2d0*kEMT)**2-eps*3d0*kf/(kf+2d0*kEMT)**2
    if(abs(df)<1d-30) exit
    dk=-f1/df; kEMT=kEMT+dk
    if(kEMT<1d-10) kEMT=1d-10
    if(abs(dk)<1d-12) exit
  enddo
  write(*,'(A)') '============================================'
  write(*,'(A)') '  POROUS MEDIA EFFECTIVE CONDUCTIVITY'
  write(*,'(A)') '============================================'
  write(*,'(A)') ''
  write(*,'(A)') '--- INPUTS ---'
  write(*,'(A,F12.4,A)') '  k_solid                 = ',ks,' W/mK'
  write(*,'(A,F12.6,A)') '  k_fluid                 = ',kf,' W/mK'
  write(*,'(A,F10.4)')    '  Porosity epsilon        = ',eps
  write(*,'(A,F10.2,A)') '  T_hot                   = ',Th,' C'
  write(*,'(A,F10.2,A)') '  T_cold                  = ',Tc,' C'
  write(*,'(A,F10.4,A)') '  Thickness L             = ',Lth,' m'
  write(*,'(A,F10.4,A)') '  Area A                  = ',A,' m2'
  write(*,'(A,F10.2)')    '  k_solid/k_fluid ratio   = ',ks/kf
  write(*,'(A)') ''
  write(*,'(A)') '--- MODEL RESULTS ---'
  write(*,'(A)') '  Model              k_eff[W/mK]   Q[W]         R_th[K/W]'
  write(*,'(A)') '  -------------------------------------------------------'
  Q=kpar*A*dT/Lth; Rth=Lth/(kpar*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Parallel (Voigt) ',kpar,Q,Rth
  Q=kser*A*dT/Lth; Rth=Lth/(kser*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Series (Reuss)   ',kser,Q,Rth
  Q=kHSu*A*dT/Lth; Rth=Lth/(kHSu*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Hashin-Shtrik Up ',kHSu,Q,Rth
  Q=kHSl*A*dT/Lth; Rth=Lth/(kHSl*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Hashin-Shtrik Lo ',kHSl,Q,Rth
  Q=kEMT*A*dT/Lth; Rth=Lth/(kEMT*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'EMT (Bruggeman)  ',kEMT,Q,Rth
  Q=kgeo*A*dT/Lth; Rth=Lth/(kgeo*A)
  write(*,'(2X,A,2X,F10.6,2X,F12.4,2X,F12.6)') 'Geometric Mean   ',kgeo,Q,Rth
  write(*,'(A)') ''
  write(*,'(A)') '--- POROSITY SWEEP ---'
  write(*,'(A)') '  eps      k_par     k_ser     k_HSu     k_HSl     k_EMT     k_geo'
  write(*,'(A)') '  ---------------------------------------------------------------------'
  do i=1,25
    es=0.01d0+0.98d0*dble(i-1)/24d0
    kpar=es*kf+(1d0-es)*ks
    kser=1d0/(es/kf+(1d0-es)/ks)
    if(abs(kf-ks)>1d-12) then
      kHSu=ks+es/(1d0/(kf-ks)+(1d0-es)/(3d0*ks))
      kHSl=kf+(1d0-es)/(1d0/(ks-kf)+es/(3d0*kf))
    else; kHSu=ks; kHSl=ks; endif
    kgeo=ks**(1d0-es)*kf**es
    keff=(kpar+kser)/2d0
    do j=1,200
      f1=(1d0-es)*(ks-keff)/(ks+2d0*keff)+es*(kf-keff)/(kf+2d0*keff)
      f2=-(1d0-es)*3d0*ks/(ks+2d0*keff)**2-es*3d0*kf/(kf+2d0*keff)**2
      if(abs(f2)<1d-30) exit; keff=keff-f1/f2
      if(keff<1d-10) keff=1d-10; if(abs(f1)<1d-12) exit
    enddo
    write(*,'(2X,F6.3,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5)') es,es*kf+(1d0-es)*ks,kser,kHSu,kHSl,keff,kgeo
  enddo
  write(*,'(A)') ''
  write(*,'(A)') '--- CORRELATIONS ---'
  write(*,'(A)') '  Parallel (Voigt): k=eps*kf+(1-eps)*ks (upper bound)'
  write(*,'(A)') '  Series (Reuss):   1/k=eps/kf+(1-eps)/ks (lower bound)'
  write(*,'(A)') '  Hashin-Shtrikman: tightest bounds for isotropic composites'
  write(*,'(A)') '  EMT/Bruggeman:    self-consistent effective medium'
  write(*,'(A)') '  Geometric:        k=ks^(1-eps)*kf^eps'
  write(*,'(A)') '  Ref: Kaviany, Heat Transfer in Porous Media, Ch.3'
  write(*,'(A)') '       Incropera Ch.3, Nield & Bejan Ch.2'
end program porous_media_cond