code:
! Generated automatically.  DO NOT EDIT!

  integer, parameter :: C_FFTW_R2R_KIND = C_INT32_T

  integer(C_INT), parameter :: FFTW_R2HC = 0
  integer(C_INT), parameter :: FFTW_HC2R = 1
  integer(C_INT), parameter :: FFTW_DHT = 2
  integer(C_INT), parameter :: FFTW_REDFT00 = 3
  integer(C_INT), parameter :: FFTW_REDFT01 = 4
  integer(C_INT), parameter :: FFTW_REDFT10 = 5
  integer(C_INT), parameter :: FFTW_REDFT11 = 6
  integer(C_INT), parameter :: FFTW_RODFT00 = 7
  integer(C_INT), parameter :: FFTW_RODFT01 = 8
  integer(C_INT), parameter :: FFTW_RODFT10 = 9
  integer(C_INT), parameter :: FFTW_RODFT11 = 10
  integer(C_INT), parameter :: FFTW_FORWARD = -1
  integer(C_INT), parameter :: FFTW_BACKWARD = +1
  integer(C_INT), parameter :: FFTW_MEASURE = 0
  integer(C_INT), parameter :: FFTW_DESTROY_INPUT = 1
  integer(C_INT), parameter :: FFTW_UNALIGNED = 2
  integer(C_INT), parameter :: FFTW_CONSERVE_MEMORY = 4
  integer(C_INT), parameter :: FFTW_EXHAUSTIVE = 8
  integer(C_INT), parameter :: FFTW_PRESERVE_INPUT = 16
  integer(C_INT), parameter :: FFTW_PATIENT = 32
  integer(C_INT), parameter :: FFTW_ESTIMATE = 64
  integer(C_INT), parameter :: FFTW_WISDOM_ONLY = 2097152
  integer(C_INT), parameter :: FFTW_ESTIMATE_PATIENT = 128
  integer(C_INT), parameter :: FFTW_BELIEVE_PCOST = 256
  integer(C_INT), parameter :: FFTW_NO_DFT_R2HC = 512
  integer(C_INT), parameter :: FFTW_NO_NONTHREADED = 1024
  integer(C_INT), parameter :: FFTW_NO_BUFFERING = 2048
  integer(C_INT), parameter :: FFTW_NO_INDIRECT_OP = 4096
  integer(C_INT), parameter :: FFTW_ALLOW_LARGE_GENERIC = 8192
  integer(C_INT), parameter :: FFTW_NO_RANK_SPLITS = 16384
  integer(C_INT), parameter :: FFTW_NO_VRANK_SPLITS = 32768
  integer(C_INT), parameter :: FFTW_NO_VRECURSE = 65536
  integer(C_INT), parameter :: FFTW_NO_SIMD = 131072
  integer(C_INT), parameter :: FFTW_NO_SLOW = 262144
  integer(C_INT), parameter :: FFTW_NO_FIXED_RADIX_LARGE_N = 524288
  integer(C_INT), parameter :: FFTW_ALLOW_PRUNING = 1048576

  type, bind(C) :: fftw_iodim
     integer(C_INT) n, is, os
  end type fftw_iodim
  type, bind(C) :: fftw_iodim64
     integer(C_INTPTR_T) n, is, os
  end type fftw_iodim64

  interface
    type(C_PTR) function fftw_plan_dft(rank,n,in,out,sign,flags) bind(C, name='fftw_plan_dft')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_dft

    type(C_PTR) function fftw_plan_dft_1d(n,in,out,sign,flags) bind(C, name='fftw_plan_dft_1d')
      import
      integer(C_INT), value :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_dft_1d

    type(C_PTR) function fftw_plan_dft_2d(n0,n1,in,out,sign,flags) bind(C, name='fftw_plan_dft_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_dft_2d

    type(C_PTR) function fftw_plan_dft_3d(n0,n1,n2,in,out,sign,flags) bind(C, name='fftw_plan_dft_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_dft_3d

    type(C_PTR) function fftw_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags) &
                         bind(C, name='fftw_plan_many_dft')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_many_dft

    type(C_PTR) function fftw_plan_guru_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) &
                         bind(C, name='fftw_plan_guru_dft')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_guru_dft

    type(C_PTR) function fftw_plan_guru_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) &
                         bind(C, name='fftw_plan_guru_split_dft')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: ri
      real(C_DOUBLE), dimension(*), intent(out) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftw_plan_guru_split_dft

    type(C_PTR) function fftw_plan_guru64_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) &
                         bind(C, name='fftw_plan_guru64_dft')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_dft

    type(C_PTR) function fftw_plan_guru64_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) &
                         bind(C, name='fftw_plan_guru64_split_dft')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: ri
      real(C_DOUBLE), dimension(*), intent(out) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_split_dft

    subroutine fftw_execute_dft(p,in,out) bind(C, name='fftw_execute_dft')
      import
      type(C_PTR), value :: p
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftw_execute_dft

    subroutine fftw_execute_split_dft(p,ri,ii,ro,io) bind(C, name='fftw_execute_split_dft')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: ri
      real(C_DOUBLE), dimension(*), intent(inout) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
    end subroutine fftw_execute_split_dft

    type(C_PTR) function fftw_plan_many_dft_r2c(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) &
                         bind(C, name='fftw_plan_many_dft_r2c')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      real(C_DOUBLE), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: flags
    end function fftw_plan_many_dft_r2c

    type(C_PTR) function fftw_plan_dft_r2c(rank,n,in,out,flags) bind(C, name='fftw_plan_dft_r2c')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_r2c

    type(C_PTR) function fftw_plan_dft_r2c_1d(n,in,out,flags) bind(C, name='fftw_plan_dft_r2c_1d')
      import
      integer(C_INT), value :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_r2c_1d

    type(C_PTR) function fftw_plan_dft_r2c_2d(n0,n1,in,out,flags) bind(C, name='fftw_plan_dft_r2c_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_r2c_2d

    type(C_PTR) function fftw_plan_dft_r2c_3d(n0,n1,n2,in,out,flags) bind(C, name='fftw_plan_dft_r2c_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_r2c_3d

    type(C_PTR) function fftw_plan_many_dft_c2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) &
                         bind(C, name='fftw_plan_many_dft_c2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: flags
    end function fftw_plan_many_dft_c2r

    type(C_PTR) function fftw_plan_dft_c2r(rank,n,in,out,flags) bind(C, name='fftw_plan_dft_c2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_c2r

    type(C_PTR) function fftw_plan_dft_c2r_1d(n,in,out,flags) bind(C, name='fftw_plan_dft_c2r_1d')
      import
      integer(C_INT), value :: n
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_c2r_1d

    type(C_PTR) function fftw_plan_dft_c2r_2d(n0,n1,in,out,flags) bind(C, name='fftw_plan_dft_c2r_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_c2r_2d

    type(C_PTR) function fftw_plan_dft_c2r_3d(n0,n1,n2,in,out,flags) bind(C, name='fftw_plan_dft_c2r_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_dft_c2r_3d

    type(C_PTR) function fftw_plan_guru_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftw_plan_guru_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru_dft_r2c

    type(C_PTR) function fftw_plan_guru_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftw_plan_guru_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru_dft_c2r

    type(C_PTR) function fftw_plan_guru_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) &
                         bind(C, name='fftw_plan_guru_split_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftw_plan_guru_split_dft_r2c

    type(C_PTR) function fftw_plan_guru_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) &
                         bind(C, name='fftw_plan_guru_split_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: ri
      real(C_DOUBLE), dimension(*), intent(out) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru_split_dft_c2r

    type(C_PTR) function fftw_plan_guru64_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftw_plan_guru64_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_dft_r2c

    type(C_PTR) function fftw_plan_guru64_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftw_plan_guru64_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_dft_c2r

    type(C_PTR) function fftw_plan_guru64_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) &
                         bind(C, name='fftw_plan_guru64_split_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_split_dft_r2c

    type(C_PTR) function fftw_plan_guru64_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) &
                         bind(C, name='fftw_plan_guru64_split_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: ri
      real(C_DOUBLE), dimension(*), intent(out) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_split_dft_c2r

    subroutine fftw_execute_dft_r2c(p,in,out) bind(C, name='fftw_execute_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftw_execute_dft_r2c

    subroutine fftw_execute_dft_c2r(p,in,out) bind(C, name='fftw_execute_dft_c2r')
      import
      type(C_PTR), value :: p
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
    end subroutine fftw_execute_dft_c2r

    subroutine fftw_execute_split_dft_r2c(p,in,ro,io) bind(C, name='fftw_execute_split_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: ro
      real(C_DOUBLE), dimension(*), intent(out) :: io
    end subroutine fftw_execute_split_dft_r2c

    subroutine fftw_execute_split_dft_c2r(p,ri,ii,out) bind(C, name='fftw_execute_split_dft_c2r')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: ri
      real(C_DOUBLE), dimension(*), intent(inout) :: ii
      real(C_DOUBLE), dimension(*), intent(out) :: out
    end subroutine fftw_execute_split_dft_c2r

    type(C_PTR) function fftw_plan_many_r2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,kind,flags) &
                         bind(C, name='fftw_plan_many_r2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      real(C_DOUBLE), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_plan_many_r2r

    type(C_PTR) function fftw_plan_r2r(rank,n,in,out,kind,flags) bind(C, name='fftw_plan_r2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_plan_r2r

    type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d')
      import
      integer(C_INT), value :: n
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind
      integer(C_INT), value :: flags
    end function fftw_plan_r2r_1d

    type(C_PTR) function fftw_plan_r2r_2d(n0,n1,in,out,kind0,kind1,flags) bind(C, name='fftw_plan_r2r_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_INT), value :: flags
    end function fftw_plan_r2r_2d

    type(C_PTR) function fftw_plan_r2r_3d(n0,n1,n2,in,out,kind0,kind1,kind2,flags) bind(C, name='fftw_plan_r2r_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_FFTW_R2R_KIND), value :: kind2
      integer(C_INT), value :: flags
    end function fftw_plan_r2r_3d

    type(C_PTR) function fftw_plan_guru_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) &
                         bind(C, name='fftw_plan_guru_r2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_plan_guru_r2r

    type(C_PTR) function fftw_plan_guru64_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) &
                         bind(C, name='fftw_plan_guru64_r2r')
      import
      integer(C_INT), value :: rank
      type(fftw_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_DOUBLE), dimension(*), intent(out) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftw_plan_guru64_r2r

    subroutine fftw_execute_r2r(p,in,out) bind(C, name='fftw_execute_r2r')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), dimension(*), intent(inout) :: in
      real(C_DOUBLE), dimension(*), intent(out) :: out
    end subroutine fftw_execute_r2r

    subroutine fftw_destroy_plan(p) bind(C, name='fftw_destroy_plan')
      import
      type(C_PTR), value :: p
    end subroutine fftw_destroy_plan

    subroutine fftw_forget_wisdom() bind(C, name='fftw_forget_wisdom')
      import
    end subroutine fftw_forget_wisdom

    subroutine fftw_cleanup() bind(C, name='fftw_cleanup')
      import
    end subroutine fftw_cleanup

    subroutine fftw_set_timelimit(t) bind(C, name='fftw_set_timelimit')
      import
      real(C_DOUBLE), value :: t
    end subroutine fftw_set_timelimit

    subroutine fftw_plan_with_nthreads(nthreads) bind(C, name='fftw_plan_with_nthreads')
      import
      integer(C_INT), value :: nthreads
    end subroutine fftw_plan_with_nthreads

    integer(C_INT) function fftw_init_threads() bind(C, name='fftw_init_threads')
      import
    end function fftw_init_threads

    subroutine fftw_cleanup_threads() bind(C, name='fftw_cleanup_threads')
      import
    end subroutine fftw_cleanup_threads

    integer(C_INT) function fftw_export_wisdom_to_filename(filename) bind(C, name='fftw_export_wisdom_to_filename')
      import
      character(C_CHAR), dimension(*), intent(in) :: filename
    end function fftw_export_wisdom_to_filename

    subroutine fftw_export_wisdom_to_file(output_file) bind(C, name='fftw_export_wisdom_to_file')
      import
      type(C_PTR), value :: output_file
    end subroutine fftw_export_wisdom_to_file

    type(C_PTR) function fftw_export_wisdom_to_string() bind(C, name='fftw_export_wisdom_to_string')
      import
    end function fftw_export_wisdom_to_string

    subroutine fftw_export_wisdom(write_char,data) bind(C, name='fftw_export_wisdom')
      import
      type(C_FUNPTR), value :: write_char
      type(C_PTR), value :: data
    end subroutine fftw_export_wisdom

    integer(C_INT) function fftw_import_system_wisdom() bind(C, name='fftw_import_system_wisdom')
      import
    end function fftw_import_system_wisdom

    integer(C_INT) function fftw_import_wisdom_from_filename(filename) bind(C, name='fftw_import_wisdom_from_filename')
      import
      character(C_CHAR), dimension(*), intent(in) :: filename
    end function fftw_import_wisdom_from_filename

    integer(C_INT) function fftw_import_wisdom_from_file(input_file) bind(C, name='fftw_import_wisdom_from_file')
      import
      type(C_PTR), value :: input_file
    end function fftw_import_wisdom_from_file

    integer(C_INT) function fftw_import_wisdom_from_string(input_string) bind(C, name='fftw_import_wisdom_from_string')
      import
      character(C_CHAR), dimension(*), intent(in) :: input_string
    end function fftw_import_wisdom_from_string

    integer(C_INT) function fftw_import_wisdom(read_char,data) bind(C, name='fftw_import_wisdom')
      import
      type(C_FUNPTR), value :: read_char
      type(C_PTR), value :: data
    end function fftw_import_wisdom

    subroutine fftw_fprint_plan(p,output_file) bind(C, name='fftw_fprint_plan')
      import
      type(C_PTR), value :: p
      type(C_PTR), value :: output_file
    end subroutine fftw_fprint_plan

    subroutine fftw_print_plan(p) bind(C, name='fftw_print_plan')
      import
      type(C_PTR), value :: p
    end subroutine fftw_print_plan

    type(C_PTR) function fftw_sprint_plan(p) bind(C, name='fftw_sprint_plan')
      import
      type(C_PTR), value :: p
    end function fftw_sprint_plan

    type(C_PTR) function fftw_malloc(n) bind(C, name='fftw_malloc')
      import
      integer(C_SIZE_T), value :: n
    end function fftw_malloc

    type(C_PTR) function fftw_alloc_real(n) bind(C, name='fftw_alloc_real')
      import
      integer(C_SIZE_T), value :: n
    end function fftw_alloc_real

    type(C_PTR) function fftw_alloc_complex(n) bind(C, name='fftw_alloc_complex')
      import
      integer(C_SIZE_T), value :: n
    end function fftw_alloc_complex

    subroutine fftw_free(p) bind(C, name='fftw_free')
      import
      type(C_PTR), value :: p
    end subroutine fftw_free

    subroutine fftw_flops(p,add,mul,fmas) bind(C, name='fftw_flops')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), intent(out) :: add
      real(C_DOUBLE), intent(out) :: mul
      real(C_DOUBLE), intent(out) :: fmas
    end subroutine fftw_flops

    real(C_DOUBLE) function fftw_estimate_cost(p) bind(C, name='fftw_estimate_cost')
      import
      type(C_PTR), value :: p
    end function fftw_estimate_cost

    real(C_DOUBLE) function fftw_cost(p) bind(C, name='fftw_cost')
      import
      type(C_PTR), value :: p
    end function fftw_cost

    integer(C_INT) function fftw_alignment_of(p) bind(C, name='fftw_alignment_of')
      import
      real(C_DOUBLE), dimension(*), intent(out) :: p
    end function fftw_alignment_of

  end interface

  type, bind(C) :: fftwf_iodim
     integer(C_INT) n, is, os
  end type fftwf_iodim
  type, bind(C) :: fftwf_iodim64
     integer(C_INTPTR_T) n, is, os
  end type fftwf_iodim64

  interface
    type(C_PTR) function fftwf_plan_dft(rank,n,in,out,sign,flags) bind(C, name='fftwf_plan_dft')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_dft

    type(C_PTR) function fftwf_plan_dft_1d(n,in,out,sign,flags) bind(C, name='fftwf_plan_dft_1d')
      import
      integer(C_INT), value :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_1d

    type(C_PTR) function fftwf_plan_dft_2d(n0,n1,in,out,sign,flags) bind(C, name='fftwf_plan_dft_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_2d

    type(C_PTR) function fftwf_plan_dft_3d(n0,n1,n2,in,out,sign,flags) bind(C, name='fftwf_plan_dft_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_3d

    type(C_PTR) function fftwf_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags) &
                         bind(C, name='fftwf_plan_many_dft')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_many_dft

    type(C_PTR) function fftwf_plan_guru_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) &
                         bind(C, name='fftwf_plan_guru_dft')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_dft

    type(C_PTR) function fftwf_plan_guru_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) &
                         bind(C, name='fftwf_plan_guru_split_dft')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: ri
      real(C_FLOAT), dimension(*), intent(out) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_split_dft

    type(C_PTR) function fftwf_plan_guru64_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) &
                         bind(C, name='fftwf_plan_guru64_dft')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: sign
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_dft

    type(C_PTR) function fftwf_plan_guru64_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) &
                         bind(C, name='fftwf_plan_guru64_split_dft')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: ri
      real(C_FLOAT), dimension(*), intent(out) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_split_dft

    subroutine fftwf_execute_dft(p,in,out) bind(C, name='fftwf_execute_dft')
      import
      type(C_PTR), value :: p
      complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftwf_execute_dft

    subroutine fftwf_execute_split_dft(p,ri,ii,ro,io) bind(C, name='fftwf_execute_split_dft')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: ri
      real(C_FLOAT), dimension(*), intent(inout) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
    end subroutine fftwf_execute_split_dft

    type(C_PTR) function fftwf_plan_many_dft_r2c(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) &
                         bind(C, name='fftwf_plan_many_dft_r2c')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      real(C_FLOAT), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: flags
    end function fftwf_plan_many_dft_r2c

    type(C_PTR) function fftwf_plan_dft_r2c(rank,n,in,out,flags) bind(C, name='fftwf_plan_dft_r2c')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_r2c

    type(C_PTR) function fftwf_plan_dft_r2c_1d(n,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_1d')
      import
      integer(C_INT), value :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_r2c_1d

    type(C_PTR) function fftwf_plan_dft_r2c_2d(n0,n1,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_r2c_2d

    type(C_PTR) function fftwf_plan_dft_r2c_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_r2c_3d

    type(C_PTR) function fftwf_plan_many_dft_c2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) &
                         bind(C, name='fftwf_plan_many_dft_c2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_INT), value :: flags
    end function fftwf_plan_many_dft_c2r

    type(C_PTR) function fftwf_plan_dft_c2r(rank,n,in,out,flags) bind(C, name='fftwf_plan_dft_c2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_c2r

    type(C_PTR) function fftwf_plan_dft_c2r_1d(n,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_1d')
      import
      integer(C_INT), value :: n
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_c2r_1d

    type(C_PTR) function fftwf_plan_dft_c2r_2d(n0,n1,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_c2r_2d

    type(C_PTR) function fftwf_plan_dft_c2r_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_dft_c2r_3d

    type(C_PTR) function fftwf_plan_guru_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftwf_plan_guru_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_dft_r2c

    type(C_PTR) function fftwf_plan_guru_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftwf_plan_guru_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_dft_c2r

    type(C_PTR) function fftwf_plan_guru_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) &
                         bind(C, name='fftwf_plan_guru_split_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_split_dft_r2c

    type(C_PTR) function fftwf_plan_guru_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) &
                         bind(C, name='fftwf_plan_guru_split_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: ri
      real(C_FLOAT), dimension(*), intent(out) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_split_dft_c2r

    type(C_PTR) function fftwf_plan_guru64_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftwf_plan_guru64_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_dft_r2c

    type(C_PTR) function fftwf_plan_guru64_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) &
                         bind(C, name='fftwf_plan_guru64_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_dft_c2r

    type(C_PTR) function fftwf_plan_guru64_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) &
                         bind(C, name='fftwf_plan_guru64_split_dft_r2c')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_split_dft_r2c

    type(C_PTR) function fftwf_plan_guru64_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) &
                         bind(C, name='fftwf_plan_guru64_split_dft_c2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: ri
      real(C_FLOAT), dimension(*), intent(out) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_split_dft_c2r

    subroutine fftwf_execute_dft_r2c(p,in,out) bind(C, name='fftwf_execute_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: in
      complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out
    end subroutine fftwf_execute_dft_r2c

    subroutine fftwf_execute_dft_c2r(p,in,out) bind(C, name='fftwf_execute_dft_c2r')
      import
      type(C_PTR), value :: p
      complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
    end subroutine fftwf_execute_dft_c2r

    subroutine fftwf_execute_split_dft_r2c(p,in,ro,io) bind(C, name='fftwf_execute_split_dft_r2c')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: in
      real(C_FLOAT), dimension(*), intent(out) :: ro
      real(C_FLOAT), dimension(*), intent(out) :: io
    end subroutine fftwf_execute_split_dft_r2c

    subroutine fftwf_execute_split_dft_c2r(p,ri,ii,out) bind(C, name='fftwf_execute_split_dft_c2r')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: ri
      real(C_FLOAT), dimension(*), intent(inout) :: ii
      real(C_FLOAT), dimension(*), intent(out) :: out
    end subroutine fftwf_execute_split_dft_c2r

    type(C_PTR) function fftwf_plan_many_r2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,kind,flags) &
                         bind(C, name='fftwf_plan_many_r2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      integer(C_INT), value :: howmany
      real(C_FLOAT), dimension(*), intent(out) :: in
      integer(C_INT), dimension(*), intent(in) :: inembed
      integer(C_INT), value :: istride
      integer(C_INT), value :: idist
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_INT), dimension(*), intent(in) :: onembed
      integer(C_INT), value :: ostride
      integer(C_INT), value :: odist
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_plan_many_r2r

    type(C_PTR) function fftwf_plan_r2r(rank,n,in,out,kind,flags) bind(C, name='fftwf_plan_r2r')
      import
      integer(C_INT), value :: rank
      integer(C_INT), dimension(*), intent(in) :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_plan_r2r

    type(C_PTR) function fftwf_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftwf_plan_r2r_1d')
      import
      integer(C_INT), value :: n
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind
      integer(C_INT), value :: flags
    end function fftwf_plan_r2r_1d

    type(C_PTR) function fftwf_plan_r2r_2d(n0,n1,in,out,kind0,kind1,flags) bind(C, name='fftwf_plan_r2r_2d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_INT), value :: flags
    end function fftwf_plan_r2r_2d

    type(C_PTR) function fftwf_plan_r2r_3d(n0,n1,n2,in,out,kind0,kind1,kind2,flags) bind(C, name='fftwf_plan_r2r_3d')
      import
      integer(C_INT), value :: n0
      integer(C_INT), value :: n1
      integer(C_INT), value :: n2
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), value :: kind0
      integer(C_FFTW_R2R_KIND), value :: kind1
      integer(C_FFTW_R2R_KIND), value :: kind2
      integer(C_INT), value :: flags
    end function fftwf_plan_r2r_3d

    type(C_PTR) function fftwf_plan_guru_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) &
                         bind(C, name='fftwf_plan_guru_r2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_plan_guru_r2r

    type(C_PTR) function fftwf_plan_guru64_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) &
                         bind(C, name='fftwf_plan_guru64_r2r')
      import
      integer(C_INT), value :: rank
      type(fftwf_iodim64), dimension(*), intent(in) :: dims
      integer(C_INT), value :: howmany_rank
      type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims
      real(C_FLOAT), dimension(*), intent(out) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
      integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind
      integer(C_INT), value :: flags
    end function fftwf_plan_guru64_r2r

    subroutine fftwf_execute_r2r(p,in,out) bind(C, name='fftwf_execute_r2r')
      import
      type(C_PTR), value :: p
      real(C_FLOAT), dimension(*), intent(inout) :: in
      real(C_FLOAT), dimension(*), intent(out) :: out
    end subroutine fftwf_execute_r2r

    subroutine fftwf_destroy_plan(p) bind(C, name='fftwf_destroy_plan')
      import
      type(C_PTR), value :: p
    end subroutine fftwf_destroy_plan

    subroutine fftwf_forget_wisdom() bind(C, name='fftwf_forget_wisdom')
      import
    end subroutine fftwf_forget_wisdom

    subroutine fftwf_cleanup() bind(C, name='fftwf_cleanup')
      import
    end subroutine fftwf_cleanup

    subroutine fftwf_set_timelimit(t) bind(C, name='fftwf_set_timelimit')
      import
      real(C_DOUBLE), value :: t
    end subroutine fftwf_set_timelimit

    subroutine fftwf_plan_with_nthreads(nthreads) bind(C, name='fftwf_plan_with_nthreads')
      import
      integer(C_INT), value :: nthreads
    end subroutine fftwf_plan_with_nthreads

    integer(C_INT) function fftwf_init_threads() bind(C, name='fftwf_init_threads')
      import
    end function fftwf_init_threads

    subroutine fftwf_cleanup_threads() bind(C, name='fftwf_cleanup_threads')
      import
    end subroutine fftwf_cleanup_threads

    integer(C_INT) function fftwf_export_wisdom_to_filename(filename) bind(C, name='fftwf_export_wisdom_to_filename')
      import
      character(C_CHAR), dimension(*), intent(in) :: filename
    end function fftwf_export_wisdom_to_filename

    subroutine fftwf_export_wisdom_to_file(output_file) bind(C, name='fftwf_export_wisdom_to_file')
      import
      type(C_PTR), value :: output_file
    end subroutine fftwf_export_wisdom_to_file

    type(C_PTR) function fftwf_export_wisdom_to_string() bind(C, name='fftwf_export_wisdom_to_string')
      import
    end function fftwf_export_wisdom_to_string

    subroutine fftwf_export_wisdom(write_char,data) bind(C, name='fftwf_export_wisdom')
      import
      type(C_FUNPTR), value :: write_char
      type(C_PTR), value :: data
    end subroutine fftwf_export_wisdom

    integer(C_INT) function fftwf_import_system_wisdom() bind(C, name='fftwf_import_system_wisdom')
      import
    end function fftwf_import_system_wisdom

    integer(C_INT) function fftwf_import_wisdom_from_filename(filename) bind(C, name='fftwf_import_wisdom_from_filename')
      import
      character(C_CHAR), dimension(*), intent(in) :: filename
    end function fftwf_import_wisdom_from_filename

    integer(C_INT) function fftwf_import_wisdom_from_file(input_file) bind(C, name='fftwf_import_wisdom_from_file')
      import
      type(C_PTR), value :: input_file
    end function fftwf_import_wisdom_from_file

    integer(C_INT) function fftwf_import_wisdom_from_string(input_string) bind(C, name='fftwf_import_wisdom_from_string')
      import
      character(C_CHAR), dimension(*), intent(in) :: input_string
    end function fftwf_import_wisdom_from_string

    integer(C_INT) function fftwf_import_wisdom(read_char,data) bind(C, name='fftwf_import_wisdom')
      import
      type(C_FUNPTR), value :: read_char
      type(C_PTR), value :: data
    end function fftwf_import_wisdom

    subroutine fftwf_fprint_plan(p,output_file) bind(C, name='fftwf_fprint_plan')
      import
      type(C_PTR), value :: p
      type(C_PTR), value :: output_file
    end subroutine fftwf_fprint_plan

    subroutine fftwf_print_plan(p) bind(C, name='fftwf_print_plan')
      import
      type(C_PTR), value :: p
    end subroutine fftwf_print_plan

    type(C_PTR) function fftwf_sprint_plan(p) bind(C, name='fftwf_sprint_plan')
      import
      type(C_PTR), value :: p
    end function fftwf_sprint_plan

    type(C_PTR) function fftwf_malloc(n) bind(C, name='fftwf_malloc')
      import
      integer(C_SIZE_T), value :: n
    end function fftwf_malloc

    type(C_PTR) function fftwf_alloc_real(n) bind(C, name='fftwf_alloc_real')
      import
      integer(C_SIZE_T), value :: n
    end function fftwf_alloc_real

    type(C_PTR) function fftwf_alloc_complex(n) bind(C, name='fftwf_alloc_complex')
      import
      integer(C_SIZE_T), value :: n
    end function fftwf_alloc_complex

    subroutine fftwf_free(p) bind(C, name='fftwf_free')
      import
      type(C_PTR), value :: p
    end subroutine fftwf_free

    subroutine fftwf_flops(p,add,mul,fmas) bind(C, name='fftwf_flops')
      import
      type(C_PTR), value :: p
      real(C_DOUBLE), intent(out) :: add
      real(C_DOUBLE), intent(out) :: mul
      real(C_DOUBLE), intent(out) :: fmas
    end subroutine fftwf_flops

    real(C_DOUBLE) function fftwf_estimate_cost(p) bind(C, name='fftwf_estimate_cost')
      import
      type(C_PTR), value :: p
    end function fftwf_estimate_cost

    real(C_DOUBLE) function fftwf_cost(p) bind(C, name='fftwf_cost')
      import
      type(C_PTR), value :: p
    end function fftwf_cost

    integer(C_INT) function fftwf_alignment_of(p) bind(C, name='fftwf_alignment_of')
      import
      real(C_FLOAT), dimension(*), intent(out) :: p
    end function fftwf_alignment_of

  end interface


code:
  INTEGER FFTW_R2HC
  PARAMETER (FFTW_R2HC=0)
  INTEGER FFTW_HC2R
  PARAMETER (FFTW_HC2R=1)
  INTEGER FFTW_DHT
  PARAMETER (FFTW_DHT=2)
  INTEGER FFTW_REDFT00
  PARAMETER (FFTW_REDFT00=3)
  INTEGER FFTW_REDFT01
  PARAMETER (FFTW_REDFT01=4)
  INTEGER FFTW_REDFT10
  PARAMETER (FFTW_REDFT10=5)
  INTEGER FFTW_REDFT11
  PARAMETER (FFTW_REDFT11=6)
  INTEGER FFTW_RODFT00
  PARAMETER (FFTW_RODFT00=7)
  INTEGER FFTW_RODFT01
  PARAMETER (FFTW_RODFT01=8)
  INTEGER FFTW_RODFT10
  PARAMETER (FFTW_RODFT10=9)
  INTEGER FFTW_RODFT11
  PARAMETER (FFTW_RODFT11=10)
  INTEGER FFTW_FORWARD
  PARAMETER (FFTW_FORWARD=-1)
  INTEGER FFTW_BACKWARD
  PARAMETER (FFTW_BACKWARD=+1)
  INTEGER FFTW_MEASURE
  PARAMETER (FFTW_MEASURE=0)
  INTEGER FFTW_DESTROY_INPUT
  PARAMETER (FFTW_DESTROY_INPUT=1)
  INTEGER FFTW_UNALIGNED
  PARAMETER (FFTW_UNALIGNED=2)
  INTEGER FFTW_CONSERVE_MEMORY
  PARAMETER (FFTW_CONSERVE_MEMORY=4)
  INTEGER FFTW_EXHAUSTIVE
  PARAMETER (FFTW_EXHAUSTIVE=8)
  INTEGER FFTW_PRESERVE_INPUT
  PARAMETER (FFTW_PRESERVE_INPUT=16)
  INTEGER FFTW_PATIENT
  PARAMETER (FFTW_PATIENT=32)
  INTEGER FFTW_ESTIMATE
  PARAMETER (FFTW_ESTIMATE=64)
  INTEGER FFTW_ESTIMATE_PATIENT
  PARAMETER (FFTW_ESTIMATE_PATIENT=128)
  INTEGER FFTW_BELIEVE_PCOST
  PARAMETER (FFTW_BELIEVE_PCOST=256)
  INTEGER FFTW_DFT_R2HC_ICKY
  PARAMETER (FFTW_DFT_R2HC_ICKY=512)
  INTEGER FFTW_NONTHREADED_ICKY
  PARAMETER (FFTW_NONTHREADED_ICKY=1024)
  INTEGER FFTW_NO_BUFFERING
  PARAMETER (FFTW_NO_BUFFERING=2048)
  INTEGER FFTW_NO_INDIRECT_OP
  PARAMETER (FFTW_NO_INDIRECT_OP=4096)
  INTEGER FFTW_ALLOW_LARGE_GENERIC
  PARAMETER (FFTW_ALLOW_LARGE_GENERIC=8192)
  INTEGER FFTW_NO_RANK_SPLITS
  PARAMETER (FFTW_NO_RANK_SPLITS=16384)
  INTEGER FFTW_NO_VRANK_SPLITS
  PARAMETER (FFTW_NO_VRANK_SPLITS=32768)
  INTEGER FFTW_NO_VRECURSE
  PARAMETER (FFTW_NO_VRECURSE=65536)
  INTEGER FFTW_NO_SIMD
  PARAMETER (FFTW_NO_SIMD=131072)

license: mit


code:
module FFTW3
  use, intrinsic :: iso_c_binding
  include 'fftw3.f03'
end module FFTW3


code:
subroutine filbig(dd,npts,f0,newdat,c4a,n4,sq0)

! Filter and downsample the real data in array dd(npts), sampled at 12000 Hz.
! Output is complex, sampled at 1378.125 Hz.

  use, intrinsic :: iso_c_binding
  use FFTW3
  use timer_module, only: timer

  parameter (NSZ=3413)
  parameter (NFFT1=672000,NFFT2=77175)
  parameter (NZ2=1000)
  real*4  dd(npts)                           !Input data
  real*4 rca(NFFT1)
  complex ca(NFFT1/2+1)                      !FFT of input
  complex c4a(NFFT2)                         !Output data
  real*4 s(NZ2)
  real*8 df
  real halfpulse(8)                 !Impulse response of filter (one sided)
  complex cfilt(NFFT2)                       !Filter (complex; imag = 0)
  real rfilt(NFFT2)                          !Filter (real)
  type(C_PTR) :: plan1,plan2,plan3           !Pointers to FFTW plans
  logical first
  equivalence (rfilt,cfilt),(rca,ca)
  data first/.true./
  data halfpulse/114.97547150,36.57879257,-20.93789101,              &
       5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
  common/refspec/dfref,ref(NSZ)
  common/patience/npatience,nthreads
  save first,plan1,plan2,plan3,rfilt,cfilt,df,ca

  if(npts.lt.0) go to 900                    !Clean up at end of program

  if(first) then
     nflags=FFTW_ESTIMATE
     if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
     if(npatience.eq.2) nflags=FFTW_MEASURE
     if(npatience.eq.3) nflags=FFTW_PATIENT
     if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE

! Plan the FFTs just once
     !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
     call fftwf_plan_with_nthreads(nthreads)
     plan1=fftwf_plan_dft_r2c_1d(nfft1,rca,ca,nflags)
     call fftwf_plan_with_nthreads(1)
     plan2=fftwf_plan_dft_1d(nfft2,c4a,c4a,-1,nflags)
     plan3=fftwf_plan_dft_1d(nfft2,cfilt,cfilt,+1,nflags)
     !$omp end critical(fftw)

! Convert impulse response to filter function
     do i=1,nfft2
        cfilt(i)=0.
     enddo
     fac=0.00625/nfft1
     cfilt(1)=fac*halfpulse(1)
     do i=2,8
        cfilt(i)=fac*halfpulse(i)
        cfilt(nfft2+2-i)=fac*halfpulse(i)
     enddo
     call fftwf_execute_dft(plan3,cfilt,cfilt)

     base=real(cfilt(nfft2/2+1))
     do i=1,nfft2
        rfilt(i)=real(cfilt(i))-base
     enddo

     df=12000.d0/nfft1
     first=.false.
  endif

! When new data comes along, we need to compute a new "big FFT"
! If we just have a new f0, continue with the existing data in ca.

  if(newdat.ne.0) then
     call timer('FFTbig  ',0)
     nz=min(npts,nfft1)
     rca(1:nz)=dd(1:nz)
     rca(nz+1:)=0.
     call fftwf_execute_dft_r2c(plan1,rca,ca)
     call timer('FFTbig  ',1)

     ib=0
     do j=1,NSZ
        ia=ib+1
        ib=nint(j*dfref/df)
        fac=sqrt(min(30.0,1.0/ref(j)))
        ca(ia:ib)=fac*conjg(ca(ia:ib))
     enddo
     newdat=0
  endif

! NB: f0 is the frequency at which we want our filter centered.
!     i0 is the bin number in ca closest to f0.
  call timer('loops   ',0)
  i0=nint(f0/df) + 1
  nh=nfft2/2
  do i=1,nh                                !Copy data into c4a and apply
     j=i0+i-1                              !the filter function
     if(j.ge.1 .and. j.le.nfft1/2+1) then
        c4a(i)=rfilt(i)*ca(j)
     else
        c4a(i)=0.
     endif
  enddo
  do i=nh+1,nfft2
     j=i0+i-1-nfft2
!     if(j.lt.1) j=j+nfft1                  !nfft1 was nfft2
     if(j.ge.1) then
        c4a(i)=rfilt(i)*ca(j)
     else
        c4a(i)=rfilt(i)*conjg(ca(2-j))
     endif
  enddo

  nadd=nfft2/NZ2
  i=0
  do j=1,NZ2
     s(j)=0.
     do n=1,nadd
        i=i+1
        s(j)=s(j) + real(c4a(i))**2 + aimag(c4a(i))**2
     enddo
  enddo
  call pctile(s,NZ2,30,sq0)
  call timer('loops   ',1)

! Do the short reverse transform, to go back to time domain.
  call timer('FFTsmall',0)
  call fftwf_execute_dft(plan2,c4a,c4a)
  call timer('FFTsmall',1)
  n4=min(npts/8,nfft2)
  return

900 continue

  !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
  call fftwf_destroy_plan(plan1)
  call fftwf_destroy_plan(plan2)
  call fftwf_destroy_plan(plan3)
  !$omp end critical(fftw)

  return
end subroutine filbig


code:
subroutine fitcal(x,y,r,iz,a,b,sigmaa,sigmab,rms)
  implicit real*8 (a-h,o-z)
  real*8 x(iz),y(iz),r(iz)

  sx=0.d0
  sy=0.d0
  sxy=0.d0
  sx2=0.d0
  do i=1,iz
     sx=sx + x(i)
     sy=sy + y(i)
     sxy=sxy + x(i)*y(i)
     sx2=sx2 + x(i)*x(i)
  enddo
  delta=iz*sx2 - sx*sx
  a=(sx2*sy - sx*sxy)/delta
  b=(iz*sxy - sx*sy)/delta

  sq=0.d0
  do i=1,iz
     r(i)=y(i) - (a + b*x(i))
     sq=sq + r(i)**2
  enddo
  rms=0.
  sigmaa=0.
  sigmab=0.
  if(iz.ge.3) then
     rms=sqrt(sq/(iz-2))
     sigmaa=sqrt(rms*rms*sx2/delta)
     sigmab=sqrt(iz*rms*rms/delta)
  endif

  return
end subroutine fitcal


code:
subroutine fix_contest_msg(mygrid,msg)

! If distance from mygrid to grid1 is more thsn 10000 km, change "grid1"
! to "R grid2" where grid2 is the antipodes of grid1.

  character*6 mygrid
  character*22 msg
  character*6 g1,g2
  logical isgrid

  isgrid(g1)=g1(1:1).ge.'A' .and. g1(1:1).le.'R' .and. g1(2:2).ge.'A' .and. &
       g1(2:2).le.'R' .and. g1(3:3).ge.'0' .and. g1(3:3).le.'9' .and.       &
       g1(4:4).ge.'0' .and. g1(4:4).le.'9' .and. g1(1:4).ne.'RR73'

  n=len(trim(msg))
  if(n.lt.4) return

  g1=msg(n-3:n)//'  '
  if(isgrid(g1)) then
     call azdist(mygrid,g1,0.d0,nAz,nEl,nDmiles,nDkm,nHotAz,nHotABetter)
     if(ndkm.gt.10000) then
        call grid2deg(g1,dlong,dlat)
        dlong=dlong+180.0
        if(dlong.gt.180.0) dlong=dlong-360.0
        dlat=-dlat
        call deg2grid(dlong,dlat,g2)
        msg=msg(1:n-4)//'R '//g2(1:4)
     endif
  endif

  return
end subroutine fix_contest_msg


code:
program fixwav

  use wavhdr
  parameter (NBANDS=23,NMODES=11)
  parameter(NMAX=120*12000)
  type(hdr) h
  integer*2 id2(NMAX),id(4)
  character*8 mode,mode0
  character*6 band
  character*12 arg
  character*80 infile
  character*1 c
  logical ok

  nargs=iargc()
  if(nargs.ne.1 .and. nargs.ne.5) then
     print*,'Usage: fixwav [fMHz mode submode TRperiod] infile'
     go to 999
  endif
  call getarg(1,infile)
  if(nargs.eq.1) go to 10
  read(infile,*) fMHz
  call getarg(2,mode0)
  call getarg(3,arg)
  nsubmode0=-1
  do i=1,8
     if(arg(1:1).eq.char(ichar('A')-1+i)) nsubmode0=i
  enddo
  if(nsubmode0.lt.0) read(arg,*) nsubmode0
  call getarg(4,arg)
  read(arg,*) ntrperiod0
  call getarg(5,infile)

10 open(10,file=infile,status='old',access='stream')
  read(10) h
  npts=h%ndata/2
  nfsample=h%nsamrate
  read(10) id2(1:npts)
  write(*,1002) h%nchan2,h%nbitsam2,h%nsamrate,npts
1002 format('Channels:',i2,'   Bits/sample:',i3,'   Sample rate:',i6,   &
          '   Npts:',i8)

  call get_wsjtx_wav_params(id2,band,mode,nsubmode,ntrperiod,ok)
  if(nfsample.ne.11025 .and. nfsample.ne.12000) ok=.false.
  if(ok) write(*,1010) band,ntrperiod,mode,char(ichar('A')-1+id2(3))
1010 format('Band: ',a6,'  T/R period:',i4,'   Mode: ',a8,1x,a1)

  if(.not.ok) write(*,'(a)') 'File has no valid WSJT-X params.'
  if(ok .and. nargs.eq.1) go to 999

  write(*,'(a)',advance='no') 'Do you want to add or change them? (Y/N): '
  read*,c
  if(c.ne.'y' .and. c.ne.'Y') go to 999

  print*,fMHz,mode0,nsubmode0,ntrperiod0
  call set_wsjtx_wav_params(fMHz,mode0,nsubmode0,ntrperiod0,id2)
  band=""
  mode=""
  nsubmode=0
  ntrperiod=0
  call get_wsjtx_wav_params(id2,band,mode,nsubmode,ntrperiod,ok)
  write(*,1010) band,ntrperiod,mode,char(ichar('A')-1+id2(3))
  rewind 10
  write(10) h,id2(1:npts)

999 end program fixwav


code:
subroutine flat1(savg,iz,nsmo,syellow)

  real savg(iz)
  real syellow(iz)
  real x(8192)

  ia=nsmo/2 + 1
  ib=iz - nsmo/2 - 1
  nstep=20
  nh=nstep/2
  do i=ia,ib,nstep
     call pctile(savg(i-nsmo/2),nsmo,50,x(i))
     x(i-nh:i+nh-1)=x(i)
  enddo
  x(1:ia-1)=x(ia)
  x(ib+1:iz)=x(ib)

  x0=0.001*maxval(x(iz/10:(9*iz)/10))
  syellow(1:iz)=savg(1:iz)/(x(1:iz)+x0)

  return
end subroutine flat1


code:
subroutine flat1a(psavg,nsmo,s2,nh,nsteps,nhmax,nsmax)

  real psavg(nh)
  real s2(nhmax,nsmax)
  real x(8192)

  ia=nsmo/2 + 1
  ib=nh - nsmo/2 - 1
  do i=ia,ib
     call pctile(psavg(i-nsmo/2),nsmo,50,x(i))
  enddo
  do i=1,ia-1
     x(i)=x(ia)
  enddo
  do i=ib+1,nh
     x(i)=x(ib)
  enddo

  do i=1,nh
     psavg(i)=psavg(i)/x(i)
     do j=1,nsteps
        s2(i,j)=s2(i,j)/x(i)
     enddo
  enddo

  return
end subroutine flat1a


code:
subroutine flat1b(psavg,nsmo,s2,nh,nsteps,nhmax,nsmax)

  real psavg(nh)
  real s2(nhmax,nsmax)
  real x(8192)

  ia=nsmo/2 + 1
  ib=nh - nsmo/2 - 1
  do i=ia,ib
     call pctile(psavg(i-nsmo/2),nsmo,50,x(i))
  enddo
  do i=1,ia-1
     x(i)=x(ia)
  enddo
  do i=ib+1,nh
     x(i)=x(ib)
  enddo

  do i=1,nh
     psavg(i)=psavg(i)/x(i)
     do j=1,nsteps
        s2(i,j)=s2(i,j)/x(i)
     enddo
  enddo

  return
end subroutine flat1b


code:
subroutine flat2(s,nz,ref)

  parameter (NSMAX=6827)
  real s(NSMAX)
  real ref(NSMAX)

  nsmo=10
  ia=nsmo+1
  ib=nz-nsmo-1
  do i=ia,ib
     call pctile(s(i-nsmo),2*nsmo+1,5,ref(i))
  enddo

  ref(:ia-1)=ref(ia)
  ref(ib+1:)=ref(ib)

  return
end subroutine flat2


code:
subroutine flat4(s,npts0,nflatten)

! Flatten a spectrum for optimum display
! Input:  s(npts)    Linear scale in power
!         nflatten   If nflatten=0, convert to dB but do not flatten
! Output: s(npts)    Flattened, with dB scale

  implicit real*8 (a-h,o-z)
  real*4 s(6827)
  real*4 base
  real*8 x(1000),y(1000),a(5)
  data nseg/10/,npct/10/

  npts=min(6827,npts0)
  if(s(1).gt.1.e29) go to 900         !Boundary between Rx intervals: do nothing
  do i=1,npts
     s(i)=10.0*log10(s(i))            !Convert to dB scale
  enddo

  if(nflatten.gt.0) then
     nterms=5
     if(nflatten.eq.2) nterms=1
     nlen=npts/nseg                   !Length of test segment
     i0=npts/2                        !Midpoint
     k=0
     do n=1,nseg                      !Skip first segment, likely rolloff here
        ib=n*nlen
        ia=ib-nlen+1
        if(n.eq.nseg) ib=npts
        call pctile(s(ia),ib-ia+1,npct,base) !Find lowest npct of points
        do i=ia,ib
           if(s(i).le.base) then
              if (k.lt.1000) k=k+1    !Save these "lower envelope" points
              x(k)=i-i0
              y(k)=s(i)
           endif
        enddo
     enddo
     kz=k
     a=0.

     call polyfit(x,y,y,kz,nterms,0,a,chisqr)  !Fit a low-order polynomial

     do i=1,npts
        t=i-i0
        yfit=a(1)+t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5)))))
        s(i)=s(i)-yfit                !Subtract the fitted baseline
     enddo
  endif

900 return
end subroutine flat4


code:
subroutine flat65(ss,nhsym,maxhsym,nsz,ref)

  real stmp(nsz)
  real ss(maxhsym,nsz)
  real ref(nsz)

  npct=28                                       !Somewhat arbitrary
  do i=1,nsz
     call pctile(ss(1,i),nhsym,npct,stmp(i))
  enddo

  nsmo=33
  ia=nsmo/2 + 1
  ib=nsz - nsmo/2 - 1
  do i=ia,ib
     call pctile(stmp(i-nsmo/2),nsmo,npct,ref(i))
  enddo
  ref(:ia-1)=ref(ia)
  ref(ib+1:)=ref(ib)
  ref=4.0*ref

  return
end subroutine flat65


code:
!-------------------------------------------------------------------------------
!
! This file is part of the WSPR application, Weak Signal Propagation Reporter
!
! File Name:    fmeasure.f90
! Description:
!
! Copyright (C) 2001-2014 Joseph Taylor, K1JT
! License: GPL-3
!
! This program is free software; you can redistribute it and/or modify it under
! the terms of the GNU General Public License as published by the Free Software
! Foundation; either version 3 of the License, or (at your option) any later
! version.
!
! This program is distributed in the hope that it will be useful, but WITHOUT
! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
! details.
!
! You should have received a copy of the GNU General Public License along with
! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
! Street, Fifth Floor, Boston, MA 02110-1301, USA.
!
!-------------------------------------------------------------------------------
program fmeasure

  parameter(NZ=1000)
  implicit real*8 (a-h,o-z)
  character infile*50
  character line*80

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage:   fmeasure <infile>'
     print*,'Example: fmeasure fmtave.out'
     go to 999
  endif
  call getarg(1,infile)

  open(10,file=infile,status='old',err=997)
  open(11,file='fcal.out',status='old',err=998)
  open(12,file='fmeasure.out',status='unknown')

  read(11,*) a,b

  write(*,1000)
  write(12,1000)
1000 format('    Freq     DF     A+B*f     Corrected  Offset'/        &
            '   (MHz)    (Hz)    (Hz)        (MHz)      (Hz)'/        &
            '-----------------------------------------------')
  i=0
  do j=1,9999
     read(10,1010,end=999) line
1010 format(a80)
     i0=index(line,' 0 ')
     if(i0.gt.0) then
        read(line,*,err=5) f,df
        dial_error=a + b*f
        fcor=f + 1.d-6*df - 1.d-6*dial_error
        offset_hz=1.d6*(fcor-f)
        write(*,1020)  f,df,dial_error,fcor,offset_hz
        write(12,1020) f,df,dial_error,fcor,offset_hz
1020    format(3f8.3,f15.9,f8.2)
     endif
5    continue
  enddo

  go to 999

997 print*,'Cannot open input file: ',infile
  go to 999
998 print*,'Cannot open fcal.out'

999 end program fmeasure


code:
program fmtave

! Average groups of frequency-calibration measurements.

  implicit real*8 (a-h,o-z)
  character infile*80
  character*8 cutc,cutc1

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage:   fmtave <infile>'
     print*,'Example: fmtave fmt.all'
     go to 999
  endif
  call getarg(1,infile)

  open(10,file=infile,status='old')
  open(12,file='fmtave.out',status='unknown')

  write(*,1000)
1000 format('    Freq     DF    CAL   N     rms    UTC     Call'/       &
            '   (kHz)    (Hz)    ?         (Hz)'/                       &
            '----------------------------------------------------')
  nkhz0=0
  sum=0.d0
  sumsq=0.d0
  n=0
  do i=1,99999
     read(10,*,end=10) cutc,nkHz,ncal,noffset,faudio,df,dblevel,snr
     if((nkHz.ne.nkHz0) .and. i.ne.1) then
        ave=sum/n
        rms=0.d0
        if(n.gt.1) then
           rms=sqrt(abs(sumsq - sum*sum/n)/(n-1.d0))
        endif
        fMHz=0.001d0*nkHz0
        write(*,1010)  fMHz,ave,ncal0,n,rms,cutc1
        write(12,1010) fMHz,ave,ncal0,n,rms,cutc1
1010    format(f8.3,f9.3,i4,i5,f8.2,2x,a8,2x,a6)
        sum=0.d0
        sumsq=0.d0
        n=0
     endif
     dial_error=faudio-noffset
     sum=sum + dial_error
     sumsq=sumsq + dial_error**2
     n=n+1
     if(n.eq.1) then
        cutc1=cutc
        ncal0=ncal
     endif
     nkHz0=nkHz
  enddo

10 ave=sum/n
  rms=0.d0
  if(n.gt.0) then
     rms=sqrt((sumsq - sum*sum/n)/(n-1.d0))
  endif
  fMHz=0.001d0*nkHz
  write(*,1010)  fMHz,ave,ncal,n,rms,cutc1
  write(12,1010) fMHz,ave,ncal,n,rms,cutc1

999 end program fmtave


code:
subroutine fmtmsg(msg,iz)

  character*22 msg

! Convert all letters to upper case
  iz=22
  do i=1,22
     if(msg(i:i).ge.'a' .and. msg(i:i).le.'z')                       &
          msg(i:i)= char(ichar(msg(i:i))+ichar('A')-ichar('a'))
     if(msg(i:i).ne.' ') iz=i
  enddo

  do iter=1,5                           !Collapse multiple blanks into one
     ib2=index(msg(1:iz),'  ')
     if(ib2.lt.1) go to 100
     msg=msg(1:ib2)//msg(ib2+2:)
     iz=iz-1
  enddo

100 return
end subroutine fmtmsg


code:
subroutine foldspec9f(s1,nq,jz,ja,jb,s2)

! Fold symbol spectra (quarter-symbol steps) from s1 into s2

  real s1(nq,jz)
  real s2(240,340)                       !340 = 4*85
  integer nsum(340)

  s2=0.
  nsum=0

  do j=ja,jb
     k=mod(j-1,340)+1
     nsum(k)=nsum(k)+1
     do i=1,NQ
        s2(i,k)=s2(i,k) + s1(i,j)
     enddo
  enddo

  do k=1,340
     fac=1.0
     if(nsum(k).gt.0) fac=1.0/nsum(k)
     s2(1:nq,k)=fac*s2(1:nq,k)
  enddo

  ave=sum(s2)/(340.0*nq)
  if(ave.gt.0.0) s2=s2/ave

  return
end subroutine foldspec9f


code:
subroutine four2a(a,nfft,ndim,isign,iform)

! IFORM = 1, 0 or -1, as data is
! complex, real, or the first half of a complex array.  Transform
! values are returned in array DATA.  They are complex, real, or
! the first half of a complex array, as IFORM = 1, -1 or 0.

! The transform of a real array (IFORM = 0) dimensioned N(1) by N(2)
! by ... will be returned in the same array, now considered to
! be complex of dimensions N(1)/2+1 by N(2) by ....  Note that if
! IFORM = 0 or -1, N(1) must be even, and enough room must be
! reserved.  The missing values may be obtained by complex conjugation.

! The reverse transformation of a half complex array dimensioned
! N(1)/2+1 by N(2) by ..., is accomplished by setting IFORM
! to -1.  In the N array, N(1) must be the true N(1), not N(1)/2+1.
! The transform will be real and returned to the input array.

! This version of four2a makes calls to the FFTW library to do the
! actual computations.

  parameter (NPMAX=2100)                 !Max numberf of stored plans
  parameter (NSMALL=16384)               !Max size of "small" FFTs
  complex a(nfft)                        !Array to be transformed
  complex aa(NSMALL)                     !Local copy of "small" a()
  integer nn(NPMAX),ns(NPMAX),nf(NPMAX)  !Params of stored plans
  integer*8 nl(NPMAX),nloc               !More params of plans
  integer*8 plan(NPMAX)                  !Pointers to stored plans
  logical found_plan
  data nplan/0/                          !Number of stored plans
  common/patience/npatience,nthreads     !Patience and threads for FFTW plans
  include 'fftw3.f90'                    !FFTW definitions
  save plan,nplan,nn,ns,nf,nl

  if(nfft.lt.0) go to 999

  nloc=loc(a)

  found_plan = .false.
  !$omp critical(four2a_setup)
  do i=1,nplan
     if(nfft.eq.nn(i) .and. isign.eq.ns(i) .and.                     &
          iform.eq.nf(i) .and. nloc.eq.nl(i)) then
        found_plan = .true.
        exit
     end if
  enddo

  if(i.ge.NPMAX) stop 'Too many FFTW plans requested.'

  if (.not. found_plan) then
     nplan=nplan+1
     i=nplan

     nn(i)=nfft
     ns(i)=isign
     nf(i)=iform
     nl(i)=nloc

! Planning: FFTW_ESTIMATE, FFTW_ESTIMATE_PATIENT, FFTW_MEASURE,
!            FFTW_PATIENT,  FFTW_EXHAUSTIVE
     nflags=FFTW_ESTIMATE
     if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
     if(npatience.eq.2) nflags=FFTW_MEASURE
     if(npatience.eq.3) nflags=FFTW_PATIENT
     if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE

     if(nfft.le.NSMALL) then
        jz=nfft
        if(iform.eq.0) jz=nfft/2
        aa(1:jz)=a(1:jz)
     endif

     !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
     if(isign.eq.-1 .and. iform.eq.1) then
        call sfftw_plan_dft_1d(plan(i),nfft,a,a,FFTW_FORWARD,nflags)
     else if(isign.eq.1 .and. iform.eq.1) then
        call sfftw_plan_dft_1d(plan(i),nfft,a,a,FFTW_BACKWARD,nflags)
     else if(isign.eq.-1 .and. iform.eq.0) then
        call sfftw_plan_dft_r2c_1d(plan(i),nfft,a,a,nflags)
     else if(isign.eq.1 .and. iform.eq.-1) then
        call sfftw_plan_dft_c2r_1d(plan(i),nfft,a,a,nflags)
     else
        stop 'Unsupported request in four2a'
     endif
     !$omp end critical(fftw)

     if(nfft.le.NSMALL) then
        jz=nfft
        if(iform.eq.0) jz=nfft/2
        a(1:jz)=aa(1:jz)
     endif
  end if
  !$omp end critical(four2a_setup)

  call sfftw_execute(plan(i))
  return

999 continue

  !$omp critical(four2a)
  do i=1,nplan
! The test is only to silence a compiler warning:
     if(ndim.ne.-999) then
        !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
        call sfftw_destroy_plan(plan(i))
        !$omp end critical(fftw)
     end if
  enddo

  nplan=0
  !$omp end critical(four2a)

  return
end subroutine four2a


code:
subroutine fox_rx(fail,called,fm,hm)

! Given fm, recently transmitted by Fox, determine hm -- the next
! message for Hound to transmit

  parameter (MAXSIG=5,NCALLS=268)
  character*6 xcall(NCALLS)
!  character*8 mycall_plus
  character*4 xgrid(NCALLS)
  integer isnr(NCALLS)
  character*32 fm
  character*22 hm
  character*6 cx,called,MyCall
  character*4 gx
  common/dxpfifo/nc,isnr,xcall,xgrid
  data MyCall/'KH1DX'/
  save

  call random_number(r)
  if(r.lt.fail) fm=''                        !Hound fails to copy
  i1=index(fm,MyCall)
  if(fm(1:3).eq.'CQ ' .and. i1.ge.4) then
     call dxped_fifo(cx,gx,isnrx)
     ntimes=1
     write(hm,1000) MyCall,cx,gx
1000 format(a6,1x,a6,1x,a4)
  endif

! Check for a "RR73" message
  ia=index(fm,trim(cx))
  ib=index(fm,';')
  ic=index(fm,trim(called))
  id=index(fm,'RR73;')
  if((ia.eq.1 .or. ic.eq.ib+2) .and. id.ge.4) then
     i1=index(fm,';')+2
     i2=index(fm,'<')-2
     cx=fm(i1:i2)                      !Callsign for next QSO
     call random_number(r)
     ireport=nint(-20+40*r)
! Send report to next caller
     write(hm,1004) MyCall,cx,ireport
1004 format(a6,1x,a6,' R',i3.2)
     if(hm(16:16).eq.' ') hm(16:16)='+'
  endif

! Check for a message with a report to Hound
  i1=index(fm,trim(called))
  i2=index(fm,MyCall)
  if(i1.eq.1 .and. i2.ge.5 .and.   &
       (index(fm,'+').ge.8 .or. index(fm,'-').ge.8)) then
! Send "R+rpt" to Fox
     write(hm,1004) MyCall,called,isnrx
     if(hm(16:16).eq.' ') hm(16:16)='+'
  endif

! Collapse multiple blanks in message
  iz=len(trim(hm))
  do iter=1,5
     ib2=index(hm(1:iz),'  ')
     if(ib2.lt.1) exit
     hm=hm(1:ib2)//hm(ib2+2:)
     iz=iz-1
  enddo

  return
end subroutine fox_rx


code:
program fox_sim

! Simulates QSO exchanges using the proposed FT8 "DXpedition" mode.
  parameter (MAXSIG=5,NCALLS=268)
  character*6 xcall(NCALLS)
  character*4 xgrid(NCALLS)
  integer isnr(NCALLS)

  character*32 fmsg(MAXSIG),fm
  character*22 hmsg(MAXSIG),hm
  character*16 log
  character*6 called(MAXSIG)
  character*4 gcalled(MAXSIG)
  character*6 MyCall
  character*4 MyGrid
  character*8 arg
  character*1 c1,c2,c3,c4
  integer ntot(MAXSIG),irate(MAXSIG),ntimes(MAXSIG)
  logical logit
  common/dxpfifo/nc,isnr,xcall,xgrid

  nargs=iargc()
  if(nargs.ne.2 .and. nargs.ne.4) then
     print*,'Usage: fox_sim  nseq maxtimes'
     print*,'       fox_sim  nseq maxtimes nsig fail'
     print*,' '
     print*,' nseq:     number of T/R sequences to execute'
     print*,' maxtimes: number of repeats of same Tx message'
     print*,' nsig:     number of simultaneous Tx sigals'
     print*,' fail:     receiving error rate'
     go to 999
  endif
  ii1=1
  ii2=5
  jj1=0
  jj2=5
  nseq=80
  if(nargs.ge.2) then
     call getarg(1,arg)
     read(arg,*) nseq
     call getarg(2,arg)
     read(arg,*) maxtimes
  endif
  if(nargs.eq.4) then
     call getarg(3,arg)
     read(arg,*) nsig
     call getarg(4,arg)
     read(arg,*) fail
     ii1=nsig
     ii2=nsig
     jj1=nint(10*fail)
     jj2=nint(10*fail)
  endif

! Read a file with calls and grids; insert random S/N values.
! This is used in place of an operator-selected FIFO
  open(10,file='xcall.txt',status='old')
  do i=1,NCALLS
     read(10,1000) xcall(i),xgrid(i)
1000 format(a6,7x,a4)
     if(i.ne.-99) cycle
     j=mod(i-1,26)
     c1=char(ichar('A')+j)
     k=mod((i-1)/26,26)
     c2=char(ichar('A')+k)
     n=mod((i-1)/260,10)
     c3=char(ichar('0')+n)
     xcall(i)='K'//c2//c3//c1//c1//c1

     j=mod(i-1,18)
     c1=char(ichar('A')+j)
     k=mod((i-1)/18,18)
     c2=char(ichar('A')+k)
     n=mod((i-1)/10,10)
     c4=char(ichar('0')+n)
     n=mod((i-1)/100,10)
     c3=char(ichar('0')+n)
     xgrid(i)=c1//c2//c3//c4

     call random_number(x)
     isnr(i)=-20+int(40*x)
  enddo
!  close(10)

! Write headings for the summary file
  minutes=nseq/4
  write(13,1002) nseq,minutes,maxtimes
1002 format(/'Nseq:',i4,'   Minutes:',i3,'   Maxtimes:',i2//                   &
    18x,'Logged QSOs',22x,'Rate (QSOs/hour)'/                                  &
    'fail Nsig: 1     2     3     4     5          1     2     3     4     5'/ &
    71('-'))

  write(*,1003)
1003 format('Seq  s n Fox messages                      Hound messages    Logged info        i Rate'/87('-'))

  ntot=0
  irate=0
  MyCall='KH1DX'
  MyGrid='AJ10'

  do jj=jj1,jj2                        !Loop over Rx failure rates
     fail=0.1*jj
     do ii=ii1,ii2                     !Loop over range of nsig
        nc=0                           !Set FIFO pointer to top
        ntimes=1
        nsig=ii
        nlogged=0
        fmsg="CQ KH1DX AJ10"
        hmsg=""
        called="      "
        do iseq=0,nseq                 !Loop over specified number of sequences
           if(iand(iseq,1).eq.0) then
              do j=1,nsig              !Loop over Fox's Tx slots
                 fm=fmsg(j)
                 hm=hmsg(j)

! Call fox_tx to determine the next Tx message for this slot
                 call fox_tx(maxtimes,fail,called(j),gcalled(j),hm,fm,    &
                      ntimes(j),log,logit)

                 fmsg(j)=fm
                 if(logit) then
! Log this QSO
                    nlogged=nlogged+1
                    nrate=0
                    if(iseq.gt.0) nrate=nint(nlogged*240.0/iseq)
                    write(*,1010) iseq,j,ntimes(j),fmsg(j),log,nlogged,nrate
1010                format(i4.4,2i2,1x,a32,20x,a16,2i4)
                    ! call log_routine()
                 else
                    write(*,1010) iseq,j,ntimes(j),fmsg(j)
                 endif
              enddo
              ! call transmit()
           endif

           if(iand(iseq,1).eq.1) then
              do j=1,nsig              !Listen for expected responses
                 fm=fmsg(j)
                 call fox_rx(fail,called(j),fm,hm)
                 if(j.ge.2) then
                    if(hm.eq.hmsg(j-1)) hm=""
                 endif
                 hmsg(j)=hm
                 write(*,1020) iseq,j,hmsg(j)
1020             format(i4.4,i2,37x,a22)
              enddo
           endif
           write(*,1021)
1021       format(87('-'))
        enddo
        ntot(ii)=nlogged
        irate(ii)=0
        if(iseq.gt.0) irate(ii)=nint(nlogged*3600.0/(15*iseq))
        write(*,1030) nsig,fail,nlogged,nc
1030    format(/'Nsig:',i3,'   Fail:',f4.1,'   Logged QSOs:',i4,      &
             '   Final nc:',i4)
     enddo

! Write the summary file
     write(13,1100) fail,ntot,irate
1100 format(f4.1,2x,5i6,5x,5i6)
  enddo

999 end program fox_sim


code:
subroutine fox_tx(maxtimes,fail,called,gcalled,hm,fm,ntimes,log,logit)

! Determine fm, the next message for Fox to transmit in this slot

  character*32 fm
  character*22 hm
  character*4 g4,MyGrid,gcalled,gx,gy
  character*6 MyCall,called,cx,cy
  character*16 log
  logical isgrid,logit
  data MyCall/'KH1DX'/,MyGrid/'AJ10'/
  save

  isgrid(g4)=g4(1:1).ge.'A' .and. g4(1:1).le.'R' .and. g4(2:2).ge.'A' .and. &
       g4(2:2).le.'R' .and. g4(3:3).ge.'0' .and. g4(3:3).le.'9' .and.       &
       g4(4:4).ge.'0' .and. g4(4:4).le.'9' .and. g4(1:4).ne.'RR73'

  logit=.false.
  n=len(trim(hm))
  g4=""
  if(n.gt.8) g4=hm(n-3:n)
  call random_number(r)
  if(r.lt.fail .and. .not.isgrid(g4)) hm=""        !Fox failed to copy

  i2=len(trim(hm))
  if(i2.gt.10) then
     i1=index(hm,' ')
     i3=index(hm(i1+1:),' ') + i1
     cx=hm(i1+1:i3)
     gx=hm(i2-3:i2)
     i4=index(hm,MyCall)

! Check for a new caller
     if(i4.eq.1 .and. isgrid(gx)) then
        call random_number(r)
        isent=nint(-20+40*r)
        write(fm,1002) cx,MyCall,isent
1002    format(a6,1x,a6,i4.2)
        if(fm(15:15).eq.' ') fm(15:15)='+'
        called=cx
        gcalled=gx
     endif
     log=''

! Check for message with R+rpt
     if(i4.eq.1 .and. cx.eq.called .and.                         &
          (index(hm,'R+').ge.8 .or. index(hm,'R-').ge.8)) then
        write(log,1006) called,gcalled,isent        !Format a log entry
1006    format(a6,2x,a4,i4.2)
        if(log(14:14).eq.' ') log(14:14)='+'
        logit=.true.
        call dxped_fifo(cy,gy,isnry)
! If FIFO is empty we should call CQ in this slot
        ntimes=1
        write(fm,1008) cx,cy,isnry
1008    format(a6,' RR73; ',a6,1x,'<KH1DX>',i4.2)
        if(fm(29:29).eq.' ') fm(29:29)='+'
        called=cy
        gcalled=gy
     endif
  endif

  if(hm.eq.'') then
     if(fm(1:3).ne.'CQ ') then
!        if(ntimes.lt.maxtimes) then
           ntimes=ntimes+1
!        else
!           ntimes=1
! If FIFO is empty we should call CQ in this slot
!           call dxped_fifo(cy,gy,isnry)
!           call random_number(r)
!           isnr=nint(-20+40*r)
!           write(fm,1010) cy,gy,isnr
           write(fm,1010) called,MyCall,isent
1010       format(a6,1x,a6,i4.2)
           if(fm(15:15).eq.' ') fm(15:15)='+'
!        endif
     endif
  endif

! Collapse multiple blanks in message
  iz=len(trim(fm))
  do iter=1,5
     ib2=index(fm(1:iz),'  ')
     if(ib2.lt.1) exit
     fm=fm(1:ib2)//fm(ib2+2:)
     iz=iz-1
  enddo

! Generate waveform for fm
  return
end subroutine fox_tx


code:
subroutine fqso_first(nfqso,ntol,ca,ncand)

! If a candidate was found within +/- ntol of nfqso, move it into ca(1).

  type candidate
     real freq
     real dt
     real sync
     real flip
  end type candidate
  type(candidate) ca(300),cb

  dmin=1.e30
  i0=0
  do i=1,ncand
     d=abs(ca(i)%freq-nfqso)
     if(d.lt.dmin) then
        i0=i
        dmin=d
     endif
  enddo

  if(dmin.lt.float(ntol)) then
     cb=ca(i0)
     do i=i0,2,-1
        ca(i)=ca(i-1)
     enddo
     ca(1)=cb
  endif

  return
end subroutine fqso_first


code:
subroutine freqcal(id2,k,nkhz,noffset,ntol,line)

  parameter (NZ=30*12000,NFFT=55296,NH=NFFT/2)
  integer*2 id2(0:NZ-1)
  complex sp,sn
  real x(0:NFFT-1)
  real xi(0:NFFT-1)
  real w(0:NFFT-1)                      !Window function
  real s(0:NH)
  character line*80,cflag*1
  logical first
  complex cx(0:NH)
  equivalence (x,cx)
  data n/0/,k0/9999999/,first/.true./
  save n,k0,w,first,pi,fs,xi

  if(first) then
     pi=4.0*atan(1.0)
     fs=12000.0
     do i=0,NFFT-1
        ww=sin(i*pi/NFFT)
        w(i)=ww*ww/NFFT
        xi(i)=2.0*pi*i
     enddo
     first=.false.
  endif

  if(k.lt.NFFT) go to 900
  if(k.lt.k0) n=0
  k0=k

  x=w*id2(k-NFFT:k-1)              !Apply window
  call four2a(x,NFFT,1,-1,0)       !Compute spectrum, r2c
  df=fs/NFFT
  if (ntol.gt.noffset) then
     ia=0
     ib=nint((noffset*2)/df)
  else
     ia=nint((noffset-ntol)/df)
     ib=nint((noffset+ntol)/df)
  endif
  smax=0.
  s=0.
  do i=ia,ib
     s(i)=real(cx(i))**2 + aimag(cx(i))**2
     if(s(i).gt.smax) then
        smax=s(i)
        ipk=i
     endif
  enddo

  call peakup(s(ipk-1),s(ipk),s(ipk+1),dx)
  fpeak=df * (ipk+dx)
  ap=(fpeak/fs+1.0/(2.0*NFFT))
  an=(fpeak/fs-1.0/(2.0*NFFT))
  sp=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*ap),-sin(xi*ap)))
  sn=sum(id2((k-NFFT):k-1)*cmplx(cos(xi*an),-sin(xi*an)))
  fpeak=fpeak+fs*(abs(sp)-abs(sn))/(abs(sp)+abs(sn))/(2*NFFT)
  xsum=0.
  nsum=0
  do i=ia,ib
     if(abs(i-ipk).gt.10) then
        xsum=xsum+s(i)
        nsum=nsum+1
     endif
  enddo
  ave=xsum/nsum
  snr=db(smax/ave)
  pave=db(ave) + 8.0
  cflag=' '
  if(snr.lt.20.0) cflag='*'
  n=n+1
  nsec=mod(time(),86400)
  nhr=nsec/3600
  nmin=mod(nsec/60,60)
  nsec=mod(nsec,60)
  ncal=1
  ferr=fpeak-noffset
  write(line,1100)  nhr,nmin,nsec,nkhz,ncal,noffset,fpeak,ferr,pave,   &
          snr,cflag,char(0)
