# TODO:
# threading? change to pp_defs
# docs improvement
# use PDL::Complex?

require 'typespec';

pp_addpm({At=>'Top'},<<'EOD');
=head1 NAME

PDL::FFTW - PDL interface to the Fastest Fourier Transform in the West v2.x

=head1 DESCRIPTION

This is a means to interface PDL with the FFTW library. It's similar to
the standard FFT routine but it's usually faster and has support for
real transforms. It works well for the types of PDLs for
which was the library was compiled (otherwise it must do conversions).

=head1 SYNOPSIS

    use PDL::FFTW

    load_wisdom("file_name");

    out_cplx_pdl =  fftw(in_cplx_pdl);
    out_cplx_pdl = ifftw(in_cplx_pdl);

    out_cplx_pdl =  rfftw(in_real_pdl);
    out_real_pdl = irfftw(in_cplx_pdl);

    cplx_pdl = nfftw(cplx_pdl);
    cplx_pdl = infftw(cplx_pdl);

    cplx_pdl = Cmul(a_cplx_pdl, b_cplx_pdl);
    cplx_pdl = Cconj(a_cplx_pdl);
    real_pdl = Cmod(a_cplx_pdl); 
    real_pdl = Cmod2(a_cplx_pdl); 

=head1 FFTW documentation

Please refer to the FFTW documentation for a better understanding of
these functions.

Note that complex numbers are represented as piddles with leading
dimension size 2 (real/imaginary pairs).

=cut
EOD

pp_addpm(<<'EOD');

=head2 load_wisdom

=for ref

Loads the wisdom from a file for better FFTW performance.

The wisdom
is automatically saved when the program ends. It will be automagically
called when the variable C<$PDL::FFT::wisdom> is set to a file name.
For example, the following is a useful idiom to have in your F<.perldlrc>
file:

  $PDL::FFT::wisdom = "$ENV{HOME}/.fftwisdom"; # save fftw wisdom in this file

Explicit usage:

=for usage

  load_wisdom($fname);

=head2 fftw

=for ref 

calculate the complex FFT of a real piddle (complex input, complex output)

=for usage

   $pdl_cplx = fftw $pdl_cplx;

=head2 ifftw

=for ref 

Complex inverse FFT (complex input, complex output).

=for usage

   $pdl_cplx = ifftw $pdl_cplx;

=head2 nfftw

=for ref 

Complex inplace FFT (complex input, complex output).

=for usage

   $pdl_cplx = nfftw $pdl_cplx;

=head2 infftw

=for ref 

Complex inplace inverse FFT (complex input, complex output).

=for usage

   $pdl_cplx = infftw $pdl_cplx;

=head2 rfftw

=for ref 

Real FFT. For an input piddle of dimensions [n1,n2,...] the 
output is [2,(n1/2)+1,n2,...] (real input, complex output).

=for usage

  $pdl_cplx = fftw $pdl_real;

=head2 irfftw

=for ref 

Real inverse FFT. Have a look at rfftw to understand the format. USE
ONLY an even n1! (complex input, real output)

=for usage

  $pdl_real = ifftw $pdl_cplx;

=head2 nrfftw

=for ref 

Real inplace FFT. If you want a transformation on a piddle
with dimensions [n1,n2,....] you MUST pass in a piddle
with dimensions [2*(n1/2+1),n2,...] (real input, complex output).

Use with care due to dimension restrictions mentioned below.
For details check the html docs that come with the fftw library.

=for usage

  $pdl_cplx = nrfftw $pdl_real;

=head2 inrfftw

=for ref 

Real inplace inverse FFT. Have a look at nrfftw to understand the format. USE
ONLY an even first dimension size! (complex input, real output)

=for usage

  $pdl_real = infftw $pdl_cplx;

=cut

EOD


