% Multirresolución mediante WENO para grado r=3 con reconstrucción AP e indicadores de Jiang y Shu. 
%       capiamincvpwenoultim(xx,ff,nivells)
%       Entradas:
%           xx      = nodos de la funcion inicial
%           ff      = f(x)
%           nivells = niveles de multiresolucion 
%       Salidas:
%           Plot de la reconstruccion(-), funcion original(o), funcion original con mas nodos (--), diferencia original y aproximacion(:) 
%
% Autor: José Jaime Noguera.   +INFO: www.uv.es/anims/noguera

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

%Discretizar por valores puntuales
for k=nivells:-1:2
   for i=1:1+2^(n-nivells+k-1)
      f(k-1,i)=f(k,2*i-1);
      x(k-1,i)=x(k,2*i-1);
   end 
end 
fr=zeros(nivells,xm); %fr significa la reconstruccion
fr(1,1:2^(n-nivells+1)+1)=f(1,1:2^(n-nivells+1)+1);

%pesos optimos
c1=405/1816;
c2=653/908;
c3=405/1816;
    
eps=10^(-6);
    
for k=2:nivells
   ultim=1+2^(n-(nivells-k+1));
   
   for i=4:ultim-2
      %indicadores de suavidad de jiang y shu
      esquerra=(1/13)*(11*fr(k-1,i)-18*fr(k-1,i-1)+9*fr(k-1,i-2)-2*fr(k-1,i-3))^2 + (13/12)*(2*fr(k-1,i)-5*fr(k-1,i-1)+4*fr(k-1,i-2)-fr(k-1,i-3))^2+(1043/960)*(fr(k-1,i)-3*fr(k-1,i-1)+3*fr(k-1,i-2)-fr(k-1,i-3))^2+(1/72)*(11*fr(k-1,i)-18*fr(k-1,i-1)+9*fr(k-1,i-2)-2*fr(k-1,i-3))*(fr(k-1,i)-3*fr(k-1,i-1)+3*fr(k-1,i-2)-fr(k-1,i-3));
      centre=(1/36)*(2*fr(k-1,i+1)+3*fr(k-1,i)-6*fr(k-1,i-1)+fr(k-1,i-2))^2 +(13/12)*(fr(k-1,i+1)-2*fr(k-1,i)+fr(k-1,i-1))^2+(1043/960)*(fr(k-1,i+1)-3*fr(k-1,i)+3*fr(k-1,i-1)-fr(k-1,i-2))^2+(1/72)*(2*fr(k-1,i+1)+3*fr(k-1,i)-6*fr(k-1,i-1)+fr(k-1,i-2))*(fr(k-1,i+1)-3*fr(k-1,i)+3*fr(k-1,i-1)-fr(k-1,i-2));
      dreta=(1/36)*(-fr(k-1,i+2)+6*fr(k-1,i+1)-3*fr(k-1,i)-2*fr(k-1,i-1))^2+(13/12)*(fr(k-1,i+2)-2*fr(k-1,i+1)+fr(k-1,i))^2+(1043/960)*(fr(k-1,i+2)-3*fr(k-1,i+1)+3*fr(k-1,i)-fr(k-1,i-1))^2+(1/72)*(-fr(k-1,i+2)+6*fr(k-1,i+1)-3*fr(k-1,i)-2*fr(k-1,i-1))*(fr(k-1,i+2)-3*fr(k-1,i+1)+3*fr(k-1,i)-fr(k-1,i-1));
      
      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);
      
      %aproximacion weno
      fr(k,2*i-2)=w1*(-(9/80)*fr(k-1,i-3)+(17/80)*fr(k-1,i-2)+(33/80)*fr(k-1,i-1)+(39/80)*fr(k-1,i))+w2*(-(1/16)*fr(k-1,i-2)+(9/16)*fr(k-1,i-1)+(9/16)*fr(k-1,i)+(-1/16)*fr(k-1,i+1))+w3*((39/80)*fr(k-1,i-1)+(33/80)*fr(k-1,i)+(17/80)*fr(k-1,i+1)-(9/80)*fr(k-1,i+2));
      fr(k,2*i-1)=fr(k-1,i);
   end
      
   %fronteras 
    fr(k,1)=fr(k-1,1); 
    fr(k,2)=(39/80)*fr(k-1,1)+(33/80)*fr(k-1,2)+(17/80)*fr(k-1,3)-(9/80)*fr(k-1,4);
    fr(k,3)=fr(k-1,2);
    fr(k,4)=-(1/16)*fr(k-1,1)+(9/16)*fr(k-1,2)+(9/16)*fr(k-1,3)-(1/16)*fr(k-1,4);
    fr(k,5)=fr(k-1,3);
    
    fr(k,2*ultim-4)=-(1/16)*fr(k-1,ultim-3)+(9/16)*fr(k-1,ultim-2)+(9/16)*fr(k-1,ultim-1)-(1/16)*fr(k-1,ultim); %i=ultim-1
    fr(k,2*ultim-3)=fr(k-1,ultim-1);
    fr(k,2*ultim-2)=(-9/80)*fr(k-1,ultim-3)+(17/80)*fr(k-1,ultim-2)+(33/80)*fr(k-1,ultim-1)+(39/80)*fr(k-1,ultim); %i=ultim
    fr(k,2*ultim-1)=fr(k-1,ultim);
  % fin fronteras 

end   

[xor,yor]=funtreballtesi(257,0); 

hold on
plot(xor, yor,'--k','LineWidth',2)                                           %original discretizada con 257 pts 
plot(x(1,1:1+2^(n-(nivells-2+1))),f(1,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


