# ---------------------------------------------------------------- #

vect  := (A,B,n) -> [seq(B[i]-A[i],i=1..n)]: 
scal  := (U,V,n) -> add( U[i]*V[i], i=1..n):
dist2 := (A,B,n) -> add( (B[i]-A[i])^2, i=1..n);
dist  := (A,B,n) -> sqrt(dist2(A,B,n)): 


CP := proc(k,i,n) local Delta;
  Delta := dist2(A,B,n)*dist2(A,P,n)-(scal(vect(A,B,n),vect(A,P,n),n))^2;
  if k = 0 then
    if i = 1 then 
       ln(( dist(B,A,n)*dist(B,P,n)+scal(vect(B,A,n),vect(B,P,n),n))
          /(dist(A,B,n)*dist(A,P,n)-scal(vect(A,B,n),vect(A,P,n),n)));
    elif i = 2 then
       ( arctan( scal(vect(B,A,n),vect(B,P,n),n) /sqrt(Delta)) 
        + arctan(scal(vect(A,B,n),vect(A,P,n),n)/sqrt(Delta))
        ) *dist(A,B,n)/sqrt(Delta);
    else 
       dist(A,B,n)/(i-2)/Delta*(
         (i-3)*dist(A,B,n)*CP(0,i-2,n) 
         + scal(vect(B,A,n),vect(B,P,n),n)/dist(B,P,n)^(i-2)
         + scal(vect(A,B,n),vect(A,P,n),n)/dist(A,P,n)^(i-2)
          );
    end if;
  elif k = 1 then
    if i=2 then 
       scal(vect(A,B,n),vect(A,P,n),n)/dist2(A,B,n)*CP(0,2,n)
      + ln(dist(B,P,n)/dist(A,P,n))/dist(A,B,n);
    else
       scal(vect(A,B,n),vect(A,P,n),n)/dist2(A,B,n)*CP(0,i,n)
      + ((dist(B,P,n))^(2-i) - (dist(A,P,n))^(2-i))/dist(A,B,n)/(2-i);
    end if;
  elif k = i-1 then 
      scal(vect(A,B,n),vect(A,P,n),n)/dist2(A,B,n)*CP(i-2,i)
     + CP(i-3,i-2,n)/dist2(A,B,n) + (dist(B,P,n))^(2-i)/(2-i)/dist(A,B,n)
  else
      (i-2*k)/(i-k-1)*scal(vect(A,B,n),vect(A,P,n),n)/dist2(A,B,n)*CP(k-1,i,n)
       + (k-1)/(i-k-1)*dist2(A,P,n)/dist2(A,B,n)*CP(k-2,i,n)
       - (dist(B,P,n))^(2-i)/dist(A,B,n)/(i-k-1) ;       
  end if;
end proc:




# ---- Checking the procedure ------  #


I[k,i] = Int( t^k/(a*t^2-2*b*t+c)^(i/2), t);
 
II := proc(k,i,a,b,c) 
   if k=0 then   
        if i=1 then ln((a*t-b)/sqrt(a)+sqrt(a*t^2-2*b*t+c))/sqrt(a)
        elif i=2 then arctan((1/2)*(2*a*t-2*b)/sqrt(c*a-b^2))/sqrt(c*a-b^2)
        elif i>2 then (  a*(i-3)*II(0,i-2,a,b,c)/(i-2)/(c*a-b^2)  + (a*t-b)/(a*t^2-2*b*t+c)^(i/2-1)/(i-2)/(c*a-b^2))
        else ( i*(b^2-c*a)*II(0,i+2,a,b,c)+ (a*t-b)/(a*t^2-2*b*t+c)^(i/2) ) / a / (1-i) ;
        end if;
   elif k=1 then     
        if i=2 then ln( a*t^2-2*b*t+c)/2/a +b/a*II(0,2,a,b,c)
        else  b/a*II(0,i,a,b,c)+ 1/(a*(2-i))*1/(a*t^2-2*b*t+c)^(i/2-1)
        end if;
   elif i=k+1 then b*II(i-2,i,a,b,c)/a+II(i-3,i-2,a,b,c)/a+ t^(i-2)/(a*t^2-2*b*t+c)^(i/2-1)/a/(2-i)
   else  ( b*(i-2*k)*II(k-1,i,a,b,c) )/(a*(i-k-1))    + (c*(k-1)*II(k-2,i,a,b,c))/ (a*(i-k-1))  - (t^(k-1)/(a*t^2-2*b*t+c)^(i/2-1) )/ (a*(i-k-1)) 
   end if;
end proc:
     
for i from -1 to 5 do 
  for k from 0 to 5 do
II(k,i,a,b,c); normal(diff(%,t));
  printf ("i=%a, k=%a - %a\n", i, k, evalb(%= t^k/(a*t^2-2*b*t+c)^(i/2))); 
end do: end do:


F[k,i] =  Int( t^k/(a*t^2-2*b*t+c)^(i/2), t=0..1);
  
F := proc(k,i) local f;
    f := unapply( II(k,i,a,b,c), t);
    unapply( normal(f(1))-normal(f(0)), [a,b,c])
end proc:


compare := proc(k,i) local f,g,flop,glop:
      f := F(k,i);
    flop :=  dist([a1,a2],[b1,b2],2)*f(dist2([a1,a2],[b1,b2],2), 
          scal(vect([a1,a2],[b1,b2],2),vect([a1,a2],[x,y],2)),  dist2([a1,a2],[x,y],2)):
      g := unapply( CP(k,i), [A,B,P]): glop := g([a1,a2], [b1,b2],[x,y]):
     normal( simplify(simplify(glop-flop,ln, symbolic),symbolic));
