/**********
Copyright 1990 Regents of the University of California.  All rights reserved.
Author: 1985 Wayne A. Christopher, U. C. Berkeley CAD Group
**********/

/*
 * Routines to do complex mathematical functions. These routines require
 * the -lm libraries. We sacrifice a lot of space to be able
 * to avoid having to do a seperate call for every vector element,
 * but it pays off in time savings.  These routines should never
 * allow FPE's to happen.
 *
 * Complex functions are called as follows:
 *  cx_something(data, type, length, &newlength, &newtype),
 *  and return a char * that is cast to complex or double.
 *
 */

#include <errno.h>

#include "ngspice/ngspice.h"
#include "ngspice/memory.h"
#include "ngspice/cpdefs.h"
#include "ngspice/dvec.h"

#include "cmath.h"
#include "cmath1.h"

#ifdef HAS_WINGUI
#define win_x_fprintf fprintf
#endif

/* This flag determines whether degrees or radians are used. The radtodeg
 * and degtorad macros are no-ops if this is FALSE.
 */

bool cx_degrees = FALSE;

void *
cx_mag(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    double *dd = (double *) data;
    ngcomplex_t *cc = (ngcomplex_t *) data;
    int i;

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_REAL)
        for (i = 0; i < length; i++)
            d[i] = fabs(dd[i]);
    else
        for (i = 0; i < length; i++)
            d[i] = cmag(cc[i]);
    return ((void *) d);
}

void *
cx_ph(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    ngcomplex_t *cc = (ngcomplex_t *) data;
    int i;

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX)
        for (i = 0; i < length; i++) {
            d[i] = radtodeg(cph(cc[i]));
        }
    /* Otherwise it is 0, but tmalloc zeros the stuff already. */
    return ((void *) d);
}

/* SJdV Modified from above to find closest from +2pi,0, -2pi */
void *
cx_cph(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    ngcomplex_t *cc = (ngcomplex_t *) data;
    int i;

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX) {
        double last_ph = cph(cc[0]);
        d[0] = radtodeg(last_ph);
        for (i = 1; i < length; i++) {
            double ph = cph(cc[i]);
            last_ph = ph - (2*M_PI) * floor((ph - last_ph)/(2*M_PI) + 0.5);
            d[i] = radtodeg(last_ph);
        }
    }
    /* Otherwise it is 0, but tmalloc zeros the stuff already. */
    return ((void *) d);
}

/* Modified from above but with real phase vector in degrees as input */
void *
cx_unwrap(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    double *dd = (double *) data;
    int i;

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_REAL) {
        double last_ph = degtorad(dd[0]);
        d[0] = last_ph;
        for (i = 1; i < length; i++) {
            double ph = degtorad(dd[i]);
            last_ph = ph - (2*M_PI) * floor((ph - last_ph)/(2*M_PI) + 0.5);
            d[i] = radtodeg(last_ph);
        }
    }
    /* Otherwise it is 0, but tmalloc zeros the stuff already. */
    return ((void *) d);
}

/* If this is pure imaginary we might get real, but never mind... */
void *cx_j(void *data, short int type, int length, int *newlength,
        short int *newtype)
{
    ngcomplex_t *c = alloc_c(length);
    *newlength = length;
    *newtype = VF_COMPLEX;

    if (type == VF_COMPLEX) {
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = -imagpart(cc[i]);
            imagpart(c[i]) = realpart(cc[i]);
        }
    }
    else {
        double *dd = (double *) data;
        int i;
        for (i = 0; i < length; i++) {
            imagpart(c[i]) = dd[i];
            /* Real part is already 0. */
        }
    }
    return (void *) c;
}

void *cx_real(void *data, short int type, int length, int *newlength,
        short int *newtype)
{
    double *d = alloc_d(length);

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX) {
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        for (i = 0; i < length; i++) {
            d[i] = realpart(cc[i]);
        }
    }
    else {
        double *dd = (double *) data;
        int i;
        for (i = 0; i < length; i++) {
            d[i] = dd[i];
        }
    }
    return (void *) d;
}

void *cx_imag(void *data, short int type, int length, int *newlength,
        short int *newtype)
{
    double *d = alloc_d(length);

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX) {
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        for (i = 0; i < length; i++) {
            d[i] = imagpart(cc[i]);
        }
    }
    else {
        double *dd = (double *) data;
        int i;
        for (i = 0; i < length; i++) {
            d[i] = dd[i];
        }
    }
    return (void *) d;
}



