matlab椭圆积分程序

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