1100 format(i2.2,':',i2.2,':',i2.2,i7,i3,i6,2f10.3,2f7.1,2x,a1,a1)

900 return
end subroutine freqcal


code:
# Compilers
CC = gcc
CXX = g++
FC = gfortran

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

# Default rules
%.o: %.c
	${CC} ${CFLAGS} -c $<
%.o: %.f
	${FC} ${FFLAGS} -c $<
%.o: %.F
	${FC} ${FFLAGS} -c $<
%.o: %.f90
	${FC} ${FFLAGS} -c $<
%.o: %.F90
	${FC} ${FFLAGS} -c $<

all:	wsprlf

OBJS0 = testpsk.o four2a.o bpfilter.o nonlinear.o tweak1.o spectrum.o smo.o
testpsk: $(OBJS0)
	$(FC) -o testpsk $(OBJS0) -lfftw3f

OBJS1 = gmsk8.o four2a.o gaussfilt.o
gmsk8: $(OBJS1)
	$(FC) -o gmsk8 $(OBJS1) -lfftw3f

OBJS2 = testfsk.o four2a.o smo.o
testfsk: $(OBJS2)
	$(FC) -o testfsk $(OBJS2) -lfftw3f

OBJS3 = fsk2sim.o four2a.o smo.o wavhdr.o gran.o
fsk2sim: $(OBJS3)
	$(FC) -o fsk2sim $(OBJS3) -lfftw3f

