kill(all)$
ratprint:false$
load("vect")$
load("eigen")$
nnorm(x):=sqrt(x.x)$
assume(h1>0)$ assume(h2>0)$ assume(b1>0)$ assume(b2>0)$
hm:(h1*b2+h2*b1)/(b1+b2)$
l1:sqrt(b1**2+(h1-hm)**2)$
l2:sqrt(b2**2+(hm-h2)**2)$

A:factor(2*h2*t/6 + l2*t/6*(2+hm/h2))$
B:factor(l2*t/6*(2+h2/hm) + l1*t/6*(2+h1/hm) + 2*hm*t/6)$
C:factor(2*h1*t/6 + l1*t/6*(2+hm/h1))$
D:C$
E:B$
F:A$
Ars:[A,B,C,D,E,F]$
nb:length(Ars)$

coords:[[b2,h2],[0,hm],[-b1,h1],[-b1,-h1],[0,-hm],[b2,-h2]]$
xs:create_list(coords[i][1], i, 1, nb)$
ys:create_list(coords[i][2], i, 1, nb)$

cents:create_list(factor((coords[mod(i,nb)+1]+coords[i])/2), i, 1, nb)$
svects:create_list(factor(uvect(coords[mod(i,nb)+1]-coords[i])), i, 1, nb)$
/* Adding extra element for panel joining B and E */
cents:append(cents,[[0,0]])$
svects:append(svects, [factor(uvect(coords[5]-coords[2]))])$

/* Length of Panel */
ws:create_list(factor(nnorm(coords[mod(i,nb)+1]-coords[i])), i, 1, nb)$
ws:append(ws, [factor(nnorm(coords[5]-coords[2]))])$

## Idealized Analysis

/* Perpendicular Distance of centroid */
ps:create_list(factor(express(cents[i] ~ svects[i])), i, 1, nb+1)$
/* Overall Properties */
Atot:factor(lsum(a, a, Ars))$
I22e:factor(lsum(y2a, y2a, ys**2*Ars))$
xC:factor(lsum(xA, xA, xs*Ars)/Atot)$
xCtrue:factor(lsum(cw, cw, cents*ws)/lsum(w,w,ws))$
/* Cell Areas */
A1:factor(ps[1]*ws[1]+ps[7]*ws[7]+ps[5]*ws[5]+ps[6]*ws[6])/2$
A2:factor(ps[2]*ws[2]+ps[3]*ws[3]+ps[4]*ws[4]-ps[7]*ws[7])/2$

qCD:qBC - factor(V3/I22*ys[3]*Ars[3])$
qDE:qCD - factor(V3/I22*ys[4]*Ars[4])$
qBE:qAB-qBC - factor(V3/I22*ys[2]*Ars[2])$
qEF:qDE+qBE - factor(V3/I22*ys[5]*Ars[5])$
qFA:qEF - factor(V3/I22*ys[6]*Ars[6])$
qs:[qAB,qBC,qCD,qDE,qEF,qFA,qBE]$

/* Twist Compatibility */ 
θ1:factor((qAB*ws[1]+qBE*ws[7]+qEF*ws[5]+qFA*ws[6])/A1)/2$
θ2:factor((qBC*ws[2]+qCD*ws[3]+qDE*ws[4]-qBE*ws[7])/A2)/2$
qBCsol:factor(solve(θ1=θ2, qBC))$

/* Moment Equality */
Mom:factor(subst(qBCsol, lsum(qpw, qpw, qs*ps*ws)))$
qABsol:factor(solve(Mom-V3*(b1-bL), qAB))$

/* Numerical Shear Flow */
pars:[h1=200,h2=100,b1=100,b2=100, bL=1.0, V3=1]$
qABn:subst(pars,qABsol)$
qBCn:subst(append(pars,qABn), qBCsol)$
q0ns:append(qABn,qBCn)$
qns:subst(append(pars,q0ns), qs);

/* Shear Center */
q00ns:solve([θ1,θ2], [qAB,qBC])$
ξid:factor(subst(I22=I22e, subst(q00ns, lsum(qpw, qpw, qs*ps*ws))/V3))$

## Exact Analysis

I22fun(p0,p1):=block(
    [ell,xy],
    ell:nnorm(p1-p0),
    xy:p0*(1-s/ell) + p1*s/ell,
    integrate(t*xy[2]**2, s, 0, ell)
    )$

I22exact:factor(sum(I22fun(coords[i], coords[mod(i,nb)+1]), i, 1, nb) + I22fun(coords[2], coords[5]))$

