program stirling_ericsson
    implicit none
    integer :: cycle_type, i, iostat_val, n_pts
    double precision :: TH, TL, Pmax, Pmin, Vmax, Vmin, R_gas, cv, cp, gamma
    double precision :: Mmol, n_mol
    double precision :: Wnet, Qin, Qout, Qregen, eta_th, eta_carnot
    double precision :: r_comp, r_press, COP_heat, COP_ref
    double precision :: W12, W23, W34, W41, Q12, Q23, Q34, Q41
    double precision :: P1, P2, P3, P4, V1, V2, V3, V4, T1, T2, T3, T4
    double precision :: theta, P_i, V_i, T_i
    double precision, parameter :: PI = 3.141592653589793d0
    double precision, parameter :: Ru = 8.314462d0
    character(len=40) :: cycle_name

    read(*,*,iostat=iostat_val) cycle_type
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Invalid cycle type input.'
        stop
    end if
    read(*,*,iostat=iostat_val) TH
    read(*,*,iostat=iostat_val) TL
    read(*,*,iostat=iostat_val) Vmin
    read(*,*,iostat=iostat_val) Vmax
    read(*,*,iostat=iostat_val) Mmol
    read(*,*,iostat=iostat_val) gamma
    read(*,*,iostat=iostat_val) n_mol
    if (iostat_val /= 0) then
        write(*,*) 'ERROR: Failed to read all cycle inputs.'
        stop
    end if
    if (TH <= TL .or. TH <= 0.0d0 .or. TL <= 0.0d0) then
        write(*,*) 'ERROR: TH must be greater than TL, both positive.'
        stop
    end if
    if (Vmin <= 0.0d0 .or. Vmax <= Vmin) then
        write(*,*) 'ERROR: Volumes must be positive, Vmax > Vmin.'
        stop
    end if
    if (Mmol <= 0.0d0) Mmol = 4.003d0   ! helium default
    if (gamma <= 1.0d0) gamma = 1.667d0
    if (n_mol <= 0.0d0) n_mol = 1.0d0

    R_gas = Ru   ! per mole
    cv = R_gas / (gamma - 1.0d0)
    cp = gamma * cv
    r_comp = Vmax / Vmin

    eta_carnot = 1.0d0 - TL / TH

    select case(cycle_type)
    case(1)
        ! ====== STIRLING CYCLE ======
        ! 1-2: Isothermal compression at TL (V1=Vmax → V2=Vmin)
        ! 2-3: Isochoric heating at Vmin  (TL → TH)
        ! 3-4: Isothermal expansion at TH  (V3=Vmin → V4=Vmax)
        ! 4-1: Isochoric cooling at Vmax   (TH → TL)
        cycle_name = 'Stirling Cycle'

        T1 = TL;  V1 = Vmax; P1 = n_mol*R_gas*T1/V1
        T2 = TL;  V2 = Vmin; P2 = n_mol*R_gas*T2/V2
        T3 = TH;  V3 = Vmin; P3 = n_mol*R_gas*T3/V3
        T4 = TH;  V4 = Vmax; P4 = n_mol*R_gas*T4/V4

        ! Work (per mole * n_mol)
        W12 = n_mol * R_gas * TL * log(V2/V1)   ! negative (compression)
        W23 = 0.0d0                               ! isochoric
        W34 = n_mol * R_gas * TH * log(V4/V3)   ! positive (expansion)
        W41 = 0.0d0                               ! isochoric

        ! Heat
        Q12 = W12                                  ! isothermal: Q = W
        Q23 = n_mol * cv * (TH - TL)             ! isochoric heating
        Q34 = W34                                  ! isothermal: Q = W
        Q41 = n_mol * cv * (TL - TH)             ! isochoric cooling

        ! With perfect regenerator: Q23 recovered from Q41
        Qregen = n_mol * cv * (TH - TL)

    case(2)
        ! ====== ERICSSON CYCLE ======
        ! 1-2: Isothermal compression at TL (P1=Pmin → P2=Pmax)
        ! 2-3: Isobaric heating at Pmax     (TL → TH)
        ! 3-4: Isothermal expansion at TH    (P3=Pmax → P4=Pmin)
        ! 4-1: Isobaric cooling at Pmin      (TH → TL)
        cycle_name = 'Ericsson Cycle'

        Pmin = n_mol * R_gas * TL / Vmax
        Pmax = n_mol * R_gas * TL / Vmin
        r_press = Pmax / Pmin

        T1 = TL; P1 = Pmin; V1 = n_mol*R_gas*T1/P1
        T2 = TL; P2 = Pmax; V2 = n_mol*R_gas*T2/P2
        T3 = TH; P3 = Pmax; V3 = n_mol*R_gas*T3/P3
        T4 = TH; P4 = Pmin; V4 = n_mol*R_gas*T4/P4

        ! Work
        W12 = n_mol * R_gas * TL * log(V2/V1)     ! negative
        W23 = P3 * (V3 - V2)                       ! isobaric expansion work (by gas during heating)
        W34 = n_mol * R_gas * TH * log(V4/V3)     ! positive
        W41 = P1 * (V1 - V4)                       ! isobaric compression work (negative)

        ! Heat
        Q12 = W12                                    ! isothermal
        Q23 = n_mol * cp * (TH - TL)               ! isobaric heating
        Q34 = W34                                    ! isothermal
        Q41 = n_mol * cp * (TL - TH)               ! isobaric cooling

        Qregen = n_mol * cp * (TH - TL)

    case default
        write(*,*) 'ERROR: Cycle type must be 1 (Stirling) or 2 (Ericsson).'
        stop
    end select

    Wnet = W12 + W23 + W34 + W41
    Qin = Q34                         ! heat added during hot isothermal
    if (cycle_type == 2) then
        ! Without regenerator, Qin includes isobaric heating too
        ! With perfect regenerator, only isothermal heat input counts
    end if
    Qout = abs(Q12)

    ! With perfect regenerator: eta = eta_Carnot
    ! Without regenerator:
    if (cycle_type == 1) then
        ! Without regen: Qin_total = Q34 + Q23
        eta_th = Wnet / (Q34 + Q23)
    else
        eta_th = Wnet / (Q34 + Q23)
    end if

    ! COP as heat pump
    COP_heat = (Q34 + Q23) / max(abs(Wnet), 1.0d-30)
    ! COP as refrigerator
    COP_ref = Qout / max(abs(Wnet), 1.0d-30)

    Pmax = max(P1, max(P2, max(P3, P4)))
    Pmin = min(P1, min(P2, min(P3, P4)))

    write(*,'(A)') '============================================================'
    write(*,'(A)') '   STIRLING & ERICSSON CYCLES ENGINE'
    write(*,'(A)') '============================================================'
    write(*,*)
    write(*,'(A)') '--- INPUTS --------------------------------------------------'
    write(*,'(A,A)')        '  Cycle Type                = ', trim(cycle_name)
    write(*,'(A,ES12.4,A)') '  Hot Temperature TH        = ', TH, ' K'
    write(*,'(A,ES12.4,A)') '  Cold Temperature TL       = ', TL, ' K'
    write(*,'(A,ES12.4,A)') '  Min Volume                = ', Vmin, ' m3'
    write(*,'(A,ES12.4,A)') '  Max Volume                = ', Vmax, ' m3'
    write(*,'(A,ES12.4)')   '  Compression Ratio         = ', r_comp
    write(*,'(A,ES12.4)')   '  Gamma                     = ', gamma
    write(*,'(A,ES12.4,A)') '  Molar Mass                = ', Mmol, ' g/mol'
    write(*,'(A,ES12.4)')   '  Moles of Gas              = ', n_mol
    write(*,*)
    write(*,'(A)') '--- STATE POINTS --------------------------------------------'
    write(*,'(A)') '  Point   T[K]          P[Pa]         V[m3]'
    write(*,'(A)') '  -------------------------------------------'
    write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') '  1     ', T1, P1, V1
    write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') '  2     ', T2, P2, V2
    write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') '  3     ', T3, P3, V3
    write(*,'(A,F12.2,2X,ES12.4,2X,ES12.4)') '  4     ', T4, P4, V4
    write(*,*)
    write(*,'(A)') '--- PROCESS ENERGY ------------------------------------------'
    write(*,'(A,ES12.4,A)') '  W12 (1→2)                 = ', W12, ' J'
    write(*,'(A,ES12.4,A)') '  Q12 (1→2)                 = ', Q12, ' J'
    write(*,'(A,ES12.4,A)') '  W23 (2→3)                 = ', W23, ' J'
    write(*,'(A,ES12.4,A)') '  Q23 (2→3)                 = ', Q23, ' J'
    write(*,'(A,ES12.4,A)') '  W34 (3→4)                 = ', W34, ' J'
    write(*,'(A,ES12.4,A)') '  Q34 (3→4)                 = ', Q34, ' J'
    write(*,'(A,ES12.4,A)') '  W41 (4→1)                 = ', W41, ' J'
    write(*,'(A,ES12.4,A)') '  Q41 (4→1)                 = ', Q41, ' J'
    write(*,*)
    write(*,'(A)') '--- CYCLE PERFORMANCE ---------------------------------------'
    write(*,'(A,ES12.4,A)') '  Net Work Wnet             = ', Wnet, ' J'
    write(*,'(A,ES12.4,A)') '  Heat Input (no regen)     = ', Q34+Q23, ' J'
    write(*,'(A,ES12.4,A)') '  Heat Input (with regen)   = ', Q34, ' J'
    write(*,'(A,ES12.4,A)') '  Heat Rejected             = ', abs(Q12)+abs(Q41), ' J'
    write(*,'(A,ES12.4,A)') '  Regenerator Heat          = ', Qregen, ' J'
    write(*,'(A,F12.4)')    '  Thermal Efficiency (no regen)  = ', eta_th
    write(*,'(A,F12.4)')    '  Thermal Efficiency (with regen)= ', eta_carnot
    write(*,'(A,F12.4)')    '  Carnot Efficiency         = ', eta_carnot
    write(*,'(A,F12.4)')    '  COP (heat pump)           = ', COP_heat
    write(*,'(A,F12.4)')    '  COP (refrigerator)        = ', COP_ref
    write(*,'(A,ES12.4,A)') '  Max Pressure              = ', Pmax, ' Pa'
    write(*,'(A,ES12.4,A)') '  Min Pressure              = ', Pmin, ' Pa'
    write(*,*)

    ! P-V diagram data
    n_pts = 50
    write(*,'(A)') '--- PV DIAGRAM DATA -----------------------------------------'
    write(*,'(A)') '  V[m3]         P[Pa]         process'
    write(*,'(A)') '  -----------------------------------------------------------'
    if (cycle_type == 1) then
        ! 1→2: isothermal TL, V1→V2
        do i = 0, n_pts
            V_i = V1 + (V2-V1)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*TL/V_i
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '1-2 iso-T(TL)'
        end do
        ! 2→3: isochoric, V=Vmin
        do i = 0, n_pts
            T_i = TL + (TH-TL)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*T_i/V2
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V2, P_i, '2-3 iso-V'
        end do
        ! 3→4: isothermal TH, V3→V4
        do i = 0, n_pts
            V_i = V3 + (V4-V3)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*TH/V_i
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '3-4 iso-T(TH)'
        end do
        ! 4→1: isochoric, V=Vmax
        do i = 0, n_pts
            T_i = TH + (TL-TH)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*T_i/V4
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V4, P_i, '4-1 iso-V'
        end do
    else
        ! Ericsson
        ! 1→2: isothermal TL
        do i = 0, n_pts
            V_i = V1 + (V2-V1)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*TL/V_i
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '1-2 iso-T(TL)'
        end do
        ! 2→3: isobaric Pmax
        do i = 0, n_pts
            T_i = TL + (TH-TL)*dble(i)/dble(n_pts)
            V_i = n_mol*R_gas*T_i/P3
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P3, '2-3 iso-P'
        end do
        ! 3→4: isothermal TH
        do i = 0, n_pts
            V_i = V3 + (V4-V3)*dble(i)/dble(n_pts)
            P_i = n_mol*R_gas*TH/V_i
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P_i, '3-4 iso-T(TH)'
        end do
        ! 4→1: isobaric Pmin
        do i = 0, n_pts
            T_i = TH + (TL-TH)*dble(i)/dble(n_pts)
            V_i = n_mol*R_gas*T_i/P1
            write(*,'(ES12.4,2X,ES12.4,2X,A)') V_i, P1, '4-1 iso-P'
        end do
    end if
    write(*,*)

    ! T-s diagram data
    write(*,'(A)') '--- TS DIAGRAM DATA -----------------------------------------'
    write(*,'(A)') '  s_rel[J/K]    T[K]          process'
    write(*,'(A)') '  -----------------------------------------------------------'
    call write_ts_data(cycle_type, T1, T2, T3, T4, V1, V2, V3, V4, &
                       P1, P2, P3, P4, n_mol, R_gas, cv, cp, n_pts)
    write(*,*)

    ! Efficiency vs temperature ratio sweep
    write(*,'(A)') '--- EFFICIENCY VS TL/TH SWEEP -------------------------------'
    write(*,'(A)') '  TL/TH         eta_Carnot    eta_no_regen'
    write(*,'(A)') '  -------------------------------------------'
    do i = 1, 40
        theta = 0.1d0 + 0.85d0*dble(i-1)/39.0d0
        write(*,'(F10.4,2X,F10.5,2X,F10.5)') theta, 1.0d0-theta, &
            (1.0d0-theta)*log(r_comp) / &
            (log(r_comp) + cv/R_gas*(1.0d0/theta - 1.0d0)/(1.0d0))
    end do
    write(*,*)
    write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
    write(*,'(A)') '  Stirling: iso-T compression/expansion + iso-V regeneration.'
    write(*,'(A)') '  Ericsson: iso-T compression/expansion + iso-P regeneration.'
    write(*,'(A)') '  With perfect regen: eta = eta_Carnot = 1 - TL/TH.'

