/* This program computes cross-correlations */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "nr.h"
#include "crc.h"

int main(int argc, char*argv[])
{
  crc_par crcpar;
  fl_st fl;
  double **datar, *corar;
  int ic;
  char suffix[3]="a1", file_name[36], fnm[8];

  if (argc >= 2) strcpy(suffix,argv[1]);
  fl.in  = fopen("crc.n", "r");
  fl.dat = fopen(strcat(strcpy(file_name,"tc.col."),suffix), "r");
  fl.cor = fopen(strcat(strcpy(file_name,"crc.cor."),suffix), "w");
  fl.out = fopen(strcat(strcpy(file_name,"crc.out."),suffix), "w");
  fl.res = fopen(strcat(strcpy(file_name,"crc.res."),suffix), "w");

  read_input(&crcpar, fl);
  write_input(&crcpar, fl);

  datar  = dmatrix(1, 2, 1, crcpar.n_line_read);
  corar  = dvector(1 - crcpar.n_line_read, crcpar.n_line_read - 1);

  for (ic=1; ic<=crcpar.nicol2ar; ic++)
  {
    printf("ic=%d nicol2ar=%d\n", ic, crcpar.nicol2ar);
    crcpar.icol[2] = crcpar.icol2ar[ic];
    if ((ic > 1) && (!crcpar.smap)) fprintf(fl.cor, "  \n");
    one_cor_cal_quantify(&crcpar, datar, corar, fl);
  }

  free_dmatrix(datar, 1, 2, 1, crcpar.n_line_read);
  free_dvector(corar, 1 - crcpar.n_line_read, crcpar.n_line_read - 1);
  free_ivector(crcpar.icol2ar, 1, crcpar.nicol2ar);

  fclose(fl.in);
  fclose(fl.dat);
  fclose(fl.cor);
  fclose(fl.out);
  fclose(fl.res);
}

/* This function read the input parameters */
void read_input(crc_par *crcpar, fl_st fl)
{
  int ic;
  int fsc;
  crcpar->epsilon = 1.0e-10;

  fsc = fscanf(fl.in, "ncol=%d detcol=%c\n", &crcpar->ncol, &crcpar->detcol);

  if (crcpar->detcol == 'o')
  {
    fsc = fscanf(fl.in, "icol1=%d icol2=%d\n", &crcpar->icol[1],
    &crcpar->icol[2]);
    crcpar->nicol2ar = 1;
    crcpar->icol2ar = ivector(1, crcpar->nicol2ar);
    crcpar->icol2ar[1] = crcpar->icol[2];
  }
  else if (crcpar->detcol == 'l')
  {
    fsc = fscanf(fl.in, "icol1=%d nicol2ar=%d icol2ar=", &crcpar->icol[1],
    &crcpar->nicol2ar);
    crcpar->icol2ar = ivector(1, crcpar->nicol2ar);
    for (ic=1; ic<=crcpar->nicol2ar; ic++)
    {
      fsc = fscanf(fl.in, "%d", &crcpar->icol2ar[ic]);
      if (ic < crcpar->nicol2ar) fsc = fscanf(fl.in, " ");
    }
    fsc = fscanf(fl.in, "\n");
  }

  fsc = fscanf(fl.in, "n_line_skip=%d n_line_read=%d smap=%d\n", 
  &crcpar->n_line_skip, &crcpar->n_line_read, &crcpar->smap);
  fsc = fscanf(fl.in, "Tmax=%lf T_between_peaks=%lf\n", &crcpar->Tmax,
  &crcpar->T_between_peaks);
  fsc = fscanf(fl.in, "kpoints=%d base_intersect=%lf inter_points=%d\n",
  &crcpar->kpoints, &crcpar->base_intersect, &crcpar->inter_points);
}

/* This function writes the input parameters */
void write_input(crc_par *crcpar, fl_st fl)
{
  int ic;

  fprintf(fl.out, "ncol=%d detcol=%c\n", crcpar->ncol, crcpar->detcol);
  fprintf(fl.out, "icol1=%d nicol2ar=%d icol2ar=", crcpar->icol[1],
  crcpar->nicol2ar);
  
  for (ic=1; ic<=crcpar->nicol2ar; ic++)
  {
    fprintf(fl.out, "%d", crcpar->icol2ar[ic]);
    if (ic < crcpar->nicol2ar) fprintf(fl.out, " ");
  }
  fprintf(fl.out, "\n");

  fprintf(fl.out, "n_line_skip=%d n_line_read=%d smap=%d\n",
  crcpar->n_line_skip, crcpar->n_line_read, crcpar->smap);
  fprintf(fl.out, "Tmax=%lf T_between_peaks=%lf\n", crcpar->Tmax,
  crcpar->T_between_peaks);
  fprintf(fl.out, "kpoints=%d base_intersect=%lf inter_points=%d\n",
  crcpar->kpoints, crcpar->base_intersect, crcpar->inter_points);
  fprintf(fl.out, "\n");
  fflush(fl.out);
}

