/* This file contains two files (cplex.c and nrutil.c).  Use an
editor and break them up. */

//* ***********************************************************************  *//
//* ****************    cplex.c file  * ***********************************  *//
//* ***********************************************************************  *//

class complex {
 double re, im;
 public:
  void operator+=(double a);
  void operator+=(complex x);
  void operator-=(double a);
  void operator-=(complex x);
  void operator*=(double a);
  void operator*=(complex x);
  void operator/=(double a);
  void operator/=(complex x);

  friend complex operator+(double a, complex x);
  friend complex operator+(complex x,double a);
  friend complex operator+(complex a, complex x);
  friend complex operator-(double a, complex x);
  friend complex operator-(complex x,double a);
  friend complex operator-(complex a, complex x);
  friend complex operator*(double a, complex x);
  friend complex operator*(complex x,double a);
  friend complex operator*(complex a, complex x);
  friend complex operator/(complex x, double a);
  friend complex operator/(complex a, complex x);
  friend complex operator/(double a, complex x);
  friend complex operator^(complex x, double a);
  friend double  operator&&(complex a, complex x);


  friend complex cconj(complex x);
  friend complex cx(double a, double b);
  friend complex cexp(complex x);
  friend complex cexp(double a);
  friend complex clog(complex x);
  friend double cmag(complex x);
  friend void ctor(complex x, double *a, double *b);
//*  Note: when calling this, ctor(x,&a,&b) for constants a + ib. *//
  friend double real(complex x);
  friend double imag(complex x);
  friend complex I(complex x);

  friend void cprintf(complex c);
  friend void cfprintf(FILE *output, complex c);

  friend void ciprintf(complex c);
  friend void cifprintf(FILE *output, complex c);

  friend void nrerror(char error_text[]);
  friend complex *cvector(int nl,int nh);
  friend complex **cmatrix(int nrl,int nrh,int ncl,int nch);
  friend void free_cmatrix(complex **m,int nrl,int nrh,int ncl,int nch);
  friend void free_cvector(complex *v,int nl,int nh);

};

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

void complex::operator+=(double a) {
  re += a;
}

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

void complex::operator+=(complex x) {
  re += x.re;
  im += x.im;
}

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

void complex::operator-=(double a) {
  re -= a;
}

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

void complex::operator-=(complex x) {
 re -= x.re;
 im -= x.im;
}

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

void complex::operator*=(double a) {
 re *= a;
 im *= a;
}

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

void complex::operator*=(complex x) {
 re = re*x.re - im*x.im;
 im = im*x.re + re*x.im;
}

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

void complex::operator/=(double a) {
 re /= a;
 im /= a;
}

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

void complex::operator/=(complex x) {
 re = (re*x.re + im*x.im)/(x.re*x.re + x.im*x.im);
 im = (im*x.re - re*x.im)/(x.re*x.re + x.im*x.im);
}


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

complex operator+(double a, complex x) {
complex c;

c.re = a + x.re;
c.im = x.im;
return(c);
}

complex operator+(complex x,double a) {
complex c;

c.re = a + x.re;
c.im = x.im;
return(c);
}

complex operator+(complex a, complex x) {
complex c;

c.re = a.re + x.re;
c.im = a.im + x.im;
return(c);
}

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

complex operator-(double a, complex x) {
complex c;

c.re = a - x.re;
c.im = -1.*x.im;
return(c);
}

complex operator-(complex x,double a) {
complex c;

c.re = x.re - a;
c.im = x.im;
return(c);
}

complex operator-(complex a, complex x) {
complex c;

c.re = a.re - x.re;
c.im = a.im - x.im;
return(c);
}

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

complex operator*(double a, complex x) {
complex c;

c.re = a * x.re;
c.im = a * x.im;
return(c);
}

complex operator*(complex x,double a) {
complex c;

c.re = a * x.re;
c.im = a * x.im;
return(c);
}

complex operator*(complex a, complex x) {
complex c;

c.re = a.re * x.re - a.im * x.im;
c.im = a.re * x.im + a.im * x.re;
return(c);
}

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

complex operator/(complex x, double a) {
complex c;

c.re = x.re/a;
c.im = x.im/a;
return(c);
}

complex operator/(complex a, complex x) {
complex c;
double magx;

magx = x.re * x.re + x.im * x.im;
c.re = (a.re * x.re + a.im * x.im)/magx;
c.im = (a.im * x.re - a.re * x.im)/magx;
return(c);
}

complex operator/(double a, complex x) {
complex c;
double magx;

magx = x.re * x.re + x.im * x.im;
c.re = a*x.re/magx;
c.im = -a*x.im/magx;
return(c);
}

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

complex operator^(complex x, double a) {
complex c;
double magx;

magx = x.re * x.re + x.im * x.im;
c.re = pow(magx,a/2.0)*cos(a*atan2(x.im,x.re));
c.im = pow(magx,a/2.0)*sin(a*atan2(x.im,x.re));

return(c);
}

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

double operator&&(complex a, complex x)  {
  
return (a.re*x.re + a.im*x.im);
}

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

complex cconj(complex x)  {
complex c;

   c.re = x.re;
   c.im = -x.im;

return (c);
}

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

complex cexp(complex x)  {
complex c;

c.re = exp(x.re)*cos(x.im);
c.im = exp(x.re)*sin(x.im);

return(c);
}

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

complex cexp(double a)  {
complex c;
//* Note: this is when it's exp(i*a), a real     *//

c.re = cos(a);
c.im = sin(a);

return(c);
}

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

