From 81f5976992870683ddc6e7a4c49ac7bb6d181e65 Mon Sep 17 00:00:00 2001
From: Devin Kees <6895754+djkees@users.noreply.github.com>
Date: Tue, 1 Sep 2026 11:25:25 -0600
Subject: [PATCH] fix: seed FAC supersonic-area-ratio Newton iteration from
 throat gamma_s (#28)

Uses soln%throat_idx instead of hardcoded station 2, matching RP-1311
Eq. 6.22 and the pattern already used elsewhere in rocket.f90. For IAC
problems throat_idx=2 so behavior is unchanged; for FAC problems
throat_idx=4, so this now seeds from the throat instead of infinity.

Fixes #10

(cherry picked from commit f0e53a88bfbcd6fd44686e792600d97d42a66c8c)
---
 source/rocket.f90 | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/source/rocket.f90 b/source/rocket.f90
index 9e43710..dcbc9f9 100644
--- a/source/rocket.f90
+++ b/source/rocket.f90
@@ -893,7 +893,7 @@ subroutine RocketSolver_solve_supar(self, soln, idx, pc, supar, h_inf, s0, weigh
             if (supar(i) < 2.0d0) then
                 ln_pinf_pe = ln_pinf_pt + sqrt(3.294d0*(log(supar(i))**2.0d0) + 1.535d0*log(supar(i)))
             else if (supar(i) >= 2.0d0) then
-                ln_pinf_pe = soln%eq_partials(2)%gamma_s + 1.4d0*log(supar(i))
+                ln_pinf_pe = soln%eq_partials(soln%throat_idx)%gamma_s + 1.4d0*log(supar(i))
             end if
 
             do j = 1, max_iter_area
@@ -987,7 +987,7 @@ subroutine RocketSolver_solve_supar_frozen(self, soln, idx, n_frz, pc, supar, h_
             if (supar(i) < 2.0d0) then
                 ln_pinf_pe = ln_pinf_pt + sqrt(3.294d0*(log(supar(i))**2.0d0) + 1.535d0*log(supar(i)))
             else if (supar(i) >= 2.0d0) then
-                ln_pinf_pe = soln%eq_partials(2)%gamma_s + 1.4d0*log(supar(i))
+                ln_pinf_pe = soln%eq_partials(soln%throat_idx)%gamma_s + 1.4d0*log(supar(i))
             end if
 
             do j = 1, max_iter_area