1#![allow(clippy::type_complexity)]
2
3use rayon::prelude::*;
4use chrono::{Datelike, TimeZone};
5use crate::{solarposition, atmosphere, clearsky, irradiance, temperature, iam, inverter};
6
7pub 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
27pub 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
45pub fn airmass_relative_batch(zenith: &[f64]) -> Vec<f64> {
51 zenith.par_iter()
52 .map(|z| atmosphere::get_relative_airmass(*z))
53 .collect()
54}
55
56pub 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
63pub 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
86pub 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
105pub 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
123pub 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
130pub 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
152pub 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#[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#[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#[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
262pub 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
278pub 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
289pub 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
296pub 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#[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 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#[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 #[must_use]
385 pub fn total_energy_wh(&self) -> f64 {
386 self.ac_power
390 .iter()
391 .filter(|p| p.is_finite() && **p > 0.0)
392 .sum()
393 }
394
395 #[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 #[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 #[must_use]
418 pub fn nan_count(&self) -> usize {
419 self.ac_power.iter().filter(|p| !p.is_finite()).count()
420 }
421}
422
423pub 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 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 pub auto_decomposition: bool,
442 pub system_losses: f64,
444 pub bifaciality_factor: f64,
446 pub bifacial_ground_albedo: f64,
448}
449
450impl BatchModelChain {
451 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 pub fn with_gamma_pdc(mut self, gamma_pdc: f64) -> Self {
477 self.gamma_pdc = gamma_pdc;
478 self
479 }
480
481 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 pub fn with_albedo(mut self, albedo: f64) -> Self {
490 self.albedo = albedo;
491 self
492 }
493
494 pub fn with_transposition(mut self, model: irradiance::DiffuseModel) -> Self {
496 self.transposition_model = model;
497 self
498 }
499
500 pub fn with_auto_decomposition(mut self, enabled: bool) -> Self {
502 self.auto_decomposition = enabled;
503 self
504 }
505
506 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 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 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 let results: Result<Vec<_>, _> = (0..n).into_par_iter().map(|i| {
541 let solpos = solarposition::get_solarposition(&self.location, weather.times[i])?;
543
544 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 let aoi_val = irradiance::aoi(
554 self.surface_tilt, self.surface_azimuth,
555 solpos.zenith, solpos.azimuth,
556 );
557
558 let doy = weather.times[i].ordinal() as i32;
563 let dni_extra = irradiance::get_extra_radiation(doy);
564
565 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 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 let iam_val = iam::physical(aoi_val, 1.526, 4.0, 0.002);
591
592 let spectral_modifier = 1.0;
594 let eff_irrad = ((poa.poa_direct * iam_val + poa.poa_diffuse) * spectral_modifier).max(0.0);
595
596 let t_cell = weather.temp_air[i] + poa.poa_global * (45.0 - 20.0) / 800.0;
599
600 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 let pdc = pdc * (1.0 - self.system_losses);
607
608 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 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}