function [A0,A,B]=coeficientes()
%function [A0,A,B]=coeficientes()
%
%Devuelve los coeficientes de la funcion efe
%(definida en el archivo efe.m)
%en terminos de su expansion en terminos de 
%senos, coseno y funciones de Bessel

%Primero introducimos a manos los ceros de las
%funciones de Bessel. Desgraciadamente, matlab
%no tiene una funcion que calcule estos ceros
besselC0n=[2.404825557695773,5.520078110286311,8.653727912911011,11.79153443901428,14.93091770848779,18.07106396791092,21.21163662987926,24.35247153074930,27.49347913204025,30.63460646843198];
besselCmn=[3.831705970207512,7.015586669815613,10.17346813506272,13.32369193631422,16.47063005087763,19.61585851046824,22.76008438059277,25.90367208761838,29.04682853491686,32.18967991097440;
5.135622301840683,8.417244140399855,11.61984117214906,14.79595178235126,17.95981949498783,21.11699705302185,24.27011231357310,27.42057354998456,30.56920449551640,33.71651950922270;
6.380161895923984,9.761023129981667,13.01520072169843,16.22346616031877,19.40941522643501,22.58272959310444,25.74816669929498,28.90835078092176,32.06485240709771,35.21867073861011;
7.588342434503804,11.06470948850119,14.37253667161759,17.61596604980483,20.82693295696239,24.01901952477111,27.19908776598125,30.37100766711725,33.53713771181922,36.69900112874465;
8.771483815959954,12.33860419746694,15.70017407971167,18.98013387517992,22.21779989656127,25.43034115422270,28.62661830729114,31.81171672404776,34.98878129455930,38.15986856196713;
9.936109524217684,13.58929017054122,17.00381966781601,20.32078921356651,23.58608443558139,26.82015198341140,30.03372238657047,33.23304176284712,36.42201966825846,39.60323941607540; 
11.08637001924509,14.82126872701317,18.28758283248173,21.64154101984840,24.93492788767302,28.19118845948320,31.42279419226558,34.63708935206932,37.83871738285361,41.03077369158554;
12.22509226400465,16.03777419088771,19.55453643099706,22.94517313187462,26.26681464117664,29.54565967099855,32.79580003734146,36.02561506386957,39.24044799517814,42.44388774327356;
13.35430047743524,17.24122038248913,20.80704778926411,24.23388525775055,27.58374896357301,30.88537896769667,34.15437792385510,37.40009997715659,40.62855371896453,43.84380142033735;
14.47550068655467,18.43346366696658,22.04698536469780,25.50945055418283,28.88737506353046,32.21185619971273,35.49990920537385,38.76180701788165,42.00419023667181,45.23157410353504];

%Numero de intervalos en que descomponemos
%los intervalos [0,1] y [0,2*pi]
M=50;

pasor=1/M;
r=0:pasor:1;
pasotheta=2*pi/M;
theta=0:pasotheta:2*pi;
[R,THETA]=meshgrid(r,theta);

%Primero calculamos los coeficientes A0
A0=zeros(10,1);
for n=1:10
    Cosenos=cos(0*THETA);
    Bessel=besselj(0,besselC0n(n)*R);
    integrando=efe(r,theta).*Cosenos.*Bessel.*R;
    %Integral por la regla del trapecio
    integral=0;
    for k=1:M
        for l=1:M
            integral = integral + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                integrando(k,l+1)+integrando(k+1,l+1))/2;
        end
    end
    %Ahora por la integral en el dividendo
    %Valores de la funcion (cos(m*theta)*besselj(m,r))^2
    integrando=((Cosenos.*Bessel).^2).*R;
    %Integral por la regla del trapecio
    integral2=0;
    for k=1:M
        for l=1:M
            integral2 = integral2 + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                integrando(k,l+1)+integrando(k+1,l+1))/2;
        end
    end
    A0(n)=integral/integral2;
end
%Ahora calculamos los coeficientes A
A=zeros(10,10);
for m=1:10
    for n=1:10
        Cosenos=cos(m*THETA);
        Bessel=besselj(m,besselCmn(m,n)*R);
        integrando=efe(r,theta).*Cosenos.*Bessel.*R;
        %Integral por la regla del trapecio
        integral=0;
        for k=1:M
            for l=1:M
                integral = integral + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                    integrando(k,l+1)+integrando(k+1,l+1))/2;
            end
        end
        %Ahora por la integral en el dividendo
        %Valores de la funcion (cos(m*theta)*besselj(m,r))^2
        integrando=((Cosenos.*Bessel).^2).*R;
        %Integral por la regla del trapecio
        integral2=0;
        for k=1:M
            for l=1:M
                integral2 = integral2 + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                    integrando(k,l+1)+integrando(k+1,l+1))/2;
            end
        end
        A(m,n)=integral/integral2;
    end
end

%Finalmente, los coeficientes B
B=zeros(10,10);
for m=1:10
    for n=1:10
        Senos=sin(m*THETA);
        Bessel=besselj(m,besselCmn(m,n)*R);
        integrando=efe(r,theta).*Senos.*Bessel.*R;
        %Integral por la regla del trapecio
        integral=0;
        for k=1:M
            for l=1:M
                integral = integral + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                    integrando(k,l+1)+integrando(k+1,l+1))/2;
            end
        end
        %Ahora por la integral en el dividendo
        %Valores de la funcion (cos(m*theta)*besselj(m,r))^2
        integrando=((Senos.*Bessel).^2).*R;
        %Integral por la regla del trapecio
        integral2=0;
        for k=1:M
            for l=1:M
                integral2 = integral2 + pasor*pasotheta*(integrando(k,l)+integrando(k+1,l)+...
                    integrando(k,l+1)+integrando(k+1,l+1))/2;
            end
        end
        B(m,n)=integral/integral2;
    end
end

