Skip to main content

pvlib/
batch.rs

1#![allow(clippy::type_complexity)]
2
3use rayon::prelude::*;
4use chrono::{Datelike, TimeZone};
5use crate::{solarposition, atmosphere, clearsky, irradiance, temperature, iam, inverter};
6
7// ---------------------------------------------------------------------------
8// Solar Position Batch
9// ---------------------------------------------------------------------------
10
11/// Batch solar position calculation for multiple timestamps.
12/// Returns (zenith_vec, azimuth_vec, elevation_vec).
13pub fn solar_position_batch(
14    location: &crate::location::Location,
15    times: &[chrono::DateTime<chrono_tz::Tz>],
16) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>), spa::SpaError> {
17    let results: Result<Vec<_>, _> = times.par_iter()
18        .map(|t| solarposition::get_solarposition(location, *t))
19        .collect();
20    let results = results?;
21    let zenith = results.iter().map(|r| r.zenith).collect();
22    let azimuth = results.iter().map(|r| r.azimuth).collect();
23    let elevation = results.iter().map(|r| r.elevation).collect();
24    Ok((zenith, azimuth, elevation))
25}
26
27/// Convenience batch solar position for UTC `NaiveDateTime` timestamps.
28///
29/// Internally creates a `Location` with the UTC timezone, converts each
30/// `NaiveDateTime` to `DateTime<Tz>`, and delegates to [`solar_position_batch`].
31pub fn solar_position_batch_utc(
32    latitude: f64,
33    longitude: f64,
34    altitude: f64,
35    times: &[chrono::NaiveDateTime],
36) -> Result<(Vec<f64>, Vec<f64>, Vec<f64>), spa::SpaError> {
37    let location = crate::location::Location::new(latitude, longitude, chrono_tz::UTC, altitude, "UTC");
38    let datetimes: Vec<chrono::DateTime<chrono_tz::Tz>> = times
39        .iter()
40        .map(|ndt| chrono::Utc.from_utc_datetime(ndt).with_timezone(&chrono_tz::UTC))
41        .collect();
42    solar_position_batch(&location, &datetimes)
43}
44
45// ---------------------------------------------------------------------------
46// Atmosphere Batch
47// ---------------------------------------------------------------------------
48
49/// Batch relative airmass for an array of zenith angles.
50pub fn airmass_relative_batch(zenith: &[f64]) -> Vec<f64> {
51    zenith.par_iter()
52        .map(|z| atmosphere::get_relative_airmass(*z))
53        .collect()
54}
55
56/// Batch absolute airmass.
57pub fn airmass_absolute_batch(airmass_relative: &[f64], pressure: f64) -> Vec<f64> {
58    airmass_relative.par_iter()
59        .map(|am| atmosphere::get_absolute_airmass(*am, pressure))
60        .collect()
61}
62
63// ---------------------------------------------------------------------------
64// Clear Sky Batch
65// ---------------------------------------------------------------------------
66
67/// Batch Ineichen clear sky model.
68/// Returns (ghi_vec, dni_vec, dhi_vec).
69pub fn ineichen_batch(
70    zenith: &[f64],
71    airmass_absolute: &[f64],
72    linke_turbidity: f64,
73    altitude: f64,
74) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
75    assert_eq!(zenith.len(), airmass_absolute.len(), "zenith and airmass_absolute must have the same length");
76    let results: Vec<_> = zenith.par_iter()
77        .zip(airmass_absolute.par_iter())
78        .map(|(z, am)| clearsky::ineichen(*z, *am, linke_turbidity, altitude, 1364.0))
79        .collect();
80    let ghi = results.iter().map(|r| r.ghi).collect();
81    let dni = results.iter().map(|r| r.dni).collect();
82    let dhi = results.iter().map(|r| r.dhi).collect();
83    (ghi, dni, dhi)
84}
85
86/// Batch Bird clear sky model.
87pub fn bird_batch(
88    zenith: &[f64],
89    airmass_relative: &[f64],
90    aod380: f64,
91    aod500: f64,
92    precipitable_water: f64,
93) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
94    assert_eq!(zenith.len(), airmass_relative.len(), "zenith and airmass_relative must have the same length");
95    let results: Vec<_> = zenith.par_iter()
96        .zip(airmass_relative.par_iter())
97        .map(|(z, am)| clearsky::bird_default(*z, *am, aod380, aod500, precipitable_water))
98        .collect();
99    let ghi = results.iter().map(|r| r.ghi).collect();
100    let dni = results.iter().map(|r| r.dni).collect();
101    let dhi = results.iter().map(|r| r.dhi).collect();
102    (ghi, dni, dhi)
103}
104
105// ---------------------------------------------------------------------------
106// Irradiance Batch
107// ---------------------------------------------------------------------------
108
109/// Batch AOI calculation.
110pub fn aoi_batch(
111    surface_tilt: f64,
112    surface_azimuth: f64,
113    solar_zenith: &[f64],
114    solar_azimuth: &[f64],
115) -> Vec<f64> {
116    assert_eq!(solar_zenith.len(), solar_azimuth.len(), "solar_zenith and solar_azimuth must have the same length");
117    solar_zenith.par_iter()
118        .zip(solar_azimuth.par_iter())
119        .map(|(z, a)| irradiance::aoi(surface_tilt, surface_azimuth, *z, *a))
120        .collect()
121}
122
123/// Batch extraterrestrial radiation.
124pub fn extra_radiation_batch(day_of_year: &[i32]) -> Vec<f64> {
125    day_of_year.par_iter()
126        .map(|d| irradiance::get_extra_radiation(*d))
127        .collect()
128}
129
130/// Batch Erbs decomposition. Returns (dni_vec, dhi_vec).
131pub fn erbs_batch(
132    ghi: &[f64],
133    zenith: &[f64],
134    day_of_year: &[u32],
135    dni_extra: &[f64],
136) -> (Vec<f64>, Vec<f64>) {
137    let n = ghi.len();
138    assert_eq!(zenith.len(), n, "zenith len mismatch");
139    assert_eq!(day_of_year.len(), n, "day_of_year len mismatch");
140    assert_eq!(dni_extra.len(), n, "dni_extra len mismatch");
141    let results: Vec<_> = ghi.par_iter()
142        .zip(zenith.par_iter())
143        .zip(day_of_year.par_iter())
144        .zip(dni_extra.par_iter())
145        .map(|(((g, z), d), e)| irradiance::erbs(*g, *z, *d, *e))
146        .collect();
147    let dni = results.iter().map(|r| r.0).collect();
148    let dhi = results.iter().map(|r| r.1).collect();
149    (dni, dhi)
150}
151
152/// Batch DISC decomposition. Returns (dni_vec, kt_vec, airmass_vec).
153pub fn disc_batch(
154    ghi: &[f64],
155    solar_zenith: &[f64],
156    day_of_year: &[i32],
157    pressure: Option<f64>,
158) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
159    let n = ghi.len();
160    assert_eq!(solar_zenith.len(), n, "solar_zenith len mismatch");
161    assert_eq!(day_of_year.len(), n, "day_of_year len mismatch");
162    let results: Vec<_> = ghi.par_iter()
163        .zip(solar_zenith.par_iter())
164        .zip(day_of_year.par_iter())
165        .map(|((g, z), d)| irradiance::disc(*g, *z, *d, pressure))
166        .collect();
167    let dni = results.iter().map(|r| r.dni).collect();
168    let kt = results.iter().map(|r| r.kt).collect();
169    let am = results.iter().map(|r| r.airmass).collect();
170    (dni, kt, am)
171}
172
173/// Batch Perez transposition model.
174#[allow(clippy::too_many_arguments)]
175pub fn perez_batch(
176    surface_tilt: f64,
177    surface_azimuth: f64,
178    dhi: &[f64],
179    dni: &[f64],
180    dni_extra: &[f64],
181    solar_zenith: &[f64],
182    solar_azimuth: &[f64],
183    airmass: &[f64],
184    aoi_vals: &[f64],
185) -> Vec<f64> {
186    let n = dhi.len();
187    assert_eq!(dni.len(), n, "dni len mismatch");
188    assert_eq!(dni_extra.len(), n, "dni_extra len mismatch");
189    assert_eq!(solar_zenith.len(), n, "solar_zenith len mismatch");
190    assert_eq!(solar_azimuth.len(), n, "solar_azimuth len mismatch");
191    assert_eq!(airmass.len(), n, "airmass len mismatch");
192    assert_eq!(aoi_vals.len(), n, "aoi_vals len mismatch");
193    (0..n).into_par_iter()
194        .map(|i| irradiance::perez(
195            surface_tilt, surface_azimuth,
196            dhi[i], dni[i], dni_extra[i],
197            solar_zenith[i], solar_azimuth[i],
198            airmass[i], aoi_vals[i],
199        ))
200        .collect()
201}
202
203/// Batch get_total_irradiance. Returns PoaComponents for each timestep.
204#[allow(clippy::too_many_arguments)]
205pub fn total_irradiance_batch(
206    surface_tilt: f64,
207    surface_azimuth: f64,
208    solar_zenith: &[f64],
209    solar_azimuth: &[f64],
210    dni: &[f64],
211    ghi: &[f64],
212    dhi: &[f64],
213    albedo: f64,
214    model: irradiance::DiffuseModel,
215    dni_extra: &[f64],
216    airmass: &[f64],
217) -> Vec<irradiance::PoaComponents> {
218    let n = solar_zenith.len();
219    assert_eq!(solar_azimuth.len(), n, "solar_azimuth len mismatch");
220    assert_eq!(dni.len(), n, "dni len mismatch");
221    assert_eq!(ghi.len(), n, "ghi len mismatch");
222    assert_eq!(dhi.len(), n, "dhi len mismatch");
223    assert_eq!(dni_extra.len(), n, "dni_extra len mismatch");
224    assert_eq!(airmass.len(), n, "airmass len mismatch");
225    (0..n).into_par_iter()
226        .map(|i| irradiance::get_total_irradiance(
227            surface_tilt, surface_azimuth,
228            solar_zenith[i], solar_azimuth[i],
229            dni[i], ghi[i], dhi[i],
230            albedo, model,
231            Some(dni_extra[i]),
232            Some(airmass[i]),
233        ))
234        .collect()
235}
236
237// ---------------------------------------------------------------------------
238// Temperature Batch
239// ---------------------------------------------------------------------------
240
241/// Batch SAPM cell temperature. Returns (cell_temp_vec, module_temp_vec).
242#[allow(clippy::too_many_arguments)]
243pub fn sapm_cell_temperature_batch(
244    poa_global: &[f64],
245    temp_air: &[f64],
246    wind_speed: &[f64],
247    a: f64, b: f64, delta_t: f64, irrad_ref: f64,
248) -> (Vec<f64>, Vec<f64>) {
249    let n = poa_global.len();
250    assert_eq!(temp_air.len(), n, "temp_air len mismatch");
251    assert_eq!(wind_speed.len(), n, "wind_speed len mismatch");
252    let results: Vec<_> = (0..n).into_par_iter()
253        .map(|i| temperature::sapm_cell_temperature(
254            poa_global[i], temp_air[i], wind_speed[i], a, b, delta_t, irrad_ref,
255        ))
256        .collect();
257    let cell = results.iter().map(|r| r.0).collect();
258    let module = results.iter().map(|r| r.1).collect();
259    (cell, module)
260}
261
262/// Batch Faiman cell temperature.
263pub fn faiman_batch(
264    poa_global: &[f64],
265    temp_air: &[f64],
266    wind_speed: &[f64],
267    u0: f64,
268    u1: f64,
269) -> Vec<f64> {
270    let n = poa_global.len();
271    assert_eq!(temp_air.len(), n, "temp_air len mismatch");
272    assert_eq!(wind_speed.len(), n, "wind_speed len mismatch");
273    (0..n).into_par_iter()
274        .map(|i| temperature::faiman(poa_global[i], temp_air[i], wind_speed[i], u0, u1))
275        .collect()
276}
277
278// ---------------------------------------------------------------------------
279// IAM Batch
280// ---------------------------------------------------------------------------
281
282/// Batch Physical IAM.
283pub fn iam_physical_batch(aoi: &[f64], n: f64, k: f64, l: f64) -> Vec<f64> {
284    aoi.par_iter()
285        .map(|a| iam::physical(*a, n, k, l))
286        .collect()
287}
288
289/// Batch ASHRAE IAM.
290pub fn iam_ashrae_batch(aoi: &[f64], b0: f64) -> Vec<f64> {
291    aoi.par_iter()
292        .map(|a| iam::ashrae(*a, b0))
293        .collect()
294}
295
296// ---------------------------------------------------------------------------
297// Inverter Batch
298// ---------------------------------------------------------------------------
299
300/// Batch PVWatts AC power.
301pub fn pvwatts_ac_batch(
302    pdc: &[f64],
303    pdc0: f64,
304    eta_inv_nom: f64,
305    eta_inv_ref: f64,
306) -> Vec<f64> {
307    pdc.par_iter()
308        .map(|p| inverter::pvwatts_ac(*p, pdc0, eta_inv_nom, eta_inv_ref))
309        .collect()
310}
311
312// ---------------------------------------------------------------------------
313// Full Pipeline Batch
314// ---------------------------------------------------------------------------
315
316/// Input data for batch simulation -- one value per timestep.
317#[derive(Debug, Clone)]
318pub struct WeatherSeries {
319    pub times: Vec<chrono::DateTime<chrono_tz::Tz>>,
320    pub ghi: Vec<f64>,
321    pub dni: Vec<f64>,
322    pub dhi: Vec<f64>,
323    pub temp_air: Vec<f64>,
324    pub wind_speed: Vec<f64>,
325    pub albedo: Option<Vec<f64>>,
326}
327
328impl WeatherSeries {
329    /// Construct a `WeatherSeries` from UTC `NaiveDateTime` timestamps and a
330    /// timezone name (e.g. `"US/Eastern"`, `"UTC"`, `"Europe/Berlin"`).
331    ///
332    /// Each `NaiveDateTime` is interpreted as UTC and converted to the target
333    /// timezone.  Returns `Err` if `tz_name` cannot be parsed.
334    pub fn from_utc(
335        timestamps: &[chrono::NaiveDateTime],
336        tz_name: &str,
337        ghi: Vec<f64>,
338        dni: Vec<f64>,
339        dhi: Vec<f64>,
340        temp_air: Vec<f64>,
341        wind_speed: Vec<f64>,
342    ) -> Result<Self, String> {
343        let tz: chrono_tz::Tz = tz_name
344            .parse()
345            .map_err(|_| format!("Unknown timezone: {}", tz_name))?;
346        let times: Vec<chrono::DateTime<chrono_tz::Tz>> = timestamps
347            .iter()
348            .map(|ndt| chrono::Utc.from_utc_datetime(ndt).with_timezone(&tz))
349            .collect();
350        Ok(Self {
351            times,
352            ghi,
353            dni,
354            dhi,
355            temp_air,
356            wind_speed,
357            albedo: None,
358        })
359    }
360}
361
362/// Output from batch simulation -- one value per timestep.
363#[derive(Debug, Clone)]
364pub struct SimulationSeries {
365    pub solar_zenith: Vec<f64>,
366    pub solar_elevation: Vec<f64>,
367    pub solar_azimuth: Vec<f64>,
368    pub airmass: Vec<f64>,
369    pub aoi: Vec<f64>,
370    pub poa_global: Vec<f64>,
371    pub poa_direct: Vec<f64>,
372    pub poa_diffuse: Vec<f64>,
373    pub cell_temperature: Vec<f64>,
374    pub effective_irradiance: Vec<f64>,
375    pub dc_power: Vec<f64>,
376    pub ac_power: Vec<f64>,
377}
378
379impl SimulationSeries {
380    /// Total energy produced in Wh (assuming 1-hour timesteps).
381    ///
382    /// Negative and non-finite AC-power values (e.g. NaN from upstream weather
383    /// data) are excluded from the sum.
384    #[must_use]
385    pub fn total_energy_wh(&self) -> f64 {
386        // `map + max(0.0)` keeps NaN out of the sum (NaN.max(0.0) is NaN, but
387        // a subsequent filter removes it) while being branchless enough for
388        // autovectorisation.
389        self.ac_power
390            .iter()
391            .filter(|p| p.is_finite() && **p > 0.0)
392            .sum()
393    }
394
395    /// Peak AC power in W. Ignores non-finite values.
396    #[must_use]
397    pub fn peak_power(&self) -> f64 {
398        self.ac_power
399            .iter()
400            .copied()
401            .filter(|p| p.is_finite())
402            .fold(0.0_f64, f64::max)
403    }
404
405    /// Capacity factor (ratio of actual energy to theoretical maximum).
406    #[must_use]
407    pub fn capacity_factor(&self, system_capacity_w: f64) -> f64 {
408        let hours = self.ac_power.len() as f64;
409        if hours == 0.0 || system_capacity_w == 0.0 { return 0.0; }
410        self.total_energy_wh() / (system_capacity_w * hours)
411    }
412
413    /// Number of timesteps in this series where the AC-power value is not
414    /// finite (typically NaN propagated from NaN-valued upstream weather).
415    /// Aggregate accessors above silently drop these; call this method if
416    /// you need to detect them explicitly.
417    #[must_use]
418    pub fn nan_count(&self) -> usize {
419        self.ac_power.iter().filter(|p| !p.is_finite()).count()
420    }
421}
422
423/// Batch ModelChain -- runs the full PV simulation pipeline on a time series
424/// using rayon for parallel processing.
425///
426/// This is the main entry point for production batch simulations.
427/// A typical TMY year (8760 hourly timesteps) completes in milliseconds.
428pub struct BatchModelChain {
429    pub location: crate::location::Location,
430    pub surface_tilt: f64,
431    pub surface_azimuth: f64,
432    pub system_capacity_dc: f64,
433    /// Temperature coefficient, e.g. -0.004
434    pub gamma_pdc: f64,
435    pub inverter_capacity: f64,
436    pub inverter_efficiency: f64,
437    pub albedo: f64,
438    pub transposition_model: irradiance::DiffuseModel,
439    /// When true, automatically decompose GHI into DNI/DHI using the Erbs model
440    /// if DNI and DHI are both near zero but GHI is positive.
441    pub auto_decomposition: bool,
442    /// System losses (wiring, connections, mismatch, soiling): 0.0 = no losses, 0.14 = 14% losses.
443    pub system_losses: f64,
444    /// Bifaciality factor: 0.0 = monofacial, 0.65-0.85 typical for bifacial modules.
445    pub bifaciality_factor: f64,
446    /// Ground albedo used for rear-side irradiance calculation.
447    pub bifacial_ground_albedo: f64,
448}
449
450impl BatchModelChain {
451    /// Create a new BatchModelChain with PVWatts-style defaults.
452    pub fn pvwatts(
453        location: crate::location::Location,
454        surface_tilt: f64,
455        surface_azimuth: f64,
456        system_capacity_dc: f64,
457    ) -> Self {
458        Self {
459            location,
460            surface_tilt,
461            surface_azimuth,
462            system_capacity_dc,
463            gamma_pdc: -0.004,
464            inverter_capacity: system_capacity_dc,
465            inverter_efficiency: 0.96,
466            albedo: 0.2,
467            transposition_model: irradiance::DiffuseModel::Perez,
468            auto_decomposition: false,
469            system_losses: 0.0,
470            bifaciality_factor: 0.0,
471            bifacial_ground_albedo: 0.2,
472        }
473    }
474
475    /// Builder: set temperature coefficient.
476    pub fn with_gamma_pdc(mut self, gamma_pdc: f64) -> Self {
477        self.gamma_pdc = gamma_pdc;
478        self
479    }
480
481    /// Builder: set inverter parameters.
482    pub fn with_inverter(mut self, capacity: f64, efficiency: f64) -> Self {
483        self.inverter_capacity = capacity;
484        self.inverter_efficiency = efficiency;
485        self
486    }
487
488    /// Builder: set ground albedo.
489    pub fn with_albedo(mut self, albedo: f64) -> Self {
490        self.albedo = albedo;
491        self
492    }
493
494    /// Builder: set transposition model.
495    pub fn with_transposition(mut self, model: irradiance::DiffuseModel) -> Self {
496        self.transposition_model = model;
497        self
498    }
499
500    /// Builder: enable/disable automatic GHI to DNI/DHI decomposition via Erbs model.
501    pub fn with_auto_decomposition(mut self, enabled: bool) -> Self {
502        self.auto_decomposition = enabled;
503        self
504    }
505
506    /// Builder: set system losses (0.0 = no losses, 0.14 = 14% losses). Clamped to [0.0, 1.0].
507    pub fn with_system_losses(mut self, losses: f64) -> Self {
508        self.system_losses = losses.clamp(0.0, 1.0);
509        self
510    }
511
512    /// Builder: set bifacial parameters.
513    ///
514    /// `bifaciality_factor` is typically 0.65-0.85 for bifacial modules (0.0 = monofacial).
515    /// `ground_albedo` is the albedo used for rear-side irradiance calculation.
516    pub fn with_bifacial(mut self, bifaciality_factor: f64, ground_albedo: f64) -> Self {
517        self.bifaciality_factor = bifaciality_factor;
518        self.bifacial_ground_albedo = ground_albedo;
519        self
520    }
521
522    /// Run the full simulation on a weather time series.
523    ///
524    /// Uses rayon for automatic parallelization across CPU cores.
525    pub fn run(&self, weather: &WeatherSeries) -> Result<SimulationSeries, spa::SpaError> {
526        let n = weather.times.len();
527        assert_eq!(weather.ghi.len(), n);
528        assert_eq!(weather.dni.len(), n, "dni len mismatch");
529        assert_eq!(weather.dhi.len(), n, "dhi len mismatch");
530        assert_eq!(weather.temp_air.len(), n, "temp_air len mismatch");
531        assert_eq!(weather.wind_speed.len(), n, "wind_speed len mismatch");
532        if let Some(albedos) = &weather.albedo {
533            assert_eq!(albedos.len(), n, "albedo len mismatch");
534        }
535
536        let pressure = atmosphere::alt2pres(self.location.altitude);
537        let albedo_default = self.albedo;
538
539        // Run entire pipeline in parallel for each timestep
540        let results: Result<Vec<_>, _> = (0..n).into_par_iter().map(|i| {
541            // 1. Solar position
542            let solpos = solarposition::get_solarposition(&self.location, weather.times[i])?;
543
544            // 2. Airmass
545            let am_rel = atmosphere::get_relative_airmass(solpos.zenith);
546            let am_abs = if am_rel.is_nan() || am_rel <= 0.0 {
547                0.0
548            } else {
549                atmosphere::get_absolute_airmass(am_rel, pressure)
550            };
551
552            // 3. AOI
553            let aoi_val = irradiance::aoi(
554                self.surface_tilt, self.surface_azimuth,
555                solpos.zenith, solpos.azimuth,
556            );
557
558            // 4. Extraterrestrial irradiance. `ordinal()` is O(1) and
559            // allocation-free; the prior `format("%j").parse()` approach
560            // allocated two `String`s per timestep — a measurable cost on
561            // an 8760-row TMY run.
562            let doy = weather.times[i].ordinal() as i32;
563            let dni_extra = irradiance::get_extra_radiation(doy);
564
565            // 4b. Auto-decompose GHI → DNI/DHI if enabled and needed
566            let (dni_i, dhi_i) = if self.auto_decomposition
567                && (weather.dni[i].abs() < 1.0 || weather.dni[i].is_nan())
568                && (weather.dhi[i].abs() < 1.0 || weather.dhi[i].is_nan())
569                && weather.ghi[i] > 0.0
570            {
571                let dni_extra_val = dni_extra;
572                irradiance::erbs(weather.ghi[i], solpos.zenith, doy as u32, dni_extra_val)
573            } else {
574                (weather.dni[i], weather.dhi[i])
575            };
576
577            // 5. Transposition
578            let poa = irradiance::get_total_irradiance(
579                self.surface_tilt, self.surface_azimuth,
580                solpos.zenith, solpos.azimuth,
581                dni_i, weather.ghi[i], dhi_i,
582                weather.albedo.as_ref().map_or(albedo_default, |a| a[i]),
583                self.transposition_model,
584                Some(dni_extra),
585                if am_rel.is_nan() { None } else { Some(am_rel) },
586            );
587
588            // 6. IAM
589            // Note: Uses Physical IAM model (Fresnel/Snell's law) matching the PVWatts configuration.
590            let iam_val = iam::physical(aoi_val, 1.526, 4.0, 0.002);
591
592            // 7. Effective irradiance (spectral modifier = 1.0 for NoLoss)
593            let spectral_modifier = 1.0;
594            let eff_irrad = ((poa.poa_direct * iam_val + poa.poa_diffuse) * spectral_modifier).max(0.0);
595
596            // 8. Cell temperature
597            // Note: Uses PVWatts temperature model (T_cell = T_air + POA*(NOCT-20)/800).
598            let t_cell = weather.temp_air[i] + poa.poa_global * (45.0 - 20.0) / 800.0;
599
600            // 9. DC power
601            let pdc = self.system_capacity_dc * (eff_irrad / 1000.0)
602                * (1.0 + self.gamma_pdc * (t_cell - 25.0));
603            let pdc = pdc.max(0.0);
604
605            // 9b. System losses
606            let pdc = pdc * (1.0 - self.system_losses);
607
608            // 9c. Bifacial rear-side gain (applied at DC level before inverter)
609            let pdc = if self.bifaciality_factor > 0.0 && poa.poa_global > 10.0 {
610                let rear_gain = (self.bifaciality_factor * self.bifacial_ground_albedo
611                    * weather.ghi[i] / poa.poa_global).min(0.25);
612                pdc * (1.0 + rear_gain)
613            } else {
614                pdc
615            };
616
617            // 10. AC power
618            let pac = inverter::pvwatts_ac(
619                pdc, self.system_capacity_dc,
620                self.inverter_efficiency, 0.9637,
621            );
622
623            Ok((solpos.zenith, solpos.azimuth, am_abs, aoi_val,
624                poa.poa_global, poa.poa_direct, poa.poa_diffuse,
625                t_cell, eff_irrad, pdc, pac))
626        }).collect();
627
628        let results = results?;
629
630        Ok(SimulationSeries {
631            solar_zenith: results.iter().map(|r| r.0).collect(),
632            solar_elevation: results.iter().map(|r| 90.0 - r.0).collect(),
633            solar_azimuth: results.iter().map(|r| r.1).collect(),
634            airmass: results.iter().map(|r| r.2).collect(),
635            aoi: results.iter().map(|r| r.3).collect(),
636            poa_global: results.iter().map(|r| r.4).collect(),
637            poa_direct: results.iter().map(|r| r.5).collect(),
638            poa_diffuse: results.iter().map(|r| r.6).collect(),
639            cell_temperature: results.iter().map(|r| r.7).collect(),
640            effective_irradiance: results.iter().map(|r| r.8).collect(),
641            dc_power: results.iter().map(|r| r.9).collect(),
642            ac_power: results.iter().map(|r| r.10).collect(),
643        })
644    }
645}