OBJS4 = fsk4sim.o four2a.o wavhdr.o gran.o tweak1.o
fsk4sim: $(OBJS4)
	$(FC) -o fsk4sim $(OBJS4) -lfftw3f

OBJS5 = wsprlf.o four2a.o
wsprlf:	$(OBJS5)
	$(FC) -o wsprlf $(OBJS5) -lfftw3f

.PHONY : clean

clean:
	$(RM) *.o testpsk testfsk fsk2sim fsk4sim wsprlf


code:
# Compilers
CC = gcc
CXX = g++
FC = gfortran

FFLAGS = -O2 -fbounds-check -Wall -Wno-conversion
CFLAGS = -O2 -I.

# Default rules
%.o: %.c
	${CC} ${CFLAGS} -c $<
%.o: %.f
	${FC} ${FFLAGS} -c $<
%.o: %.F
	${FC} ${FFLAGS} -c $<
%.o: %.f90
	${FC} ${FFLAGS} -c $<
%.o: %.F90
	${FC} ${FFLAGS} -c $<

all:	dbpsksim.exe

OBJS0 = testpsk.o four2a.o bpfilter.o nonlinear.o tweak1.o spectrum.o smo.o
testpsk: $(OBJS0)
	$(FC) -o testpsk $(OBJS0) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS1 = gmsk8.o four2a.o gaussfilt.o