/* Create complex conjugate of data */
void *cx_conj(void *data, short int type, int length,
        int *p_newlength, short int *p_newtype)
{
    /* Length and type do not change */
    *p_newlength = length;
    *p_newtype = type;

    /* For complex, copy with conjugation */
    if (type == VF_COMPLEX) {
        ngcomplex_t * const c_dst = alloc_c(length);
        ngcomplex_t *c_dst_cur = c_dst;
        ngcomplex_t *c_src_cur = (ngcomplex_t *) data;
        ngcomplex_t * const c_src_end = c_src_cur + length;
        for ( ; c_src_cur < c_src_end;  c_src_cur++, c_dst_cur++) {
            c_dst_cur->cx_real = c_src_cur->cx_real;
            c_dst_cur->cx_imag = -c_src_cur->cx_imag;
        }
        return (void *) c_dst;
    }

    /* Else real, so just copy */
    return memcpy(alloc_d(length), data, (unsigned int) length * sizeof(double));
} /* end of function cx_conj */



/* This is obsolete... */

void *
cx_pos(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    double *dd = (double *) data;
    ngcomplex_t *cc = (ngcomplex_t *) data;
    int i;

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX)
        for (i = 0; i < length; i++)
            d[i] = ((realpart(cc[i]) > 0.0) ? 1.0 : 0.0);
    else
        for (i = 0; i < length; i++)
            d[i] = ((dd[i] > 0.0) ? 1.0 : 0.0);
    return ((void *) d);
}

void *cx_db(void *data, short int type, int length,
        int *newlength, short int *newtype)
{
    int xrc = 0;
    double *d = alloc_d(length);

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_COMPLEX) {
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        for (i = 0; i < length; i++) {
            const double tt = cmag(cc[i]);
            rcheck(tt > 0, "db");
            /*
                if (tt == 0.0)
                    d[i] = 20.0 * - log(HUGE);
                else
            */
            d[i] = 20.0 * log10(tt);
        }
    }
    else {
        double *dd = (double *) data;
        int i;
        for (i = 0; i < length; i++) {
            const double tt = dd[i];
            rcheck(tt > 0, "db");
            /*
                if (dd[i] == 0.0)
                    d[i] = 20.0 * - log(HUGE);
                else
            */
            d[i] = 20.0 * log10(tt);
        }
    }

EXITPOINT:
    if (xrc != 0) {
        txfree(d);
        d = (double *) NULL;
    }
    return ((void *) d);
} /* end of function cx_db */



void *cx_log10(void *data, short int type, int length,
        int *newlength, short int *newtype)
{
    int xrc = 0;
    void *rv;

    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        rv = c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            double td;

            td = cmag(cc[i]);
            /* Perhaps we should trap when td = 0.0, but Ken wants
             * this to be possible...
             */
            rcheck(td >= 0, "log10");
            if (td == 0.0) {
                realpart(c[i]) = - log10(HUGE);
                imagpart(c[i]) = 0.0;
            }
            else {
                realpart(c[i]) = log10(td);
                imagpart(c[i]) = atan2(imagpart(cc[i]), realpart(cc[i]));
            }
        }
    }
    else {
        double *d;
        double *dd = (double *) data;
        int i;

        rv = d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++) {
            rcheck(dd[i] >= 0, "log10");
            if (dd[i] == 0.0) {
                d[i] = - log10(HUGE);
            }
            else {
                d[i] = log10(dd[i]);
            }
        }
    }

    *newlength = length;

EXITPOINT:
    if (xrc != 0) { /* Free resources on error */
        txfree(rv);
        rv = NULL;
    }

    return rv;
} /* end of function cx_log10 */



