text:
### A Pluto.jl notebook ###
# v0.14.3

using Markdown
using InteractiveUtils

# ╔═╡ c695e3e0-2f3b-11eb-027f-ed0248d96007
using Plots

# ╔═╡ bd1b2930-7458-11eb-3483-71e2935f5945
begin
	using PlutoUI
	almost(text) = Markdown.MD(Markdown.Admonition("warning", "Almost there!", [text]))
	keep_working(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("danger", "Keep working on it!", [text]));
	alert(text=md"Alert!") = Markdown.MD(Markdown.Admonition("danger", "Alert!", [text]));
	definition(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("warning", "Definition", [text]));
	exercise(text=md"The answer is not quite right.") = Markdown.MD(Markdown.Admonition("correct", "Exercise", [text]));
	correct(text=md"Great! You got the right answer! Let's move on to the next section.") = Markdown.MD(Markdown.Admonition("correct", "Got it!", [text]));
	hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));
	warning(text=md"Warning.") = Markdown.MD(Markdown.Admonition("warning", "Warning!", [text]));
	tip(text=md"Tip.") = Markdown.MD(Markdown.Admonition("warning", "Tip", [text]));
end;

# ╔═╡ 1b0fda70-7433-11eb-3fd1-5bde27831603
md"""
# Milestone 2: Matching Chart + Drag Polar estimation
"""

# ╔═╡ 4450d52a-bcc4-4a75-a56d-479e81ad53c8
plotlyjs()

# ╔═╡ 59bc0250-7463-11eb-0880-7348a2b01ad8
md"### Aircraft Parameters
The wing span length is set to $b = 8.1m$ and the wing area $S = 7m^2$.
Here we consider an aircraft with aspect ratio $AR = 8.67$, and an Oswald span
efficiency factor $e = 0.8$.
The zero-lift drag coefficient $C_{D_0}$ can be calculated by various methods, such
as regression, wetted area estimation, drag build-up, etc. Here we use a value
$C_{D_0} = 0.013$."

# ╔═╡ dca970c0-2f3b-11eb-21c6-9b2ceb7b62ca
begin #clean
CD0 = 13e-3
AR = 7
e = 0.8
end;

# ╔═╡ 6e4f4e40-7457-11eb-3561-f7458f6c733a
md"""
### Induced Drag Coefficient

```math
k = \frac{1}{\pi e AR}
```
"""

# ╔═╡ cb4cb52e-2f3b-11eb-0121-c92bc4b2b950
induced_drag_coefficient(e, AR) = 1 / (π * e * AR)

# ╔═╡ e19a16c0-2f3b-11eb-204b-a1d79f7a31ad
k = induced_drag_coefficient(e, AR)

# ╔═╡ 80ca4b10-7457-11eb-27d2-dd5fd4e77383
md"""
### Parabolic Drag Polar

```math
C_D = C_{D_0} + k\left(C_L - C_{L_{\alpha = 0}}\right)^2, \quad k = \frac{1}{\pi e AR}
```
"""

# ╔═╡ d701fde0-2f3b-11eb-345f-955af341b009
drag_polar(CD0, k, CL, CL0 = 0.) = CD0 + k * (CL - CL0)^2

# ╔═╡ 3aaa03d0-746d-11eb-2b50-dd705284a42f
md"Let's define a range of lift coefficients over which we would like to see the variation of the drag polar."

# ╔═╡ 46d54fa0-2f3c-11eb-2c81-c972cc10fc41
cls = -1.5:0.05:1.5

# ╔═╡ dfe1e4c0-746e-11eb-267e-3546c9416205
md"Here we consider an cambered airfoil with NACA $2412$, which has $C_{L_{\alpha = 0}} = 0.3$."

# ╔═╡ 19df47e0-746e-11eb-2b1b-a9b4a09db0a5
CL0_1 = 0.3

# ╔═╡ 3707e260-7463-11eb-17f7-073926ca3929
cds = drag_polar.(CD0, k, cls, CL0_1);

# ╔═╡ 17dc7840-2f3c-11eb-26d4-67e3a57e9730
test_plot = plot(cds, cls,
	 			 label = "Clean",
				 xlabel = "CD", ylabel = "CL", title = "Drag Polar")

# ╔═╡ 45a4c000-7476-11eb-3b15-7795df62140c
# savefig(test_plot, "Drag Polar")

# ╔═╡ 06d40dd1-1244-40df-9ac4-e81ea7f37c1a
begin #for cruise
CD0_2 = 0.014
CL0_2 = 0.093
cds_2 = drag_polar.(CD0_2, k, cls, CL0_2)
end;

# ╔═╡ 75c452c3-cad2-472e-9949-108aa62b53e5
begin #landing
CD0_3 = 0.013
CL0_3 = 0.04
cds_3 = drag_polar.(CD0_3, k, cls, CL0_3)
end;

# ╔═╡ 3ebcb070-746e-11eb-2e5c-8129abd1179e
begin
	plot(cds, cls,
	label = "Clean",
	xlabel = "CD", ylabel = "CL", title = "Drag Polar")
	plot!(cds_2, cls, marker = :dot,
	label = "CD0 = $CD0_2, CL0 = $CL0_2")
	plot!(cds_3, cls, marker = :dot,
	label = "CD0 = $CD0_3, CL0 = $CL0_3")

end

# ╔═╡ 4eb37530-7433-11eb-3149-afb0f1d19b2c
md"""
## Constraints Analysis
The constraints of the initial sizing procedure are determined by equating the thrust-to-weight ratio $(T/W)$ to a function of the wing loading $(W/S)$ for different flight conditions.

```math
(T/W) = f(W/S)
```

The objective is to obtain the maximum $(W/S)$ and the minimum $(T/W)$ to obtain the optimum size and propulsion. These will allow determination of the wing area $S_w$ and required thrust, once the MTOW has been determined.

*Note*: We use $S$ and $S_w$ interchangeably for brevity, differentiating in context.

Here we define a range of wing loading values we consider for our aircraft.
"""

# ╔═╡ 614a1050-7433-11eb-0821-41c5a3438a0f
wing_loadings = 5:1:2000;

# ╔═╡ 9831b280-7433-11eb-30f9-7b6994bc6653
md"""
### Power Loading

For engines or motors with propellers, the relevant conversion accounting for speed and efficiency of the propeller is as follows:

\begin{align}
	(P/W) = \frac{(T/W) V}{\eta_\text{prop}}
\end{align}
"""

# ╔═╡ 9be5449e-7433-11eb-103d-8f4adda19729
η_prop_FW = 0.8;

# ╔═╡ 93acedb0-7433-11eb-23cf-e76771b143e2
power_loading(tbW, V, η_p) = tbW * V / η_p

# ╔═╡ 31f691a0-743a-11eb-2d56-a5a00e193c49
md"""### Stall Speed

Dynamic pressure:

```math
q = \frac{1}{2} \rho V^2
```

Stall speed:
```math
(W/S)^{FW}_\text{stall} = \frac{1}{2} \rho V^2_\text{stall} C_{L_\max}
```
"""

# ╔═╡ b979b830-2f3b-11eb-2c1f-95035ea2a15b
dynamic_pressure(ρ, V) = 1/2 * ρ * V^2

# ╔═╡ 3d9c49f2-743a-11eb-0cd9-a32988a847f4
wing_loading_stall_speed(V_stall, CL_max, ρ = 1.225) = dynamic_pressure(ρ, V_stall) * CL_max

# ╔═╡ 338b05c0-7465-11eb-1ef1-75134695ef8e
V_stall = 31 #why did they put 97.5?, 61 knots in prev. FAR23 req.

# ╔═╡ 34f56b30-7465-11eb-35b7-b3147697cc9d
CL_max = 1.8

# ╔═╡ 45ed0ae0-743a-11eb-34b1-8dc823f25349
n = 20

# ╔═╡ 426a3e30-7465-11eb-3956-f7d1b0e8630f
stalls = fill(wing_loading_stall_speed(V_stall, CL_max), n);

# ╔═╡ 1c9f3d20-746c-11eb-15c5-5f0e9ada6740
md"For the small-scale aircraft, the $(T/W)$ ratios may exceed $1$ with oversized
motors."

# ╔═╡ 30610740-746b-11eb-103a-97a1a677aa08
tbWs_stall = range(0, 5.; length = n)

# ╔═╡ 426cd640-7465-11eb-13b6-cb0a8c77df0a
pbWs_stall = power_loading.(tbWs_stall, V_stall, CL_max);

# ╔═╡ bc98aa10-8bc2-11eb-3761-275f8c315941
pbWs_stall