/* This function calculate the cross-correlation and quantify it  */
/* for two specific columns.                                      */
void one_cor_cal_quantify(crc_par *crcpar, double **datar, double *corar,
     fl_st fl)
{
  read_data(crcpar, datar, fl);
  calculate_correlation(crcpar, datar, 2, 2, corar, fl);
  /* amp_ph_cal(crcpar, 2, 2, corar, fl); */
  calculate_correlation(crcpar, datar, 1, 2, corar, fl);
  /* amp_ph_cal(crcpar, 1, 2, corar, fl); */
}

/* This function read the data file */
void read_data(crc_par *crcpar, double **datar, fl_st fl)
{
  double *vec, time1, time2;
  int iread, icol;
  char line[Mline], *p1;

  rewind(fl.dat);
  vec = dvector(1, crcpar->ncol);

  iread = 1;
  while (iread <= crcpar->n_line_skip) 
  {
    if (fgets(line, Mline, fl.dat) != NULL)
    {
      iread++;
    }
    else
    {
      break;
    }
  }

  if (iread != crcpar->n_line_skip + 1)
  {
    printf("iread=%d != n_line_skip=%d + 1\n", iread, crcpar->n_line_skip);
    exit(0);
  }

  iread  = 1;
  while (iread <= crcpar->n_line_read) 
  {
    if (fgets(line, Mline, fl.dat) != NULL)
    {
      p1= strtok(line, " ");
      if (iread == 1) sscanf(p1, "%lf\n", &time1);
      else if (iread == 2) 
      {
        sscanf(p1, "%lf\n", &time2);
        crcpar->deltat = time2 - time1;
        fprintf(fl.out, "deltat=%lf\n", crcpar->deltat);
      }

      sscanf(p1, "%lf\n", &vec[1]);
      for(icol=2; icol<=crcpar->ncol; icol++)
      {
        p1= strtok(NULL, " ");
        sscanf(p1, "%lf\n", &vec[icol]);
      }

      datar[1][iread] = vec[crcpar->icol[1]];
      datar[2][iread] = vec[crcpar->icol[2]];
      iread++;
    }
    else
    {
      break;
    }
  }

  if (iread != crcpar->n_line_read + 1)
  {
    printf("iread=%d != n_line_read=%d + 1\n", iread, crcpar->n_line_read);
    exit(0);
  }

  /*
  for (iread=1; iread<=crcpar->n_line_read; iread++)
    fprintf(fl.out, "%d %lf %lf\n", iread, datar[1][iread], datar[2][iread]);
  fprintf(fl.out, "\n");
  */
  free_dvector(vec, 1, crcpar->ncol);
}

/* This function calculates the cross-correlation */
void calculate_correlation(crc_par *crcpar, double **datar, int idat1,
     int idat2, double *corar, fl_st fl)
{
  double **var, *ans, V1av, V2av, xcc;
  int npt, ipt, jpt;
  int ido;

  npt = crcpar->n_line_read;
  var = dmatrix(1, 2, 1, 2*npt);
  ans = dvector(1, 4*npt);

  /*  calculating the non-connected correlation <V_1 V_2> */

  for (ipt=1; ipt<=npt; ipt++) var[1][ipt] = datar[idat1][ipt];
  for (ipt=1; ipt<=npt; ipt++) var[1][npt+ipt] = 0.0;
  for (ipt=1; ipt<=npt; ipt++) var[2][ipt] = datar[idat2][ipt];
  for (ipt=1; ipt<=npt; ipt++) var[2][npt+ipt] = 0.0;


  fprintf(fl.out, "var\n");
  for (ipt=1; ipt<=2*npt; ipt++)
    fprintf(fl.out, "%d %lf %lf\n", ipt, var[1][ipt], var[2][ipt]);
  fprintf(fl.out, "\n");

  xcorrel(var[1], var[2], (unsigned long) (2*npt), ans);

  for (ipt=1-npt; ipt<=-1; ipt++)
      corar[ipt] = ans[2*npt+ipt+1]; 
  for (ipt=0; ipt<=npt-1; ipt++)
    corar[ipt] = ans[ipt+1];

  /* calculating the average <V> */

  V1av = 0.0;
  for (ipt=1; ipt<=npt; ipt++) V1av += datar[idat1][ipt];
  V1av /= npt;

  V2av = 0.0;
  for (ipt=1; ipt<=npt; ipt++) V2av += datar[idat2][ipt];
  V2av /= npt;

  fprintf(fl.out, "V1av=%lf V2av=%lf\n",  V1av, V2av);

  fprintf(fl.out, "cc npt=%d smap=%d\n", npt, crcpar->smap);
  for (ipt=1-npt; ipt<=npt-1; ipt++)
  {
    corar[ipt] /= npt - abs(ipt);
    corar[ipt] -= V1av * V2av; 
    if (!crcpar->smap)
    {
      fprintf(fl.cor, "%lf %lf\n", ipt * crcpar->deltat, corar[ipt]);
    }
  }
  fprintf(fl.cor, "  \n");

  free_dmatrix(var, 1, 2, 1, 2*npt);
  free_dvector(ans, 1, 4*npt);
}

