# The Base expansion
restart;
de := diff(rho(h), h) = c(h)*k/rho(h); 
dsolve(de, rho(h));
DiBa := solve(%[1], rho(h)):
DiBa := subs(Int(c(h), h) = ICh, DiBa):
E1 := subs(ICh = NH, DiBa) = Ba(phi):
_C1:=solve(E1,_C1);
Eq := h*Int(subs(k = 1, DiBa)^2/2, phi = -Pi..Pi) 
        = H*Int(Ba(phi)^2/2, phi = -Pi..Pi):
EI1 := ICh = value(solve(Eq, ICh));
EI2 := subs(int(Ba(phi)^2, phi = -Pi..Pi) = 2*So, EI1):
EI3 := subs(ICh = int(c(h), h), EI2);
EI4 := normal(diff(EI3, h));
NH := subs(h = H, rhs(EI4));
DiBa:=normal(subs(ICh = rhs(EI3), DiBa));
de_test := simplify(subs(rho(h) = DiBa, c(h) = rhs(EI4), de)):
evalb(%);
DiBa:=unapply(subsop(1 = collect(expand(op(1, DiBa)), [So, Pi, k]), 
                  DiBa),phi);
C_t := evalb(r*sqrt(n) = simplify(subs(So = Pi*r^2, 
               Ba(phi) = r, h = H/n, k = 1, DiBa(phi)),symbolic));
save(DiBa, `DiBa.sav`);


# Base Described by One and Several Functions
restart;
with(linalg): with(plots):
e1 := cx^2 + cy^2 = b^2:
e2 := (cx - b)^2 + cy^2 = a^2:
allvalues(solve({e1, e2}, {cx, cy}));
assign(%[1]);
A := [0, 0]:  B := [c, 0]:  C := [cx, cy]:
T := (A + B+ C)/3:
A := normal(A - T):  B := normal(B - T):  C := normal(C - T): 
Aa := abs(arctan(A[2]/A[1])):
Ab := abs(arctan(B[2]/B[1])):
Ac := Pi/2 - arctan(C[1]/C[2]):
rho := (Q[y]*P[x] - Q[x]*P[y])/
       ((P[x] - Q[x])*sin(phi) -(P[y] - Q[y])*cos(phi));
r1 := normal(subs(P[x] = A[1], Q[x] = B[1], P[y] = A[2], 
                  Q[y] = B[2], rho));
r2 := normal(subs(P[x]  =B[1], Q[x] = C[1], P[y] = B[2], 
                  Q[y] = C[2], rho));
r3 := normal(subs(P[x] = C[1], Q[x] = A[1], P[y] = C[2], 
                  Q[y] = A[2], rho));
o := (a + b + c)/2:
Sb := sqrt(o*(o - a)*(o - b)*(o - c)):
save(r1, r2, r3, Sb, Aa, Ab, Ac, `Tri.sav`);
read `DiBa.sav`:
Sus := a = 8, b = 7, c = 6:  Suf := h = N*H, So = Sb:
Aa := evalf(subs(Sus, Aa)):
Ab := evalf(subs(Sus, Ab)):
Ac := evalf(subs(Sus, Ac)):
R1 := evalf(subs(Suf, Ba(phi) = r1, Sus, DiBa(phi))):
R2 := evalf(subs(Suf, Ba(phi) = r2, Sus, DiBa(phi))):
R3 := evalf(subs(Suf, Ba(phi) = r3, Sus, DiBa(phi))):
RE := subs(h = N*H, Ba(phi) = a*b/sqrt(a^2*sin(phi)^2
      + b^2*cos(phi)^2), So = Pi*a*b, a = 3.5, b = 1.85, DiBa(phi)):
SqN := [seq(1 - i*0.1, i = 0..8)]:
with(plots):
for k from 1/3 by 1/3 to 1 do:
  P1 := polarplot({seq(R1, N = SqN)}, phi= - Pi + Aa..-Ab):
  P2 := polarplot({seq(R2, N = SqN)}, phi = -Ab..Ac):
  P3 := polarplot({seq(R3, N = SqN)}, phi = Ac..Pi + Aa):
  print(display({P1, P2 ,P3}, scaling = constrained, 
          view = [-6.5..6.5, -6..7]));
  print(polarplot({seq(RE, N = SqN)}, phi = -Pi..Pi, 
             scaling = constrained, view = [-6.5..6.5, -6..7]));
od:
k := 'k':


