with(plots):
with(plottools):
with(Student[Calculus1]):

# Define the 1-soliton solution
u := (x, t, k, delta) -> (k^2/2)*sech(1/2*(k*x - k^3*t + delta))^2:

# Interactive plot with sliders for t, k, delta
interactive_plot := Explore(
    plot(u(x, t, k, delta), x = -10 .. 10, color = blue, thickness = 2, axes = boxed, labels = ["x", "u(x,t)"]),
    t = -5 .. 5,
    k = 0.5 .. 3,
    delta = -5 .. 5
);

BlaschkeAnalyse := proc(xfun, yfun, tinterval := 0..2*Pi, naam := "Curve", kleur := blue, nulsplits := 3)
local dx, dy, ddx, ddy, kappa, f, i, sol, t1, t2, x1, y1, x2, y2, k1, k2,
      curveplot, tangent_plot, tangents, points, labels, table, solutions, L, norm, vx0, vy0, x0, y0, r, p;

  uses plots, plottools ;
  
  print("�� Analysis of ", naam);

  # Derivatives
  dx := unapply(diff(xfun(t), t), t):
  dy := unapply(diff(yfun(t), t), t):
  ddx := unapply(diff(dx(t), t), t):
  ddy := unapply(diff(dy(t), t), t):

  # Curvature
  kappa := t -> abs(dx(t)*ddy(t) - dy(t)*ddx(t)) / ((dx(t)^2 + dy(t)^2)^(3/2)):

  # Plot the curve
  curveplot := plot([xfun(t), yfun(t), t = tinterval], color = kleur, thickness=2, labels = ["x", "y"]):

  # Difference function f(t) = kappa(t) - kappa(t+Pi)
  f := t -> evalf(kappa(t) - kappa(t + Pi)):

  # Search for roots numerically in split intervals
  solutions := []:
  for i from 0 to nulsplits-1 do
    sol := fsolve(f(t) = 0, t = op(1, tinterval) + i*Pi/nulsplits + 0.01 .. op(1, tinterval) + (i+1)*Pi/nulsplits - 0.01):
    if type(sol, numeric) then
      solutions := [op(solutions), sol]:
    end if;
  end do:

  print("✅ Blaschke pairs found: ", nops(solutions));

  # Visualization elements
  tangents := []:
  points := []:
  labels := []:
  table := []:
  L := 3:

  tangent_plot := proc(t0)
    x0 := evalf(xfun(t0)):
    y0 := evalf(yfun(t0)):
    vx0 := evalf(dx(t0)):
    vy0 := evalf(dy(t0)):
    norm := sqrt(vx0^2 + vy0^2):
    vx0 := vx0 / norm:
    vy0 := vy0 / norm:
    return line([x0 - L*vx0, y0 - L*vy0], [x0 + L*vx0, y0 + L*vy0], color=red, thickness=2):
  end proc:

  for i from 1 to nops(solutions) do
    t1 := evalf(solutions[i]):
    t2 := evalf(t1 + Pi):
    x1 := evalf(xfun(t1)):
    y1 := evalf(yfun(t1)):
    x2 := evalf(xfun(t2)):
    y2 := evalf(yfun(t2)):
    k1 := evalf(kappa(t1)):
    k2 := evalf(kappa(t2)):

    # Tangents
    tangents := [op(tangents), tangent_plot(t1), tangent_plot(t2)]:

    # Points
    points := [op(points),
      pointplot([[x1, y1]], symbol=solidcircle, symbolsize=14, color=red),
      pointplot([[x2, y2]], symbol=solidcircle, symbolsize=14, color=red)
    ]:

    # Labels
    labels := [op(labels),
      textplot([[x1, y1, cat("P", i)]], font=[Times, Bold, 14], align=above, color=red),
      textplot([[x2, y2, cat("P", i, "'")]], font=[Times, Bold, 14], align=above, color=red)
    ]:

    # Table rows
    table := [op(table),
      sprintf("P%d:    t = %.4f   (x,y) = (%.4f, %.4f)   κ = %.4f", i, t1, x1, y1, k1),
      sprintf("P%d':   t = %.4f   (x,y) = (%.4f, %.4f)   κ = %.4f", i, t2, x2, y2, k2),
      "------------------------------------------------------------"
    ]:
  end do:

  # Display table
  print("�� Table of Blaschke pairs:");
  for r in table do
    print(r);
  end do:

  # Show everything
  return display([curveplot, op(tangents), op(points), op(labels)], scaling=constrained,
    title=cat("Blaschke Analysis of ", naam));

