; IDL Version 6.0, Mac OS X (darwin ppc m32)
; Journal File for bob@Bobs-Computer.local
; Working directory: /Users/Teach/AST802
; Date: Tue Jan 13 13:11:58 2004
 
;Spatial Grid
nx=100
x=findgen(nx)
help,x
dx=x(1)-x(0)
print,dx
;      1.00000
;Time step
dt=0.4

;Initial condition
uorig=fltarr(nx)
amp=1
uorig=amp*sin(2.*!pi*x/50.)+1.
plot,uorig

;Velocity array
nt=50
u=fltarr(nx,nt)
;Initial condition everywhere
u(*,0)=uorig
;Boundary condition
u(0,*)=u(0,0)

;Implicit Scheme (Note there was a sign error in front of the 2nd term 
;in the square root which made it negative.  It now works
for t=1,49 do begin $
for j=1,99 do u[j,t]=-(dx/dt)+sqrt((dx/dt)^2 + $
((2*dx/dt)*u(j,t-1)+u(j-1,t)^2))
plot,u(*,0),line=1
oplot,u(*,20),psym=-1
oplot,u(*,40),psym=-1

;Shift operator
y=[1,2,3,4,5]
print,shift(y,+1)
;       5       1       2       3       4
print,shift(y,-1)
;       2       3       4       5       1

;Explicit Scheme
;Upwind
u(*,0)=uorig
for n=1,nt-1 do u(*,n)=u(*,n-1)-(dt/(2*dx))*(u[*,n-1]^2-shift(u[*,n-1]^2,+1))
plot,u(*,0),line=1
oplot,u(*,20),psym=-1
oplot,u(*,40),psym=-1

;Centered
u(*,0)=uorig
for n=1,nt-1 do u(*,n)=u(*,n-1) - (dt/(2*dx))*(shift(u[*,n-1]^2,-1)- $
shift(u[*,n-1]^2,+1))
plot,u(*,10),psym=-1
oplot,u(*,5),psym=-4
oplot,u(*,0),line=1

;Downwind
u(*,0)=uorig
for n=1,nt-1 do u(*,n)=u(*,n-1) - (dt/(2*dx))*(shift(u[*,n-1]^2,-1)- $
shift(u[*,n-1]^2,0))
plot,u(*,15),psym=-1
oplot,u(*,5),line=2
oplot,u(*,0),line=1