1use earths::biosphere::ecosystems::{borealforest, coralreef, tropicalrainforest};
2use earths::biosphere::fauna::{PredatorPrey, africanelephant, bluewhale};
3use earths::biosphere::vegetation::{
4 beerlambertlightextinction, coniferousforest, temperategrassland, tropicalbroadleaf,
5};
6use earths::geodata::bathymetry::{BathymetryData, OceanBasin};
7use earths::geodata::coordinates::{EARTHFLATTENING, EARTHSEMIMAJORM, EARTHSEMIMINORM, LatLon};
8use earths::geodata::elevation::ElevationProvider;
9use earths::geodata::regions::{RegionDatabase, RegionType};
10use earths::lighting::day_night::{DayNightCycle, DaylightState};
11use earths::lighting::seasons::{
12 AXIALTILTDEG as SEASONSTILT, Season, TROPICALYEARDAYS, VERNALEQUINOXJD, seasonat, subsolarpoint,
13};
14use earths::lighting::solar_position::{EARTHAXIALTILTDEG, SolarPosition};
15use earths::physics::collisions::{chicxulubequivalent, tunguskaequivalent};
16use earths::physics::orbit::{
17 ARGUMENTPERIHELIONDEG, ECCENTRICITY, EarthOrbit, INCLINATIONDEG, LONGITUDEASCENDINGNODEDEG,
18 SEMIMAJORAXIS,
19};
20use earths::physics::rotation::{
21 AXIALTILTDEG, AXIALTILTRAD, EarthRotation, PRECESSIONPERIODYEARS, SIDEREALDAYS, SOLARDAYS,
22};
23use earths::physics::tides::{
24 TidalForce, lunartosolartideratio, neaptideamplitude, springtideamplitude,
25};
26use earths::rendering::atmosphere_scattering::AtmosphereEndpoint;
27use earths::rendering::clouds::{CloudLayer, CloudSystemEndpoint, CloudType};
28use earths::rendering::materials::PbrMaterial;
29use earths::rendering::ocean_rendering::OceanEndpoint;
30use earths::rendering::shaders::{ShaderEndpoint, UniformValue};
31use earths::satellites::artificial::{ArtificialSatellite, Constellation, OrbitType};
32use earths::satellites::moon::{
33 EARTHMOONDISTANCE, LUNARMASS, LUNARORBITALPERIOD, LUNARRADIUS, MoonSource, MoonState,
34};
35use earths::temporal::calendar::{DateTime, J2000JD, SECONDSPERDAY, UNIXEPOCHJD};
36use earths::temporal::epoch::{Epoch, J1950EPOCH, J2000EPOCH, MJDOFFSET};
37use earths::temporal::time_scale::TimeScale;
38use earths::terrain::heightmap::Heightmap;
39use earths::terrain::lod::{Face, LodConfig, LodTerrain};
40use earths::terrain::mesh::TerrainMesh;
41use earths::terrain::texturing::{Biome, BiomeClassifier};
42
43fn main() {
44 let orbit = EarthOrbit::new();
45 let perioddays = orbit.orbitalperioddays();
46 let vperihelion = orbit.velocityatdistance(orbit.perihelionm());
47 let vaphelion = orbit.velocityatdistance(orbit.aphelionm());
48 let energy = orbit.specific_orbital_energy();
49 let angmom = orbit.specific_angular_momentum();
50 let vesc = EarthOrbit::escape_velocity_at_surface();
51 let fsun = orbit.gravitational_force_sun();
52 let r = orbit.current_radius();
53 let vmean = orbit.mean_orbital_velocity();
54 println!(
55 "Orbit: P={:.2}d vperi={:.0}m/s vaph={:.0}m/s vmean={:.0}m/s",
56 perioddays, vperihelion, vaphelion, vmean
57 );
58 println!(
59 " E={:.2e}J/kg L={:.2e}m²/s vesc={:.0}m/s Fsun={:.2e}N r={:.3e}m",
60 energy, angmom, vesc, fsun, r
61 );
62 println!(
63 " a={:.3e}m e={} i={}° Ω={}° ω={}°",
64 SEMIMAJORAXIS,
65 ECCENTRICITY,
66 INCLINATIONDEG,
67 LONGITUDEASCENDINGNODEDEG,
68 ARGUMENTPERIHELIONDEG
69 );
70
71 let rot = EarthRotation::new();
72 let vequator = rot.surfacevelocityatlatitude(0.0);
73 let v45 = rot.surfacevelocityatlatitude(45.0);
74 let accel = rot.centripetalaccelerationatlatitude(0.0);
75 let coriolis = rot.coriolisparameter(45.0);
76 let moi = rot.momentofinertia();
77 let ke = rot.rotationalkineticenergy();
78 let lrot = rot.angular_momentum();
79 let prec = rot.precession_rate_rad_per_year();
80 let daylensummer = rot.day_length_variation_due_to_tilt(172, 60.0);
81 let daylenwinter = rot.day_length_variation_due_to_tilt(355, 60.0);
82 println!(
83 "Rotation: veq={:.0}m/s v45={:.0}m/s ac={:.4}m/s² f45={:.5e}",
84 vequator, v45, accel, coriolis
85 );
86 println!(
87 " I={:.3e}kg·m² KE={:.3e}J L={:.3e}kg·m²/s prec={:.2e}rad/yr",
88 moi, ke, lrot, prec
89 );
90 println!(
91 " day@60°N: summer={:.1}h winter={:.1}h",
92 daylensummer, daylenwinter
93 );
94 println!(
95 " sidereal={:.4}s solar={}s tilt={}°={:.4}rad precperiod={}yr",
96 SIDEREALDAYS, SOLARDAYS, AXIALTILTDEG, AXIALTILTRAD, PRECESSIONPERIODYEARS
97 );
98
99 let moontide = TidalForce::frommoon();
100 let suntide = TidalForce::fromsun();
101 let amoon = moontide.tidalacceleration();
102 let asun = suntide.tidalacceleration();
103 let pot = moontide.tidalpotential(0.0);
104 let bulgemoon = moontide.tidalbulgeheight();
105 let bulgesun = suntide.tidalbulgeheight();
106 let grav = moontide.gravitationalattraction();
107 let spring = springtideamplitude();
108 let neap = neaptideamplitude();
109 let ratio = lunartosolartideratio();
110 println!(
111 "Tides: amoon={:.2e} asun={:.2e} ratio={:.2}",
112 amoon, asun, ratio
113 );
114 println!(
115 " bulgemoon={:.3}m bulgesun={:.3}m spring={:.3}m neap={:.3}m",
116 bulgemoon, bulgesun, spring, neap
117 );
118 println!(" potential(0)={:.2e} Fgrav={:.2e}N", pot, grav);
119
120 let chix = chicxulubequivalent();
121 let tung = tunguskaequivalent();
122 let kechix = chix.kineticenergymt();
123 let crater = chix.craterdiameterm(2700.0);
124 let fireball = chix.fireballradiusm();
125 let ejecta = chix.ejectavolumem3(2700.0);
126 let vimp = chix.impactvelocity();
127 let ketung = tung.kineticenergymt();
128 println!(
129 "Chicxulub: {:.2e}Mt crater={:.0}m fireball={:.0}m ejecta={:.2e}m³ vimp={:.0}m/s",
130 kechix, crater, fireball, ejecta, vimp
131 );
132 println!("Tunguska: {:.2e}Mt", ketung);
133
134 let mut moon = MoonState::new();
135 let source = match moon.source.get() {
136 MoonSource::Binary => "binary",
137 MoonSource::Simulation => "simulated",
138 };
139 let pos0 = moon.position();
140 moon.step(3600.0);
141 let pos1 = moon.position();
142 let gmoon = moon.gravityat(EARTHMOONDISTANCE);
143 println!(
144 "Moon({}): pos0=({:.0},{:.0}) pos1h=({:.0},{:.0}) g@dist={:.4e}m/s²",
145 source, pos0.0, pos0.1, pos1.0, pos1.1, gmoon
146 );
147 println!(
148 " M={:.3e}kg R={:.4e}m d={:.3e}m P={:.0}s",
149 LUNARMASS, LUNARRADIUS, EARTHMOONDISTANCE, LUNARORBITALPERIOD
150 );
151
152 let iss = ArtificialSatellite::leo("ISS", 420000.0, 408000.0);
153 let geo = ArtificialSatellite::geo("GEO-Sat", 300.0);
154 let custom = ArtificialSatellite::new("Molniya", 1500.0, 500000.0, 0.7, 1.1);
155 let piss = iss.orbitalperiods();
156 let viss = iss.orbitalvelocityms();
157 let pgeo = geo.orbitalperiods();
158 let gsurf = iss.gravityatsurface();
159 let posiss = iss.position();
160 let posgeo = geo.position();
161
162 let orbittype = |sat: &ArtificialSatellite| -> &str {
163 match sat.orbittype {
164 OrbitType::LEO => "LEO",
165 OrbitType::MEO => "MEO",
166 OrbitType::GEO => "GEO",
167 OrbitType::HEO => "HEO",
168 OrbitType::Custom => "Custom",
169 }
170 };
171
172 println!(
173 "ISS({}): P={:.0}s v={:.0}m/s gsurf={:.2} pos=({:.0},{:.0},{:.0})",
174 orbittype(&iss),
175 piss,
176 viss,
177 gsurf,
178 posiss.0,
179 posiss.1,
180 posiss.2
181 );
182 println!(
183 "GEO({}): P={:.0}s pos=({:.0},{:.0},{:.0})",
184 orbittype(&geo),
185 pgeo,
186 posgeo.0,
187 posgeo.1,
188 posgeo.2
189 );
190 println!(
191 "Custom({}): e={:.1} i={:.1}rad",
192 orbittype(&custom),
193 custom.eccentricity,
194 custom.inclinationrad
195 );
196
197 let mut constellation = Constellation::new("Starlink-test");
198 constellation.add(ArtificialSatellite::leo("S1", 260.0, 550000.0));
199 constellation.add(ArtificialSatellite::leo("S2", 260.0, 550000.0));
200 let mut satstepped = iss;
201 satstepped.step(60.0);
202 constellation.stepall(60.0);
203 let positions = constellation.positions();
204 println!(
205 "Constellation: {} sats positions={}",
206 positions.len(),
207 positions.len()
208 );
209
210 let dt = DateTime::new(2024, 6, 21, 12, 0, 0.0);
211 let jd = dt.tojuliandate();
212 let dtback = DateTime::fromjuliandate(jd);
213 let unix = dt.tounixtimestamp();
214 let dtunix = DateTime::fromunixtimestamp(unix);
215 println!(
216 "DateTime: {}-{}-{} {}:{}:{:.0} JD={:.4} unix={:.0}",
217 dt.year, dt.month, dt.day, dt.hour, dt.minute, dt.second, jd, unix
218 );
219 println!(
220 " roundtrip JD: {}-{}-{} roundtrip unix: {}-{}-{}",
221 dtback.year, dtback.month, dtback.day, dtunix.year, dtunix.month, dtunix.day
222 );
223 println!(
224 " J2000JD={} UNIXJD={} SPD={}",
225 J2000JD, UNIXEPOCHJD, SECONDSPERDAY
226 );
227
228 let mut epoch = Epoch::j2000();
229 let epochmjd = Epoch::frommjd(51544.5);
230 let centuries = epoch.centuriessincej2000();
231 let days = epoch.dayssincej2000();
232 let gmst = epoch.gmstdegrees();
233 let mjd = epoch.tomjd();
234 epoch.advancedays(365.25);
235 let centuries1yr = epoch.centuriessincej2000();
236 let mut epoch2 = Epoch::fromjd(J2000JD);
237 epoch2.advanceseconds(86400.0);
238 println!(
239 "Epoch: T0={:.4} days={:.1} GMST={:.2}° MJD={:.1}",
240 centuries, days, gmst, mjd
241 );
242 println!(
243 " +1yr: T={:.6} frommjd.jd={:.1} epoch2.days={:.1}",
244 centuries1yr,
245 epochmjd.juliandate,
246 epoch2.dayssincej2000()
247 );
248 println!(
249 " J2000={} J1950={} MJDOFF={}",
250 J2000EPOCH, J1950EPOCH, MJDOFFSET
251 );
252
253 let mut ts = TimeScale::realtime();
254 ts.step(1.0);
255 let simdt = ts.simulationdt(1.0);
256 let hours = ts.simhours();
257 let daysts = ts.simdays();
258 let years = ts.simyears();
259 ts.pause();
260 ts.resume();
261 ts.togglepause();
262 ts.setspeed(100.0);
263 let mut ff = TimeScale::fastforward(10.0);
264 ff.step(1.0);
265 let sm = TimeScale::slowmotion(2.0);
266 println!(
267 "TimeScale: dt={:.1} h={:.6} d={:.8} yr={:.10} ffsim={:.1}s smspeed={:.1}",
268 simdt, hours, daysts, years, ff.simulationtimes, sm.speedmultiplier
269 );
270
271 let sunpos = SolarPosition::compute(jd, 48.8, 2.3);
272 let above = sunpos.isabovehorizon();
273 let distsun = sunpos.distancem();
274 println!(
275 "Solar: el={:.1}° az={:.1}° above={} dist={:.3e}m tilt={}°",
276 sunpos.elevationdeg, sunpos.azimuthdeg, above, distsun, EARTHAXIALTILTDEG
277 );
278
279 let dnc = DayNightCycle::new(jd);
280 let stateparis = dnc.stateat(48.8, 2.3);
281 let statetext = match stateparis {
282 DaylightState::Day => "day",
283 DaylightState::Night => "night",
284 DaylightState::CivilTwilight => "civiltwilight",
285 DaylightState::NauticalTwilight => "nauticaltwilight",
286 DaylightState::AstronomicalTwilight => "astrotwilight",
287 };
288 let terminator = dnc.terminatorpoints(36);
289 let ambient = dnc.ambientlight(48.8, 2.3);
290 println!(
291 "DayNight: Paris={} ambient={:.2} terminatorpts={}",
292 statetext,
293 ambient,
294 terminator.len()
295 );
296
297 let seasonstate = seasonat(jd, 48.8);
298 let seasonname = match seasonstate.seasonnorth {
299 Season::Spring => "spring",
300 Season::Summer => "summer",
301 Season::Autumn => "autumn",
302 Season::Winter => "winter",
303 };
304 let southname = match seasonstate.seasonsouth {
305 Season::Spring => "spring",
306 Season::Summer => "summer",
307 Season::Autumn => "autumn",
308 Season::Winter => "winter",
309 };
310 let (sublat, sublon) = subsolarpoint(jd);
311 println!(
312 "Season: N={} S={} decl={:.2}° dayh={:.2} subsolar=({:.1},{:.1})",
313 seasonname,
314 southname,
315 seasonstate.solardeclinationdeg,
316 seasonstate.daylengthhours,
317 sublat,
318 sublon
319 );
320 println!(
321 " tilt={}° tropyear={}d vernaljd={}",
322 SEASONSTILT, TROPICALYEARDAYS, VERNALEQUINOXJD
323 );
324
325 let paris = LatLon::new(48.8566, 2.3522, 35.0);
326 let tokyo = LatLon::new(35.6762, 139.6503, 40.0);
327 let ecefparis = paris.toecef();
328 let cart = paris.tocartesiansimple();
329 let distpt = paris.distanceto(&tokyo);
330 let latlonback = ecefparis.tolatlon();
331 println!(
332 "Paris->ECEF: ({:.0},{:.0},{:.0}) cart=({:.0},{:.0},{:.0})",
333 ecefparis.x, ecefparis.y, ecefparis.z, cart[0], cart[1], cart[2]
334 );
335 println!(
336 " dist(Paris-Tokyo)={:.0}m roundtriplat={:.4}°",
337 distpt, latlonback.latdeg
338 );
339 println!(
340 " f={:.9} a={:.0}m b={:.3}m",
341 EARTHFLATTENING, EARTHSEMIMAJORM, EARTHSEMIMINORM
342 );
343
344 let regions = RegionDatabase::continents();
345 let region = regions.pointinregion(48.8, 2.3);
346 if let Some(r) = region {
347 let rtype = match r.regiontype {
348 RegionType::Continent => "continent",
349 RegionType::Ocean => "ocean",
350 RegionType::Sea => "sea",
351 RegionType::Country => "country",
352 RegionType::Island => "island",
353 };
354 println!(
355 "Region(48.8,2.3): {} ({}) area={:.0}km²",
356 r.name, rtype, r.areakm2
357 );
358 }
359
360 let elevprov = ElevationProvider::global(30.0);
361 let eleveverest = elevprov.sample(27.988, 86.925);
362 let notables = ElevationProvider::notableelevations();
363 println!(
364 "Elevation(Everest)={:.0}m notables={}",
365 eleveverest,
366 notables.len()
367 );
368
369 let bathy = BathymetryData::global(60.0);
370 let depthmariana = bathy.sample(11.35, 142.2);
371 let isocean = bathy.isocean(0.0, -30.0);
372 let stats = BathymetryData::oceanstats();
373 let basins = OceanBasin::majorbasins();
374 println!(
375 "Bathy(Mariana)={:.0}m isocean(0,-30)={} basins={} avgdepth={:.0}m",
376 depthmariana,
377 isocean,
378 basins.len(),
379 stats.avgdepthm
380 );
381
382 let mut hmap = Heightmap::new(64);
383 hmap.generateprocedural(42, 6, 0.5, 2.0);
384 let hsample = hmap.sample(45.0, 10.0);
385 let rat = hmap.radiusat(45.0, 10.0);
386 println!(
387 "Heightmap(64): sample(45,10)={:.1}m radius={:.0}m",
388 hsample, rat
389 );
390
391 let config = LodConfig::default();
392 let mut lod = LodTerrain::new(config);
393 lod.update([6371000.0 + 1000.0, 0.0, 0.0]);
394
395 let mesh = TerrainMesh::fromregion(0.0, 10.0, 0.0, 10.0, 8, &|lat, lon| hmap.sample(lat, lon));
396 println!(
397 "Mesh: {} vertices {} triangles",
398 mesh.vertexcount(),
399 mesh.trianglecount()
400 );
401
402 let classifier = BiomeClassifier::default();
403 let biome = classifier.classify(2000.0, 45.0, 0.6);
404 let biomename = match biome {
405 Biome::Ocean => "ocean",
406 Biome::Beach => "beach",
407 Biome::Desert => "desert",
408 Biome::Grassland => "grassland",
409 Biome::Forest => "forest",
410 Biome::Tundra => "tundra",
411 Biome::Snow => "snow",
412 Biome::Mountain => "mountain",
413 Biome::Volcanic => "volcanic",
414 Biome::Taiga => "taiga",
415 };
416 let splat = classifier.splat(2000.0, 45.0, 0.6);
417 println!(
418 "Biome(2000m,45°,0.6)={} splattop={:.2}",
419 biomename, splat.weights[0].1
420 );
421
422 let face = Face::PosZ;
423 println!("LOD face={:?}", face);
424
425 let all_mats = PbrMaterial::all_earth();
426 let total_subsurface: f32 = all_mats.iter().map(|m| m.subsurface).sum();
427 let avg_ior: f32 = all_mats.iter().map(|m| m.ior).sum::<f32>() / all_mats.len() as f32;
428 let avg_normal: f32 =
429 all_mats.iter().map(|m| m.normal_strength).sum::<f32>() / all_mats.len() as f32;
430 let emissive_count = all_mats
431 .iter()
432 .filter(|m| m.emissive[0] > 0.0 || m.emissive[1] > 0.0 || m.emissive[2] > 0.0)
433 .count();
434 println!(
435 "Materials({}) avg_ior={:.3} avg_normal={:.2} total_sss={:.2} emissive={}",
436 all_mats.len(),
437 avg_ior,
438 avg_normal,
439 total_subsurface,
440 emissive_count
441 );
442 for m in &all_mats {
443 let luminance = m.albedo[0] * 0.3 + m.albedo[1] * 0.59 + m.albedo[2] * 0.11;
444 let fresnel_r0 = ((m.ior - 1.0) / (m.ior + 1.0)).powi(2);
445 println!(
446 " {}: L={:.3} R={:.2} M={:.2} F₀={:.4} IOR={:.3} SSS={:.2} N={:.2} α={:.2} E=({:.2},{:.2},{:.2})",
447 m.name,
448 luminance,
449 m.roughness,
450 m.metallic,
451 fresnel_r0,
452 m.ior,
453 m.subsurface,
454 m.normal_strength,
455 m.albedo[3],
456 m.emissive[0],
457 m.emissive[1],
458 m.emissive[2]
459 );
460 }
461
462 let shterrain = ShaderEndpoint::terrain();
463 let shatm = ShaderEndpoint::atmosphere();
464 let shocean = ShaderEndpoint::ocean();
465 let shnights = ShaderEndpoint::night_lights();
466 let total_uniform_bytes: usize = [&shterrain, &shatm, &shocean, &shnights]
467 .iter()
468 .map(|sh| {
469 sh.name.len()
470 + sh.uniforms
471 .iter()
472 .map(|u| match &u.value {
473 UniformValue::Float(_) => 4,
474 UniformValue::Vec3(_) => 12,
475 UniformValue::Vec4(_) => 16,
476 UniformValue::Int(_) => 4,
477 })
478 .sum::<usize>()
479 })
480 .sum();
481 println!(
482 "Shaders: terrain={} atm={} ocean={} night={} total_bytes={}",
483 shterrain.uniforms.len(),
484 shatm.uniforms.len(),
485 shocean.uniforms.len(),
486 shnights.uniforms.len(),
487 total_uniform_bytes
488 );
489 for u in &shterrain.uniforms {
490 let val = match &u.value {
491 UniformValue::Float(f) => format!("f{:.2}", f),
492 UniformValue::Vec3(v) => format!("v3({:.1},{:.1},{:.1})", v[0], v[1], v[2]),
493 UniformValue::Vec4(v) => format!("v4({:.1},{:.1},{:.1},{:.1})", v[0], v[1], v[2], v[3]),
494 UniformValue::Int(i) => format!("i{}", i),
495 };
496 println!(" {}={}", u.name, val);
497 }
498
499 let atm = AtmosphereEndpoint::earth();
500 let scale_height_ratio = atm.rayleigh_scale_height_m / atm.mie_scale_height_m;
501 let atmo_shell_vol = 4.0 / 3.0
502 * std::f64::consts::PI
503 * ((atm.planet_radius_m + atm.atmosphere_height_m).powi(3) - atm.planet_radius_m.powi(3));
504 let air_density_sl = atm.sea_level_pressure_pa
505 / (8.314_462_618 / atm.mean_molar_mass_kg_mol * atm.sea_level_temperature_k);
506 let total_atmo_mass_kg = air_density_sl * atmo_shell_vol * 0.5;
507 let total_rayleigh_scatter: f64 = atm
508 .species
509 .iter()
510 .map(|s| {
511 let n_s = atm.sea_level_number_density_m3 * s.volume_fraction;
512 let king =
513 (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
514 let dn = s.refractive_index_stp - 1.0;
515 8.0 * std::f64::consts::PI.powi(3) * (2.0 * dn).powi(2) * king
516 / (3.0 * n_s * (550e-9_f64).powi(4))
517 })
518 .sum();
519 let ozone_optical_depth = atm.ozone_column_du
520 * 1e-5
521 * (-((25000.0 - atm.ozone_peak_altitude_m) / 5000.0).powi(2)).exp();
522 let g = atm.mie_asymmetry_g;
523 let hg_forward = (1.0 - g * g) / (1.0 + g * g - 2.0 * g).powf(1.5);
524 println!(
525 "Atmosphere: M_atm={:.3e}kg H_ray/H_mie={:.2} σ_ray550={:.4e} O₃τ={:.4e} HG_fwd={:.3} S₀={:.0}W/m²",
526 total_atmo_mass_kg,
527 scale_height_ratio,
528 total_rayleigh_scatter,
529 ozone_optical_depth,
530 hg_forward,
531 atm.sun_irradiance_w_m2
532 );
533 for s in &atm.species {
534 let king = (6.0 + 3.0 * s.depolarization_factor) / (6.0 - 7.0 * s.depolarization_factor);
535 println!(
536 " {} ({}): f={:.6} M={:.4e}kg/mol n={:.6} δ={:.3} K={:.4}",
537 s.name,
538 s.symbol,
539 s.volume_fraction,
540 s.molar_mass_kg_mol,
541 s.refractive_index_stp,
542 s.depolarization_factor,
543 king
544 );
545 }
546 println!(
547 " β_rgb=({:.4e},{:.4e},{:.4e}) mie_coeff={:.2e}",
548 atm.rayleigh_coefficients_rgb[0],
549 atm.rayleigh_coefficients_rgb[1],
550 atm.rayleigh_coefficients_rgb[2],
551 atm.mie_coefficient
552 );
553
554 let atlantic = OceanEndpoint::earth_atlantic();
555 let pacific = OceanEndpoint::earth_pacific();
556 let arctic = OceanEndpoint::earth_arctic();
557 for ocean in [&atlantic, &pacific, &arctic] {
558 let fresnel_r0 = ((ocean.ior - 1.0) / (ocean.ior + 1.0)).powi(2);
559 let fetch_m = ocean.fetch_km * 1000.0;
560 let grav = ocean.surface_gravity_m_s2;
561 let ws = ocean.wind_speed_m_s.max(0.1);
562 let omega_p = grav / ws * (0.74 * (grav * fetch_m / (ws * ws)).powf(-0.33));
563 let sig_wave_h = 4.0 * (8.1e-3 * grav / omega_p.powi(2)).sqrt();
564 let foam_fraction = if ws > ocean.foam_threshold_wind_m_s {
565 3.84e-6 * (ws - ocean.foam_threshold_wind_m_s).powf(3.41)
566 } else {
567 0.0
568 };
569 let dx = ocean.patch_size_m / ocean.grid_size as f64;
570 let penetration_rgb = [
571 1.0 / ocean.absorption_coefficients_rgb_m[0],
572 1.0 / ocean.absorption_coefficients_rgb_m[1],
573 1.0 / ocean.absorption_coefficients_rgb_m[2],
574 ];
575 let wind_bearing = ocean.wind_direction[1]
576 .atan2(ocean.wind_direction[0])
577 .to_degrees();
578 println!(
579 "Ocean(d={:.0}m,S={:.1}psu,T={:.1}K): R₀={:.4} Hs={:.2}m foam={:.4} dx={:.1}m pen_g={:.1}m wind={:.0}°",
580 ocean.mean_depth_m,
581 ocean.salinity_psu,
582 ocean.surface_temperature_k,
583 fresnel_r0,
584 sig_wave_h,
585 foam_fraction,
586 dx,
587 penetration_rgb[1],
588 wind_bearing
589 );
590 }
591
592 let cumulus = CloudLayer::cumulus();
593 let stratus = CloudLayer::stratus();
594 let cirrus = CloudLayer::cirrus();
595 let cb = CloudLayer::cumulonimbus();
596 let cloudtype = |ct: &CloudType| match ct {
597 CloudType::Cumulus => "Cu",
598 CloudType::Stratus => "St",
599 CloudType::Cirrus => "Ci",
600 CloudType::Cumulonimbus => "Cb",
601 CloudType::Altocumulus => "Ac",
602 CloudType::Stratocumulus => "Sc",
603 CloudType::Nimbostratus => "Ns",
604 };
605 for l in [&cumulus, &stratus, &cirrus, &cb] {
606 let r_eff = l.droplet_radius_um * 1e-6;
607 let lwc = 0.3 * l.density;
608 let tau = 1.5 * 2.0 * lwc * l.thickness_m * l.coverage / (1000.0 * r_eff);
609 let ice_corr = 1.0 - 0.15 * l.ice_fraction;
610 let transmittance = (-(tau * ice_corr)).exp();
611 let wind_dir = l.wind_direction[1].atan2(l.wind_direction[0]).to_degrees();
612 println!(
613 " {} base={:.0}m thick={:.0}m cov={:.2} ρ={:.2} τ={:.2} T={:.4} r={:.0}µm ice={:.1} wind={:.1}m/s@{:.0}°",
614 cloudtype(&l.cloud_type),
615 l.base_altitude_m,
616 l.thickness_m,
617 l.coverage,
618 l.density,
619 tau * ice_corr,
620 transmittance,
621 l.droplet_radius_um,
622 l.ice_fraction,
623 l.wind_speed_m_s,
624 wind_dir
625 );
626 }
627
628 let csys = CloudSystemEndpoint::earth_default();
629 let stormy = CloudSystemEndpoint::earth_stormy();
630 let total_od_default: f64 = csys
631 .layers
632 .iter()
633 .map(|l| {
634 let r = l.droplet_radius_um * 1e-6;
635 1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
636 * (1.0 - 0.15 * l.ice_fraction)
637 })
638 .sum();
639 let total_od_stormy: f64 = stormy
640 .layers
641 .iter()
642 .map(|l| {
643 let r = l.droplet_radius_um * 1e-6;
644 1.5 * 2.0 * 0.3 * l.density * l.thickness_m * l.coverage / (1000.0 * r)
645 * (1.0 - 0.15 * l.ice_fraction)
646 })
647 .sum();
648 println!(
649 "CloudSystem: default {} layers τ_total={:.2} stormy {} layers τ_total={:.2}",
650 csys.layers.len(),
651 total_od_default,
652 stormy.layers.len(),
653 total_od_stormy
654 );
655
656 let rainforest = tropicalrainforest();
657 let boreal = borealforest();
658 let reef = coralreef();
659 println!("Ecosystems:");
660 for eco in [&rainforest, &boreal, &reef] {
661 let sp = eco.totalspecies();
662 let ind = eco.totalindividuals();
663 let shannon = eco.shannonindex();
664 let simpson = eco.simpsondiversity();
665 let areasp = eco.expectedspeciesfromarea(0.25, 100.0);
666 let npp = eco.totalnppgcyr();
667 let turnover = eco.biomassturnovertimeyr(15000.0);
668 println!(
669 " {} species {} ind H'={:.2} D={:.4} SAR={:.0} NPP={:.0} τ={:.1}yr",
670 sp, ind, shannon, simpson, areasp, npp, turnover
671 );
672 }
673
674 let tb = tropicalbroadleaf();
675 let tg = temperategrassland();
676 let cf = coniferousforest();
677 let photo = tb.photosynthesisrate(400.0, 1500.0, 28.0);
678 let canopy = tb.canopyphotosynthesis(photo);
679 let transp = tg.transpirationmmday(1.5, 0.3);
680 let nppveg = cf.nppkgcm2yr(15.0);
681 let crt = cf.carbonresidencetimeyr(nppveg);
682 let light = beerlambertlightextinction(2000.0, 4.0, 0.5);
683 println!(
684 "Vegetation: photo={:.2} canopy={:.2} transp={:.2}mm/d NPP={:.3} CRT={:.1}yr light={:.1}",
685 photo, canopy, transp, nppveg, crt, light
686 );
687
688 let elephant = africanelephant();
689 let whale = bluewhale();
690 let grel = elephant.growthrate();
691 let proj = elephant.projectforward(10.0);
692 let metab = whale.metabolicratew();
693 let range = elephant.homerangekm2();
694 let gentime = whale.generationtimeyears();
695 let lifespan = elephant.maxlifespanyears();
696 println!(
697 "Elephant: r={:.3} N(+10yr)={:.0} range={:.1}km² lifespan={:.0}yr",
698 grel, proj, range, lifespan
699 );
700 println!("Whale: metab={:.0}W gen={:.1}yr", metab, gentime);
701
702 let mut pp = PredatorPrey {
703 prey: africanelephant(),
704 predator: bluewhale(),
705 attackrate: 0.01,
706 conversionefficiency: 0.1,
707 predatordeathrate: 0.05,
708 };
709 let preyr = pp.preygrowthrate();
710 let predr = pp.predatorgrowthrate();
711 pp.step(0.1);
712 println!(
713 "Lotka-Volterra: preyr={:.1} predr={:.1} -> prey={:.1} pred={:.2}",
714 preyr, predr, pp.prey.count, pp.predator.count
715 );
716}