diff --git a/examples/ram_air_kite.jl b/examples/ram_air_kite.jl index 1460e45b..a8d14055 100644 --- a/examples/ram_air_kite.jl +++ b/examples/ram_air_kite.jl @@ -20,10 +20,6 @@ VSM = Solver( aerodynamic_model_type=:VSM, is_with_artificial_damping=false ) -VSM_with_stall_correction = Solver( - aerodynamic_model_type=:VSM, - is_with_artificial_damping=true -) # Setting velocity conditions v_a = 15.0 diff --git a/scripts/polars.jl b/scripts/polars.jl index bf7c5d1e..7f9de6ca 100644 --- a/scripts/polars.jl +++ b/scripts/polars.jl @@ -1,4 +1,4 @@ -using Distributed, Timers, Serialization, SharedArrays +using Distributed, Timers, Serialization, SharedArrays, StaticArrays using Interpolations using Xfoil using ControlPlots @@ -204,6 +204,10 @@ try cl_matrix[:, j], cd_matrix[:, j], cm_matrix[:, j] = run_solve_alpha(alphas, d_trailing_edge_angles[j], reynolds_number, x, y, lower, upper, kite_speed, SPEED_OF_SOUND, x_turn) end + cl_matrix = Matrix{Float64}(cl_matrix) + cd_matrix = Matrix{Float64}(cd_matrix) + cm_matrix = Matrix{Float64}(cm_matrix) + display(cl_matrix) function plot_values(alphas, d_trailing_edge_angles, matrix, interp, name) @@ -213,7 +217,8 @@ try X_data = collect(d_trailing_edge_angles) .+ zeros(length(alphas))' Y_data = collect(alphas)' .+ zeros(length(d_trailing_edge_angles)) - interp_matrix = similar(matrix) + matrix = Matrix{Float64}(matrix) + interp_matrix = zeros(size(matrix)...) int_alphas, int_d_trailing_edge_angles = alphas .+ deg2rad(0.5), d_trailing_edge_angles .+ deg2rad(0.5) interp_matrix .= [interp(alpha, d_trailing_edge_angle) for alpha in int_alphas, d_trailing_edge_angle in int_d_trailing_edge_angles] X_int = collect(int_d_trailing_edge_angles) .+ zeros(length(int_alphas))' @@ -241,7 +246,7 @@ try @info "Relative trailing_edge height: $(upper - lower)" @info "Reynolds number for flying speed of $kite_speed is $reynolds_number" - serialize(polar_path, (cl_interp, cd_interp, cm_interp)) + serialize(polar_path, (alphas, d_trailing_edge_angles, cl_matrix, cd_matrix, cm_matrix)) catch e @info "Removing processes" diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 27483333..48e77cee 100644 --- a/src/VortexStepMethod.jl +++ b/src/VortexStepMethod.jl @@ -10,7 +10,8 @@ using ControlPlots using Measures using LaTeXStrings using NonlinearSolve -using Interpolations: linear_interpolation, Line, Extrapolation +using Interpolations +using Interpolations: linear_interpolation, Line, Extrapolation, FilledExtrapolation using Serialization using SharedArrays diff --git a/src/kite_geometry.jl b/src/kite_geometry.jl index 51cb6d9f..5c5c0665 100644 --- a/src/kite_geometry.jl +++ b/src/kite_geometry.jl @@ -241,7 +241,7 @@ mutable struct KiteWing <: AbstractWing end @info "Loding polars and kite info from $polar_path and $info_path" - (cl_interp, cd_interp, cm_interp) = deserialize(polar_path) + (alpha_range, beta_range, cl_matrix::Matrix, cd_matrix::Matrix, cm_matrix::Matrix) = deserialize(polar_path) (center_of_mass, inertia_tensor, circle_center_z, radius, gamma_tip, le_interp, te_interp, area_interp) = deserialize(info_path) @@ -249,7 +249,7 @@ mutable struct KiteWing <: AbstractWing # Create sections sections = Section[] for gamma in range(-gamma_tip, gamma_tip, n_sections) - aero_input = (:interpolations, (cl_interp, cd_interp, cm_interp)) + aero_input = (:polar_data, (alpha_range, beta_range, cl_matrix, cd_matrix, cm_matrix)) LE_point = [0.0, 0.0, circle_center_z] .+ [le_interp(gamma), sin(gamma) * radius, cos(gamma) * radius] if !isapprox(alpha, 0.0) local_y_vec = [0.0, sin(-gamma), cos(gamma)] × [1.0, 0.0, 0.0] diff --git a/src/panel.jl b/src/panel.jl index f21715d7..87db1ec9 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -1,5 +1,9 @@ using LinearAlgebra +# static types for interpolations +const I1 = Interpolations.FilledExtrapolation{Float64, 1, Interpolations.GriddedInterpolation{Float64, 1, Vector{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64} +const I2 = Interpolations.FilledExtrapolation{Float64, 2, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64} + """ Panel @@ -36,9 +40,9 @@ mutable struct Panel cl_coefficients::Vector{Float64} cd_coefficients::Vector{Float64} cm_coefficients::Vector{Float64} - cl_interp::Function - cd_interp::Function - cm_interp::Function + cl_interp::Union{Nothing, I1, I2} + cd_interp::Union{Nothing, I1, I2} + cm_interp::Union{Nothing, I1, I2} aero_center::MVec3 control_point::MVec3 bound_point_1::MVec3 @@ -82,45 +86,65 @@ mutable struct Panel # Initialize aerodynamic properties cl_coeffs, cd_coeffs, cm_coeffs = zeros(3), zeros(3), zeros(3) - cl_interp, cd_interp, cm_interp = (α::Float64)->0.0, (α::Float64)->0.0, (α::Float64)->0.0 + # Calculate width + width = norm(bound_point_2 - bound_point_1) + + # Initialize filaments + filaments = [ + BoundFilament(bound_point_2, bound_point_1), + BoundFilament(bound_point_1, TE_point_1), + BoundFilament(TE_point_2, bound_point_2) + ] + + cl_interp, cd_interp, cm_interp = nothing, nothing, nothing + if aero_model === :lei_airfoil_breukels cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2) + elseif aero_model === :polar_data aero_1 = section_1.aero_input[2] aero_2 = section_2.aero_input[2] - if size(aero_1) != size(aero_2) + if !all(size.(aero_1) .== size.(aero_2)) throw(ArgumentError("Polar data must have same shape")) end - polar_data = (aero_1 + aero_2) / 2 - cl_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,2]; extrapolation_bc=NaN) - cd_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,3]; extrapolation_bc=NaN) - cm_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,4]; extrapolation_bc=NaN) - cl_interp = (α::Float64) -> cl_interp_(α) - cd_interp = (α::Float64) -> cd_interp_(α) - cm_interp = (α::Float64) -> cm_interp_(α) - elseif aero_model === :interpolations - cl_left, cd_left, cm_left = section_1.aero_input[2] - cl_right, cd_right, cm_right = section_2.aero_input[2] - cl_interp = cl_left === cl_right ? (α::Float64) -> cl_left(α, 0.0) : - (α::Float64) -> 0.5cl_left(α, 0.0) + 0.5cl_right(α, 0.0) # TODO: add trailing edge deflection - cd_interp = cd_left === cd_right ? (α::Float64) -> cd_left(α, 0.0) : - (α::Float64) -> 0.5cd_left(α, 0.0) + 0.5cd_right(α, 0.0) - cm_interp = cm_left === cm_right ? (α::Float64) -> cm_left(α, 0.0) : - (α::Float64) -> 0.5cm_left(α, 0.0) + 0.5cm_right(α, 0.0) + + if length(aero_1) == 4 + !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations." + + polar_data = ( + Vector{Float64}((aero_1[2] + aero_2[2]) / 2), + Vector{Float64}((aero_1[3] + aero_2[3]) / 2), + Vector{Float64}((aero_1[4] + aero_2[4]) / 2) + ) + alphas = Vector{Float64}(aero_1[1]) + + cl_interp = linear_interpolation(alphas, polar_data[1]; extrapolation_bc=NaN) + cd_interp = linear_interpolation(alphas, polar_data[2]; extrapolation_bc=NaN) + cm_interp = linear_interpolation(alphas, polar_data[3]; extrapolation_bc=NaN) + + elseif length(aero_1) == 5 + !all(isapprox.(aero_1[1], aero_2[1])) && @error "Make sure you use the same alpha range for all your interpolations." + !all(isapprox.(aero_1[2], aero_2[2])) && @error "Make sure you use the same beta range for all your interpolations." + + polar_data = ( + Matrix{Float64}((aero_1[3] + aero_2[3]) / 2), + Matrix{Float64}((aero_1[4] + aero_2[4]) / 2), + Matrix{Float64}((aero_1[5] + aero_2[5]) / 2) + ) + alphas = Vector{Float64}(aero_1[1]) + betas = Vector{Float64}(aero_1[2]) + + cl_interp = linear_interpolation((alphas, betas), polar_data[1]; extrapolation_bc=NaN) + cd_interp = linear_interpolation((alphas, betas), polar_data[2]; extrapolation_bc=NaN) + cm_interp = linear_interpolation((alphas, betas), polar_data[3]; extrapolation_bc=NaN) + else + throw(ArgumentError("Polar data in wrong format: $aero_1")) + end + elseif !(aero_model === :inviscid) throw(ArgumentError("Unsupported aero model: $aero_model")) end - - # Calculate width - width = norm(bound_point_2 - bound_point_1) - - # Initialize filaments - filaments = [ - BoundFilament(bound_point_2, bound_point_1), - BoundFilament(bound_point_1, TE_point_1), - BoundFilament(TE_point_2, bound_point_2) - ] new( TE_point_1, LE_point_1, TE_point_2, LE_point_2, @@ -262,8 +286,14 @@ function calculate_cl(panel::Panel, alpha::Float64)::Float64 end elseif panel.aero_model === :inviscid cl = 2π * alpha - elseif panel.aero_model === :polar_data || panel.aero_model === :interpolations - cl = panel.cl_interp(alpha)::Float64 + elseif panel.aero_model === :polar_data + if isa(panel.cl_interp, I1) + cl = panel.cl_interp(alpha)::Float64 + elseif isa(panel.cl_interp, I2) + cl = panel.cl_interp(alpha, 0.0)::Float64 + else + throw(ArgumentError("cl_interp is $(panel.cl_interp)")) + end else throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end @@ -284,9 +314,14 @@ function calculate_cd_cm(panel::Panel, alpha::Float64) if abs(alpha) > (π/9) # Outside ±20 degrees cd = 2 * sin(alpha)^3 end - elseif panel.aero_model === :polar_data || panel.aero_model === :interpolations - cd = panel.cd_interp(alpha)::Float64 - cm = panel.cm_interp(alpha)::Float64 + elseif panel.aero_model === :polar_data + if isa(panel.cd_interp, I1) + cd = panel.cd_interp(alpha)::Float64 + cm = panel.cm_interp(alpha)::Float64 + elseif isa(panel.cd_interp, I2) + cd = panel.cd_interp(alpha, 0.0)::Float64 + cm = panel.cm_interp(alpha, 0.0)::Float64 + end elseif !(panel.aero_model === :inviscid) throw(ArgumentError("Unsupported aero model: $(panel.aero_model)")) end diff --git a/src/solver.jl b/src/solver.jl index 0973d702..4257c228 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -147,7 +147,7 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n ) # Try again with reduced relaxation factor if not converged if !converged && relaxation_factor > 1e-3 - @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)" + log && @warn "Running again with half the relaxation_factor = $(relaxation_factor/2)" converged, gamma_new, alpha_array, v_a_array = gamma_loop( solver, body_aero, @@ -158,7 +158,8 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n y_airf_array, z_airf_array, panels, - relaxation_factor/2 + relaxation_factor/2; + log ) end diff --git a/src/wing_geometry.jl b/src/wing_geometry.jl index fa21d422..c185897a 100644 --- a/src/wing_geometry.jl +++ b/src/wing_geometry.jl @@ -68,7 +68,8 @@ Add a new section to the wing. - aero_input: Can be: - :inviscid - :`lei_airfoil_breukels` - - (:interpolations, (`cl_interp`, `cd_interp`, `cm_interp`)) + - (:polar_data, (`alpha_range`, `cl_vector`, `cd_vector`, `cm_vector`)) + - (:polar_data, (`alpha_range`, `beta_range`, `cl_matrix`, `cd_matrix`, `cm_matrix`)) """ function add_section!(wing::Wing, LE_point::Vector{Float64}, TE_point::Vector{Float64}, aero_input) @@ -157,34 +158,6 @@ function refine_aerodynamic_mesh(wing::AbstractWing) end end -""" - interpolate_to_common_alpha(alpha_common::Vector{Float64}, - alpha_orig::Vector{Float64}, - CL_orig::Vector{Float64}, - CD_orig::Vector{Float64}, - CM_orig::Vector{Float64}) - -Interpolate aerodynamic coefficients to a common alpha range. -""" -function interpolate_to_common_alpha(alpha_common, - alpha_orig, - CL_orig, - CD_orig, - CM_orig) - - # Create interpolation objects - itp_CL = linear_interpolation(alpha_orig, CL_orig) - itp_CD = linear_interpolation(alpha_orig, CD_orig) - itp_CM = linear_interpolation(alpha_orig, CM_orig) - - # Evaluate at common alpha points - CL_common = itp_CL.(alpha_common) - CD_common = itp_CD.(alpha_common) - CM_common = itp_CM.(alpha_common) - - return CL_common, CD_common, CM_common -end - """ calculate_new_aero_input(aero_input, section_index::Int, @@ -212,46 +185,38 @@ function calculate_new_aero_input(aero_input, polar_right = aero_input[section_index + 1][2] # Unpack polar data - @views begin - alpha_left, CL_left, CD_left, CM_left = ( - polar_left[:, i] for i in 1:4 - ) - alpha_right, CL_right, CD_right, CM_right = ( - polar_right[:, i] for i in 1:4 - ) - end - - # Create common alpha array - alpha_common = sort(unique(vcat(alpha_left, alpha_right))) - - # Interpolate both polars - CL_left_common, CD_left_common, CM_left_common = interpolate_to_common_alpha( - alpha_common, alpha_left, CL_left, CD_left, CM_left - ) - CL_right_common, CD_right_common, CM_right_common = interpolate_to_common_alpha( - alpha_common, alpha_right, CL_right, CD_right, CM_right - ) - - # Weighted interpolation - CL_interp = CL_left_common .* left_weight .+ CL_right_common .* right_weight - CD_interp = CD_left_common .* left_weight .+ CD_right_common .* right_weight - CM_interp = CM_left_common .* left_weight .+ CM_right_common .* right_weight - - return (:polar_data, hcat(alpha_common, CD_interp, CL_interp, CM_interp)) + if length(polar_left) == 4 + alpha_left, CL_left, CD_left, CM_left = polar_left + alpha_right, CL_right, CD_right, CM_right = polar_right - elseif model_type === :interpolations - cl_left, cd_left, cm_left = aero_input[section_index][2] - cl_right, cd_right, cm_right = aero_input[section_index + 1][2] - - cl_interp = cl_left === cl_right ? (α, β) -> cl_left[α, β] : - (α, β) -> left_weight * cl_left[α, β] + right_weight * cl_right[α, β] - cd_interp = cd_left === cd_right ? (α, β) -> cd_left[α, β] : - (α, β) -> left_weight * cd_left[α, β] + right_weight * cd_right[α, β] - cm_interp = cm_left === cm_right ? (α, β) -> cm_left[α, β] : - (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β] + # Create common alpha array + !all(isapprox.(alpha_left, alpha_right)) && @error "Make sure you use the same alpha range for all your interpolations." + !isa(CL_right, AbstractVector) && @error "Provide polar data in the correct format: (alpha, cl, cd, cm)" + + # Weighted interpolation + CL_data = CL_left .* left_weight .+ CL_right .* right_weight + CD_data = CD_left .* left_weight .+ CD_right .* right_weight + CM_data = CM_left .* left_weight .+ CM_right .* right_weight + + return (:polar_data, (alpha_left, CL_data, CD_data, CM_data)) + + elseif length(polar_left) == 5 + alpha_left, beta_left, CL_left, CD_left, CM_left = polar_left + alpha_right, beta_right, CL_right, CD_right, CM_right = polar_right + + # Create common alpha array + !all(isapprox.(alpha_left, alpha_right)) && @error "Make sure you use the same alpha range for all your interpolations." + !all(isapprox.(beta_left, beta_right)) && @error "Make sure you use the same alpha range for all your interpolations." + !isa(CL_right, AbstractMatrix) && @error "Provide polar data in the correct format: (alpha, beta, cl, cd, cm)" + + # Weighted interpolation + CL_data = CL_left .* left_weight .+ CL_right .* right_weight + CD_data = CD_left .* left_weight .+ CD_right .* right_weight + CM_data = CM_left .* left_weight .+ CM_right .* right_weight + + return (:polar_data, (alpha_left, beta_left, CL_data, CD_data, CM_data)) + end - return (:interpolations, (cl_interp, cd_interp, cm_interp)) - elseif model_type === :lei_airfoil_breukels tube_diameter_left = aero_input[section_index][2][1] tube_diameter_right = aero_input[section_index + 1][2][1] diff --git a/test/bench.jl b/test/bench.jl index aa8d4840..fdf3a459 100644 --- a/test/bench.jl +++ b/test/bench.jl @@ -86,25 +86,45 @@ using LinearAlgebra y_airf_array[i, :] .= panel.y_airf z_airf_array[i, :] .= panel.z_airf end + + n_angles = 5 + alphas = range(-deg2rad(10), deg2rad(10), n_angles) + cls = [2π * α for α in alphas] + cds = [0.01 + 0.05 * α^2 for α in alphas] + cms = [-0.1 * α for α in alphas] + for model in models - solver = Solver( - aerodynamic_model_type=model - ) - result = @benchmark gamma_loop( - $solver, - $body_aero, - $gamma_new, - $va_array, - $chord_array, - $x_airf_array, - $y_airf_array, - $z_airf_array, - $body_aero.panels, - 0.5; - log = false - ) samples = 1 evals = 1 - @test result.allocs ≤ 100 - @info "Model: $model \t Allocations: $(result.allocs) Memory: $(result.memory)" + for aero_model in [:inviscid, (:polar_data, (alphas, cls, cds, cms))] + wing = Wing(n_panels, spanwise_panel_distribution=:linear) + add_section!(wing, + [0.0, span/2, 0.0], # Left tip LE + [chord, span/2, 0.0], # Left tip TE + aero_model) + add_section!(wing, + [0.0, -span/2, 0.0], # Right tip LE + [chord, -span/2, 0.0], # Right tip TE + aero_model) + body_aero = BodyAerodynamics([wing]) + + solver = Solver( + aerodynamic_model_type=model + ) + result = @benchmark gamma_loop( + $solver, + $body_aero, + $gamma_new, + $va_array, + $chord_array, + $x_airf_array, + $y_airf_array, + $z_airf_array, + $body_aero.panels, + 0.5; + log = false + ) samples = 1 evals = 1 + @test result.allocs ≤ 100 + @info "Model: $model \t Aero_model: $aero_model \t Allocations: $(result.allocs) Memory: $(result.memory)" + end end end @@ -121,16 +141,18 @@ using LinearAlgebra va_unit_array = zeros(n_panels, 3) reference_point = zeros(3) - # Fill arrays with data - for (i, panel) in enumerate(body_aero.panels) - chord_array[i] = panel.chord - x_airf_array[i, :] .= panel.x_airf - y_airf_array[i, :] .= panel.y_airf - z_airf_array[i, :] .= panel.z_airf - va_array[i, :] .= panel.va - end - + + # # Fill arrays with data + # for (i, panel) in enumerate(body_aero.panels) + # chord_array[i] = panel.chord + # x_airf_array[i, :] .= panel.x_airf + # y_airf_array[i, :] .= panel.y_airf + # z_airf_array[i, :] .= panel.z_airf + # va_array[i, :] .= panel.va + # end + set_va!(body_aero, vel_app) results = @MVector zeros(3) + result = @benchmark calculate_results( $body_aero, $gamma, diff --git a/test/test_kite_geometry.jl b/test/test_kite_geometry.jl index bf962c2c..3b131853 100644 --- a/test/test_kite_geometry.jl +++ b/test/test_kite_geometry.jl @@ -100,14 +100,9 @@ using Serialization end end - # Create interpolations - cl_interp = extrapolate(scale(interpolate(cl_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN) - cd_interp = extrapolate(scale(interpolate(cd_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN) - cm_interp = extrapolate(scale(interpolate(cm_matrix, BSpline(Linear())), alphas, d_trailing_edge_angles), NaN) - # Serialize polar data polar_path = test_dat_path[1:end-4] * "_polar.bin" - serialize(polar_path, (cl_interp, cd_interp, cm_interp)) + serialize(polar_path, (alphas, d_trailing_edge_angles, cl_matrix, cd_matrix, cm_matrix)) # Create and serialize obj file faces = [[i, i+1, i+2] for i in 1:2:length(vertices)-2] diff --git a/test/test_panel.jl b/test/test_panel.jl index 2b06f308..404a743f 100644 --- a/test/test_panel.jl +++ b/test/test_panel.jl @@ -115,16 +115,22 @@ end # Generate mock polar data alphas = deg2rad.(-10:1:25) n_points = length(alphas) - polar_data = zeros(n_points, 4) - big_polar_data = zeros(n_points, 4) + polar_data = [zeros(n_points) for _ in 1:4] + big_polar_data = [zeros(n_points) for _ in 1:4] # Fill polar data with realistic values for (i, alpha) in enumerate(alphas) cd = 0.015 + 0.015 * abs(alpha/10)^1.5 # Drag increases with angle cl = 0.1 * alpha + 0.01 * alpha^2 * exp(-alpha/20) # Lift with stall behavior cm = -0.02 * alpha # Linear pitching moment - polar_data[i,:] = [alpha, cl, cd, cm] - big_polar_data[i,:] = [alpha, 1.1cl, 1.1cd, 1.1cm] + polar_data[1][i] = alpha + polar_data[2][i] = cl + polar_data[3][i] = cd + polar_data[4][i] = cm + big_polar_data[1][i] = alpha + big_polar_data[2][i] = 1.1cl + big_polar_data[3][i] = 1.1cd + big_polar_data[4][i] = 1.1cm end # Create two sections with slightly different polar data @@ -142,9 +148,9 @@ end cd, cm = calculate_cd_cm(panel, alpha) # Calculate expected values through interpolation - expected_cl = linear_interpolation(polar_data[:,1], polar_data[:,2], extrapolation_bc=Line())(alpha) - expected_cd = linear_interpolation(polar_data[:,1], polar_data[:,3], extrapolation_bc=Line())(alpha) - expected_cm = linear_interpolation(polar_data[:,1], polar_data[:,4], extrapolation_bc=Line())(alpha) + expected_cl = linear_interpolation(polar_data[1], polar_data[2], extrapolation_bc=Line())(alpha) + expected_cd = linear_interpolation(polar_data[1], polar_data[3], extrapolation_bc=Line())(alpha) + expected_cm = linear_interpolation(polar_data[1], polar_data[4], extrapolation_bc=Line())(alpha) # Average with slightly different data (1.1 factor) expected_cl = (expected_cl + 1.1 * expected_cl) / 2