function [b,alpha,s,e,map,condM] = aproximacao(parametros,NP,...
    p,np_comp,CD,LD,base,A,f,M,pasta,plota,estagio,grava_tex)

fprintf(2,'--- inicio de aproximacao -----------------------\n');
nb = size(base.elementos,2);
np = size(p,1);
nf = size(f,2);
fprintf(2,'%d elementos na base\n',nb);
fprintf(2,'%d pontos de amostragem\n',np);
fprintf(2,'%d funcoes a aproximar\n',nf);
assert(size(A,1) == np)
assert(size(A,2) == nb)
assert(size(M,1) == nb)
assert(size(M,2) == nb)
assert(size(f,1) == np)
assert(size(p,2) == 2)
fprintf(2,'analizando as matrizes A e M\n');
nao_nulos_A = nnz(A)/size(A,2); %nao nulos por coluna
nao_nulos_M = nnz(M)/size(M,1); %nao nulos por coluna
fprintf(2,'nao nulos por coluna: A = %.1f M = %.1f\n',...
    nao_nulos_A,nao_nulos_M);

% figure(8); clf;
% analisa_esparsidade(M);
% if plota.salva_eps
%     saveas(gcf,sprintf('%s/histM.eps',pasta));
% end
if nb <= 1600
    fprintf(2,'calculando condicao de M\n');
    rng(0,'twister')
    condM = condest(M);
    fprintf(2,'condM = %8.6f\n',condM);
    %    auto_valor = eigs(M);  %eigs autovalor da matriz esparca
    %    maior_autovalor = max(auto_valor);
    %    menor_autovalor = min(auto_valor);
    %    fprintf(2,'condM = %.1e min_av = %.3e max_av = %.3e \n',condM,...
    %    menor_autovalor,maior_autovalor);
    
else
    fprintf(2,'condicao de M nao calculada, muito grande\n');
    condM = 0;
    maior_autovalor = 0;
    menor_autovalor = 0;
end

if grava_tex
    nomearq = sprintf('%s/parametros.tex',pasta);
    fprintf(2,'gravando parametros em %s\n',nomearq);
    arq = fopen(nomearq,'wt');
    %fprintf(arq,'\\newcommand{\\PARMmaiorautovalor}{%s}\n',...
    %   converte_erro_tex(maior_autovalor));
    %fprintf(arq,'\\newcommand{\\PARMmenorautovalor}{%s}\n',...
    %   converte_erro_tex(menor_autovalor));
    fprintf(arq,'\\newcommand{\\PARMcondicao}{%s}\n',...
        converte_erro_tex(condM));
    fprintf(arq,'\\newcommand{\\PARMnnA}{%d}\n',nnz(A));
    fprintf(arq,'\\newcommand{\\PARMnnM}{%d}\n',nnz(M));
    fprintf(arq,'\\newcommand{\\PARMnnpcA}{%s}\n',...
        converte_nnz_tex(nao_nulos_A));
    fprintf(arq,'\\newcommand{\\PARMnnpcM}{%s}\n',...
        converte_nnz_tex(nao_nulos_M));
end

%determina a aproximacao com esta base
Minv = inverte_matriz_M(M);
s = zeros(np,nf);
if (parametros.metodo == 1) || (parametros.metodo == 2)
    niter = 1;
else
    niter = 6;