pp_addpm('

use PDL;
use PDL::ImageND qw/kernctr/;
use strict;

my ($wisdom_fname, $wisdom_loaded, %plan_cache, $datatype, $warn_conv);

#
# package vriables: $wisdom_fname, $wisdom_loaded, %plan_cache, $datatype, $warn_conv, 
#           $COMPILED_TYPE: (Type "double" or "float") that PDL::FFTW was compiled/linked to
#

$wisdom_loaded = 0;

$PDL::FFTW::COMPILED_TYPE = "'.$fftwtype.'";

sub load_wisdom {
    $wisdom_fname = shift;
    
    if (!$wisdom_loaded) {
	my ($wisdom, $ok);
	$wisdom_loaded = 1;
	if (!open(TMPFILE,$wisdom_fname)) {
	    print STDERR "Warning: couldn\'t find wisdom file\n";
	    return;
	}
	$wisdom = <TMPFILE>;
	close(TMPFILE);
	$ok = PDL_fftw_import_wisdom_from_string($wisdom);
	if (!$ok) {
	    print STDERR "Warning: couldn\'t import wisdom\n";
	}
    }
}

sub save_wisdom {
  if ($wisdom_loaded) {
    my $wisdom = PDL_fftw_export_wisdom_to_string();

    if (!open(TMPFILE,">$wisdom_fname")) {
      print STDERR "Warning: couldn\'t write wisdom file\n";
      return;
    }
    print TMPFILE $wisdom;
    close(TMPFILE);
  }
}

sub BEGIN {
  my $HOME=$ENV{HOME};
  my $test_type = ' . $fftwtype  . ' 1;

  $datatype = $test_type->get_datatype;
  $warn_conv = 0;
#  if ( -e "$HOME/.wisdom" ) {
#    print STDERR "Autoloading wisdom\n";
#    load_wisdom("$HOME/.wisdom");
#  }
}

sub END {
  save_wisdom();
}

if (defined $PDL::FFT::wisdom) {
    load_wisdom($PDL::FFT::wisdom);
    print "loading wisdom from $PDL::FFT::wisdom\n" if $PDL::debug;
}

# real fftw
*PDL::rfftw = \&rfftw;
sub rfftw {
  my $in = PDL::Core::topdl(shift);
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i);

 
  if ($datatype != $in->get_datatype) {
    $in = convert $in, ' . $fftwtype . ';
    if ($warn_conv != 1) {
      print STDERR "PDL::FFTW Warning! doing conversion.\n";
      $warn_conv = 1;
    }
  }
  foreach $i (@dims) {
    unshift @newdims, $i;
  }
  $name="r_@{newdims}_f";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_rfftwnd_create_plan($pdldims,0,$wisdom_loaded);
  }
  $dims[0]= int($dims[0]/2)+1;
  unshift @dims, 2;
  $out = zeroes ' . $fftwtype . ',@dims;
  $in->make_physical;
  PDL_rfftwnd_one_real_to_complex($plan_cache{$name},$in,$out);
  return $out;
}

# this assumes even last input dimension!

*PDL::irfftw = \&irfftw;
sub irfftw {
  my $in = PDL::Core::topdl(shift);
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i);

  $i = shift @dims;
  if ($i != 2) {
     barf("Working only on complex numbers!");
  } 
  if ($datatype != $in->get_datatype) {
    $in = convert $in, ' . $fftwtype . ';
    if ($warn_conv != 1) {
      print STDERR "PDL::FFTW Warning! doing conversion.\n";
      $warn_conv = 1;
    }
  }
  $i= shift @dims;
  unshift @newdims, 2*($i-1);
  foreach $i (@dims) {
    unshift @newdims, $i;
  } 
  $name="r_@{newdims}_b";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_rfftwnd_create_plan($pdldims,1,$wisdom_loaded);
  }
  unshift @dims, 2*($i-1);
  $out = zeroes ' . $fftwtype . ',@dims;
  $in->make_physical;
  PDL_rfftwnd_one_complex_to_real($plan_cache{$name},$in,$out);
  return $out;
}

# real inplace fftw

*PDL::nrfftw = \&nrfftw;
sub nrfftw {
  my ($in)=@_;
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i,$ndim);

 
  if ($datatype != $in->get_datatype) {
	barf("Cannot do inplace fftw: compiled for ' . $fftwtype . ' !");
  }
  foreach $i (@dims) {
    unshift @newdims, $i;
  }
  $ndim = $newdims[$#newdims] = 2*(int($newdims[$#newdims]/2) - 1);
  $name="nr_@{newdims}_f";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_rfftwnd_create_plan($pdldims,0,$wisdom_loaded+2);
  }
  $dims[0]= int($ndim/2)+1;
  unshift @dims, 2;
  $in->make_physical;
  PDL_inplace_rfftwnd_one_real_to_complex($plan_cache{$name},$in);
  $in->reshape(@dims);
  return $in;
}

# this assumes even last input dimension!

