program economic_pipe
    implicit none
    integer::iostat_val,i,ns,best_i
    double precision::Q,rho,mu,cp_fac,ec,hrs,eta,f_in,L,rate,life
    double precision::pi_v,D,V,Re,f_sw,dP,Pw,E_cost,C_cost,T_cost
    double precision::CRF,D_opt,V_opt,Re_opt,T_min,E_opt,C_opt
    double precision::D_mm,f_auto
    pi_v=3.14159265358979d0
    read(*,*,iostat=iostat_val) Q;read(*,*,iostat=iostat_val) rho
    read(*,*,iostat=iostat_val) mu;read(*,*,iostat=iostat_val) cp_fac
    read(*,*,iostat=iostat_val) ec;read(*,*,iostat=iostat_val) hrs
    read(*,*,iostat=iostat_val) eta;read(*,*,iostat=iostat_val) f_in
    read(*,*,iostat=iostat_val) L;read(*,*,iostat=iostat_val) rate
    read(*,*,iostat=iostat_val) life
    if(iostat_val/=0)then;write(*,*)'ERROR: Bad input.';stop;end if
    if(Q<=0)Q=0.01d0;if(rho<=0)rho=998.0d0;if(mu<=0)mu=0.001d0
    if(cp_fac<=0)cp_fac=5.0d0;if(ec<=0)ec=0.10d0;if(hrs<=0)hrs=8760.0d0
    if(eta<=0.or.eta>1)eta=0.70d0;if(L<=0)L=100.0d0
    if(rate<=0)rate=5.0d0;if(life<=0)life=20.0d0
    CRF=(rate/100.0d0)*(1.0d0+rate/100.0d0)**life/((1.0d0+rate/100.0d0)**life-1.0d0)
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   ECONOMIC PIPE DIAMETER OPTIMIZER'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,ES10.3,A)') '  Volume Flow Q             = ',Q,' m3/s'
    write(*,'(A,F10.2,A)') '  Fluid Density rho         = ',rho,' kg/m3'
    write(*,'(A,ES10.3,A)') '  Viscosity mu              = ',mu,' Pa.s'
    write(*,'(A,F10.2,A)') '  Pipe Cost Factor          = ',cp_fac,' $/m per mm'
    write(*,'(A,F10.4,A)') '  Energy Cost               = ',ec,' $/kWh'
    write(*,'(A,F10.0,A)') '  Operating Hours           = ',hrs,' h/yr'
    write(*,'(A,F10.2)')   '  Pump Efficiency           = ',eta
    write(*,'(A,F10.2,A)') '  Pipe Length                = ',L,' m'
    write(*,'(A,F10.2,A)') '  Interest Rate             = ',rate,' percent'
    write(*,'(A,F10.0,A)') '  Lifetime                  = ',life,' years'
    write(*,'(A,F10.6)')   '  Capital Recovery Factor   = ',CRF
    write(*,*)
    T_min=1.0d30;best_i=0;D_opt=0;V_opt=0;Re_opt=0;E_opt=0;C_opt=0
    ns=50
    write(*,'(A)')'--- COST VS DIAMETER ----------------------------------------'
    write(*,'(A)')'  D[mm]   V[m/s]  Re          f       E_cost[$/yr] C_cost[$/yr] Total[$/yr]'
    write(*,'(A)')'  --------------------------------------------------------------------------'
    do i=1,ns
      D_mm=15.0d0+dble(i-1)*(500.0d0-15.0d0)/dble(ns-1)
      D=D_mm/1000.0d0
      V=4.0d0*Q/(pi_v*D**2)
      Re=rho*V*D/mu
      if(f_in>0)then;f_sw=f_in
      else
        if(Re>0)then;f_auto=0.25d0/(log10(5.74d0/Re**0.9d0))**2
          if(f_auto<0.005d0)f_auto=0.005d0;if(f_auto>0.1d0)f_auto=0.1d0
          f_sw=f_auto
        else;f_sw=0.02d0;end if
      end if
      dP=f_sw*L/D*rho*V**2/2.0d0
      Pw=Q*dP/eta
      E_cost=Pw*hrs*ec/1000.0d0
      C_cost=cp_fac*D_mm*L*CRF
      T_cost=E_cost+C_cost
      if(T_cost<T_min)then;T_min=T_cost;best_i=i;D_opt=D_mm;V_opt=V;Re_opt=Re
        E_opt=E_cost;C_opt=C_cost;end if
      write(*,'(F7.1,2X,F6.2,2X,ES10.3,2X,F7.5,2X,F10.2,4X,F10.2,4X,F10.2)') &
        D_mm,V,Re,f_sw,E_cost,C_cost,T_cost
    end do
    write(*,*)
    write(*,'(A)')'--- OPTIMAL RESULT ------------------------------------------'
    write(*,'(A,F10.1,A)') '  Optimal Diameter D_opt    = ',D_opt,' mm'
    write(*,'(A,F10.3,A)') '  Optimal Velocity V_opt    = ',V_opt,' m/s'
    write(*,'(A,ES12.3)')  '  Reynolds at D_opt         = ',Re_opt
    write(*,'(A,F10.2,A)') '  Annual Energy Cost        = ',E_opt,' $/yr'
    write(*,'(A,F10.2,A)') '  Annual Capital Cost       = ',C_opt,' $/yr'
    write(*,'(A,F10.2,A)') '  Min Total Annual Cost     = ',T_min,' $/yr'
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  Swamee-Jain: f = 0.25/[log10(5.74/Re^0.9)]^2'
    write(*,'(A)')'  E_cost = (Q*dP/eta)*hours*price/1000'
    write(*,'(A)')'  C_cost = cost_per_mm*D*L*CRF'
    write(*,'(A)')'  CRF = i(1+i)^n / ((1+i)^n - 1)'
end program economic_pipe