program adsorption_isotherms
    implicit none
    integer :: model, i, n_points, iostat_val
    double precision :: C, Prel, T, qmax, KL, KF, nF, qmBET, CBET, areaMolecule, MWads, massAds
    double precision :: q, theta, surfaceArea, totalAds, x, qL, qF, qB, betlin
    double precision, parameter :: NAvo = 6.02214076d23
    character(len=80) :: modelName

    read(*,*,iostat=iostat_val) model
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid model input.'
        stop
    end if
    read(*,*,iostat=iostat_val) C
    read(*,*,iostat=iostat_val) Prel
    read(*,*,iostat=iostat_val) T
    read(*,*,iostat=iostat_val) qmax
    read(*,*,iostat=iostat_val) KL
    read(*,*,iostat=iostat_val) KF
    read(*,*,iostat=iostat_val) nF
    read(*,*,iostat=iostat_val) qmBET
    read(*,*,iostat=iostat_val) CBET
    read(*,*,iostat=iostat_val) areaMolecule
    read(*,*,iostat=iostat_val) MWads
    read(*,*,iostat=iostat_val) massAds
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all adsorption inputs.'
        stop
    end if
    if (C < 0.0d0 .or. Prel < 0.0d0 .or. Prel >= 1.0d0) then
        write(*,*) 'ERROR: Concentration must be non-negative and relative pressure must be 0 <= P/P0 < 1.'
        stop
    end if
    if (qmax < 0.0d0 .or. KL < 0.0d0 .or. KF < 0.0d0 .or. nF <= 0.0d0 .or. qmBET < 0.0d0 .or. CBET < 0.0d0) then
        write(*,*) 'ERROR: Isotherm parameters must be positive or non-negative as applicable.'
        stop
    end if
    if (areaMolecule < 0.0d0 .or. MWads <= 0.0d0 .or. massAds <= 0.0d0) then
        write(*,*) 'ERROR: Surface-area inputs and mass basis must be valid.'
        stop
    end if

    select case(model)
    case(1)
        modelName = 'Langmuir Isotherm'
        q = langmuir(C, qmax, KL)
    case(2)
        modelName = 'Freundlich Isotherm'
        q = freundlich(C, KF, nF)
    case(3)
        modelName = 'BET Isotherm'
        q = bet(Prel, qmBET, CBET)
    case default
        write(*,*) 'ERROR: Invalid model. Use 1 Langmuir, 2 Freundlich, 3 BET.'
        stop
    end select

    if (qmax > 0.0d0) then
        theta = min(1.0d0, langmuir(C, qmax, KL)/qmax)
    else
        theta = 0.0d0
    end if
    ! qmBET is mmol/g. area molecule is nm2/molecule. Result: m2/g.
    surfaceArea = qmBET*1.0d-3 * NAvo * areaMolecule*1.0d-18
    totalAds = q*massAds

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   ADSORPTION ISOTHERMS ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A,A)')        '  Model Name              = ', trim(modelName)
    write(*,'(A,ES12.4)')   '  Concentration (C)       = ', C
    write(*,'(A,ES12.4)')   '  Relative Pressure       = ', Prel
    write(*,'(A,ES12.4,A)') '  Temperature             = ', T, ' K'
    write(*,*)
    write(*,'(A)') '--- ISOTHERM PARAMETERS -------------------------------------'
    write(*,'(A,ES12.4)')   '  Saturation Capacity (qmax) = ', qmax
    write(*,'(A,ES12.4)')   '  Langmuir Constant (K_L)  = ', KL
    write(*,'(A,ES12.4)')   '  Freundlich K_F          = ', KF
    write(*,'(A,ES12.4)')   '  Freundlich n            = ', nF
    write(*,'(A,ES12.4)')   '  BET Monolayer q_m       = ', qmBET
    write(*,'(A,ES12.4)')   '  BET Constant (C_BET)    = ', CBET
    write(*,'(A,ES12.4,A)') '  Molecular Area          = ', areaMolecule, ' nm2/molecule'
    write(*,*)
    write(*,'(A)') '--- ADSORPTION RESULTS --------------------------------------'
    write(*,'(A,ES12.4)')   '  Adsorbed Capacity (q)   = ', q
    write(*,'(A,ES12.4)')   '  Surface Coverage        = ', theta
    write(*,'(A,ES12.4,A)') '  Surface Area            = ', surfaceArea, ' m2/g'
    write(*,'(A,ES12.4)')   '  Total Adsorbed Amount   = ', totalAds
    write(*,*)

    write(*,'(A)') '--- ISOTHERM PROFILE ----------------------------------------'
    write(*,'(A)') '  X             q_Langmuir    q_Freundlich  q_BET        BET_linear'
    write(*,'(A)') '  ----------------------------------------------------------------'
    n_points = 70
    do i=1,n_points
        x = 0.001d0 + 0.95d0*dble(i-1)/dble(n_points-1)
        qL = langmuir(x*max(C,1.0d0), qmax, KL)
        qF = freundlich(x*max(C,1.0d0), KF, nF)
        qB = bet(x, qmBET, CBET)
        if (qB > 0.0d0 .and. x < 1.0d0) then
            betlin = x/(qB*(1.0d0-x))
        else
            betlin = 0.0d0
        end if
        write(*,'(F10.6,2X,ES12.4,2X,ES12.4,2X,ES12.4,2X,ES12.4)') x, qL, qF, qB, betlin
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Langmuir: q = qmax*K_L*C/(1+K_L*C).'
    write(*,'(A)') '  Freundlich: q = K_F*C^(1/n).'
    write(*,'(A)') '  BET: q = qm*C*x/((1-x)*(1+(C-1)*x)).'
    write(*,'(A)') '  Surface area: S = qm(mol/g)*N_A*area_molecule.'

contains
    double precision function langmuir(Ce, qm, K)
        double precision, intent(in) :: Ce, qm, K
        langmuir = qm*K*Ce/(1.0d0 + K*Ce)
    end function langmuir
    double precision function freundlich(Ce, Kf, nf)
        double precision, intent(in) :: Ce, Kf, nf
        if (Ce <= 0.0d0) then
            freundlich = 0.0d0
        else
            freundlich = Kf*Ce**(1.0d0/nf)
        end if
    end function freundlich
    double precision function bet(x, qm, Cb)
        double precision, intent(in) :: x, qm, Cb
        if (x <= 0.0d0 .or. x >= 1.0d0) then
            bet = 0.0d0
        else
            bet = qm*Cb*x/((1.0d0-x)*(1.0d0 + (Cb-1.0d0)*x))
        end if
    end function bet
end program adsorption_isotherms