PROGRAM d12r3(input,output);
(* driver for routine REALFT *)
LABEL 1,99;
CONST
   eps=1.0e-3;
   np=32;
   np2=34;   (* np+2 *)
   width=50.0;
   pi=3.1415926;
TYPE
   gldarray = ARRAY [1..np2] OF real;
VAR
   big,per,scal,small : real;
   i,j,n,nlim : integer;
   data,size : gldarray;

(*$I MODFILE.PAS *)
(*$I FOUR1.PAS *)

(*$I REALFT.PAS *)

BEGIN
   n := np DIV 2;
1:   writeln('Period of sinusoid in channels (2-',np:2,')');
   readln(per);
   IF (per <= 0.0) THEN GOTO 99;
   FOR i := 1 to np DO BEGIN
      data[i] := cos(2.0*pi*(i-1)/per)
   END;
   realft(data,n,+1);
   big := -1.0e10;
   FOR i := 1 to n DO BEGIN
      size[i] := sqrt(sqr(data[2*i-1])+sqr(data[2*i]));
      IF (size[i] > big) THEN big := size[i]
   END;
   size[1] := data[1];
   IF (size[1] > big) THEN big := size[1];
   scal := width/big;
   FOR i := 1 to n DO BEGIN
      nlim := round(scal*size[i]+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   writeln('press RETURN to continue ...');
   readln;
   realft(data,n,-1);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF (data[i] < small) THEN small := data[i];
      IF (data[i] > big) THEN big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   GOTO 1;
99:
END.