#include "NR_F3Tensor.h"
#include <stdlib.h>
extern "C" {
#include "nrutil.h"
}
#include "Msg.h"

// Normally in nr.h
extern "C" {
void rlft3(float ***data, float **speq, unsigned long nn1,
	   unsigned long nn2, unsigned long nn3, int isign);
}

//
NR_F3Tensor::NR_F3Tensor()
{
  data = NULL;
  spec = NULL;
  domain = Spatial;
  xsize=0;
  ysize=0;
  
}

NR_F3Tensor::~NR_F3Tensor()
{
  if (data!=NULL)
    free_f3tensor(data,1,1,1,ysize,1,xsize);
  if (spec!=NULL)
    free_matrix(spec,1,1,1,2*ysize);
}

void NR_F3Tensor::fft()
{
  if (domain!=Spatial)
    msg.err("NR_F3Tensor::fft() - Not Spatial");
  
  rlft3(data,spec,1,ysize,xsize,1);
  domain=Frequency;
  
}

void NR_F3Tensor::ifft()
{
  if (domain!=Frequency)
    msg.err("NR_F3Tensor::ifft() - Not Frequency");
  
  rlft3(data,spec,1,ysize,xsize,-1);
  domain=Spatial;
  
}

void NR_F3Tensor::ifftScaled()
{
  // IFFT
  ifft();
  
  // Fix up a constant scaling
  float constant = 1.0/(xsize*ysize/2.0);
  float **data2d = data[1];
  
  for (int  j=1 ; j<=ysize; j++)
    for (int i=1; i<=xsize; i++)
      {
	data2d[j][i] *= constant;
      }
}

void NR_F3Tensor::setSize(int x, int y)
{
  if (xsize==x && ysize==y)
    return;
  
  if (data!=NULL)
    free_f3tensor(data,1,1,1,ysize,1,xsize);
  if (spec!=NULL)
    free_matrix(spec,1,1,1,2*ysize);
  
  xsize=x;
  ysize=y;
  
  data = f3tensor(1,1,1,ysize,1,xsize);
  spec = matrix(1,1,1,2*ysize);
}