*PDL::inrfftw = \&inrfftw;
sub inrfftw {
  my ($in)=@_;
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i,$ndim);

  $i = shift @dims;
  if ($i != 2) {
     barf("Working only on complex numbers!");
  } 
  if ($datatype != $in->get_datatype) {
	barf("Cannot do inplace fftw: compiled for ' . $fftwtype . ' !");
  }
  $i= shift @dims;
  $ndim =  2*($i-1);
  unshift @newdims, $ndim;
  foreach $i (@dims) {
    unshift @newdims, $i;
  } 
  $name="nr_@{newdims}_b";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_rfftwnd_create_plan($pdldims,1,$wisdom_loaded+2);
  }
  unshift @dims, 2*(int($ndim/2)+1);
  $in->make_physical;
  PDL_inplace_rfftwnd_one_complex_to_real($plan_cache{$name},$in);
  $in->reshape(@dims);
  return $in;
}

# complex fftw

*PDL::fftw = \&fftw;
sub fftw {
  my $in = PDL::Core::topdl(shift);
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i);

  if ($datatype != $in->get_datatype) {
    $in = convert $in, ' . $fftwtype . ';
    if ($warn_conv != 1) {
      print STDERR "PDL::FFTW Warning! doing conversion.\n";
      $warn_conv = 1;
    }
  }
  $i=shift @dims;
  if ($i!=2) {
    barf("It works only on complex!");
  }
  foreach $i (@dims) {
    unshift @newdims, $i;
  }
  $name="c_@{newdims}_f";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_fftwnd_create_plan($pdldims,0,$wisdom_loaded);
  }
  unshift @dims, 2;
  $out = zeroes ' . $fftwtype . ',@dims;
  $in->make_physical;
  PDL_fftwnd_one($plan_cache{$name},$in,$out);
  return $out;
}

*PDL::ifftw = \&ifftw;
sub ifftw {
  my $in = PDL::Core::topdl(shift);
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i);

  if ($datatype != $in->get_datatype) {
    $in = convert $in, ' . $fftwtype . ';
    if ($warn_conv != 1) {
      print STDERR "PDL::FFTW Warning! doing conversion.\n";
      $warn_conv = 1;
    }
  } 
  $i = shift @dims;
  if ($i != 2) {
     barf("Working only on complex numbers!");
  }
  foreach $i (@dims) {
    unshift @newdims, $i;
  } 
  $name="c_@{newdims}_b";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_fftwnd_create_plan($pdldims,1,$wisdom_loaded);
  }
  unshift @dims, 2;
  $out = zeroes ' . $fftwtype . ',@dims;
  $in->make_physical;
  PDL_fftwnd_one($plan_cache{$name},$in,$out);
  return $out;
}

# complex inplace fftw

*PDL::nfftw = \&nfftw;
sub nfftw {
  my ($in)=@_;
  my @dims = $in->dims;
  my @newdims;
  my ($name,$i);

  if ($datatype != $in->get_datatype) {
	barf("Cannot do inplace fftw: compiled for ' . $fftwtype . ' !");
  } 
  $i=shift @dims;
  if ($i!=2) {
    barf("It works only on complex!");
  }
  foreach $i (@dims) {
    unshift @newdims, $i;
  }
  $name="n_@{newdims}_f";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_fftwnd_create_plan($pdldims,0,$wisdom_loaded+2);
  }
  $in->make_physical;
  PDL_inplace_fftwnd_one($plan_cache{$name},$in);
  return $in;
}

*PDL::infftw = \&infftw;
sub infftw {
  my ($in)=@_;
  my @dims = $in->dims;
  my @newdims;
  my ($out,$name,$i);

  if ($datatype != $in->get_datatype) {
	barf("Cannot do inplace fftw: compiled for ' . $fftwtype . q| !");
  } 
  $i = shift @dims;
  if ($i != 2) {
     barf("Working only on complex numbers!");
  } 
  foreach $i (@dims) {
    unshift @newdims, $i;
  } 
  $name="n_@{newdims}_b";
  if ( !defined($plan_cache{$name}) ) {
    my $pdldims = long [@newdims];
    $plan_cache{$name} = PDL_fftwnd_create_plan($pdldims,1,$wisdom_loaded+2);
  }
  $in->make_physical;
  PDL_inplace_fftwnd_one($plan_cache{$name},$in);
  return $in;
}

=head2 rfftwconv

=for ref

ND convolution using real ffts from the FFTW library

=for example

  $conv = rfftwconv $im, kernctr $im, $k; # kernctr is from PDL::FFT

=cut

