;lecture_advection_draft
;Written by Thomas Wiegelmann
;11.09.2008
;Email: wiegelmann@mps.mpg.de
;This program is a draft for use in PDE-exercise III

;
FUNCTION rho_0, x
;Gauss profile
x0=3.0
rho_bg=0.1
rho0=rho_bg+1.0*exp(-(x-x0)^2/2)
return,rho0
END
;
FUNCTION rho_1, x
;Rectangular profile
;Square wave test for advection equation
x0=5
rho0=0.0*x+0.1
u=where((x ge x0-1.) and (x le x0+1.))
rho0[u]=1.1
return,rho0
END
;
FUNCTION compute_F, utest
;Flux F for Advection equation
v0=0.1
F0=v0*utest
return,F0
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Main program
prof=1
print,'Initial profile: 1=Gauss,  2=Square box'
read,prof
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Initial conditions
nx=100
xmax=10.0
xmin=0.0
x=xmax*findgen(nx)/(nx)-xmin
rho0=rho_0(x) ; Gauss-shape
if prof eq 2 then $
rho0=rho_1(x) ; Square box
dt=0.5; iteration time step; vary dt for experiments!
;Compute the Courant-number and choose time-step dt accordingly
dx=xmax/nx; equidistant grid
rho=rho0 & rho1=rho0
;
for it=0, 200 do begin
t=it*dt
F=compute_F(rho)
rho1=0.5*(shift(rho,1)+shift(rho,-1))-(dt/(2*dx))*(shift(F,-1)-shift(F,1))
;Update density
rho=rho1
;plot solution
plot,x,rho,xstyle=1,ystyle=1,xrange=[xmin,xmax],$
yrange=[-0.2,1.2],xtitle='x',ytitle='rho(x)',thick=2,$
charsize=1.5,title='Lax-Method'

oplot,x,rho0,linestyle=1
xyouts,3,-0.05,'t='+strcompress(string(t)),size=2
;

wait,0.05
endfor
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




end