void *cx_log(void *data, short int type, int length,
        int *newlength, short int *newtype)
{
    int xrc = 0;
    void *rv;

    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;

        rv = c = alloc_c(length);
        *newtype = VF_COMPLEX;

        int i;
        for (i = 0; i < length; i++) {
            double td;

            td = cmag(cc[i]);
            rcheck(td >= 0, "log");
            if (td == 0.0) {
                realpart(c[i]) = - log(HUGE);
                imagpart(c[i]) = 0.0;
            }
            else {
                realpart(c[i]) = log(td);
                imagpart(c[i]) = atan2(imagpart(cc[i]), realpart(cc[i]));
            }
        }
    }
    else {
        double *d;
        double *dd = (double *) data;

        rv = d = alloc_d(length);
        *newtype = VF_REAL;

        int i;
        for (i = 0; i < length; i++) {
            rcheck(dd[i] >= 0, "log");
            if (dd[i] == 0.0)
                d[i] = - log(HUGE);
            else
                d[i] = log(dd[i]);
        }
    }

    *newlength = length;

EXITPOINT:
    if (xrc != 0) { /* Free resources on error */
        txfree(rv);
        rv = NULL;
    }

    return rv;
} /* end of function cx_log */



void *
cx_exp(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            double td;

            td = exp(realpart(cc[i]));
            realpart(c[i]) = td * cos(imagpart(cc[i]));
            imagpart(c[i]) = td * sin(imagpart(cc[i]));
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = exp(dd[i]);
        return ((void *) d);
    }
}

void *
cx_sqrt(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = NULL;
    ngcomplex_t *c = NULL;
    double *dd = (double *) data;
    ngcomplex_t *cc = (ngcomplex_t *) data;
    int i, cres = (type == VF_REAL) ? 0 : 1;

    if (type == VF_REAL)
        for (i = 0; i < length; i++)
            if (dd[i] < 0.0)
                cres = 1;
    if (cres) {
        c = alloc_c(length);
        *newtype = VF_COMPLEX;
    } else {
        d = alloc_d(length);
        *newtype = VF_REAL;
    }
    *newlength = length;
    if (type == VF_COMPLEX) {
        for (i = 0; i < length; i++) {
            if (realpart(cc[i]) == 0.0) {
                if (imagpart(cc[i]) == 0.0) {
                    realpart(c[i]) = 0.0;
                    imagpart(c[i]) = 0.0;
                } else if (imagpart(cc[i]) > 0.0) {
                    realpart(c[i]) = sqrt (0.5 *
                                            imagpart(cc[i]));
                    imagpart(c[i]) = realpart(c[i]);
                } else {
                    imagpart(c[i]) = sqrt( -0.5 *
                                            imagpart(cc[i]));
                    realpart(c[i]) = - imagpart(c[i]);
                }
            } else if (realpart(cc[i]) > 0.0) {
                if (imagpart(cc[i]) == 0.0) {
                    realpart(c[i]) =
                        sqrt(realpart(cc[i]));
                    imagpart(c[i]) = 0.0;
                } else if (imagpart(cc[i]) < 0.0) {
                    realpart(c[i]) = -sqrt(0.5 *
                                            (cmag(cc[i]) + realpart(cc[i])));
                } else {
                    realpart(c[i]) = sqrt(0.5 *
                                           (cmag(cc[i]) + realpart(cc[i])));
                }
                imagpart(c[i]) = imagpart(cc[i]) / (2.0 *
                                                      realpart(c[i]));
            } else { /* realpart(cc[i]) < 0.0) */
                if (imagpart(cc[i]) == 0.0) {
                    realpart(c[i]) = 0.0;
                    imagpart(c[i]) =
                        sqrt(- realpart(cc[i]));
                } else {
                    if (imagpart(cc[i]) < 0.0)
                        imagpart(c[i]) = - sqrt(0.5 *
                                                 (cmag(cc[i]) -
                                                  realpart(cc[i])));
                    else
                        imagpart(c[i]) = sqrt(0.5 *
                                               (cmag(cc[i]) -
                                                realpart(cc[i])));
                    realpart(c[i]) = imagpart(cc[i]) /
                                      (2.0 * imagpart(c[i]));
                }
            }
        }
        return ((void *) c);
    } else if (cres) {
        for (i = 0; i < length; i++)
            if (dd[i] < 0.0)
                imagpart(c[i]) = sqrt(- dd[i]);
            else
                realpart(c[i]) = sqrt(dd[i]);
        return ((void *) c);
    } else {
        for (i = 0; i < length; i++)
            d[i] = sqrt(dd[i]);
        return ((void *) d);
    }
}