/* ------------------------------------------------------------ */

/* This function computes the amplitudes and the phases of the  */
/* cross-correlation                                            */
/* Adapted from ampph.c                                         */
void amp_ph_cal (crc_par *crcpar, int idat1, int idat2, double *corar,
     fl_st fl)
{
  peak_str peakstr;

  crcpar->ntc = (crcpar->Tmax + crcpar->epsilon) / crcpar->deltat;
  fprintf(fl.out, "ntc=%d\n", crcpar->ntc);

  peakstr.mp = (int) (2.0 * crcpar->Tmax / crcpar->T_between_peaks);
  peakstr.iarmx = ivector(1, peakstr.mp);
  peakstr.np = 0;
  peakstr.iarmn = ivector(1, 4);
  peakstr.t_isect = dmatrix(1, 3, 1, 3);

  find_local_max(crcpar, &peakstr, corar, fl);
  find_peak_near_zero(crcpar, &peakstr, corar, fl);
  find_local_min(crcpar, &peakstr, corar, fl);
  find_intersections_with_zero(crcpar, &peakstr, corar, fl);

  free_ivector(peakstr.iarmx, 1, peakstr.mp);
  free_ivector(peakstr.iarmn, 1, 4);
  free_dmatrix(peakstr.t_isect, 1, 3, 1, 3);

  fprintf(fl.res, "%d %d %lf %lf %lf %lf\n", crcpar->icol[idat1],
  crcpar->icol[idat2], peakstr.amplitude, peakstr.t_period, peakstr.t_shift,
  peakstr.phase);
  fflush(fl.res);
}

/* This function calculates the local maxima, within time interval */
/* T_between_peaks                                                 */
void find_local_max(crc_par *crcpar, peak_str *peakstr, double *corar, 
     fl_st fl)
{
  int ii;
  int ipeak;

  for (ii=crcpar->kpoints - crcpar->ntc; ii<=crcpar->ntc - crcpar->kpoints; 
  ii++)
  {
    if (find_very_local_max(ii, crcpar, corar, fl))
    {
      if (!crcpar->smap)
        fprintf(fl.out, "%lf ii=%d\n", corar[ii], ii);
      if (peakstr->np == 0)
      {
        peakstr->np++;
        peakstr->iarmx[peakstr->np] = ii;
      }
      else 
      {
        if (ii * crcpar->deltat -  
            peakstr->iarmx[peakstr->np]* crcpar->deltat < 
            crcpar->T_between_peaks)
        {
          if (corar[ii] >  corar[peakstr->iarmx[peakstr->np]])
          { 
   	    peakstr->iarmx[peakstr->np] = ii;         
	  }
	}
        else
	{
          peakstr->np++;
          if (peakstr->np > peakstr->mp)
          {
            printf("np=%d > mp=%d\n", peakstr->np, peakstr->mp);
            exit(0);
	  }
          peakstr->iarmx[peakstr->np] = ii;
	}
      }
    }
  }

  if (!crcpar->smap)
  {
    fprintf(fl.out, "\n");
    for (ipeak=1; ipeak<=peakstr->np; ipeak++)
    {
      ii = peakstr->iarmx[ipeak];
      fprintf(fl.out, "%lf %lf %d %d\n", ii * crcpar->deltat, corar[ii], ii,
      ipeak);
    }
  }
}

