close,/all
lambdau = 0.210
;
; create the y scale (500 points in the interval [-.84,.84])
npts = 501
y=findgen(npts)/float(npts-1)
y=y*2.0*0.84
y=y-0.84
;
; read harmonics from second column in ex3_h.dat
;
a1 = rascii('ex3_h.dat')
hmax = max(a1[0,*])
h = Reform(a1[1,*])
itmp = Dialog_Message(/Info,'Harmonics: '+Vect2String(h))
phase = -!pi + 2.0D0*!pi*y/lambdau
help,phase

; first harmonic
b = h[0] * cos (phase)
help,b
; the others
FOR i=1,hmax-1 DO BEGIN
  n = i+1
  b = b + h[i] * cos(phase*n)
ENDFOR
;itmp = Dialog_Message(/Info,'Harmonics: '+Vect2String(y))
xplot,y,b