nx_atmosphere.nx source

1// nx_atmosphere.nx -- atmospheric scattering math primitives. 2// 3// Per honest audit 2026-05-16: visual quality of a world depends on 4// atmospheric depth, fog distribution, and sky color. This primitive 5// ships substrate-native math for: 6// - Per-body atmospheric coefficients (Rayleigh + Mie scaling) 7// - Beer-Lambert extinction along a viewing distance 8// - Sky color as a function of zenith + sun angles (Preetham-style) 9// - Fog density as a function of altitude 10// 11// Used by the rasterizer (queued) and by atmospheric quality graders. 12// 13// genealogy_id: rayleigh_1871 + mie_1908 + preetham_1999_practical_sky + 14// hosek_wilkie_2012_analytic_skylight 15// lineage_id: nx_atmosphere_scattering_q14_v1 16 17// nx_safety_envelope: 18// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 19// sil_target: SIL1 20// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 21// verdict: NOT_YET_EVALUATED 22 23import "nx_syscalls.nx" 24import "nx_tier.nx" 25const NX_MAGIC_6029: i64 = 6029 26const NX_MAGIC_2212: i64 = 2212 27const NX_MAGIC_5800: i64 = 5800 28const NX_MAGIC_6200: i64 = 6200 29 30const NX_ATM_Q: nx_int = 16384 31 32// ===== Per-body atmospheric coefficients =========================== 33// These match the body_id space from nx_solar_system_catalog. 34// Coefficients are normalised so Earth = 1.0Q at sea level. 35// 36// Rayleigh: predominantly blue scattering (1/lambda^4 wavelength dep). 37// Mie: predominantly white/aerosol scattering (wavelength-independent). 38// Density scale: relative to Earth sea level. 39// 40// Earth: R=Q, M=Q, density=Q 41// Venus: R=Q*4, M=Q*2, density=Q*92 (dense CO2) 42// Mars: R=Q/100, M=Q/200, density=Q/100 (thin) 43// Titan: R=Q*2, M=Q*3, density=Q (thick orange haze) 44// Airless: all 0 45const NX_BODY_SUN_LOCAL: nx_int = 0 // local alias for SUN 46const NX_BODY_MERCURY_LOCAL: nx_int = 1 47const NX_BODY_VENUS_LOCAL: nx_int = 2 48const NX_BODY_EARTH_LOCAL: nx_int = 3 49const NX_BODY_MARS_LOCAL: nx_int = 4 50const NX_BODY_TITAN_LOCAL: nx_int = 16 51 52func nx_atmosphere_rayleigh_q14(body_id: nx_int) -> nx_int { 53 let q: nx_int = NX_ATM_Q 54 if body_id == NX_BODY_EARTH_LOCAL { return q } 55 if body_id == NX_BODY_VENUS_LOCAL { return q * 4 } 56 if body_id == NX_BODY_MARS_LOCAL { return q / 100 } 57 if body_id == NX_BODY_TITAN_LOCAL { return q * 2 } 58 return 0 59} 60 61func nx_atmosphere_mie_q14(body_id: nx_int) -> nx_int { 62 let q: nx_int = NX_ATM_Q 63 if body_id == NX_BODY_EARTH_LOCAL { return q } 64 if body_id == NX_BODY_VENUS_LOCAL { return q * 2 } 65 if body_id == NX_BODY_MARS_LOCAL { return q / 200 } 66 if body_id == NX_BODY_TITAN_LOCAL { return q * 3 } 67 return 0 68} 69 70// ===== Beer-Lambert extinction ===================================== 71// I = I_0 * exp(-tau * d) 72// Approximate exp(-x) for x in [0, 4Q] via piecewise-linear table. 73// Returns attenuation in Q14: 1.0 = no attenuation, 0 = fully extincted. 74func _atm_exp_neg_q14(x_q14: nx_int) -> nx_int { 75 let q: nx_int = NX_ATM_Q 76 if x_q14 <= 0 { return q } 77 if x_q14 >= 4 * q { return 0 } 78 // Anchors: exp(-0)=1.0, exp(-1)=0.368, exp(-2)=0.135, exp(-3)=0.050, exp(-4)=0.018 79 // In Q14: 16384 / 6029 / 2212 / 819 / 295. 80 if x_q14 <= q { 81 // Linear from 16384 to 6029 over [0, Q]. 82 return q - ((q - NX_MAGIC_6029) * x_q14) / q 83 } 84 if x_q14 <= 2 * q { 85 let off: nx_int = x_q14 - q 86 return NX_MAGIC_6029 - ((NX_MAGIC_6029 - NX_MAGIC_2212) * off) / q 87 } 88 if x_q14 <= 3 * q { 89 let off: nx_int = x_q14 - 2 * q 90 return NX_MAGIC_2212 - ((NX_MAGIC_2212 - 819) * off) / q 91 } 92 let off: nx_int = x_q14 - 3 * q 93 return 819 - ((819 - 295) * off) / q 94} 95 96// Beer-Lambert extinction at a given distance for a given body. 97// distance_q14_km: viewing distance in Q14 km 98// body_id : the body whose atmosphere we're looking through 99// Returns transmittance in Q14 [0, Q]. 100func nx_atmosphere_extinction_q14( 101 distance_q14_km: nx_int, body_id: nx_int 102) -> nx_int { 103 let q: nx_int = NX_ATM_Q 104 let ray: nx_int = nx_atmosphere_rayleigh_q14(body_id) 105 let mie: nx_int = nx_atmosphere_mie_q14(body_id) 106 let tau: nx_int = ray + mie 107 if tau <= 0 { return q } 108 // Optical depth = tau * distance / scale_height (~ 8 km for Earth). 109 // tau_eff_q14 = tau * distance / (8 * q) 110 let scale_h_q: nx_int = 8 * q 111 let tau_eff: nx_int = (tau * distance_q14_km) / scale_h_q 112 return _atm_exp_neg_q14(tau_eff) 113} 114 115// ===== Fog density at altitude ===================================== 116// Fog density decreases exponentially with altitude (Beer-Lambert). 117// fog(h) = base * exp(-h / scale_h) 118// Returns Q14 fog density. 119func nx_atmosphere_fog_density_q14( 120 altitude_q14_km: nx_int, body_id: nx_int 121) -> nx_int { 122 let q: nx_int = NX_ATM_Q 123 let base_density: nx_int = nx_atmosphere_mie_q14(body_id) 124 if base_density <= 0 { return 0 } 125 let scale_h: nx_int = 2 * q // 2 km scale height for fog (vs 8 km for full atm) 126 let frac: nx_int = (altitude_q14_km * q) / scale_h 127 return (base_density * _atm_exp_neg_q14(frac)) / q 128} 129 130// ===== Sky color (simplified Preetham) ============================ 131// Sky color = base_horizon_color * (1 - extinction) + zenith_color * extinction 132// modulated by sun-angle proximity. 133// 134// Writes (R, G, B) Q14 [0, Q] to out_rgb. 135// zenith_angle_q14: 0 = zenith, Q = horizon (90 deg) 136// sun_zenith_q14: sun's zenith angle (0 = noon, Q = sunset) 137// body_id: atmospheric body 138func nx_atmosphere_sky_color_q14( 139 zenith_angle_q14: nx_int, sun_zenith_q14: nx_int, 140 body_id: nx_int, out_rgb: *i64 141) { 142 let q: nx_int = NX_ATM_Q 143 // Earth default colors: zenith deep blue (R=70, G=130, B=240), horizon 144 // pale (R=200, G=220, B=250). 145 var zR: nx_int = 70 * q / 255 146 var zG: nx_int = 130 * q / 255 147 var zB: nx_int = 240 * q / 255 148 var hR: nx_int = 200 * q / 255 149 var hG: nx_int = 220 * q / 255 150 var hB: nx_int = 250 * q / 255 151 152 // Mars: pinkish/orange sky. 153 if body_id == NX_BODY_MARS_LOCAL { 154 zR = 200 * q / 255; zG = 130 * q / 255; zB = 100 * q / 255 155 hR = 240 * q / 255; hG = 180 * q / 255; hB = 130 * q / 255 156 } 157 // Titan: orange haze. 158 if body_id == NX_BODY_TITAN_LOCAL { 159 zR = 220 * q / 255; zG = 120 * q / 255; zB = 30 * q / 255 160 hR = 240 * q / 255; hG = 160 * q / 255; hB = 60 * q / 255 161 } 162 // Venus: yellow-white. 163 if body_id == NX_BODY_VENUS_LOCAL { 164 zR = 220 * q / 255; zG = 200 * q / 255; zB = 100 * q / 255 165 hR = 250 * q / 255; hG = 230 * q / 255; hB = 150 * q / 255 166 } 167 // Mercury / airless: black. 168 if body_id == NX_BODY_MERCURY_LOCAL { 169 zR = 0; zG = 0; zB = 0 170 hR = 0; hG = 0; hB = 0 171 } 172 173 // Mix zenith vs horizon by zenith_angle. 174 let t: nx_int = zenith_angle_q14 175 let one_minus_t: nx_int = q - t 176 var R: nx_int = (zR * one_minus_t + hR * t) / q 177 var G: nx_int = (zG * one_minus_t + hG * t) / q 178 var B: nx_int = (zB * one_minus_t + hB * t) / q 179 180 // Sunset reddening: when sun is near horizon (sun_zenith_q14 ~ Q), 181 // warm the color toward red/orange. Apply if sun_zenith >= 0.7Q. 182 let sunset_band: nx_int = q * 7 / 10 183 if sun_zenith_q14 >= sunset_band { 184 let intensity: nx_int = ((sun_zenith_q14 - sunset_band) * q) / (q - sunset_band) 185 // R bumps up; B drops. 186 R = R + (intensity * (q - R)) / (q * 2) 187 B = B - (intensity * B) / (q * 2) 188 if B < 0 { B = 0 } 189 } 190 191 out_rgb[0] = R 192 out_rgb[1] = G 193 out_rgb[2] = B 194} 195 196// ===== Self-test ==================================================== 197func main() -> i64 { 198 let q: nx_int = NX_ATM_Q 199 200 // T1: Per-body coefficients. 201 if nx_atmosphere_rayleigh_q14(NX_BODY_EARTH_LOCAL) != q { return __syscall(93, 1, 0, 0, 0, 0, 0) } 202 if nx_atmosphere_rayleigh_q14(NX_BODY_MERCURY_LOCAL) != 0 { return __syscall(93, 2, 0, 0, 0, 0, 0) } 203 if nx_atmosphere_mie_q14(NX_BODY_VENUS_LOCAL) != 2 * q { return __syscall(93, 3, 0, 0, 0, 0, 0) } 204 205 // T2: exp(-0) = Q; exp(-1) ~ 6029; exp(-4) ~ 295. 206 if _atm_exp_neg_q14(0) != q { return __syscall(93, 10, 0, 0, 0, 0, 0) } 207 let e1: nx_int = _atm_exp_neg_q14(q) 208 if e1 < NX_MAGIC_5800 { return __syscall(93, 11, 0, 0, 0, 0, 0) } 209 if e1 > NX_MAGIC_6200 { return __syscall(93, 12, 0, 0, 0, 0, 0) } 210 if _atm_exp_neg_q14(5 * q) != 0 { return __syscall(93, 13, 0, 0, 0, 0, 0) } 211 212 // T3: Extinction. Earth at 0 distance -> full transmittance. 213 if nx_atmosphere_extinction_q14(0, NX_BODY_EARTH_LOCAL) != q { 214 return __syscall(93, 20, 0, 0, 0, 0, 0) 215 } 216 // Earth at very long distance -> fully extincted. 217 let ext_far: nx_int = nx_atmosphere_extinction_q14(1000 * q, NX_BODY_EARTH_LOCAL) 218 if ext_far >= q / 4 { return __syscall(93, 21, 0, 0, 0, 0, 0) } 219 // Mercury (no atm) at any distance -> full transmittance. 220 if nx_atmosphere_extinction_q14(1000 * q, NX_BODY_MERCURY_LOCAL) != q { 221 return __syscall(93, 22, 0, 0, 0, 0, 0) 222 } 223 224 // T4: Fog density: high at sea level, low at altitude. 225 let fog_low: nx_int = nx_atmosphere_fog_density_q14(0, NX_BODY_EARTH_LOCAL) 226 let fog_high: nx_int = nx_atmosphere_fog_density_q14(10 * q, NX_BODY_EARTH_LOCAL) 227 if fog_low <= fog_high { return __syscall(93, 30, 0, 0, 0, 0, 0) } 228 229 // T5: Sky color. 230 let rgb: *i64 = (sys_mmap(3 * NX_SIZEOF_NX_INT)) as *i64 231 // Earth, zenith=0 (looking straight up), noon -> deep blue 232 nx_atmosphere_sky_color_q14(0, 0, NX_BODY_EARTH_LOCAL, rgb) 233 // B should be the highest of the three (Earth daytime zenith). 234 if rgb[2] < rgb[0] { return __syscall(93, 40, 0, 0, 0, 0, 0) } 235 if rgb[2] < rgb[1] { return __syscall(93, 41, 0, 0, 0, 0, 0) } 236 237 // Mars zenith -> R highest (orange sky). 238 nx_atmosphere_sky_color_q14(0, 0, NX_BODY_MARS_LOCAL, rgb) 239 if rgb[0] < rgb[2] { return __syscall(93, 50, 0, 0, 0, 0, 0) } 240 241 // Mercury airless -> black. 242 nx_atmosphere_sky_color_q14(0, 0, NX_BODY_MERCURY_LOCAL, rgb) 243 if rgb[0] != 0 { return __syscall(93, 60, 0, 0, 0, 0, 0) } 244 if rgb[1] != 0 { return __syscall(93, 61, 0, 0, 0, 0, 0) } 245 if rgb[2] != 0 { return __syscall(93, 62, 0, 0, 0, 0, 0) } 246 247 // Sunset reddening: Earth, zenith=horizon (Q), sun_zenith near horizon (0.9Q). 248 nx_atmosphere_sky_color_q14(q, q * 9 / 10, NX_BODY_EARTH_LOCAL, rgb) 249 // R should be higher than the noon-Q-zenith earth sky. 250 let rgb_noon: *i64 = (sys_mmap(3 * NX_SIZEOF_NX_INT)) as *i64 251 nx_atmosphere_sky_color_q14(q, 0, NX_BODY_EARTH_LOCAL, rgb_noon) 252 if rgb[0] <= rgb_noon[0] { return __syscall(93, 70, 0, 0, 0, 0, 0) } 253 254 return 0 255}