function?Fun=EllipticInt（stringphikerrtol）
phi=double（phi）;k=double（k）;
if?strcmp（string‘E‘）
????tmp1=RF（cos（phi）^2/sin（phi）^21/sin（phi）^2-k*k1/sin（phi）^2）;
????tmp2=RD（cos（phi）^2/sin（phi）^21/sin（phi）^2-k*k1/sin（phi）^2）;
????Fun=tmp1-k*k*tmp2/3;
????if?errtol<=1e-15
????????tmp=mfun（‘EllipticE‘sin（phi）k）;
????????if?max（abs（tmp-Fun））<1e-6
????????????Fun=tmp;
????????end
????end
elseif?strcmp（string‘F‘）
????Fun=RF（cos（phi）^2/sin（phi）^21/sin（phi）^2-k*k1/sin（phi）^2）;
????if?errtol<=1e-15
????????tmp=mfun（‘EllipticF‘sin（phi）k）;
????????if?max（abs（tmp-Fun））<1e-6
????????????Fun=tmp;
????????end
????end
elseif?strcmp（string‘EF‘）
????tmp1=RF（cos（phi）^2/sin（phi）^21/sin（phi）^2-k*k1/sin（phi）^2）;
????tmp2=RD（cos（phi）^2/sin（phi）^21/sin（phi）^2-k*k1/sin（phi）^2）;
????Fun=[tmp1-k*k*tmp2/3?tmp1];
????if?errtol<=1e-15
????????tmp=[mfun（‘EllipticE‘sin（phi）k）?mfun（‘EllipticF‘sin（phi）k）];
????????if?max（abs（tmp-Fun））<1e-6
????????????Fun=tmp;
????????end
????end
end

????function?Vrd=RD（xyz）
????????x=double（x）;y=double（y）;z=double（z）;
????????x_old=x;
????????y_old=y;
????????z_old=z;
????????Vrd=0;
????????for?count=1:1:1000
????????????lambda=sqrt（x_old*y_old）+sqrt（x_old*z_old）+sqrt（y_old*z_old）;
????????????u=（x_old+y_old+3*z_old）/5;
????????????x_new=（x_old+lambda）/4;
????????????y_new=（y_old+lambda）/4;
????????????z_new=（z_old+lambda）/4;
????????????X=1-x_o