end proc:

for i from 1 to 5 do
for k from 0 to 3 do 
   print(k,i, compare(k,i));
end do:
end do:

for i from 1 to 4 do
for k from 0 to 4 do 
   print(k,i, compare(k,i));
end do:
end do:

# ---Illustration ----------------------------------------------- #

Skel := [ [[0,0],[0,-0.7]], [[0,0],[0,1.5]], [[0,0],[3/4,1/2]], [[3/4,1/2],[1.,0.1]], [[0,0],[-3/4,1/2]],
      [[-3/4,1/2],[-1.,0.1]] , [[0.2,-0.2], [0.6,-0.6]],  [[-0.2,-0.2], [-0.6,-0.6]]]:

i:= 3:
p := t->b*t^2+c*t+d:
solve({p(0)=0.01,  p(1/3)=0.008, p(1)=0.2}, {a,b,c,d,f}):
p1 := unapply(expand(subs(%, p(t))),t);plot(p1(s), s=0..1);
f1 := unapply(unapply( add( coeff(p1(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[1]),[x,y]), [x,y]):
p := t->f*t^4+a*t^3+b*t^2+c*t+d:
solve({p(0)=0.01,  D(p)(0)=0, p(1/3)=0.03, p(2/3)=0.2,p(1)=0}, {a,b,c,d,f}):
p2 := unapply(expand(subs(%, p(t))),t);plot(p2(s), s=0..1);
f2:= unapply(unapply( add( coeff(p2(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[2]),[x,y]), [x,y]):
p := t->b*t^2+c*t+d:
solve({p(0)=0.01, D(p)(0)=0, p(1)=0.1}, {a,b,c,d,f}):
p3 := unapply(expand(subs(%, p(t))),t);plot(p3(s), s=0..1);
P5 := p3:
f3:= unapply(unapply( add( coeff(p3(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[3]),[x,y]), [x,y]):
f5:= unapply(unapply( add( coeff(p3(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[5]),[x,y]), [x,y]):
p := t->c*t+d:
solve({p(0)=0.1,  p(1)=0.005}, {a,b,c,d,f}):
p4 := unapply(expand(subs(%, p(t))),t);
p6 := p4: plot(p4(s), s=0..1);
f4:= unapply(unapply( add( coeff(p4(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[4]),[x,y]), [x,y]):
f6:= unapply(unapply( add( coeff(p4(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[6]),[x,y]), [x,y]):
p := t->a*t^3+b*t^2+c*t+d:
solve({p(0)=0.001, D(p)(0)=0.1,  p(2/3)=0.05,p(1)=0.2}, {a,b,c,d,f}):
p7 := unapply(expand(subs(%, p(t))),t);plot(p7(s), s=0..1);
p8 := p7:
f7:= unapply(unapply( add( coeff(p7(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[7]),[x,y]), [x,y]):
f8:= unapply(unapply( add( coeff(p7(t),t,k)*CP(k,i),k=0..degree(p(t),t)), [A,B,P])(op(Skel[8]),[x,y]), [x,y]):

# for i=3
implicitplot( [seq(f1(x,y)+f2(x,y)+f3(x,y)+f4(x,y)+f5(x,y)+f6(x,y) +f7(x,y)+f8(x,y)= c, c=[2.5,5,15])], x=-1.5..1.5,y=-1..2, grid=[250,250],scaling=constrained, color=[ COLOR(RGB, 0.26,0.9,0.1),COLOR(RGB, 0.26,0.65,0.95),COLOR(RGB, 0.65,0.25,0.95),COLOR(RGB, 0.95,0.25,0.9), COLOR(RGB, 0.9,0.35,0.65), COLOR(RGB, 0.9,0.35,0.15) , COLOR(RGB, 0.65,0.35,0.15)],view=[-1.2..1.2,-0.9..1.7], axes=none);

# for i=5
implicitplot( [seq(f1(x,y)+f2(x,y)+f3(x,y)+f4(x,y)+f5(x,y)+f6(x,y) +f7(x,y)+f8(x,y)= c, c=[50,150,350])], x=-1.5..1.5,y=-1..2, grid=[250,250],scaling=constrained, color=[ COLOR(RGB, 0.26,0.9,0.1),COLOR(RGB, 0.26,0.65,0.95),COLOR(RGB, 0.65,0.25,0.95),COLOR(RGB, 0.95,0.25,0.9), COLOR(RGB, 0.9,0.35,0.65), COLOR(RGB, 0.9,0.35,0.15) , COLOR(RGB, 0.65,0.35,0.15)],view=[-1.2..1.2,-0.9..1.7],axes=none);

perp := proc(A,B) local C; 
   C := [B[2]-A[2], -B[1]+A[1]]:
   [C[1]/sqrt(C[1]^2+C[2]^2),C[2]/sqrt(C[1]^2+C[2]^2)];
end proc:
pg := polygonplot(Skel, thickness=3,color=blue):
for i from 1 to nops(Skel) do 
  Skel[i][1]+ map(`*`,Skel[i][2]-Skel[i][1],t)+ map(`*`,perp(op(Skel[i])),p||i(t)); 
d1 := plot([op(%), t=0..1]):
Skel[i][1]+ map(`*`,Skel[i][2]-Skel[i][1],t)+ map(`*`,perp(op(Skel[i])),-p||i(t)); 
d2 := plot([op(%), t=0..1]):
g||i := display({d1,d2}):
end do:
display({pg, seq(g||i,i=1..nops(Skel))},axes=none);