gmsk8: $(OBJS1)
	$(FC) -o gmsk8 $(OBJS1) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS2 = testfsk.o four2a.o smo.o
testfsk: $(OBJS2)
	$(FC) -o testfsk $(OBJS2) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS3 = fsk2sim.o four2a.o smo.o wavhdr.o gran.o
fsk2sim: $(OBJS3)
	$(FC) -o fsk2sim $(OBJS3) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS4 = fsk4sim.o four2a.o gran.o genfsk4.o smo.o getsnr.o spec4.o \
      watterson.o db.o snr2_wsprlf.o pctile.o shell.o snr_wsprlf.o
fsk4sim.exe: $(OBJS4)
	$(FC) -o fsk4sim.exe $(OBJS4) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS5 = wsprlf.o four2a.o downsample.o
wsprlf.exe:	$(OBJS5)
	$(FC) -o wsprlf.exe $(OBJS5) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS6 = wspr_gmsk.o four2a.o gaussfilt.o
wspr_gmsk.exe: $(OBJS6)
	$(FC) -o wspr_gmsk.exe $(OBJS6) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS7 = wspr_msk.o four2a.o bpfilter.o
wspr_msk.exe: $(OBJS7)
	$(FC) -o wspr_msk.exe $(OBJS7) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS8 = dbpsksim.o four2a.o gran.o genbpsk.o watterson.o db.o \
      encode120.o bpdecode120.o platanh.o
