program hx_pinch
    implicit none

    integer, parameter :: max_streams = 100
    integer, parameter :: max_temps = 200

    ! Inputs
    double precision :: dT_min
    integer :: N
    integer :: stream_type(max_streams) ! 1 = Hot, 2 = Cold
    double precision :: Ts(max_streams), Tt(max_streams), CP(max_streams)

    ! Shifted temps
    double precision :: ts_shift(max_streams), tt_shift(max_streams)
    double precision :: shifted_temps(max_temps)
    double precision :: unique_temps(max_temps)
    integer :: m, i, j, k, iostat_val

    ! Intervals
    double precision :: dH(max_temps)
    double precision :: t_high, t_low, sum_cp_hot, sum_cp_cold
    double precision :: s_min, s_max

    ! Cascades
    double precision :: cascade_I(0:max_temps)
    double precision :: cascade_I_corr(0:max_temps)
    double precision :: cascade_I_min, Q_H_min, Q_C_min, Q_recup_max
    double precision :: T_pinch_shifted, T_pinch_hot, T_pinch_cold

    ! Hot Composite
    double precision :: hot_temps(max_temps)
    double precision :: unique_hot_temps(max_temps)
    integer :: num_hot_temps, num_unique_hot
    double precision :: hot_dH(max_temps)
    double precision :: hot_H(max_temps)
    double precision :: Q_hot_total

    ! Cold Composite
    double precision :: cold_temps(max_temps)
    double precision :: unique_cold_temps(max_temps)
    integer :: num_cold_temps, num_unique_cold
    double precision :: cold_dH(max_temps)
    double precision :: cold_H(max_temps)
    double precision :: Q_cold_total

    double precision :: val, temp_val
    logical :: found

    ! Read inputs
    read(*,*,iostat=iostat_val) dT_min
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid dT_min'
        stop
    end if
    read(*,*,iostat=iostat_val) N
    if (iostat_val /= 0 .or. N <= 0 .or. N > max_streams) then
        write(*,*) 'ERROR: Invalid number of streams'
        stop
    end if

    do i = 1, N
        read(*,*,iostat=iostat_val) stream_type(i), Ts(i), Tt(i), CP(i)
        if (iostat_val /= 0) then
            write(*,*) 'ERROR: Invalid stream data at index ', i
            stop
        end if
        if (stream_type(i) /= 1 .and. stream_type(i) /= 2) then
            write(*,*) 'ERROR: Stream type must be 1 (Hot) or 2 (Cold) at index ', i
            stop
        end if
        if (CP(i) < 0.0d0) then
            write(*,*) 'ERROR: CP must be positive or zero at index ', i
            stop
        end if
    end do

    ! 1. Shift Temperatures
    do i = 1, N
        if (stream_type(i) == 1) then
            ts_shift(i) = Ts(i) - dT_min / 2.0d0
            tt_shift(i) = Tt(i) - dT_min / 2.0d0
        else
            ts_shift(i) = Ts(i) + dT_min / 2.0d0
            tt_shift(i) = Tt(i) + dT_min / 2.0d0
        end if
    end do

    ! 2. Gather unique shifted temperatures
    do i = 1, N
        shifted_temps(2*i - 1) = ts_shift(i)
        shifted_temps(2*i)     = tt_shift(i)
    end do

    m = 0
    do i = 1, 2*N
        val = shifted_temps(i)
        found = .false.
        do j = 1, m
            if (abs(unique_temps(j) - val) < 1.0d-5) then
                found = .true.
                exit
            end if
        end do
        if (.not. found) then
            m = m + 1
            unique_temps(m) = val
        end if
    end do

    ! Sort unique shifted temperatures descending
    do i = 1, m-1
        do j = i+1, m
            if (unique_temps(j) > unique_temps(i)) then
                temp_val = unique_temps(i)
                unique_temps(i) = unique_temps(j)
                unique_temps(j) = temp_val
            end if
        end do
    end do

    ! 3. Interval heat balances
    do j = 1, m-1
        t_high = unique_temps(j)
        t_low = unique_temps(j+1)
        sum_cp_hot = 0.0d0
        sum_cp_cold = 0.0d0
        do i = 1, N
            s_min = min(ts_shift(i), tt_shift(i))
            s_max = max(ts_shift(i), tt_shift(i))
            if (s_min <= t_low + 1.0d-5 .and. s_max >= t_high - 1.0d-5) then
                if (stream_type(i) == 1) then
                    sum_cp_hot = sum_cp_hot + CP(i)
                else
                    sum_cp_cold = sum_cp_cold + CP(i)
                end if
            end if
        end do
        dH(j) = (t_high - t_low) * (sum_cp_hot - sum_cp_cold)
    end do

    ! 4. Heat Cascade
    cascade_I(0) = 0.0d0
    do j = 1, m-1
        cascade_I(j) = cascade_I(j-1) + dH(j)
    end do

    cascade_I_min = 0.0d0
    do j = 0, m-1
        if (cascade_I(j) < cascade_I_min) then
            cascade_I_min = cascade_I(j)
        end if
    end do

    Q_H_min = 0.0d0
    if (cascade_I_min < 0.0d0) then
        Q_H_min = -cascade_I_min
    end if

    cascade_I_corr(0) = Q_H_min
    do j = 1, m-1
        cascade_I_corr(j) = cascade_I_corr(j-1) + dH(j)
    end do
    Q_C_min = cascade_I_corr(m-1)

    ! 5. Pinch point
    T_pinch_shifted = -999.0d0
    do j = 0, m-1
        if (abs(cascade_I_corr(j)) < 1.0d-5) then
            T_pinch_shifted = unique_temps(j+1)
            exit
        end if
    end do

    if (T_pinch_shifted == -999.0d0) then
        ! Fallback if no exact zero, use the minimum of cascade
        T_pinch_shifted = unique_temps(m)
    end if

    T_pinch_hot = T_pinch_shifted + dT_min / 2.0d0
    T_pinch_cold = T_pinch_shifted - dT_min / 2.0d0

    ! Total heat duty changes
    Q_hot_total = 0.0d0
    Q_cold_total = 0.0d0
    do i = 1, N
        if (stream_type(i) == 1) then
            Q_hot_total = Q_hot_total + CP(i) * abs(Ts(i) - Tt(i))
        else
            Q_cold_total = Q_cold_total + CP(i) * abs(Ts(i) - Tt(i))
        end if
    end do

    Q_recup_max = Q_hot_total - Q_C_min

    ! 6. Hot Composite Curve construction
    num_hot_temps = 0
    do i = 1, N
        if (stream_type(i) == 1) then
            num_hot_temps = num_hot_temps + 1
            hot_temps(num_hot_temps) = Ts(i)
            num_hot_temps = num_hot_temps + 1
            hot_temps(num_hot_temps) = Tt(i)
        end if
    end do

    num_unique_hot = 0
    do i = 1, num_hot_temps
        val = hot_temps(i)
        found = .false.
        do j = 1, num_unique_hot
            if (abs(unique_hot_temps(j) - val) < 1.0d-5) then
                found = .true.
                exit
            end if
        end do
        if (.not. found) then
            num_unique_hot = num_unique_hot + 1
            unique_hot_temps(num_unique_hot) = val
        end if
    end do

    ! Sort unique hot temps descending
    do i = 1, num_unique_hot-1
        do j = i+1, num_unique_hot
            if (unique_hot_temps(j) > unique_hot_temps(i)) then
                temp_val = unique_hot_temps(i)
                unique_hot_temps(i) = unique_hot_temps(j)
                unique_hot_temps(j) = temp_val
            end if
        end do
    end do

    ! Integrate Hot Composite Curve Enthalpy
    ! We start at lowest temperature (unique_hot_temps(num_unique_hot)) with H = 0.0
    hot_H(num_unique_hot) = 0.0d0
    do j = num_unique_hot-1, 1, -1
        t_low = unique_hot_temps(j+1)
        t_high = unique_hot_temps(j)
        sum_cp_hot = 0.0d0
        do i = 1, N
            if (stream_type(i) == 1) then
                s_min = min(Ts(i), Tt(i))
                s_max = max(Ts(i), Tt(i))
                if (s_min <= t_low + 1.0d-5 .and. s_max >= t_high - 1.0d-5) then
                    sum_cp_hot = sum_cp_hot + CP(i)
                end if
            end if
        end do
        hot_dH(j) = (t_high - t_low) * sum_cp_hot
        hot_H(j) = hot_H(j+1) + hot_dH(j)
    end do


    ! 7. Cold Composite Curve construction
    num_cold_temps = 0
    do i = 1, N
        if (stream_type(i) == 2) then
            num_cold_temps = num_cold_temps + 1
            cold_temps(num_cold_temps) = Ts(i)
            num_cold_temps = num_cold_temps + 1
            cold_temps(num_cold_temps) = Tt(i)
        end if
    end do

    num_unique_cold = 0
    do i = 1, num_cold_temps
        val = cold_temps(i)
        found = .false.
        do j = 1, num_unique_cold
            if (abs(unique_cold_temps(j) - val) < 1.0d-5) then
                found = .true.
                exit
            end if
        end do
        if (.not. found) then
            num_unique_cold = num_unique_cold + 1
            unique_cold_temps(num_unique_cold) = val
        end if
    end do

    ! Sort unique cold temps descending
    do i = 1, num_unique_cold-1
        do j = i+1, num_unique_cold
            if (unique_cold_temps(j) > unique_cold_temps(i)) then
                temp_val = unique_cold_temps(i)
                unique_cold_temps(i) = unique_cold_temps(j)
                unique_cold_temps(j) = temp_val
            end if
        end do
    end do

    ! Integrate Cold Composite Curve Enthalpy
    ! We start at lowest temperature (unique_cold_temps(num_unique_cold)) with H = 0.0,
    ! then in output we shift by Q_C_min.
    cold_H(num_unique_cold) = 0.0d0
    do j = num_unique_cold-1, 1, -1
        t_low = unique_cold_temps(j+1)
        t_high = unique_cold_temps(j)
        sum_cp_cold = 0.0d0
        do i = 1, N
            if (stream_type(i) == 2) then
                s_min = min(Ts(i), Tt(i))
                s_max = max(Ts(i), Tt(i))
                if (s_min <= t_low + 1.0d-5 .and. s_max >= t_high - 1.0d-5) then
                    sum_cp_cold = sum_cp_cold + CP(i)
                end if
            end if
        end do
        cold_dH(j) = (t_high - t_low) * sum_cp_cold
        cold_H(j) = cold_H(j+1) + cold_dH(j)
    end do


    ! Output results to stdout in KEY=value format
    write(*,'(A,F18.4)') 'QH_MIN=', Q_H_min
    write(*,'(A,F18.4)') 'QC_MIN=', Q_C_min
    write(*,'(A,F18.4)') 'Q_RECUP_MAX=', Q_recup_max
    write(*,'(A,F18.4)') 'T_PINCH_HOT=', T_pinch_hot
    write(*,'(A,F18.4)') 'T_PINCH_COLD=', T_pinch_cold

    ! Output Hot Composite coordinates (sorted descending)
    ! They are already stored with unique_hot_temps descending, hot_H matching them.
    write(*,'(A,I4)') 'HOT_COMP_PTS_COUNT=', num_unique_hot
    write(*,'(A)',advance='no') 'HOT_COMP_H='
    do j = 1, num_unique_hot
        if (j > 1) write(*,'(A)',advance='no') ','
        write(*,'(F18.4)',advance='no') hot_H(j)
    end do
    write(*,*) ''

    write(*,'(A)',advance='no') 'HOT_COMP_T='
    do j = 1, num_unique_hot
        if (j > 1) write(*,'(A)',advance='no') ','
        write(*,'(F18.4)',advance='no') unique_hot_temps(j)
    end do
    write(*,*) ''

    ! Output Cold Composite coordinates (sorted descending)
    ! They are already stored with unique_cold_temps descending, cold_H matching them.
    ! We shift the H value by Q_C_min.
    write(*,'(A,I4)') 'COLD_COMP_PTS_COUNT=', num_unique_cold
    write(*,'(A)',advance='no') 'COLD_COMP_H='
    do j = 1, num_unique_cold
        if (j > 1) write(*,'(A)',advance='no') ','
        write(*,'(F18.4)',advance='no') cold_H(j) + Q_C_min
    end do
    write(*,*) ''

    write(*,'(A)',advance='no') 'COLD_COMP_T='
    do j = 1, num_unique_cold
        if (j > 1) write(*,'(A)',advance='no') ','
        write(*,'(F18.4)',advance='no') unique_cold_temps(j)
    end do
    write(*,*) ''

end program hx_pinch
