function PlotFunctionFGen(w)

%%%% Plot FUNCTION F (Formula 0.1 and 0.3 of the corrigendum: corrected version of Formula B.1 and B.3 of the JFM
%%%% (2013) paper) versus r for a fixed value of the (positive) frequency w 



if w<2
r1=(1-w^2/4)^(1/2)
r1m=[0.001:0.001:(r1-0.001)];
r1p=[r1+0.001:0.001:1];
Fp1=FunctionF(r1p,w);
Fm1=FunctionFwsmall(r1m,w); 
figure(1)
hold off
plot(r1m,Fm1,'k-','linewidth',2)
hold on
plot(r1p,Fp1,'k-','linewidth',2)
axis([0 1 -1 1])

else 
r=[0:0.001:0.999];
F=FunctionF(r,w);
figure(1)
hold off
plot(r,F,'k-','linewidth',2)
axis([0 1 -1 1])
end
xlabel('r')
ylabel('F')

end


function F=FunctionF(r,w)


%%%% FUNCTION F FOR OMEGA>2

theta=asin(r);

sp=sqrt(w./2+cos(theta));
sm=sqrt(w./2-cos(theta));

%F1a=-sin(theta).^2.*(sp-sqrt(cos(theta))).*((sp.^3+sm.^3).*(1-cos(theta).^2)+4.*sm.^2.*sp.^2.*(sm-sp).*cos(theta))./(32.*sm.^3.*sp.^3.*(sp.^2-cos(theta)));
%F1b=sin(theta).^2.*(sm-sqrt(cos(theta))).*((sp.^3+sm.^3).*(1-cos(theta).^2)+4.*sm.^2.*sp.^2.*(sm-sp).*cos(theta))./(32.*(sm.^2+cos(theta)).*sp.^3.*sm.^3);
FCor=-((sp.^3+sm.^3).*sin(theta).^4./(32.*sp.^3.*sm.^3)+(sm-sp).*cos(theta).*sin(theta).^2./(8.*sp.*sm)).*((sp-sqrt(cos(theta)))./(sp.^2-cos(theta))-(sm-sqrt(cos(theta)))./(sm.^2+cos(theta)));
F1c=sin(theta).^2.*(sp+sqrt(cos(theta))).*(2.*sp.^2-2.*sp.*sqrt(cos(theta))+cos(theta)).*(1+4.*sp.^2.*cos(theta)-cos(theta).^2)./(32.*sp.^3.*(cos(theta).^2+4.*sp.^4));
F1d=sin(theta).^2.*(sm+sqrt(cos(theta))).*(2.*sm.^2-2.*sm.*sqrt(cos(theta))+cos(theta)).*(1-4.*sm.^2.*cos(theta)-cos(theta).^2)./(32.*sm.^3.*(cos(theta).^2+4.*sm.^4));
F1e=-((sm.^2+sp.^2).*(sm.^6+7.*sp.^2.*sm.^4+8.*sm.^7.*sp.^3-7.*sm.^2.*sp.^4+16.*sm.^5.*sp.^5-sp.^6+4.*sp.^3.*sm.^3.*(-3+2.*sp.^4))+(4.*sm.^6.*sp-sm.^4.*sp.^3-16.*sm.^2.*sp.^5-3.*sp.^7+4.*sm.^5.*sp.^2.*(-3+16.*sp.^4)+sp.^4.*sm.^3.*(-19+32.*sp.^4)+sm.^7.*(-1+32.*sp.^4)).*sqrt(cos(theta))+(4.*sm.^8.*sp.^2-3.*sp.^6+sm.^6.*(1-20.*sp.^4)+sm.^2.*sp.^4.*(-11+4.*sp.^4)+4.*sp.^3.*sm.^3.*(-3+16.*sp.^4)+sm.^5.*(-4.*sp+64.*sp.^5)+sm.^4.*(sp.^2-20.*sp.^6)).*cos(theta)-(4.*sm.^7.*sp.^2+4.*sm.^6.*sp.^3+sp.^5+56.*sm.^4.*sp.^5+sm.^3.*sp.^2.*(3-76.*sp.^4)+sm.^5.*(1-8.*sp.^4)+3.*sp.^3.*sm.^2.*(1-4.*sp.^4)).*cos(theta).^(3./2)+(-sm.^8-4.*sm.^6.*sp.^2-4.*sp.^3.*sm.^5-48.*sm.^4.*sp.^4+60.*sp.^5.*sm.^3+20.*sm.^2.*sp.^6+sp.^8).*cos(theta).^2+(sm.^7-4.*sm.^6.*sp+8.*sm.^5.*sp.^2-11.*sm.^4.*sp.^3+31.*sp.^4.*sm.^3+20.*sm.^2.*sp.^5+3.*sp.^7).*cos(theta).^(5./2)+(-sm.^6+4.*sm.^5.*sp-sp.^2.*sm.^4+12.*sp.^3.*sm.^3+11.*sm.^2.*sp.^4+3.*sp.^6).*cos(theta).^3+(sm.^5+3.*sm.^3.*sp.^2+3.*sp.^3.*sm.^2+sp.^5).*cos(theta).^(7./2)).*sin(theta).^2./(32.*sm.^3.*sp.^3.*(sm.^2+sp.^2).*(sm.^2+sp.^2+2.*sp.*sqrt(cos(theta))+cos(theta)).^2);
F1f= (1/4)*sin(theta).^2;
%%FCor=F1a+F1b

