function FF=fun(s,t)
%--------------------------------------------------
%Compute the exact value of the integral for the function #s 
% on the 
%  t=1 cone    
%  t=2 cylinder
%  t=3 torus
%  t=4 sphere
%--------------------------------------------------
fprintf('Do you prefer to compute the integral in this session? 1=Compute/0=NotCompute \n');
K=input('Y/N?');

if(K==0)
    switch t
        case 1
            switch s
                case 1
                    FF=0.637848242027904;
                case 2
                    FF=1.297528862046102;
                case 3
                    FF=1.159852041305101;
                case 4
                    FF=0.145771246135679;
            end
        case 2
            switch s
                case 1
                    FF=0.712504034681439;
                case 2
                    FF=0.578397783759171;
                case 3
                    FF=1.318195439838512;
                case 4
                    FF=-0.016702772737230;
            end
        case 3
            switch s
                case 1
                    FF=0.343528272832302;
                case 2
                    FF=4.700517385919818;
                case 3
                    FF=1.306729852997971;                        
                case 4
                    FF=-0.012692747954474;
            end
        case 4
            switch s
                case 1
                    FF=0.729465700516434;
                case 2
                    FF=0.280864805817563;
                case 3
                    FF=1.340544583085376;
                case 4
                    FF=0.089499074722825;
            end
    end
