/* Steven Andrews, 5/99 */
/* See documentation called Cn doc */
/* Copyright 2003 by Steven Andrews.  Permission is granted
   for non-commercial use of and modifications to the code. */

#include "Rn.h"
#include "Cn.h"
#include "RnSort.h"
#include "math2.h"

float *makecmplx(float *ar,float *ai,float *c,int n) {
	int i;
	
	if(ar&&ai)
		for(i=0;i<n;i++)	{c[2*i]=ar[i];c[2*i+1]=ai[i];}
	else if(ar)
		for(i=0;i<n;i++)	{c[2*i]=ar[i];c[2*i+1]=0;}
	else if(ai)
		for(i=0;i<n;i++)	{c[2*i]=0;c[2*i+1]=ai[i];}
	else
		for(i=0;i<n;i++)	{c[2*i]=c[2*i+1]=0;}
	return c; }

float *real(float *a,float *cr,int n) {
	int i;
	
	for(i=0;i<n;i++)	cr[i]=a[2*i];
	return cr; }

float *imag(float *a,float *ci,int n) {
	int i;

	for(i=0;i<n;i++)	ci[i]=a[2*i+1];
	return ci; }

float *magnitude(float *a,float *cr,int n) {
	int i;
	
	for(i=0;i<n;i++) cr[i]=sqrt(a[2*i]*a[2*i]+a[2*i+1]+a[2*i+1]);
	return cr; }

float *cmplxphase(float *a,float *cr,int n) {
	int i;
	
	for(i=0;i<n;i++) cr[i]=a[2*i]?atan(a[2*i+1]/a[2*i]):0;
	return cr; }

float *CompConj(float *a,float *c,int n) {
	int i;
	
	for(i=0;i<2*n;i+=2)	c[i]=a[i];
	for(i=1;i<2*n;i+=2)	c[i]=-a[i];
	return c; }

float *rotateCV(float *a,float *c,int n,int p) {
	int i;
	float temp;
	
	if(p<0) p+=4*(-p/4+1);
	else p=p%4;
	if(p==1)
		for(i=0;i<2*n;i+=2)	{temp=-a[i+1];c[i+1]=a[i];c[i]=temp;}
	else if(p==2)
		for(i=0;i<2*n;i+=2)	{c[i]=-a[i];c[i+1]=-a[i+1];}
	else if(p==3)
		for(i=0;i<2*n;i+=2)	{temp=a[i+1];c[i+1]=-a[i];c[i]=temp;}
	else
		for(i=0;i<2*n;i+=2)	{c[i]=a[i];c[i+1]=a[i+1];}
	return c; }

float *rotate2CV(float *a,float *c,int n,float phi) {
	int i;
	float cp,sp,temp;
	
	cp=cos(phi);
	sp=sin(phi);
	for(i=0;i<2*n;i+=2)	{
		temp=cp*a[i]-sp*a[i+1];
		c[i+1]=sp*a[i]+cp*a[i+1];
		c[i]=temp;	}
	return c; }

float *multeikx(float *xr,float *a,float *c,int n,float k) {
	int i;
	float th,temp;
	
	for(i=0;i<n;i++)	{
		th=k*xr[i];
		temp=a[2*i]*cos(th)-a[2*i+1]*sin(th);
		c[2*i+1]=a[2*i]*sin(th)+a[2*i+1]*cos(th);
		c[2*i]=temp; }
	return c; }

float *multCV(float *a,float *b,float *c,int n) {
	int i;
	float temp;
	
	for(i=0;i<n;i++)	{
		temp=a[2*i]*b[2*i]-a[2*i+1]*b[2*i+1];
		c[2*i+1]=a[2*i]*b[2*i+1]+a[2*i+1]*b[2*i];
		c[2*i]=temp; }
	return c; }