end proc:

# Voorbeeldgebruik: superellips
sgn := t -> piecewise(t < 0, -1, t = 0, 0, 1):
a := 3: b := 2: n := 2.5:
x := t -> a * sgn(cos(t)) * abs(cos(t))^n:
y := t -> b * sgn(sin(t)) * abs(sin(t))^n:

BlaschkeAnalyse(x, y, 0..2*Pi, naam="Superellips (n=2.5)", magenta);

BlaschkeAnalyse(
    t -> 5*cos(t),
    t -> (5 + sin(6*t))*sin(t),
    0..2*Pi,
    "Voorbeeldkromme",
    blue,
    3
);

restart;
with(plots):

### Voorwaarden ###
assume(x::real, y::real);
additionally(cosh(2*y) - cos(2*x) <> 0);  # Voorkom singulariteit

### Stap 1: Definieer u(x,y) ###
u := sin(2*x)/(cosh(2*y) - cos(2*x));
print("Stap 1: Gegeven reëel deel u(x,y):");
print(u);

### Stap 2: Controleer harmoniciteit ###
laplacian_u := diff(u, x$2) + diff(u, y$2);
simplified_laplacian := simplify(laplacian_u) assuming cosh(2*y) - cos(2*x) <> 0;
print("Stap 2: Laplaciaan van u (moet 0 zijn):");
print(simplified_laplacian);

if simplified_laplacian = 0 then
    print("u is harmonisch ✅");
else
    print("u is NIET harmonisch ❌");
    return;
end if;

### Stap 3: Cauchy-Riemann - bereken afgeleiden ###
ux := diff(u, x):
uy := diff(u, y):
print("Stap 3: Partiële afgeleiden van u:");
print("∂u/∂x =", simplify(ux));
print("∂u/∂y =", simplify(uy));

### Stap 4: Integreer voor v(x,y) ###
v_x := -uy;
v := int(v_x, x) + C1;  # C1 is integratieconstante
print("Stap 4: v(x,y) (inclusief C1):");
print(simplify(v));

### Stap 5: Harmonisch geconjugeerde functie ###
print("Stap 5: C1 = constante, dus behouden");
f_xy := u + I*v;
print("Stap 6: Analytische functie f(x,y):");
print(simplify(f_xy));

### Stap 6: Vergelijk met cot(z) ###
f_xy := subs(C1 = 0, f_xy);  # Constante = 0 voor vergelijking
cot_z := sin(2*x)/(cosh(2*y) - cos(2*x)) - I*sinh(2*y)/(cosh(2*y) - cos(2*x));
difference := simplify(f_xy - cot_z);
print("Stap 7: Verschil tussen f(x,y) en cot(z):");
print(difference);

if difference = 0 then
    print("✅ f(z) = cot(z)");
else
    print("❌ f(z) ≠ cot(z), verschil ≠ 0 (controle C1?)");
end if;

### Stap 8: Plot ter verificatie ###
plot3d(u, x = -Pi .. Pi, y = -1 .. 1, title = "Reëel deel u(x,y)");
plot3d(Im(subs(C1=0, f_xy)), x = -Pi .. Pi, y = -1 .. 1, title = "Imaginaire deel v(x,y) (C1=0)");
 

SEIRTP_R0_sensitivity_report := proc(params::list, normalize::boolean := true)
description "R₀ Sensitivity Analysis – Table + Barplot + Dynamic Policy Advice for SEIRTP model with compartment effects";

global comp_matrix_export, comp_matrix_labels;

