; 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.
;Sine wave
uorig=amp*sin(2.*!pi*x/50.)+1.
;Shock wave
uorig(0:50)=amp
uorig(51:*)=0.
;Rarefaction wave
;uorig(0:50)=0.
;uorig(51:*)=amp
;plot,uorig,line=1
tmp=' '
;read,tmp

;Velocity array
nt=50
u=fltarr(4,nx,nt)
;Initial condition everywhere
u(0,*,0)=uorig
;Boundary condition
u(0,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
;Backwards time difference
;(u[j,n+1]-u[j,n])/dt + (u[j,n+1]^2-u[j-1,n+1]^2)/2dx = 0
for t=1,49 do begin $
for j=1,99 do u[0,j,t]=-(dx/dt)+sqrt((dx/dt)^2 + $
((2*dx/dt)*u[0,j,t-1]+u[0,j-1,t]^2))
plot,u(0,*,40),line=0,psym=-1,title='Implicit',ystyle=3
oplot,shift(uorig,amp*40*dt/2),line=1
;plot,u(0,*,0),line=1,yrange=[0,3]
;oplot,u(0,*,20),psym=-1
;oplot,u(0,*,40),psym=-1
read,tmp

;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
;Forward time difference, Upwind space difference (assuming u>0)
;(u[j,n+1]-u[j,n])/dt+((u[j,n]^2-u[j-1,n]^2)/dx = 0
u(1,*,0)=uorig
for n=1,nt-1 do u(1,*,n)=u(1,*,n-1)-(dt/(2*dx))*(u[1,*,n-1]^2-shift(u[1,*,n-1]^2,+1))
;plot,u(1,*,0),line=1,title='Explicit: upwind',ystyle=3
;oplot,u(1,*,20),psym=-1
oplot,u(1,*,40),line=2,psym=-4
oplot,shift(uorig,amp*40*dt/2),line=1
read,tmp

;Centered
;Forward time difference, centered space difference
;(u[j,n+1]-u[j,n])/dt+((u[j+1,n]^2-u[j-1,n]^2)/dx = 0
u(2,*,0)=uorig
for n=1,nt-1 do u(2,*,n)=u(2,*,n-1) - (dt/(2*dx))*(shift(u[2,*,n-1]^2,-1)- $
shift(u[2,*,n-1]^2,+1))
plot,u(2,*,10),psym=-1,title='Explicit: centered'
;oplot,u(2,*,5),psym=-4
;oplot,u(2,*,10),line=4,psym=-4
;oplot,u(2,*,0),line=1
oplot,shift(uorig,10*dt/2),lin=1
read,tmp

;Downwind
;Forward time difference, downwind space difference (assuming u>0)
;(u[j,n+1]-u[j,n])/dt+((u[j+1,n]^2-u[j,n]^2)/dx = 0
u(3,*,0)=uorig
for n=1,nt-1 do u(3,*,n)=u(3,*,n-1) - (dt/(2*dx))*(shift(u[3,*,n-1]^2,-1)- $
shift(u[3,*,n-1]^2,0))
plot,u(3,*,15),psym=-1,title='Explicit: downwind',ystyle=3
;oplot,u(3,*,5),line=2
;oplot,u(3,*,15),line=5,psym=-5
;plot,u(3,*,0),line=1
oplot,shift(uorig,15*dt/2),lin=1

;label,x=.5,y=.9,'IC',line=1
;label,'Backward Time (implicit)',psym=-1
;label,'Upwind (explicit)',lin=2,psym=-2
;label,'Centered (explicit)',lin=4,psym=-4
;label,'Downwind (explicit)',lin=5,psym=-5