# Finding the Minimal Illumination Point on a Road
restart;
S[1] := P[1]*h[1]/(h[1]^2 + x^2)^(3/2):
S[2] := P[2]*h[2]/(h[2]^2 + (s - x)^2)^(3/2):
C :=  S[1] + S[2]:
dC := diff(C, x);
# solve(dC=0,x):
eq := diff(S[1], x)^2 - diff(S[2], x)^2; 
eq := collect(primpart(numer(eq)), x);
P[1] := 2000:  P[2] := 3000: s := 20: h[1] := 5: h[2] := 6: 
plot({C, S[1], S[2], dC}, x = -s/2..s*3/2);
Xm := fsolve(dC = 0, x, 5..10);
Cmin := subs(x = Xm, C);
Xe := fsolve(S[1] = S[2], x, 0..s);
dX := Xm - Xe;
sol:=[solve(eq, x)]:# 
solf:=select(t -type(t, numeric) and t 0 and t < s, evalf(sol));
map(t -> if subs(x = t, diff(dC, x)) < 0 then max else min fi, solf);
map(t -> subs(x = t, C), solf);


# Varying h2 to Maximize the Illumination
h[2] := 'h[2]':
H2 := array(0..30):  # array for the values of h[2] 
X  := array(0..30):  # array for the values of x(h[2])
for i from 0 to 30 do
   H2[i] := 3 + 6*i/30:
   X[i] := fsolve(subs(h[2]=H2[i], dC), x, s/4..3*s/4):
od:
H2 := convert(H2, list):
X  := convert(X, list):
plot(zip((h2, x) -[h2, x], H2, X));# 
f := unapply(C, x, h[2]):
Cu := [seq([X[i], H2[i], f(X[i], H2[i])], i = 1..31)]:
with(plots):
PL1 := spacecurve(Cu, thickness = 4):
PL2 := plot3d(subs(h[2] = h2 ,C), x = -s/2..3*s/2, h2 = 3..9, 
              style = wireframe):
display({PL1, PL2});# 
with(linalg):
g := grad(C, [x, h[2]]);
Sol := fsolve({g[1] = 0, g[2] = 0}, {x, h[2]}, 
              {x = 0..s, h[2] = 3..9});
H := subs(Sol, hessian(C, [x, h[2]]));
eigenvals(H); 
{solve(g[2] = 0, h[2])};
tan(alpha[2]) = normal(%[2]/(s-x));
evalf(arctan(rhs(%)));
evalf(convert(%, degrees));

# Optimal Illumination
P[1] := 'P[1]': P[2] := 'P[2]': h[1] := 'h[1]': 
h[2] := 'h[2]': s := 's': 
g := grad(C, [x, h[1], h[2]]);
g:=subs(s=sigma+x,convert(g,list));
sh1 := {solve(g[2] = 0, h[1])};
sh2 := expand({solve(g[3] = 0, h[2])});
ho[1]:=remove(has,sh1,-1/2)[]; #ho[1] := sh1[2];
ho[2]:=remove(has,sh2,-1/2)[]; #ho[2]:= sh2[1];
Q := simplify(subs(h[1] = ho[1], h[2] = ho[2], g[1]),symbolic);
q:=subs(P[1]=p[1]^3,P[2]=p[2]^3,select(has,numer(Q),P[1]))=0;# 
q:=map(u->simplify((u+p[2]^3*x^3)^(1/3),symbolic),q);# 
q:=subs(sigma=s-x,p[1]=P[1]^(1/3),p[2]=P[2]^(1/3),q);
Xo:=solve(q,x);
s := 20: 
plot3d(Xo, P[1] = 0..2000, P[2] = 0..3000, orientation = [110, 60],
axes=BOXED);# 
P[1] := 2000: P[2] := 3000: x := Xo: 
h[1] := evalf(ho[1]);
h[2]:=evalf(subs(sigma=s-x,ho[2]));
evalf( (Xm - Xo)/Xm );
evalf( (C - Cmin)/Cmin );
Cop := evalf(C);
x := 'x': Cop := C: h[1] := 5: h[2] := 6: 
plot({C, Cop}, x = -s/2..3*s/2);#