float *deriv2CV(float *a,float *c,int n,int p) {
	int i;
	
	for(i=1;i<2*(n-1);i++)
		c[i]=a[i-2]+a[i+2]-2*a[i];
	if(p)	{														// set p=1 for periodic boundary conditions
		c[0]=a[2*n-2]+a[2]-2*a[0];
		c[1]=a[2*n-1]+a[3]-2*a[1];
		c[2*n-2]=a[2*n-4]+a[0]-2*a[2*n-2];
		c[2*n-1]=a[2*n-3]+a[1]-2*a[2*n-1]; }
	else {
		c[0]=a[0]+a[4]-2*a[2];
		c[1]=a[1]+a[5]-2*a[3];
		c[2*n-2]=a[2*n-6]+a[2*n-2]-2*a[2*n-4];
		c[2*n-1]=a[2*n-5]+a[2*n-1]-2*a[2*n-3]; }
	return c; }

float *integCV(float *a,float *c,int n) {
	int i;
	
	c[0]=a[0]/2;
	c[1]=a[1]/2;
	for(i=2;i<2*n;i++)
		c[i]=c[i-2]+(a[i-2]+a[i])/2;
	return c; }

float FTStartDflt(float *xr,int n) {
	if(n%2) return -3.14159265358*(n-1.0)/n*(n-1.0)/(xr[n-1]-xr[0]);
	else return -3.14159265358*(n-1.0)/(xr[n-1]-xr[0]); }

float *fourier(float *xr,float *a,float *kr,float *c,int nx,int nk,int isign) {
	int i,j;
	double sumr,sumi,dx,cp,sp;
	
	dx=(xr[nx-1]-xr[0])/(nx-1)/sqrt(2*PI);
	for(i=0;i<nk;i++)	{
		sumr=sumi=0;
		for(j=0;j<nx;j++)	{		// sum+=exp(isign*kr*xr)*a[j]
			cp=cos(kr[i]*xr[j]);
			sp=isign*sin(kr[i]*xr[j]);
			sumr+=cp*a[2*j]-sp*a[2*j+1];
			sumi+=cp*a[2*j+1]+sp*a[2*j]; }
		c[2*i]=sumr*dx;
		c[2*i+1]=sumi*dx; }
	return c; }

float *realcosft(float *xr,float *ar,float *kr,float *cr,int nx,int nk) {
	int i,j;
	double sumr,dx;

	dx=(xr[nx-1]-xr[0])/(nx-1)/sqrt(PI/2);
	for(i=0;i<nk;i++)	{
		sumr=0;
		for(j=0;j<nx;j++)
			sumr+=cos(kr[i]*xr[j])*ar[j];
		cr[i]=sumr*dx; }
	return cr; }

float *hankel(float *xr,float *ar,float *kr,float *cr,int nx,int nk,int samp) {
	int i,j;
	double sum,dx,x;

	dx=(xr[nx-1]-xr[0])/(nx-1);
	dx/=samp;
	for(i=0;i<nk;i++) {
		sum=0;
		j=-2;
		for(x=xr[0];x<=xr[nx-1];x+=dx)
			sum+=interpolate1(xr,ar,nx,&j,x)*bessj0(kr[i]*x)*x;
		cr[i]=sum*dx; }
	return cr; }

float *fft(float *xr,float *a,float *kr,float *c,int nn,int isign) {
	int ofst;
	double dx,dk;
	
	float *data;
	int n,mmax,m,j,istep,i;
	double wtemp,wr,wpr,wpi,wi,theta;
	float tempr,tempi;
	
	copyV(a,c,2*nn);
	data=c-1;
	// start of Numerical Recipes procedure	
	n=nn<<1;
	j=1;
	for(i=1;i<n;i+=2)	{
		if(j>i)	{
			tempr=data[j];data[j]=data[i];data[i]=tempr;
			tempi=data[j+1];data[j+1]=data[i+1];data[i+1]=tempi;}
		m=n>>1;
		while(m>=2&&j>m)	{
			j-=m;
			m>>=1; }
		j+=m; }
	mmax=2;
	while(n>mmax)	{
		istep=2*mmax;
		theta=6.28318530717959/(isign*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; }
	// end of Numerical Recipes procedure
	dx=(xr[nn-1]-xr[0])/(nn-1);
	dk=2*PI/(nn*dx);
	ofst=floor(kr[0]/dk+0.5);
	for(j=0;j<nn;j++) kr[j]=(j+ofst)*dk;
	leftrotV(c,c,2*nn,2*ofst);
	multKV(dx/sqrt(2*PI),c,c,2*nn);
	multeikx(kr,c,c,nn,-xr[0]);
	return c; }