;demo_wave_sep
;Written by Thomas Wiegelmann
;10.09.2008
;Email: wiegelmann@mps.mpg.de
;Solve wave-equation by separation-ansatz
;Here Boundary condition: Phi(0,t)=Phi(Lx,t)=0
;Intial profile Phi(x,0) approximated by
;Fourier-series (sin-functions)
;See PDE-lecture for details
;
FUNCTION rho_0, x
;Gauss profile
x0=5.0
rho_bg=0.0
rho0=rho_bg+1.0*exp(-(x-x0)^2/2)
return,rho0
END
;
FUNCTION rho_1, x
;Rectangular profile
;Square wave test
x0=5
rho0=0.0*x+0.0
u=where((x ge x0-1.) and (x le x0+1.))
rho0[u]=1.1
return,rho0
END

;PRO demo_wave_sep
;Solve wave-equation by separation
;Some definitions
Lx=10.0 & nx=101
tend=120.0
nt=201
tstep=tend/(nt-1)

c=0.5; wave velocty

;Define grid and time discretisation
x=Lx*findgen(nx)/(nx-1) ; spatial grid
dx=Lx/(nx-1)
tvek=tend*findgen(nt)/(nt-1)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Compute coefficients ;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
prof=1
print,'Wave equation solved by variable separation'
print,' '
print,'Initial profile: 1=Gauss,  2=Square box'
read,prof
fx=rho_0(x)
if prof eq 2 then $
fx=rho_1(x)
;
boundary=1
print,'Which boundary condition ?'
print,'1 : Dirichlet, Phi(0,t) = Phi(Lx,t)=0'
print,'2 : Von Neumann, d Phi(0,t)/dx = d Phi(Lx,t)/dx=0'
print,'3 : Periodic'
read, boundary
;
maxi=1.2*max(fx)
plot,x,fx,yrange=[-maxi,maxi],linestyle=1,xtitle='x',$
ytitle='Phi(x,0)', charsize=1.5,title='Intitial condition'
xyouts,1,-0.8,'dotted line = Original',size=2
print,'How many Fourier coefficients?'
fouri=-1
read, fouri
if fouri le 0 then fouri=nx-1
maxi=1.2*max(abs(fx)); plotting range
ak_vec=fltarr(fouri+1)
fx0=fltarr(nx)
if boundary eq 1 then begin
for k=1,fouri do begin
dummy=(2.0/Lx)*dx*total(fx*sin(k*!pi*x/Lx))
ak_vec[k]=dummy
fx0=fx0+ak_vec[k]*sin(k*!pi*x/Lx)
endfor
endif; boundary 1
if (boundary eq 2) or (boundary eq 3) then begin
ak_vec[0]=(1.0/Lx)*dx*total(fx)
fx0=ak_vec[0]
for k=1, fouri do begin
dummy=(2.0/Lx)*dx*total(fx*cos(k*!pi*x/Lx))
ak_vec[k]=dummy
fx0=fx0+ak_vec[k]*cos(k*!pi*x/Lx)
endfor
endif; boundary 2, 3
;wait,2
oplot,x,fx0,linestyle=0
xyouts,1,-1.2,'solid line = Fourier Series',size=2
;wait,5
dummy='y'
print,'Satisfied? (y or n)'
read,dummy
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Compute solution of wave equation ;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
if dummy eq 'y' then begin
for it=0,nt-1 do begin
t=tvek[it]
phi=fltarr(nx)
for k=0, n_elements(ak_vec)-1 do begin
CASE boundary of
   1: phi=phi+ak_vec[k]*sin(k*!pi*x/Lx)*cos(c*k*!pi*t/Lx)
   2: phi=phi+ak_vec[k]*cos(k*!pi*x/Lx)*cos(c*k*!pi*t/Lx)
   3: phi=phi+ak_vec[k]*cos(k*!pi*(x-c*t)/Lx)

ENDCASE
endfor
plot,x,phi,yrange=[-maxi,maxi],xtitle='x',$
ytitle='Phi(x,t)', charsize=1.5
oplot,x,fx,linestyle=1
xyouts,1,1.2,'t='+strcompress(string(t)),size=2
wait,0.1

endfor
endif

END