void *
cx_sin(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = sin(degtorad(realpart(cc[i]))) *
                              cosh(degtorad(imagpart(cc[i])));
            imagpart(c[i]) = cos(degtorad(realpart(cc[i]))) *
                              sinh(degtorad(imagpart(cc[i])));
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = sin(degtorad(dd[i]));
        return ((void *) d);
    }
}

void *
cx_sinh(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        double u, v;
        c = alloc_c(length);

        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            /* sinh(x+iy) = sinh(x)*cos(y) + i * cosh(x)*sin(y) */
            u = degtorad(realpart(cc[i]));
            v = degtorad(imagpart(cc[i]));
            realpart(c[i]) = sinh(u)*cos(v);
            imagpart(c[i]) = cosh(u)*sin(v);
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = sinh(degtorad(dd[i]));
        return ((void *) d);
    }
}


void *
cx_cos(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        c = alloc_c(length);
        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            realpart(c[i]) = cos(degtorad(realpart(cc[i]))) *
                              cosh(degtorad(imagpart(cc[i])));
            imagpart(c[i]) = - sin(degtorad(realpart(cc[i]))) *
                              sinh(degtorad(imagpart(cc[i])));
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = cos(degtorad(dd[i]));
        return ((void *) d);
    }
}


void *
cx_cosh(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *c;
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;
        double u, v;

        c = alloc_c(length);

        *newtype = VF_COMPLEX;
        for (i = 0; i < length; i++) {
            /* cosh(x+iy) = cosh(x)*cos(y) + i * sinh(x)*sin(y) */
            u = degtorad(realpart(cc[i]));
            v = degtorad(imagpart(cc[i]));
            realpart(c[i]) = cosh(u)*cos(v);
            imagpart(c[i]) = sinh(u)*sin(v);
        }
        return ((void *) c);
    } else {
        double *d;
        double *dd = (double *) data;
        int i;

        d = alloc_d(length);
        *newtype = VF_REAL;
        for (i = 0; i < length; i++)
            d[i] = cosh(degtorad(dd[i]));
        return ((void *) d);
    }
}



static double *d_tan(double *dd, int length)
{
    int xrc = 0;
    double *d = alloc_d(length);

    int i;
    for (i = 0; i < length; i++) {
        rcheck(tan(degtorad(dd[i])) != 0, "tan");
        d[i] = tan(degtorad(dd[i]));
    }

EXITPOINT:
    if (xrc != 0) { /* Free resources on error */
        txfree(d);
        d = (double *) NULL;
    }

    return d;
} /* end of function d_tan */



static double *
d_tanh(double *dd, int length)
{
    double *d;
    int i;

    d = alloc_d(length);
    for (i = 0; i < length; i++) {
        d[i] = tanh(degtorad(dd[i]));
    }
    return d;
}

/* See https://proofwiki.org/wiki/Tangent_of_Complex_Number (formulation 4) among
 * others for the tangent formula:

 * sin z = sin(x + iy) = sin x cos(iy) + cos x sin(iy) = sin x cosh y + i cos x sinh y
 * cos z = cos(x + iy) = cos x cos(iy) + sin x sin(iy) = cos x cosh y - i sin x sinh y
 * tan z = ((sin x cosh y + i cos x sinh y) / (cos x cosh y - i sin x sinh y)) *
            (cos x cosh y + isin x sinh y) / (cos x cosh y + i sin x sinh y)
      = ...
 *
 *
 * tan(a + bi) = (sin(2a) + i * sinh(2b)) / (cos(2a) + cosh(2b))
 */
static ngcomplex_t *c_tan(ngcomplex_t *cc, int length)
{
    ngcomplex_t * const c = alloc_c(length);

    int i;
    for (i = 0; i < length; i++) {
        errno = 0;
        ngcomplex_t *p_dst = c + i;
        ngcomplex_t *p_src = cc + i;
        const double a = p_src->cx_real;
        const double b = p_src->cx_imag;
        const double u = 2 * degtorad(a);
        const double v = 2 * degtorad(b);
        const double n_r = sin(u);
        const double n_i = sinh(v);
        const double d1 = cos(u);
        const double d2 = cosh(v);
        const double d = d1 + d2;
        if (errno != 0 || d == 0.0) {
            (void) fprintf(cp_err,
                    "Invalid argument %lf + %lf i for compex tangent", a, b);
            txfree(c);
            return (ngcomplex_t *) NULL;
        }
        p_dst->cx_real = n_r / d;
        p_dst->cx_imag = n_i / d;
    } /* end of loop over elements in array */
    return c;
} /* end of function c_tan */



