/******************************************************************************

   This file is a part of the

             COMPLEX PSEUDO-RAW SAR IMAGE PROCESSING TOOLBOX

   Copyright (C) 2017 Rémy Abergel (Remy.Abergel AT parisdescartes.fr), Loïc
   Denis (Loic.Denis AT univ-st-etienne.fr), Saïd Ladjal (Said.Ladjal AT
   telecom-paristech.fr) and Florence Tupin (Florence.Tupin AT
   telecom-paristech.fr).

   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, see <http://www.gnu.org/licenses/>.

******************************************************************************/

#include <math.h>
#include <fftw3.h>

/* horizontal translation using Shannon interpolation */
void horizontal_translation(dft_out,dft_in,nx,ny,tx,iplan_out)
     fftw_complex *dft_in,*dft_out;
     fftw_plan iplan_out;
     int nx,ny;
     double tx;
{
  int a,b,adr;
  double c,s,re,im,nrm,coef;

  nrm = 1./((double)nx*(double)ny);
  coef = -2.*M_PI*tx/(double)nx;

  // fill dft_out //
  for(a=-(nx/2);a<=(nx-1)/2;a++) {
    c = cos(coef*a);
    s = sin(coef*a);
    for(b=-(ny/2);b<=(ny-1)/2;b++) {
      adr = ((ny+b)%ny)*nx+(nx+a)%nx;
      re = dft_in[adr][0];
      im = dft_in[adr][1];
      dft_out[adr][0] = (re*c-im*s) * nrm;
      dft_out[adr][1] = (re*s+im*c) * nrm;
    }
  }

  // apply inverse DFT //
  fftw_execute(iplan_out);

}


/* vertical translation using Shannon interpolation */
void vertical_translation(dft_out,dft_in,nx,ny,ty,iplan_out)
     fftw_complex *dft_in,*dft_out;
     fftw_plan iplan_out;
     int nx,ny;
     double ty;
{
  int a,b,adr;
  double c,s,re,im,nrm,coef;

  nrm = 1./((double)nx*(double)ny);
  coef = -2.*M_PI*ty/(double)ny;

  // fill dft_out //
  for(b=-(ny/2);b<=(ny-1)/2;b++) {
    c = cos(coef*b);
    s = sin(coef*b);
    for(a=-(nx/2);a<=(nx-1)/2;a++) {
      adr = ((ny+b)%ny)*nx+(nx+a)%nx;
      re = dft_in[adr][0];
      im = dft_in[adr][1];
      dft_out[adr][0] = (re*c-im*s) * nrm;
      dft_out[adr][1] = (re*s+im*c) * nrm;
    }
  }

  // apply inverse DFT //
  fftw_execute(iplan_out);

}

/* cardinal sine function */
double sinc(double x)
{
  double px = M_PI*x;
  return (x==0) ? 1. : sin(px)/px;
}