else
    switch t
        case 1
            switch s
                case 1
                    %Function sqrt((z+1)(1-z))cos(x/2+y/3+z/3) 
                    fun1=@(u,v) sqrt(-(-1+u.^2).*(abs(u).^2+abs(u.*cos(v)).^2+ ...
                        abs(u.*sin(v)).^2)).*cos(1/30.*u.*(6+15.*cos(v)+10.*sin(v)));
                    Ff=dblquad(fun1,0,1,0,2*pi,1e-11);
                case 2
                    %Function cos(30*x*y*z) for z<1/2 and sqrt(x^2+y^2+z^2)^3 for z=>1/2
                    fun2=@(u,v) sqrt(abs(u).^2+abs(u.*cos(v)).^2+abs(u.*sin(v)).^2).*(cos(30.*u.^3.*...
                        cos(v).*sin(v)).*(u<1/2)+sqrt(2.*u.^2).^3.*(u>=1/2));   
                    Ff=dblquad(fun2,0,1,0,2*pi,1e-11);
                case 3
                    %Function exp(-sin(2x^2+3y^2+5z^2))
                    fun3=@(u,v) exp(-sin(2.*(u.*cos(v)).^2+3.*(u.*sin(v)).^2+5.*u.^2)).*sqrt(abs(u).^2+...
                        abs(u.*cos(v)).^2+abs(u.*sin(v)).^2);     
                    Ff=dblquad(fun3,0,1,0,2*pi,1e-11);
                case 4
                    %Function exp(-sqrt(x^2+y^2+z^2))*1/(1+x.^2)*cos(1+x^2)*sin(1-y^2)*exp(abs(z))
                    fun4=@(u,v) exp(-sqrt(2.*u.^2))./(1+(u.*cos(v)).^2).*cos(1+(u.*cos(v)).^2).*sin(1-(u.*sin(v)).^2).*...
                        exp(abs(u)).*sqrt(abs(u).^2+abs(u.*cos(v)).^2+abs(u.*sin(v)).^2);             
                    Ff=dblquad(fun4,0,1,0,2*pi,1e-11);
            end
                Area=@(u,v) sqrt(abs(u).^2+abs(u.*cos(v)).^2+abs(u.*sin(v)).^2);
                Ar=dblquad(Area,0,1,0,2*pi,1e-11);
            case 2
                switch s
                    case 1
                        %Function sqrt((z+1)(1-z))cos(x/2+y/3+z/3) 
                        fun1=@(u,v) sqrt((1-v.^2)).*cos(cos(u)./2+sin(u)./3+v./5).*sqrt(abs(cos(u)).^2+abs(sin(u)).^2);
                        Ff=dblquad(fun1,0,2*pi,-1,1,1e-11);
                    case 2
                        %Function cos(30*x*y*z) for z<1/2 and sqrt(x^2+y^2+z^2)^3 for z=>1/2
                        fun2=@(u,v) sqrt(abs(cos(u)).^2+abs(sin(u)).^2).*(cos(30.*cos(u).*...
                            sin(u).*v).*(v<1/2)+sqrt(1+v.^2).^3.*(v>=1/2));
                        Ff=dblquad(fun2,0,2*pi,-1,1,1e-11);
                    case 3
                        %Function exp(-sin(2x^2+3y^2+5z^2))
                        fun3=@(u,v) exp(-sin(2.*cos(u).^2+3.*sin(u).^2+5.*v.^2)).*sqrt(abs(cos(u)).^2+abs(sin(u)).^2);
                        Ff=dblquad(fun3,0,2*pi,-1,1,1e-11);
                    case 4
                        %Function exp(-sqrt(x^2+y^2+z^2))*1/(1+x^2)*cos(1+x)*sin(1-y)*exp(abs(z))
                        fun4=@(u,v) exp(-sqrt(1+v.^2))./(1+cos(u).^2).*cos(1+cos(u).^2).*sin(1-sin(u).^2).*...
                            exp(abs(v)).*sqrt(abs(cos(u)).^2+abs(sin(u)).^2);
                        Ff=dblquad(fun4,0,2*pi,-1,1,1e-11);
                end
                Area=@(u,v) sqrt(abs(cos(u)).^2+abs(sin(u)).^2);
                Ar=dblquad(Area,0,2*pi,-1,1,1e-11);
        case 3
            switch s
                case 1
                    %Function sqrt((z+1)(1-z))cos(x/2+y/3+z/5) 
                    fun1=@(u,v) sqrt(1-sin(u).^2).*sqrt(abs(cos(u).*(2+cos(u)).*cos(v)).^2+... 
                        abs((2+cos(u)).*sin(u)).^2+abs(cos(u).*(2+cos(u)).*sin(v)).^2).*cos((2+cos(u)).*cos(v)./...
                        2+(2+cos(u)).*sin(v)./3+sin(u)./5);
                    Ff=dblquad(fun1,0,2*pi,0,2*pi,1e-11);
                case 2
                    %Function cos(30*x*y*z) for z<1/2 and sqrt(x^2+y^2+z^2)^3 for z=>1/2
                    fun2=@(u,v) (cos(30.*(2+cos(u)).^2.*cos(v).*sin(v).*sin(u)).*...
                        (sin(u)<1/2)+sqrt((2+cos(u)).^2+sin(u).^2).^3.*(sin(u)>=1/2)).*...
                        sqrt(abs(cos(u).*(2+cos(u)).*cos(v)).^2+abs((2+cos(u)).*sin(u)).^2+...
                        abs(cos(u).*(2+cos(u)).*sin(v)).^2);
                    Ff=dblquad(fun2,0,2*pi,0,2*pi,1e-11);
                case 3
                    %Function exp(-sin(2x^2+3y^2+5z^2))
                    fun3=@(u,v) exp(-sin(2.*(2+cos(u)).^2.*cos(v).^2+3.*(2+cos(u)).^2.*...
                        sin(v).^2+5.*sin(u).^2)).*sqrt(abs(cos(u).*(2+cos(u)).*...
                        cos(v)).^2+abs((2+cos(u)).*sin(u)).^2+abs(cos(u).*(2+cos(u)).*sin(v)).^2);
                    Ff=dblquad(fun3,0,2*pi,0,2*pi,1e-11);
                case 4
                    %Function exp(-sqrt(x^2+y^2+z^2))*1/(1+x.^2)*cos(1+x^2)*sin(1-y^2)*exp(abs(z))
                    fun4=@(u,v) exp(-sqrt((2+cos(u)).^2+sin(u).^2))./(1+(2+cos(u)).^2.*...
                        cos(v).^2).*cos(1+(2+cos(u)).^2.*cos(v).^2).*sin(1-(2+cos(u)).^2.*...
                        sin(v).^2).*exp(abs(sin(u))).*sqrt(abs(cos(u).*(2+cos(u)).*cos(v)).^2+...
                        abs((2+cos(u)).*sin(u)).^2+abs(cos(u).*(2+cos(u)).*sin(v)).^2);
                    Ff=dblquad(fun4,0,2*pi,0,2*pi,1e-11);
            end
            Area=@(u,v) sqrt(abs(cos(u).*(2+cos(u)).*cos(v)).^2+...
                abs((2+cos(u)).*sin(u)).^2+abs(cos(u).*(2+cos(u)).*sin(v)).^2);
            Ar=dblquad(Area,0,2*pi,0,2*pi,1e-11);    
    case 4
            switch s
                case 1
                    %Function sqrt((z+1)(1-z))cos(x/2+y/3+z/3) 
                    fun1=@(u,v) sqrt((1-cos(u)).*(1+cos(u))).*cos(cos(u)./5+...
                        1/2.*cos(v).*sin(u)+1/3.*sin(u).*sin(v)).*sqrt(abs(cos(u).*sin(u)).^2+...
                        abs(cos(v).*sin(u).^2).^2+abs(sin(u).^2.*sin(v)).^2);
                    Ff=dblquad(fun1,0,pi,0,2*pi,1e-11);
                case 2
                    %Function cos(30*x*y*z) for z<1/2 and sqrt(x^2+y^2+z^2)^3 for z=>1/2
                    fun2=@(u,v) sqrt(abs(cos(u).*sin(u)).^2+abs(cos(v).*sin(u).^2).^2+...
                    abs(sin(u).^2.*sin(v)).^2).*((cos(30.*sin(u).^2.*cos(u).*cos(v).*sin(v))).*...
                    (cos(u)<1/2)+1.*(cos(u)>=1/2));   
                    Ff=dblquad(fun2,0,pi,0,2*pi,1e-11);
                case 3
                    %Function exp(-sin(2x^2+3y^2+5z^2))
                    fun3=@(u,v) sqrt(abs(cos(u).*sin(u)).^2+abs(cos(v).*sin(u).^2).^2+...
                    abs(sin(u).^2.*sin(v)).^2).*(exp(-sin(2.*(sin(u).*cos(v)).^2+...
                    3.*(sin(u).*sin(v)).^2+5.*(cos(u)).^2)));     
                    Ff=dblquad(fun3,0,pi,0,2*pi,1e-11);
                case 4
                    %Function exp(-sqrt(x^2+y^2+z^2))*1/(1+x.^2)*cos(1+x^2)*sin(1-y^2)*exp(abs(z))
                    fun4=@(u,v) sqrt(abs(cos(u).*sin(u)).^2+abs(cos(v).*sin(u).^2).^2+...
                    abs(sin(u).^2.*sin(v)).^2).*(exp(-1)./(1+(sin(u).*cos(v)).^2).*cos(1+...
                    (sin(u).*cos(v)).^2).*sin(1-(sin(u).*sin(v)).^2).*exp(abs(cos(u))));             
                    Ff=dblquad(fun4,0,pi,0,2*pi,1e-11);
            end
                Area=@(u,v) sqrt(abs(cos(u).*sin(u)).^2+abs(cos(v).*sin(u).^2).^2+...
                    abs(sin(u).^2.*sin(v)).^2);
                Ar=dblquad(Area,0,pi,0,2*pi,1e-11);
    end
FF=Ff./Ar;
end

end