F=FCor+F1c+F1d+F1e+F1f;

end 

function F=FunctionFwsmall(r,w)




theta=asin(r);
sp=sqrt(w./2+cos(theta));
sm=sqrt(-w./2+cos(theta));


%%%% Relation (0.3) 

F1=sin(theta).^2./32.*8;
F2= sin(theta).^2./32.*((sp+sqrt(cos(theta))).*(4.*sp.^2.*cos(theta)+sin(theta).^2)./(sp.^3.*(2.*sp.^2+2.*sp.*sqrt(cos(theta))+cos(theta))));
F3=sin(theta).^2./32.*((sm+sqrt(cos(theta))).*(4.*sm.^2.*cos(theta)+sin(theta).^2)./(sm.^3.*(2.*sm.^2+2.*sm.*sqrt(cos(theta))+cos(theta))));
F4=sin(theta).^2./32.*(-4.*(sm.^3+sm.^2.*sqrt(cos(theta))+(sm+sp).*cos(theta))./(sp.*sm.*(sm+sqrt(cos(theta)))));
F5=sin(theta).^2./32.*(-sm.^4.*sin(theta).^2.*(sm+sqrt(cos(theta))+cos(theta).*(sp.^3+sm.^3)./sm.^4)./(cos(theta).*sm.^3.*sp.^3.*(sm+sqrt(cos(theta)))));
F6=sin(theta).^2./32.*(-4.*(sp.^3+sp.^2.*sqrt(cos(theta))+(sm+sp).*cos(theta))./(sp.*sm.*(sp+sqrt(cos(theta)))));
F7=sin(theta).^2./32.*(-sp.^4.*sin(theta).^2.*(sp+sqrt(cos(theta))+cos(theta).*(sp.^3+sm.^3)./sp.^4)./(cos(theta).*sm.^3.*sp.^3.*(sp+sqrt(cos(theta)))));
F8=sin(theta).^2./32.*(1./(cos(theta).*sm.^3.*sp.^3.*(sm+sp).*(sm+sp+sqrt(cos(theta))).^2).*(-8.*sm.^2.*sp.^2.*(sm+sp).*cos(theta).^2.*(sp-sm).^2.*(cos(theta)-1)+4.*sm.^2.*sp.^2.*(sm-sp).^2.*cos(theta).^(5./2).*(1+2.*(sm+sp).*sqrt(cos(theta)))+cos(theta).*(sm.*sp.*(sm+sp).*(2.*(2.*sm.^2-3.*sm.*sp+2.*sp.^2).*sin(theta).^2+4.*sm.^5.*sp-8.*sm.^3.*sp.^3+4.*sm.*sp.^5)+2.*(sm.^5+sp.^5).*sin(theta).^2)+(sm+sp).^2.*(sm.^4+sp.^4).*sin(theta).^2.*((sm+sp)+2.*sqrt(cos(theta)))+cos(theta).^(3./2).*(8.*sp.^2.*sm.^2.*(sm.^2-sp.^2).^2+(sm.^4+sp.^4+sm.*sp.*(sm-sp).^2).*sin(theta).^2)));

F=F1+F2+F3+F4+F5+F6+F7+F8;

end