local Lambda, alpha, beta, r, xi, K, eta, phi, mu, sigma, tau, theta, delta, gamma_;
local Ro, Ro_num, dp, Sp, p, i, paramlist, param_names, values, label_list, maxval;
local param_values, n_params, barplot, interpretation, recommendations, threshold;
local dir, advice, v, comp_labels, comp_effects, all_effects, j;

# Full symbolic parameter list
paramlist := [Lambda, alpha, beta, r, xi, K, eta, phi, mu, sigma, tau, theta, delta, gamma_];
param_names := ["Λ", "α", "β", "r", "ξ", "K", "η", "φ", "μ", "σ", "τ", "θ", "δ", "γ"];
n_params := nops(paramlist);
param_values := [seq(rhs(params[i]), i=1..n_params)];

# R₀ formula
Ro := (Lambda*alpha*beta*r*xi*(-K*eta*phi + K*eta*r - r)) /
      (mu*(K*phi*sigma - K*r*sigma - r*tau)*
      (K*beta*eta*phi - K*beta*eta*r + beta*r + mu*r + r*theta)*(delta + gamma_ + mu));

printf("\n�� SEIRTP-Parameterverklaring / Parameter Descriptions:\n");
printf("╔════════════╦════════════════════════════════════════════════════════════╗\n");
printf("║ Parameter  ║ Beschrijving (NL) / Description (EN)                      ║\n");
printf("╠════════════╬════════════════════════════════════════════════════════════╣\n");
printf("║ Λ (Lambda) ║ Instroom / geboorteratio – inflow rate / birth rate       ║\n");
printf("║ α (alpha)  ║ Infectieratio bij contact – infection probability         ║\n");
printf("║ β (beta)   ║ Contactfrequentie – contact rate                          ║\n");
printf("║ r          ║ Detectieratio / testintensiteit – detection rate          ║\n");
printf("║ ξ (xi)     ║ Tracing efficiëntie – tracing effectiveness               ║\n");
printf("║ K          ║ Zorgcapaciteit – healthcare system capacity               ║\n");
printf("║ η (eta)    ║ Snelheid van ziek worden – progression to infectious      ║\n");
printf("║ φ (phi)    ║ Behandelingsimpact – treatment effect on transmission     ║\n");
printf("║ μ (mu)     ║ Natuurlijke sterfte – natural death rate                  ║\n");
printf("║ σ (sigma)  ║ Sterfte bij behandeling – death during treatment          ║\n");
printf("║ τ (tau)    ║ Behandelduur – duration of treatment                      ║\n");
printf("║ θ (theta)  ║ Terugval / herinfectie – relapse or reinfection rate      ║\n");
printf("║ δ (delta)  ║ Infectie-specifieke sterfte – infection-induced mortality ║\n");
printf("║ γ (gamma)  ║ Herstelkans – recovery rate                               ║\n");
printf("╚════════════╩════════════════════════════════════════════════════════════╝\n\n");

Ro_num := evalf(subs(paramlist =~ param_values, Ro));
if Ro_num = 0 then error "R₀ = 0, normalization not possible"; end if;

printf("�� R₀ computed = %.5f\n", Ro_num);
if abs(Ro_num - 1.0) < 0.01 then
    printf("⚖️  R₀ ≈ 1 → Epidemic threshold.\n");
elif Ro_num < 1.0 then
    printf("✅ R₀ < 1 → Outbreak likely to die out.\n");
else
    printf("�� R₀ > 1 → Epidemic spread possible.\n");
end if;

# Sensitivity values
values := [];
for i from 1 to n_params do
    p := paramlist[i];
    dp := evalf(subs(paramlist =~ param_values, diff(Ro, p)));
    Sp := dp * param_values[i] / Ro_num;
    values := [op(values), Sp];
end do;

if normalize then
    maxval := max(seq(abs(v), v in values));
    values := map(v -> v / maxval, values);
end if;

# Interpretation
interpretation := proc(val)
    if abs(val) >= 0.9 then return "�� Strong effect";
    elif abs(val) >= 0.4 then return "�� Moderate influence";
    elif abs(val) >= 0.01 then return "�� Small effect";
    else return "�� Negligible"; end if;
