gravity

代码说明：

  function dg=gplus(x,y,z,x1,x2,y1,y2,z1,z2,p) G=6.67e-11 t222=myfun(x,y,z,x2,y2,z2) t122=myfun(x,y,z,x1,y2,z2) t212=myfun(x,y,z,x2,y1,z2) t221=myfun(x,y,z,x2,y2,z1) t211=myfun(x,y,z,x2,y1,z1) t121=myfun(x,y,z,x1,y2,z1) t112=myfun(x,y,z,x1,y1,z2) t111=myfun(x,y,z,x1,y1,z1) dg=-G*p*1e6*(t222-t122-t212-t221+t211+t121+t112-t111) 转化为g.u. function y=myfun(x,y,z,xi,eta,zeta) R=sqrt((xi-x).^2+(eta-y).^2+(zeta-z).^2) temp1=(xi-x).*log((eta-y)+R) temp2=(eta-y).*log((xi-x)+R) temp3=(zeta-z).*atan((zeta-z).*R/((xi-x).*(eta-y))) 书中所给公式 temp3=-(zeta-z).*atan((xi-x).*(eta-y)./(zeta-z)./R) 上式变换 y=temp1+temp2+temp3 (function dg = gplus (x, y, z, x1, x2, y1, y2, z1, z2, p) G = 6.67e-11 t222 = myfun (x, y, z, x2, y2, z2) t122 = myfun (x, y, z, x1, y2, z2) t212 = myfun (x, y, z, x2, y1, z2) t221 = myfun (x, y, z, x2, y2, z1) t211 = myfun (x, y, z, x2, y1, z1) t121 = myfun (x, y, z, x1, y2, z1) t112 = myfun (x, y, z, x1, y1, z2) t111 = myfun (x, y, z, x1, y1, z1) dg =-G* p* 1e6* (t222-t122-t212-t221+t211+t121+t112-t111) converted the Gu function y = myfun, (x y, z, xi, eta, zeta) R = sqrt ((xi-x). ^ 2+ (eta-y). ^ 2+ (zeta-z). ^ 2) temp1 = (xi-x).* log ((eta-y)+R) temp2 = (eta-y).* log ((xi-x)+R) temp3 = (zeta-z).* atan ((zeta-z).* R/( (xi-x).* (eta-y))) book given formula Temp3 =- (Zeta-z)* ATAN ((Xi-X)* (ETA-y)/(Zeta- Z)/R) the above formula transform Y = temp1+temp2+temp3 )