*PDL::rfftwconv = \&rfftwconv;
sub rfftwconv {
  my ($a,$b) = @_;
  my ($ca,$cb) = (rfftw($a), rfftw($b));
  my $cmul = Cmul($ca,$cb);
  my $ret = irfftw($cmul);
  my $mul = 1; for (0..$a->getndims-1) {$mul *= $a->getdim($_)}
  $ret /= $mul;
  return $ret;
}

=head2 fftwconv

=for ref

ND convolution using complex ffts from the FFTW library

Assumes real input!

=for example

  $conv = fftwconv $im, kernctr $im, $k; # kernctr is from PDL::FFT

=cut

*PDL::fftwconv = \&fftconv;
sub fftwconv {
  my ($a,$b) = @_;
  my ($ca,$cb) = (PDL->zeroes(2,$a->dims),PDL->zeroes(2,$b->dims));
  $ca->slice('(0)') .= $a;
  $cb->slice('(0)') .= $b;
  nfftw($ca); nfftw($cb);
  my $cmul = Cmul($ca,$cb);
  infftw $cmul;  # transfer back inplace
  # and return real part only
  my $ret = $cmul->slice('(0)')->sever;
  my $mul = 1; for (0..$a->getndims-1) {$mul *= $a->getdim($_)}
  $ret /= $mul;
  return $ret;
}

|);

pp_def('Cmul',
       Pars => 'a(n); b(n); [o]c(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c(n => 0)= $a(n => 0)*$b(n => 0) - $a(n => 1)*$b(n => 1);
$c(n => 1)= $a(n => 0)*$b(n => 1) + $a(n => 1)*$b(n => 0);
%}
',
	Doc => <<'EOD',
=for ref

Complex multiplication

=for usage

   $out_pdl_cplx = Cmul($a_pdl_cplx,$b_pdl_cplx);

EOD
      );

pp_def('Cscale',
       Pars => 'a(n); b(); [o]c(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c(n => 0)= $a(n => 0)*$b();
$c(n => 1)= $a(n => 1)*$b();
%}
',
	Doc => <<'EOD',
=for ref

Complex by real multiplation.

=for usage

  Cscale($a_pdl_cplx,$b_pdl_real);
EOD
      );

pp_def('Cdiv',
       Pars => 'a(n); b(n); [o]c(n);',
       Code => '
$GENERIC() divi;

if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
divi = $b(n => 0)*$b(n => 0) + $b(n => 1)*$b(n => 1);

$c(n => 0)= ( $a(n => 0)*$b(n => 0) + $a(n => 1)*$b(n => 1) ) / divi;
$c(n => 1)= ( $a(n => 1)*$b(n => 0) - $a(n => 0)*$b(n => 1) ) / divi;
%}
',
	Doc => <<'EOD',
=for ref

Complex division.

=for usage

  $out_pdl_cplx = Cdiv($a_pdl_cplx,$b_pdl_cplx);

EOD
      );

pp_def('Cbmul',
       Pars => 'a(n); b(n);',
       Code => '
$GENERIC() tmp1,tmp;

if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
tmp = $a(n => 0);
tmp1 = $a(n => 1);
$a(n => 0)= tmp * $b(n => 0) - tmp1 * $b(n => 1);
$a(n => 1)= tmp * $b(n => 1) + tmp1 * $b(n => 0);
%}
',
	Doc => <<'EOD',
=for ref

Complex inplace multiplication.

=for usage

   Cbmul($a_pdl_cplx,$b_pdl_cplx);
EOD
      );

pp_def('Cbscale',
       Pars => 'a(n); b();',
       Code => '
$GENERIC() tmp1,tmp,divi;

if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$a(n => 0) *= $b();
$a(n => 1) *= $b();
%}
',
	Doc => <<'EOD',
=for ref

Complex inplace multiplaction by real.

=for usage

   Cbscale($a_pdl_cplx,$b_pdl_real);
EOD
      );

pp_def('Cconj',
       Pars => 'a(n); [o]c(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c(n => 0)= $a(n => 0);
$c(n => 1)= -$a(n => 1);
%}
',
	Doc => <<'EOD',
=for ref

Complex conjugate.

=for usage

   $out_pdl_cplx = Cconj($a_pdl_cplx);
EOD
      );

pp_def('Cbconj',
       Pars => 'a(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$a(n => 1)= -$a(n => 1);
%}
',
	Doc => <<'EOD',
=for ref

Complex inplace conjugate.

