Add the passes of XSTAR's npass

docs/usage.md

@@ -66,4 +66,4 @@ What differs:
 - The atomic database of the XSTAR package (`$HEADAS/refdata/atdb.fits`) is not the one of its source tree: hydrogen has other records, and the temperatures of the tree are 3-6% lower than the package's. Use the one that you want to compare with.
 - With `vturbi > 0` XSTAR 2.59j smooths its continuum so that it is not absorbed below 20 keV (`gsmooth2`): Radix does not, and the comparisons with it use `vturbi = 0`.
 - To reproduce the rounded constants of XSTAR (`0.861707 eV` for k × 10⁴ K, 12.56 for 4π, ...) use `set_constants!(ucalc_constants())` first; the default is CODATA 2022.
-- Not done: more than one pass (`npass`).
+- `passes=3` (XSTAR's `npass`) repeats the march with the optical depths of the lines and edges between each zone and the outer edge, as the escape probabilities of the outer side. It changes the ion fractions by up to 2.5% on the thick slab and follows XSTAR's change to 0.15%. `passes` must be odd.

src/transfer.jl

@@ -297,7 +297,7 @@ its own radius with the radiation of the incident spectrum attenuated by the con
 `exp(-dpthc)` (XSTAR's `trnfrc`) and mapped (`map_spectrum`), and with the escape probabilities of the lines and recombination
 edges that the optical depths of those zones give. With `equilibrium=true` the temperature of each is the thermal equilibrium
 (`thermal_equilibrium`, from `T` and the `xee` of the previous zone), otherwise `T` is kept and the electron fraction iterated
-(`ionization_balance`). `attenuate=false` leaves the spectrum as it is; `lines=false` leaves the lines out of the continuum opacity (`add_line!`), `luminous=false` does not add up the luminosities of the lines and edges (`Luminosities`), with the turbulent speed `vturb` (km/s). `processes` is a collection of `AbstractContinuum` processes (`standard_processes(compton)`) whose heating, cooling and
+(`ionization_balance`). `attenuate=false` leaves the spectrum as it is; `lines=false` leaves the lines out of the continuum opacity (`add_line!`), `luminous=false` does not add up the luminosities of the lines and edges (`Luminosities`), with the turbulent speed `vturb` (km/s). `passes` (odd, XSTAR's `npass`) repeats the march over the same zones: the even passes go from the outer edge inwards at the temperatures of the zones to give the optical depths of the lines and edges beyond each zone, which the odd passes use as the outward depths of the escape probabilities. `processes` is a collection of `AbstractContinuum` processes (`standard_processes(compton)`) whose heating, cooling and
 opacity are those of the zones; `cfrac` is the covering fraction of the escape probabilities (give `Thomson` the same). Further keywords go to those functions.
 
 After each zone the opacities are added to the depths (`stpcut`): the lines and edges to the `OpticalDepths`, the continuum
@@ -309,8 +309,50 @@ fields `r`, `Δr`, `radiation`, `escape`, `opacity`), the final `OpticalDepths`
 not added to the radiation.
 """
 function march_zones(mixture::Mixture, ntot, processes, E, L, zones; T=100.0, equilibrium=false, xee=1.0, cfrac=0.0,
-        attenuate=true, diffuse=true, lines=true, luminous=true, vturb=default_turbulence, lfast=photoionization_lfast, kw...)
+        attenuate=true, diffuse=true, lines=true, luminous=true, vturb=default_turbulence, lfast=photoionization_lfast, passes=1, kw...)
+    isodd(passes) || throw(ArgumentError("passes must be odd: the last pass goes in the same direction as the first"))
+    result = sweep(mixture, ntot, processes, E, L, zones; T, equilibrium, xee, cfrac, attenuate, diffuse, lines, luminous, vturb, lfast, kw...)
+    for pass in 2:passes
+        grid = [(zone.r, zone.Δr) for zone in result.zones]
+        if iseven(pass)
+            outer = backward(mixture, result, processes, E, L; cfrac, lfast, vturb, kw...)
+            result = (; result..., outer)
+        else
+            result = sweep(mixture, ntot, processes, E, L, grid; T, equilibrium, xee, cfrac, attenuate, diffuse, lines, luminous, vturb, lfast,
+                outer=result.outer, kw...)
+        end
+    end
+    (; result..., passes)
+end
+
+# backward: the sweep from the outer zone inwards at the temperatures and densities of the zones of `result` (XSTAR's second pass), with the
+# radiation of the incident spectrum attenuated by the continuum depth before each zone, to give the optical depths of the lines and edges
+# between it and the outer edge: the `outward` depths of each zone (before it is added), as the vectors of `OpticalDepths`
+function backward(mixture::Mixture, result, processes, E, L; cfrac, lfast, vturb, kw...)
+    depths = OpticalDepths(mixture)
+    dpthc = zeros(length(E))
+    before = map(result.zones) do zone
+        previous = copy(dpthc)
+        dpthc .+= zone.opacity.total .* zone.Δr
+        previous
+    end
+    outer = Vector{Vector{Vector{Float64}}}(undef, length(result.zones))
+    for i in reverse(eachindex(result.zones))
+        zone = result.zones[i]
+        outer[i] = map(copy, depths.outward)
+        radiation = map_spectrum(point_source(E, L .* exp.(-before[i]), zone.r))
+        escape = escape_probabilities(mixture, OpticalDepths(result.inner[i], depths.outward); cfrac)
+        n = zone.ntot
+        balance = ionization_balance(mixture, zone.T, n; radiation, escape, xee=zone.xee, lfast, kw...)
+        edges = record_opacities(mixture, balance, zone.T, n; radiation, lfast, vturb)
+        add_zone!(depths, edges, zone.Δr; direction=:outward)
+    end
+    outer
+end
+
+function sweep(mixture::Mixture, ntot, processes, E, L, zones; T, equilibrium, xee, cfrac, attenuate, diffuse, lines, luminous, vturb, lfast, outer=nothing, kw...)
     depths = OpticalDepths(mixture)
+    inner = Vector{Vector{Float64}}[]
     continuum_gas = Continuum(mixture, processes; lfast, vturb)
     dpthc = zeros(length(E))
     dpthcont = zeros(length(E))
@@ -324,7 +366,8 @@ function march_zones(mixture::Mixture, ntot, processes, E, L, zones; T=100.0, eq
         r, Δr = step
         incident = point_source(E, spectrum.L, r)
         radiation = map_spectrum(incident)
-        escape = escape_probabilities(mixture, depths; cfrac)
+        push!(inner, map(copy, depths.inward))
+        escape = escape_probabilities(mixture, outer === nothing ? depths : OpticalDepths(depths.inward, outer[length(results)+1]); cfrac)
         zone = if equilibrium
             thermal_equilibrium(mixture, t -> gas_density(ntot, t, r), processes; radiation, escape, T, xee, lfast, kw...)
         else
@@ -345,7 +388,7 @@ function march_zones(mixture::Mixture, ntot, processes, E, L, zones; T=100.0, eq
         dpthc .+= continuum.total .* Δr
         dpthcont .+= continuum.continuum .* Δr
     end
-    (; zones=results, depths, dpthc, dpthcont, spectrum, luminosities, vturb)
+    (; zones=results, depths, dpthc, dpthcont, spectrum, luminosities, vturb, inner, outer)
 end
 
 """

test/reference/xstar_thick_slab/README.md

@@ -16,3 +16,6 @@ incident one between 0.15 eV and 20 keV (and Thomson scattering outside it), and
 999 points (and than the 0.16% bins of 9999 for T below 4×10⁵ K), the loop ends at the first neighbour with `exp(-earg) = 0`
 and never adds the bin itself, so the smoothed opacity and emissivity arrays are 0 below 20 keV. The continuum is then not absorbed by the
 slab at all. Radix does not do it.
+
+`xout_abund1_npass3.fits` is the same run with `npass=3` (it takes four times as long): the ion fractions of the zones
+after the pass back from the outer edge, which `passes=3` of `slab_model` reproduces.

test/transfer_tests.jl

@@ -402,6 +402,18 @@ function transfer_balance_tests(db)
             @test [z.r - auto.r for z in auto.zones][1:6] ≈ [0; depth[1:5]] atol=3e-4*1.8e16 rtol=1e-3
             @test round.(log10.(1e4 .* cumsum([z.Δr for z in auto.zones])[2:end]); digits=2) == [20.26, 20.56, 20.74, 20.86, 20.95, 21.0]
             @test sum(z.Δr for z in auto.zones) ≈ 1e17
+            # three passes (`npass=3`: the outward optical depths of the lines and edges of the pass back from the outer edge): the ion fractions change by
+            # 2.5% (He II in the first zone) to 0.2% and follow those of XSTAR in the change, which the escape probabilities of the outer side make
+            auto3 = Radix.slab_model(mixture, processes; density=1e4, column=1e21, logξ=1.0, luminosity=1e8*1e38, α=-1.0, emult=1.0, taumax=5.0, steps=3, T=10.0, iterate=false, vturb=0.0, luminous=false, passes=3)
+            pass3 = fits(joinpath(dir, "xout_abund1_npass3.fits"))[2].data
+            @test auto3.passes == 3 && length(auto3.zones) == 7
+            for k in 1:7
+                shift = auto3.zones[k].fractions[2][2]/auto.zones[k].fractions[2][2]
+                @test shift ≈ Float64(pass3.he_ii[k])/Float64(slabab.he_ii[k]) atol=1.5e-3
+                @test auto3.zones[k].fractions[2][2] ≈ Float64(pass3.he_ii[k]) rtol=4e-3
+            end
+            @test auto3.zones[1].fractions[2][2] < 0.98*auto.zones[1].fractions[2][2]
+            @test_throws ArgumentError Radix.slab_model(mixture, processes; density=1e4, column=1e21, logξ=1.0, luminosity=1e8*1e38, passes=2)
             for k in 2:length(zones)
                 @test thickrun.zones[k].fractions[2][2] ≈ slabab.he_ii[k + 1] rtol=1.5e-3
                 @test thickrun.zones[k].fractions[findfirst(==(8), mixture.Z)][8] ≈ slabab.o_viii[k + 1] rtol=3e-4