/*
 *    CCL (Complex C Library) - part of Mathlib
 */

#include "mathlib.constants"
#include "mathlib.errors"

typedef struct complex {
    double real;
    double imag;
} COMPLEX;

/*************************************************************************
 *    Cadd(z1,z2) = z1 + z2
 */
COMPLEX *Cadd(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;

    ztmp.real = z1->real + z2->real;
    ztmp.imag = z1->imag + z2->imag;
    return(&ztmp);
}

/*************************************************************************
 *    Csub(z1,z2) = z1 - z2
 */
COMPLEX *Csub(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;

    ztmp.real = z1->real - z2->real;
    ztmp.imag = z1->imag - z2->imag;
    return(&ztmp);
}

/*************************************************************************
 *    Cmult(z1,z2) = z1 * z2
 */
COMPLEX *Cmult(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;

    ztmp.real = z1->real * z2->real - z1->imag * z2->imag;
    ztmp.imag = z1->real * z2->imag + z2->real * z1->imag;
    return(&ztmp);
}

/*************************************************************************
 *    Cmultd(z,d) = z * d
 */
COMPLEX *Cmultd(z,d)
COMPLEX *z;
double d;
{
    static COMPLEX ztmp;

    ztmp.real = z->real * d;
    ztmp.imag = z->imag * d;
    return(&ztmp);
}

/*************************************************************************
 *    Cscale(z,d1,d2) = real * d1 + j imag * d2
 *    where z = real + j imag
 */
COMPLEX *Cscale(z,d1,d2)
COMPLEX *z;
double d1,d2;
{
    static COMPLEX ztmp;

    ztmp.real = z->real * d1;
    ztmp.imag = z->imag * d2;
    return(&ztmp);
}

/*************************************************************************
 *    Cdiv(z1,z2) = z1 / z2
 */