dbpsksim.exe: $(OBJS8)
	$(FC) -o dbpsksim.exe $(OBJS8) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS9 = fsk4a.o four2a.o gran.o genfsk4a.o spec4.o \
      watterson.o db.o
fsk4a.exe: $(OBJS9)
	$(FC) -o fsk4a.exe $(OBJS9) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS10 = gmsk8.o gaussfilt.o four2a.o
gmsk8.exe: $(OBJS10)
	   $(FC) -o gmsk8.exe $(OBJS10) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS11 = gmsksim.o four2a.o gran.o gengmsk.o genbpsk.o watterson.o db.o \
      encode168.o bpdecode168.o platanh.o gaussfilt.o tweak1.o smo121.o
gmsksim.exe: $(OBJS11)
	$(FC) -o gmsksim.exe $(OBJS11) C:\JTSDK\fftw3f\libfftw3f-3.dll

OBJS12 = mskhfsim.o four2a.o gran.o genmskhf.o watterson.o db.o \
      encode168.o bpdecode168.o platanh.o twkfreq1.o smo121.o \
      polyfit4.o
mskhfsim.exe: $(OBJS12)
	$(FC) -o mskhfsim.exe $(OBJS12) C:\JTSDK\fftw3f\libfftw3f-3.dll

.PHONY : clean

