To perform a quad-summation

Does anyone have idea how to perform the following summation?

where,

I used the following code-

Phi01=0.4*pi;
Phi02=0.8*pi;
N1=10;
s=0.4;
P1=linspace(-4,4,N1);
i1=0;
for x=P1
    i1=i1+1;
    tic
    syms  v1 v m1 m
    Bmv=(2*m)-v;
    Bm1v1=(2*m1)-v1;
    psimvm1v1=((pi*0.5)*(m-m1))-((2*(Bmv-Bm1v1)*(Bmv+Bm1v1+1)*x*(x^2+1)*log(1-s))...
        *(((x^2)+(2*Bmv+1)^2)*((x^2)+(2*Bm1v1+1)^2))^(-1));
    lambdamvm1v1=((Bmv+Bm1v1+1)*(x^2+1)*(x^2+(2*Bmv+1)*(2*Bm1v1+1)))*((x^2+(2*Bmv+1)^2)*(x^2+(2*Bm1v1+1)^2))^(-1);
    
    
    Fmvm1v1=((Phi01)^(v+v1)*(Phi02)^(m+m1-v-v1))...
        *(factorial(v)*factorial(v1)*factorial(m-v)*factorial(m1-v1)*(Bmv+Bm1v1+1)*(x^2+1)^(Bmv+Bm1v1))^(-1);
    S1 = symsum((Fmvm1v1*(1-s)^(lambdamvm1v1)*cos(psimvm1v1)),v1,0,m1);
    S2 = symsum(S1,v,0,m);
    S3 = symsum(S2,m1,0,17);
    S4 = symsum(S3,m,0,17);
    T=(s)^(-1)*(1-S4);
    T1(i1)=T;
    toc
end

plot(P1(1:i1),T1(1:i1),'b','LineWidth',3);

However, this does not work!

[Walter Roberson]

No definition for s ?

You can compute with a symbolic s. The steps are pretty slow, but it can be done. But you cannot plot() the result. fplot() over a range of s perhaps.

[Walter Roberson]

Make that "the steps are quite slow", at least for the last one. Hours. And that is just for one particular x value.

[Subhajit Bej]

Hi, s is a constant and can take any value between 0 and 1. I had assumed s=0.4

[Asvin Kumar]

My observations are same as Walter's. Nothing seems to be wrong with the code. It's just very slow and takes more than 12 hours for one particular x value or at least that's when I stopped the function.

[Walter Roberson]

The graphs look pretty different for fairly small values of s by the way.

[Subhajit Bej]

Yes, they should look diffrent as s is related to the size of a circular aperture which is responsible for diffraction of light.

Answer 1

[Subhajit Bej]

Hello all,

Thank you for your responses. BTW, I modified the code to make it faster. Here it is-

Phi01=-0.4*pi;
Phi02=0.8*pi;
s=0.8;
N1=300;
P1=linspace(-4,4,N1);
SF5=zeros(size(P1));
SF6=zeros(size(P1));
tic
i1=0;
for x=P1
    i1=i1+1;
    SF4=0;
    for m=0:15
        SF3=0;
        for m1=0:15
            SF2=0;
            for v=0:m
                SF1=0;
                for v1=0:m1
                    Bmv=(2*m)-v;
                    Bm1v1=(2*m1)-v1;
                    ps1mvm1v1=((pi*0.5)*(m-m1))-((2*(Bmv-Bm1v1)*(Bmv+Bm1v1+1)*x*(x^2+1)*log(1-s))...
                        *(((x^2)+(2*Bmv+1)^2)*((x^2)+(2*Bm1v1+1)^2))^(-1));
                    lambdamvm1v1=((Bmv+Bm1v1+1)*(x^2+1)*(x^2+(2*Bmv+1)*(2*Bm1v1+1)))*((x^2+(2*Bmv+1)^2)*(x^2+(2*Bm1v1+1)^2))^(-1);
                    
                    
                    Fmvm1v1=((Phi01)^(v+v1)*(Phi02)^(m+m1-v-v1))...
                        *(factorial(v)*factorial(v1)*factorial(m-v)*factorial(m1-v1)*(Bmv+Bm1v1+1)*(x^2+1)^(Bmv+Bm1v1))^(-1);
                    SF1=SF1+(Fmvm1v1*((1-s)^(lambdamvm1v1))*cos(ps1mvm1v1));
                end
                SF2=SF2+SF1;
            end
            SF3=SF3+SF2+SF1;
        end
        SF4=SF4+SF3+SF2+SF1;
    end
    SF5(i1)=SF4;
    SF6(i1)=(s^(-1))*(1-(SF5(i1)))
    toc
end
divd1=(SF5(1)+SF5(end))*0.5;
divd2=(SF6(1)+SF6(end))*0.5;
plot(P1(1:i1),(s^(-1))*(1-(SF5(1:i1))),'g--o','LineWidth',2)
xlabel('x=z/z_r')
ylabel('\Delta T/T_0')

Though the code works and is much faster with m=m1=15 which is enough to have acurate results, I do not repeat the results from this article- Fig. 6. Don't know what is wrong in the implementation!