# The Lateral Side Distortion
restart;
LaSi := Q2*z^2 + Q1*z + Q0:
E1 := subs(z = 0, LaSi) = DiBa(phi);
assign(%);
E2 := subs(z = h, LaSi) = DiBa(phi):
Q1 := solve(E2, Q1);
E3 := Int(Int(LaSi^2/2, z = 0..h), phi = -Pi..Pi) = H*So;
I1 := normal(int(LaSi^2/2, z = 0..h));
I2:=map(u->simplify(int(u,phi=-Pi..Pi)),I1);
read(`DiBa.sav`):
I3 := subs(int(DiBa(phi)^2,phi = -Pi .. Pi) = -So*(k*h - k*H - h)*2/h, I2);
E3:=subs(select(has,select(has,I3,int),int)=IDiBa,I3)=H*So;
Sol := solve(E3, Q2):
T1 := Sol[1]: T2 := Sol[2]:
LaSi_1:=[subs(Q2 = T1, LaSi), subs(Q2 = T2, LaSi)]:
Cylinder_Subs := IDiBa = int(DiBa(phi), phi = -Pi..Pi), 
    So = Pi*r^2, Ba(phi) = r, r=1 , H = 4, k = kappa, h = 1:
LaSi_2 := simplify(subs(Cylinder_Subs, LaSi_1)):
Kappa := [seq(1 - i/2, i = 0..2)]:
with(plots):
P1 := plot({seq([LaSi_2[1], z, z = 0..1], kappa = Kappa)}, 
             thickness = 1):
kappa := 'kappa':
P2 := plot({seq([LaSi_2[2], z, z = 0..1], kappa = Kappa)}, 
             thickness = 3):
kappa := 'kappa':
display({P1, P2});
a1 := expand(subs(Q2 = T2, LaSi)):
a2:=factor(remove(has,a1,Ba(phi))):
a3:=select(has,select(has,a2,So),So)/sqrt(5):
a4:=collect(a3^2,[So,k,Pi]):
a5:=subs(24*h-24*H=B,-24*h+24*H=-B,a4):
a6:=sqrt(map(u->factor(u/h^2),a5)):
LaSi:=subs(a3=a6,B=24*h*(1-H/h),IDiBa=Int('DiBa(phi)',phi=A..B),a2*h+'DiBa(phi)');
save(LaSi, `LaSi.sav`):


# Non-centered Bases
restart;
Tx := T*cos(theta):  Ty := T*sin(theta):
r1x := rho(phi)*cos(phi) + Tx = r(psi)*cos(psi):
r1y := rho(phi)*sin(phi) + Ty = r(psi)*sin(psi):
r2x := map(t -> t - Tx, r1x): r2y := map(t -> t - Ty, r1y):
R1 := combine(map(t -> t^2, r2x) + map(t -> t^2, r2y)):
sol := rho(phi) = combine(solve(R1, rho(phi))[1]);
R2 := phi = solve(lhs(r2y)/lhs(r2x) = rhs(r2y)/rhs(r2x), phi):
dR2 := D(phi) = normal(diff(rhs(R2),psi)):
df := factor(combine(rhs(dR2)));
Cf := solve(r2x, cos(phi)):  Sf := solve(r2y, sin(phi)):
Cf := subs(sol, Cf);  Sf := subs(sol, Sf);
read(`LaSi.sav`):  read(`DiBa.sav`):
FRzS := subs(Int(DiBa, phi = A .. B) = 
     Int(DiBa*df, psi = A ..B), Ba(phi) = Ba(psi), LaSi):
with(codegen,optimize,makeproc,cost):
FO := codegen[optimize]([W1 = FRzS,  W2 = Cf, W3 = Sf]):
FI := codegen[optimize](op(1,rhs(op(1, select(has, [FO], Int))))):
EIco := NDigits:
EIco := NDigits, _NumMethod;
#?evalf[int]


# 3D Graphics / Centered Base
restart;
read(`Tri.sav`):
a := 8:  b := 7:  c := 6:
Tr := [unapply(r1, phi), unapply(r2, phi), unapply(r3, phi)];
 TrBo := [-Pi + Aa, -Ab, Ac, Pi + Aa];