=for usage

   Cbconj($a_pdl_cplx);
EOD
      );

pp_def('Cexp',
       Pars => 'a(n); [o]c(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c(n => 0)= exp($a(n => 0)) * cos($a(n => 1));
$c(n => 1)= exp($a(n => 0)) * sin($a(n => 1));
%}
',
	Doc => <<'EOD',
=for ref

Complex exponentation.

=for usage

   $out_pdl_cplx = Cexp($a_pdl_cplx);
EOD
      );

pp_def('Cbexp',
       Pars => 'a(n);',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
double re = $a(n => 0);
double im = $a(n => 1);

$a(n => 0)= exp(re) * cos(im);
$a(n => 1)= exp(re) * sin(im);
%}
',
	Doc => <<'EOD',
=for ref

Complex inplace exponentation.

=for usage

   $out_pdl_cplx = Cexp($a_pdl_cplx);
EOD
      );

pp_def('Cmod',
       Pars => 'a(n); [o]c();',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c()= sqrt ( $a(n => 0)*$a(n => 0) + $a(n => 1)*$a(n => 1) );
%}
',
	Doc => <<'EOD',
=for ref

modulus of a complex piddle.

=for usage

  $out_pdl_real = Cmod($a_pdl_cplx);
EOD
      );

pp_def('Carg',
       Pars => 'a(n); [o]c();',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c()= atan2($a(n=>1),$a(n=>0));
%}
',
	Doc => <<'EOD',
=for ref

argument of a complex number.

=for usage

   $out_pdl_real = Carg($a_pdl_cplx);
EOD
      );

pp_def('Cmod2',
       Pars => 'a(n); [o]c();',
       Code => '
if ($SIZE(n)!=2) barf("This function works only on complex\n");
threadloop %{
$c()= ( $a(n => 0)*$a(n => 0) + $a(n => 1)*$a(n => 1) );
%}
',
	Doc => <<'EOD',
=for ref

Returns squared modulus of a complex number.

=for usage

   $out_pdl_real = Cmod2($a_pdl_cplx);
EOD
      );


### real fftw

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int
PDL_rfftwnd_create_plan(dims, dir, flag)
  pdl* dims
  int dir
  int flag
CODE:
 {
  fftw_direction fdir=0;
  int fflag=FFTW_USE_WISDOM;

  if (dims->ndims != 1) {barf("Only 1d input dimesions make sense");} 
  if (dims->data == NULL) {barf("input piddles must be physical");} 
  if (dims->datatype != PDL_L) {barf("Only integers please");} 

  if (dir) {
    fdir=FFTW_COMPLEX_TO_REAL;
  }
  else {
    fdir=FFTW_REAL_TO_COMPLEX;
  }
  if (flag & 1 ) { 
    fflag |= FFTW_ESTIMATE;
  }
  else {
    fflag |= FFTW_MEASURE;
  }
  if (flag & 2 ) { 
    fflag |= FFTW_IN_PLACE;
  }
  else {
    fflag |= FFTW_OUT_OF_PLACE;
  }
 
  RETVAL = 
  (int) rfftwnd_create_plan( dims->dims[0], 
			    ( int *) dims->data,
			    fdir,
			    fflag);

 }
OUTPUT:
 RETVAL
'
);

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_rfftwnd_one_real_to_complex(plan, in, out)
  int plan
  pdl* in
  pdl* out
 CODE:
  if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad Type");}
  rfftwnd_one_real_to_complex( (rfftwnd_plan) plan, (fftw_real *) in->data, (fftw_complex *) out->data);
');

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_rfftwnd_one_complex_to_real(plan, in, out)
  int plan
  pdl* in
  pdl* out
 CODE:
  if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad type");} 
  rfftwnd_one_complex_to_real( (rfftwnd_plan) plan, (fftw_complex *) in->data, (fftw_real *) out->data);


');

# NOTE: BUG! the inplace code below will not work with slices!
# no backpropagation of results!
pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_inplace_rfftwnd_one_real_to_complex(plan, in)
  int plan
  pdl* in
 CODE:
  if (in->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE) {barf("Bad Type");}
  PDL->children_changesoon(in, PDL_PARENTDATACHANGED);
  rfftwnd_one_real_to_complex( (rfftwnd_plan) plan, (fftw_real *) in->data, NULL);
  /* this call is crucial to propagate changes back if slices are given as arguments
   * Note: must not used vaffinechanged (!) since any slice has physical data
   */
  PDL->changed( in , PDL_PARENTDATACHANGED , 0 );