/* This function calculates a very local maximum */
int find_very_local_max(int ii, crc_par *crcpar, double *corar, fl_st fl)
{
  int ll1, ll2, vlm, jj;

  ll1 = max(ii - crcpar->kpoints, -crcpar->ntc);
  ll2 = min(ii + crcpar->kpoints, crcpar->ntc);

  vlm = 1;

  for (jj=ll1; jj<=ii-1; jj++)
  {
    if (corar[jj] > corar[ii])
    {
      vlm=0;
      break;
    }
  }

  if (vlm)
  {
    for(jj=ii+1; jj<=ll2; jj++)
    {
      if (corar[jj] > corar[ii])
      {
        vlm=0;
        break;
      }
    }
  }

  return(vlm);
}

/* This function find the peak that is closest to tau=0.                */
void find_peak_near_zero(crc_par *crcpar, peak_str *peakstr, double *corar,
     fl_st fl)
{
  int ipeak, iold, inew;

  peakstr->zp = 1;

  for (ipeak=2; ipeak<=peakstr->np; ipeak++)
  {
    iold = peakstr->iarmx[peakstr->zp];
    inew = peakstr->iarmx[ipeak];
    if (fabs(inew) < fabs(iold))
    {
      peakstr->zp = ipeak;
    }
  }

  if (!crcpar->smap)
    fprintf(fl.out, "\nzp=%d ii=%d c=%lf %lf\n", peakstr->zp,
    peakstr->iarmx[peakstr->zp], peakstr->iarmx[peakstr->zp] * crcpar->deltat,
    corar[peakstr->iarmx[peakstr->zp]]);

  if (peakstr->zp < 3)
  {
    printf("zp=%d < 3\n", peakstr->zp);
    exit(0);
  }
}

/* This function calculates the four centeral local manima */
void find_local_min(crc_par *crcpar, peak_str *peakstr, double *corar,
     fl_st fl)
{
  int ipeak, jpeak, ileft, iright, ii;

  if (!crcpar->smap) fprintf(fl.out, "\n");
  for (ipeak=1; ipeak<=4; ipeak++)
  {
    jpeak = ipeak + peakstr->zp - 3;
    ileft = peakstr->iarmx[jpeak]+1;
    iright = peakstr->iarmx[jpeak+1]-1;
    
    if (!crcpar->smap)
      fprintf(fl.out, "ipeak=%d jpeak=%d ileft=%d iright=%d\n", ipeak, jpeak, 
      ileft, iright);

    peakstr->iarmn[ipeak] = ileft;
    for (ii=ileft+1; ii<=iright; ii++)
    {
      if (corar[ii] < corar[peakstr->iarmn[ipeak]])
      {
        peakstr->iarmn[ipeak] = ii;
      }
    }

    if (!crcpar->smap)
      fprintf(fl.out, "%lf %lf %d %d\n", peakstr->iarmn[ipeak]*crcpar->deltat,
      corar[peakstr->iarmn[ipeak]], ipeak, peakstr->iarmn[ipeak]);
  }

  if (!crcpar->smap)
  {
    fprintf(fl.out, "\nmax for amp\n");
    fprintf(fl.out, "%lf %lf %d\n", peakstr->iarmx[peakstr->zp-1] *
    crcpar->deltat, corar[peakstr->iarmx[peakstr->zp-1]], peakstr->zp-1);
    fprintf(fl.out, "%lf %lf %d\n", peakstr->iarmx[peakstr->zp+1] *
    crcpar->deltat, corar[peakstr->iarmx[peakstr->zp+1]], peakstr->zp-1);
  }

  peakstr->amplitude = 0.0;
  for (ipeak=-1; ipeak<=1; ipeak+=2)
  {
    peakstr->amplitude += 2.0 * corar[peakstr->iarmx[peakstr->zp+ipeak]];
  }
  for (ipeak=1; ipeak<=4; ipeak++)
  {
    peakstr->amplitude -= corar[peakstr->iarmn[ipeak]];
  }
  peakstr->amplitude /= 4.0;
  if (!crcpar->smap) fprintf(fl.out, "amplitude=%lf\n", peakstr->amplitude);
}

