bodyx = -pi*E*((3*(ni - 1/2)*(-(4*cos(X*pi/Lx)^6*sin(Y*pi/(2*(Ly/2)))^7*a^4*(Ly/2)^6*pi^4)/3 - 8*cos(X*pi/Lx)^5*sin(Y*pi/(2*(Ly/2)))^6*a^3*Lx*(Ly/2)^6*pi^3 + Lx^2*pi^2*cos(X*pi/Lx)^4*(a^2*pi^2*(cos(X*pi/Lx)^2 + 1)*cos(Y*pi/(2*(Ly/2)))^2 + (2*a^2*pi^2*cos(X*pi/Lx)^2)/3 - 16*(Ly/2)^2)*((Ly/2)/2)^4*a^2*sin(Y*pi/(2*(Ly/2)))^5 + (16*Lx^3*pi*cos(X*pi/Lx)^3*(Ly/2)^4*(a^2*pi^2*(cos(X*pi/Lx)^2 + 3/4)*cos(Y*pi/(2*(Ly/2)))^2 + (3*a^2*pi^2*cos(X*pi/Lx)^2)/4 - 2*(Ly/2)^2)*a*sin(Y*pi/(2*(Ly/2)))^4)/3 + Lx^4*pi^2*cos(X*pi/Lx)^2*(sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/6 + (5*(Ly/2)^2)/3)*cos(Y*pi/(2*(Ly/2)))^4 + 104*(cos(X*pi/Lx)^2 + 10/13)*(Ly/2)^4*cos(Y*pi/(2*(Ly/2)))^2/3 + (104*cos(X*pi/Lx)^2*(Ly/2)^4)/3)*a^2*sin(Y*pi/(2*(Ly/2)))^3/4 + (7*Lx^5*pi*cos(X*pi/Lx)*(sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/14 + (5*(Ly/2)^2)/7)*cos(Y*pi/(2*(Ly/2)))^4 + (32*(Ly/2)^4*(cos(X*pi/Lx)^2 + 1)*cos(Y*pi/(2*(Ly/2)))^2)/7 + (48*cos(X*pi/Lx)^2*(Ly/2)^4)/7)*a*sin(Y*pi/(2*(Ly/2)))^2)/6 + (4*Lx^6*(sin(X*pi/Lx)^6*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4/64 + sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/32 + (Ly/2)^2/2)*cos(Y*pi/(2*(Ly/2)))^4 + 2*cos(Y*pi/(2*(Ly/2)))^2*(Ly/2)^4 + 2*cos(X*pi/Lx)^2*(Ly/2)^4)*sin(Y*pi/(2*(Ly/2))))/3 + (2*cos(X*pi/Lx)*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^4*a*Lx^7*(Ly/2)^2*pi)/3)*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2)/32 - pi*sin(Y*pi/(2*(Ly/2)))*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*a^2*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 16*(Ly/2)^2)/4)^2/4)*log((a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 8*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 8*Lx^2*(Ly/2)^2 - pi*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/(Lx^2*(Ly/2)^2)) + (-3*(ni - 1/2)*(-(4*cos(X*pi/Lx)^6*sin(Y*pi/(2*(Ly/2)))^7*a^4*(Ly/2)^6*pi^4)/3 - 8*cos(X*pi/Lx)^5*sin(Y*pi/(2*(Ly/2)))^6*a^3*Lx*(Ly/2)^6*pi^3 + Lx^2*pi^2*cos(X*pi/Lx)^4*(a^2*pi^2*(cos(X*pi/Lx)^2 + 1)*cos(Y*pi/(2*(Ly/2)))^2 + (2*a^2*pi^2*cos(X*pi/Lx)^2)/3 - 16*(Ly/2)^2)*(Ly/2)^4*a^2*sin(Y*pi/(2*(Ly/2)))^5 + (16*Lx^3*pi*cos(X*pi/Lx)^3*(Ly/2)^4*(a^2*pi^2*(cos(X*pi/Lx)^2 + 3/4)*cos(Y*pi/(2*(Ly/2)))^2 + (3*a^2*pi^2*cos(X*pi/Lx)^2)/4 - 2*(Ly/2)^2)*a*sin(Y*pi/(2*(Ly/2)))^4)/3 + Lx^4*pi^2*cos(X*pi/Lx)^2*(sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/6 + (5*(Ly/2)^2)/3)*cos(Y*pi/(2*(Ly/2)))^4 + 104*(cos(X*pi/Lx)^2 + 10/13)*(Ly/2)^4*cos(Y*pi/(2*(Ly/2)))^2/3 + (104*cos(X*pi/Lx)^2*(Ly/2)^4)/3)*a^2*sin(Y*pi/(2*(Ly/2)))^3/4 + (7*Lx^5*pi*cos(X*pi/Lx)*(sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/14 + (5*(Ly/2)^2)/7)*cos(Y*pi/(2*(Ly/2)))^4 + (32*(Ly/2)^4*(cos(X*pi/Lx)^2 + 1)*cos(Y*pi/(2*(Ly/2)))^2)/7 + (48*cos(X*pi/Lx)^2*(Ly/2)^4)/7)*a*sin(Y*pi/(2*(Ly/2)))^2)/6 + (4*Lx^6*(sin(X*pi/Lx)^6*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4/64 + sin(X*pi/Lx)^2*a^2*pi^2*((Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/32 + (Ly/2)^2/2)*cos(Y*pi/(2*(Ly/2)))^4 + 2*cos(Y*pi/(2*(Ly/2)))^2*(Ly/2)^4 + 2*cos(X*pi/Lx)^2*(Ly/2)^4)*sin(Y*pi/(2*(Ly/2))))/3 + (2*cos(X*pi/Lx)*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^4*a*Lx^7*(Ly/2)^2*pi)/3)*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2)/32 - pi*sin(Y*pi/(2*(Ly/2)))*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*a^2*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 16*(Ly/2)^2)/4)^2/4)*log((a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 8*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 8*Lx^2*(Ly/2)^2 + pi*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/(Lx^2*(Ly/2)^2)) + pi*(-pi^4*cos(X*pi/Lx)^8*(-1 - (3*log(2))/2 + ni)*(Ly/2)^8*a^4*sin(Y*pi/(2*(Ly/2)))^9 - 8*Lx*pi^3*cos(X*pi/Lx)^7*(-1 - (3*log(2))/2 + ni)*(Ly/2)^8*a^3*sin(Y*pi/(2*(Ly/2)))^8 + Lx^2*pi^2*cos(X*pi/Lx)^6*(pi^2*((-1 + ni)*cos(X*pi/Lx)^2 + 1/2 + (3*sin(X*pi/Lx)^2*log(2))/2)*a^2*cos(Y*pi/(2*(Ly/2)))^2 - 24*(-1 - (3*log(2))/2 + ni)*(Ly/2)^2)*(Ly/2)^6*a^2*sin(Y*pi/(2*(Ly/2)))^7 + 7*Lx^3*pi*cos(X*pi/Lx)^5*(Ly/2)^6*(pi^2*((ni - 13/14)*cos(X*pi/Lx)^2 + (9*sin(X*pi/Lx)^2*log(2))/7 + 3/7)*a^2*cos(Y*pi/(2*(Ly/2)))^2 - 32*(-1 - (3*log(2))/2 + ni)*(Ly/2)^2/7)*a*sin(Y*pi/(2*(Ly/2)))^6 + (3*Lx^4*cos(X*pi/Lx)^4*(pi^4*((Ly/2)^2*(-1 + ni)*cos(X*pi/Lx)^2 + ((3*(Ly/2)^2*log(2))/2 + Lx^2/6)*sin(X*pi/Lx)^2 + (Ly/2)^2*(1 + ni)/3)*sin(X*pi/Lx)^2*a^4*cos(Y*pi/(2*(Ly/2)))^4 - (2*pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 + (-76*ni + 66)*(Ly/2)^2)*cos(X*pi/Lx)^2 - 84*(sin(X*pi/Lx)^2*log(2) - ni/21 + 5/14)*(Ly/2)^2)*(Ly/2)^2*a^2*cos(Y*pi/(2*(Ly/2)))^2)/3 - 128*(-1 - (3*log(2))/2 + ni)*(Ly/2)^6/3)*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^5)/8 + (9*Lx^5*pi*cos(X*pi/Lx)^3*(Ly/2)^4*cos(Y*pi/(2*(Ly/2)))^2*(pi^2*((ni - 5/6)*cos(X*pi/Lx)^2 + sin(X*pi/Lx)^2*log(2) + (2*ni)/9 + 2/9)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^2 + (-(2*a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2)/3 + (32*(Ly/2)^2*(ni - 5/6))/3)*cos(X*pi/Lx)^2 + 32*(sin(X*pi/Lx)^2*log(2) - ni/6 + 5/12)*(Ly/2)^2/3)*a*sin(Y*pi/(2*(Ly/2)))^4)/4 + Lx^6*(sin(X*pi/Lx)^4*a^4*pi^4*(cos(X*pi/Lx)^2*ni + (3*sin(X*pi/Lx)^2*log(2))/2 + 1/2)*cos(Y*pi/(2*(Ly/2)))^4 - 2*pi^2*sin(X*pi/Lx)^2*((sin(X*pi/Lx)^2*pi^2*a^2*ni + (-42*ni + 31)*(Ly/2)^2)*cos(X*pi/Lx)^2 - 30*(sin(X*pi/Lx)^2*log(2) + (4*ni)/15 + 1/5)*(Ly/2)^2)*a^2*cos(Y*pi/(2*(Ly/2)))^2 - 52*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (48*(Ly/2)^2*(ni - 5/6))/13)*cos(X*pi/Lx)^2 - 48*(sin(X*pi/Lx)^2*log(2) - ni/3 + 1/2)*(Ly/2)^2/13)*(Ly/2)^2)*cos(X*pi/Lx)^2*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2*sin(Y*pi/(2*(Ly/2)))^3/16 + (5*Lx^7*pi*cos(X*pi/Lx)*(pi^2*((ni - 7/10)*cos(X*pi/Lx)^2 + (3*sin(X*pi/Lx)^2*log(2))/5 + 1/5)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 + ((-(8*a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2)/5 + (96*(Ly/2)^2*(ni - 2/3))/5)*cos(X*pi/Lx)^2 + 48*(sin(X*pi/Lx)^2*log(2) + ni/3 + 1/6)*(Ly/2)^2/5)*cos(Y*pi/(2*(Ly/2)))^2 - (48*cos(X*pi/Lx)^2*(ni - 1/2)*(Ly/2)^2)/5)*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*a*sin(Y*pi/(2*(Ly/2)))^2)/16 + Lx^8*(pi^4*((-1 + ni)*cos(X*pi/Lx)^2 + (3*sin(X*pi/Lx)^2*log(2))/2 - ni + 1)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^6 - 4*pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 + (-36*ni + 22)*(Ly/2)^2)*cos(X*pi/Lx)^2 + (-12*(Ly/2)^2*log(2) + Lx^2*(ni - 1/2))*sin(X*pi/Lx)^2 + 4*(ni - 3/2)*(Ly/2)^2)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 - 224*(Ly/2)^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (24*(Ly/2)^2*(ni - 2/3))/7)*cos(X*pi/Lx)^2 + (-(12*(Ly/2)^2*log(2))/7 + (2*Lx^2*(ni - 1/2))/7)*sin(X*pi/Lx)^2 - (4*(Ly/2)^2)/7)*cos(Y*pi/(2*(Ly/2)))^2 - 256*cos(X*pi/Lx)^2*(ni - 1/2)*(Ly/2)^4)*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*sin(Y*pi/(2*(Ly/2)))/256 + Lx^9*pi*cos(X*pi/Lx)*(sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^4*a^2*pi^2 + (-2*sin(X*pi/Lx)^2*a^2*pi^2 + 32*(Ly/2)^2)*cos(Y*pi/(2*(Ly/2)))^2 - 48*(Ly/2)^2)*(ni - 1/2)*cos(Y*pi/(2*(Ly/2)))^4*sin(X*pi/Lx)^4*a/64)*a^2)*sin(X*pi/Lx)/(2*rho*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 16*(Ly/2)^2)/4)^2*(a*sin(Y*pi/(2*(Ly/2)))*pi*cos(X*pi/Lx) + Lx)^2*(1 + ni)*(ni - 1/2)*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*a);