end proc:

# Sensitivity Table
printf("\n�� %s Sensitivity Index of R₀:\n", `if`(normalize, "Normalized", "Original"));
printf("--------------------------------------------------------------------------\n");
printf("%-13s | %-10s | %-30s\n", "Parameter", "Sₚ", "Interpretation");
printf("--------------------------------------------------------------------------\n");
for i from 1 to n_params do
    printf("%-13s | %+10.5f | %-30s\n", param_names[i], values[i], interpretation(values[i]));
end do;
printf("--------------------------------------------------------------------------\n");

# Policy Advice
threshold := 0.4;
printf("\n�� Dynamic Policy Implications:\n");
printf("--------------------------------------------------------------------------\n");
printf("%-13s | %-10s | %-30s\n", "Parameter", "Direction", "Policy Advice");
printf("--------------------------------------------------------------------------\n");
for i from 1 to n_params do
    if abs(values[i]) >= threshold then
        if values[i] > 0 then
            dir := "↑ Increase"; advice := "⚠️ Reduce to lower R₀";
        else
            dir := "↓ Decrease"; advice := "✅ Increase to lower R₀";
        end if;
        printf("%-13s | %-10s | %-30s\n", param_names[i], dir, advice);
    end if;
end do;
printf("--------------------------------------------------------------------------\n");

# Compartment Effects Setup
comp_labels := ["S", "E", "I", "R", "P", "T"];
all_effects := table([
    Lambda = [1, 0, 0, 0, 0, 0],
    alpha  = [-1, 1, 1, 0, 0, 0],
    beta   = [0, -1, 1, 0, 0, 0],
    r      = [0, 0, 0, 0, 0, 1],
    xi     = [0, 0, 0, 0, 1, 0],
    K      = [0, 0, 0, 0, 0, 1],
    eta    = [0, -1, -1, 0, 0, 0],
    phi    = [0, 0, 0, 0, -1, 0],
    mu     = [-1, -1, -1, -1, -1, -1],
    sigma  = [0, 0, 0, 0, -1, 0],
    tau    = [0, 0, 0, 0, -1, 0],
    theta  = [0, -1, 0, 1, 0, 0],
    delta  = [0, 0, -1, 0, 0, 0],
    gamma_ = [0, 0, -1, 1, 0, 0]
]);

# Corrected matrix build using paramlist symbols
comp_effects := [seq(all_effects[paramlist[i]], i = 1..n_params)];
comp_matrix_labels := param_names;

# Export matrix
comp_matrix_export := comp_effects;

# Display Matrix
printf("\n�� Compartmental Effect Matrix (Qualitative):\n");
printf("--------------------------------------------------------------------------\n");
printf("%-10s | %3s %3s %3s %3s %3s %3s\n", "Parameter", op(comp_labels));
printf("--------------------------------------------------------------------------\n");
for i from 1 to nops(comp_effects) do
    printf("%-10s |", comp_matrix_labels[i]);
    for j to 6 do
        if comp_effects[i][j] = 0 then printf(" %3s", ".");
        elif comp_effects[i][j] > 0 then printf(" %3s", "+");
        else printf(" %3s", "-"); end if;
    end do;
    printf("\n");
end do;
printf("--------------------------------------------------------------------------\n");

# Optional Plot: UITGEZET werkt nog niet
##label_list := [seq(sprintf("%s (%s)", param_names[i], paramlist[i]), i=1..n_params)];
#barplot := plots:-barplot(label_list, values, orientation=horizontal,
             #title=sprintf("Sensitivity of R₀ (%s)", `if`(normalize, "Normalized", "Original")),
             #width=0.6, color="CornflowerBlue", axes=boxed,
            # labels=["Parameter", "Sensitivity Sₚ"]);

#return barplot;

end proc:

params := [
  Lambda = 0.0133, alpha = 0.07954551, beta = 0.567, r = 0.076, xi = 0.0487,
  K = 300, eta = 0.009, phi = 0.09, mu = 0.00056, sigma = 0.000456,
  tau = 0.0900982, theta = 0.009, delta = 0.9, gamma_ = 0.00917
];
 

