% WENO para IC de grado 3 para F y 2 para f
%       preliipmcwenou(xx,ff,nivells)
%       Entradas:
%           x       = nodos de la funcion inicial
%           f       = f(x)
%           nivells = niveles de multirresolucion 
%       Salidas
%           Plot de la reconstruccion(-), funcion en el nivel mas grosero(o), funcion original con mas nodos (--), diferencia original y aproximacion(:) 
%
% Autor: José Jaime Noguera.  +INFO: www.uv.es/anims/noguera

function preliipmcwenou(xx,ff,nivells)
grau=3; %para F

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);
r=grau;

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 optimos
c1=3/16;
c2=10/16;
c3=3/16;
eps=10^(-6);
    
for k=2:nivells
  [xn,xm]=size(x);
  nn=log(xm)/log(2);
  [m1,n1]=size(x(k-1,1:2^(nn-(nivells-k+1))));
   
   for i=3:n1-2
      %indicador de suavidad de Jiang y Shu 
      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/8)*fr(k-1,i-2)+(1/2)*fr(k-1,i-1)+(5/8)*fr(k-1,i))+w2*((1/8)*fr(k-1,i-1)+fr(k-1,i)-(1/8)*fr(k-1,i+1))+w3*((11/8)*fr(k-1,i)-(1/2)*fr(k-1,i+1)+(1/8)*fr(k-1,i+2));     
      fr(k,2*i)=2*fr(k-1,i)-fr(k,2*i-1);
   end
   
   %fronteras
   for i=1:2
      bizq=coeficientesipmc(0,r-1);
      fr(k,2*i-1)=bizq*fr(k-1,i:i+r-1)';
      fr(k,2*i)=2*fr(k-1,i)-fr(k,2*i-1);
   end
   for i=n1-1:n1
      bder=coeficientesipmc(-(r-1),r-1);
      fr(k,2*i-1)=bder*fr(k-1,i-r+1:i)';
      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