# ╔═╡ 20ee4160-746b-11eb-0707-61ac484435ea
plot(stalls, tbWs_stall, label = "Stall Speed",
	 ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²")

# ╔═╡ 2b7cae61-0f9d-4f4a-94e3-318a679d6d26
md"""
### Take-off
"""

# ╔═╡ 00ef41ef-9241-4331-971e-db20e6c68467
takeoff_condition(CL_max,TOP,wing_load)=(1/(CL_max*TOP))*wing_load

# ╔═╡ 401133b0-81db-4419-bf38-78dd836c23d0
begin
	CL_takeoff=2
	TOP=650
	V_to=V_stall*1.1
end;

# ╔═╡ 0c098d7d-a67e-4053-8f07-738e29eee056
tbw_takeoff=takeoff_condition.(CL_takeoff,TOP,wing_loadings)

# ╔═╡ 1177e156-a184-4dda-bb63-81b76c097548
pbw_takeoff=power_loading.(tbw_takeoff, V_to, η_prop_FW)

# ╔═╡ 72b35b03-ff38-4b0d-ac51-65ef7d45cb7f
begin
	Stall_takeoff_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart (Fixed Wing)")
	plot!(stalls, tbWs_stall, label = "Stall Speed")
	#plot!(wing_loadings, tbw_takeoff, label = "Takeoff")
	plot!(wing_loadings,tbw_takeoff,label="takeoff power")
end

# ╔═╡ b5aa6a47-1e0b-4a8f-877b-7e69a998f679
md"""
### Landing
"""

# ╔═╡ e1341b5a-9751-4cb9-954d-947959795225
begin
	CLmax_land=2.5
	safety_marg=0.6
	S_a=183
	g=9.18
	MTOW=2541.309934319981;
	Mland=MTOW-(MTOW*0.31135388287326593/1.1)
	WLW0=Mland/MTOW
	V_land=V_stall*1.3
	S_fl=750
end;

# ╔═╡ 9fe7102e-0c47-4e3e-98ec-c0fa209e1bbb
landing_condition(g, CLmax_land, S_fl, S_a, safety_marg)=((g*CLmax_land)/5)*((S_fl*safety_marg)-S_a)

# ╔═╡ 1970ac8e-2eb8-463f-93d6-ff9f826ef178
landing_condition(g, CLmax_land, S_fl, S_a, safety_marg)

# ╔═╡ 7ed46d06-3118-4ef0-b76f-44e5b1e0857b
landing_sls=(1/WLW0)*landing_condition(g, CLmax_land, S_fl, S_a, safety_marg)

# ╔═╡ 6cb0fe34-ad41-472b-9847-1ae42c3679c9
wbs_landing=fill(landing_sls,n);

# ╔═╡ 82d3f856-20ef-42ef-afa5-1fb6a2f7781c
tbWs_landing = range(0, 5.; length = n);

# ╔═╡ 0c8b36b7-dd69-4be2-ac26-a89b53a7bd30
pbWs_landing = power_loading.(tbWs_landing, V_land, η_prop_FW);

# ╔═╡ 027f830f-9591-4024-81c4-15d3df9a40f4
begin
	Stall_land_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart (Fixed Wing)")
	plot!(stalls, tbWs_stall, label = "Stall Speed")
	plot!(wing_loadings, tbw_takeoff, label = "Takeoff")
	plot!(wbs_landing, tbWs_landing, label = "Landing")
end

# ╔═╡ 13c57067-22a6-4f40-a52b-3c97ca09c22c
stalls;

# ╔═╡ 752bf020-7433-11eb-0393-d3ab1d104f47
md"""
### Cruise

The matching equation for cruise is obtained from the cruise conditions $T = D,~L = W$ and substituting into the drag polar relation.

```math
\begin{align}
	(T/W)^\text{FW}_\text{cruise} = \frac{q C_{D_0}}{(W/S)} + \frac{k}{q}(W/S)
\end{align}
```
"""

# ╔═╡ 7f11a8f0-7433-11eb-20e8-3b8ddaee4060
thrust_to_weight_fw_cruise(q, k, CD0, wing_loading) = ((1/(0.909/1.225))^0.6)*0.990016*(q * CD0 / (wing_loading) + k / q * wing_loading)

# ╔═╡ 85c91220-7468-11eb-114d-51f4af036910
begin
	V_cruise = 77 #my input:77m/s
	ρ = 0.905 #1.31
	density_ratio_cruise=ρ/1.225
	q = dynamic_pressure(ρ, V_cruise)
	k_cruise=induced_drag_coefficient(0.65,AR)
	weightFrac_cruise=0.990016
end;

# ╔═╡ 839d4b90-7433-11eb-2273-79686fbd7d7e
begin
	tbWs_cruise = thrust_to_weight_fw_cruise.(q, k, CD0_2, wing_loadings)
	pbWs_cruise = power_loading.(tbWs_cruise, V_cruise, η_prop_FW)
end;

# ╔═╡ dc76e9af-c578-4f8c-970c-800e2a7877fa
cruise_condition(wing_load,q,CD0,k)=(q*CD0/wing_load)+(k*wing_load/q)

# ╔═╡ cef05bca-8384-4bb0-8df8-3a8b2cdfd3bf
tbw_cruise=1/density_ratio_cruise^0.6 * weightFrac_cruise* cruise_condition.(wing_loadings,q,CD0,k_cruise)

# ╔═╡ 39e74b9b-abe4-4c13-9cb0-38851886907f
pbw_cruise=power_loading.(tbw_cruise, V_cruise, η_prop_FW);

# ╔═╡ 1cc66961-edd1-480b-951d-2eb96da19856
plot(wing_loadings, tbw_cruise, label = "Cruise",
	 ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²")

# ╔═╡ 59cc1dfd-2394-448d-8120-1e0b702715eb
begin
	Stall_land_cruise_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart (Fixed Wing)")
	plot!(stalls, tbWs_stall, label = "Stall Speed")
	plot!(wing_loadings, tbw_takeoff, label = "Takeoff")
	plot!(wbs_landing, tbWs_landing, label = "Landing")
	plot!(wing_loadings,tbw_cruise,label="Cruise")
end

# ╔═╡ efbbc3f0-7439-11eb-161d-3166990ea8f1
md"""
### Climb
"""

# ╔═╡ fcaae67d-f693-42e9-a4d9-8f1cf9b45a6e
climb_condition(k_s, CD0, CL_max, K, G)=((k_s^2)*CD0/CL_max)+((CL_max/(k_s^2))*k)+G;

# ╔═╡ 00a688f5-fca3-4b6a-967d-0245917cc0da
function thrust_corr_climb(k_s, CD0, CL_max, K, G, weight_factor, MCT=false)

	MCT_factor=ifelse(MCT, 1/0.94, 1)

	(1/0.8)*MCT_factor*weight_factor*climb_condition(k_s, CD0, CL_max, K, G)
end

# ╔═╡ 44b52b65-ff5f-48df-89d5-e40b368deece
begin
	CD0_climb=CD0
	CL_max_climb=2.2
end;

# ╔═╡ 26b2336a-7a75-4ff6-842b-83c0e7d257cb
md"""
#### Takeoff climb
"""

# ╔═╡ 533111bd-b958-43c6-83ba-b3c19973e5a5
begin
	k_s_takeoff=1.1
	G_takeoff=0.083
	CD0_takeoff=CD0
	k_takeoff=induced_drag_coefficient(0.85, AR)
	weight_fact1=1
	MCT_takeoff=false

	takeoff_climbs=fill(thrust_corr_climb(k_s_takeoff,CD0_takeoff,CL_max_climb,k_takeoff,G_takeoff,weight_fact1,MCT_takeoff),length(wing_loadings))
end;

# ╔═╡ ef4a6911-1ed0-4ed3-b6ea-16d9ca443da1
md"""
#### Transition climb
"""

# ╔═╡ e3eb9744-6ae1-4a7b-ae13-cea30c5d5d48
begin
	k_s_transition=1.1
	G_transition=0
	CD0_transition=CD0
	k_transition=induced_drag_coefficient(0.75,AR)
	weight_fact2=0.99
	MCT_transition=false

	transition_climb=fill(thrust_corr_climb(k_s_transition,CD0_transition,CL_max_climb,k_transition,G_transition,weight_fact2,MCT_transition),length(wing_loadings))
end;

# ╔═╡ 5ca80935-8f66-45c7-ab31-b61dadc0168d
md"""
#### Second climb
"""

# ╔═╡ 30570a80-6945-4b42-bf61-e35d4cde78a3
begin
	k_s_second=1.1
	G_second=0.083
	CD0_second=CD0
	k_second=induced_drag_coefficient(0.85,AR)
	weight_fact3=0.98
	MCT_second=false

	second_climb=fill(thrust_corr_climb(k_s_second,CD0_second,CL_max_climb,k_second,G_second,weight_fact3,MCT_second),length(wing_loadings))
end;

# ╔═╡ 34eb5b06-8be0-4c71-8d65-be7f4d9eaff9
md"""
#### Enroute climb
"""

# ╔═╡ 673cec13-25aa-41be-a98e-80ff4cd15d6f
begin
	k_s_enroute=1.2
	G_enroute=0.083
	CD0_enroute=CD0
	k_enroute=induced_drag_coefficient(0.85,AR)
	weight_fact4=1
	MCT_enroute=true

	enroute_climb=fill(thrust_corr_climb(k_s_enroute,CD0_enroute,CL_max_climb,k_enroute,G_enroute,weight_fact4,MCT_enroute),length(wing_loadings))
end;

# ╔═╡ 71780241-ffb3-4115-8821-4f4843e372ec
md"""
#### Balked landing climb
"""

# ╔═╡ d0caf29c-a624-4612-8369-3175be159e98
begin
	k_s_balked=1.5
	G_balked=0.03
	CD0_balked=CD0
	k_balked=induced_drag_coefficient(0.75,AR)
	weight_fact5=1
	MCT_balked=false

	balked_climb=fill(thrust_corr_climb(k_s_balked,CD0_balked,CL_max_climb,k_balked,G_balked,weight_fact5,MCT_balked),length(wing_loadings))
end;

# ╔═╡ d27f356b-f728-40f6-9af1-0ad6e3943c8e
begin
	climb_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart (climb)")
	plot!(wing_loadings,takeoff_climbs, label="Climb 1")
	plot!(wing_loadings,transition_climb, label="Climb 2")
	plot!(wing_loadings,second_climb, label="Climb 3")
	plot!(wing_loadings,enroute_climb, label="Climb 4")
	plot!(wing_loadings,balked_climb, label="Climb 5")
end

# ╔═╡ 6baa7798-7b9f-48b1-8dc3-3bbc4b82980d
begin
	match_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart (Fixed Wing)")
	plot!(stalls, tbWs_stall, label = "Stall Speed")
	plot!(wing_loadings, tbw_takeoff, label = "Takeoff")
	plot!(wbs_landing, tbWs_landing, label = "Landing")
	plot!(wing_loadings,tbw_cruise,label="Cruise")
	plot!(wing_loadings,takeoff_climbs, label="Climb 1")
	plot!(wing_loadings,transition_climb, label="Climb 2")
	plot!(wing_loadings,second_climb, label="Climb 3")
	plot!(wing_loadings,enroute_climb, label="Climb 4")
	plot!(wing_loadings,balked_climb, label="Climb 5")
end

# ╔═╡ 4f6a3de0-743a-11eb-109b-63f93165bd36
md"""
### Service Ceiling
The service ceiling is the altitude at which the maximum rate of climb is $0.001~m/s$, so we re-use the climb equation with this value.
"""

# ╔═╡ bbe43fc0-746c-11eb-1435-d31af1391f71
begin
	RoC_ceiling = 0.001
	density_ceiling=0.89469763258
	weight_ratio_ceiling=0.718339816991928
	k_ceiling=k_cruise
end;

# ╔═╡ 8bd45fb1-9056-4a68-914c-14cabf344e2e
ceiling_condition(RoC_ceiling,CD0_2, k)=RoC_ceiling*sqrt(2*CD0_2*k)

# ╔═╡ dc2d5cf4-dc53-43bd-a188-ad06c064b270
best_climb_rate(density,wing_load,k,CD0)=sqrt((2/density)*wing_load*sqrt(k/(3*CD0)))

# ╔═╡ c8d0da8e-5fcf-4bed-a5af-3c4fa03ea177
best_VRoC=best_climb_rate.(density_ceiling,wing_loadings,k_ceiling,CD0_2)

# ╔═╡ ed7768ff-e274-4f0e-912e-353c35a720ed
begin
	tbWs_ceiling=fill((((1.225/density_ceiling)^0.6))*weight_ratio_ceiling*ceiling_condition(RoC_ceiling,CD0_2, k_ceiling),length(wing_loadings))
	pbWs_Ceiling=power_loading.(tbWs_ceiling, best_VRoC, η_prop_FW)
end;

# ╔═╡ 0962b6c0-af94-4582-8986-92cea23c030a
begin
	match_chart2 = plot(ylabel = "(T/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart")
	plot!(stalls, tbWs_stall, label = "Stall Speed")
	plot!(wing_loadings, tbw_takeoff, label = "Takeoff")
	plot!(wbs_landing, tbWs_landing, label = "Landing")
	plot!(wing_loadings,tbw_cruise,label="Cruise")
	plot!(wing_loadings,takeoff_climbs, label="Climb 1")
	plot!(wing_loadings,transition_climb, label="Climb 2")
	plot!(wing_loadings,second_climb, label="Climb 3")
	plot!(wing_loadings,enroute_climb, label="Climb 4")
	plot!(wing_loadings,balked_climb, label="Climb 5")
	plot!(wing_loadings,tbWs_ceiling,label="Ceiling")
end

# ╔═╡ 624bebc0-743a-11eb-1832-ddee60ff70b5
md"""
## Matching Charts
"""

# ╔═╡ 66e263d0-743a-11eb-0259-bd958c5897d1
begin
	fw_chart = plot(ylabel = "(P/W), W/N", xlabel = "(W/S), N/m²",
		 			title = "Matching Chart")
	plot!(stalls, pbWs_stall, label = "Stall Speed")
	plot!(wing_loadings, tbw_cruise, label = "Cruise")
	#plot!(wing_loadings, pbWs_fw_climb, label = "Fixed-Wing Climb")
	plot!(wing_loadings, pbWs_Ceiling, label = "Ceiling")
	plot!(wbs_landing, pbWs_landing, label = "Landing")
	plot!(wing_loadings,takeoff_climbs, label="Climb 1")
	plot!(wing_loadings,transition_climb, label="Climb 2")
	plot!(wing_loadings,second_climb, label="Climb 3")
	plot!(wing_loadings,enroute_climb, label="Climb 4")
	plot!(wing_loadings,balked_climb, label="Climb 5")
	plot!(wing_loadings, pbw_takeoff, label = "Takeoff")
end

# ╔═╡ 43e98480-75c0-11eb-1807-f5c025c1b476
md"Now that we have the $(W/S)$ and $(T/W)$ or $(P/W)$ values, we can obtain the reference wing area $S_w$ and required thrust $T$ or power $P$ by using the maximum takeoff weight $W_0$ obtained from the initial weight estimation procedure."

# ╔═╡ Cell order:
# ╟─1b0fda70-7433-11eb-3fd1-5bde27831603
# ╠═4450d52a-bcc4-4a75-a56d-479e81ad53c8
# ╟─59bc0250-7463-11eb-0880-7348a2b01ad8
# ╠═dca970c0-2f3b-11eb-21c6-9b2ceb7b62ca
# ╟─6e4f4e40-7457-11eb-3561-f7458f6c733a
# ╠═cb4cb52e-2f3b-11eb-0121-c92bc4b2b950
# ╠═e19a16c0-2f3b-11eb-204b-a1d79f7a31ad
# ╟─80ca4b10-7457-11eb-27d2-dd5fd4e77383
# ╠═d701fde0-2f3b-11eb-345f-955af341b009
# ╟─3aaa03d0-746d-11eb-2b50-dd705284a42f
# ╠═46d54fa0-2f3c-11eb-2c81-c972cc10fc41
# ╟─dfe1e4c0-746e-11eb-267e-3546c9416205
# ╠═19df47e0-746e-11eb-2b1b-a9b4a09db0a5
# ╠═3707e260-7463-11eb-17f7-073926ca3929
# ╠═c695e3e0-2f3b-11eb-027f-ed0248d96007
# ╠═17dc7840-2f3c-11eb-26d4-67e3a57e9730
# ╠═45a4c000-7476-11eb-3b15-7795df62140c
# ╠═06d40dd1-1244-40df-9ac4-e81ea7f37c1a
# ╠═75c452c3-cad2-472e-9949-108aa62b53e5
# ╠═3ebcb070-746e-11eb-2e5c-8129abd1179e
# ╟─4eb37530-7433-11eb-3149-afb0f1d19b2c
# ╠═614a1050-7433-11eb-0821-41c5a3438a0f
# ╟─9831b280-7433-11eb-30f9-7b6994bc6653
# ╠═9be5449e-7433-11eb-103d-8f4adda19729
# ╠═93acedb0-7433-11eb-23cf-e76771b143e2
# ╟─31f691a0-743a-11eb-2d56-a5a00e193c49
# ╠═b979b830-2f3b-11eb-2c1f-95035ea2a15b
# ╠═3d9c49f2-743a-11eb-0cd9-a32988a847f4
# ╠═338b05c0-7465-11eb-1ef1-75134695ef8e
# ╠═34f56b30-7465-11eb-35b7-b3147697cc9d
# ╠═45ed0ae0-743a-11eb-34b1-8dc823f25349
# ╠═426a3e30-7465-11eb-3956-f7d1b0e8630f
# ╟─1c9f3d20-746c-11eb-15c5-5f0e9ada6740
# ╠═30610740-746b-11eb-103a-97a1a677aa08
# ╠═426cd640-7465-11eb-13b6-cb0a8c77df0a
# ╠═bc98aa10-8bc2-11eb-3761-275f8c315941
# ╠═20ee4160-746b-11eb-0707-61ac484435ea
# ╟─2b7cae61-0f9d-4f4a-94e3-318a679d6d26
# ╠═00ef41ef-9241-4331-971e-db20e6c68467
# ╠═401133b0-81db-4419-bf38-78dd836c23d0
# ╠═0c098d7d-a67e-4053-8f07-738e29eee056
# ╠═1177e156-a184-4dda-bb63-81b76c097548
# ╠═72b35b03-ff38-4b0d-ac51-65ef7d45cb7f
# ╟─b5aa6a47-1e0b-4a8f-877b-7e69a998f679
# ╠═e1341b5a-9751-4cb9-954d-947959795225
# ╠═9fe7102e-0c47-4e3e-98ec-c0fa209e1bbb
# ╠═1970ac8e-2eb8-463f-93d6-ff9f826ef178
# ╠═7ed46d06-3118-4ef0-b76f-44e5b1e0857b
# ╠═6cb0fe34-ad41-472b-9847-1ae42c3679c9
# ╠═82d3f856-20ef-42ef-afa5-1fb6a2f7781c
# ╠═0c8b36b7-dd69-4be2-ac26-a89b53a7bd30
# ╠═027f830f-9591-4024-81c4-15d3df9a40f4
# ╠═13c57067-22a6-4f40-a52b-3c97ca09c22c
# ╟─752bf020-7433-11eb-0393-d3ab1d104f47
# ╠═7f11a8f0-7433-11eb-20e8-3b8ddaee4060
# ╠═85c91220-7468-11eb-114d-51f4af036910
# ╠═839d4b90-7433-11eb-2273-79686fbd7d7e
# ╠═dc76e9af-c578-4f8c-970c-800e2a7877fa
# ╠═cef05bca-8384-4bb0-8df8-3a8b2cdfd3bf
# ╠═39e74b9b-abe4-4c13-9cb0-38851886907f
# ╠═1cc66961-edd1-480b-951d-2eb96da19856
# ╠═59cc1dfd-2394-448d-8120-1e0b702715eb
# ╟─efbbc3f0-7439-11eb-161d-3166990ea8f1
# ╠═fcaae67d-f693-42e9-a4d9-8f1cf9b45a6e
# ╠═00a688f5-fca3-4b6a-967d-0245917cc0da
# ╠═44b52b65-ff5f-48df-89d5-e40b368deece
# ╟─26b2336a-7a75-4ff6-842b-83c0e7d257cb
# ╠═533111bd-b958-43c6-83ba-b3c19973e5a5
# ╟─ef4a6911-1ed0-4ed3-b6ea-16d9ca443da1
# ╠═e3eb9744-6ae1-4a7b-ae13-cea30c5d5d48
# ╟─5ca80935-8f66-45c7-ab31-b61dadc0168d
# ╠═30570a80-6945-4b42-bf61-e35d4cde78a3
# ╟─34eb5b06-8be0-4c71-8d65-be7f4d9eaff9
# ╠═673cec13-25aa-41be-a98e-80ff4cd15d6f
# ╟─71780241-ffb3-4115-8821-4f4843e372ec
# ╠═d0caf29c-a624-4612-8369-3175be159e98
# ╠═d27f356b-f728-40f6-9af1-0ad6e3943c8e
# ╠═6baa7798-7b9f-48b1-8dc3-3bbc4b82980d
# ╟─4f6a3de0-743a-11eb-109b-63f93165bd36
# ╠═bbe43fc0-746c-11eb-1435-d31af1391f71
# ╠═8bd45fb1-9056-4a68-914c-14cabf344e2e
# ╠═dc2d5cf4-dc53-43bd-a188-ad06c064b270
# ╠═c8d0da8e-5fcf-4bed-a5af-3c4fa03ea177
# ╠═ed7768ff-e274-4f0e-912e-353c35a720ed
# ╠═0962b6c0-af94-4582-8986-92cea23c030a
# ╟─624bebc0-743a-11eb-1832-ddee60ff70b5
# ╠═66e263d0-743a-11eb-0259-bd958c5897d1
# ╟─43e98480-75c0-11eb-1807-f5c025c1b476
# ╟─bd1b2930-7458-11eb-3483-71e2935f5945

max_stars_repo_path: Matching Chart.jl

max_stars_repo_name: rahmanfahd/aircraft-design

max_issues_repo_path: Matching Chart.jl

max_issues_repo_name: rahmanfahd/aircraft-design

max_forks_repo_path: Matching Chart.jl

max_forks_repo_name: rahmanfahd/aircraft-design


text:
# Autogenerated wrapper script for Shoco_jll for aarch64-apple-darwin
export libshoco

JLLWrappers.@generate_wrapper_header("Shoco")
JLLWrappers.@declare_library_product(libshoco, "@rpath/libshoco.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libshoco,
        "lib/libshoco.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
module HWfuncapp

	using FastGaussQuadrature # to get chebyshevnodes

	# you dont' have to use PyPlot. I did much of it in PyPlot, hence
	# you will see me often qualify code with `PyPlot.plot` or so.
	#using PyPlot
	using LinearAlgebra
	import ApproXD: getBasis, BSpline
	using Distributions
	using ApproxFun
	using Plots

	export scalemap, unitmap, TT, ChebyT, q1, q2, q3, q7, runall

	ChebyT(x,deg) = cos(acos(x)*deg)
	unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1	#[a,b] -> [-1,1]
	scalemap(x,lb,ub) = 0.5 .* (lb .+ ub) .+ 0.5 .* (ub .- lb) .* x #scaling to [a,b]

	function TT(x,J,lb,ub)
		M = length(x)
		Phi = zeros(M,J)
		for i in 1:M
			for j in 1:J
				Phi[i,j] = ChebyT(unitmap(x[i],lb,ub),j-1)
			end
		end
		return Phi
	end

	#function IPmatCh(n) #Chebyshev interpolation matrix
	#	M = zeros(n,n)
	#	for i in 1:n
	#		for j in 1:n
	#			M[i,j] = cos(((n - i + 0.5)*(j - 1)*(pi))/(n))
	#		end
	#	end
	#end

	#function cHat(n)
	#	return inv(IPmatCh(n)) * f(gch(-3,3,n))
	#end
	#f(x) = x .+ 2x.^2 - exp.(-x)

	#function gch(lb, ub, nnodes)
	#	return scalemap(gausschebyshev(nnodes)[1], lb, ub)
	#end

	function q1(n=15)
		lb = -3
		ub = 3
		n_new = 100
		f(x) = x .+ 2x.^2 - exp.(-x)
		chNodes = scalemap(gausschebyshev(n)[1], lb, ub)
		cIP = inv(TT(chNodes,n,lb,ub))*f(chNodes)
		grid = Vector(range(lb,length = n_new,stop = ub))
		true_new = f(grid)
		IP_new = (TT(grid,n,lb,ub))*cIP
		err = true_new .- IP_new
		p = Plots.plot(layout = 2, dpi = 400)
		Plots.plot!(p[1],grid,true_new,label = "True",lw = 1,linecolor = "black")
		Plots.plot!(p[1],grid,IP_new, label = "IP Approx", lw = 2,linestyle = :dot, linecolor = "red")
		Plots.plot!(p[2],grid,err, label = "error", lw = 1, linecolor = "green")
		# without using PyPlot, just erase the `PyPlot.` part
		Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q1.png"))
		return Dict(:error=>maximum(abs,err))
	end

	function q2(b::Number = 4)
		@assert b > 0
		# use ApproxFun.jl to do the same:
		S = Chebyshev(-b..b)
		ft = Fun(x -> x .+ 2x.^2 - exp.(-x),-b..b)
		n,m = 100,50
		grid = Vector(range(-b, stop = b, length = n))
		true_v = ft.(grid)
		V = zeros(n,m)
		for k = 1:m
    		V[:,k] = Fun(S,[zeros(k-1);1]).(grid)
		end
		fa = Fun(S,V\true_v)
		approx_v = fa.(grid)
		err = approx_v .- true_v
		p = Plots.plot(layout = 2, dpi = 400)
		Plots.plot!(p[1],grid,true_v,label = "True",lw = 1,linecolor = "black")
		Plots.plot!(p[1],grid,approx_v, label = "IP Approx", lw = 2,linestyle = :dot, linecolor = "red")
		Plots.plot!(p[2],grid,err, label = "error", lw = 1, linecolor = "green")
		Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q2.png"))
	end

	function q3(b::Number = 10)
		@assert b > 0
		x = Fun(identity,-b..b)
		f = sin(x^2)
		g = cos(x)
		h = f - g
		rh = roots(h)
		p = Plots.plot(h, label = "h(x)", dpi = 400, lw = 1, linecolor = "black")
		Plots.scatter!(p,rh,h.(rh), label = "roots of h(x)", markercolor = "green", markersize = 5)
		Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q3.png"))
		y = Fun(identity,-10..0)
		hh = sin(y^2) - cos(y)
		hh_int1 = cumsum(hh)
		hh_int = hh_int1 + hh(-10)
		integral = ApproxFun.norm(hh_int(0) - hh_int(-10))
		return (p,integral)
	end

	# optinal
	function q4()

		return fig
	end

	# I found having those useful for q5
	mutable struct ChebyType
		f::Function # fuction to approximate
		nodes::Union{Vector,LinRange} # evaluation points
		basis::Matrix # basis evaluated at nodes
		coefs::Vector # estimated coefficients

		deg::Int 	# degree of chebypolynomial
		lb::Float64 # bounds
		ub::Float64

		# constructor
		function ChebyType(_nodes::Union{Vector,LinRange},_deg,_lb,_ub,_f::Function)
			n = length(_nodes)
			y = _f(_nodes)
			_basis = Float64[ChebyT(unitmap(_nodes[i],_lb,_ub),j) for i=1:n,j=0:_deg]
			_coefs = _basis \ y  # type `?\` to find out more about the backslash operator. depending the args given, it performs a different operation
			# create a ChebyComparer with those values
			new(_f,_nodes,_basis,_coefs,_deg,_lb,_ub)
		end
	end

	# function to predict points using info stored in ChebyType
	function predict(Ch::ChebyType,x_new)

		true_new = Ch.f(x_new)
		basis_new = Float64[ChebyT(unitmap(x_new[i],Ch.lb,Ch.ub),j) for i=1:length(x_new),j=0:Ch.deg]
		basis_nodes = Float64[ChebyT(unitmap(Ch.nodes[i],Ch.lb,Ch.ub),j) for i=1:length(Ch.nodes),j=0:Ch.deg]
		preds = basis_new * Ch.coefs
		preds_nodes = basis_nodes * Ch.coefs

		return Dict("x"=> x_new,"truth"=>true_new, "preds"=>preds, "preds_nodes" => preds_nodes)
	end

	function q5(deg=(5,9,15),lb=-1.0,ub=1.0)

		runge(x) = 1.0 ./ (1 .+ 25 .* x.^2)

		PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q5.png"))

	end

	function q6()

		# compare 2 knot vectors with runge's function

		PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q6.png"))

	end

	function q7(pr = 0.025)
		deg = 3
		lb = -1
		ub = 1
		nk = 13
		nev = 65
		f(x) = (abs.(x)).^0.5
		bs1 = BSpline(nk,deg,lb,ub)
		knots2 = vcat(range(lb,stop = -pr, length = 5),0,0,0,range(pr,stop = ub, length = 5))
		bs2 = BSpline(knots2,deg)
		eval = Vector(range(lb,stop = ub, length = nev))
		c1 = getBasis(eval,bs1) \ f(eval)
		c2 = getBasis(eval,bs2) \ f(eval)
		grid = Vector(range(lb, stop = ub, length = 1000))
		tr = f(grid)
		a1 = getBasis(grid, bs1) * c1
		a2 = getBasis(grid, bs2) * c2
		err1 = a1 .- tr
		err2 = a2 .- tr
		p = Plots.plot(layout = 3, dpi = 400)
		Plots.plot!(p[1],grid,tr,label = "True",lw = 1.5,linecolor = "black")
		Plots.plot!(p[2],grid,a1, label = "Approx1", lw = 1.5, linecolor = "blue")
		Plots.plot!(p[2],grid,a2, label = "Approx2", lw = 1.5, linecolor = "red")
		Plots.plot!(p[3],grid,err1, label = "error1", lw = 1.5, linecolor = "blue")
		Plots.plot!(p[3],grid,err2, label = "error2", lw = 1.5, linecolor = "red")
		Plots.savefig(p,joinpath(dirname(@__FILE__),"..","q7.png"))
		return (p,maximum(abs,err2))
	end

	#function q7bis()
	#	G = Vector(range(0.001, stop = 0.2, length = 1000))
	#	return Plots.plot(G,q7.(G))
	#end

	# function to run all questions
	function runall()
		@info("running all questions of HW-funcapprox:")
		q1(15)
		q2(3)
		q3(10)
		#q4()
		#q5()
		#q6()
		q7()
	end

end # module


text:
#=
    This file includes all the funtions needed to compute the techno-economic
    indicators
 =#

 mutable struct NPV{T <: Array{Float64}}
      capex::T
      opex::T
      salvage::T
      cf::T
      cumulative_npv::T
      total::T
 end

 # Compute costs
NPV(mg::Microgrid, designer::AbstractDesigner) = NPV(1:mg.parameters.ns, mg, designer)
# Compute costs for a given scenario s
function NPV(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, designer::AbstractDesigner)

    # Discount factor
    γ = repeat(1. ./ (1. + mg.parameters.τ) .^ range(0, length = mg.parameters.ny, step = mg.parameters.Δy), 1, length(s))

    # Discounted capex
    capexx = γ .* capex(s, mg, designer)

    # Discounted opex
    opex = γ .* dropdims(baseline_cost(s, mg) .-  grid_cost(s, mg), dims=1)

    # Discounted salvage value
    salvage = γ .*  salvage_value(s, mg)

    # Discounted cash flow
    cf = - capexx .+ opex .+ salvage

    # Discounted NPV each year
    cumulative = cumsum(cf, dims = 1)

    # Discounted NPV
    total = sum(cf, dims=1)

    return NPV(capexx, opex, salvage, cf, cumulative, total)
end

# Baseline cost
baseline_cost(mg::Microgrid) = baseline_cost(1:mg.parameters.ns, mg)
baseline_cost(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid) = baseline_cost(1:mg.parameters.ny, s, mg)
function baseline_cost(y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid)
    # TODO : compute baseline cost whatever the microgrid...
    total = 0.
    if !isempty(mg.demands)
        for (k,a) in enumerate(mg.demands)
            if a.carrier isa Electricity
                total = total .+ sum(a.carrier.power[:,y,s] .* mg.grids[1].cost_in[:,y,s] * mg.parameters.Δh, dims = 1)
            elseif a.carrier isa Heat
                total = total .+ sum(a.carrier.power[:,y,s] / mg.converters[isin(mg.converters, Heater)[2]].η_E_H .* mg.grids[1].cost_in[:,y,s] * mg.parameters.Δh, dims = 1)
            end
        end
    end
    return total
end
# Grid cost
grid_cost(mg::Microgrid) = grid_cost(1:mg.parameters.ns, mg)
grid_cost(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid) = grid_cost(1:mg.parameters.ny, s, mg)
grid_cost(y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid) = sum(sum(max.(0., a.carrier.power[:,y,s]) .* a.cost_in[:,y,s] .- min.(0., a.carrier.power[:,y,s]) .* a.cost_out[:,y,s], dims = 1) * mg.parameters.Δh for a in mg.grids)

# CAPEX
capex(mg::Microgrid, designer::AbstractDesigner) = capex(1:mg.parameters.ns, mg, designer)
# CAPEX for a given scenario s
function capex(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, designer::AbstractDesigner)
    capex = 0.
    # Generations
    for (k, a) in enumerate(mg.generations)
        capex = capex .+ designer.decisions.generations[k][:,s] .* a.cost[:,s]
    end
    # Storages
    for (k, a) in enumerate(mg.storages)
        capex = capex .+ designer.decisions.storages[k][:,s] .* a.cost[:,s]
    end
    # Converters
    for (k, a) in enumerate(mg.converters)
        capex = capex .+ designer.decisions.converters[k][:,s] .* a.cost[:,s]
    end
    return capex
end

# Salvage value
salvage_value(mg::Microgrid) = salvage_value(1:mg.parameters.ns, mg)
# Salvage value for a given scenario s
function salvage_value(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid)
    # TODO salvage as a function of SoH
    # Linear depreciation of components
    nh, ny = mg.parameters.nh, mg.parameters.ny
    salvage = zeros(mg.parameters.ny, length(s))
    salvage[ny,:] .= 0.
    # Generations
    for a in mg.generations
        salvage[ny,:] = salvage[ny,:] .+ (a.lifetime .- ny) ./ a.lifetime .* a.cost[ny, s]
    end
    # Storages
    for a in mg.storages
        salvage[ny,:] = salvage[ny,:] .+ (a.lifetime .- ny) ./ a.lifetime .* a.cost[ny, s]
    end
    # Converters
    for a in mg.converters
        salvage[ny,:] = salvage[ny,:] .+ (a.lifetime .- ny) ./ a.lifetime .* a.cost[ny, s]
    end
    return salvage
end

mutable struct EAC{T <: Array{Float64}}
     capex::T
     opex::T
     total::T
end

function EAC(y::Int64, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, designer::AbstractDesigner)
    # Annualised capex
    capex = Genesys.annualised_capex(y:y, s, mg, designer)
    # opex
    opex = grid_cost(y+1, s, mg)

    return EAC(capex, opex, capex .+ opex)
end
# Annualised CAPEX
function annualised_capex(y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, designer::AbstractDesigner)
    # Preallocation
    capex = 0.
    # Generations
    for (k, a) in enumerate(mg.generations)
        Γ = (mg.parameters.τ * (mg.parameters.τ + 1.) ^ a.lifetime) / ((mg.parameters.τ + 1.) ^ a.lifetime - 1.)
        capex = capex .+ Γ .* designer.decisions.generations[k][y,s] .* a.cost[y,s]
    end
    # Storages
    for (k, a) in enumerate(mg.storages)
        Γ = (mg.parameters.τ * (mg.parameters.τ + 1.) ^ a.lifetime) / ((mg.parameters.τ + 1.) ^ a.lifetime - 1.)
        capex = capex .+ Γ .* designer.decisions.storages[k][y,s] .* a.cost[y,s]
    end
    # Converters
    for (k, a) in enumerate(mg.converters)
        Γ = (mg.parameters.τ * (mg.parameters.τ + 1.) ^ a.lifetime) / ((mg.parameters.τ + 1.) ^ a.lifetime - 1.)
        capex = capex .+ Γ .* designer.decisions.converters[k][y,s] .* a.cost[y,s]
    end
    return capex
end

# Share of renewables
renewable_share(mg::Microgrid) = renewable_share(1:mg.parameters.ns, mg)
# Share of renewables for a given scenario s
renewable_share(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid) = renewable_share(1:mg.parameters.ny, s, mg)
# Share of renewables for a given year y of a givn scenario s
function renewable_share(y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid)
    # TODO to be changed if there is no grid...
    total = 0.
    for (k,a) in enumerate(mg.demands)
        if a.carrier isa Electricity
            total = total .+ sum(a.carrier.power[:,y,s], dims = 1)
        elseif a.carrier isa Heat
            total = total .+ sum(a.carrier.power[:,y,s], dims = 1) ./ mg.converters[isin(mg.converters, Heater)[2]].η_E_H
        end
    end
    for (k,a) in enumerate(mg.grids)
        if a.carrier isa Electricity
            return share = dropdims(1. .- sum(max.(0., a.carrier.power[:,y,s]), dims = 1) ./ total, dims=1)
        else
            println("Renewable share not yet defined!")
            return nothing
        end
    end
end

# LPSP
mutable struct LPSP{T}
    elec::Union{Nothing, T}
    heat::Union{Nothing, T}
    hydrogen::Union{Nothing, T}
end

LPSP(mg::Microgrid) = LPSP(1:mg.parameters.ns, mg)
# LPSP for a given scenario s
LPSP(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid) = LPSP(1:mg.parameters.ny, s, mg)
# LPSP for a given scenario s and year y
function LPSP(y::Union{Int64, UnitRange{Int64}}, s::Union{Int64, UnitRange{Int64}}, mg::Microgrid)
    # Initialization
    elec, heat, hydrogen = nothing, nothing, nothing
    # Computation
    for a in mg.demands
        if a.carrier isa Electricity
            elec = sum(max.(0., power_balance(1:mg.parameters.nh, y, s, mg, Electricity)), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
            for aa in mg.grids
                if aa.carrier isa Electricity
                    elec = elec .- sum(max.(0., aa.carrier.power[:, y, s]), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
                end
            end
            elec = dropdims(elec,dims=1)
        elseif a.carrier isa Heat
            heat = sum(max.(0., power_balance(1:mg.parameters.nh, y, s, mg, Heat)), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
            for aa in mg.grids
                if aa.carrier isa Heat
                    heat = heat .- sum(max.(0., aa.carrier.power[:, y, s]), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
                end
            end
            heat = dropdims(heat,dims=1)
        elseif a.carrier  isa Hydrogen
            hydrogen = sum(max.(0., power_balance(1:mg.parameters.nh, y, s, mg, Hydrogen)), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
            for aa in mg.grids
                if aa.carrier isa Electricity
                    hydrogen = hydrogen .- sum(max.(0., aa.carrier.power[:, y, s]), dims=1) ./ sum(a.carrier.power[:, y, s], dims = 1)
                end
            end
            hydrogen = dropdims(hydrogen,dims=1)
        end
    end
    return LPSP(elec, heat, hydrogen)
end

mutable struct Metrics{T}
    baseline::T
    npv::NPV{T}
    eac::EAC{T}
    renewable_share::T
    lpsp::LPSP{T}
end

# Compute indicators
Metrics(mg::Microgrid, designer::AbstractDesigner) = Metrics(1:mg.parameters.ns, mg, designer)
# Compute indicators for a given scenario s
function Metrics(s::Union{Int64, UnitRange{Int64}}, mg::Microgrid, designer::AbstractDesigner)
    # Baseline cost
    baseline = dropdims(baseline_cost(mg), dims = 1)
    # NPV
    npv = NPV(s, mg, designer)
    # EAC
    eac = EAC(1, s, mg, designer)
    # Share of renewables
    share = renewable_share(s, mg)
    # LPSP
    lpsp = LPSP(s, mg)

    return Metrics(baseline, npv, eac, share, lpsp)
end


text:
function get_a(obj::ButcherArrayProvider)
    return jcall(obj, "getA", Vector{Vector{jdouble}}, ())
end

function get_b(obj::ButcherArrayProvider)
    return jcall(obj, "getB", Vector{jdouble}, ())
end

function get_c(obj::ButcherArrayProvider)
    return jcall(obj, "getC", Vector{jdouble}, ())
end


text:
developing = false; # make it false for CI testing
if developing
	# ] add Revise # if you don't have it
	path = pwd()
	if path[end-6:end] != "MatLang"
		error("cd(to the MatLang path)")
	end
	push!(LOAD_PATH, ".")
	using Revise
	using MatLang
	using Test

else
	using MatLang
	using Test
end
################################################################
# These tests don't test every single method. For more examples, refer to examples folder.

@testset "MatLang.j - Entering_Commands" begin

end
@testset "MatLang.jl - Matrices_and_Arrays" begin
	@test zerosM(:mat, 2) == zeros(2, 2)

	@test zerosM(:mat, 2, like = ones(Int32, 2, 2)) == zeros(Int32, 2, 2)

	@test zerosM([2, 3]) == zerosM((2, 3))

	@test zerosM(Int8, [2, 3]) == zerosM(Int8, (2, 3))

	@test onesM(:mat, 2) == ones(2, 2)

	@test onesM([2, 3]) == onesM((2, 3))

	@test onesM(Int8, [2, 3]) == onesM(Int8, (2, 3))

	@test falseM([2, 2]) == 0*trueM(:mat,2)

	@test eyeM(2, 3) == [1 0 0; 0 1 0]

	@test eyeM(Int32, 2, 3) == [1 0 0; 0 1 0]

	@test eyeM(2, 3, like = ones(Int32, 2, 2)) == [1 0 0; 0 1 0]

	@test diagM([2, 3, 4], 1) == [0 2 0 0; 0 0 3 0; 0 0 0 4; 0 0 0 0]

	@test diagM([1 2 3; 4 5 6; 7 8 9]) == [1; 5; 9]

	@test blkdiagM(ones(2, 2), 2 * ones(2, 2)) == [ones(2, 2) zeros(2, 2); zeros(2, 2)  2 * ones(2, 2)]

	@test blkdiagM(ones(2, 4), 2 * ones(3, 2)) == [ones(2, 4) zeros(2, 2); zeros(3, 4)  2 * ones(3, 2)]

	@test catM(1, ones(3, 3), zeros(3, 3), 2 * ones(3, 3)) == [ones(3, 3); zeros(3, 3); 2 * ones(3, 3)]

	@test horzcatM(ones(3, 3), zeros(3, 3)) == [ones(3, 3) zeros(3, 3)]

	@test repelemM([1 2 3 4], 3) == [1; 1; 1; 2; 2; 2; 3; 3; 3; 4; 4; 4]

	@test repmatM([1; 2], 3, 2) == [[1; 2] [1; 2]; [1; 2] [1; 2]; [1; 2] [1; 2]]

	@test linspaceM(1, 10, 2) == 1:2:10

	@test logspaceM(1, 5, :equal) == 10.0.^(1:5)

	@test meshgridM(1:2:5, 1:2:5)==([1 3 5;1 3 5; 1 3 5],[1 1 1; 3 3 3; 5 5 5])

	@test freqspaceM(10, dim = 1) == 0:0.2:1

	@test ndgridM(1:2:5, 1:2:5)==([1 1 1; 3 3 3; 5 5 5], [1 3 5;1 3 5; 1 3 5])

	@test lengthM([2 3 4 5 6 7; 1 2 3 4 5 6]) == 6

	@test sizeM(:arr, [2 3 4 5 6 7; 1 2 3 4 5 6]) == [2; 6]

	@test sizeM(rand(3, 5, 4), 2, 3) == (5, 4)

	@test numelM([1, 2, 3, 4]) == 4

	@test [
	    isscalarM(:mat, 1)
	    isscalarM(:mat, 5 * ones(1, 1, 1))
	    isscalarM(:mat, "H")
	    isscalarM(:mat, "Hi")
	    isscalarM(:mat, ["Hi"])
	    isscalarM(:mat, ["Hi" "Bye"])
	    isscalarM(:mat, 'H')
	] == [true; true; true; false; true; false; true]

	@test [
	    isvectorM(:mat, [1; 2; 3])
	    isvectorM(:mat, [1 2 3])
	    isvectorM(:mat, ones(3, 1))
	    isvectorM(:mat, ones(1, 3))
	    isvectorM(:mat, ones(3))
	    isvectorM(:mat, 1)
	    isvectorM(:mat, "Hi")
	    isvectorM(:mat, ["Hi", "Bye"])
	] == [true; true; true; true; true; true; true; true]

	@test [
	ismatrixM(:mat, [1 2 3; 3 5 6])
	ismatrixM(:mat, [1;2;3])
	ismatrixM(:mat, ones(3, 1))
	ismatrixM(:mat, ones(1, 3))
	ismatrixM(:mat, ones(3))
	ismatrixM(:mat, "Hi")
	ismatrixM(:mat, ["Hi", "Bye"])
	] == [true; true; true; true; true; true; true]

	@test [
	isrowM([1 2 3]) # true
	isrowM(ones(3, 1)) # false
	isrowM(ones(1, 3)) # true
	isrowM(ones(3)) # false
	isrowM(1) # false
	isrowM("Hi") # false
	isrowM(["Hi", "Bye"]) # false
	isrowM(["Hi" "Bye"]) # true
	] == [true; false; true; false; false; false; false; true]

	@test [

	iscolumnM([1; 2; 3]) # false
	iscolumnM([1 2 3]) # false
	iscolumnM(ones(3, 1)) # tue
	iscolumnM(ones(1, 3)) #false
	iscolumnM(ones(3)) # false
	iscolumnM(1) # false
	iscolumnM("Hi") # false
	iscolumnM(["Hi", "Bye"]) # false
	iscolumnM(["Hi" "Bye"]) # true

	] == [false; false; true; false; false; false; false; false; false]
	@test sortM([1 5 3; 4 1 10], 2, :ascend) == [1 3 5; 1 4 10]

	B = [1 5 3; 4 1 10]
	sorted, isort = sortM(B, 2, :ascend, I = true)
	@test 	B[isort] == sorted

	@test flipM([1 2; 3 4]) == [3 4; 1 2]

	@test  rot90M([1 2; 3 4], 3) == [3 1; 4 2]

	@test transposeM([2 3 4 5 6 7; 1 2 3 4 5 6]) == transposeM(:arr, [2 3 4 5 6 7; 1 2 3 4 5 6])

	@test squeezeM(ones(2,1,2)) == [1 1; 1 1]
end

@testset "MatLang.j - Numeric_Types" begin

	@test typeof(doubleM(10.0)) == Float64

	@test typeof(singleM(10.0)) == Float32

	@test typeof(int8M(-10)) == Int8

	@test typeof(int16M(-10)) == Int16

	@test typeof(int32M(-10)) == Int32

	@test typeof(int64M(-10)) == Int64

	@test typeof(uint8M(10)) == UInt8

	@test typeof(uint16M(10)) == UInt16

	@test typeof(uint32M(10)) == UInt32

	@test typeof(uint64M(10)) == UInt64

end


text:
using Printf

include("INC.jl")

function create_lattice_vars()
    LatVecs = zeros(3,3)
    A = 5.13
    LatVecs[1,:] = [A, A, 0.0]
    LatVecs[2,:] = [A, 0.0, A]
    LatVecs[3,:] = [0.0, A, A]
    lattice_vars = LatticeVars( LatVecs )
    return lattice_vars
end

function create_atomic_vars(lattice_vars)

    maxatoms = 200
    maxspecies = 8
    atposl = zeros(3,maxatoms,maxspecies)

    nspecies = 2
    natoms = [1,1]

    # species 1, atom 1
    atposl[:,1,1] = [0.0, 0.0, 0.0]
    # species 2, atom 1
    atposl[:,1,2] = [0.25, 0.25, 0.25]

    atomic_vars = AtomicVars(nspecies, natoms, atposl, lattice_vars)
end

function main()

    latt_vars = create_lattice_vars()
    atm_vars = create_atomic_vars(latt_vars)

    Nspecies = 2
    atsp_vars = AtomicSpeciesVars(Nspecies)
    mtr_vars = MuffinTinRadialVars(Nspecies)
    apwlo_vars = APWLOVars(Nspecies, mtr_vars.lmaxapw)

    readspecies!(1, "DATA_species/Si.in", atsp_vars, mtr_vars, apwlo_vars)
    readspecies!(2, "DATA_species/Pt.in", atsp_vars, mtr_vars, apwlo_vars)

    checkmt!( latt_vars, atm_vars, atsp_vars.spsymb, mtr_vars )

    println("Pass here")

end

main()


text:
using MPIFiles, PyPlot

filename = "../../../test/measurement_V2.mdf"
f = MPIFile(filename)

u = getMeasurementsFD(f, frames=1:100, numAverages=100)
figure(6, figsize=(6,4))
semilogy(abs.(u[1:400,1,1,1]))
savefig("../assets/spectrum1.png")


text:
# This file was generated, do not modify it. # hide
r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
using Documenter
using TensorKit

makedocs(modules=[TensorKit],
            sitename="TensorKit.jl",
            authors = "Jutho Haegeman",
            format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true",
                                        mathengine = MathJax()),
            pages = [
                "Home" => "index.md",
                "Manual" => ["man/intro.md", "man/tutorial.md", "man/categories.md",
                                "man/spaces.md", "man/sectors.md", "man/tensors.md"],
                "Library" => ["lib/sectors.md","lib/spaces.md","lib/tensors.md"],
                "Index" => ["index/index.md"]
            ])

deploydocs(repo = "github.com/Jutho/TensorKit.jl.git")


text:
include("../src/Cinta.jl")

using .Cinta

my_parameters = Parameters(10.5, .2)
my_line_retriever = LineRetrieverFromFile("src/datos_cinta.txt")

# println("Ruta\tInc\tMin\tVel ofrecida\t\tBotón\tSubir o bajar\tVel")
println(get_output(my_parameters, my_line_retriever))

max_stars_repo_name: AndresParraSilva/julia-start

max_issues_repo_name: AndresParraSilva/julia-start

max_forks_repo_name: AndresParraSilva/julia-start


text:
module RoesserModels

VERSION < v"0.4-" && using Docile

@doc doc"""
The Roesser model is a 2D state-space model, useful for 2D spatial dynamics
or a combination of 1D spatial dynamics and time.

References:
  \cite{Roesser:1975}
"""->
type RoesserModel
  A1::Array{Float64}
  A2::Array{Float64}
  A3::Array{Float64}
  A4::Array{Float64}
  B1::Array{Float64}
  B2::Array{Float64}
  C1::Array{Float64}
  C2::Array{Float64}
  D::Array{Float64}
  in_dim::Int
  out_dim::Int
  x_h_dim::Int
  x_v_dim::Int
  function RoesserModel(A1, A2, A3, A4, B1, B2, C1, C2, D)
    x_h_dim = size(A1, 1)
    x_v_dim = size(A4, 1)
    out_dim = size(C1, 1)
    in_dim = size(B1, 2)

    if ((size(A1, 2) != x_h_dim) || (size(A3,2) != x_h_dim)
        || (size(B1,1) != x_h_dim) || (size(C1,2) != x_h_dim)
        || (size(A2,1) != x_h_dim))
      error("Matrix dimensions for Horizontal states do not agree");
    end
    if  ((size(A2,2) != x_v_dim) || (size(A4,1) != x_v_dim)
         || (size(A4,2) != x_v_dim) || (size(B2,1) != x_v_dim)
         || (size(C2,2) != x_v_dim))
      error("Matrix dimensions for Vertical states do not agree");
    end
    if ((size(C2,1) != out_dim) || (size(D,1) != out_dim)
        || (size(B2,2) != in_dim) || (size(D,2) != in_dim))
      error("Matrix dimensions for Input/Output do not agree")
    end
    new(A1, A2, A3, A4, B1, B2, C1, C2, D, in_dim, out_dim, x_h_dim, x_v_dim)
  end
end

@doc doc"""
The CRSD (Causal Recursive Separable in Denominator) model is a specific type
of Roesser Models.

References:
  \cite{Roesser:1975}
"""->
function CrsdRoesserModel(A1, A2, A4, B1, B2, C1, C2, D)
  x_h_dim = size(A1, 1)
  x_v_dim = size(A4, 1)
  out_dim = size(C1, 1)
  in_dim = size(B1, 2)
  A3 = zeros(x_v_dim, x_h_dim);
  return RoesserModel(A1, A2, A3, A4, B1, B2, C1, C2, D)
end

@doc doc"""
Simulates the Roesser model in an iterative approach.

"""->
function simulate(system::RoesserModel, inputs)
  if system.in_dim > 1 && size(inputs,3) != system.in_dim
    error("Input dimension mismatch")
  end

  M = size(inputs,1)
  N = size(inputs,2)
  output = zeros(M, N, system.out_dim)
  x_h = zeros(size(inputs,1)+1, size(inputs,2)+1, system.x_h_dim)
  x_v = zeros(size(inputs,1)+1, size(inputs,2)+1, system.x_v_dim)

  A1 = system.A1
  A2 = system.A2
  A3 = system.A3
  A4 = system.A4
  B1 = system.B1
  B2 = system.B2
  C1 = system.C1
  C2 = system.C2
  D = system.D
  for i = 1:M
    for j = 1:N
      input_ij = vec(inputs[i,j,1:end])';
      x_h[i+1,j, 1:end] = A1*vec(x_h[i,j, 1:end]) + A2*vec(x_v[i,j,1:end]) + B1 * input_ij
      x_v[i,j+1, 1:end] = A3*vec(x_h[i,j, 1:end]) + A4*vec(x_v[i,j,1:end]) + B2 * input_ij
      output[i,j, 1:end] = C1 * vec(x_h[i,j, 1:end]) + C2 * vec(x_v[i,j,1:end]) + D * input_ij
    end
  end
  return output
end

@doc doc"""
Generate a completely random system. No check is performed for observability
or controllability.

"""->
function randomsystem(nh::Integer, nv::Integer, outputs::Integer, inputs::Integer)
  A1 = randn(nh, nh);
  A2 = randn(nh, nv);
  A3 = randn(nv, nh);
  A4 = randn(nv, nv);
  C1 = randn(outputs, nh);
  C2 = randn(outputs, nv);
  B1 = randn(nh, inputs);
  B2 = randn(nv, inputs);
  D = randn(outputs, inputs);

  return RoesserModel(A1, A2, A3, A4, B1, B2, C1, C2, D);
end

@doc doc"""
Tests the model for global controllability, according to the criteria
formulated by Kung in 1976.

References:
  \cite{Kung:1976}
"""->
function is_global_controllable(system::RoesserModel, M::Integer, N::Integer)

end

@doc doc"""
Tests the model for global controllability, according to the criteria
formulated by Kung in 1976.

References:
  \cite{Kung:1976}
"""->
function is_global_observable(system::RoesserModel, M::Integer, N::Integer)

end

@doc doc"""
Tests the model for local controllability, according to the criteria
formulated by Roesser in 1975.

References:
  \cite{Roesser:1975}
  \cite{Kung:1976}
"""->
function is_local_controllable(system::RoesserModel)

end

@doc doc"""
Tests the model for local observability, according to the criteria
formulated by Roesser in 1975.

References:
  \cite{Roesser:1975}
  \cite{Kung:1976}
"""->
function is_local_observable(system::RoesserModel)

end

@doc doc"""
Tests the model for modal controllability, according to the criteria
formulated by Kung in 1976.

References:
  \cite{Kung:1976}
"""->
function is_modal_controllable(system::RoesserModel)

end

@doc doc"""
Tests the model for modal observability, according to the criteria
formulated by Kung in 1976.

References:
  \cite{Kung:1976}
"""->
function is_modal_observable(system::RoesserModel)

end

end # module


text:
"""
    function plot_PCR_point_est_fit(θ_fitted::NamedTuple,model::CoVModel,k::Integer)

Plot the fit of a point estimate with and without vaccination. Incidence is not adjusted for number of test and probability of test positive among cases

"""
function plot_PCR_point_prediction(params,init,model,model_novacc,PCR_array,θ_fitted,all_IHRIFR_fits,p_ID,p_IH,scenario,startdate,name)

    #Test per day to scale prediction to match numbers of tests
    tests_per_day = sum(model.case_data[:,:,2],dims = 2)[:] #sum(model.case_data[:,:,:],dims = [2,3])[:]
    mean_tests_per_day_after_sept = mean(tests_per_day[300:end]) #+2000
    tests_per_day = tests_per_day[(startdate-Date(2020,2,20)).value:end]  # for when data is available after start date
    tests_per_day_new =  vcat(tests_per_day,ones(366-length(tests_per_day)).*mean_tests_per_day_after_sept)
    #tests_per_day_new = ones(366).* mean(tests_per_day)  # use average number of tests in the past for the year 2021

    immune_escape_variant_time = Float64((Date(2021,9,1) - Date(2021,1,31)).value) #  September 1st, time to implement variant immunity escape
    function immune_escape_variant!(integrator)
        # Reduce vaccine efficacies by 50%
        integrator.p[45:59] *= 0.5

        # Start with 'variant_E' exposed number of individuals
        variant_E = 50
        prop_vacc = sum(integrator.u[:,1:7,:],dims=2)[:,1,:] ./ sum(integrator.u[:,1:7,:],dims=[2,3])[:]
        integrator.u[:,2,:] .= variant_E .* prop_vacc # number exposed in each age group and vaccinatin status
        integrator.u[:,1,:] .= integrator.u[:,1,:] .- integrator.u[:,2,:] #substract E from S

        # move 50% from R -> S
        loss_susceptible = 0.5 .* integrator.u[:,7,:]
        integrator.u[:,1,:] .+= loss_susceptible
        integrator.u[:,7,:] .-= loss_susceptible
    end
    immune_escape_variant_cb = PresetTimeCallback([immune_escape_variant_time],immune_escape_variant!,save_positions=(false,false))

   # Update callbacks based on scenario
    if  scenario>2  # For Scenario 3 and 4
     model.cb = CallbackSet(model.cb, immune_escape_variant_cb)
     model_novacc.cb = CallbackSet(model_novacc.cb, immune_escape_variant_cb)
    end

    #Solve dynamical system
    p = copy(params)
    sol= solve(model.prob, Vern9(); callback = model.cb, u0=init, p=p, saveat = 1)
    p = copy(params)
    sol_novacc= solve(model_novacc.prob, Vern9(); callback = model_novacc.cb, u0=init, p=p, saveat = 1)

    incidence = KenyaCoVaccines.get_incidence_time_array(sol)
    PCR_pred = VectorOfArray([KenyaCoVaccines.simple_conv(incidence[:,j],PCR_array) for j = 1:size(incidence,2)])
    total_PCR_pred = sum(PCR_pred[:,:],dims = 2)[:]
    beta_bin_pred = θ_fitted.χ_PCR_array.*total_PCR_pred./((θ_fitted.χ_PCR_array-1).*total_PCR_pred .+ sum(model.N))
    pred = beta_bin_pred.*tests_per_day_new[1:size(PCR_pred,1)]

    incidence = KenyaCoVaccines.get_incidence_time_array(sol_novacc)
    PCR_pred = VectorOfArray([KenyaCoVaccines.simple_conv(incidence[:,j],PCR_array) for j = 1:size(incidence,2)])
    total_PCR_pred= sum(PCR_pred[:,:],dims = 2)[:]
    beta_bin_pred = θ_fitted.χ_PCR_array.*total_PCR_pred./((θ_fitted.χ_PCR_array-1).*total_PCR_pred .+ sum(model.N))
    pred_novacc = beta_bin_pred.*tests_per_day_new[1:size(PCR_pred,1)]

    cum_inc_vac_children = VectorOfArray([sum(u[1:4,8,:]) for u in sol.u])
    cum_inc_novac_children = VectorOfArray([sum(u[1:4,8,:]) for u in sol_novacc.u])
    cum_inc_vac_adults = VectorOfArray([sum(u[5:14,8,:]) for u in sol.u])
    cum_inc_novac_adults = VectorOfArray([sum(u[5:14,8,:]) for u in sol_novacc.u])
    cum_inc_vac_elderly = VectorOfArray([sum(u[15:end,8,:]) for u in sol.u])
    cum_inc_novac_elderly = VectorOfArray([sum(u[15:end,8,:]) for u in sol_novacc.u])

    #xticks_monthly = [(Date(y,m,1) - Date(2020,2,20)).value for m = 1:12,y = 2020:2021][1:end]
    #xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2020:2021][1:end]
    xticks_monthly = [(Date(y,m,1) - Date(2021,1,1)).value for m = 1:12,y = 2021][1:end]
    xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2021][1:end]

    pcr_pos = sum(model.case_data[:,:,1],dims = 2)[(Date(2021,1,31)-Date(2020,2,20)).value : end]
    pcr_pos_ma = KenyaCoVaccines.sma(pcr_pos,7)  # add 7-day movinf average to the plot
    plt = plot(pred,
            lab = "With vaccination",
            lw=1,
            legend = :topright,
            title = "PCR prevalence $(name)",
            xticks = (xticks_monthly,xticklabs_monthly),
            ylabel = "Daily cases",
            color=:blue,
            size = (800,600),dpi = 300);
    plot!(plt,pred_novacc,color=:red,lab = "Without vaccination",lw=1);
    plot!(plt,pcr_pos_ma,lw=1,color=:green,lab="7-day MA PCR prevalence");
    scatter!(plt,pcr_pos,color=:green,lab="Observed PCR prevalence")

    plt2=plot(cum_inc_vac_children,lab = "vaccines: 0-19 yrs old",legend = :bottomright,title = "Cumulative incidence $(name)",xlabel = "days",lw=3);
    plot!(plt2,cum_inc_novac_children,lab = "no vaccines: 0-19 yrs old",lw=3)

    plt3=plot(cum_inc_vac_adults,lab = "vaccines: 20-64 yrs old",legend = :bottomright,title = "Cumulative incidence $(name)",xlabel = "days",lw=3);
    plot!(plt3,cum_inc_novac_adults,lab = "no vaccines: 20-64 yrs old",lw=3)

    plt4=plot(cum_inc_vac_elderly,lab = "vaccines: 65+ yrs old",legend = :bottomright,title = "Cumulative incidence $(name)",xlabel = "days",lw=3);
    plot!(plt4,cum_inc_novac_elderly,lab = "no vaccines: 65+ yrs old",lw=3)

    p = copy(params)
    vacc_efficacy_severe_disease_risk = p[55:59]

    IHR = all_IHRIFR_fits[[nm[1] .== name for nm in all_IHRIFR_fits]][1].IHR

    υ = (1 .- vacc_efficacy_severe_disease_risk') .* IHR # efficacy multiplied by risk of severe disease by age
    υ_nv = ones(5)' .* IHR # No vaccine efficacy

    incidence = diff([u[:,8,:] for u in sol.u])
    dis_inc =  [inc .* υ for inc in incidence]
    disease_inc =  [sum(inc) for inc in dis_inc]  #sum over age and dose
    disease_inc =  KenyaCoVaccines.simple_conv(disease_inc,p_IH)  # convolve with time to hsop function to get right timing of Hospitalisation

    incidence_nv = diff([u[:,8,:] for u in sol_novacc.u])
    dis_inc_nv =  [inc .* υ_nv for inc in incidence_nv]
    disease_inc_nv = [sum(inc) for inc in dis_inc_nv]  #sum over age and dose
    disease_inc_nv =  KenyaCoVaccines.simple_conv(disease_inc_nv,p_IH)  # convolve with time to hsop function to get right timing of Hospitalisation

    l = min(length(disease_inc_nv),length(disease_inc))
    severe_disease_averted = cumsum(disease_inc_nv[1:l] .- disease_inc[1:l])  # cummulative

    cum_incidence =  VectorOfArray([sum(u[:,8,:]) for u in sol.u])
    cum_incidence_nv =  VectorOfArray([sum(u[:,8,:]) for u in sol_novacc.u])
    cases_averted  = cum_incidence_nv[1:(l-1)] .- cum_incidence[1:(l-1)] # cummulative

    # deaths averted
    IFR = all_IHRIFR_fits[[nm[1] .== name for nm in all_IHRIFR_fits]][1].IFR
    υ = (1 .- vacc_efficacy_severe_disease_risk') .* IFR # efficacy multiplied by risk of death by age, assumes only those with severe disease die
    υ_nv = ones(5)' .* IFR # No vaccine efficacy

    incidence = diff([u[:,8,:] for u in sol.u])
    d_inc =  [inc .* υ for inc in incidence]
    death_inc =  [sum(inc) for inc in d_inc]  #sum over age and dose
    deaths =  KenyaCoVaccines.simple_conv(death_inc,p_ID)  # convolve with time to death function to get right timing of deaths
    cumm_deaths = cumsum(deaths)

    incidence = diff([u[:,8,:] for u in sol_novacc.u])
    d_inc =  [inc .* υ_nv for inc in incidence]
    death_inc =  [sum(inc) for inc in d_inc]  #sum over age and dose
    deaths =  KenyaCoVaccines.simple_conv(death_inc,p_ID)  # convolve with time to death function to get right timing of deaths
    cumm_deaths_novacc = cumsum(deaths)

    deaths_averted = cumm_deaths_novacc .- cumm_deaths

    plot(severe_disease_averted,dpi = 300,lab="Severe disease", title = "Cummulative Cases/Disease/Deaths averted $(name)",
     lw=3,color=:red,legend =:topleft,ylabel= "Disease/Deaths",xticks=nothing)
    plot!(cases_averted,lab="Cases", ylabel= "Cases",xticks=nothing,legend =:topleft, lw=3,color=:blue)
    plt5=plot!(twinx(),deaths_averted,lab="Deaths",size=(800,600),xticks=(xticks_monthly,xticklabs_monthly), lw=3, color=:black,legend =:topright)

    @save("scenario_$(scenario)/predicted_states/solutions_$(name).jld2",sol)
    @save("scenario_$(scenario)/deaths_averted/deaths_averted_$(name).jld2",deaths_averted)
    @save("scenario_$(scenario)/cases_averted/cases_averted_$(name).jld2", cases_averted)
    @save("scenario_$(scenario)/severe_disease_averted/severe_disease_averted_$(name).jld2",severe_disease_averted)

    @save("scenario_$(scenario)/deaths/deaths_$(name).jld2",cumm_deaths)
    @save("scenario_$(scenario)/cases/cases_$(name).jld2", cum_incidence)
    @save("scenario_$(scenario)/severe_disease/severe_disease_$(name).jld2",disease_inc)

    @save("scenario_$(scenario)/deaths/deaths_$(name).jld2",cumm_deaths_novacc)
    @save("scenario_$(scenario)/cases/cases_$(name).jld2", cum_incidence_nv)
    @save("scenario_$(scenario)/severe_disease/severe_disease_$(name).jld2",disease_inc_nv)

    return plt, plt2, plt3, plt4, plt5
end

"""
Function to plot vaccine uptake for a given county
    Plots vaccine doses over time (daily) and the cumulative number of doses at a given time by age group
"""

function plot_vaccination_uptake(time,dose1_data,dose2_data)
    xticks_monthly = [(Date(y,m,1) - Date(2021,1,31)).value for m = 1:12,y = 2021][1:end]
    xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2021][1:end]

    vacc_data = dose1_data .+ dose2_data

    plt1= plot(sum(vacc_data, dims=2), lw=3, lab="Dose 1 + Dose 2",xticks = (xticks_monthly,xticklabs_monthly),legend = :topleft, ylims=[0,maximum(sum(vacc_data, dims=2)[:])+500])
    plot!(plt1,sum(dose1_data,dims=2),ylabel="Number of doses", xlabel = "Time",lab="Dose 1", lw=3)

    cum_doses = zeros(size(dose1_data))
    for  age in 1:17
     cum_doses[:,age] = cumsum(vacc_data[:,age])
    end

    xticks_age = 1:17
    xticklabs_age= vcat([string(a*5)*"-"*string(a*5 + 4) for a = 0:15],"80+")

    plt2=bar(1:17,cum_doses[time,:], xticks = (xticks_age,xticklabs_age),
        legend= :none,
        xlabel = "Age group",
        title = "Number of doses administered as at $(Date(2021,1,31) + Dates.Day(time))",
        ylabel = "Cummulative number of doses",
        size = (800,600))

    return plt1, plt2
end

"""
    function gather_monthly_sero_data(sero_data,years,date_first_day::Date)
Calculate the all age group sero-positivity with Jeffery 95% CIs for the point estimate. X tick positions are given
relative to the start date of the data set: `date_first_day`.
"""

function gather_monthly_sero_data(sero_data,years,date_first_day::Date)
        data_dates = [date_first_day + Day(k-1) for k=1:size(sero_data,1)]
        monthly_x_pos = Float64[]
        monthly_seropos = Float64[]
        lb_seropos = Float64[]
        ub_seropos = Float64[]

        for yr in years
                for mnth in 1:12
                        daysinmonthyear = findall((month.(data_dates).== mnth).&(year.(data_dates).== yr) )
                        pos_sero = sum(sero_data[daysinmonthyear,:,1])
                        neg_sero = sum(sero_data[daysinmonthyear,:,2])
                        if pos_sero + neg_sero > 0
                                p_sero = pos_sero/(pos_sero+neg_sero)
                                push!(monthly_x_pos,mean(daysinmonthyear))
                                push!(monthly_seropos,p_sero)
                                d = Beta(pos_sero +0.5,neg_sero +0.5)
                                push!(lb_seropos,p_sero-invlogcdf(d,log(0.025)))
                                push!(ub_seropos,invlogcdf(d,log(0.975)) - p_sero )
                        end
                end
        end
        return (monthly_x_pos=monthly_x_pos,monthly_seropos=monthly_seropos,lb_seropos=lb_seropos,ub_seropos=ub_seropos)
end

"""
    function plot_PCR_point_est_fit(θ_fitted::NamedTuple,model::CoVModel,k::Integer)
Plot the fit of a point estimate `θ_fitted` of the underlying parameters to the PCR data for the `k`th
area of the `CoVModel` type `model`.
"""
function plot_PCR_point_est_fit(θ_fitted::NamedTuple,model,k;startdate = Date(2020,2,20),enddate = Date(2021,3,1))
        #  Get parameters
        u0 = KenyaCoVaccines.create_initial_conditions(θ_fitted.E₀_array[k],model.subareas[k].N)
        p = deepcopy(model.subareas[k].p)
        p[1] = θ_fitted.β₀
        p[41] = θ_fitted.β_home
        p[42] = θ_fitted.β_school
        p[43] = θ_fitted.β_other
        p[44] = θ_fitted.β_work

        #Solve dynamical system
        l_solve = (enddate - startdate).value
        sol = solve(model.subareas[k].prob, BS3();callback = model.subareas[k].cb, u0=u0, p=p, saveat = 1,tspan = (0,l_solve))
        incidence = KenyaCoVaccines.get_incidence_time_array(sol)
        PCR_pred = VectorOfArray([KenyaCoVaccines.simple_conv(incidence[:,j],model.PCR_array) for j = 1:size(incidence,2)])
        total_PCR_pred = sum(PCR_pred[:,:],dims = 2)[:]

        #Find days with and without neg. PCR test result data and scale prediction to match numbers of tests
        total_days_with_neg_test = trues(length(total_PCR_pred))
        days_without_neg_test = findall(model.subareas[k].case_data[:,1,2] .== -1)
        days_with_neg_test = findall(model.subareas[k].case_data[:,1,2] .>= 0)
        total_days_with_neg_test[days_without_neg_test] .= false
        tests_per_day_from_data = sum(model.subareas[k].case_data[:,:,2],dims = 2)[:]
        mean_tests_per_day_after_sept = mean(tests_per_day_from_data[194:end])

        tests_per_day_for_pred = fill(mean_tests_per_day_after_sept,length(total_PCR_pred))
        tests_per_day_for_pred[days_with_neg_test] .= tests_per_day_from_data[days_with_neg_test]

        neg_bin_pred = total_PCR_pred.*θ_fitted.p_test_array[k].*model.relative_testing_rate[1:size(PCR_pred,1)]
        beta_bin_pred = θ_fitted.χ_PCR_array[k].*total_PCR_pred./((θ_fitted.χ_PCR_array[k]-1).*total_PCR_pred .+ sum(model.subareas[k].N))
        beta_bin_pred = beta_bin_pred.*tests_per_day_for_pred
        pred = neg_bin_pred.*(.~total_days_with_neg_test) .+ beta_bin_pred.*total_days_with_neg_test

        #Plot

        xticks_monthly = [(Date(y,m,1) - Date(2020,2,20)).value for m = 1:12,y = 2020:2021][3:end]
        xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2020:2021][3:end]
        plt = scatter(sum(model.subareas[k].case_data[:,:,1],dims = 2),
                lab = "Data",
                legend = :topleft,
                title = "KenyaCoV: fit to $(model.subareas[k].name) PCR data",
                xticks = (xticks_monthly,xticklabs_monthly),
                ylabel = "Daily cases",
                size = (800,600),dpi = 250);
        plot!(plt,pred,lab = "KenyaCoV fit",lw=2);
        return plt
end

# function plot_PCR_point_est_fit(E₀_var,θ_fitted::NamedTuple,model,k;startdate = Date(2020,2,20),enddate = Date(2021,3,1))
#         #  Get parameters
#         u0 = KenyaCoVaccines.create_initial_conditions(θ_fitted.E₀_array[k],model.subareas[k].N)
#         p = deepcopy(model.subareas[k].p)
#         p[1] = θ_fitted.β₀
#         p[41] = θ_fitted.β_home
#         p[42] = θ_fitted.β_school
#         p[43] = θ_fitted.β_other
#         p[44] = θ_fitted.β_work
#         #New variant callback
#         newvariant_time = Float64((Date(2021,1,31) - Date(2020,2,20)).value)
#     function affect_newvariant!(integrator)
#         loss_susceptible = 0.75.*integrator.u[:,7,1]
#         integrator.u[:,1,1] .+= loss_susceptible
#         integrator.u[:,7,1] .-= loss_susceptible
#         u0_nv = KenyaCoVaccines.create_initial_conditions(E₀_var,model.subareas[k].N)
#         integrator.u[:,2,1] .= u0_nv[:,2,1]
#         integrator.u[:,3:6,1] .= 0.
#     end
#     newvariant_cb = PresetTimeCallback([newvariant_time],affect_newvariant!,save_positions=(false,false))

#     cb = CallbackSet(model.subareas[k].cb,newvariant_cb)

#         #Solve dynamical system
#         l_solve = (enddate - startdate).value
#         sol = solve(model.subareas[k].prob, BS3();callback = cb, u0=u0, p=p, saveat = 1,tspan = (0,l_solve))
#         incidence = KenyaCoVaccines.get_incidence_time_array(sol)
#         PCR_pred = VectorOfArray([KenyaCoVaccines.simple_conv(incidence[:,j],model.PCR_array) for j = 1:size(incidence,2)])
#         total_PCR_pred = sum(PCR_pred[:,:],dims = 2)[:]

#         #Find days with and without neg. PCR test result data and scale prediction to match numbers of tests
#         total_days_with_neg_test = trues(length(total_PCR_pred))
#         days_without_neg_test = findall(model.subareas[k].case_data[:,1,2] .== -1)
#         days_with_neg_test = findall(model.subareas[k].case_data[:,1,2] .>= 0)
#         total_days_with_neg_test[days_without_neg_test] .= false
#         tests_per_day_from_data = sum(model.subareas[k].case_data[:,:,2],dims = 2)[:]
#         mean_tests_per_day_after_sept = mean(tests_per_day_from_data[194:end])

#         tests_per_day_for_pred = fill(mean_tests_per_day_after_sept,length(total_PCR_pred))
#         tests_per_day_for_pred[days_with_neg_test] .= tests_per_day_from_data[days_with_neg_test]

#         neg_bin_pred = total_PCR_pred.*θ_fitted.p_test_array[k].*model.relative_testing_rate[1:size(PCR_pred,1)]
#         beta_bin_pred = θ_fitted.χ_PCR_array[k].*total_PCR_pred./((θ_fitted.χ_PCR_array[k]-1).*total_PCR_pred .+ sum(model.subareas[k].N))
#         beta_bin_pred = beta_bin_pred.*tests_per_day_for_pred
#         pred = neg_bin_pred.*(.~total_days_with_neg_test) .+ beta_bin_pred.*total_days_with_neg_test

#         #Plot

#         xticks_monthly = [(Date(y,m,1) - Date(2020,2,20)).value for m = 1:12,y = 2020:2021][3:end]
#         xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2020:2021][3:end]
#         plt = scatter(sum(model.subareas[k].case_data[:,:,1],dims = 2),
#                 lab = "Data",
#                 legend = :topleft,
#                 title = "KenyaCoV: fit to $(model.subareas[k].name) PCR data",
#                 xticks = (xticks_monthly,xticklabs_monthly),
#                 ylabel = "Daily cases");
#         plot!(plt,pred,lab = "KenyaCoV fit",lw=2);
#         return plt
# end

"""
    function plot_sero_point_est_fit(θ_fitted::NamedTuple,model::CoVModel,k::Integer)
Plot the fit of a point estimate `θ_fitted` of the underlying parameters to the serological data for the `k`th
area of the `CoVModel` type `model`.
"""

function plot_sero_point_est_fit(θ_fitted::NamedTuple,model,k;startdate = Date(2020,2,20),enddate = Date(2021,3,1))
        #  Get parameters
        u0 = KenyaCoVaccines.create_initial_conditions(θ_fitted.E₀_array[k],model.subareas[k].N)
        p = deepcopy(model.subareas[k].p)
        p[1] = θ_fitted.β₀
        p[41] = θ_fitted.β_home
        p[42] = θ_fitted.β_school
        p[43] = θ_fitted.β_other
        p[44] = θ_fitted.β_work

        #Find the best fit sero-detection curve
        sero_array = vcat(model.baseline_sero_array[1:30],[(1-θ_fitted.seroreversionrate)^day for day in 1:length(model.baseline_sero_array[31:end])])

        #Solve dynamical system
        l_solve = (enddate - startdate).value
        sol = solve(model.subareas[k].prob, BS3();callback = model.subareas[k].cb, u0=u0, p=p, saveat = 1,tspan = (0,l_solve))
        incidence = KenyaCoVaccines.get_incidence_time_array(sol)
        sero_pred = VectorOfArray(model.sero_sensitivity*[KenyaCoVaccines.simple_conv(incidence[:,j],sero_array) for j = 1:size(incidence,2)])
        total_prop_inf = [sum(u[:,8,:]) for u in sol.u]./sum(model.subareas[k].N)

        #Rescale sero_pred numbers by pop. size in each age group
        pred_sero_prop_by_age = sero_pred./repeat(model.subareas[k].N',size(sero_pred,1),1)

        if sum(model.subareas[k].sero_data) > 0
                #Rescale the sero proportion prediction by the observed testing rate per age group
                num_tests_by_age = sum(model.subareas[k].sero_data[:,:,:],dims = [1,3])[:]
                num_tests = sum(model.subareas[k].sero_data[:,:,:])
                pred_sero_prop_overall = pred_sero_prop_by_age*(num_tests_by_age./num_tests)

                #Gather monthly sero-pos data
                monthly_sero_data = gather_monthly_sero_data(model.subareas[k].sero_data,[2020,2021],Date(2020,2,21))

                xticks_monthly = [(Date(y,m,1) - Date(2020,2,20)).value for m = 1:12,y = 2020:2021][3:end]
                xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2020:2021][3:end]

                plt = scatter(monthly_sero_data.monthly_x_pos,monthly_sero_data.monthly_seropos,
                                yerror = (monthly_sero_data.lb_seropos,monthly_sero_data.ub_seropos),
                                ms = 8,
                                color = :green,
                                xticks = (xticks_monthly,xticklabs_monthly),
                                lab = "Data",
                                legend = :topleft,
                                ylims = (-0.02,1),
                                ylabel = "Prop. of population",
                                title = "KenyaCoV: fit to $(model.subareas[k].name) sero data",
                                size = (800,600),dpi = 250);
                plot!(plt,pred_sero_prop_overall,
                        lab = "KenyaCoV sero. pred.",
                        lw = 2,
                        color = :green);

                plot!(plt,total_prop_inf,
                lab = "KenyaCoV cum. inf. pred.",
                lw = 2,
                color = :red);

                return plt
        else
                xticks_monthly = [(Date(y,m,1) - Date(2020,2,20)).value for m = 1:12,y = 2020:2021][3:end]
                xticklabs_monthly = [monthname(m)[1:3]*"-"*string(y-2000) for m = 1:12,y = 2020:2021][3:end]

                plt = plot(total_prop_inf,
                        xticks = (xticks_monthly,xticklabs_monthly),
                        lab = "KenyaCoV cum. inf. pred.",
                        legend = :topleft,
                        ylabel = "Prop. of population",
                        lw = 2,
                        title = "KenyaCoV: fit to $(model.subareas[k].name), no sero data",
                        color = :red,
                        size = (800,600),dpi = 250);

                return plt
        end
end

"""
    function plot_sero_age_for_monthyear_point_est(θ_fitted::NamedTuple,model,k,mnth,yr,date_first_day::Date)
Plot the fit of a point estimate `θ_fitted` of the underlying parameters to the age profile of serological data for the `k`th
area of the `CoVModel` type `model` in month `mnth`, year `yr`.
"""

function plot_sero_age_for_monthyear_point_est(θ_fitted::NamedTuple,model,k,mnth,yr,date_first_day::Date;startdate = Date(2020,2,20),enddate = Date(2021,3,1))
        #  Get parameters
        u0 = KenyaCoVaccines.create_initial_conditions(θ_fitted.E₀_array[k],model.subareas[k].N)
        p = deepcopy(model.subareas[k].p)
        p[1] = θ_fitted.β₀
        p[41] = θ_fitted.β_home
        p[42] = θ_fitted.β_school
        p[43] = θ_fitted.β_other
        p[44] = θ_fitted.β_work

        #Find the best fit sero-detection curve
        sero_array = vcat(model.baseline_sero_array[1:30],[(1-θ_fitted.seroreversionrate)^k for k in 1:length(model.baseline_sero_array[31:end])])

        #Solve dynamical system
        l_solve = (enddate - startdate).value
        sol = solve(model.subareas[k].prob, BS3();callback = model.subareas[k].cb, u0=u0, p=p, saveat = 1,tspan = (0,l_solve))
        incidence = KenyaCoVaccines.get_incidence_time_array(sol)
        sero_pred = VectorOfArray(model.sero_sensitivity*[KenyaCoVaccines.simple_conv(incidence[:,j],sero_array) for j = 1:size(incidence,2)])
        total_prop_inf = [sum(u[:,8,1]) for u in sol.u]./sum(model.subareas[k].N)

        #Rescale sero_pred numbers by pop. size in each age group and select the month-year
        data_dates = [date_first_day + Day(k-1) for k=1:size(model.subareas[k].sero_data,1)]
        daysinmonthyear = findall((month.(data_dates).== mnth).&(year.(data_dates).== yr) )

        pred_sero_prop_by_age = sero_pred./repeat(model.subareas[k].N',size(sero_pred,1),1)
        #Find average predicted sero-pos by age group in month
        pred_sero_prop_by_age = mean(pred_sero_prop_by_age[daysinmonthyear,:],dims = 1)[:]
        #Gather data for selected month-year

        pos_sero_by_age = sum(model.subareas[k].sero_data[daysinmonthyear,:,1],dims = 1)[:]
        neg_sero_by_age = sum(model.subareas[k].sero_data[daysinmonthyear,:,2],dims = 1)[:]
        prop_sero_by_age = pos_sero_by_age./(pos_sero_by_age .+ neg_sero_by_age)

        data_err = zeros(17)
        for a in 1:length(data_err)
                if pos_sero_by_age[a] + neg_sero_by_age[a] > 0
                        d = Beta(pos_sero_by_age[a]+0.5,neg_sero_by_age[a]+0.5)
                        P_a = pos_sero_by_age[a]/(pos_sero_by_age[a]+neg_sero_by_age[a])
                        data_err[a] = P_a - invlogcdf(d,log(0.025))
                end
        end

        plt = groupedbar(hcat(prop_sero_by_age,pred_sero_prop_by_age),
                        yerr = hcat(data_err,zeros(17)),
                        lab = ["KNBS data" "KenyaCoV Pred."],
                        xticks = (1:17,vcat([string((a-1)*5)*"-"*string(a*5-1) for a = 1:16 ],"80+") ),
                        legend = :topleft,
                        size = (700,400),dpi =250,
                        ylabel = "Prop. of age group",
                        title = "KenyaCoV: fit to $(model.subareas[k].name) sero data in $(monthname(mnth)), $(yr)");
        return plt
end

"""
function plot_sero_age_for_monthyear_point_est_fit_with_comparison(θ_fitted::NamedTuple,model,k,mnth,yr,date_first_day::Date;
                                                    comparison_sero_array,comparison_sero_label::String)
Plot the fit of a point estimate `θ_fitted` of the underlying parameters to the age profile of serological data for the `k`th
area of the `CoVModel` type `model` in month `mnth`, year `yr`. Also, point an additional comparison array `comparison_sero_array`.
"""
function plot_sero_age_for_monthyear_point_est_fit_with_comparison(θ_fitted::NamedTuple,model,k,mnth,yr,date_first_day::Date;
                                                    comparison_sero_array,comparison_sero_label::String,startdate = Date(2020,2,20),enddate = Date(2021,3,1))
        #  Get parameters
        u0 = KenyaCoVaccines.create_initial_conditions(θ_fitted.E₀_array[k],model.subareas[k].N)
        p = deepcopy(model.subareas[k].p)
        p[1] = θ_fitted.β₀
        p[41] = θ_fitted.β_home
        p[42] = θ_fitted.β_school
        p[43] = θ_fitted.β_other
        p[44] = θ_fitted.β_work

        #Find the best fit sero-detection curve
        sero_array = vcat(model.baseline_sero_array[1:30],[(1-θ_fitted.seroreversionrate)^k for k in 1:length(model.baseline_sero_array[31:end])])

        #Solve dynamical system
        l_solve = (enddate - startdate).value
        sol = solve(model.subareas[k].prob, BS3();callback = model.subareas[k].cb, u0=u0, p=p, saveat = 1, tspan = (0,l_solve))
        incidence = KenyaCoVaccines.get_incidence_time_array(sol)
        sero_pred = VectorOfArray(model.sero_sensitivity*[KenyaCoVaccines.simple_conv(incidence[:,j],sero_array) for j = 1:size(incidence,2)])
        total_prop_inf = [sum(u[:,8,:]) for u in sol.u]./sum(model.subareas[k].N)

        #Rescale sero_pred numbers by pop. size in each age group and select the month-year
        data_dates = [date_first_day + Day(k-1) for k=1:size(model.subareas[k].sero_data,1)]
        daysinmonthyear = findall((month.(data_dates).== mnth).&(year.(data_dates).== yr) )

        pred_sero_prop_by_age = sero_pred./repeat(model.subareas[k].N',size(sero_pred,1),1)
        #Find average predicted sero-pos by age group in month
        pred_sero_prop_by_age = mean(pred_sero_prop_by_age[daysinmonthyear,:],dims = 1)[:]

        plt = groupedbar(hcat(pred_sero_prop_by_age,comparison_sero_array),
                        lab = ["KenyaCoV Pred." comparison_sero_label],
                        xticks = (1:17,vcat([string((a-1)*5)*"-"*string(a*5-1) for a = 1:16 ],"80+") ),
                        legend = :topleft,
                        size = (700,400),dpi =250,                        ylabel = "Prop. of age group",
                        title = "KenyaCoV: fit to $(model.subareas[k].name) sero data in $(monthname(mnth)), $(yr)");
        return plt
end

function movingaverage(X::Vector,numofele::Int)
    BackDelta = div(numofele,2)
    ForwardDelta = isodd(numofele) ? div(numofele,2) : div(numofele,2) - 1
    len = length(X)
    Y = similar(X)
    for n = 1:len
        lo = max(1,n - BackDelta)
        hi = min(len,n + ForwardDelta)
        Y[n] = mean(X[lo:hi])
    end
    return Y
end


text: struct NotImplementedError <:Exception end


text:
#
# This file is part of the Simulate.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

using Simulate, Printf, Random

function simple(c::Clock, input::Channel, output::Channel, name, id, op)
    token = take!(input)         # take something from the input
    now!(c, SF(println, @sprintf("%5.2f: %s %d took token %d", tau(c), name, id, token)))
    d = delay!(c, rand())        # after a delay
    put!(output, op(token, id))  # put it out with some op applied
end

clk = Clock()      # create a clock
Random.seed!(123)  # seed the random number generator

ch1 = Channel(32)  # create two channels
ch2 = Channel(32)

for i in 1:2:8    # create and register 8 SimProcesses SP
    process!(clk, SP(i, simple, ch1, ch2, "foo", i, +))
    process!(clk, SP(i+1, simple, ch2, ch1, "bar", i+1, *))
end

put!(ch1, 1)    # put first token into channel 1
yield()         # let the first task take it
run!(clk, 10)   # and run for 10 time units


text:
"""
    struct DemoPage <: Any
    DemoPage(root::String)

Constructs a demo page object.

# Fields

Besides the root path to the demo page folder `root`, this struct has some other fields:

* `title`: page title
* `template`: template content of the demo page.
* `sections`: demo sections found in `root`

# Configuration

You can manage an extra `config.json` file to customize rendering of a demo page.
Supported items are:

* `order`: specify the sections order. By default, it's case-insensitive alphabetic order.
* `template`: path to template filename. By default, it's `"index.md"`. The content of the template file should has one and only one `{{{democards}}}`.
* `theme`: specify which card theme should be used to generate the index page. If not specified, it
  will default to `nothing`.
* `title`: specify the title of this demo page. By default, it's the folder name of `root`. Will be override by `template`.
* `properties`: a dictionary of properties that can be propagated to its children items. The same properties in
  the children items, if exist, have higher priority.

The following is an example of `config.json`:

```json
{
    "template": "template.md",
    "theme": "grid",
    "order": [
        "basic",
        "advanced"
    ],
    "properties": {
        "notebook": "false",
        "julia": "1.6",
        "author": "Johnny Chen"
    }
}
```

# Examples

The following is the simplest folder structure of a `DemoPage`:

```text
demos
└── basic
    ├── demo_1.md
    ├── demo_2.md
    ├── demo_3.md
    ├── demo_4.md
    ├── demo_5.md
    ├── demo_6.md
    ├── demo_7.md
    └── demo_8.md
```

!!! note

    A `DemoPage` doesn't manage demo files directly, so here you'll need
    a DemoSection `basic` to manage them.

The following is a typical folder structure of a `DemoPage`:

```text
demos
├── advanced
│   ├── advanced_demo_1.md
│   └── advanced_demo_2.md
├── basic
│   ├── part1
│   │   ├── basic_demo_1_1.md
│   │   └── basic_demo_1_2.md
│   └── part2
│       ├── config.json
│       ├── basic_demo_2_1.md
│       └── basic_demo_2_2.md
├── config.json
└── template.md
```

!!! warning

    A section should only hold either subsections or demo files. A folder that has both subfolders and demo files (e.g., `*.md`) is invalid.

See also: [`MarkdownDemoCard`](@ref DemoCards.MarkdownDemoCard), [`DemoSection`](@ref DemoCards.DemoSection)
"""
mutable struct DemoPage
    root::String
    sections::Vector{DemoSection}
    template::String
    theme::Union{Nothing, String}
    title::String
    # These properties will be shared by all children of it during build time
    properties::Dict{String, Any}
end

basename(page::DemoPage) = basename(page.root)

function DemoPage(root::String)::DemoPage
    root = replace(root, r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
    isdir(root) || throw(ArgumentError("page root does not exist: $(root)"))
    root = rstrip(root, '/') # otherwise basename(root) will returns `""`

    path = joinpath.(root, filter(x->!startswith(x, "."), readdir(root))) # filter out hidden files
    section_paths = filter(x->isdir(x)&&!(basename(x) in ignored_dirnames),
                           path)

    if isempty(section_paths)
        # This is an edge and trivial case that we need to disallow; we want the following
        # folder structure be uniquely parsed as a demo section and not a demo page without
        # ambiguity.
        # ```
        # demos/
        # └── demo1.jl
        # ```
        msg = "can not find a valid page structure in page dir \"$root\"\nit should have at least one folder as section inside it"
        throw(ArgumentError(msg))
    end

    json_path = joinpath(root, config_filename)
    json_config = isfile(json_path) ? JSON.parsefile(json_path) : Dict()

    template_file = joinpath(root, get(json_config, "template", template_filename))
    config = parse(Val(:Markdown), template_file)
    config = merge(json_config, config) # template has higher priority over config file

    sections = filter(map(DemoSection, section_paths)) do sec
        empty_section = isempty(sec.cards) && isempty(sec.subsections)
        if empty_section
            @warn "Empty section detected, remove from the demo page tree." section=relpath(sec.root, root)
            return false
        else
            return true
        end
    end
    isempty(sections) && error("Empty demo page, you have to add something.")

    page = DemoPage(root, sections, "", nothing, "", Dict{String, Any}())
    page.theme = load_config(page, "theme"; config=config)

    section_orders = load_config(page, "order"; config=config)
    section_orders = map(sections) do sec
        findfirst(x-> x == basename(sec.root), section_orders)
    end
    ordered_sections = sections[section_orders]

    title = load_config(page, "title"; config=config)
    page.sections = ordered_sections
    page.title = title

    # default template requires a title
    template = load_config(page, "template"; config=config)
    page.template = template

    if haskey(config, "properties")
        properties = config["properties"]::Dict
        merge!(page.properties, properties)
    end
    return page
end

function load_config(page::DemoPage, key; config=Dict())
    if isempty(config)
        json_path = joinpath(page.root, config_filename)
        json_config = isfile(json_path) ? JSON.parsefile(json_path) : Dict()

        template_file = joinpath(page.root, get(json_config, "template", template_filename))
        config = parse(Val(:Markdown), template_file)
        config = merge(json_config, config) # template has higher priority over config file
    end
    # config could still be an empty dictionary

    if key == "order"
        haskey(config, key) || return get_default_order(page)

        order = config[key]
        validate_order(order, page)
        return order
    elseif key == "title"
        return get(config, key, get_default_title(page))
    elseif key == "template"
        return load_template(page, config)
    elseif key == "theme"
        theme = get(config, key, nothing)
        if !isnothing(theme) && lowercase(theme) == "nothing"
            theme = nothing
        end
        return theme
    else
        throw(ArgumentError("Unrecognized key $(key) for DemoPage"))
    end
end

function load_template(page, config)
    key = "template"
    if haskey(config, key)
        # windows compatibility
        template_name = replace(config[key],
                            r"[/\\]" => Base.Filesystem.path_separator)
        template_file = joinpath(page.root, template_name)
    else
        template_file = joinpath(page.root, template_filename)
    end

    isfile(template_file) || return get_default_template(page)

    check_ext(template_file, :markdown)
    if 1 != sum(occursin.("{{{democards}}}", readlines(template_file)))
        throw(ArgumentError("invalid template file $(template_file): it should has one and only one {{{democards}}}"))
    end
    return read(template_file, String)
end

get_default_order(page::DemoPage) =
    sort(basename.(page.sections), by = x->lowercase(x))

function get_default_template(page::DemoPage)
    header = "# $(page.title)\n\n"
    # TODO: by doing this we loss the control on section-level template
    content = "{{{democards}}}" # render by Mustache
    footer = ""
    return header * content * footer
end

# page utils

"""
    infer_pagedir(card_path; rootdir="")

Given a demo card path, infer the *outmost* dir path that makes it a valid demo page. If it fails to
find such dir path, return `nothing`.

The inference is done recursively, `rootdir` sets a stop condition for the inference process.

!!! warning
    This inference may not be the exact page dir in trivial cases, for example:

    ```
    testdir/
    └── examples
        └── sections
            └── demo1.md
    ```

    Both `testdir` and `examples` can be valid dir path for a demo page, this function would
    just return `testdir` as it is the outmost match.

"""
function infer_pagedir(path; rootdir::String=first(splitdrive(path)))
    # the last root path that `DemoPage(root)` successfully parses
    if is_demopage(path)
        root, dir = splitdir(path) # we need the outmost folder, so not return right here
    elseif is_demosection(path)
        root, dir = splitdir(path)
    elseif is_democard(path)
        root, dir = splitdir(dirname(path))
    else
        return nothing
    end

    while true
        if is_demopage(root)
            root, dir = splitdir(root)
        else
            pagedir = joinpath(root, dir)
            return is_demopage(pagedir) ? pagedir : nothing # in case it fails at the first try
        end

        root in push!(["/", ""], rootdir) && return nothing
    end
end

function is_demopage(dir)
    try
        # if fails to parse, then it is not a valid demo page
        @suppress_err DemoPage(dir)
        return true
    catch err
        @debug err
        return false
    end
end


text:
# Equation 33 of https://www.mdpi.com/2227-7390/6/2/16/htm

# Loaded Packages
using FdeSolver
using Plots
using SpecialFunctions

# Parameters
tSpan = [0, 5]         # Time Span
y0 = [1, 0.5, 0.3]     # Initial values
β = [0.5, 0.2, 0.6]    # Order of derivation

# Definition of the System
function F(t, n, β, y)

    F1 = 1 / sqrt(pi) * (((y[n, 2] - 0.5) * (y[n, 3] - 0.3))^(1 / 6) + t[n]^(1 / 2))
    F2 = gamma(2.2) * (y[n, 1] - 1)
    F3 = gamma(2.8) / gamma(2.2) * (y[n, 2] - 0.5)

    return [F1, F2, F3]

end

# Numerical Solution
t, Yapp = FDEsolver(F, tSpan, y0, β, nothing, nc = 5)

# Plot
plot(t, Yapp, linewidth = 5, title = "Solution of system 33",
     xaxis = "Time (t)", yaxis = "y(t)", label = "Approximation")

plot!(t, t -> (t .+ 1), lw = 3, ls = :dash, color= "red", label = "Exact solution")
plot!(t, t -> (t.^1.2 .+ 0.5), lw = 3, ls = :dash, color= "cyan", label = "Exact solution")
plot!(t, t -> (t.^1.8 .+ 0.3), lw = 3, ls = :dash, color= "black" ,label = "Exact solution")


text:
using Documenter, SDWBA

makedocs(
  sitename = "SDWBA.jl",
  authors = "Sam Urmy",
  # assets = ["assets/style.css"],
  pages = [
    "Introduction" => "index.md",
    "API reference" => "API.md",
    "Included Scatterering Models" => "Models.md"
  ]
)

deploydocs(
    repo = "github.com/ElOceanografo/SDWBA.jl.git",
)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

import Base.show

##### Observatory
"""
Type holding information about an observing site.  Its fields are:

* `name`: the name of the site
* `latitude`: North-ward latitude of the site in degrees
* `longitude`: East-ward longitude of the site in degrees
* `altitude`: altitude of the site in meters
* `tz`: the number of hours of offset from UTC
"""
struct Observatory
    name::String
    latitude::Float64
    longitude::Float64
    altitude::Float64
    tz::Float64 # There are non-integer time zones
    # Define constructor that automatically converts longitude and latitude with
    # "ten", for convenience.
    Observatory(name, lat, long, alt, tz) =
        new(String(name), Float64(ten(lat)), Float64(ten(long)),
            Float64(float(alt)), Float64(ten(tz)))
end

# New type representation
function show(io::IO, obs::Observatory)
    println(io, "Observatory: ", obs.name)
    println(io, "latitude:    ", obs.latitude, "°N")
    println(io, "longitude:   ", obs.longitude, "°E")
    println(io, "altitude:    ", obs.altitude, " m")
    tzdec, tzint = modf(obs.tz)
    print(io,   "time zone:   ", @sprintf("UTC%+d", tzint),
          tzdec == 0 ? "" : @sprintf(":%d", abs(tzdec*60)))
end

##### Planet
"""
Type holding information about a planet.  Its fields are:

Designation:

* `name`: the name

Physical characteristics:

* `radius`: mean radius in meters
* `eqradius`: equatorial radius in meters
* `polradius`: polar radius in meters
* `mass`: mass in kilogram

Orbital characteristics (epoch J2000):

* `ecc`: eccentricity of the orbit
* `axis`: semi-major axis of the orbit in meters
* `period`: sidereal orbital period in seconds

Position characteristics (epoch J2000):

* `inc`: inclination in degrees
* `asc_long`: longitude of the ascending node in degrees
* `per_long`: longitude of perihelion in degrees
* `mean_long`: mean longitude in degrees
"""
struct Planet
    name::String
    radius::Float64
    eqradius::Float64
    polradius::Float64
    mass::Float64
    ecc::Float64
    axis::Float64
    period::Float64
    inc::Float64
    asc_long::Float64
    per_long::Float64
    mean_long::Float64
end

# New type representation
function show(io::IO, pl::Planet)
    println(io, "Planet:                    ", uppercasefirst(pl.name))
    println(io, "mean radius:               ", pl.radius, " m")
    println(io, "equatorial radius:         ", pl.eqradius, " m")
    println(io, "polar radius:              ", pl.polradius, " m")
    println(io, "mass:                      ", pl.mass, " kg")
    println(io, "eccentricity:              ", pl.ecc)
    println(io, "semi-major axis:           ", pl.axis, " m")
    println(io, "period:                    ", pl.period, " s")
    println(io, "inclination                ", pl.inc, " °")
    println(io, "longitude of ascending node", pl.asc_long, " °")
    println(io, "longitude of perihelion    ", pl.per_long, " °")
    print(io,   "mean longitude             ", pl.mean_long, " °")
end

export Observatory, Planet


text:
module TestSimple

include("preamble.jl")

# Test that task group are synchronized and the task context can be
# passed around across threads/tasks.
@testset "wait" begin
    results = SyncSeq()

    bg(context, id) = @in context @go @go push!(results, id)

    @taskgroup begin
        @go bg(_, 1)
        @go bg(_, 2)
        @go bg(_, 3)
    end

    @test sort(collect(results)) == 1:3
end

# Test cancellation on error.
@testset "cancel on error" begin
    results = SyncSeq()

    function bg(context, id)
        @in context begin
            @go begin
                push!(results, ("before sleep", id))
                @await tickingfor(_, 1)
                push!(results, ("after sleep", id))
            end
            if id == 1
                @await ErrorException("terminate")  # "throw" an error
            end
        end
    end

    @taskgroup begin
        @go bg(_, 1)
        @go bg(_, 2)
    end

    @test sort(collect(results)) == [("before sleep", 1), ("before sleep", 2)]
end

# Test manual cancellation
@testset "manual cancel" begin
    results = SyncSeq()

    function bg(context, id)
        @in context begin
            @go begin
                push!(results, ("before sleep", id))
                @await tickingfor(_, 1)
                push!(results, ("after sleep", id))
            end
            if id == 1
                cancel!(context)
            end
        end
    end

    @taskgroup begin
        @go bg(_, 1)
        @go bg(_, 2)
    end

    @test sort(collect(results)) == [("before sleep", 1), ("before sleep", 2)]
end

# Test that cancellation of `@cancelscope` does not leak out.
@testset "`@cancelscope`" begin
    results = SyncSeq()

    function bg(context, id)
        @in context begin
            subctx = @cancelscope begin
                @go begin
                    push!(results, ("before sleep", id))
                    @await tickingfor(_, 1)
                    push!(results, ("after sleep", id))
                end
            end
            if id == 1
                # Manually cancel only this scope:
                cancel!(subctx)
                # Note: `@await ErrorException("terminate")` inside
                # `@cancelscope` does not work because it cancels
                # everything; i.e., the error bubbles up until the
                # root `@taskgroup` which in turn cancels everything,
                # not just this `@cancelscope`.
            end
        end
    end

    @taskgroup begin
        @go bg(_, 1)
        @go bg(_, 2)
    end

    @test sort(collect(results)) ==
        [("after sleep", 2), ("before sleep", 1), ("before sleep", 2)]
end

end  # module


text:
function write_doc(docformat::LaTeX2PDF, doc, rendered, out_path)
    cd_back = let d = pwd(); () -> cd(d); end
    cd(doc.cwd)
    try
        tex_path = basename(out_path)
        write(tex_path, rendered)
        cmds = copy(docformat.latex_cmd)
        push!(cmds, tex_path)
        cmd = Cmd(cmds)
        run(cmd); run(cmd) # XXX: is twice enough for every case ?
    catch
        @warn "Error converting document to pdf. Try running latex manually"
        rethrow()
    finally
        cd_back()
    end

    return get_out_path(doc, out_path, "pdf")
end