
% Método de las características
% Este programa no considera la aparición de cavitación

clc
clear all
close all
 
% Datos
Lgp = 3*154.8 ;     % longitud de la galería de presión(m)
Ltf = 1197.25-Lgp ; % longitud de la tubería forzada(m)
Ltdif = 20 ;        % longitud del tubo difusor

L= Lgp + Ltf ;

landa = 0.02 ;      % coeficiente de fricción del conducto
D = 4.4 ;           % diámetro del conducto (m)
v_0 = 1 ;           % m/s
c_0 = 1300 ;        % m/s

N = 50 ;
inc_t = 1.0/N ;     % Incremento de tiempo
tfin = 50 ;
t_p = (0:1.0/N:tfin);
lt=length(t_p);     % = tfin*N + 1

x_p = (0:1.0/N:1);  
lx=length(x_p);     % = N + 1

% Condiciones de iniciales para t=0 (primera fila) y para toda x (todas las columnas)

p_p(1,1:lx) = 0 ;
v_p(1,1:lx) = 1 ;

% A continuación, se resuelven el resto de los puntos

A = [1 1 ; 1 -1];

for i = 2:lt
    j=1;
    p_p(i,j) = 0 ;
    v_p(i,j) = v_p(i-1,j+1) - p_p(i-1,j+1) ;
    j=lx;
    v_p(i,j) = 0 ;
    p_p(i,j) = v_p(i-1,j-1) + p_p(i-1,j-1) ;
    for j = 2:lx-1
        b(1,1) = v_p(i-1,j-1) + p_p(i-1,j-1);
        b(2,1) = v_p(i-1,j+1) - p_p(i-1,j+1);
        x=A\b;
        v_p(i,j) = x(1);
        p_p(i,j) = x(2);
    end
end
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
% % Representación gráfica (animación)
% 
%  m = size(v_p,1) ;
%  figure(1)
%  clf
% % hv = plot(x,v(1:n),'EraseMode','background');
%  hp = plot(x_p,p_p(1,:),'EraseMode','xor');
%  hold on
%  hv = plot(x_p,v_p(1,:),'r','EraseMode','xor');
%  axis([0 1 -1.2 1.2])
% % axis square
% % box on
%  set(gca,'xtick',[0 0.2 0.4 0.6 0.8 1])
%  set(gca,'ytick',[-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 1.2])
% % grid
%  xlabel('{\it x/L}','Fontsize',15)
%  ylabel('{\it v/v}_o ,   {\it p/p}_o','Fontsize',15)
%  title('Presión y velocidad a lo largo del conducto')
% 
%  for j=2:m
%    pause(0.01);
%    % set(hv,'XData',x,'YData',v( n*(j-1)+1 : n*j ))
%    set(hp,'YData',p_p(j,:))
%    set(hv,'YData',v_p(j,:))
%    drawnow
%  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(2)
  axis([0 1 -2 2])
%  plot(t_p,p_p(:,12)','r')
%  hold on 
 plot(t_p,p_p(:,51)','r')
 xlabel('{\it t/t}_i_v','Fontsize',15)
 ylabel('{\it p/p}_o','Fontsize',15)
 hold on
 

 