');

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_inplace_rfftwnd_one_complex_to_real(plan, in)
  int plan
  pdl* in
 CODE:
  if (in->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE) {barf("Bad type");} 
  PDL->children_changesoon(in, PDL_PARENTDATACHANGED);
  rfftwnd_one_complex_to_real( (rfftwnd_plan) plan, (fftw_complex *) in->data, NULL);
  /* this call is crucial to propagate changes back if slices are given as arguments
   * Note: must not used vaffinechanged (!) since any slice has physical data
   */
  PDL->changed( in , PDL_PARENTDATACHANGED , 0 );

');

### complex fftw

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int
PDL_fftwnd_create_plan(dims, dir, flag)
  pdl* dims
  int dir
  int flag
CODE:
 {
  fftw_direction fdir=0;
  int fflag=FFTW_USE_WISDOM;

  if (dims->ndims != 1) {barf("Only 1d input dimesions make sense");} 
  if (dims->data == NULL) {barf("input piddles must be physical");} 
  if (dims->datatype != PDL_L) {barf("Only integers please");} 

  if (dir) {
    fdir=FFTW_BACKWARD;
  }
  else {
    fdir=FFTW_FORWARD;
  }
  if (flag & 1 ) { 
    fflag |= FFTW_ESTIMATE;
  }
  else {
    fflag |= FFTW_MEASURE;
  }
  if (flag & 2 ) { 
    fflag |= FFTW_IN_PLACE;
  }
  else {
    fflag |= FFTW_OUT_OF_PLACE;
  }
 
  RETVAL = 
  (int) fftwnd_create_plan( dims->dims[0], 
			    ( int *) dims->data,
			    fdir,
			    fflag);

 }
OUTPUT:
 RETVAL
'
);

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_fftwnd_one(plan, in, out)
  int plan
  pdl* in
  pdl* out
 CODE:
  if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad type!");} 
  fftwnd_one( (fftwnd_plan) plan, (fftw_complex *) in->data, (fftw_complex *) out->data);
');

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void
PDL_inplace_fftwnd_one(plan, in)
  int plan
  pdl* in
 CODE:
  if (in->data==NULL) {barf("Need a physical pdl!");}
  if (in->datatype != PDL_MYTYPE) {barf("Only float please");} 
  PDL->children_changesoon(in, PDL_PARENTDATACHANGED);
  fftwnd_one( (fftwnd_plan) plan, (fftw_complex *) in->data, NULL);
  /* this call is crucial to propagate changes back if slices are given as arguments
   * Note: must not used vaffinechanged (!) since any slice has physical data
   */
  PDL->changed( in , PDL_PARENTDATACHANGED , 0 );
');

### wisdom stuff

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int
PDL_fftw_import_wisdom_from_string (wisdom)
 char* wisdom
 CODE:
  RETVAL = ( fftw_import_wisdom_from_string(wisdom) == FFTW_SUCCESS );
 OUTPUT:
  RETVAL
');

pp_addxs('','
MODULE = PDL::FFTW PACKAGE = PDL::FFTW

char* 
PDL_fftw_export_wisdom_to_string ()
 CODE:
  RETVAL = fftw_export_wisdom_to_string();
 OUTPUT:
  RETVAL
');


pp_add_exported('','load_wisdom save_wisdom rfftw irfftw fftw ifftw
nfftw infftw nrfftw inrfftw fftwconv rfftwconv kernctr');

# I don't see a point in exporting these (CS)
# 'PDL_rfftwnd_create_plan PDL_rfftwnd_one_real_to_complex PDL_rfftwnd_one_complex_to_real PDL_fftw_export_wisdom_to_string PDL_fftw_import_wisdom_from_string PDL_inplace_fftwnd_one PDL_fftwnd_one PDL_fftwnd_create_plan PDL_inplace_rfftwnd_one_real_to_complex PDL_inplace_rfftwnd_one_complex_to_real';

pp_addpm({At => 'Bot'},<< 'EOD');
=head1 AUTHOR

Copyright (C) 1999 Christian Pellegrin, 2000 Christian Soeller.
All rights reserved. There is no warranty. You are allowed
to redistribute this software / documentation under certain
conditions. For details, see the file COPYING in the PDL
distribution. If this file is separated from the PDL distribution,
the copyright notice should be included in the file.

=cut
EOD

pp_done();