/* This function calculates intersections of the cross-correlation with */
/* zero                                                                 */
void find_intersections_with_zero(crc_par *crcpar, peak_str *peakstr, 
     double *corar, fl_st fl)
{
  int ipeak, pleft, pright;

  /* Finds the closest point to zero */
  if (!crcpar->smap) fprintf(fl.out, "\n");
  for (ipeak=1; ipeak<=3; ipeak++)
  {
    pleft = peakstr->iarmn[ipeak] + 1;
    pright = peakstr->iarmx[peakstr->zp + ipeak - 2] - 1;
    if (!crcpar->smap)
      fprintf(fl.out, "\npleft=%d nleft=%d\n", pleft, pright);
    find_one_intersection(pleft, pright, crcpar, corar,
    &peakstr->t_isect[ipeak][1],fl);
    if (!crcpar->smap) 
      fprintf(fl.out, "t_isect=%lf\n", peakstr->t_isect[ipeak][1]); 
  }

  for (ipeak=1; ipeak<=3; ipeak++)
  {
    pleft = peakstr->iarmx[peakstr->zp + ipeak - 2] + 1;
    pright = peakstr->iarmn[ipeak+1] - 1;
    if (!crcpar->smap)
      fprintf(fl.out, "\npleft=%d nleft=%d\n", pleft, pright);
    find_one_intersection(pleft, pright, crcpar, corar,
    &peakstr->t_isect[ipeak][2],fl);
    if (!crcpar->smap) 
      fprintf(fl.out, "t_isect=%lf\n", peakstr->t_isect[ipeak][2]);
  }

  if (!crcpar->smap) fprintf(fl.out, "\n");
  for (ipeak=1; ipeak<=3; ipeak++)
  {
    peakstr->t_isect[ipeak][3] = (peakstr->t_isect[ipeak][1] +
    peakstr->t_isect[ipeak][2]) / 2.0;
    if (!crcpar->smap)
      fprintf(fl.out, "%lf %lf %lf\n", peakstr->t_isect[ipeak][1], 
      peakstr->t_isect[ipeak][2], peakstr->t_isect[ipeak][3]);
  }

  peakstr->t_period = (peakstr->t_isect[3][3] - peakstr->t_isect[1][3]) / 2.0;
  peakstr->t_shift = (peakstr->t_isect[1][3] + peakstr->t_isect[2][3] +
    peakstr->t_isect[3][3]) / 3.0;
  peakstr->phase = peakstr->t_shift / peakstr->t_period;
  if (!crcpar->smap)
    fprintf(fl.out, "\nt_period=%lf t_shift=%lf phase=%lf\n",
    peakstr->t_period, peakstr->t_shift, peakstr->phase);
}

/* This function finds one intersection of the cross-correlation with zero. */
void find_one_intersection(int pleft, int pright, crc_par *crcpar, 
     double *corar, double *t_isect, fl_st fl)
{
  int ipl, ipr, find_isect_l, find_isect_r;

  /* Finds left intersection point */
  ipl = pleft;
  while (!(find_isect_l = cond_intersect(ipl, crcpar, corar, fl)) && 
         ipl < pright)
  {
   ipl++;
  }

  ipr = pright;
  while (!(find_isect_r = cond_intersect(ipr-1, crcpar, corar, fl)) &&
         ipr > pleft)
  {
   ipr--;
  }

  if (!crcpar->smap)
    fprintf(fl.out, "find_isect_l=%d ipl=%d find_isect_r=%d ipr=%d\n",
  find_isect_l, ipl, find_isect_r, ipr);

  if (find_isect_l && find_isect_r)
  {
    linear_regression(ipl, ipr, crcpar, corar, t_isect, fl);
  }

}

/* This function formulates the condition for intersection with zero  */
/* between two adjacent points.                                       */
int cond_intersect(int ip, crc_par *crcpar, double *corar, fl_st fl)
{
  int cond;

  cond = 0;

  if (fabs(corar[ip]) < crcpar->epsilon)
  {
    cond = 1;
  }
  else if (corar[ip] * corar[ip+1] < 0.0)
  {
    cond = 1;
  }

  return(cond);
}