contains

    subroutine write_ts_data(ctype, t1i, t2i, t3i, t4i, v1i, v2i, v3i, v4i, &
                             p1i, p2i, p3i, p4i, nm, Rg, cvi, cpi, npts)
        implicit none
        integer, intent(in) :: ctype, npts
        double precision, intent(in) :: t1i,t2i,t3i,t4i,v1i,v2i,v3i,v4i
        double precision, intent(in) :: p1i,p2i,p3i,p4i,nm,Rg,cvi,cpi
        double precision :: s_acc, ds_i, T_loc, frac
        integer :: j
        s_acc = 0.0d0
        if (ctype == 1) then
            ! 1→2: isothermal at TL
            do j = 0, npts
                frac = dble(j)/dble(npts)
                ds_i = nm*Rg*log(v1i + (v2i-v1i)*frac) - nm*Rg*log(v1i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t1i, '1-2'
            end do
            s_acc = s_acc + nm*Rg*log(v2i/v1i)
            ! 2→3: isochoric
            do j = 0, npts
                frac = dble(j)/dble(npts)
                T_loc = t2i + (t3i-t2i)*frac
                ds_i = nm*cvi*log(T_loc/t2i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '2-3'
            end do
            s_acc = s_acc + nm*cvi*log(t3i/t2i)
            ! 3→4: isothermal at TH
            do j = 0, npts
                frac = dble(j)/dble(npts)
                ds_i = nm*Rg*log(v3i + (v4i-v3i)*frac) - nm*Rg*log(v3i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t3i, '3-4'
            end do
            s_acc = s_acc + nm*Rg*log(v4i/v3i)
            ! 4→1: isochoric
            do j = 0, npts
                frac = dble(j)/dble(npts)
                T_loc = t4i + (t1i-t4i)*frac
                ds_i = nm*cvi*log(T_loc/t4i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '4-1'
            end do
        else
            ! Ericsson
            ! 1→2: isothermal
            do j = 0, npts
                frac = dble(j)/dble(npts)
                ds_i = nm*Rg*log(v1i + (v2i-v1i)*frac) - nm*Rg*log(v1i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t1i, '1-2'
            end do
            s_acc = s_acc + nm*Rg*log(v2i/v1i)
            ! 2→3: isobaric
            do j = 0, npts
                frac = dble(j)/dble(npts)
                T_loc = t2i + (t3i-t2i)*frac
                ds_i = nm*cpi*log(T_loc/t2i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '2-3'
            end do
            s_acc = s_acc + nm*cpi*log(t3i/t2i)
            ! 3→4: isothermal
            do j = 0, npts
                frac = dble(j)/dble(npts)
                ds_i = nm*Rg*log(v3i + (v4i-v3i)*frac) - nm*Rg*log(v3i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, t3i, '3-4'
            end do
            s_acc = s_acc + nm*Rg*log(v4i/v3i)
            ! 4→1: isobaric
            do j = 0, npts
                frac = dble(j)/dble(npts)
                T_loc = t4i + (t1i-t4i)*frac
                ds_i = nm*cpi*log(T_loc/t4i)
                write(*,'(ES12.4,2X,F12.2,2X,A)') s_acc+ds_i, T_loc, '4-1'
            end do
        end if
    end subroutine write_ts_data

end program stirling_ericsson