#
# fftl02.rb
#
# $Id: fftl02.rb,v 1.2 2000/11/27 01:53:56 keiko Exp $
#

require "narray"
require "numru/dcl"

include NumRu
include Math


#-----------------------------------------------------------------------
module Tfftp
#
#
# PURPOSE                TO DEMONSTRATE THE USE OF FFTPACK, AND TO
#                        TEST THE PERFORMANCE OF FFTPACK ON ONE
#                        WELL-CONDITIONED PROBLEM.
#
# USAGE                  CALL TFFTPK (IERROR)
#
# ARGUMENTS
#
# ON OUTPUT              IERROR
#                          INTEGER VARIABLE SET TO ZERO IF FFTPACK
#                          CORRECTLY SOLVED THE TEST PROBLEM, AND
#                          ONE IF FFTPACK FAILED.
#
# I/O                    IF THE TEST SUCCEEDS(FAILS), THE MESSAGE
#
#                           FFTPACK TEST SUCCESSFUL (UNSUCCESSFUL)
#
#                        IS WRITTEN ON UNIT 6. IN THE CASE OF FAILURE,
#                        ADDITIONAL MESSAGES ARE WRITTEN IDENTIFYING THE
#                        FAILURE MORE EXPLICITLY.
#
# PRECISION              SINGLE
#
# REQUIRED LIBRARY       NONE
# FILES
#
# LANGUAGE               FORTRAN
#
# HISTORY                WRITTEN BY MEMBERS OF THE SCIENTIFIC
#                        COMPUTING DIVISION OF NCAR,
#                        BOULDER COLORADO.
#
# ALGORITHM              FOR EACH OF THE ROUTINES, RFFTF, RFFTB, EZFFTF,
#                        AND EZFFTB IN THE FFTPACK PACKAGE A SIMILAR
#                        TEST IS RUN. AN APPROPIATE VECTOR, FOR WHICH
#                        THE EXACT TRANSFORM IS KNOWN IS USED AS THE
#                        INPUT VECTOR. THE ROUTINE IS CALLED TO PERFORM
#                        THE TRANSFORM. THE CALCULATED TRANSFORM VECTOR
#                        IS COMPARED WITH THE EXACT TRANSFORM TO SEE
#                        WHETHER THE PERFORMANCE CRITERION IS MET WITHIN
#                        THE SPECIFED TOLERANCE.
#
#                        FOR RFFTF AND EZFFTF, A REAL VECTOR, THE ELEMEN
#                        WHICH ARE EQUAL TO ONE, IS USED AS INPUT. THE
#                        TRANSFORMED VECTOR HAS THE FIRST ELEMENT EQUAL
#                        TO THE LENGTH OF THE INPUT VECTOR. ALL OTHER
#                        ELEMENTS ARE EQUAL TO ZERO.
#
#                        FOR RFFTB AND EZFFTB, THE INPUT VECTOR HAS FIRS
#                        ELEMENT EQUAL TO ONE AND ALL THE OTHER ELEMENTS
#                        EQUAL TO ZERO. THE TRANSFORMED VECTOR HAS ALL
#                        COMPONENTS EQUAL TO ONE.
#
# PORTABILITY            ANSI STANDARD
#
#
  N = 36
  DIM1 = 2*N+15
  DIM2 = 3*N+15
  ND2 = N/2
  TOL = 0.01

  M998  = " FFTPACK TEST SUCCESSFUL\n"
  M999  = " FFTPACK TEST UNSUCCESSFUL\n"
  M1001 = " IN FFTPACK, ENTRY RFFTF RESULTS IN ERROR\n"
  M1002 = " IN FFTPACK, ENTRY RFFTB RESULTS IN ERROR\n"
  M1003 = " IN FFTPACK, ENTRY EZFFTF RESULTS IN ERROR\n"
  M1004 = " IN FFTPACK, ENTRY EZFFTB RESULTS IN ERROR\n"