clean:
	$(RM) *.o testpsk.exe testfsk.exe fsk2sim.exe fsk4sim.exe wsprlf.exe


code:
subroutine chkcrc10(decoded,nbadcrc)

  use crc
  integer*1 decoded(60)
  integer*1, target:: i1Dec8BitBytes(9)

! Check the CRC
  do ibyte=1,6
     itmp=0
     do ibit=1,8
        itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
     enddo
     i1Dec8BitBytes(ibyte)=itmp
  enddo
  i1Dec8BitBytes(7)=decoded(49)*128+decoded(50)*64

! Pack received CRC into bytes 8 and 9 for crc10_check
  i1Dec8BitBytes(8)=decoded(51)*2+decoded(52)
  i1Dec8BitBytes(9)=decoded(53)*128 + decoded(54)*64+decoded(55)*32 +   &
       decoded(56)*16
  i1Dec8BitBytes(9)=i1Dec8BitBytes(9) + decoded(57)*8+decoded(58)*4 +   &
       decoded(59)*2+decoded(60)*1
  nbadcrc=1
  if(crc10_check(c_loc(i1Dec8BitBytes),9)) nbadcrc=0

  return
end subroutine chkcrc10


code:
subroutine chkcrc12(decoded,nbadcrc)

  use crc
  integer*1 decoded(84)
  integer*1, target:: i1Dec8BitBytes(11)

