#dimensioning
n = 100;
qa = ones(n,2);
q = qa;
f = zeros(n+1,2);
x = zeros(n,1);
#dados
t = 0;
dt = 0.00125;
nt = 1000;
dx = 1./n;
ui = 0.;
hi = 0.;
g=1.;
omega=2*pi/0.2;
for i=1:n;
 x(i)=(i-0.4)*dx;
 eta=2.*exp(-(x(i)-0.4)^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
#code----------------------------------
plot(x,q(:,1),"-"); pause
for m=1:nt
# update
        qa = q;
# fluxes
for i=2:n
xx = 0.5*(x(i-1)+x(i));
R=zeros(2,2);
h0=2;
u0=0;
c0=sqrt(g*h0);
A=[0 1;-u0*u0+g*h0 2*u0];
R=[1 1;u0-c0 u0+c0];
D=[u0-c0 0;0 u0+c0];
Dmenos=min(D,0);
Imenos=abs(sign(Dmenos));
Mmenos=R*Imenos*inv(R);
qint_godunov=qa(i-1,:)'+Mmenos*(qa(i,:)'-qa(i-1,:)');
f(i,:)=(A*qint_godunov);
endfor
#boundary conditions
w11=1/(2*c0)*((u0+c0)*qa(1,1)-qa(1,2));
#closed left
hezq=2*c0*w11/(u0+c0);
q2ezq=0;
qezq=[hezq; q2ezq];
f(1,:)=A*qezq;

##closed right
#w2n=1/(2*c0)*(-(u0-c0)*qa(n,1)+qa(n,2));
#hder=2*c0*w2n/(c0-u0);
#qder=[hder;0];

#open right
w2n=1/(2*c0)*(-(u0-c0)*qa(n,1)+qa(n,2));
hder=w2n;
qder=[hder;(u0+c0)*hder];
f(n+1,:)=A*qder;

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