complex clog(complex x)  {
complex c;

c.re = log(sqrt(x.re*x.re + x.im*x.im));
c.im = atan2(x.im,x.re);

return(c);
}

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

double cmag(complex x)  {
    return  (sqrt(x.re*x.re + x.im*x.im));
}

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

complex cx(double a, double b) {
complex c;

c.re = a;
c.im = b;

return (c);
}


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

void ctor(complex x, double *a, double *b)  {

*a = x.re;
*b = x.im;

}


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

double real(complex x)   {
   return (x.re);
}

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

double imag(complex x)   {
   return (x.im);
}

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

complex I(complex x)   {
complex c;
  
c.re = -1.*x.im;
c.im = x.re;

return(c);;
}


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

void cprintf(complex c)   {

if (c.im < 0.0)
   printf ("%-10lf   -%-8lf", c.re, fabs(c.im));
else
   printf ("%-10lf   +%-8lf", c.re, c.im);
   
}
 

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

void cfprintf(FILE *output, complex c)   {

if (c.im < 0.0)
   fprintf (output, "%-10lf   -%-8lf", c.re, fabs(c.im));
else
   fprintf (output, "%-10lf   +%-8lf", c.re, c.im);
   
}


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

void ciprintf(complex c)   {

if (c.im < 0.0)
   printf ("%-10lf   -%-8lfi", c.re, fabs(c.im));
else
   printf ("%-10lf   +%-8lfi", c.re, c.im);

}


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

void cifprintf(FILE *output, complex c)   {

if (c.im < 0.0)
   fprintf (output, "%-10lf   -%-8lfi", c.re, fabs(c.im));
else
   fprintf (output, "%-10lf   +%-8lfi", c.re, c.im);

}


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


void nrerror(char error_text[])
{
	void exit();

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}


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

complex *cvector(int nl,int nh)
{
        complex *v;

        v=(complex *)malloc((unsigned) (nh-nl+1)*sizeof(complex));
        if (!v) nrerror("allocation failure in vector()");
        return v-nl;
}

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

complex **cmatrix(int nrl,int nrh,int ncl,int nch)
{
        int i;
        complex **m;

        m=(complex **) malloc((unsigned) (nrh-nrl+1)*sizeof(complex*));
        if (!m) nrerror("allocation failure 1 in matrix()");
        m -= nrl;

        for(i=nrl;i<=nrh;i++) {
                m[i]=(complex *) malloc((unsigned) (nch-ncl+1)*sizeof(complex));
                if (!m[i]) nrerror("allocation failure 2 in matrix()");
                m[i] -= ncl;
        }
        return m;
}

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

void free_cvector(complex *v,int nl,int nh)
{
        free((char*) (v+nl));
}

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

void free_cmatrix(complex **m,int nrl,int nrh,int ncl,int nch)
{
        int i;

        for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
        free((char*) (m+nrl));
}


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

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

//* ********************************************************************   *//
//* ************    nrutil.c file   ************************************   *//
//* ********************************************************************   *//

double step(double a)   //* Not in nrutil.c *//
{
 if (a>=0.0) return(1.0);
 else return(0.0);
}


float *vector(int nl,int nh)
{
	float *v;

	v=(float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl;
}

int *ivector(int nl,int nh)
{
	int *v;

	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl;
}

double *dvector(int nl,int nh)
{
	double *v;

	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl;
}



float **matrix(int nrl,int nrh,int ncl,int nch)
{
	int i;
	float **m;

	m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
		if (!m[i]) nrerror("allocation failure 2 in matrix()");
		m[i] -= ncl;
	}
	return m;
}

double **dmatrix(int nrl,int nrh,int ncl,int nch)
{
	int i;
	double **m;

	m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
	if (!m) nrerror("allocation failure 1 in dmatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
		if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
		m[i] -= ncl;
	}
	return m;
}

int **imatrix(int nrl,int nrh,int ncl,int nch)
{
	int i,**m;

	m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
	if (!m) nrerror("allocation failure 1 in imatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
		if (!m[i]) nrerror("allocation failure 2 in imatrix()");
		m[i] -= ncl;
	}
	return m;
}



float **submatrix(float **a,int oldrl,int oldrh,int oldcl,int oldch,
                  int newrl,int newcl)
{
	int i,j;
	float **m;

	m=(float **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(float*));
	if (!m) nrerror("allocation failure in submatrix()");
	m -= newrl;

	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;

	return m;
}



void free_vector(float *v,int nl,int nh)
{
	free((char*) (v+nl));
}

void free_ivector(int *v,int nl,int nh)
{
	free((char*) (v+nl));
}

void free_dvector(double *v,int nl,int nh)
{
	free((char*) (v+nl));
}



void free_matrix(float **m,int nrl,int nrh,int ncl,int nch)
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_dmatrix(double **m,int nrl,int nrh,int ncl,int nch)
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_imatrix(int **m,int nrl,int nrh,int ncl,int nch)
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}



void free_submatrix(float **b,int nrl,int nrh,int ncl,int nch)
{
	free((char*) (b+nrl));
}



float **convert_matrix(float *a,int nrl,int nrh,int ncl,int nch)
{
	int i,j,nrow,ncol;
	float **m;

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	m = (float **) malloc((unsigned) (nrow)*sizeof(float*));
	if (!m) nrerror("allocation failure in convert_matrix()");
	m -= nrl;
	for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
	return m;
}



void free_convert_matrix(float **b,int nrl,int nrh,int ncl,int nch)
{
	free((char*) (b+nrl));
}