! Check the CRC
! Collapse 84 decoded bits to 11 bytes. Bytes 1-9 are the message,
! byte 10 and first half of byte 11 is the crc
  do ibyte=1,9
      itmp=0
    do ibit=1,8
      itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
    enddo
    i1Dec8BitBytes(ibyte)=itmp
  enddo

! Pack the crc into bytes 10 and 11 for crc12_check
  i1Dec8BitBytes(10)=decoded(73)*8 + decoded(74)*4 + decoded(75)*2 + decoded(76)
  i1Dec8BitBytes(11)=decoded(77)*128 + decoded(78)*64 +                      &
       decoded(79)*32 + decoded(80)*16 + decoded(81)*8 + decoded(82)*4 +   &
       decoded(83)*2 + decoded(84)
  nbadcrc=1
  if( crc12_check(c_loc (i1Dec8BitBytes), 11) ) nbadcrc=0

  return
end subroutine chkcrc12


code:
subroutine cpolyfit(c,pp,id,maxn,aa,bb,zz,nhardsync)

  parameter (KK=84)                     !Information bits (72 + CRC12)
  parameter (ND=168)                    !Data symbols: LDPC (168,84), r=1/2
  parameter (NS=65)                     !Sync symbols (2 x 26 + Barker 13)
  parameter (NR=3)                      !Ramp up/down
  parameter (NN=NR+NS+ND)               !Total symbols (236)
  parameter (NSPS=16)                   !Samples per MSK symbol (16)
  parameter (N2=2*NSPS)                 !Samples per OQPSK symbol (32)
  parameter (N13=13*N2)                 !Samples in central sync vector (416)
  parameter (NZ=NSPS*NN)                !Samples in baseband waveform (3760)
  parameter (NFFT1=4*NSPS,NH1=NFFT1/2)

  complex c(0:NZ-1)                     !Complex waveform
  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z,z0
  real x(NS),yi(NS),yq(NS)              !For complex polyfit
  real pp(2*NSPS)                       !Shaped pulse for OQPSK
  real aa(20),bb(20)                    !Fitted polyco's
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data

  ib=NSPS-1
  ib2=N2-1
  n=0
  do j=1,117                                !First-pass demodulation
     ia=ib+1
     ib=ia+N2-1
     zz(j)=sum(pp*c(ia:ib))/NSPS
     if(abs(id(j)).eq.2) then               !Save all sync symbols
        n=n+1
        x(n)=float(ia+ib)/NZ - 1.0
        yi(n)=real(zz(j))*0.5*id(j)
        yq(n)=aimag(zz(j))*0.5*id(j)