SEIRTP_R0_sensitivity_report(params, true);   # Genormaliseerde sensitiviteit

�� SEIRTP-Parameterverklaring / Parameter Descriptions:
╔════════════╦════════════════════════════════════════════════════════════╗
║ Parameter  ║ Beschrijving (NL) / Description (EN)                      ║
╠════════════╬════════════════════════════════════════════════════════════╣
║ Λ (Lambda) ║ Instroom / geboorteratio – inflow rate / birth rate       ║
║ α (alpha)  ║ Infectieratio bij contact – infection probability         ║
║ β (beta)   ║ Contactfrequentie – contact rate                          ║
║ r          ║ Detectieratio / testintensiteit – detection rate          ║
║ ξ (xi)     ║ Tracing efficiëntie – tracing effectiveness               ║
║ K          ║ Zorgcapaciteit – healthcare system capacity               ║
║ η (eta)    ║ Snelheid van ziek worden – progression to infectious      ║
║ φ (phi)    ║ Behandelingsimpact – treatment effect on transmission     ║
║ μ (mu)     ║ Natuurlijke sterfte – natural death rate                  ║
║ σ (sigma)  ║ Sterfte bij behandeling – death during treatment          ║
║ τ (tau)    ║ Behandelduur – duration of treatment                      ║
║ θ (theta)  ║ Terugval / herinfectie – relapse or reinfection rate      ║
║ δ (delta)  ║ Infectie-specifieke sterfte – infection-induced mortality ║
║ γ (gamma)  ║ Herstelkans – recovery rate                               ║
╚════════════╩════════════════════════════════════════════════════════════╝

�� R₀ computed = 1.54099
�� R₀ > 1 → Epidemic spread possible.

�� Normalized Sensitivity Index of R₀:
--------------------------------------------------------------------------
Parameter     | Sₚ       | Interpretation                
--------------------------------------------------------------------------
Λ            |   +0.39683 | �� Small effect             
α            |   +0.39683 | �� Small effect             
β            |   +0.00442 | �� Negligible               
r             |   -1.00000 | �� Strong effect            
ξ            |   +0.39683 | �� Small effect             
K             |   +0.15556 | �� Small effect             
η            |   +0.00147 | �� Negligible               
φ            |   +1.00000 | �� Strong effect            
μ            |   -0.39733 | �� Small effect             
σ            |   +0.15409 | �� Small effect             
τ            |   -0.55091 | �� Moderate influence       
θ            |   -0.00416 | �� Negligible               
δ            |   -0.39258 | �� Small effect             
γ            |   -0.00400 | �� Negligible               
--------------------------------------------------------------------------

�� Dynamic Policy Implications:
--------------------------------------------------------------------------
Parameter     | Direction  | Policy Advice                 
--------------------------------------------------------------------------
r             | ↓ Decrease | ✅ Increase to lower R₀    
φ            | ↑ Increase | ⚠️ Reduce to lower R₀   
τ            | ↓ Decrease | ✅ Increase to lower R₀    
--------------------------------------------------------------------------

�� Compartmental Effect Matrix (Qualitative):
--------------------------------------------------------------------------
Parameter  |   S   E   I   R   P   T
--------------------------------------------------------------------------
Λ         |   +   .   .   .   .   .
α         |   -   +   +   .   .   .
β         |   .   -   +   .   .   .
r          |   .   .   .   .   .   +
ξ         |   .   .   .   .   +   .
K          |   .   .   .   .   .   +
η         |   .   -   -   .   .   .
φ         |   .   .   .   .   -   .
μ         |   -   -   -   -   -   -
σ         |   .   .   .   .   -   .
τ         |   .   .   .   .   -   .
θ         |   .   -   .   +   .   .
δ         |   .   .   -   .   .   .
γ         |   .   .   -   +   .   .
--------------------------------------------------------------------------

SEIRTP_R0_sensitivity_report(params, false);  # Originele sensitiviteit