qfun(p0,p1,q0):=block(
    [ell,p,xy,q,qintS,pqintS,qend],
    ell:nnorm(p1-p0),
    xy:p0*(1-s/ell) + p1*s/ell,
    p:factor(express(xy ~ uvect(p1-p0))),
    q:factor(q0 - t*V3/I22*integrate(xy[2], s, 0, s)),
    qintS:factor(integrate(q, s, 0, ell)),
    pqintS:factor(integrate(factor(p*q), s, 0, ell)),
    qend:factor(subst([s=ell], q)),
    [q, qend, qintS, pqintS]
    )$

/* Left Cell (Cell 1) */
qBCall:qfun(coords[2], coords[3], q0B1)$
q0C:qBCall[2]$
qCDall:qfun(coords[3], coords[4], q0C)$
q0D:qCDall[2]$
qDEall:qfun(coords[4], coords[5], q0D)$
q0E1:qDEall[2];

/* Right Cell (Cell 2) */
qABall:qfun(coords[1], coords[2], q0A)$
q0B2:qABall[2]$
qBEall:qfun(coords[2], coords[5], q0B2-q0B1)$
q0Em:qBEall[2]$
qEFall:qfun(coords[5], coords[6], q0E1+q0Em)$
q0F:qEFall[2]$
qFAall:qfun(coords[6], coords[1], q0F)$
q0A1:qFAall[2]$

/* Twist Rates */
θ1ex:factor((qBCall[3]+qCDall[3]+qDEall[3]-qBEall[3])/2/A1)$
θ2ex:factor((qABall[3]+qBEall[3]+qEFall[3]+qFAall[3])/2/A2)$
q0AB1:solve([θ1ex=0,θ2ex=0], [q0A,q0B1])[1]$

/* Twisting Moment */
Mtot:factor(subst(q0AB1, qBCall[4]+qCDall[4]+qDEall[4]+qBEall[4]+qABall[4]+qEFall[4]+qFAall[4]))$
ξex:factor(subst(I22=I22exact, Mtot/V3))$

## Visualization

wxanimate_autoplay:false$
parsh1:[h1=h1h,h2=10,b1=100,b2=100, bL=50]$
hmin:0.1*100$ hmax:3*100$ nh:10$
with_slider(
    h1h,
    makelist(hmin + (hmax-hmin)*i/nh, i, 0, nh), 
    subst([h1=h1h,h2=100,b1=100,b2=100, bL=50], 
        [
            [discrete, map(lambda([i], coords[i]), [1,2,3,4,5,6,1,2,5])], 
            [discrete, coords],
            [discrete, [[xCtrue[1],0]]],
            [discrete, [[xC,0]]],
            [discrete, [[ξex,0]]],
            [discrete, [[ξid,0]]]
        ]),
    [style, lines, points, points, points, points, points],
    [legend, "Section", "Booms", "True Centroid", "Idealized Centroid", "True Shear Center", "Idealized Shear Center"],
    [y, -hmax, hmax], [x, -hmax/2, hmax/2], [gnuplot_preamble, "set key outside right center"],
    [same_xy]
), wxplot_size=[700,400]$

parsh1:[h1=h1h,h2=100,b1=100,b2=100, bL=1.0]$
hmin:0.1*100$ hmax:25*100$ nh:100$
h1s:makelist(hmin+(hmax-hmin)*i/nh, i, 0,nh)$
ξis:subst(h1h=h1s, subst(parsh1, ξid))$
ξes:subst(h1h=h1s, subst(parsh1, ξex))$
xCs:subst(h1h=h1s, subst(parsh1,xC))$
xCts:subst(h1h=h1s, subst(parsh1,xCtrue[1]))$

wxplot2d([
        [discrete, makelist([h1s[i], xCts[i]], i, 1, nh+1)],
        [discrete, makelist([h1s[i], xCs[i]], i, 1, nh+1)], 
        [discrete, makelist([h1s[i], ξes[i]], i, 1, nh+1)],
        [discrete, makelist([h1s[i], ξis[i]], i, 1, nh+1)]
    ], 
    [style, [lines,2,0], [lines, 2,1], [lines, 2,2], [lines, 2,3]],
    [legend, "True Centroid", "Idealized Centroid", "Exact Shear Center", "Idealized Shear Center"],
    [xlabel, "Length h_1"], [ylabel, "X_2 Coordinate"], [gnuplot_preamble, "set key left top"]), wxplot_size=[500,400];

parsh1;

float(subst(h1h=200, subst(parsh1, ξex)));