bodyy = -pi*E*cos(Y*pi/(2*(Ly/2)))*((-((ni - 1/2)*(cos(X*pi/Lx)^6*sin(Y*pi/(2*(Ly/2)))^7*a^5*(Ly/2)^6*pi^5 + cos(X*pi/Lx)^5*a^4*Lx*(Ly/2)^6*pi^4*(cos(X*pi/Lx)^2 + 6)*sin(Y*pi/(2*(Ly/2)))^6 + 2*Lx^2*pi^3*cos(X*pi/Lx)^4*(Ly/2)^4*((3*a^2*pi^2*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2)/8 + (Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/2 + 9*(Ly/2)^2)*a^3*sin(Y*pi/(2*(Ly/2)))^5 + (3*Lx^3*pi^2*cos(X*pi/Lx)^3*(sin(X*pi/Lx)^2*a^2*pi^2*(cos(X*pi/Lx)^2 + 4)*cos(Y*pi/(2*(Ly/2)))^2 - (8*(Ly/2)^2*cos(X*pi/Lx)^2)/3 + (128*(Ly/2)^2)/3)*(Ly/2)^4*a^2*sin(Y*pi/(2*(Ly/2)))^4)/4 - 4*Lx^4*pi*(-(3*sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^4*a^4*pi^4)/64 - sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2*a^2*(Ly/2)^2*pi^2 + (Ly/2)^4*(cos(X*pi/Lx)^2 - 7))*cos(X*pi/Lx)^2*(Ly/2)^2*a*sin(Y*pi/(2*(Ly/2)))^3 + (3*Lx^5*(sin(X*pi/Lx)^4*a^4*pi^4*(cos(X*pi/Lx)^2 + 2)*cos(Y*pi/(2*(Ly/2)))^4 - (16*sin(X*pi/Lx)^2*a^2*(Ly/2)^2*pi^2*(cos(X*pi/Lx)^2 - 2)*cos(Y*pi/(2*(Ly/2)))^2)/3 - (64*(Ly/2)^2*(cos(X*pi/Lx)^2*sin(X*pi/Lx)^2*a^2*pi^2 - 2*(Ly/2)^2))/3)*cos(X*pi/Lx)*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)/16 + Lx^6*pi*(sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4/8 + sin(X*pi/Lx)^2*a^2*(Ly/2)^2*pi^2*(cos(X*pi/Lx)^2 + 3)*cos(Y*pi/(2*(Ly/2)))^4 + 2*((sin(X*pi/Lx)^2*a^2*pi^2 - 8*(Ly/2)^2)*cos(X*pi/Lx)^2 + sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2)*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2 - 40*cos(X*pi/Lx)^2*(Ly/2)^4)*sin(X*pi/Lx)^2*a*sin(Y*pi/(2*(Ly/2)))/8 + cos(X*pi/Lx)*sin(X*pi/Lx)^2*Lx^7*(sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4 - 8*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^4*a^2*(Ly/2)^2*pi^2 + 32*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2*a^2*(Ly/2)^2*pi^2 - 128*(Ly/2)^4)/64)*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/2 + pi*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*(sin(Y*pi/(2*(Ly/2)))*a*pi + cos(X*pi/Lx)*Lx)*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2 + 16*(Ly/2)^2)/4)^2*a^2)*log((a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 8*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 8*Lx^2*(Ly/2)^2 - pi*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/(Lx^2*(Ly/2)^2)) + (((ni - 1/2)*(cos(X*pi/Lx)^6*sin(Y*pi/(2*(Ly/2)))^7*a^5*(Ly/2)^6*pi^5 + cos(X*pi/Lx)^5*a^4*Lx*(Ly/2)^6*pi^4*(cos(X*pi/Lx)^2 + 6)*sin(Y*pi/(2*(Ly/2)))^6 + 2*Lx^2*pi^3*cos(X*pi/Lx)^4*(Ly/2)^4*((3*a^2*pi^2*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2)/8 + (Ly/2)^2*cos(X*pi/Lx)^2 - sin(X*pi/Lx)^2*Lx^2/2 + 9*(Ly/2)^2)*a^3*sin(Y*pi/(2*(Ly/2)))^5 + (3*Lx^3*pi^2*cos(X*pi/Lx)^3*(sin(X*pi/Lx)^2*a^2*pi^2*(cos(X*pi/Lx)^2 + 4)*cos(Y*pi/(2*(Ly/2)))^2 - (8*(Ly/2)^2*cos(X*pi/Lx)^2)/3 + (128*(Ly/2)^2)/3)*(Ly/2)^4*a^2*sin(Y*pi/(2*(Ly/2)))^4)/4 - 4*Lx^4*pi*(-(3*sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^4*a^4*pi^4)/64 - sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2*a^2*(Ly/2)^2*pi^2 + (Ly/2)^4*(cos(X*pi/Lx)^2 - 7))*cos(X*pi/Lx)^2*(Ly/2)^2*a*sin(Y*pi/(2*(Ly/2)))^3 + (3*Lx^5*(sin(X*pi/Lx)^4*a^4*pi^4*(cos(X*pi/Lx)^2 + 2)*cos(Y*pi/(2*(Ly/2)))^4 - (16*sin(X*pi/Lx)^2*a^2*(Ly/2)^2*pi^2*(cos(X*pi/Lx)^2 - 2)*cos(Y*pi/(2*(Ly/2)))^2)/3 - (64*(Ly/2)^2*(cos(X*pi/Lx)^2*sin(X*pi/Lx)^2*a^2*pi^2 - 2*(Ly/2)^2))/3)*cos(X*pi/Lx)*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)/16 + Lx^6*pi*(sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4/8 + sin(X*pi/Lx)^2*a^2*(Ly/2)^2*pi^2*(cos(X*pi/Lx)^2 + 3)*cos(Y*pi/(2*(Ly/2)))^4 + 2*((sin(X*pi/Lx)^2*a^2*pi^2 - 8*(Ly/2)^2)*cos(X*pi/Lx)^2 + sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2)*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2 - 40*cos(X*pi/Lx)^2*(Ly/2)^4)*sin(X*pi/Lx)^2*a*sin(Y*pi/(2*(Ly/2)))/8 + cos(X*pi/Lx)*sin(X*pi/Lx)^2*Lx^7*(sin(X*pi/Lx)^4*cos(Y*pi/(2*(Ly/2)))^6*a^4*pi^4 - 8*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^4*a^2*(Ly/2)^2*pi^2 + 32*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2*a^2*(Ly/2)^2*pi^2 - 128*(Ly/2)^4)/64)*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/2 + pi*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*(sin(Y*pi/(2*(Ly/2)))*a*pi + cos(X*pi/Lx)*Lx)*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2 + 16*(Ly/2)^2)/4)^2*a^2)*log((a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 8*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 8*Lx^2*(Ly/2)^2 + pi*sqrt((Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2 + 4*cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2)*(a^2*pi^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2*Lx^2 + 4*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 16*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + 16*Lx^2*(Ly/2)^2)*a^2))/(Lx^2*(Ly/2)^2)) - 4*pi*(-pi^5*cos(X*pi/Lx)^8*(sin(X*pi/Lx)^2*Lx^2*ni - 2*(ni + (3*log(2))/2)*(Ly/2)^2)*(Ly/2)^6*a^5*sin(Y*pi/(2*(Ly/2)))^9/2 + Lx*pi^4*cos(X*pi/Lx)^7*((ni + (3*log(2))/2)*(Ly/2)^2*cos(X*pi/Lx)^2 - 4*sin(X*pi/Lx)^2*Lx^2*ni + 8*(ni + (3*log(2))/2)*(Ly/2)^2)*(Ly/2)^6*a^4*sin(Y*pi/(2*(Ly/2)))^8 + Lx^2*(2*pi^2*(sin(X*pi/Lx)^2*ni + 3*log(2) + 1/2)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^2 + (sin(X*pi/Lx)^2*a^2*pi^2 + 24*(Ly/2)^2*(ni + 2*log(2)))*cos(X*pi/Lx)^2 + (-50*Lx^2*ni - 4*(Ly/2)^2)*sin(X*pi/Lx)^2 + 104*(Ly/2)^2*(ni + (18*log(2))/13))*pi^3*cos(X*pi/Lx)^6*(Ly/2)^6*a^3*sin(Y*pi/(2*(Ly/2)))^7/4 + (3*Lx^3*(sin(X*pi/Lx)^2*a^2*pi^2*(cos(X*pi/Lx)^2*log(2) + 2*sin(X*pi/Lx)^2*ni + 6*log(2) + 1)*cos(Y*pi/(2*(Ly/2)))^2 + ((4*sin(X*pi/Lx)^2*a^2*pi^2)/3 + 8*(Ly/2)^2*(ni + 3*log(2)))*cos(X*pi/Lx)^2 + (-(38*Lx^2*ni)/3 - 4*(Ly/2)^2)*sin(X*pi/Lx)^2 + 88*(ni + (12*log(2))/11)*(Ly/2)^2/3)*pi^2*cos(X*pi/Lx)^5*(Ly/2)^6*a^2*sin(Y*pi/(2*(Ly/2)))^6)/2 - (3*Lx^4*(-(3*pi^4*(log(2) + 1/3)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^4)/2 + pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 24*(Ly/2)^2*(log(2) + 1/18))*cos(X*pi/Lx)^2 - 20*(sin(X*pi/Lx)^2*ni + (14*log(2))/5 + 13/30)*(Ly/2)^2)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^2 - 50*((sin(X*pi/Lx)^2*a^2*pi^2 + (32*(Ly/2)^2*(ni + 6*log(2)))/25)*cos(X*pi/Lx)^2 + (-(56*Lx^2*ni)/25 - (48*(Ly/2)^2)/25)*sin(X*pi/Lx)^2 + (32*(Ly/2)^2*(ni + (3*log(2))/5))/5)*(Ly/2)^2/3)*pi*cos(X*pi/Lx)^4*(Ly/2)^4*a*sin(Y*pi/(2*(Ly/2)))^5)/8 + (3*Lx^5*cos(X*pi/Lx)^3*(pi^4*((ni + (3*log(2))/2)*cos(X*pi/Lx)^2 + 6*log(2) + 2)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^4 - 6*pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (28*(Ly/2)^2*(log(2) + 5/42))/3)*cos(X*pi/Lx)^2 - 40*(sin(X*pi/Lx)^2*ni + (12*log(2))/5 + 3/10)*(Ly/2)^2/9)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^2 + (76*(Ly/2)^2*((sin(X*pi/Lx)^2*a^2*pi^2 + (48*(Ly/2)^2*log(2))/19)*cos(X*pi/Lx)^2 + (-(8*Lx^2*ni)/19 - (16*(Ly/2)^2)/19)*sin(X*pi/Lx)^2 + (32*(Ly/2)^2*ni)/19))/3)*(Ly/2)^4*sin(Y*pi/(2*(Ly/2)))^4)/8 - (3*Lx^6*pi*cos(X*pi/Lx)^2*(Ly/2)^2*(pi^4*(ni - 3*log(2) - 3/2)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^6/3 + pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 16*(ni + (3*log(2))/2)*(Ly/2)^2)*cos(X*pi/Lx)^2 + (4*(Ly/2)^2*(ni - 30*log(2) - 21/2))/3)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 + (172*(Ly/2)^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 192*(log(2) + 5/24)*(Ly/2)^2/43)*cos(X*pi/Lx)^2 - 64*(sin(X*pi/Lx)^2*ni + (3*log(2))/2)*(Ly/2)^2/43)*cos(Y*pi/(2*(Ly/2)))^2)/3 - (224*cos(X*pi/Lx)^2*(Ly/2)^4)/3)*sin(X*pi/Lx)^2*a*sin(Y*pi/(2*(Ly/2)))^3)/32 + Lx^7*cos(X*pi/Lx)*(Ly/2)^2*(pi^4*((ni + (3*log(2))/2)*cos(X*pi/Lx)^2 - ni + 3*log(2) + 3/2)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^6 - 6*pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 6*(ni + (5*log(2))/3 + 1/18)*(Ly/2)^2)*cos(X*pi/Lx)^2 + (2*(Ly/2)^2*(ni - 12*log(2) - 9/2))/3)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 - 104*(Ly/2)^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (24*(Ly/2)^2*(log(2) + 1/3))/13)*cos(X*pi/Lx)^2 - (8*(Ly/2)^2*(sin(X*pi/Lx)^2*ni - 1/2))/13)*cos(Y*pi/(2*(Ly/2)))^2 + 32*cos(X*pi/Lx)^2*(Ly/2)^4)*sin(X*pi/Lx)^2*sin(Y*pi/(2*(Ly/2)))^2/16 - Lx^8*pi*(pi^4*(-1 - (3*log(2))/2 + ni)*sin(X*pi/Lx)^4*a^4*cos(Y*pi/(2*(Ly/2)))^6/2 + pi^2*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 24*(ni + log(2) - 1/6)*(Ly/2)^2)*cos(X*pi/Lx)^2 + Lx^2*(ni - 1/2)*sin(X*pi/Lx)^2 + 12*(Ly/2)^2*(ni - 2*log(2) - 7/6))*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 + 76*((a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (48*(Ly/2)^2*(ni + 2*log(2) + 1/6))/19)*cos(X*pi/Lx)^2 + (4*Lx^2*(ni - 1/2)*sin(X*pi/Lx)^2)/19 + (16*(Ly/2)^2*(ni - 3*log(2) - 3/2))/19)*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2 + 512*cos(X*pi/Lx)^2*(ni - 1/2)*(Ly/2)^4)*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^4*a*sin(Y*pi/(2*(Ly/2)))/128 + Lx^9*cos(X*pi/Lx)*(sin(X*pi/Lx)^4*a^4*pi^4*(ni + (3*log(2))/2)*cos(Y*pi/(2*(Ly/2)))^6 - 4*pi^2*(a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - 12*(ni + log(2) - 1/6)*(Ly/2)^2)*sin(X*pi/Lx)^2*a^2*cos(Y*pi/(2*(Ly/2)))^4 - 112*(a^2*pi^2*(ni - 1/2)*sin(X*pi/Lx)^2 - (24*(Ly/2)^2*(log(2) + 1/3))/7)*(Ly/2)^2*cos(Y*pi/(2*(Ly/2)))^2 - 256*(ni - 1/2)*(Ly/2)^4)*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^4/256)*a^2)/(16*rho*((Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2*cos(X*pi/Lx)^2*pi^2*a^2 + 4*Lx*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))*cos(X*pi/Lx)*pi*a + Lx^2*(a^2*pi^2*sin(X*pi/Lx)^2*cos(Y*pi/(2*(Ly/2)))^2 + 16*(Ly/2)^2)/4)^2*(cos(X*pi/Lx)^2*(Ly/2)^2*sin(Y*pi/(2*(Ly/2)))^2 + Lx^2*cos(Y*pi/(2*(Ly/2)))^2*sin(X*pi/Lx)^2/4)^2*(Ly/2)*(a*sin(Y*pi/(2*(Ly/2)))*pi*cos(X*pi/Lx) + Lx)^2*(ni - 1/2)*(1 + ni)*a);