So := simplify(Sb);
T := 0:  theta := 0:  n:= nops(Tr): EIco := 12:
read(`Body.map`):
with(plots):
for i from 0 to 3 do;
  h := 16 - 4*i;
  Body(Tr, TrBo, 16, h, 0.4);
  for j from 1 to n do:
    PL[j] := cylinderplot(Rho(u, z)[j], u = TrBo[j]..TrBo[j + 1],
             z = 0..h, grid=[3^i + 1, 3^i + 1]); 
  od:
  SPL[i] := display({seq(PL[j], j = 1..n)}, color = black,
            style = hidden):
od:
display({seq(SPL[i], i = 0..3)}, axes=boxed, 
         scaling = constrained);


# Non-centered, Segmented Base
restart;
EL := 18/sqrt(36*sin(u)^2 + 9*cos(u)^2):
p1 := -6/cos(u):
p2 := -4/sin(u):
k2 := solve(subs(x = r*cos(u), y = r*sin(u), 
           (x+1)^2 + y^2 = 16) ,r)[1]:
CC := [u -> 3, unapply(EL, u), unapply(p1, u), unapply(p2, u),
         unapply(k2, u)];
alpha[3] := Pi + arctan(4/6):
alpha[4] := 3/2*Pi - arctan(1/4):
CB := [0, Pi/2, Pi, alpha[3], alpha[4], 2*Pi];
save(CC,CB, `crazy.sav`);
EIco := 10, _CCquad:
n := nops(CC): 
read(`MsCe.map`):
MsCe(CC, CB, 0);
read(`Body.map`):
Body(CC, CB, 20, 8, 0.2);
with(plots):
for i from 1 to n do:
  pl[i] := plot3d([X(u, z)[i], Y(u, z)[i], z], 
           u = CB[i]..CB[i + 1], z = 0..8);
od:
display({seq(pl[i], i = 1..n)}, color = black,
        style = hidden, orientation = [-100, 15], axes = normal);


# Convex Polygon Base
restart;  
PNTi := [[1,1], [4,3], [3,4], [2,4], [0,2]]:
PNTo := [PNTi[], PNTi[1]]:
n := nops(PNTi):
Tx := sum(op(i, PNTi)[1], i=1..n)/n; 
Ty := sum(op(i, PNTi)[2], i= 1..n)/n;
readlib(polar):  read(`MsCe.map`):   EIco := 10, _CCquad:  
T := op(1, polar(Tx + I*Ty)):  theta := op(2, polar(Tx + I*Ty)):
while evalf(T) > 0 do;
  Tx := T*cos(theta):
  Ty := T*sin(theta):
  PNTo := [seq([op(i, PNTo)[1] - Tx, op(i, PNTo)[2] - Ty], 
            i = 1..n+1)];
  i := 'i';  j := i + 1;
  r := [seq(unapply(normal((PNTo[j, 1]*PNTo[i, 2]
        -PNTo[i, 1]*PNTo[j, 2])/
         (sin(u)*(PNTo[j, 1] - PNTo[i, 1]) + 
         cos(u)*(PNTo[i, 2]  - PNTo[j, 2]))), u), i = 1..n)];
  bo := [seq(op(2, polar(PNTo[i, 1] + I*PNTo[i,2 ])), 
          i = 1..n + 1)]:
  for i from 1 to n do; 
     if evalf(bo[j]) < evalf(bo[i]) then 
        bo[j] := bo[j] + 2*Pi; 
     fi; od:
  MsCe(r, bo, 0);
od:
read(`Body.map`):
Body(r, bo, 4, 2, 0.3):
with(plots):
for i from 1 to n do;
  pl[i] := cylinderplot(Rho(u, z)[i], u=bo[i]..bo[j], z = 0..2); 
od:
display({seq(pl[i], i = 1..n)});


# 3 D Animation
restart; with(plots):
Grid := [5, 5, 3, 3, 5]:
Frames := 20: Hi := 20: hf := 6: k := 1/3:
read(`crazy.sav`): read(`MsCe.map`):
n := nops(CC): EIco := 10,_Dexp:
MsCe(CC,CB,1):
read(`Body.map`):
for j from 0 to Frames do;
  h := Hi - (Hi - hf)*j/Frames;
  Body(CC, CB, Hi, h, k);
  An[j] := [seq(op(plot3d([X(u, z)[i], Y(u, z)[i], z], 
            u = CB[i]..CB[i+1], z = 0..h, grid = [Grid[i], 15])), 
            i = 1..n)]; 
  print(j);
od:
display(seq(display(An[i]),i=0..20),insequence=true,title=`Animation`,axes=boxed,scaling=constrained);
display(subs(TITLE(Animation)=NULL,%));