WRITEFILE("quadvel1.txt");

/* Roots of Legendre polynomial of specified order */
lr2(i):=
    block([],
        if i = 0 then return(-1/sqrt(3)),
    return(1/sqrt(3))
    )$

lr3(i):=
    block([],
        if i = 0 then return(0),
        if i = 1 then return(-sqrt(3/5)),
        return(sqrt(3/5))
    )$

lr4(i):=
    block([],
        if i = 0 then return(-sqrt((15-sqrt(120))/35)),
        if i = 1 then return( sqrt((15-sqrt(120))/35)),
        if i = 2 then return(-sqrt((15+sqrt(120))/35)),
        return(sqrt((15+sqrt(120))/35))
    )$

/* Associated weights for Gauss Legendre quadrature */
glw2(i):=1;

glw3(i):=
    block([],
        if i = 0 then return(8/9),
        return(5/9)
    )$

glw4(i):=
    block([],
        if i = 0 or i = 1 then return(0.6521451549),
        return(0.3478548451)
    )$


/* Wrappers for supported polynomials */
legendre_root(order,i):=
    block([],
        if order = 2 then return(lr2(i)),
        if order = 3 then return(lr3(i)),
        if order = 4 then return(lr4(i))
    )$

gauss_legendre_weight(order,i):=
    block([],
        if order = 2 then return(glw2(i)),
        if order = 3 then return(glw3(i)),
        if order = 4 then return(glw4(i))
    )$

/* Function to perform the numerical integration */
quad(order,func):=
    block([sum, x, w],
        sum:0,
        for i:0 thru order-1 do (
            x : legendre_root(order,i),
            w : gauss_legendre_weight(order,i),
            sum : sum + (w * func(x))),
        return(sum)
    )$


phi1(XI):=1/2*(1-XI);
phi2(XI):=1/2*(1+XI);

AA(t):=alpha(t)^2*A0A;
dAAdt(t):=diff(AA(t),t);

A(xi,t):=AA(t)+(AB-AA(t))*(1/lAtoB)*(s0+((le/2)*(1+xi)));
dAdt(xi,t):=dAAdt(t)*(1-((1/lAtoB)*(s0+(le/2)*(1+xi))));

/* integrands */
Ai11(xi):=phi1(xi)*V1*diff(phi1(xi)*A(xi,t),xi);

/* integrate numerically */
A11q:quad(4,lambda([xi],Ai11(xi)));

/* Test it out
yi1(xi,t):=phi1(xi) * (le/2) * dAdt(xi,t);
yi2(xi,t):=phi2(xi) * (le/2) * dAdt(xi,t);

Ai11(xi):=phi1(xi)*V1*diff(phi1(xi)*A(xi,t),xi);
Ai12(xi):=phi1(xi)*V2*diff(phi2(xi)*A(xi,t),xi);
Ai21(xi):=phi2(xi)*V1*diff(phi1(xi)*A(xi,t),xi);
Ai22(xi):=phi2(xi)*V2*diff(phi2(xi)*A(xi,t),xi);

y1q:ratsimp(quad(3,lambda([xi],yi1(xi,t))),lAtoB,A0A,AB);
y1e:ratsimp(integrate(yi1(xi,t), xi,-1,1),lAtoB,A0A,AB);

A11q:ratsimp(quad(4,lambda([xi],Ai11(xi))),V1,lAtoB,alpha(t)^2,AB);
A11e:ratsimp(integrate(Ai11(xi), xi,-1,1),V1,lAtoB,alpha(t)^2,AB);

tex(y1q);
tex(y1e);

tex(A11q);
tex(A11e);
*/

CLOSEFILE();