PROGRAM synthesize
CALL plotf(0,0.5,0.5,1)
CALL plotf(0,0.5,0,0.5)
CALL plotf(0.5,1,0.5,1)
CALL plotf(0.5,1,0,0.5)
END

SUB plotf(xmin,xmax,ymin,ymax)
    OPEN #1: screen xmin,xmax,ymin,ymax
    SET WINDOW -4,4,-2,2
    PLOT LINES: -pi,0;pi,0
    PLOT LINES: 0,-1.5;0,1.5
    INPUT prompt "number of modes = ": N
    SET COLOR "red"
    CALL fourier(N)
    CLOSE #1
END SUB

SUB fourier(N)
    ! compute Fourier series and plot function
    DIM a(0 to 1000),b(1000)
    CALL coefficients(N,a(),b())
    LET nplot = 100
    LET t = -pi
    LET dt = pi/100
    DO while t <= pi
       LET f = a(0)/2
       FOR k = 1 to N
           IF a(k) <> 0 then LET f = f + a(k)*cos(k*t)
           IF b(k) <> 0 then LET f = f + b(k)*sin(k*t)
       NEXT k
       PLOT LINES: t,f;
       LET t = t + dt
    LOOP
END SUB

SUB coefficients(N,a(),b())
    ! generate Fourier coefficients for special case 
    LET a(0) = 0
    FOR k = 1 to N
        LET a(k) = 0
        IF mod(k,2) <> 0 then
           LET b(k) = 2/(k*pi)
        ELSE
           LET b(k) = 0
        END IF
    NEXT k
END SUB