#dimensionamentos
n = 100;
qa = ones(n,2);
q = qa;
f = zeros(n+1,2);
x = zeros(n,1);
#dados
t = 0;
dt = 0.001;
nt = 1000;
dx = 1./n;
ui = 0.;
hi = 1.;
g=1.;
omega=2*pi/0.2;
for i=1:n;
 x(i)=(i-0.5)*dx;
 eta=0.1*exp(-(x(i)-0.5)^2/0.002);
 q(i,1)=hi+eta;
# q(i,2)=hi*ui;
 q(i,2)=hi*ui+eta*(ui+sqrt(g*hi));
# if (x(i)<0.3) q(i,1)=12.;
# else q(i,1)=1.; endif
# q(i,2)=0.;
endfor
#codigo

plot(x,q(:,1),"-"); pause
for m=1:nt
# atualizamos
        qa = q;
# calculo dos fluxos
for i=2:n
xx = 0.5*(x(i-1)+x(i));
fezq(1)=q(i-1,2);
fezq(2)=q(i-1,2)^2/q(i-1,1)+0.5*g*q(i-1,1)^2;
fder(1)=q(i,2);
fder(2)=q(i,2)^2/q(i,1)+0.5*g*q(i,1)^2;

qint_lw=(q(i-1,:)+q(i,:))/2-dt/2/dx*(fder-fezq);

f(i,1)=qint_lw(2);
f(i,2)=qint_lw(2)^2/qint_lw(1)+0.5*g*qint_lw(1)^2;
endfor
#condicoes de contorno
#ezquerda
#media de Roe
#h0=0.5*(qa(i-1,1)+qa(i,1));
#u0=(qa(i-1,2)/sqrt(qa(i-1,1))+qa(i,2)/sqrt(qa(i,1)))/(sqrt(qa(i-1,1))+sqrt(qa(i,1)));
#
h0=hi;
u0=ui;
c0=sqrt(g*h0);
A(1,1)=0;A(1,2)=1;A(2,1)=-u0*u0+g*h0;A(2,2)=2*u0;
R=[1 1;u0-c0 u0+c0];
D=[u0-c0 0;0 u0+c0];

z11=1/(2*c0)*((u0+c0)*qa(1,1)-qa(1,2));
#fechado
hezq=2*c0*z11/(u0+c0);
q2ezq=0;
#onda
#eta=0.0*sin(omega*t);
#if t>0.2 eta=0; endif
#z2ezq=1/2*(h0+2*eta);
#zezq=[z11,z2ezq]';
#qezq=R*zezq;
qezq=[hezq; q2ezq];
f(1,1)=qezq(2);
f(1,2)=qezq(2)^2/qezq(1)+0.5*g*qezq(1)^2;
#f(1,:)=A*qezq;
#direita
z2n=1/(2*c0)*(-(u0-c0)*qa(n,1)+qa(n,2));
#lado direito fechado
hder=2*c0*z2n/(c0-u0);
#lado direito aberto
#z1n=0;
#qder=R*[z1n,z2n]';
qder=[hder; 0];
f(n+1,1)=qder(2);
f(n+1,2)=qder(2)^2/qder(1)+0.5*g*qder(1)^2;
#f(n+1,:)=A*qder;

# calculo do q
for i=1:n
	q(i,:)=qa(i,:)-dt/dx*(f(i+1,:)-f(i,:));
endfor
# graficamos
plot(x,q(:,1),"-+"); pause
t=t+dt;
endfor