�� SEIRTP-Parameterverklaring / Parameter Descriptions:
╔════════════╦════════════════════════════════════════════════════════════╗
║ Parameter  ║ Beschrijving (NL) / Description (EN)                      ║
╠════════════╬════════════════════════════════════════════════════════════╣
║ Λ (Lambda) ║ Instroom / geboorteratio – inflow rate / birth rate       ║
║ α (alpha)  ║ Infectieratio bij contact – infection probability         ║
║ β (beta)   ║ Contactfrequentie – contact rate                          ║
║ r          ║ Detectieratio / testintensiteit – detection rate          ║
║ ξ (xi)     ║ Tracing efficiëntie – tracing effectiveness               ║
║ K          ║ Zorgcapaciteit – healthcare system capacity               ║
║ η (eta)    ║ Snelheid van ziek worden – progression to infectious      ║
║ φ (phi)    ║ Behandelingsimpact – treatment effect on transmission     ║
║ μ (mu)     ║ Natuurlijke sterfte – natural death rate                  ║
║ σ (sigma)  ║ Sterfte bij behandeling – death during treatment          ║
║ τ (tau)    ║ Behandelduur – duration of treatment                      ║
║ θ (theta)  ║ Terugval / herinfectie – relapse or reinfection rate      ║
║ δ (delta)  ║ Infectie-specifieke sterfte – infection-induced mortality ║
║ γ (gamma)  ║ Herstelkans – recovery rate                               ║
╚════════════╩════════════════════════════════════════════════════════════╝

�� R₀ computed = 1.54099
�� R₀ > 1 → Epidemic spread possible.

�� Original Sensitivity Index of R₀:
--------------------------------------------------------------------------
Parameter     | Sₚ       | Interpretation                
--------------------------------------------------------------------------
Λ            |   +1.00000 | �� Strong effect            
α            |   +1.00000 | �� Strong effect            
β            |   +0.01113 | �� Small effect             
r             |   -2.51999 | �� Strong effect            
ξ            |   +1.00000 | �� Strong effect            
K             |   +0.39200 | �� Small effect             
η            |   +0.00370 | �� Negligible               
φ            |   +2.51999 | �� Strong effect            
μ            |   -1.00127 | �� Strong effect            
σ            |   +0.38830 | �� Small effect             
τ            |   -1.38830 | �� Strong effect            
θ            |   -0.01048 | �� Small effect             
δ            |   -0.98930 | �� Strong effect            
γ            |   -0.01008 | �� Small effect             
--------------------------------------------------------------------------

�� Dynamic Policy Implications:
--------------------------------------------------------------------------
Parameter     | Direction  | Policy Advice                 
--------------------------------------------------------------------------
Λ            | ↑ Increase | ⚠️ Reduce to lower R₀   
α            | ↑ Increase | ⚠️ Reduce to lower R₀   
r             | ↓ Decrease | ✅ Increase to lower R₀    
ξ            | ↑ Increase | ⚠️ Reduce to lower R₀   
φ            | ↑ Increase | ⚠️ Reduce to lower R₀   
μ            | ↓ Decrease | ✅ Increase to lower R₀    
τ            | ↓ Decrease | ✅ Increase to lower R₀    
δ            | ↓ Decrease | ✅ Increase to lower R₀    
--------------------------------------------------------------------------

�� Compartmental Effect Matrix (Qualitative):
--------------------------------------------------------------------------
Parameter  |   S   E   I   R   P   T
--------------------------------------------------------------------------
Λ         |   +   .   .   .   .   .
α         |   -   +   +   .   .   .
β         |   .   -   +   .   .   .
r          |   .   .   .   .   .   +
ξ         |   .   .   .   .   +   .
K          |   .   .   .   .   .   +
η         |   .   -   -   .   .   .
φ         |   .   .   .   .   -   .
μ         |   -   -   -   -   -   -
σ         |   .   .   .   .   -   .
τ         |   .   .   .   .   -   .
θ         |   .   -   .   +   .   .
δ         |   .   .   -   .   .   .
γ         |   .   .   -   +   .   .
--------------------------------------------------------------------------