% Multirresolución WENO para grado  r=3  para F (grado 2 para f) utilizando la recostrucción IAC e indicadores de Jiang y Shu
%       capiaiamcwenou(xx,ff,nivells)
%       Entradas:
%           xx      = nodos de la funcion inicial
%           ff      = f(x)
%           nivells = niveles de multiresolucion 
%       Salidas:
%          plot de la funcion original(--), de la recosntruccion(-),de los nodos del nivel mas grosero (o) y de la diferencia entre la original y la reconstruccion (:)


function capiaiamcwenou(xx,ff,nivells)

nivells=nivells+1;
[xn,xm]=size(xx);
n=log(xm)/log(2);
f=zeros(nivells,xm);
x=zeros(nivells,xm);
f(nivells,1:xm)=ff(1:xm);
x(nivells,1:xm)=xx(1:xm);

%BAJAMOS NIVELES
for k=nivells:-1:2   
   for i=1:2^(n-(nivells-k+1))
      f(k-1,i)=(1/2)*(f(k,2*i-1)+f(k,2*i));
      x(k-1,i)=(x(k,2*i)+x(k,2*i-1))/2;
   end 
end
    
  fr=zeros(nivells,xm);
  fr(1,1:2^(n-(nivells-1)))=f(1,1:2^(n-(nivells-1)));

%Pesos óptimos    
c1=1/2;
c2=0;
c3=1/2;

eps=10^(-6);

 
%SUBIMOS NIVELES
for k=2:nivells
   [m1,n1]=size(x(k-1,1:2^(n-(nivells-k+1))));  
   for i=3:n1-2
      esquerra=(13/12)*(fr(k-1,i)-2*fr(k-1,i-1)+fr(k-1,i-2))^2 + (1/4)*(3*fr(k-1,i)-4*fr(k-1,i-1)+fr(k-1,i-2))^2;
      centre=(13/12)*(fr(k-1,i+1)-2*fr(k-1,i)+fr(k-1,i-1))^2 + (1/4)*(fr(k-1,i+1)-fr(k-1,i-1))^2;
      dreta=(13/12)*(fr(k-1,i+2)-2*fr(k-1,i+1)+fr(k-1,i))^2 + (1/4)*(-fr(k-1,i+2)+4*fr(k-1,i+1)-3*fr(k-1,i))^2;
      
      alpha1=c1/(eps+esquerra);
      alpha2=c2/(eps+centre);
      alpha3=c3/(eps+dreta);
      
      w1=alpha1/(alpha1+alpha2+alpha3);
      w2=alpha2/(alpha1+alpha2+alpha3);
      w3=alpha3/(alpha1+alpha2+alpha3);
      
      fr(k,2*i-1)=w1*((1/10)*fr(k-1,i-2)+(1/20)*fr(k-1,i-1)+(17/20)*fr(k-1,i))+w2*((1/8)*fr(k-1,i-1)+fr(k-1,i)-(1/8)*fr(k-1,i+1))+w3*((23/20)*fr(k-1,i)-(1/20)*fr(k-1,i+1)-(1/10)*fr(k-1,i+2));      
      fr(k,2*i)=2*fr(k-1,i)-fr(k,2*i-1);
   end
     
  %fronteras 
  m=4; % numero de nodos del stencil para F
  nn=3;   
  for i=1:2
     fr(k,2*i-1)=penomcrn(0,i,nn,k,m,fr);
     fr(k,2*i)=2*fr(k-1,i)-fr(k,2*i-1);   
  end   
  for i=n1-1:n1
     fr(k,2*i-1)=penomcrn(m,i,nn,k,m,fr);
     fr(k,2*i)=2*fr(k-1,i)-fr(k,2*i-1);
  end  
  %fin fronteras 

end   
  
[xor,yor]=funtesi(256,0); 

hold on
plot(xor, yor,'--k','LineWidth',2)                                           %original discretizada con 257 pts 
plot(x(1,1:2^(n-(nivells-2+1))),f(1,1:2^(n-(nivells-2+1))),'ko')             %original en el nivel mas grosero
plot(xx,(ff-fr(nivells,1:xm))-(0.1)*ones(1,xm),':k','LineWidth',3)           %original-reconstruccion
plot(xx,fr(nivells,1:xm),'-k','LineWidth',3)                                 %reconstruccion
axis([0 1 -0.2 1.1 ]);
set(gca,'Xtick',[0 .5 1], 'FontSize',15)
set(gca,'Ytick',[0 .5 1], 'FontSize',15)
box
hold off
