; This program shows how a triangle wave on a string evolves in time.
; RWB, October 18, 2000
PRO square
   char1=''; input buffer
   ncmax=100; number of coefficients
   coef=fltarr(ncmax); array for coefficients, initialized to zero.
   ;*** Calculate the coefficients for the sine-wave series.
   ;*** This calculation is for a triangle wave.
   for n = 1,ncmax-1,2 do begin
      coef(n)=4./n/!pi
   ENDFOR
   ;*** Now set up the calculation of the string displacement.
   npt=100; number of points on the string.
   fltot=14.; length of string, in meters.
   v=14.; speed of waves, in m/s.
   k0=2.*!pi/(2.*fltot); wave number for lowest mode
   omega0=k0*v; angular frequency for the lowest mode
   x=findgen(npt)*fltot/npt; array of x values.
   nc=10; number of coefficients to use in reconstructing
   ;*** In this while loop, the waveform will be calculated using
   ;*** nc coefficients.  The first time through, nc = 10; after that,
   ;*** nc is entered from the console.
   WHILE nc ne 0 DO BEGIN
      f=fltarr(npt)
      FOR ic=0,nc-1 DO BEGIN
	 f=f+coef(ic)*sin(ic*k0*x)
      ENDFOR
      plot,x,f
      ncsave=nc; save last value for next phase.
      read,'maximum order = '+string(nc-1)+' (return to exit loop) ',char1
      IF strlen(char1) GT 0 THEN nc = fix(char1)+1
      IF strlen(char1) EQ 0 THEN nc = 0
      nc=min([nc,ncmax])
   ENDWHILE
   ;*** Now the number of terms has been decided.  Calculate the time
   ;***  evolution.
   tmax=10.
   dt=.0205
   istep=1; step flag; 1 means one step at a time.
   plot,x,f,yrange=[-1,1]
   oplot,[0.,fltot],[0.,0.],linestyle=5
   FOR t=0.,tmax,dt DO BEGIN
      y=fltarr(npt)
      FOR ic=0,ncsave-1 DO BEGIN
	 y=y+coef(ic)*sin(ic*k0*x)*cos(ic*omega0*t)
      ENDFOR
      oplot,x,y
      IF istep eq 1 THEN BEGIN
	 READ,'Return for one step, 0 to finish fast ',char1
	 IF strlen(char1) ne 0 THEN istep=0
      ENDIF ELSE BEGIN
	 ;print,t
      ENDELSE
   ENDFOR

   SAVE
END