program gas_mixture
    implicit none
    integer::nsp,iostat_val,i,j,ns
    double precision::y(10),MW(10),Cp(10),mu(10),kk(10)
    double precision::y_sum,Mmix,Cp_mol,Cp_mass,mu_mix,k_mix
    double precision::phi_ij,num,den,R_mix,gamma_mix,rho_STP
    double precision::R_gas,y1_sw,mu_sw,k_sw,y_tmp(10)
    R_gas=8.314462d0
    read(*,*,iostat=iostat_val) nsp
    if(iostat_val/=0.or.nsp<2)then;write(*,*)'ERROR: Bad input.';stop;end if
    if(nsp>10)nsp=10
    do i=1,nsp
      read(*,*,iostat=iostat_val) y(i),MW(i),Cp(i),mu(i),kk(i)
      if(iostat_val/=0)then;write(*,*)'ERROR: Bad input line',i;stop;end if
    end do
    y_sum=0.0d0;do i=1,nsp;y_sum=y_sum+y(i);end do
    if(y_sum>0)then;do i=1,nsp;y(i)=y(i)/y_sum;end do;end if
    Mmix=0;Cp_mol=0
    do i=1,nsp;Mmix=Mmix+y(i)*MW(i);Cp_mol=Cp_mol+y(i)*Cp(i);end do
    Cp_mass=Cp_mol/Mmix*1000.0d0
    R_mix=R_gas/Mmix*1000.0d0
    gamma_mix=Cp_mass/(Cp_mass-R_mix)
    rho_STP=101325.0d0*Mmix/(R_gas*273.15d0*1000.0d0)
    mu_mix=0;k_mix=0
    do i=1,nsp
      num=y(i)*mu(i); den=0
      do j=1,nsp
        phi_ij=(1.0d0+sqrt(mu(i)/mu(j))*(MW(j)/MW(i))**0.25d0)**2 &
               /sqrt(8.0d0*(1.0d0+MW(i)/MW(j)))
        den=den+y(j)*phi_ij
      end do
      mu_mix=mu_mix+num/max(den,1.0d-30)
    end do
    do i=1,nsp
      num=y(i)*kk(i); den=0
      do j=1,nsp
        phi_ij=(1.0d0+sqrt(mu(i)/mu(j))*(MW(j)/MW(i))**0.25d0)**2 &
               /sqrt(8.0d0*(1.0d0+MW(i)/MW(j)))
        den=den+y(j)*phi_ij
      end do
      k_mix=k_mix+num/max(den,1.0d-30)
    end do
    write(*,'(A)')'============================================================'
    write(*,'(A)')'   GAS MIXTURE PROPERTIES (Wilke / Wassilijewa)'
    write(*,'(A)')'============================================================'
    write(*,*)
    write(*,'(A)')'--- INPUTS --------------------------------------------------'
    write(*,'(A,I4)')      '  Number of Species         = ',nsp
    write(*,'(A)')'  #   y_i       MW        Cp[J/molK]  mu[uPa.s]  k[W/mK]'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,nsp
      write(*,'(I4,2X,F8.4,2X,F8.2,4X,F8.2,4X,F8.2,4X,F8.4)') i,y(i),MW(i),Cp(i),mu(i),kk(i)
    end do
    write(*,*)
    write(*,'(A)')'--- MIXTURE PROPERTIES --------------------------------------'
    write(*,'(A,F12.4,A)') '  Molar Mass M_mix         = ',Mmix,' g/mol'
    write(*,'(A,F12.4,A)') '  Cp_mix (molar)            = ',Cp_mol,' J/(mol.K)'
    write(*,'(A,F12.2,A)') '  Cp_mix (mass)             = ',Cp_mass,' J/(kg.K)'
    write(*,'(A,F12.4,A)') '  Viscosity mu_mix          = ',mu_mix,' uPa.s'
    write(*,'(A,F12.6,A)') '  Conductivity k_mix        = ',k_mix,' W/(m.K)'
    write(*,'(A,F12.4,A)') '  Gas Constant R_mix        = ',R_mix,' J/(kg.K)'
    write(*,'(A,F12.4)')   '  Gamma (Cp/Cv) approx      = ',gamma_mix
    write(*,'(A,F12.4,A)') '  Density at STP            = ',rho_STP,' kg/m3'
    write(*,*)
    ns=20
    write(*,'(A)')'--- SWEEP: y(1) VARIATION -----------------------------------'
    write(*,'(A)')'  y(1)        mu_mix      k_mix'
    write(*,'(A)')'  -----------------------------------------------------------'
    do i=1,ns
      y1_sw=dble(i-1)/dble(ns-1)
      y_tmp=0;y_tmp(1)=y1_sw
      if(nsp>1)then
        den=0;do j=2,nsp;den=den+y(j);end do
        if(den>0)then;do j=2,nsp;y_tmp(j)=(1.0d0-y1_sw)*y(j)/den;end do
        else;y_tmp(2)=1.0d0-y1_sw;end if
      end if
      mu_sw=0;k_sw=0
      do j=1,nsp
        num=y_tmp(j)*mu(j);den=0
        do ns=1,nsp
          phi_ij=(1.0d0+sqrt(mu(j)/mu(ns))*(MW(ns)/MW(j))**0.25d0)**2 &
                 /sqrt(8.0d0*(1.0d0+MW(j)/MW(ns)))
          den=den+y_tmp(ns)*phi_ij
        end do
        mu_sw=mu_sw+num/max(den,1.0d-30)
      end do
      do j=1,nsp
        num=y_tmp(j)*kk(j);den=0
        do ns=1,nsp
          phi_ij=(1.0d0+sqrt(mu(j)/mu(ns))*(MW(ns)/MW(j))**0.25d0)**2 &
                 /sqrt(8.0d0*(1.0d0+MW(j)/MW(ns)))
          den=den+y_tmp(ns)*phi_ij
        end do
        k_sw=k_sw+num/max(den,1.0d-30)
      end do
      write(*,'(F8.4,4X,F10.4,4X,F10.6)') y1_sw,mu_sw,k_sw
    end do
    write(*,*)
    write(*,'(A)')'--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)')'  Wilke (1950): mu_mix = sum(yi*mui / sum(yj*phi_ij))'
    write(*,'(A)')'  Wassilijewa (1904): k_mix same form with conductivity.'
    write(*,'(A)')'  phi_ij = [1+sqrt(mui/muj)*(Mj/Mi)^0.25]^2 / sqrt(8(1+Mi/Mj))'
end program gas_mixture