COMPLEX *Cdiv(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;

    double den,r;
    double absr,absi;
    double Dabs();

    absr = Dabs(z2->real);
    absi = Dabs(z2->imag);

    if (z1->real == 0.0 && z1->imag == 0.0) {
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else if (z2->real == 0.0 && z2->imag == 0.0) {
        ccl_error = ZERO_DIVIDE;
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else if (absr >= absi) {
        r = z2->imag / z2->real;
        den = z2->real + r * z2->imag;
        ztmp.real = (z1->real + z1->imag * r) / den;
        ztmp.imag = (z1->imag - z1->real * r) / den;
    } else {
        r = z2->real / z2->imag;
        den = z2->imag + r * z2->real;
        ztmp.real = (z1->real * r + z1->imag) / den;
        ztmp.imag = (z1->imag * r - z1->real) / den;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cdivd(z1,d) = z / d
 */
COMPLEX *Cdivd(z,d)
COMPLEX *z;
double d;
{
    static COMPLEX ztmp;

    if (d == 0.0) {
        ccl_error = ZERO_DIVIDE;
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else if (z->real == 0.0 && z->imag == 0.0) {
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else {
        ztmp.real = z->real / d;
        ztmp.imag = z->imag / d;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Crcp(z) = 1 / z
 */
COMPLEX *Crcp(z)
COMPLEX *z;
{
    static COMPLEX ztmp;

    double denom;

    if (z->real == 0.0 && z->imag == 0.0) {
        ccl_error = ZERO_DIVIDE;
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else if (z->imag == 0.0) {
        ztmp.real = 1.0 / z->real;
        ztmp.imag = 0.0;
    } else if (z->real == 0.0) {
        ztmp.real = 0.0;
        ztmp.imag = -1.0 / z->imag;
    } else {
        denom = z->real * z->real + z->imag * z->imag;
        ztmp.real =   z->real / denom;
        ztmp.imag = - z->imag / denom;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cpow(z1,z2) = z1 ^ z2
 */
COMPLEX *Cpow(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;
    double Datan(),Dln(),Dsin(),Dcos(),Dexp(),Dsign();
    double p,r,v,w,sign;

    sign = Dsign(z1->imag);

    if (z1->real == 0.0) {
        if (z1->imag == 0.0) {
            ztmp.real = 0.0;
            ztmp.imag = 0.0;
            return(&ztmp);
        } else
            p = PID2 * sign;
    }
    else {
        p = Datan(z1->imag/z1->real);
        if (z1->real < 0.0) {
            if (z1->imag >= 0.0)
                p += PI;
            else
                p -= PI;
        }
    }
    r = 0.5 * Dln(z1->real * z1->real + z1->imag * z1->imag);
    v = z2->real * p + z2->imag * r;
    w = Dexp(z2->real * r - z2->imag * p);

    ztmp.real = w * Dcos(v);
    ztmp.imag = w * Dsin(v);
    return(&ztmp);
}

/*************************************************************************
 *    Cpowd(z,d) = z ^ d
 */
COMPLEX *Cpowd(z,d)
COMPLEX *z;
double d;
{
    static COMPLEX ztmp;
    double Cabs();
    double Datan(),Dln(),Dcos(),Dsin(),Dexp(),Dabs();

    double phi,r,absr,absi,sqrr,sqri;

    if (z->real > 0.0)
        phi = Datan(z->imag/z->real);
    else if (z->real < 0.0 && z->imag >= 0.0)
        phi = Datan(z->imag/z->real) + PI;
    else if (z->real < 0.0 && z->imag < 0.0)
        phi = Datan(z->imag/z->real) - PI;
    else if (z->real == 0.0 && z->imag == 0.0) {
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
        return(&ztmp);
    } else if (z->real == 0.0 && z->imag > 0.0)
        phi = PID2;
    else if (z->real == 0.0 && z->imag < 0.0)
        phi = -PID2;

    absr = Dabs(z->real);
    absi = Dabs(z->imag);

    r = Cabs(z);
    r = Dexp(d * Dln(r));
    phi = d * phi;

    ztmp.real = r * Dcos(phi);
    ztmp.imag = r * Dsin(phi);
    return(&ztmp);
}

/*************************************************************************
 *    Csqrt(z) = u (+-) jv
 *    where u = sqrt(0.5*(root+x)) and v = sqrt(0.5*(root-x))
 *    and root is the magnitude of z
 *    the sign is the same as that of the imaginary part of z
 */
COMPLEX *Csqrt(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Cabs();
    double Dsqrt();
    double root;

    if (z->real == 0.0 && z->imag == 0.0) {
        ztmp.real = 0.0;
        ztmp.imag = 0.0;
    } else if (z->imag == 0.0) {
        ztmp.real = Dsqrt(z->real);
        ztmp.imag = 0.0;
    } else {
        root = Cabs(z);
        ztmp.real = Dsqrt(0.5 * (root + z->real));
        ztmp.imag = Dsqrt(0.5 * (root - z->real));
        if (z->imag < 0.0)
            ztmp.imag = - ztmp.imag;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cln(z) = ln( Cabs(z) ) + j Carg(z)
 *    where z = real + j imag
 */
COMPLEX *Cln(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Cabs(),Carg();
    double Dln();

    if (z->imag == 0.0 && z->real > 0.0) {
        ztmp.real = Dln(z->real);
        ztmp.imag = 0.0;
    } else if (z->real == 0.0) {
        if (z->imag > 0.0) {
            ztmp.real = Dln(z->imag);
            ztmp.imag = PID2;
        } else {
            ztmp.real = Dln(-(z->imag));
            ztmp.imag = - PID2;
        }
    } else {
        ztmp.real = Dln(Cabs(z));
        ztmp.imag = Carg(z);
    }
    return(&ztmp);
}

/*************************************************************************
 *    Clogd(z,d) = ln(z) / ln(d)
 */
COMPLEX *Clogd(z,d)
COMPLEX *z;
double d;
{
    static COMPLEX ztmp;
    COMPLEX *Cln();
    double Dln();
    double lnd;

    lnd = Dln(d);
    ztmp = *Cln(z);
    ztmp.real /= lnd;
    ztmp.imag /= lnd;
    return(&ztmp);
}

/*************************************************************************
 *    Cexp(z) = exp(real) cos(imag) + j( exp(real) sin(imag) )
 *    where z = real + j imag
 */
COMPLEX *Cexp(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dsin(),Dcos(),Dexp();
    double r;

    if (z->real == 0.0 && z->imag == 0.0) {
        ztmp.real = 1.0;
        ztmp.imag = 0.0;
    } else if (z->imag == 0.0) {
        ztmp.real = Dexp(z->real);
        ztmp.imag = 0.0;
    } else if (z->real == 0.0) {
        ztmp.real = Dcos(z->imag);
        ztmp.imag = Dsin(z->imag);
    } else {
        r = Dexp(z->real);
        ztmp.real = r * Dcos(z->imag);
        ztmp.imag = r * Dsin(z->imag);
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cconj(z) = real - j imag
 *    where z = real + j imag
 */
COMPLEX *Cconj(z)
COMPLEX *z;
{
    static COMPLEX ztmp;

    ztmp.real = z->real;
    ztmp.imag = - z->imag;
    return(&ztmp);
}

/*************************************************************************
 *    Cxj(z) = j z
 */
COMPLEX *Cxj(z)
COMPLEX *z;
{
    static COMPLEX ztmp;

    ztmp.real = - z->imag;
    ztmp.imag = z->real;
    return(&ztmp);
}

/*************************************************************************
 *    Crotate(z,d) = r exp[ j (phi + d) ]
 */
COMPLEX *Crotate(z,d)
COMPLEX *z;
double d;
{
    static COMPLEX ztmp;
    double Dsin(),Dcos();
    double sind,cosd;

    sind = Dsin(d);
    cosd = Dcos(d);
    ztmp.real = z->real * cosd - z->imag * sind;
    ztmp.imag = z->real * sind + z->imag * cosd;
    return(&ztmp);
}

/*************************************************************************
 *    Cmax(z1,z2) = max(z1,z2)
 */
COMPLEX *Cmax(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;
    double abs1,abs2;

    abs1 = z1->real * z1->real + z1->imag * z1->imag;
    abs2 = z2->real * z2->real + z2->imag * z2->imag;
    if (abs1 <= abs2) {
        ztmp.real = z2->real;
        ztmp.imag = z2->imag;
    } else {
        ztmp.real = z1->real;
        ztmp.imag = z1->imag;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cmin(z1,z2) = min(z1,z2)
 */
COMPLEX *Cmin(z1,z2)
COMPLEX *z1,*z2;
{
    static COMPLEX ztmp;
    double abs1,abs2;

    abs1 = z1->real * z1->real + z1->imag * z1->imag;
    abs2 = z2->real * z2->real + z2->imag * z2->imag;
    if (abs1 >= abs2) {
        ztmp.real = z2->real;
        ztmp.imag = z2->imag;
    } else {
        ztmp.real = z1->real;
        ztmp.imag = z1->imag;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cswap(z1,z2), swaps z1 and z2
 */
int Cswap(z1,z2)
COMPLEX *z1,*z2;
{
    COMPLEX temp;

    temp.real = z2->real;
    temp.imag = z2->imag;
    z2->real = z1->real;
    z2->imag = z1->imag;
    z1->real = temp.real;
    z1->imag = temp.imag;
}

/*************************************************************************
 *    Creal(z) = real
 *    where z = real + j imag
 */
double Creal(z)
COMPLEX *z;
{
    return(z->real);
}

/*************************************************************************
 *    Cimag(z) = imag
 *    where z = real + j imag
 */
double Cimag(z)
COMPLEX *z;
{
    return(z->imag);
}

/*************************************************************************
 *    Cabs(z) = sqrt( real^2 + imag^2 )
 *    where z = real + j imag
 */
double Cabs(z)
COMPLEX *z;
{
    double Dsqrt(),Dabs();
    double absr,absi,sqrr,sqri;

    absr = Dabs(z->real);
    absi = Dabs(z->imag);

    if (absr == 0.0)
        return(absi);
    else if (absi == 0.0)
        return(absr);
    else if (absr > absi) {
        sqrr = absr * absr;
        sqri = absi * absi;
        return(absr * Dsqrt(1 + sqri/sqrr));
    } else {
        sqrr = absr * absr;
        sqri = absi * absi;
        return(absi * Dsqrt(1 + sqrr/sqri));
    }
}

/*************************************************************************
 *    Carg(z) Return the angle of z to the positive real axis
 */
double Carg(z)
COMPLEX *z;
{
    double Datan(),Dsign();
    double sign;

    sign = Dsign(z->imag);

    if (z->real > 0.0)
        return(Datan(z->imag / z->real));
    else if (z->real < 0.0) {
        if (z->imag != 0.0)
            return(Datan(z->imag / z->real) + sign * PI);
        else
            return(-PI);
    } else {
        if (z->imag != 0.0)
            return(sign * PID2);
        else
            return(0.0);
    }
}

/*************************************************************************
 *    Cmplx(d1,d2) = z
 *    where z = d1 + j d2
 */
COMPLEX *Cmplx(d1,d2)
double d1,d2;
{
    static COMPLEX ztmp;

    ztmp.real = d1;
    ztmp.imag = d2;
    return(&ztmp);
}

/*************************************************************************
 *    Cmake(d1,d2) = z
 *    where z + d1 ( cos(d2) + j sin(d2) )
 */
COMPLEX *Cmake(d1,d2)
double d1,d2;
{
    double Dsin(),Dcos();
    static COMPLEX ztmp;

    ztmp.real = d1 * Dcos(d2);
    ztmp.imag = d1 * Dsin(d2);
    return(&ztmp);
}

/*************************************************************************
 *    Csin(z) = sin(real) cosh(imag) + j cos(real) sinh(imag)
 *    where z = real + j imag
 */
COMPLEX *Csin(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dsin(),Dcos(),Dsinh(),Dcosh();

    if (z->imag == 0.0) {
        ztmp.real = Dsin(z->real);
        ztmp.imag = 0.0;
    } else {
        ztmp.real = Dsin(z->real) * Dcosh(z->imag);
        ztmp.imag = Dcos(z->real) * Dsinh(z->imag);
    }
    return(&ztmp);
}

/*************************************************************************
 *    Ccos(z) = cos(real) cosh(imag) - j sin(real) sinh(imag)
 */
COMPLEX *Ccos(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dsin(),Dcos(),Dsinh(),Dcosh();

    if (z->imag == 0.0) {
        ztmp.real = Dcos(z->real);
        ztmp.imag = 0.0;
    } else {
        ztmp.real = Dcos(z->real) * Dcosh(z->imag);
        ztmp.imag = - Dsin(z->real) * Dsinh(z->imag);
    }
    return(&ztmp);
}

/*************************************************************************
 *    Ctan(z) = ( sin(2*real) + jsinh(2*imag) )
 *              -------------------------------
 *              ( cos(2*real) + cosh(2*imag) )
 *    where z = real + j imag
 */
COMPLEX *Ctan(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dtan(),Dcos(),Dsin(),Dsinh(),Dcosh();
    double denom,real2,imag2;

    if (z->imag == 0.0) {
        ztmp.real = Dtan(z->real);
        ztmp.imag = 0.0;
    } else {
        real2 = 2.0 * z->real;
        imag2 = 2.0 * z->imag;
        denom = Dcos(real2) + Dcosh(imag2);
        ztmp.real = Dsin(real2) / denom;
        ztmp.imag = Dsinh(imag2) / denom;
    }
    return(&ztmp);
}

/*************************************************************************
 *    Casin(z) = k*pi + (-1)^k asin(b)
 *                    + j (-1)^k ln(a + sqrt(a^2 - 1))
 *    where a = 0.5 sqrt((x+1)^2 + y^2) + 0.5 sqrt((x-1)^2 + y^2)
 *    and   b = 0.5 sqrt((x+1)^2 + y^2) - 0.5 sqrt((x-1)^2 + y^2)
 *    and z = x + jy, k an integer
 */
COMPLEX *Casin(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dsqrt(),Dasin(),Dln();
    double a,b;
    double xm1,xp1,x2,y2;
    double part1,part2;

    if (z->imag == 0.0) {
        ztmp.real = Dasin(z->real);
        ztmp.imag = 0.0;
    } else {
        x2 = z->real * z->real;
        y2 = z->imag * z->imag;
        xp1 = x2 + 2.0 * z->real + 1.0;
        xm1 = x2 - 2.0 * z->real + 1.0;
        part1 = 0.5 * Dsqrt(xp1 + y2);
        part2 = 0.5 * Dsqrt(xm1 + y2);
        a = part1 + part2;
        b = part1 - part2;
        ztmp.real = Dasin(b);
        ztmp.imag = Dln( a + Dsqrt(a * a - 1.0) );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Cacos(z) = 2*k*pi (+-) [ acos(b)
 *                        - j ln(a + sqrt(a^2 - 1))
 *    where a = 0.5 sqrt((x+1)^2 + y^2) + 0.5 sqrt((x-1)^2 + y^2)
 *    and   b = 0.5 sqrt((x+1)^2 + y^2) - 0.5 sqrt((x-1)^2 + y^2)
 *    and   z = x + jy, K an integer.
 */
COMPLEX *Cacos(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Dsqrt(),Dacos(),Dln(),Dsqrt();
    double a,b;
    double xm1,xp1,x2,y2;
    double part1,part2;

    if (z->imag == 0.0) {
        ztmp.real = Dacos(z->real);
        ztmp.imag = 0.0;
    } else {
        x2 = z->real * z->real;
        y2 = z->imag * z->imag;
        xp1 = x2 + 2.0 * z->real + 1.0;
        xm1 = x2 - 2.0 * z->real + 1.0;
        part1 = 0.5 * Dsqrt(xp1 + y2);
        part2 = 0.5 * Dsqrt(xm1 + y2);
        a = part1 + part2;
        b = part1 - part2;
        ztmp.real = Dacos(b);
        ztmp.imag = - Dln( a + Dsqrt(a*a - 1.0) );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Catan(z) = k*pi + 0.5 * atan(2x/(1-x^2-y^2))
 *                    + j/4 * ln (( x^2+(y+1)^2) / ( x^2+(y-1)^2))
 *    where z = x + jy, K an integer.
 */
COMPLEX *Catan(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    double Datan(),Dln();
    double ym1,yp1,x2,y2,denom;

    if (z->imag == 0.0) {
        ztmp.real = Datan(z->real);
        ztmp.imag = 0.0;
    } else {
        x2 = z->real * z->real;
        y2 = z->imag * z->imag;
        denom = 1.0 - x2 - y2;
        yp1 = x2 + y2 + 2.0 * z->imag + 1.0;
        ym1 = x2 + y2 - 2.0 * z->imag + 1.0;
        ztmp.real = 0.5 * Datan( 2.0 * z->real / denom );
        ztmp.imag = 0.25 * Dln( yp1 / ym1 );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Csinh(z) = 0.5 ( Cexp(z) - Cexp(-z) )
 */
COMPLEX *Csinh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX *Cexp();
    double Dsinh();
    COMPLEX mz,zt1,zt2;

    if (z->imag == 0.0) {
        ztmp.real = Dsinh(z->real);
        ztmp.imag = 0.0;
    } else {
        mz.real = - z->real;
        mz.imag = - z->imag;
        zt1 = *Cexp(z);
        zt2 = *Cexp(&mz);
        ztmp.real = 0.5 * (zt1.real - zt2.real );
        ztmp.imag = 0.5 * (zt1.imag - zt2.imag );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Ccosh(z) = 0.5 ( Cexp(z) + Cexp(-z) )
 */
COMPLEX *Ccosh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX *Cexp();
    COMPLEX mz,zt1,zt2;
    double Dcosh();

    if (z->imag == 0.0) {
        ztmp.real = Dcosh(z->real);
        ztmp.imag = 0.0;
    } else {
        mz.real = - z->real;
        mz.imag = - z->imag;
        zt1 = *Cexp(z);
        zt2 = *Cexp(&mz);
        ztmp.real = 0.5 * ( zt1.real + zt2.real );
        ztmp.imag = 0.5 * ( zt1.imag + zt2.imag );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Ctanh(z) = ( 1 - Cexp(-2z) ) / ( 1 + Cexp(-2z) )
 */
COMPLEX *Ctanh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX zt1,zt2,num,denom;
    COMPLEX *Cdiv(),*Cexp();
    double Dtanh();

    if (z->imag == 0.0) {
        ztmp.real = Dtanh(z->real);
        ztmp.imag = 0.0;
    } else {
        zt1.real = -2.0 * z->real;
        zt1.imag = -2.0 * z->imag;
        zt2 = *Cexp(&zt1);
        num.real = 1.0 - zt2.real;
        num.imag = - zt2.imag;
        denom.real = 1.0 + zt2.real;
        denom.imag = zt2.imag;
        ztmp = *Cdiv( &num , &denom );
    }
    return(&ztmp);
}

/*************************************************************************
 *    Casinh(z) = Cln( z + Csqrt( z^2 + 1 ))
 */
COMPLEX *Casinh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX zt1,zt2;
    COMPLEX *Cln(),*Csqrt();

    zt1.real = z->real * z->real - z->imag * z->imag + 1.0;
    zt1.imag = 2.0 * z->real * z->imag;
    zt2 = *Csqrt(&zt1);
    zt2.real += z->real;
    zt2.imag += z->imag;
    ztmp = *Cln(&zt2);
    return(&ztmp);
}

/*************************************************************************
 *    Cacosh(z) = Cln ( z + Csqrt(z^2 - 1) )
 */
COMPLEX *Cacosh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX zt1,zt2;
    COMPLEX *Csqrt(),*Cln();

    zt1.real = z->real * z->real - z->imag * z->imag - 1.0;
    zt1.imag = 2.0 * z->real * z->imag;
    zt2 = *Csqrt(&zt1);
    zt2.real += z->real;
    zt2.imag += z->imag;
    ztmp = *Cln(&zt2);
    return(&ztmp);
}

/*************************************************************************
 *    Catanh(z) = 0.5 * Cln( (1+z) / (1-z) )
 */
COMPLEX *Catanh(z)
COMPLEX *z;
{
    static COMPLEX ztmp;
    COMPLEX zp1,zm1,zt1;
    COMPLEX *Cln(),*Cdiv();

    zp1.real = 1.0 + z->real;
    zp1.imag = z->imag;
    zm1.real = 1.0 - z->real;
    zm1.imag = - (z->imag);
    zt1 = *Cln( Cdiv( &zp1 , &zm1 ) );
    ztmp.real = zt1.real * 0.5;
    ztmp.real = zt1.imag * 0.5;
    return(&ztmp);
}

/*************************************************************************
 *    Cuniform() = z
 *    where the magnitude and phase of z are uniformly distributed
 *    between (0,sqrt 2) and (0,2 pi) respectively
 */
COMPLEX *Cuniform()
{
    static COMPLEX ztmp;
    double Dsin(),Dcos(),Duniform();
    double d1,d2;

    d1 = SQRT2 * Duniform();
    d2 = TWOPI * Duniform();
    ztmp.real = d1 * Dcos(d2);
    ztmp.imag = d1 * Dsin(d2);
    return(&ztmp);
}

/*************************************************************************
 *    Cnormal() = z
 *    where the real and imaginary parts of z are normally
 *    distributed with mean 0 and variance 1, ie: N(0,1)
 *
 *    Given two samples u1 and u2 from a uniform distribution
 *    on (0,1)
 *    then    n1 = sqrt(2 ln(u1)) cos(2 pi u2)
 *    and     n2 = sqrt(2 ln(u1)) sin(2 pi u2)
 *    are two independent normally distributed samples ie:  N(0,1)
 */
COMPLEX *Cnormal()
{
    static COMPLEX ztmp;
    double Dln(),Dsqrt(),Dcos(),Dsin(),Duniform();
    double d1,d2,rand1;

    while ( (rand1 = Duniform()) == 0.0) ;
    d1 = Dsqrt( 2.0 * Dln(rand1) );
    d2 = TWOPI * ran();
    ztmp.real = d1 * Dcos(d2);
    ztmp.imag = d1 * Dsin(d2);
    return(&ztmp);
}

/*************************************************************************
 *    Crayleigh() = z
 *    where the real and imaginary parts are uniformly
 *    distributed on (-1,1)
 */
COMPLEX *Crayleigh()
{
    static COMPLEX ztmp;
    double Duniform();

    ztmp.real = 2.0 * Duniform() - 1.0;
    ztmp.imag = 2.0 * Duniform() - 1.0;
    return(&ztmp);
}

/*************************************************************************
 *    Cpoly(z1,n,z) evaluates a polynomial at z1 of degree n
 *    whose coefficients are in the complex array z.
 *
 *    f(z) = a0 + a1 z + a2 z^2 + .... + an z^n
 *         = a0 + z ( a1 + z ( a2 + .... + z an) ) ... )
 *    The a's are complex coeficients in increasing order.
 *    This function evaluates f(z1)
 */
COMPLEX *Cpoly(z1,n,z)
COMPLEX *z1;
int n;
COMPLEX z[];
{
    static COMPLEX ztmp;
    COMPLEX zp;
    int i;

    ztmp.real = z1->real;
    ztmp.imag = z1->imag;
    for (i = n; i > 0; i--) {
        zp.real = ztmp.real * z[i].real - ztmp.imag * z[i].imag;
        zp.imag = ztmp.real * z[i].imag + ztmp.imag * z[i].real;
        ztmp.real += z[i-1].real;
        ztmp.imag += z[i-1].imag;
    }
    return(&ztmp);
}