void linear_regression(int ipl, int ipr, crc_par *crcpar, double *corar,
     double *t_isect, fl_st fl)
{
  double *xx, *yy;
  double sumx, sumxx, sumy, sumxy, denom, aa, bb;
  int lenv, ii;

  lenv = ipr - ipl + 2 * crcpar->inter_points - 1; 
  if (!crcpar->smap) fprintf(fl.out, "lenv=%d\n", lenv);

  xx = dvector(1, lenv);
  yy = dvector(1, lenv);

  for (ii=1; ii<=lenv; ii++)
  {
    xx[ii] = (ipl - crcpar->inter_points + ii) * crcpar->deltat;
    yy[ii] = corar[ipl - crcpar->inter_points + ii];
    if (!crcpar->smap)
      fprintf(fl.out, "ii=%d jj=%d xx=%lf yy=%lf\n", ii, 
      ipl - crcpar->inter_points + ii, xx[ii], yy[ii]);
  }

  sumx  = 0.0;
  sumxx = 0.0;
  sumy  = 0.0;
  sumxy = 0.0;

  for (ii=1; ii<=lenv; ii++)
  {
    sumx  += xx[ii];
    sumxx += xx[ii] * xx[ii];
    sumy  += yy[ii];
    sumxy += xx[ii] * yy[ii];
  }

  sumx  /= lenv;
  sumxx /= lenv;
  sumy  /= lenv;
  sumxy /= lenv;

  denom = sumxx - sumx * sumx;
  if (fabs(denom) < crcpar->epsilon)
  {
    *t_isect = -999.0;
  }
  else
  {
    aa = (-sumxy * sumx + sumxx * sumy) / denom;
    if (fabs(aa) < crcpar->epsilon)
    {
      *t_isect = -998.0;
    }
    else
    {
      bb = (sumxy - sumx * sumy) / denom;
      *t_isect = (crcpar->base_intersect - aa) / bb;
      if (!crcpar->smap)
        fprintf(fl.out, "denom=%lf aa=%lf bb=%lf ti=%lf\n", denom, aa, bb,
        *t_isect);
    }
  }

  free_dvector(xx, 1, lenv);
  free_dvector(yy, 1, lenv);
}

/* CORREL */
void xcorrel(double data1[], double data2[], unsigned long n, double ans[])
{
        unsigned long no2,i;
        double dum,*fft;

        fft=dvector(1,n<<2);
        xtwofft(data1,data2,fft,ans,n);

        no2=n>>1;
        for (i=2;i<=n+2;i+=2) {
            ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/no2;
            ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/no2;
	}
        ans[2]=ans[n+1];
        realft(ans,n,-1);
        free_dvector(fft,1,n<<2);
}

void xtwofft(double data1[], double data2[], double fft1[], double fft2[],
	unsigned long n)
{
	unsigned long nn3,nn2,jj,j;
	double rep,rem,aip,aim;

	nn3=1+(nn2=2+n+n);

	for (j=1,jj=2;j<=n;j++,jj+=2) {
		fft1[jj-1]=data1[j];
		fft1[jj]=data2[j];
	}

	xfour1(fft1,n,1);

	fft2[1]=fft1[2];
	fft1[2]=fft2[2]=0.0;

	for (j=3;j<=n+1;j+=2) {
		rep=0.5*(fft1[j]+fft1[nn2-j]);
		rem=0.5*(fft1[j]-fft1[nn2-j]);
		aip=0.5*(fft1[j+1]+fft1[nn3-j]);
		aim=0.5*(fft1[j+1]-fft1[nn3-j]);
		fft1[j]=rep;
		fft1[j+1]=aim;
		fft1[nn2-j]=rep;
		fft1[nn3-j] = -aim;
		fft2[j]=aip;
		fft2[j+1] = -rem;
		fft2[nn2-j]=aip;
		fft2[nn3-j]=rem;
	}
}

void xfour1(double data[], unsigned long nn, int isign)
{
	unsigned long n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	double tempr,tempi,temp;

	n=nn << 1;
	j=1;
	for (i=1;i<n;i+=2) {
		if (j > i) {
			SWAP(data[j],data[i]);
			SWAP(data[j+1],data[i+1]);
		}
		m=n >> 1;
		while (m >= 2 && j > m) {
			j -= m;
			m >>= 1;
		}
		j += m;
	}
	mmax=2;

	while (n > mmax) {
		istep=mmax << 1;
		theta=isign*(6.28318530717959/mmax);
		wtemp=sin(0.5*theta);
		wpr = -2.0*wtemp*wtemp;
		wpi=sin(theta);
		wr=1.0;
		wi=0.0;

		for (m=1;m<mmax;m+=2) {

			for (i=m;i<=n;i+=istep) {
				j=i+mmax;
				tempr=wr*data[j]-wi*data[j+1];
				tempi=wr*data[j+1]+wi*data[j];
				data[j]=data[i]-tempr;
				data[j+1]=data[i+1]-tempi;
				data[i] += tempr;
				data[i+1] += tempi;
			}

			wr=(wtemp=wr)*wpr-wi*wpi+wr;
			wi=wi*wpr+wtemp*wpi+wi;
		}

		mmax=istep;
	}

}