/* complex tanh function, uses tanh(z) = -i * tan(i * z).
 * Unfortunately, it did so very inefficiently (a memory allocations per
 * value calculated!) and leaked memory badly (all of the aforementioned
 * allocations.) These issues were fixed 2020-03-04 */
static ngcomplex_t *c_tanh(ngcomplex_t *cc, int length)
{
    ngcomplex_t * const tmp = alloc_c(length); /* i * z */

    /* Build the i * z array to allow tan() to be called */
    {
        int i;
        for (i = 0; i < length; ++i) {
            ngcomplex_t *p_dst = tmp + i;
            ngcomplex_t *p_src = cc + i;

            /* multiply by i */
            p_dst->cx_real = -p_src->cx_imag;
            p_dst->cx_imag = p_src->cx_real;
        }
    }

   /* Calculat tan(i * z), exiting on failure */
    ngcomplex_t *const c = c_tan(tmp, length);
    if (c == (ngcomplex_t *) NULL) {
        txfree(tmp);
        return (ngcomplex_t *) NULL;
    }

    /* Multiply by -i to find final result */
    {
        int i;
        for (i = 0; i < length; ++i) {
            ngcomplex_t *p_cur = c + i;
            const double cx_real = p_cur->cx_real;
            p_cur->cx_real = p_cur->cx_imag;
            p_cur->cx_imag = -cx_real;
        }
    }

    return c;
} /* end of function c_tanh */



void *
cx_tan(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_REAL) {
        *newtype = VF_REAL;
        return (void *) d_tan((double *) data, length);
    } else {
        *newtype = VF_COMPLEX;
        return (void *) c_tan((ngcomplex_t *) data, length);
    }
}

void *
cx_tanh(void *data, short int type, int length, int *newlength, short int *newtype)
{
    *newlength = length;
    if (type == VF_REAL) {
        *newtype = VF_REAL;
        return (void *) d_tanh((double *) data, length);
    } else {
        *newtype = VF_COMPLEX;
        return (void *) c_tanh((ngcomplex_t *) data, length);
    }
}

void *
cx_atan(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d;

    d = alloc_d(length);
    *newtype = VF_REAL;
    *newlength = length;
    if (type == VF_COMPLEX) {
        ngcomplex_t *cc = (ngcomplex_t *) data;
        int i;

        for (i = 0; i < length; i++)
            d[i] = radtodeg(atan(realpart(cc[i])));
    } else {
        double *dd = (double *) data;
        int i;

        for (i = 0; i < length; i++)
            d[i] = radtodeg(atan(dd[i]));
    }
    return ((void *) d);
}

/* Struct to store and order the values of the amplitudes preserving the index in the original array */
typedef struct {
    double amplitude;
    int index;
} amplitude_index_t;

static int compare_structs (const void *a, const void *b);

/*
 *  Returns the positions of the elements in a real vector
 *  after they have been sorted into increasing order using a stable method (qsort).
 */

void *
cx_sortorder(void *data, short int type, int length, int *newlength, short int *newtype)
{
    double *d = alloc_d(length);
    double *dd = (double *) data;
    int i;

    amplitude_index_t * const array_amplitudes = (amplitude_index_t *)
            tmalloc(sizeof(amplitude_index_t) * (size_t) length);

    *newlength = length;
    *newtype = VF_REAL;
    if (type == VF_REAL) {

        for(i = 0; i < length; i++){
            array_amplitudes[i].amplitude = dd[i];
            array_amplitudes[i].index = i;
        }

        qsort(array_amplitudes, (size_t) length, sizeof(array_amplitudes[0]), compare_structs);

        for(i = 0; i < length; i++)
            d[i] = array_amplitudes[i].index;
    }

    txfree(array_amplitudes);

    /* Otherwise it is 0, but tmalloc zeros the stuff already. */
    return ((void *) d);
}

static int
compare_structs(const void *a, const void *b)
{
    amplitude_index_t *aa = (amplitude_index_t *) a;
    amplitude_index_t *bb = (amplitude_index_t *) b;

    if (aa->amplitude > bb->amplitude)
        return 1;
    else if (aa->amplitude == bb->amplitude)
        return 0;
    else
        return -1;
}