#  def stores(x)
#
# FORCES THE ARGUMENT VALUE X TO BE STORED IN MEMORY LOCATION V.
#
#    common /value/ v
#    v=x
#  end
  def trunc(x)
#
# TRUNC IS A PORTABLE FORTRAN FUNCTION WHICH TRUNCATES A VALUE TO THE
# MACHINE SINGLE PRECISION WORD SIZE, REGARDLESS OF WHETHER LONGER
# PRECISION INTERNAL REGISTERS ARE USED FOR FLOATING POINT ARITHMETIC IN
# COMPUTING THE VALUE INPUT TO TRUNC.  THE METHOD USED IS TO FORCE A
# STORE INTO MEMORY BY USING A COMMON BLOCK IN ANOTHER SUBROUTINE.
#
#    common /value/ v
#    stores(x)
#    trunc=v
    x.to_f  #????
  end
#
# STATEMENT FUNCTION SMALL(EX) IS FOR TESTING WHETHER X IS CLOSE TO ZERO
# INDEPENDENTLY OF MACHINE WORD SIZE. SMALL(EX) IS EXACTLY ZERO ONLY IF
# ABS(X) .LT. EPS/TOL, WHERE EPS IS THE MACHINE PRECESION AND TOL IS A
# TOLERANCE FACTOR USED TO CONTROL THE STRICTNESS OF THE TEST.
#
  def small(ex)
    trunc(1.0+TOL*ex.abs)-1.0
  end

  def fft_test
    rldat = NArray.sfloat(N)
#
# CALL INITIALIZATION ROUTINE FOR RFFTF AND RFFTB.
#
    wrfft = DCL::rffti(N)
#
# TEST OF RFFTF.
#
    rldat.fill!(1.0)
#
    result = DCL::rfftf(rldat, wrfft)
#
# TEST RESULTS OF RFFTF
#
    error = (N.to_f - result[0]).abs
    for i in 1..N-1
      error = [error, (result[i]).abs].max
    end
    if (small(error) == 0)
      ier1 = 0
    else
      ier1 = 1
      print M1001
    end
#
# TEST OF RFFTB.
#
    rldat.fill!(0.0)
    rldat[0] = 1.0
#
    result = DCL::rfftb(rldat, wrfft)
#
# TEST RESULTS OF RFFTB
#
    error = 0.0
    for i in 0..N-1
      error = [error, (1.0 - result[i]).abs].max
    end
    if (small(error) == 0)
      ier2 = 0
    else
      ier2 = 1
      print M1002
    end
#
# CALL INITIALIZATION ROUTINE EZFFTI FOR EZFFTF AND EZFFTB
#
    wezfft = DCL::ezffti(N)
#
# TEST OF EZFFTF.
#
    rldat.fill!(1.0)
#
#
    azero, a, b = DCL::ezfftf(rldat, wezfft)
#
# TEST RESULTS OF EZFFTF
#
    error = (azero - 1.0).abs
    for i in 0..ND2-1
      error = [(a[i]).abs + (b[i]).abs, error].max
    end
    if (small(error) == 0)
      ier3 = 0
    else
      ier3 = 1
      print M1003
    end
#
# TEST OF EZFFTB.
#
    azero = 1.0
    a[0..ND2-1] = 0.0
    b[0..ND2-1] = 0.0
#
    result = DCL::ezfftb(azero, a, b, wezfft)
#
# TEST RESULTS OF EZFFTB
#
    error = 0.0
    for i in 0..N-1
      error = [(1.0 - result[i]).abs, error].max
    end
    if (small(error) == 0)
      ier4 = 0
    else
      ier4 = 1
      print M1004
    end
#
#
    ierror = ier1 + ier2 + ier3 + ier4
    if (ierror == 0)
      print M998
    else
      ierror = 1
      print M999
    end
    ierror
  end
end


include Tfftp

ier = Tfftp::fft_test