!        write(54,1225) n,x(n),yi(n),yq(n)
!1225    format(i5,3f12.4)
     endif
     if(j.le.116) then
        zz(j+117)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS
     endif
  enddo

  aa=0.
  bb=0.
  nterms=0
  chisqa=0.
  chisqb=0.
  if(maxn.gt.0) then
     npts=n
     mode=0
     nterms=maxn
     call polyfit4(x,yi,yi,npts,nterms,mode,aa,chisqa)
     call polyfit4(x,yq,yq,npts,nterms,mode,bb,chisqb)
  endif

  nhardsync=0
  do j=1,117
     if(abs(id(j)).ne.2) cycle
     xx=j*2.0/117.0 - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=real(z)
     if(p*id(j).lt.0) nhardsync=nhardsync+1
  enddo

  return
end subroutine cpolyfit


code:
subroutine cpolyfitw(c,pp,id,maxn,aa,bb,zz,nhardsync)

  include 'wsprlf_params.f90'

  complex c(0:NZ-1)                     !Complex waveform
  complex zz(NS+ND)                     !Complex symbol values (intermediate)
  complex z,z0
  real x(NS),yi(NS),yq(NS)              !For complex polyfit
  real pp(2*NSPS)                       !Shaped pulse for OQPSK
  real aa(20),bb(20)                    !Fitted polyco's
  integer id(NS+ND)                     !NRZ values (+/-1) for Sync and Data

  ib=NSPS-1
  ib2=N2-1
  n=0
  jz=(NS+ND+1)/2
  do j=1,jz                                !First-pass demodulation
     ia=ib+1
     ib=ia+N2-1
     zz(j)=sum(pp*c(ia:ib))/NSPS
     if(abs(id(j)).eq.2) then               !Save all sync symbols
        n=n+1
        x(n)=float(ia+ib)/NZ - 1.0
        yi(n)=real(zz(j))*0.5*id(j)
        yq(n)=aimag(zz(j))*0.5*id(j)
!        write(54,1225) n,x(n),yi(n),yq(n)
!1225    format(i5,3f12.4)
     endif
     if(j.lt.jz) then
        zz(j+jz)=sum(pp*c(ia+NSPS:ib+NSPS))/NSPS
     endif
  enddo

  aa=0.
  bb=0.
  nterms=0
  chisqa=0.
  chisqb=0.
  if(maxn.gt.0) then
     npts=n
     mode=0
     nterms=maxn
     call polyfit4(x,yi,yi,npts,nterms,mode,aa,chisqa)
     call polyfit4(x,yq,yq,npts,nterms,mode,bb,chisqb)
  endif

  nhardsync=0
  do j=1,205
     if(abs(id(j)).ne.2) cycle
     xx=j*2.0/205.0 - 1.0
     yii=1.
     yqq=0.
     if(nterms.gt.0) then
        yii=aa(1)
        yqq=bb(1)
        do i=2,nterms
           yii=yii + aa(i)*xx**(i-1)
           yqq=yqq + bb(i)*xx**(i-1)
        enddo
     endif
     z0=cmplx(yii,yqq)
     z=zz(j)*conjg(z0)
     p=real(z)
     if(p*id(j).lt.0) nhardsync=nhardsync+1
  enddo

  return
end subroutine cpolyfitw