end
iter = 0;
while true
    iter = iter + 1;
    if iter == 1
        if parametros.metodo == 2
            fprintf(2,'calculando o lado direito do sistema (IC)\n');
            b = f;
        else
            fprintf(2,'calculando o lado direito do sistema (MQ)\n');
            b = calcula_vetor_de_aprox_MQ(A,f,np_comp);
        end
    else
        sigma = ((niter + 1 - iter)^2)*parametros.tol;
        if parametros.metodo == 4
            alpha = 0.25;
            fprintf(2,'calculando os pesos dos pontos\n');
            w = recalcula_peso_WQ(f,s,alpha,sigma);
            fprintf(2,'recalculando a matriz de momentos\n');
            M = (A'*diag(w)*A)/np_comp;
            Minv = inverte_matriz_M(M);
            fprintf(2,'calculando o lado direito do sistema (WQ)\n');
            b = calcula_vetor_de_aprox_WQ(w,A,f,np_comp);
        else
            assert(parametros.metodo == 3)
            fprintf(2,'calculando o lado direito do sistema (SO)\n');
            b = calcula_vetor_de_aprox_SO(A,f,s,np_comp,...
                sigma);
        end
    end
    alpha = calcula_coefs_aproximacao(M,Minv,b);
    s = A*alpha; %aproximacao nos pontos de amostragem
    convergiu = false;
    if (parametros.metodo == 1) || (convergiu) || (iter >= niter)
        break
    end
end
e = f - s;
map = zeros(np,1);
for i = 1:np
    %verificar se eh norma 2 ou norma_erro
    map(i,1) = norm(e(i,:),2);
end
erro_L2 = norma_erro_L2(map,np_comp);
erro_Linf = norma_erro_Linf(map);
fprintf(2,'norma do erro = %9.6f (L2) %9.6f (Linf)\n', erro_L2, erro_Linf);
if grava_tex
    fprintf(arq,'\\newcommand{\\PARMerroHilb}{%s}\n',...
        converte_erro_tex(erro_L2));
    fprintf(arq,'\\newcommand{\\PARMerroLinf}{%s}\n',...
        converte_erro_tex(erro_Linf));
    fclose(arq);
    fprintf(2,'arquivo %s fechado\n',nomearq);
end
for j = 1:nf
    if parametros.funcao == 0
        nomef = sprintf('hartley(%d)',j);
    else
        %nomef = sprintf('funcao %d',parametros.funcao);
        nomef = 'funcao';
    end
    
    if plota.componentes
        for k = 1:nb
            figure(4+plota.janela_base);
            if isnan(NP(1))
                plota_funcao_irreg(p,A(:,k)*alpha(k,j),CD,LD,NaN,NaN);
            else
                plota_funcao(p,A(:,k)*alpha(k,j),true,CD,LD,...
                    NP,NaN,NaN);
            end
            define_zlim2(alpha,f,s);
            title(sprintf('elemento escalonado %d da %s %s',...
                k,nomef,estagio),'FontSize',20);
            hold off
            saveas(gcf,sprintf('%s/f%04d_alpha_psi%04d.eps',...
                pasta,j,k));
        end
    end
    
    if plota.coeficientes
        %plota os coeficientes alpha na grade de centros
        assert(false,'plota.coeficietes nao funciona')
        figure(14+plota.janela_base);
        C = centros_da_base(base);
        plota_funcao(C,alpha(:,j),true,CD,LD,NP,NaN,NaN);
        title(sprintf('coeficientes para %s %s',...
            nomef,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/alpha%04d.eps',pasta,j));
        end
    end
    
    if plota.aproximacao
        figure(5+plota.janela_base);
        if isnan(NP(1))
            plota_funcao_irreg(p,s(:,j),CD,LD,NaN,NaN);
        else
            plota_funcao(p,s(:,j),true,CD,LD,NP,NaN,NaN);
        end
        title(sprintf('aproximacao da %s %s',...
            nomef,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/s%04d.eps',pasta,j));
        end
    end
    
    if plota.erro
        figure(6+plota.janela_base);
        if isnan(NP(1))
            plota_funcao_irreg(p,e(:,j),CD,LD,NaN,NaN);
        elseif plota.erro_curva_nivel
            plot_contour_funcao(parametros,p,e(:,j),CD,LD,NP);
        else
            emax = max([plota.emax, max(abs(e(:,j)))]);
            plota_funcao(p,e(:,j),false,CD,LD,NP,-emax,emax);
        end
        title(sprintf('erro da %s %s',...
            nomef,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/e%04d.eps',pasta,j));
            
        end
    end
end


if plota.base_dual
    %determina a base lagrangiana/Dual
    fprintf(2,'calculando a base dual \n');
    L = A*Minv;
    
    %intervalo de elementos da a plotar
    ib = escolhe_elemento_central(base,CD);
    fb = ib;
    
    for k = ib:fb
        figure(3+plota.janela_base);
        if isnan(NP(1))
            plota_funcao_irreg(p,L(:,k),CD,LD,NaN,NaN);
        elseif plota.curva_nivel
            plot_contour_funcao(parametros,p,L(:,k),CD,LD,NP);
        else
            plota_funcao(p,L(:,k),true,CD,LD,NP,NaN,NaN);
        end
        title(sprintf('elemento dual %d %s',...
            k,estagio),'FontSize',20);
        hold off
        if plota.salva_eps
            saveas(gcf,sprintf('%s/theta%04d.eps',pasta,k));
        end
        %pause
    end
end

fprintf(2,'--- fim de aproximacao -----------------------\n');

