text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <stdbool.h>
#include <png.h>
#include "beachball.h"
#include "parmt_config.h"
#ifdef PARMT_USE_INTEL
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
//#include "pixmap.h"
//#include "pixmap_png.h"
//#include "pixmap_jpg.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661923
#endif
#ifndef MAX
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
#endif
#ifndef MIN
#define MIN(x,y) (((x) < (y)) ? (x) : (y))
#endif
#define PIXMAP_COLORS 3

static void vecstrd(const double str, const double dip, double r[3]);
static void strdv(const double r[3], double *str, double *dip);
static void crossp(const double *__restrict__ u,
                   const double *__restrict__ v,
                   double *__restrict__ n);
static void strdvec(double r[3], double *str, double *dip);
static void rotmat(const int naxis, const double theta,
                   double *__restrict__ r,
                   double *__restrict__ rt);
static void isItNear(const double *__restrict__ pt,
                     double *howClose,
                     double *__restrict__ ptold);
static double dot3(const double *__restrict__ x,
                   const double *__restrict__ y);
static void gemv3(const double *__restrict__ A, const double *__restrict__ x,
                  double *__restrict__ y);
static double transform(const double a, const double b,
                        const double c, const double d,
                        const double x);

/*
int beachball_setMomentTensor(const double m11, const double m22,
                              const double m33, const double m12,
                              const double m13, const double m23,
                              const enum cmopad_basis_enum basis,
                              struct beachballPlot_struct *beachball)
{
    struct cmopad_struct cmt;
    double mt[3][3], m9[9], eigs[3];
    int ierr, info;
    bool lexplosion, limplosion;
    memset(&cmt, 0, sizeof(struct cmopad_struct));
    // the edge case is for purely isotropic sources
    lexplosion = false;
    limplosion = false;
    m9[0] = m11;
    m9[1] = m12;
    m9[2] = m13;
    m9[3] = m12;
    m9[4] = m22;
    m9[5] = m23;
    m9[6] = m13;
    m9[7] = m23;
    m9[8] = m33;
    info = LAPACKE_dsyev(LAPACK_COL_MAJOR, 'N', 'L', 3, m9, 3, eigs);
    if (info != 0)
    {
        printf("Failed classifying eigenvalues\n");
    }
    if ((eigs[0] > 0.0 && eigs[1] > 0.0 && eigs[2] > 0.0))
    {
        //printf("Source is pure explosion\n");
        lexplosion = true;
    }
    if ((eigs[0] < 0.0 && eigs[1] < 0.0 && eigs[2] < 0.0))
    {
        printf("Source is pure implosion\n");
        limplosion = true;
    }
    // copy the input moment tensor
    mt[0][0] = m11;
    mt[1][1] = m22;
    mt[2][2] = m33;
    mt[0][1] = mt[1][0] = m12;
    mt[0][2] = mt[2][0] = m13;
    mt[1][2] = mt[2][1] = m23;
    ierr = cmopad_basis_transformMatrixM33(mt, basis, NED);
    if (ierr != 0)
    {
        printf("Error transforming matrix\n");
        return -1;
    }
    // decompose the moment tensor
    ierr = cmopad_standardDecomposition(mt, &cmt);
    if (ierr != 0)
    {
        printf("Failed initializing moment tensor from strike, dip, rake\n");
        return -1;
    }
    // get the principal tension, null, and plunge axes
    ierr = cmopad_MT2PrincipalAxisSystem(0, &cmt);
    if (ierr != 0)
    {
        printf("Error computing principal axes\n");
        return -1;
    }
    if (lexplosion || limplosion)
    {
        cmt.p_axis[2] = cmt.eigenvalues[0];
        cmt.null_axis[2] = cmt.eigenvalues[1];
        cmt.t_axis[2] = cmt.eigenvalues[2];
    }
    else
    {
        // Convert pressure, null, and tension principal axes to
        // azimuth, plunge, and length
        ierr += cmopad_Eigenvector2PrincipalAxis(NED, cmt.eig_pnt[0],
                                                 cmt.p_axis,    cmt.p_axis);
        ierr += cmopad_Eigenvector2PrincipalAxis(NED, cmt.eig_pnt[1],
                                                 cmt.null_axis, cmt.null_axis);
        ierr += cmopad_Eigenvector2PrincipalAxis(NED, cmt.eig_pnt[2],
                                                 cmt.t_axis,    cmt.t_axis);
        if (ierr != 0)
        {
            printf("Error converting eigenvector to principal axis!\n");
            return -1;
        }
    }
    memcpy(&beachball->cmt, &cmt, sizeof(struct cmopad_struct));
    beachball->lhaveMT = true;
    return 0;
}
*/

int beachball_writePNG(const char *fileName,
                       size_t height, size_t width, unsigned char *bytes)
{
    FILE *file = NULL;
    png_structp png = NULL;
    png_infop info = NULL;
    png_bytepp rows = NULL;
    int h;

    file = fopen(fileName, "wb");
    if (file == NULL){printf("failed to open file\n");}
    png = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
    if (png == NULL){printf("error making png struct\n");}
    info = png_create_info_struct(png);
    if (info == NULL){printf("error making png info struct\n");}
    png_init_io(png, file);
    png_set_IHDR(png, info, width, height, 8,
                 PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE,
                 PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE);
    png_write_info(png, info);
    rows = calloc(height, sizeof(png_bytep));
    for (h=0; h<(int) height; h++)
    {
        rows[h] = bytes + (PIXMAP_COLORS*h*(int) width);
    }
    png_write_rows(png, rows, height);
    png_write_end(png, NULL); //info);
    png_destroy_write_struct(&png, &info);
    if (rows != NULL){free(rows);}
    fclose(file);
    return 0;
}

int beachball_draw(struct beachballPlot_struct beachball)
{
    const char *fcnm = "beachball_draw\0";
    struct beachBall_struct beachImage;
    unsigned char *image;
    double baz, bdp, bev, paz, pdp, pev, taz, tdp, tev, x, xnorm, y;
    int i, i1, i2, indx, ix, iy, jy, k;
    size_t npixel;
    taz = beachball.tAxis[0];
    tdp = beachball.tAxis[1];
    tev = beachball.tAxis[2];
    baz = beachball.nAxis[0];
    bdp = beachball.nAxis[1];
    bev = beachball.nAxis[2];
    paz = beachball.pAxis[0];
    pdp = beachball.pAxis[1];
    pev = beachball.pAxis[2];
    // normalize the moment tensor for my sanity
    xnorm = MAX(fabs(tev), MAX(fabs(bev), fabs(pev)));
    if (xnorm == 0.0)
    {
        printf("%s: Error normalization is zero %f %f %f\n", fcnm, tev, bev, pev);
        return -1;
    }
    // compute everything describing the beachball
    memset(&beachImage, 0, sizeof(struct beachBall_struct));
    cliffsNodes_beachBallP(beachball.xc, beachball.yc, beachball.rad,
                           taz, tdp, tev,
                           baz, bdp, bev,
                           paz, pdp, pev,
                           beachball.nNodalLine,
                           beachball.nFocalSphere,
                           beachball.nxPixel, beachball.nyPixel,
                           true, true, true,
                           &beachImage);
    // draw the image
    npixel = (size_t) (beachball.nxPixel*beachball.nyPixel);
    image = (unsigned char *)
            calloc(PIXMAP_COLORS*npixel, sizeof(unsigned char));
    // Set a white background
    for (indx=0; indx<3*beachball.nxPixel*beachball.nyPixel; indx++)
    {
        image[indx] = 255;
    }
    // Draw the polarity
    if (beachball.lwantPolarity)
    {
        for (i=0; i<beachImage.npolarity; i++)
        {
            x = transform(beachball.xc - beachball.rad,
                          beachball.xc + beachball.rad,
                          0, (double) (beachball.nxPixel-1),
                          beachImage.xPolarity[i]);
            y = transform(beachball.yc - beachball.rad,
                          beachball.yc + beachball.rad,
                          0, (double) (beachball.nyPixel-1),
                          beachImage.yPolarity[i]);
            ix = (int) (x + 0.5);
            iy = (int) (y + 0.5);
            ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
            iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
            jy = beachball.nyPixel - 1 - iy;
            if (beachImage.polarity[i] ==-1)
            {
                // white
                image[3*(jy*beachball.nxPixel+ix)+0] = 255;
                image[3*(jy*beachball.nxPixel+ix)+1] = 255;
                image[3*(jy*beachball.nxPixel+ix)+2] = 255;
            }
            else
            {
                // red
                image[3*(jy*beachball.nxPixel+ix)+0] = 255;
                image[3*(jy*beachball.nxPixel+ix)+1] = 0;
                image[3*(jy*beachball.nxPixel+ix)+2] = 0;
            }
        }
    }
    // Add the nodal lines
    if (beachball.lwantNodalLines)
    {
        // draw first nodal line
        for (k=0; k<beachImage.nNodalLineSegs1; k++)
        {
            i1 = beachImage.nodalLine1Ptr[k];
            i2 = beachImage.nodalLine1Ptr[k+1];
            for (i=i1; i<i2; i++)
            {
                // start point in vector (line)
                x = transform(beachball.xc - beachball.rad,
                              beachball.xc + beachball.rad,
                              0, (double) (beachball.nxPixel-1),
                              beachImage.xNodalLine1[2*i]);
                y = transform(beachball.yc - beachball.rad,
                              beachball.yc + beachball.rad,
                              0, (double) (beachball.nyPixel-1),
                              beachImage.yNodalLine1[2*i]);
                ix = (int) (x + 0.5);
                iy = (int) (y + 0.5);
                ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
                iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
                jy = beachball.nyPixel - 1 - iy;
                // black
                image[3*(jy*beachball.nxPixel+ix)+0] = 0;
                image[3*(jy*beachball.nxPixel+ix)+1] = 0;
                image[3*(jy*beachball.nxPixel+ix)+2] = 0;
                // end point in vector (line)
                x = transform(beachball.xc - beachball.rad,
                              beachball.xc + beachball.rad,
                              0, (double) (beachball.nxPixel-1),
                              beachImage.xNodalLine1[2*i+1]);
                y = transform(beachball.yc - beachball.rad,
                              beachball.yc + beachball.rad,
                              0, (double) (beachball.nyPixel-1),
                              beachImage.yNodalLine1[2*i+1]);
                ix = (int) (x + 0.5);
                iy = (int) (y + 0.5);
                ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
                iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
                jy = beachball.nyPixel - 1 - iy;
                // black
                image[3*(jy*beachball.nxPixel+ix)+0] = 0;
                image[3*(jy*beachball.nxPixel+ix)+1] = 0;
                image[3*(jy*beachball.nxPixel+ix)+2] = 0;
            }
        }
        // draw second nodal line
        for (k=0; k<beachImage.nNodalLineSegs2; k++)
        {
            i1 = beachImage.nodalLine2Ptr[k];
            i2 = beachImage.nodalLine2Ptr[k+1];
            for (i=i1; i<i2; i++)
            {
                // start point in vector (line)
                x = transform(beachball.xc - beachball.rad,
                              beachball.xc + beachball.rad,
                              0, (double) (beachball.nxPixel-1),
                              beachImage.xNodalLine2[2*i]);
                y = transform(beachball.yc - beachball.rad,
                              beachball.yc + beachball.rad,
                              0, (double) (beachball.nyPixel-1),
                              beachImage.yNodalLine2[2*i]);
                ix = (int) (x + 0.5);
                iy = (int) (y + 0.5);
                ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
                iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
                jy = beachball.nyPixel - 1 - iy;
                // black
                image[3*(jy*beachball.nxPixel+ix)+0] = 0;
                image[3*(jy*beachball.nxPixel+ix)+1] = 0;
                image[3*(jy*beachball.nxPixel+ix)+2] = 0;
                // end point in vector (line)
                x = transform(beachball.xc - beachball.rad,
                              beachball.xc + beachball.rad,
                              0, (double) (beachball.nxPixel-1),
                              beachImage.xNodalLine2[2*i+1]);
                y = transform(beachball.yc - beachball.rad,
                              beachball.yc + beachball.rad,
                              0, (double) (beachball.nyPixel-1),
                              beachImage.yNodalLine2[2*i+1]);
                ix = (int) (x + 0.5);
                iy = (int) (y + 0.5);
                ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
                iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
                jy = beachball.nyPixel - 1 - iy;
                // black
                image[3*(jy*beachball.nxPixel+ix)+0] = 0;
                image[3*(jy*beachball.nxPixel+ix)+1] = 0;
                image[3*(jy*beachball.nxPixel+ix)+2] = 0;
            }
        }
    }
    // Draw an outline
    for (i=0; i<beachImage.nFocalSphere; i++)
    {
        x = transform(beachball.xc - beachball.rad,
                      beachball.xc + beachball.rad,
                      0, (double) (beachball.nxPixel-1),
                      beachImage.xFocalSphere[i]);
        y = transform(beachball.yc - beachball.rad,
                      beachball.yc + beachball.rad,
                      0, (double) (beachball.nyPixel-1),
                      beachImage.yFocalSphere[i]);
        ix = (int) (x + 0.5);
        iy = (int) (y + 0.5);
        ix = MAX(0, MIN(ix, beachball.nxPixel - 1));
        iy = MAX(0, MIN(iy, beachball.nyPixel - 1));
        jy = beachball.nyPixel - 1 - iy;
        // black
        image[3*(jy*beachball.nxPixel+ix)+0] = 0;
        image[3*(jy*beachball.nxPixel+ix)+1] = 0;
        image[3*(jy*beachball.nxPixel+ix)+2] = 0;
    }
    // Write it
    beachball_writePNG("temp2.png",
                       (size_t) beachball.nxPixel, (size_t) beachball.nyPixel,
                       image);
    free(image);
    cliffsNodes_free(&beachImage);
    return 0;
}

int beachball_plotDefaults(struct beachballPlot_struct *beachball)
{
    memset(beachball, 0, sizeof(struct beachballPlot_struct));
    beachball->rad = 1.0;
    beachball->xc = 1.5;
    beachball->yc = 1.5;
    beachball->nxPixel = 251;
    beachball->nyPixel = 251;
    beachball->nNodalLine = 3*beachball->nxPixel;
    beachball->nFocalSphere = 4*beachball->nxPixel;
    beachball->lwantPolarity = true;
    beachball->lwantNodalLines = true;
    beachball->lhaveMT = false;
    return 0;
}

static double transform(const double a, const double b,
                        const double c, const double d,
                        const double x)
{
    double det, c1, c2, xi;
    det = 1.0/(b - a);
    c1 = det*(b*c - a*d);
    c2 = det*(-c + d);
    xi = c1 + x*c2;
    return xi;
}
//============================================================================//
int cliffsNodes_plotStations(const double xc, const double yc, const double rad,
                             const int nstat,
                             const double *__restrict__ az,
                             const double *__restrict__ toa,
                             double *__restrict__ xp,
                             double *__restrict__ yp);

/*!
 * @brief Computes the station locations for the given azimuth
 *        and take-off angle on a focal sphere with radius rad
 *        centered at (xc, yc).  This is based on Cliff Frohlich's
 *        stnPlt but now processes multiple stations and does no plotting.
 *
 * @param[in] xc     x center of focal sphere
 * @param[in] yc     y center of focal sphere
 * @param[in] rad    radius of sphere
 * @param[in] nstat  number of stations
 * @param[in] az     source station azimuths (degrees) [nstat]
 * @param[in] toa    take off angle (degrees) from source to stations [nstat]
 *
 * @param[out] xp    station's x location on focal sphere [nstat]
 * @param[out] yp    station's y location on focal sphere [nstat]
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
int cliffsNodes_plotStations(const double xc, const double yc, const double rad,
                             const int nstat,
                             const double *__restrict__ az,
                             const double *__restrict__ toa,
                             double *__restrict__ xp,
                             double *__restrict__ yp)
{
    const char *fcnm = "cliffsNodes_plotStations\0";
    double azrad, rp, sign, toaN;
    int i;
    const double sqrt2 = 1.4142135623730951;
    const double degrad = M_PI/180.0;
    //------------------------------------------------------------------------//
    //
    // error handling
    if (nstat < 1 || az == NULL || toa == NULL || xp == NULL || yp == NULL)
    {
        if (nstat < 1){printf("%s: No stations\n", fcnm);}
        if (az == NULL){printf("%s: az is NULL\n", fcnm);}
        if (toa == NULL){printf("%s: toa is NULL\n", fcnm);}
        if (xp == NULL){printf("%s: xp is NULL\n", fcnm);}
        if (yp == NULL){printf("%s: yp is NULL\n", fcnm);}
        return -1;
    }
    // compute point on focal shere for each station
    for (i=0; i<nstat; i++)
    {
        // ray normally goes down
        sign = 1.0;
        toaN = toa[i];
        // ray is going up
        if (toaN > 90.)
        {
            toaN = 180. - toa[i];
            sign =-1.;
        }
        rp = sign*rad*sqrt2*sin(0.5*degrad*toaN);
        azrad = az[i]*degrad;
        xp[i] = xc + rp*sin(azrad);
        yp[i] = yc + rp*cos(azrad);
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Given the azimuth and dip of the tension and pressure eigenvectors
 *        in geographic coordinates along with normalized eigenvalues
 *        this computes the moment tensort in geographic coordinates as
 *        well as the rotation and transpose rotation matrices which rotate
 *        from a system where t=1-axis and p=2-axis.  This is based on
 *        Cliff Frohlich's mtensor.
 *
 * @param[in] taz    tension axis azimuth (degrees)
 * @param[in] tdp    tension axis dip (degrees)
 * @param[in] paz    pressure axis azimuth (degrees)
 * @param[in] pdp    pressure axis dip
 * @param[in] tev1   normalized eigenvalue for tension eigenvector
 * @param[in] bev1   normalized eigenvalue for null eigenvector
 * @param[in] pev1   normalized eigenvalue for pressure eigenvector
 *
 * @param[out] mt    [3 x 3] moment tensor (column major)
 * @param[out] rm    [3 x 3] rotation matrix as described above (column major)
 * @param[out] rmt   [3 x 3] transpose rotation matrix described above
 *                   (column major)
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
void cliffsNodes_mtensor(const double taz, const double tdp,
                         const double paz, const double pdp,
                         const double tev1, const double bev1,
                         const double pev1,
                         double *__restrict__ mt,
                         double *__restrict__ rm,
                         double *__restrict__ rmt,
                         double *ampMax)
{
    double mpt[9], sm1[9], smt1[9], sm2[9], smt2[9],  sm3[9],
           smt3[9], work[9], p[3], pnew[3], pxt[3], t[3], pang;
    const double degrad = M_PI/180.0;
    // set up moment tensor in t-p system
    memset(mpt, 0, 9*sizeof(double));
    mpt[0] = tev1;
    mpt[4] = pev1;
    mpt[8] = bev1;
    // condition p-axis so that it really is perpendicular to t-axis
    vecstrd(taz, tdp, t);
    vecstrd(paz, pdp, p);
    crossp(p, t, pxt);
    crossp(t, pxt, p);
    // Step 1: rotate 1-axis about 3-axis so that 1-axis lies along t azimuth
    rotmat(3, taz-90., sm1, smt1);
    // Step 2: rotate 1-axis about 2-axis so that 1-axis liesalong t
    rotmat(2, tdp, sm2, smt2);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, sm2, 3, sm1, 3, 0.0, sm3, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, smt1, 3, smt2, 3, 0.0, smt3, 3);
    // in this system, p should lie in 2-3 plane
    cblas_dgemv(CblasColMajor, CblasNoTrans, 3, 3,
                1.0, sm3, 3, p, 1, 0.0, pnew, 1);
    // Step 3: rotate about 1 axis so that 2 axis lies along p
    pang = atan2(pnew[2], pnew[1])/degrad; //(3.14159/180.)
    rotmat(1, -pang, sm1, smt1);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, sm1, 3, sm3, 3, 0.0, rm, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, smt3, 3, smt1, 3, 0.0, rmt, 3);
    // Finally: find moment tensor in space coordinates.
    // This is because Mpt=RM*Mspace*RMT, thus Mspace=RMT*Mpt*RM
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, rmt, 3, mpt, 3, 0.0, work, 3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 3, 3, 3,
                1.0, work, 3, rm, 3, 0.0, mt, 3);
    // ampMax is the absolute maximum amplitude
    *ampMax = MAX(fabs(pev1), MAX(fabs(tev1),  fabs(bev1)));
    return;
}
//============================================================================//
/*!
 * @brief Computes the lower-hemisphere projection `beachball' diagram for a
 *        focal mechanism.  This is based on Cliff Frohlich's bballP and
 *        has been translated to C and does no plotting.
 *
 * @param[in] xc               plot coordinates (e.g. pixels) of center of
 *                             beachball
 * @param[in] yc               plot coordinates (e.g. pixels) of center of
 *                             beachball
 * @param[in] rad              beachball radius (e.g. pixels)
 * @param[in] taz1             azimuth (degrees) of tension eigenvector
 * @param[in] tdp1             dip (degrees) of tension eigenvector
 * @param[in] tevI             eigenvalue corresponding to tension eigenvector
 *                             (e.g. Newton meters)
 * @param[in] baz1             azimuth (degrees) of null eigenvector
 * @param[in] bdp1             dip (degrees) of null eigenvector
 * @param[in] bevI             eigenvalue corresponding to null eigenvector
 *                             (e.g. Newton meters)
 * @param[in] paz1             azimuth (degrees) of pressure eigenvector
 * @param[in] pdp1             dip (degrees) of pressure eigenvector
 * @param[in] pevI             eigenvalue corresponding to pressure eigenvector
 *                             (e.g. Newton meters)
 * @param[in] nNodalLine       number of points to draw in nodal line.
 * @param[in] nFocalSphere     number of points to draw in focal sphere.
 * @param[in] nxPixel          number of pixels in x when drawing polarity
 * @param[in] nyPixel          number of pixels in y when drawing polairty
 * @param[in] lwantNodalLines  if true then compute the nodal lines.
 * @param[in] lwantPTBAxes     if true then compute the PTB axis points.
 * @param[in] lwantPolarity    if true then compute the polarity.
 *
 * @param[out] beachBall       container for the pressure, tension and null axes,
 *                             the nodal lines, and the polarity
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 * @bug this routine fails if rad is not = 1.
 *
 */
void cliffsNodes_beachBallP(const double xc, const double yc, const double rad,
                            const double taz1, const double tdp1,
                            const double tevI,
                            const double baz1, const double bdp1,
                            const double bevI,
                            const double paz1, const double pdp1,
                            const double pevI,
                            const int nNodalLine,
                            const int nFocalSphere,
                            const int nxPixel, const int nyPixel,
                            const bool lwantNodalLines,
                            const bool lwantPTBAxes,
                            const bool lwantPolarity,
                            struct beachBall_struct *beachBall)
{
    double *xw1, *xw2, *yw1, *yw2;
    double rm[9] __attribute__ ((aligned(64)));
    double rmt[9] __attribute__ ((aligned(64)));
    double tm[9] __attribute__ ((aligned(64)));
    double pt[3] __attribute__ ((aligned(64)));
    double ptnew[3] __attribute__ ((aligned(64)));
    double ax1[3] __attribute__ ((aligned(64)));
    double ax2[3] __attribute__ ((aligned(64)));
    double ax3[3] __attribute__ ((aligned(64)));
    double ptLine1[3] __attribute__ ((aligned(64)));
    double ptLine2[3] __attribute__ ((aligned(64)));
    double ptLine1Old[3] __attribute__ ((aligned(64)));
    double ptLine2Old[3] __attribute__ ((aligned(64)));
    double ampMax, az, baz, bdp, bev, cosPsi, cosTheta,
           dip, dp, dpsi, dx, dy, evI, fclvd, fIso, howClose, rad2,
           rb, rp, rpt, rpt2, rt, paz, pdp, pev, plusMin, psi,
           sinsq, sinPsi, sinTheta, str,
           taz, tdp, tev, theta, x, x0, xpt, y, y0, ypt;
    int *pn1, *pn2, ix, ixy, iy, k, nlines, nnl1, nnl2, npts, nseg,
        nx, ny, penDown, psave1, psave2;
    size_t nwork;
    const double sqrt2 = 1.4142135623730951;
    const double sqrt2i = 1.0/sqrt2;
    const double degrad = M_PI/180.0;
    const double degradi = 180.0/M_PI;
    const double third = 1.0/3.0;
    const double tol = 0.95;
    // null out result
    memset(beachBall, 0, sizeof(struct beachBall_struct));
    // Avoid overwriting memory
    taz = taz1;
    tdp = tdp1;
    baz = baz1;
    bdp = bdp1;
    paz = paz1;
    pdp = pdp1;
    // Step 1: set up axes:
    //  - D axis ('dominant' 1-axis - with largest absolute eigenvalue)
    //  - M axis ('minor'    3-axis)
    //  - B axis (           2-axis)
    evI = (tevI + bevI + pevI)*third;
    tev = tevI - evI;
    bev = bevI - evI;
    pev = pevI - evI;
    vecstrd(paz, pdp, ax1);
    vecstrd(taz, tdp, ax3);
    fclvd = fabs(bev/pev);
    fIso = evI/pev;
    if (fabs(tev) > fabs(pev))
    {
        vecstrd(taz, tdp, ax1);
        vecstrd(paz, pdp, ax3);
        fclvd = fabs(bev/tev);
        fIso = evI/tev;
    }
    //if (fabs(evI) > .03)
    //{
    //    write(6,901) evI
    //}
    crossp(ax1, ax3, ax2);
    strdvec(ax2, &baz, &bdp);
    cliffsNodes_mtensor(taz, tdp, paz, pdp, tevI, bevI, pevI,
                        tm, rm, rmt, &ampMax);
    // Step 2: Now: plot upper half nodal line:
    // point on nodal line is represented by vector pt;
    //  pt is defined by two angles
    //  - psi is azimuth with respect to 3-axis (M-axis)
    //    in 2-3 plane (in B-M plane)
    //  - theta is angle with respect to 1 axis.
    // It can be shown (i. e., cliff did it once) that:
    //   sin(theta)**2=(2+2*fIso)/(3+(1-2*f)*cos(2*psi))
    // Additionally, a constraint corresponding to Cliff's Rule 1
    // for existance of P nodal planes is applied
    //psi = 0.;
    //dpsi = 2.;
    //npts  = nNodalLine1; //1 + (int)(360.0/dpsi + 0.5); //nint((360./dpsi))
    //dpsi = dpsi*degrad;
    psave1 =-1;
    psave2 =-1;
    memset(ptLine1Old, 0, 3*sizeof(double));
    memset(ptLine2Old, 0, 3*sizeof(double));
    if (lwantNodalLines && nNodalLine > 1 &&
        !(tevI > 0.0 && bevI > 0.0 && pevI > 0.0) &&
        !(tevI < 0.0 && bevI < 0.0 && pevI < 0.0))
    {
        npts = nNodalLine;
        psi = 0.0;
        dpsi = 360.0/(double) (nNodalLine - 1);
        dpsi = dpsi*degrad;
        xw1 = (double *) calloc((size_t) npts, sizeof(double));
        yw1 = (double *) calloc((size_t) npts, sizeof(double));
        pn1 = (int *) calloc((size_t) npts+1, sizeof(int));
        xw2 = (double *) calloc((size_t) npts, sizeof(double));
        yw2 = (double *) calloc((size_t) npts, sizeof(double));
        pn2 = (int *) calloc((size_t) npts+1, sizeof(int));
        nnl1 = 0;
        nnl2 = 0;
        for (k=0; k<npts; k++)
        {
            psi = k*dpsi;
            theta = 0.0;
            sinsq = (2. + 2.*fIso)/(3. + (1. - 2.*fclvd)*cos(2.*psi));
            if (sinsq > 0 && sinsq <= 1.0)
            {
                theta = asin(sqrt(sinsq));
                cosTheta = cos(theta);
                sinTheta = sin(theta);
                cosPsi = cos(psi);
                sinPsi = sin(psi);
                // upper half nodal line
                ptLine1[0] = cosTheta*ax1[0]
                           + sinTheta*(sinPsi*ax2[0] + cosPsi*ax3[0]);
                ptLine1[1] = cosTheta*ax1[1]
                           + sinTheta*(sinPsi*ax2[1] + cosPsi*ax3[1]);
                ptLine1[2] = cosTheta*ax1[2]
                           + sinTheta*(sinPsi*ax2[2] + cosPsi*ax3[2]);
                // lower half nodal line
                ptLine2[0] =-cosTheta*ax1[0]
                           - sinTheta*(sinPsi*ax2[0] + cosPsi*ax3[0]);
                ptLine2[1] =-cosTheta*ax1[1]
                           - sinTheta*(sinPsi*ax2[1] + cosPsi*ax3[1]);
                ptLine2[2] =-cosTheta*ax1[2]
                           - sinTheta*(sinPsi*ax2[2] + cosPsi*ax3[2]);
                isItNear(ptLine1, &howClose, ptLine1Old);
                // upper half nodal line
                penDown =-1;
                if (ptLine1[2] < 0.0){penDown = 1;}
                if (howClose < tol){penDown =-1;}
                strdv(ptLine1, &str, &dip);
                rpt = rad*sqrt2*sin(0.5*degrad*(90.-dip));
                xpt = xc + rpt*sin(str*degrad);
                ypt = yc + rpt*cos(str*degrad);
                xw1[k] = xpt;
                yw1[k] = ypt;
                if (k == 0){pn1[k] = 3;}
                if (penDown == 1 && psave1 == 1)
                {
                    pn1[k] = 2;
                    //printf("%f %f\n", xpt, ypt);
                }
                if (penDown ==-1){pn1[k] = 3;}
                nnl1 = nnl1 + 1;
                psave1 = penDown;
                // lower half nodal line
                penDown =-1;
                if (ptLine2[2] < 0.0){penDown = 1;}
                if (howClose < tol){penDown =-1;}
                strdv(ptLine2, &str, &dip);
                rpt = rad*sqrt2*sin(0.5*degrad*(90.-dip));
                xpt = xc + rpt*sin(str*degrad);
                ypt = yc + rpt*cos(str*degrad);
                xw2[k] = xpt;
                yw2[k] = ypt;
                if (k == 0){pn2[k] = 3;} // pen down
                if (penDown == 1 && psave2 == 1)
                {
                    pn2[k] = 2;
                    //xptsLine2[k] = xpt;
                    //yptsLine2[k] = ypt;
                    //penLine2[k] = penDwn;
                }
                if (penDown ==-1){pn2[k] = 3;}
                nnl2 = nnl2 + 1;
                psave2 = penDown;
            }
        } // loop on points
        nwork = (size_t) (MAX(2, 2*nnl1));
        beachBall->xNodalLine1 = (double *) calloc(nwork, sizeof(double));
        beachBall->yNodalLine1 = (double *) calloc(nwork, sizeof(double));
        nwork = (size_t) (nnl1 + 1);
        beachBall->nodalLine1Ptr = (int *) calloc(nwork, sizeof(int));
        nwork = (size_t) (MAX(2, 2*nnl2));
        beachBall->xNodalLine2 = (double *) calloc(nwork, sizeof(double));
        beachBall->yNodalLine2 = (double *) calloc(nwork, sizeof(double));
        nwork = (size_t) (nnl2 + 1);
        beachBall->nodalLine2Ptr = (int *) calloc(nwork, sizeof(int));
        // make sure the lines end
        pn1[nnl1] = 3;
        pn2[nnl2] = 3;
        // compute the nodal line 1 line segments
        nlines = 0;
        nseg = 0;
        for (k=0; k<nnl1; k++)
        {
            // line begins
            if (pn1[k] == 3 && pn1[k+1] == 2)
            {
                beachBall->xNodalLine1[2*nlines  ] = xw1[k];
                beachBall->xNodalLine1[2*nlines+1] = xw1[k+1];
                beachBall->yNodalLine1[2*nlines  ] = yw1[k];
                beachBall->yNodalLine1[2*nlines+1] = yw1[k+1];
                nlines = nlines + 1;
            }
            // line continues
            if (pn1[k] == 2 && pn1[k+1] == 2)
            {
                beachBall->xNodalLine1[2*nlines  ] = xw1[k];
                beachBall->xNodalLine1[2*nlines+1] = xw1[k+1];
                beachBall->yNodalLine1[2*nlines  ] = yw1[k];
                beachBall->yNodalLine1[2*nlines+1] = yw1[k+1];
                nlines = nlines + 1;
            }
            // line ends
            if (pn1[k] == 2 && pn1[k+1] == 3)
            {
/*
                beachBall->xNodalLine1[2*nlines  ] = xw1[k];
                beachBall->xNodalLine1[2*nlines+1] = xw1[k+1];
                beachBall->yNodalLine1[2*nlines  ] = yw1[k];
                beachBall->yNodalLine1[2*nlines+1] = yw1[k+1];
                nlines = nlines + 1;
*/
                nseg = nseg + 1;
                beachBall->nodalLine1Ptr[nseg] = nlines;
            }
        }
        beachBall->nNodalLineSegs1 = nseg;
        // compute the nodal line 2 line segments
        nlines = 0;
        nseg = 0;
        for (k=0; k<nnl2; k++)
        {
            // line begins
            if (pn2[k] == 3 && pn2[k+1] == 2)
            {
                beachBall->xNodalLine2[2*nlines  ] = xw2[k];
                beachBall->xNodalLine2[2*nlines+1] = xw2[k+1];
                beachBall->yNodalLine2[2*nlines  ] = yw2[k];
                beachBall->yNodalLine2[2*nlines+1] = yw2[k+1];
                nlines = nlines + 1;
            }
            // line continues
            if (pn2[k] == 2 && pn2[k+1] == 2)
            {
                beachBall->xNodalLine2[2*nlines  ] = xw2[k];
                beachBall->xNodalLine2[2*nlines+1] = xw2[k+1];
                beachBall->yNodalLine2[2*nlines  ] = yw2[k];
                beachBall->yNodalLine2[2*nlines+1] = yw2[k+1];
                nlines = nlines + 1;
            }
            // line ends
            if (pn2[k] == 2 && pn2[k+1] == 3)
            {
                beachBall->xNodalLine2[2*nlines  ] = xw2[k];
                beachBall->xNodalLine2[2*nlines+1] = xw2[k+1];
                beachBall->yNodalLine2[2*nlines  ] = yw2[k];
                beachBall->yNodalLine2[2*nlines+1] = yw2[k+1];
                nlines = nlines + 1;
                nseg = nseg + 1;
                beachBall->nodalLine2Ptr[nseg] = nlines;
            }
        }
        beachBall->nNodalLineSegs2 = nseg;
//printf("%d %d\n", nlines, nseg);
        free(pn1);
        free(pn2);
        free(xw1);
        free(xw2);
        free(yw1);
        free(yw2);
    }
    // Step 4: Plot the, P, T, and B axes
    if (lwantPTBAxes)
    {
        rp = rad*sqrt2*sin(0.5*degrad*(90.-pdp));
        rt = rad*sqrt2*sin(0.5*degrad*(90.-tdp));
        rb = rad*sqrt2*sin(0.5*degrad*(90.-bdp));
        beachBall->xp = xc + rp*sin(paz*degrad);
        beachBall->yp = yc + rp*cos(paz*degrad);
        beachBall->xt = xc + rt*sin(taz*degrad);
        beachBall->yt = yc + rt*cos(taz*degrad);
        beachBall->xb = xc + rb*sin(baz*degrad);
        beachBall->yb = yc + rb*cos(baz*degrad);
    }
    // Step 5: Plot the focal sphere
    if (nFocalSphere > 1)
    {
        beachBall->nFocalSphere = nFocalSphere;
        beachBall->xFocalSphere = (double *)
                                  calloc((size_t) nFocalSphere, sizeof(double));
        beachBall->yFocalSphere = (double *)
                                  calloc((size_t) nFocalSphere, sizeof(double));
        dpsi = 2.0*M_PI/(double) (nFocalSphere - 1);
#ifdef _OPENMP
        #pragma omp simd
#endif
        for (k=0; k<nFocalSphere; k++)
        {
            psi = (double) k*dpsi;
            beachBall->xFocalSphere[k] = xc + rad*sin(psi);
            beachBall->yFocalSphere[k] = yc + rad*cos(psi);
        }
    }
    // Step 6: Plot the + and -
    if (lwantPolarity && nxPixel > 0 && nyPixel > 0)
    {
        // begin in lower left corner
        nx = nxPixel;
        ny = nyPixel;
        dx = 2.0*rad/(double) (nx - 1);
        dy = 2.0*rad/(double) (ny - 1);
        x0 = xc - rad;
        y0 = yc - rad;
        rad2 = rad*rad; // saves a sqrt computation when comparing distances
        xw1 = (double *)calloc((size_t) (nx*ny), sizeof(double));
        yw1 = (double *)calloc((size_t) (nx*ny), sizeof(double));
        pn1 = (int *)calloc((size_t) (nx*ny), sizeof(int));
        nwork = 0;
        // loop on pixel grid
        for (iy=0; iy<ny; iy++)
        {
            for (ix=0; ix<nx; ix++)
            {
                x = x0 + (double) ix*dx;
                y = y0 + (double) iy*dy;
                k = iy*nx + ix;
                pn1[k] = 0;
                xw1[k] = x;
                yw1[k] = y;
                // require (x,y) be in the focal sphere
                rpt2 = pow(x - xc, 2) + pow(y - yc, 2);
                if (rpt2 < rad2)
                {
                    nwork = nwork + 1;
                    rpt = sqrt(rpt2);
                    // compute the azimuth
                    theta = atan2(y - yc, x - xc);  // measure from x axis
                    az = (M_PI_2 - theta)*degradi;  // az is measured from north
                    if (az < 0.0){az = az + 360.0;} // convention is [0,360]
                    // compute the dip
                    dp = (M_PI_2 - 2.0*asin(rpt*sqrt2i))*degradi;
                    vecstrd(az, dp, pt);
                    //cblas_dgemv(CblasColMajor, CblasNoTrans, 3, 3,
                    //            1.0, tm, 3, pt, 1, 0.0, ptnew, 1);
                    gemv3(tm, pt, ptnew);
                    plusMin = dot3(pt, ptnew);
                    pn1[k] =-1;
                    if (plusMin > 100.0*DBL_EPSILON*ampMax){pn1[k] = 1;}
                } // end check on location
            } // loop on x
        } // loop on y
        // Copy the points to be plotted
        if (nwork > 0)
        {
            beachBall->npolarity = (int) nwork;
            beachBall->xPolarity = (double *)
                                   calloc((size_t) nwork, sizeof(double));
            beachBall->yPolarity = (double *)
                                   calloc((size_t) nwork, sizeof(double));
            beachBall->polarity = (int *)
                                  calloc((size_t) nwork, sizeof(double));
            k = 0;
            for (ixy=0; ixy<nx*ny;  ixy++)
            {
                if (pn1[ixy] != 0)
                {
                    beachBall->xPolarity[k] = xw1[ixy];
                    beachBall->yPolarity[k] = yw1[ixy];
                    beachBall->polarity[k] = pn1[ixy];
                    k = k + 1;
                }
            }
        }
        free(pn1);
        free(xw1);
        free(yw1);
    }
    return;
}

void cliffsNodes_free(struct beachBall_struct *beachBall)
{
    if (beachBall->xNodalLine1 != NULL){free(beachBall->xNodalLine1);}
    if (beachBall->yNodalLine1 != NULL){free(beachBall->yNodalLine1);}
    if (beachBall->xNodalLine2 != NULL){free(beachBall->xNodalLine2);}
    if (beachBall->yNodalLine2 != NULL){free(beachBall->yNodalLine2);}
    if (beachBall->xFocalSphere != NULL){free(beachBall->xFocalSphere);}
    if (beachBall->yFocalSphere != NULL){free(beachBall->yFocalSphere);}
    if (beachBall->xPolarity != NULL){free(beachBall->xPolarity);}
    if (beachBall->yPolarity != NULL){free(beachBall->yPolarity);}
    if (beachBall->polarity != NULL){free(beachBall->polarity);}
    if (beachBall->nodalLine1Ptr != NULL){free(beachBall->nodalLine1Ptr);}
    if (beachBall->nodalLine2Ptr != NULL){free(beachBall->nodalLine2Ptr);}
    memset(beachBall, 0, sizeof(struct beachBall_struct));
    return;
}
//============================================================================//
/*!
 * @brief Finds simple rotation matrix R and transpose R^T for rotation
 *        of a vector by angle theta about x, y, or z.  This is based on
 *        Cliff Rohlich's rotmat.
 *
 * @param[in] naxis    rotation axis.  if 1 then x.  if 2 then y.  if 3 then z.
 * @param[in] theta    rotation angle (degrees)
 * @param[out] r       rotation matrix [3 x 3] in column major format.
 * @param[out] rt      tranpsose rotation matrix [3 x 3] in column major format.
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void rotmat(const int naxis, const double theta,
                   double *__restrict__ r,
                   double *__restrict__ rt)
{
    double cosThetar, sinThetar, thetar;
    int i, j;
    const double degrad = M_PI/180.0;
    memset(r, 0, 9*sizeof(double));
    thetar = theta*degrad;
    cosThetar = cos(thetar);
    sinThetar = sin(thetar);
    r[0] = cosThetar;
    r[4] = cosThetar;
    r[8] = cosThetar;
    // naxis = x
    if (naxis == 1)
    {
        r[0] = 1.0;       // (1, 1)
        r[5] = sinThetar; // (3, 2)
        r[7] =-sinThetar; // (2, 3)
    }
    // naxis = y
    else if (naxis == 2)
    {
        r[2] = sinThetar; // (3, 1)
        r[4] = 1.0;       // (2, 2)
        r[6] =-sinThetar; // (1, 3)
    }
    // naxis = z
    else if (naxis == 3)
    {
        r[1] = sinThetar; // (2, 1)
        r[3] =-sinThetar; // (1, 2)
        r[8] = 1.0;       // (3, 3)
    }
    // Compute the transpose
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            rt[3*i+j] = r[3*j+i];
        }
    }
    return;
}
//============================================================================//
/*!
 * @brief Find components of downward pointing unit vector pole from the
 *        strike and dip angles.  This is based on Cliff Frolich's vecstrd.
 *
 * @param[in] str    strike angle of pole (degrees)
 * @param[in] dip    dip angle of pole (degrees)
 *
 * @param[out] r     corresponding unit vector pole (3)
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void vecstrd(const double str, const double dip, double r[3])
{
    double cosdip, cosstr, diprad, sindip, sinstr, strrad;
    const double degrad = M_PI/180.0;
    // convert to radians
    strrad = str*degrad;
    diprad = dip*degrad;
    // compute cosines and sines
    sinstr = sin(strrad);
    cosstr = cos(strrad);
    cosdip = cos(diprad);
    sindip = sin(diprad);
    // compute (x,y,z) from strike and dip for unit radius
    r[0] = sinstr*cosdip;
    r[1] = cosstr*cosdip;
    r[2] =-sindip;
    return;
}
//============================================================================//
static void isItNear(const double *__restrict__ pt,
                     double *howClose,
                     double *__restrict__ ptold)
{
    *howClose = dot3(ptold, pt); //cblas_ddot(3, ptold, 1, pt, 1);
    ptold[0] = pt[0];
    ptold[1] = pt[1];
    ptold[2] = pt[2];
    return;
}
/*!
 * @brief Finds the strike and dip of downward pointing unit vector.  This
 *        is based on Cliff Frolich's strdvec.
 *
 * @param[in,out] r      on input contains the downward pointing unit
 *                       vector.
 *                       if r[2] is > 0 then the unit vector is pointing
 *                       updward.  in this case then on output r3 will
 *                       be negated.
 *
 * @param[out] str       strike angle (degrees) \f$ s.t. \phi \in [0,360] \f$.
 * @param[out] dip       dip angle (degrees)
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void strdvec(double r[3], double *str, double *dip)
{
    double rs;
    const double degradi = 180.0/M_PI;
    if (r[2] > 0.0)
    {
        r[0] =-r[0];
        r[1] =-r[1];
        r[2] =-r[2];
    }
    *str = atan2(r[0], r[1])*degradi;
    if (*str < 0.){*str = *str + 360.;}
    rs = sqrt(r[0]*r[0] + r[1]*r[1]);
    *dip=atan2(-r[2], rs)*degradi;
    return;
}
/*!
 * @brief Cross product of 2 length 3 vectors u and v.  Returns normal vector n.
 *
 * @param[in] u    first vector in cross product u x v [3]
 * @param[in] v    second vector in cross product u x v [3]
 *
 * @param[out] n   normal vector n = u x v [3]
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void crossp(const double *__restrict__ u,
                   const double *__restrict__ v,
                   double *__restrict__ n)
{
    n[0] = u[1]*v[2] - u[2]*v[1];
    n[1] = u[2]*v[0] - u[0]*v[2];
    n[2] = u[0]*v[1] - u[1]*v[0];
    return;
}
/*!
 * @brief Finds strike and dip given components of unit vector.  This is
 *        based on Cliff Frohlich's strdv.
 *
 * @param[in] r      unit vector
 * @param[out] str   strike angle (degrees)
 * @param[out] dip   dip angle (degrees) s.t. dip down is positive
 *
 * @author Ben Baker
 *
 * @copyright BSD
 *
 */
static void strdv(const double r[3], double *str, double *dip)
{
    double rs;
    const double degradi = 180.0/M_PI;
    *str = atan2(r[0], r[1])*degradi;
    rs = sqrt(r[0]*r[0] + r[1]*r[1]);
    *dip = atan2(-r[2], rs)*degradi;
    return;
}

/*!
 * @brief Utility function for computing dot product of length 3 vector
 */
static double dot3(const double *__restrict__ x,
                   const double *__restrict__ y)
{
    double dot;
    dot = x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
    return dot;
}

/*!
 * @brief Utility function for computing y = Ax for 3 x 3 matrix.
 */
static void gemv3(const double *__restrict__ A, const double *__restrict__ x,
                  double *__restrict__ y)
{
    y[0] = A[0]*x[0] + A[3]*x[1] + A[6]*x[2];
    y[1] = A[1]*x[0] + A[4]*x[1] + A[7]*x[2];
    y[2] = A[2]*x[0] + A[5]*x[1] + A[8]*x[2];
    return;
}


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "parmt_polarity.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "ttimes.h"
#include "iscl/array/array.h"
#include "iscl/geodetic/geodetic.h"
#include "iscl/memory/memory.h"

#define LDG 8
#ifndef MAX
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
#endif

static void fillBasis(const double i, const double phi,
                      double *__restrict__ gam,
                      double *__restrict__ phat,
                      double *__restrict__ phihat);
static double computeContraction3x3(const double *__restrict__ a,
                                    const double *__restrict__ M,
                                    const double *__restrict__ b);
static void setM3x3(const int k, double *__restrict__ M);
static int computePadding64f(const int n);
static int performPolaritySearch64f(const int nmt, const int ldm,
                                    const int nobs,
                                    const int blockSize, const int mblock,
                                    const int Mrows, const int Kcols,
                                    const double *__restrict__ Dmat,
                                    const double *__restrict__ G,
                                    const double *__restrict__ Sigma,
                                    const double *__restrict__ mts,
                                    double *__restrict__ phi);

/*!
 * @brief Driver routine for computing the polarity Green's functions from the
 *        ttimes ak135 global travel time table.
 *
 * @param[in] globalComm     Global MPI communicator.
 * @param[in] parms          Contains the polarity modeling parameters.
 * @param[in] data           The SAC data whose header information will define
 *                           the polarity and channel information.
 *
 * @param[out] polarityData  On exit contains the corresponding polarity Green's
 *                           functions that can be used by the grid-search to
 *                           estimate polarities from a given moment tensor.
 *
 * @result 0 indicates success.
 */
int parmt_polarity_computeTTimesGreens(
    const MPI_Comm globalComm,
    const struct parmtPolarityParms_struct parms,
    const struct parmtData_struct data,
    struct polarityData_struct *polarityData)
{
    char kt0[8], kcmpnm[8], stat[8];
    double G6[6], *cmpazs, *cmpincs, *deps, *evlas, *evlos,
           *GxxBuf, *GyyBuf, *GzzBuf, *GxyBuf, *GxzBuf, *GyzBuf,
           *stlas, *stlos, cmpinc, cmpincSAC, cmpaz, stla, stlo;
    int *icomps, *observation, *polarity, *waveType, icomp, ierr, ierrAll,
        iloc, iobs, ipol, it, iwav, jloc, k, kt, myid, nPolarity, nprocs;
    size_t lenos;
    const int master = 0;
    const int nTimeVars = 11;
    const enum sacHeader_enum timeVarNames[11]
       = {SAC_CHAR_KA,
          SAC_CHAR_KT0, SAC_CHAR_KT1, SAC_CHAR_KT2, SAC_CHAR_KT3,
          SAC_CHAR_KT4, SAC_CHAR_KT5, SAC_CHAR_KT6, SAC_CHAR_KT7,
          SAC_CHAR_KT8, SAC_CHAR_KT9};
    //------------------------------------------------------------------------//
    //
    // Initialize
    ierr = 0;
    cmpazs = NULL;
    cmpincs = NULL;
    deps = NULL;
    evlas = NULL;
    evlos = NULL;
    GxxBuf = NULL;
    GyyBuf = NULL;
    GzzBuf = NULL;
    GxyBuf = NULL;
    GxzBuf = NULL;
    GyzBuf = NULL;
    icomps = NULL;
    observation = NULL;
    polarity = NULL;
    stlas = NULL;
    stlos = NULL;
    waveType = NULL;
    memset(polarityData, 0, sizeof(struct polarityData_struct));
    // ttimes uses read-only file-io and i don't want to fix it
    MPI_Comm_rank(globalComm, &myid);
    MPI_Comm_size(globalComm, &nprocs);
    //if (myid != master){goto WAIT;}
    // Verify there is a chance for something to do
    if (data.nobs < 1 || data.nlocs < 1)
    {
        if (data.nobs < 1){fprintf(stderr, "%s: No observations\n", __func__);}
        if (data.nlocs < 1){fprintf(stderr, "%s: No locations\n", __func__);}
        return 0;
    }
    // Extract the depths in the grid-search from the first waveform
    iobs = 0;
    deps = memory_calloc64f(data.nlocs);
    evlas = memory_calloc64f(data.nlocs);
    evlos = memory_calloc64f(data.nlocs);
    for (iloc=0; iloc<data.nlocs; iloc++)
    {
        k = iobs*data.nlocs + iloc;
        ierr  = sacio_getFloatHeader(SAC_FLOAT_EVDP, data.sacGxx[k].header,
                                     &deps[iloc]);
        ierr += sacio_getFloatHeader(SAC_FLOAT_EVLA, data.sacGxx[k].header,
                                     &evlas[iloc]);
        ierr += sacio_getFloatHeader(SAC_FLOAT_EVLO, data.sacGxx[k].header,
                                     &evlos[iloc]);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Unable to get event coordinates on %d\n",
                    __func__, myid);
            goto ERROR;
        }
    }
    // Compute the number of usable polarities (i.e. observations with P picks
    // or S picks with first motions)
    nPolarity = 0;
    if (myid == master)
    {
        polarity = memory_calloc32i(data.nobs);
        waveType = memory_calloc32i(data.nobs);
        icomps = memory_calloc32i(data.nobs);
        stlas = memory_calloc64f(data.nobs);
        stlos = memory_calloc64f(data.nobs);
        cmpazs = memory_calloc64f(data.nobs);
        cmpincs = memory_calloc64f(data.nobs);
        observation = array_set32i(data.nobs, -1, &ierr);
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            // Ensure the essential preliminary information is defined
            ierr  = sacio_getFloatHeader(SAC_FLOAT_CMPINC,
                                         data.data[iobs].header, &cmpincSAC);
            ierr += sacio_getFloatHeader(SAC_FLOAT_CMPAZ,
                                         data.data[iobs].header, &cmpaz);
            ierr += sacio_getFloatHeader(SAC_FLOAT_STLA,
                                         data.data[iobs].header, &stla);
            ierr += sacio_getFloatHeader(SAC_FLOAT_STLO,
                                         data.data[iobs].header, &stlo);
            ierr += sacio_getCharacterHeader(SAC_CHAR_KCMPNM,
                                             data.data[iobs].header, kcmpnm);
            if (ierr != 0)
            {
                fprintf(stderr, "%s: Failed to get header information\n",
                        __func__);
                continue;
            }
            // Figure out the component
            ierr = parmt_utils_getComponent(kcmpnm, cmpincSAC, &icomp);
            if (ierr != 0)
            {
                fprintf(stderr, "%s: Failed to classify component\n", __func__);
            }
            // SAC to SEED convention
            cmpinc = cmpincSAC - 90.0;
            // Figure out the component
/*
            lenos = MAX(1, strlen(kcmpnm));
            icomp = 1;
            if (kcmpnm[lenos-1] == 'Z' || kcmpnm[lenos-1] == '1')
            {
                icomp = 1;
            }
            else if (kcmpnm[lenos-1] == 'N' || kcmpnm[lenos-1] == '2')
            {
                icomp = 2;
            }
            else if (kcmpnm[lenos-1] == 'E' || kcmpnm[lenos-1] == '3')
            {
                icomp = 3;
            }
            else
            {
                fprintf(stderr, "%s: Cannot classify component %s\n", __func__, kcmpnm);
                continue;
            }
*/
            // Get the primary pick
            for (it=0; it<nTimeVars; it++)
            {
                ierr = sacio_getCharacterHeader(timeVarNames[it],
                                                data.data[iobs].header,
                                                kt0);
                if (ierr == 0){break;}
            }
            lenos = strlen(kt0);
            // Ensure the pick is defined
            if (ierr == 0 && lenos > 1)
            {
                iwav = 1;
                if (kt0[0] == 'P' || kt0[0] == 'p')
                {
                    iwav = 1;
                }
                else if (kt0[0] == 'S' || kt0[0] == 's')
                {
                    iwav = 2;
                }
                else
                {
                    // surface waves will commonly be processed and not
                    // have polarities
                    if (kt0[0] == 'R' || kt0[0] == 'r' ||
                        kt0[0] == 'L' || kt0[0] == 'l')
                    {
                        continue;
                    }
                    fprintf(stderr, "%s: t0 phase is not a P or S phase %s\n",
                            __func__, kt0);
                    continue;
                }
                if (kt0[lenos-1] == '+')
                {
                    ipol = 1;
                }
                else if (kt0[lenos-1] == '-')
                {
                    ipol =-1;
                }
                else
                {
                    fprintf(stderr, "%s: could not classify polarity %s\n",
                            __func__, kt0);
                    continue;
                }
                // let user know something happened
                sacio_getCharacterHeader(SAC_CHAR_KSTNM,
                                         data.data[iobs].header, stat);
                fprintf(stdout, "%s: Polarity for %s is %d\n",
                        __func__, stat, ipol);
//printf("%f %f %f %f %d %d %d %d\n", stla, stlo, cmpinc, cmpaz, icomp, iobs, iwav, ipol);
                // save information for modeling
                stlas[nPolarity] = stla;
                stlos[nPolarity] = stlo;
                cmpincs[nPolarity] = cmpinc;
                cmpazs[nPolarity] = cmpaz;
                icomps[nPolarity] = icomp;
                observation[nPolarity] = iobs;
                waveType[nPolarity] = iwav;
                polarity[nPolarity] = ipol;
                nPolarity = nPolarity + 1;
                //printf("%d %d %f %f\n", iwav, ipol, cmpinc, cmpaz);
           } // End basic header info check
        } // Loop on observations
    } // End check on myid == master
    // Tell all other processes about the forward modeling information
    MPI_Bcast(&nPolarity, 1, MPI_INT, master, globalComm);
    if (nPolarity == 0)
    {
        if (myid == master)
        {
            fprintf(stdout, "%s: There are no polarities\n", __func__);
        }
        ierr = 0;
        goto ERROR;
    }
    if (myid != master)
    {
        stlas = memory_calloc64f(nPolarity);
        stlos = memory_calloc64f(nPolarity);
        cmpincs = memory_calloc64f(nPolarity);
        cmpazs = memory_calloc64f(nPolarity);
        icomps = memory_calloc32i(nPolarity);
        observation = memory_calloc32i(nPolarity);
        waveType = memory_calloc32i(nPolarity);
        polarity = memory_calloc32i(nPolarity);
    }
    else
    {
        fprintf(stdout, "%s: Warning - i'm setting wts to unity for now\n",
                __func__);
    }
    MPI_Bcast(stlas,   nPolarity, MPI_DOUBLE, master, globalComm);
    MPI_Bcast(stlos,   nPolarity, MPI_DOUBLE, master, globalComm);
    MPI_Bcast(cmpincs, nPolarity, MPI_DOUBLE, master, globalComm);
    MPI_Bcast(cmpazs,  nPolarity, MPI_DOUBLE, master, globalComm);
    MPI_Bcast(icomps,      nPolarity, MPI_INT, master, globalComm);
    MPI_Bcast(observation, nPolarity, MPI_INT, master, globalComm);
    MPI_Bcast(waveType,    nPolarity, MPI_INT, master, globalComm);
    MPI_Bcast(polarity,    nPolarity, MPI_INT, master, globalComm);
    // Set space
    ierrAll = 0;
    GxxBuf = memory_calloc64f(data.nlocs*nPolarity);
    GyyBuf = memory_calloc64f(data.nlocs*nPolarity);
    GzzBuf = memory_calloc64f(data.nlocs*nPolarity);
    GxyBuf = memory_calloc64f(data.nlocs*nPolarity);
    GxzBuf = memory_calloc64f(data.nlocs*nPolarity);
    GyzBuf = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gxx = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gyy = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gzz = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gxy = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gxz = memory_calloc64f(data.nlocs*nPolarity);
    polarityData->Gyz = memory_calloc64f(data.nlocs*nPolarity);
    // Set the data
    polarityData->polarity = memory_calloc64f(nPolarity);
    for (ipol=0; ipol<nPolarity; ipol++)
    {
        polarityData->polarity[ipol] = (double) polarity[ipol];
    }
    // Set the weights
    // TODO fix me
    polarityData->wts = array_set64f(nPolarity, 1.0, &ierr);
    // Compute the forward modeling matrix columns
    for (jloc=0; jloc<data.nlocs; jloc=jloc+nprocs)
    {
        iloc = jloc + myid;
        if (iloc >= data.nlocs){continue;}
        // Compute the forward modeling matrix for this observation group
        for (ipol=0; ipol<nPolarity; ipol++)
        {
            kt = iloc*nPolarity + ipol;
            ierr = parmt_polarity_computeGreensRowFromTtimes(
                                waveType[ipol], icomps[ipol],
                                evlas[iloc], evlos[iloc], deps[iloc],
                                stlas[ipol], stlos[ipol],
                                cmpincs[ipol], cmpazs[ipol],
                                parms.ttimesTablesDir,
                                parms.ttimesModel,
                                G6);
            if (ierr != 0)
            {
                fprintf(stderr,
                        "%s: Error computing polarities %d %d on PID %d\n",
                       __func__, iloc, ipol, myid);
                ierrAll = ierrAll + 1;
                continue;
            }
            // Save it
            GxxBuf[kt] = G6[0];
            GyyBuf[kt] = G6[1];
            GzzBuf[kt] = G6[2];
            GxyBuf[kt] = G6[3];
            GxzBuf[kt] = G6[4];
            GyzBuf[kt] = G6[5];
            //printf("%f %f %f %e %e %e %e %e %e\n", stlas[ipol], stlos[ipol], deps[iloc], G6[0], G6[1], G6[2], G6[3], G6[4], G6[5]);
        } // Loop on the polarities
        //printf("\n");
    } // Loop on the locations
    ierr = ierrAll;
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error computing polarities from ttimes\n",
                __func__);
        ierr = 1;
        goto ERROR;
    }
    polarityData->nPolarity = nPolarity;
    polarityData->nlocs = data.nlocs;
    polarityData->nobs = data.nobs;
    polarityData->obsMap = array_copy32i(nPolarity, observation, &ierr);
    MPI_Allreduce(GxxBuf, polarityData->Gxx, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    MPI_Allreduce(GyyBuf, polarityData->Gyy, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    MPI_Allreduce(GzzBuf, polarityData->Gzz, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    MPI_Allreduce(GxyBuf, polarityData->Gxy, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    MPI_Allreduce(GxzBuf, polarityData->Gxz, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    MPI_Allreduce(GyzBuf, polarityData->Gyz, data.nlocs*nPolarity,
                  MPI_DOUBLE, MPI_SUM, globalComm);
    // Finally assemble Green's functions into modeling matrix
//WAIT:;
ERROR:;
    MPI_Barrier(globalComm);
    // free memory
    memory_free64f(&GxxBuf);
    memory_free64f(&GyyBuf);
    memory_free64f(&GzzBuf);
    memory_free64f(&GxyBuf);
    memory_free64f(&GxzBuf);
    memory_free64f(&GyzBuf);
    memory_free64f(&deps);
    memory_free64f(&evlas);
    memory_free64f(&evlos);
    memory_free64f(&stlas);
    memory_free64f(&stlos);
    memory_free64f(&cmpincs);
    memory_free64f(&cmpazs);
    memory_free32i(&observation);
    memory_free32i(&polarity);
    memory_free32i(&icomps);
    memory_free32i(&waveType);
    return ierr;
}
//============================================================================//
/*!
 * @brief Computes the row of the Green's functions matrix for modeling
 *        a polarity from the ttimes model.
 *
 * @param[in] data      Contains the pick type and channel information.
 * @param[in] wavetype  P_WAVE (1) indicates a P-wave.
 * @param[in] wavetype  S_WAVE (2) indicates an S-wave.
 * @param[in] evdp      Depth of the event in kilometers.
 * @param[in] dirnm     Directory containing the iasp-tau binary files.
 * @param[in] model     Name of the model, e.g., ak135.
 *
 * @param[out] G        On exit contains the Green's functions so that for
 *                      the i'th row of G, \f$ G_i \cdot \textbf{m} \f$
 *                      computes an estimate of polarity.  This is packed
 *                      in order
 *                      \f$
 *                        \{m_{xx}, m_{yy}, m_{zz}, m_{xy}, m_{xz}, m_{yz} \}
 *                      \f$
 *                      and must have a dimension of at least [6].
 *
 * @result 0 indicates success.
 *
 * @copyright ISTI distributed under the Apache 2 license.
 *
 */
int parmt_polarity_computeGreensRowFromData(const struct sacData_struct data,
                                            const int wavetype,
                                            const double evdp,
                                            const char *dirnm,
                                            const char *model,
                                            double *__restrict__ G)
{
    char kcmpnm[16];
    double cmpaz, cmpinc, cmpincSAC, evla, evlo, stla, stlo;
    int icomp, ierr;
    size_t lenos;
    memset(kcmpnm, 16, 16*sizeof(char));
    ierr  = sacio_getFloatHeader(SAC_FLOAT_CMPINC, data.header, &cmpincSAC);
    ierr += sacio_getFloatHeader(SAC_FLOAT_CMPAZ,  data.header, &cmpaz);
    ierr += sacio_getFloatHeader(SAC_FLOAT_EVLA,   data.header, &evla);
    ierr += sacio_getFloatHeader(SAC_FLOAT_EVLO,   data.header, &evlo);
    ierr += sacio_getFloatHeader(SAC_FLOAT_STLA,   data.header, &stla);
    ierr += sacio_getFloatHeader(SAC_FLOAT_STLO,   data.header, &stlo);
    ierr += sacio_getCharacterHeader(SAC_CHAR_KCMPNM, data.header, kcmpnm);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to get header information\n", __func__);
        return -1;
    }
    // Figure out the component
    ierr = parmt_utils_getComponent(kcmpnm, cmpincSAC, &icomp);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to classify component\n", __func__);
        return -1;
    }
    // SAC to SEED convention
    cmpinc = cmpincSAC - 90.0;
    // Figure out the component
    lenos = MAX(1, strlen(kcmpnm));
    icomp = 1;
    if (kcmpnm[lenos-1] == 'Z' || kcmpnm[lenos-1] == '1')
    {
        icomp = 1;
    }
    else if (kcmpnm[lenos-1] == 'N' || kcmpnm[lenos-1] == '2')
    {
        icomp = 2;
    }
    else if (kcmpnm[lenos-1] == 'E' || kcmpnm[lenos-1] == '3')
    {
        icomp = 3;
    }
    else
    {
        fprintf(stderr, "%s: Can't classify component %s\n", __func__, kcmpnm);
        return -1;
    }
    ierr = parmt_polarity_computeGreensRowFromTtimes(wavetype, icomp,
                                                     evla, evlo, evdp,
                                                     stla, stlo,
                                                     cmpinc, cmpaz,
                                                     dirnm, model,
                                                     G);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to compute G\n", __func__);
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Computes a row of the Green's functions matrix from ttimes
 *
 * @param[in] wavetype   If 1 then this is a P wave.
 * @param[in] wavetype   If 2 then this is an S wave.
 * @param[in] icomp      If 1 then this is the vertical channel.
 * @param[in] icomp      If 2 then this is the north (1 or 2) channel.
 * @param[in] icomp      If 3 then this is the east (2 or 3) channel.
 * @param[in] evla       Event latitude (degrees).
 * @param[in] evlo       Event longitude (degrees).
 * @param[in] evdp       Event depth (km).
 * @param[in] stla       Station latitude (degrees).
 * @param[in] stlo       Station longitude (degrees).
 * @param[in] cmpaz      Component azimuth (0 north, +90 east).
 * @param[in] cmpinc     Component inclination (-90 up, 0 east/north, +90 down).
 * @param[in] dirnm      Directory containing the ttimes precomputed binary
 *                       files.  If NULL then the default as dicated by the
 *                       ttimes configuration will be used.
 * @param[in] model      Model name (e.g., ak135 or iasp91)
 *
 * @param[out] G         Row of matrix s.t. G*m produces estimates the polarity
 *                       at the station.  Here m is packed
 *                       \f$ \{m_{xx}, m_{yy}, m_{zz},
 *                             m_{xy}, m_{xz}, m_{yz} \} \f$
 *                       This must have dimension of at least 6.
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_polarity_computeGreensRowFromTtimes(
    const int wavetype, const int icomp,
    const double evla, const double evlo, const double evdp,
    const double stla, const double stlo,
    const double cmpinc, const double cmpaz,
    const char *dirnm, const char *model,
    double *__restrict__ G)
{
    double aoiRec, az, azSrc, baz, bazRec, dist, delta, toaSrc;
    struct ttimesTravelTime_struct ttime;
    int ierr;
    ierr = 0;
    if (G == NULL)
    {
        fprintf(stderr, "%s: Error G is NULL\n", __func__);
        return -1;
    }
    memset(G, 0, 6*sizeof(double));
    if (evdp < 0.0 || evdp > ttimes_getMaxDepth())
    {
        fprintf(stderr, "%s: Error depth must be between [0,%f]\n", __func__,
                ttimes_getMaxDepth());
        return -1;
    }
    memset(&ttime, 0, sizeof(struct ttimesTravelTime_struct));
    geodetic_gps2distanceAzimuth(evla, evlo, stla, stlo,
                                 &dist, &delta, &az, &baz);
    if (wavetype == 1)
    {
        ierr = ttimes_getFirstPPhase(delta, evdp, dirnm, model, &ttime);
    }
    else if (wavetype == 2)
    {
        ierr = ttimes_getFirstSPhase(delta, evdp, dirnm, model, &ttime);
    }
    else
    {
        fprintf(stderr, "%s: Invalid phase type - must be 1 (P) or 2 (S)\n",
                __func__);
        return -1;
    }
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error computing theoretical traveltime info\n",
                 __func__);
        return -1;
    }
    // Compute the column in the Green's function matrix
    toaSrc = ttime.toang;
    azSrc  = az;
    bazRec = baz;
    aoiRec = ttime.aoi;
    //printf("%f %f %f %f %f %f %f\n", delta, stla, stlo, toaSrc, azSrc, aoiRec, bazRec);
    ierr = parmt_polarity_computeGreensMatrixRow(wavetype, icomp,
                                                 azSrc, toaSrc, bazRec, aoiRec,
                                                 cmpinc, cmpaz, G);

    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to compute polarity for row\n", __func__);
        memset(G, 0, 6*sizeof(double));
        return -1;
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Computes a column for the Green's function s.t. G*m produces an
 *        estimate of polarity on the icomp'th component measured in the
 *        far-field.  The moment tensor which would be applied to this row
 *        is packed:
 *          \f$ \{m_{xx}, m_{yy}, m_{zz}, m_{xy}, m_{xz}, m_{yz} \} \f$
 *        with convention North, East, Down (e.g., Jost and Herrmann).
 *        For more see Quantitative Seismology - Aki and Richards 2002,
 *        Eqn 4.96 on pg 111 and Source Mechanisms of Earthquakes: Theory
 *        and Practice - Udias et al. pg 100..
 *
 * @param[in] wavetype   =1 -> P wave
 *                       =2 -> S wave
 * @param[in] icomp      receiver component of motion:
 *                         =1 -> Vertical channel
 *                         =2 -> 1 or North channel
 *                         =3 -> 2 or East channel
 * @param[in] azSrc      source to receiver azimuth is measured positive from
 *                       north (degrees)
 * @param[in] toaSrc     take-off angle (measured positive from x3 where x3
 *                       points down) (degrees)
 * @param[in] bazRec     receiver to source back azimuth measured positive
 *                       from north (degrees)
 * @param[in] aoiRec     angle of incidence at receiver
 * @param[in] cmpaz      component azimuth (0 north, +90 east)
 * @param[in] cmpinc     component inclinantion (-90 up, 0 east/north, +90 down)
 *
 * @param[out] G         row of matrix s.t. G*m produces estimates the polarity
 *                       at the station.  Here m is packed
 *                       \f$ \{m_{xx}, m_{yy}, m_{zz},
 *                             m_{xy}, m_{xz}, m_{yz} \} \f$
 *
 * @bugs i've really only looked at the vertical p-teleseismic case - no idea
 *        about other phases, wavetypes
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_polarity_computeGreensMatrixRow(const int wavetype,
                                          const int icomp,
                                          const double azSrc,
                                          const double toaSrc,
                                          const double bazRec,
                                          const double aoiRec,
                                          const double cmpinc,
                                          const double cmpaz,
                                          double *__restrict__ G)
{
    double M[9], G63[18], up[6], ush[6], usv[6],
           gam[3], lhat[3], phat[3], phihat[3],
           cosba, cos_cmpaz, cost_rec,
           r11, r12, r13, r21, r22, r23, r31, r32, r33,
           sinba, sin_cmpaz,
           sint_rec,
           t1, t2, t3, theta, xsign, u1, u2, ue, un, uz;
    int k;
    const double pi180 = M_PI/180.0;
    const bool lrot = true;//false;
    const int P_WAVE = 1;
    const int S_WAVE = 2;
    //------------------------------------------------------------------------//
    //
    // Error check
    if (icomp < 1 || icomp > 3)
    {
        fprintf(stderr, "%s: Invalid component\n", __func__);
        return -1;
    }
    if (wavetype < P_WAVE || wavetype > S_WAVE)
    {
        fprintf(stderr, "%s: Invalid wavetype\n", __func__);
        return -1;
    }
    // Fill the basis at the source (Aki and Richards Eqn 4.88)
    fillBasis(toaSrc, azSrc, gam, phat, phihat);
    // Compute the contractions for u_p, u_sv, and u_sh (A&R Eqn 4.96)
    // for all 6 individual moment tensor terms
    for (k=0; k<6; k++)
    {
        // Set the moment tensor with the k'th mt term 1 others zero
        setM3x3(k, M);
        // Compute contraction
        up[k]  = computeContraction3x3(gam,    M, gam);
        usv[k] = computeContraction3x3(phat,   M, gam);
        ush[k] = computeContraction3x3(phihat, M, gam);
    }
    // Fill the basis at the receiver - notice the basis uses the forward
    // azimuth from the receiver to the source
    fillBasis(aoiRec, (bazRec + 180.0), lhat, phat, phihat);
    // Compute geometric factors at receiver
//printf("%f %f %f %f %f\n", toaSrc, azSrc, bazRec, aoiRec, cmpaz);
    theta = aoiRec*pi180;
    cosba     = cos(bazRec*pi180);
    cost_rec  = cos(theta);
    sinba     = sin(bazRec*pi180);
    sint_rec  = sin(theta);
    cos_cmpaz = cos(cmpaz*pi180);
    sin_cmpaz = sin(cmpaz*pi180);
    // Set 3D rotation matrix
    r11 = cost_rec;
    r21 = sint_rec;
    r31 = 0.0;
    r12 =-sint_rec*sinba;
    r22 = cost_rec*sinba;
    r32 =         -cosba;
    r13 =-sint_rec*cosba;
    r23 = cost_rec*cosba;
    r33 =          sinba;
    // Flip sign for receivers that acquire positive down
    xsign = 1.0;
    if (fabs(cmpinc - 90.0) < 1.e-4){xsign =-1.0;}
    // Compute 6 x 3 subforward modeling matrix with row for up (1 channel)
    if (wavetype == P_WAVE)
    {
        // Loop on mts terms
        for (k=0; k<6; k++)
        {
            // Extract the (north, east, down) component
            t3 =-up[k]*lhat[0]; // z-down -> z-up
            //t3 = up[k]*lhat[0]; // z-up
            t2 = up[k]*lhat[1]; // east
            t1 = up[k]*lhat[2]; // north
            // Not sure if i have to rotate LQT -> ZNE
            if (lrot)
            {
                uz = t1*r11 + t2*r21 + t3*r31; // Z
                ue = t1*r12 + t2*r22 + t3*r32; // E
                un = t1*r13 + t2*r23 + t3*r33; // N
            }
            else
            {
                uz = t3;
                ue = t2;
                un = t1;
            }
            // Rotate into (1, 2)
            u1 = un*cos_cmpaz + ue*sin_cmpaz;
            u2 =-un*sin_cmpaz + ue*cos_cmpaz;
            // Finish
            G63[k*3+0] = xsign*uz;
            //G63[k*3+0] =-xsign*uz;
            G63[k*3+1] = u1;
            G63[k*3+2] = u2;
        }
    }
    // SH wave
    else
    {
        // Loop on mts terms
        for (k=0; k<6; k++)
        {
            // Extract the (north, east, down) component
            t3 =-usv[k]*phat[0] - ush[k]*phihat[0]; // z-down -> z-up
            t2 = usv[k]*phat[1] + ush[k]*phihat[1]; // east
            t1 = usv[k]*phat[2] + ush[k]*phihat[2]; // north
            // Not sure if i have to rotate LQT -> ZNE
            if (lrot)
            {
                uz = t1*r11 + t2*r21 + t3*r31; // Z
                ue = t1*r12 + t2*r22 + t3*r32; // E
                un = t1*r13 + t2*r23 + t3*r33; // N
            }
            else
            {
                uz = t3;
                ue = t2;
                un = t1;
            }
            // Rotate into (1, 2)
            u1 = un*cos_cmpaz + ue*sin_cmpaz;
            u2 =-un*sin_cmpaz + ue*cos_cmpaz;
            // Finish
            G63[k*3+0] = xsign*uz;
            G63[k*3+1] = u1;
            G63[k*3+2] = u2;
        }
    }
    // Copy the result - vertical
    if (icomp == 1)
    {
        for (k=0; k<6; k++)
        {
            G[k] = G63[k*3+0];
        }
    }
    // 1 component
    else if (icomp == 2)
    {
        for (k=0; k<6; k++)
        {
            G[k] = G63[k*3+1];
        }
    }
    // 2 component
    else if (icomp == 3)
    {
        for (k=0; k<6; k++)
        {
            G[k] = G63[k*3+2];
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Fills in the basis vectors - Aki and Richards pg 108 Eqn 4.88.
 *
 * @param[in] i       Take-off angle (degrees).
 * @param[in] phi     Azimuth angle (degrees).
 * @param[in] gam     P-wave direction.  This has dimension [3].
 * @param[in] phat    SV-wave direction.  This has dimension [3].
 * @param[in] phihat  SH-wave direction.  This has dimension [3].
 *
 * @copyright ISTI distributed under the Apache 2 license.
 *
 */
static void fillBasis(const double i, const double phi,
                      double *__restrict__ gam,
                      double *__restrict__ phat,
                      double *__restrict__ phihat)
{
    double cosi, sini, cosp, sinp;
    const double pi180 = M_PI/180.0;
    cosi = cos(i*pi180);
    sini = sin(i*pi180);
    cosp = cos(phi*pi180);
    sinp = sin(phi*pi180);
    gam[0] = sini*cosp;
    gam[1] = sini*sinp;
    gam[2] = cosi;

    phat[0] = cosi*cosp;
    phat[1] = cosi*sinp;
    phat[2] =-sini;

    phihat[0] =-sinp;
    phihat[1] = cosp;
    phihat[2] = 0.0;
}
//============================================================================//
/*!
 * @brief Sets the moment tensor matrix where the k'th moment tensor
 *        term is 1 and others are zero.  Here moment tensor terms
 *        are counted {0,1,2,3,4,5} = {xx,yy,zz,xy,xz,yz}
 *
 * @param[in] k    Moment tensor index.  This is in the range [0,5] and follows
 *                 the mapping {0,1,2,3,4,5} = {xx,yy,zz,xy,xz,yz}.
 *
 * @param[out] M   The 3x3 NED moment tensor. This is an array of dimension [9].
 *
 */
static void setM3x3(const int k, double *__restrict__ M)
{
    memset(M, 0, 9*sizeof(double));
    // mxx (fill mtt)
    if (k == 0)
    {
        M[4] = 1.0; //M[0][0] = 1.0;
    }
    // myy (fill mpp)
    else if (k == 1)
    {
        M[8] = 1.0; //M[1][1] = 1.0;
    }
    // mzz (fill mrr)
    else if (k == 2)
    {
        M[0] = 1.0; //M[2][2] = 1.0;
    }
    // mxy and myz (fill mtp)
    else if (k == 3)
    {
        M[5] =-1.0; //M[0][1] = 1.0;
        M[7] =-1.0; //M[1][0] = 1.0;
    }
    // mxz and mzx (fill mrp)
    else if (k == 4)
    {
        M[1] = 1.0; //M[0][2] = 1.0;
        M[3] = 1.0; //M[2][0] = 1.0;
    }
    // myz and mzy (fill mrp)
    else
    {
       M[2] =-1.0; //M[1][2] = 1.0;
       M[6] =-1.0; //M[2][1] = 1.0;
    }
/*
    // mxx
    if (k == 0)
    {
        M[0] = 1.0; //M[0][0] = 1.0;
    }
    // myy
    else if (k == 1)
    {
        M[4] = 1.0; //M[1][1] = 1.0;
    }
    // mzz
    else if (k == 2)
    {
        M[8] = 1.0; //M[2][2] = 1.0;
    }
    // mxy and myz
    else if (k == 3)
    {
        M[1] = 1.0; //M[0][1] = 1.0;
        M[3] = 1.0; //M[1][0] = 1.0;
    }
    // mxz and mzx
    else if (k == 4)
    {
        M[2] = 1.0; //M[0][2] = 1.0;
        M[6] = 1.0; //M[2][0] = 1.0;
    }
    // myz and mzy
    else
    {
       M[5] = 1.0; //M[1][2] = 1.0;
       M[7] = 1.0; //M[2][1] = 1.0;
    }
*/
    return;
}
//============================================================================//
/*!
 * @brief Computes the contraction in Aki and Richards Eqn 4.96
 *        for the given basis vectors and moment tensor.
 */
static double computeContraction3x3(const double *__restrict__ a,
                                    const double *__restrict__ M,
                                    const double *__restrict__ b)
{
    double res;
    int i, ij, j;
    res = 0.0;
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            ij = 3*i + j;
            res = res + a[i]*M[ij]*b[j];
        }
    }
    return res;
}
//============================================================================//
/*!
 * @brief Tabulates the objective function over the earthquake locations.
 *
 * @param[in] locComm       Location MPI communicator.
 * @param[in] blockSize     Block-size for performing matrix-matrix
 *                          multiplications.
 * @param[in] polarityData  Contains the polarity data and the Green's functions.
 * @param[in] mtloc         Contains the local moment tensors in this rank's
 *                          grid search.
 *
 * @param[out] phi          The objective function tabulated for all locations
 *                          and moment tensors.  This is only accessed by
 *                          master process in locComm; and for this case has
 *                          dimension [nlocs x mtloc.nmtAll] with leading
 *                          dimension mtloc.nmtAll. Otherwise, this can be NULL.
 *
 * @result 0 indicates success.
 *
 * @copyright ISTI distributed under the Apache 2 license.
 *
 */
int polarity_performLocationSearch64f(const MPI_Comm locComm,
                                      const int blockSize,
                                      struct polarityData_struct polarityData,
                                      struct localMT_struct mtloc,
                                      double *__restrict__ phi)
{
    double *Dmat, *G, *Sigma, *phiLoc, *phiWork;
    int icol, ierr, ierrAll, iloc, ipol, jloc, jndx,
        kt, mylocID, nlocProcs, mblock, pad;
    const int nPolarity = polarityData.nPolarity;
    const int Mrows = nPolarity;
    const int Kcols = 6;
    const int nlocs = polarityData.nlocs;
    const int ldm = mtloc.ldm;
    const int nmt = mtloc.nmt;
    const int master = 0;
    // Initialize
    ierr = 0;
    ierrAll = 0;
    MPI_Comm_size(locComm, &nlocProcs);
    MPI_Comm_rank(locComm, &mylocID);
    // Compute padding for 64 bit data alignment in observations and mt blocks
    pad = computePadding64f(blockSize);
    mblock = blockSize + pad;
    // Set space
    phiLoc = memory_calloc64f(nmt);
    G = memory_calloc64f(nPolarity*LDG);
    Dmat = memory_calloc64f(nPolarity*mblock);
    Sigma = array_set64f(nPolarity, 1.0, &ierr); // Default to identity
    phiWork = memory_calloc64f(nlocs*mtloc.nmtAll);
    // Set the row major data matrix where each row is an observation
    for (ipol=0; ipol<nPolarity; ipol++)
    {
        for (icol=0; icol<mblock; icol++)
        {
            Dmat[ipol*mblock+icol] = polarityData.polarity[ipol];
            if (icol == 1){printf("%f\n", polarityData.polarity[ipol]);}
        }
    }
    // Input weights are given - copy them to diagonal weight matrix
    if (polarityData.wts != NULL)
    {
        ierr = array_copy64f_work(nPolarity, polarityData.wts, Sigma);
    }
    // Loop on the source locations in the grid search
    for (jloc=0; jloc<nlocs; jloc=jloc+nlocProcs)
    {
        iloc = jloc + mylocID;
        kt = iloc*nPolarity;
        // Assemble the row major Green's functions matrix
        array_zeros64f_work(LDG*nPolarity, G);
        #pragma omp simd
        for (ipol=0; ipol<nPolarity; ipol++)
        {
            G[LDG*ipol+0] = polarityData.Gxx[kt+ipol];
            G[LDG*ipol+1] = polarityData.Gyy[kt+ipol];
            G[LDG*ipol+2] = polarityData.Gzz[kt+ipol];
            G[LDG*ipol+3] = polarityData.Gxy[kt+ipol];
            G[LDG*ipol+4] = polarityData.Gxz[kt+ipol];
            G[LDG*ipol+5] = polarityData.Gyz[kt+ipol];
        }
        // Perform the polarity search for all mt's at this location
        ierr = performPolaritySearch64f(nmt, ldm,
                                        nPolarity,
                                        blockSize, mblock,
                                        Mrows, Kcols,
                                        Dmat, G, Sigma, mtloc.mts,
                                        phiLoc);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error in mt polarity search\n", __func__);
            ierrAll = ierrAll + 1;
        }
        // Gather the moment tensor search results onto the master
        jndx = 0;
        if (mtloc.myid == master){jndx = iloc*mtloc.nmtAll;}
        if (mtloc.commSize == 1)
        {
            array_copy64f_work(mtloc.nmt, phiLoc, &phiWork[jndx]);
        }
        else
        {
            MPI_Gatherv(phiLoc, mtloc.nmt, MPI_DOUBLE,
                        &phiWork[jndx], mtloc.nmtProc, mtloc.offset,
                        MPI_DOUBLE, master, mtloc.comm);
        }
        fprintf(stdout, "%s: mylocID: %d min: %f max: %f\n",
                __func__, mylocID, array_min64f(nmt, phiLoc, &ierr),
                array_max64f(nmt, phiLoc, &ierr));

    }
    // Reduce the search onto the master
    //if (linLocComm && mylocID == master)
    {
         MPI_Reduce(phiWork, phi, nlocs*mtloc.nmtAll, MPI_DOUBLE,
                    MPI_SUM, master, locComm);
    }
    // Free memory
    memory_free64f(&phiWork);
    memory_free64f(&phiLoc);
    memory_free64f(&G);
    memory_free64f(&Dmat);
    memory_free64f(&Sigma);
    return 0;
}
//============================================================================//
/*!
 * @brief Performs the polarity search for a host of moment tensors.
 *
 * @param[in] nmt        Number of moment tensors.
 * @param[in] ldm        Leading dimension of moment tensor matrix.  This must
 *                       be at least 6.
 * @param[in] nPolarity  The number of polarities (observations).
 * @param[in] blockSize  Controls the number of forward problems GM to compute
 *                       simultaneously.
 * @param[in] mblock     The largest block size.  This serves as a leading
 *                       dimension.
 * @param[in] Mrows      Number of rows in Green's functions matrix.
 * @param[in] Kcols      Number of columns in Green's functions matrix.  This
 *                       should be 6.
 * @param[in] Dmat       The polarity data replicated to the max block size.
 *                       This has dimension [nPolarity x mblock] with leading
 *                       dimension mblock.
 * @param[in] G          Forward modeling matrix.  This has dimension
 *                       [nPolarity x LDG] with leading dimension LDG.
 * @param[in] Sigma      These are the data weights and has dimension
 *                       [nPolarity].
 * @param[in] mts        Matrix of moment tensors and has dimension [nmt x ldm].
 *                       with leading dimension ldm.
 *
 * @param[out] phi       The objective function (variance reduction) tabulated
 *                       for all the moment tensors and observations.  This has
 *                       dimension [nmt].
 *
 * @result 0 indicates success.
 *
 * @copyright ISTI distributed under the Apache 2 license.
 *
 */
static int performPolaritySearch64f(const int nmt, const int ldm,
                                    const int nPolarity,
                                    const int blockSize, const int mblock,
                                    const int Mrows, const int Kcols,
                                    const double *__restrict__ Dmat,
                                    const double *__restrict__ G,
                                    const double *__restrict__ Sigma,
                                    const double *__restrict__ mts,
                                    double *__restrict__ phi)
{
    double *U, *Usign, *res2, *sqrtSigmaWt, res, traceSigma;
    int i, ic, ierr, ipol, jmt, kmt, nmtBlocks, Ncols;
    const double one = 1.0;
    const double zero = 0.0;
    ierr = 0;
    nmtBlocks = (int) ((double) (nmt)/(double) (blockSize) + 1.0);
    if (nmtBlocks*blockSize < nmt)
    {
        fprintf(stderr, "%s: Internal error - all mts wont be visited\n",
                __func__);
        return -1;
    }
#ifdef __INTEL_COMPILER
    __assume_aligned(G, 64);
    __assume_aligned(Dmat, 64);
#endif
    // The VR in Chiang 2016 is:
    //   (1 - \frac{ \sum_i w_i (Pol_{i,obs} - Pol_{i,est})^2 }
    //             { \sum_i w_i (Pol_{i,obs}^2) } )
    // Because Pol_{obs} and Pol_{est} take values of + or -1  we can
    // reduce the demoniator so that the VR is
    //   (1 - \frac{ \sum_i w_i (Pol_{i,obs} - Pol_{i,est})^2 }
    //             { \sum_i w_i })
    // Furthermore, we can incorporate this term into the weighting
    // function s.t.  Sigma/trace(Sigma) to obtain
    //   (1 - \frac{ \sum_i \hat{Sigma} (Pol_{i,obs} - Pol_{i,est})^2)
    // Finally, because we want likelihoods we want to rescale the
    // residual squared from [0,4] to [0,1].  To do this we multiply
    // the weight by 0.5 s.t. the residual [-2,2] goes to [-1,1] which,
    // when squared, fits in the desired bounds of [0,1]; i.e. multiply by
    // (1/2)^2 = 0.25
    traceSigma = array_sum64f(nPolarity, Sigma, &ierr);
    if (ierr != 0){traceSigma = 1.0;}
    sqrtSigmaWt = memory_calloc64f(nPolarity);
#ifdef __INTEL_COMPILER
    __assume_aligned(sqrtSigmaWt, ISCL_MEMORY_ALIGN);
#endif
    for (i=0; i<nPolarity; i++){sqrtSigmaWt[i] = sqrt(0.25*Sigma[i]/traceSigma);}
    #pragma omp parallel \
     private (i, ic, ipol, jmt, kmt, Ncols, res, res2, U, Usign) \
     shared (G, Dmat, mts, nmtBlocks, phi, sqrtSigmaWt) \
     default (none) reduction(+:ierr)
    {
    U = memory_calloc64f(mblock*nPolarity);
    res2 = memory_calloc64f(mblock);
    Usign = memory_calloc64f(mblock*nPolarity);
#ifdef __INTEL_COMPILER
    __assume_aligned(res2, ISCL_MEMORY_ALIGN);
    __assume_aligned(Usign, ISCL_MEMORY_ALIGN);
#endif
    #pragma omp for
    for (kmt=0; kmt<nmtBlocks; kmt++)
    {
        jmt = kmt*blockSize;
        Ncols = MIN(blockSize, nmt - jmt); // Number of columns of M
        // Compute U = GM
        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
                    Mrows, Ncols, Kcols, one, G, LDG,
                    &mts[ldm*jmt], ldm, zero, U, mblock);
        // Make the theoretical polarity +1 or -1 to match the data
        array_copysign64f_work(mblock*nPolarity, U, Usign);
        memset(res2, 0, (size_t) mblock*sizeof(double));
        // Compute the weighted residual
        for (ipol=0; ipol<nPolarity; ipol++)
        {
            #pragma omp simd aligned(Dmat, Usign, res2: ISCL_MEMORY_ALIGN)
            for (ic=0; ic<Ncols; ic++)
            {
                // Change residual range from [0,2] to [0,1] and weight
                res = sqrtSigmaWt[ipol]*( Dmat[ipol*mblock+ic]
                                        - Usign[ipol*mblock+ic]);
                res2[ic] = res2[ic] + res*res;
            }
        }
        // Compute the variance reduction and put it into objective function
        for (ic=0; ic<Ncols; ic++)
        {
            phi[jmt+ic] = 1.0 - res2[ic];
        }
    } // Loop on moment tensor blocks
    memory_free64f(&Usign);
    memory_free64f(&U);
    memory_free64f(&res2);
    } // end parallel section
    memory_free64f(&sqrtSigmaWt);
    return 0;
}

static int computePadding64f(const int n)
{
    size_t mod, pad;
    int ipad;
    // Set space and make G matrix
    pad = 0;
    mod = ((size_t) n*sizeof(double))%64;
    if (mod != 0)
    {
        pad = (64 - mod)/sizeof(double);
    }
    ipad = (int) pad;
    return ipad;
}

ext: c


text:
\documentclass{report}
\usepackage[T1]{fontenc}
\usepackage{bera}

\usepackage[pdftex,usenames,dvipsnames]{color}
\usepackage{listings}
\definecolor{mycode}{rgb}{0.9,0.9,1}
\lstset{language=c++,tabsize=1,basicstyle=\scriptsize,backgroundcolor=\color{mycode}}

\usepackage{hyperref}
\usepackage{fancyhdr}
\pagestyle{fancy}
\usepackage{verbatim}
\begin{document}

%\title{The Powell Class for Minimization of Functions}
%\author{G.A.}
%\maketitle

\begin{comment}
@o Minimizer.h -t
@{
/*
@i license.txt
*/
@}
\end{comment}

\section{A class to minimize a function without using Derivatives}
This requires the GNU Scientific Library, libraries and devel. headers.

Explain usage here. FIXME.

This is the structure of this file:

@o Minimizer.h -t
@{
#ifndef MINIMIZER_H
#define MINIMIZER_H

#include <iostream>
#include <vector>
#include <stdexcept>
extern "C" {
#include <gsl/gsl_multimin.h>
}

namespace PsimagLite {
	@<MockVector@>
	@<MyFunction@>
	@<theClassHere@>
} // namespace PsimagLite
#endif // MINIMIZER_H

@}

And this is the class here:

@d theClassHere
@{
template<typename RealType,typename FunctionType>
class Minimizer {
	@<privateTypedefsAndConstants@>
public:
	@<constructor@>
	@<destructor@>
	@<publicFunctions@>
private:
	@<privateFunctions@>
	@<privateData@>
}; // class Minimizer
@}

@d privateTypedefsAndConstants
@{
typedef typename FunctionType::FieldType FieldType;
typedef typename Vector<FieldType>::Type VectorType;
typedef Minimizer<RealType,FunctionType> ThisType;

@}

@d privateData
@{
FunctionType& function_;
SizeType maxIter_;
const gsl_multimin_fminimizer_type *gslT_;
gsl_multimin_fminimizer *gslS_;
@}

@d constructor
@{
Minimizer(FunctionType& function,SizeType maxIter)
		: function_(function),
		  maxIter_(maxIter),
		  gslT_(gsl_multimin_fminimizer_nmsimplex2),
		  gslS_(gsl_multimin_fminimizer_alloc(gslT_,function_.size()))
{
}
@}

@d publicFunctions
@{
@<simplex@>
@}

@d destructor
@{
~Minimizer()
{
	gsl_multimin_fminimizer_free (gslS_);
}
@}

@d simplex
@{
int simplex(VectorType& minVector,RealType delta=1e-3,RealType tolerance=1e-3)
{
	gsl_vector *x;
	/* Starting point,  */
	x = gsl_vector_alloc (function_.size());
	for (SizeType i=0;i<minVector.size();i++)
		gsl_vector_set (x, i, minVector[i]);

	gsl_vector *xs;
	xs = gsl_vector_alloc (function_.size());
	for (SizeType i=0;i<minVector.size();i++)
		gsl_vector_set (xs, i, delta);

	gsl_multimin_function func;
	func.f= myFunction<FunctionType>;
	func.n = function_.size();
	func.params = &function_;
	gsl_multimin_fminimizer_set (gslS_, &func, x, xs);

	for (SizeType iter=0;iter<maxIter_;iter++) {
		int status = gsl_multimin_fminimizer_iterate (gslS_);

		if (status) throw RuntimeError("Minimizer::simplex(...): Error encountered\n");

		RealType size = gsl_multimin_fminimizer_size(gslS_);
		status = gsl_multimin_test_size(size, tolerance);

		if (status == GSL_SUCCESS) {
			found(minVector,gslS_->x,iter);
			gsl_vector_free (x);
			gsl_vector_free (xs);
			return iter;
		}
	}
	gsl_vector_free (x);
	gsl_vector_free (xs);
	return -1;
}
@}

@d privateFunctions
@{
@<found@>
@}

@d found
@{
void found(VectorType& minVector,gsl_vector* x,SizeType iter)
{
	for (SizeType i=0;i<minVector.size();i++)
		minVector[i] = gsl_vector_get(x,i);
}
@}

@d MockVector
@{
template<typename FieldType>
class MockVector {
public:
	MockVector(const gsl_vector *v) : v_(v)
	{
	}
	const FieldType& operator[](SizeType i) const
	{
		return v_->data[i];
	}
	SizeType size() const { return v_->size; }
private:
	const gsl_vector *v_;
}; // class MockVector
@}

@d MyFunction
@{
template<typename FunctionType>
typename FunctionType::FieldType myFunction(const gsl_vector *v, void *params)
{
	MockVector<typename FunctionType::FieldType> mv(v);
	FunctionType* ft = (FunctionType *)params;
	return ft->operator()(mv);
}
@}

\end{document}

ext: w


text:
#ifndef LPC_H
#define LPC_H

#include <cmath>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_blas.h>

gsl_vector *hanning(size_t L);

double lpcCoeffs(double *lpc, gsl_vector *x, size_t order);

#endif // LPC_H

ext: h

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
#pragma once

#include <cstdint>
#include <string>
#include <vector>

#include <gsl/span>

namespace terminal_editor {

/// Returns name of passed control character, or nullptr if given byte was not recognized.
/// ISO 30112 defines POSIX control characters as Unicode characters U+0000..U+001F, U+007F..U+009F, U+2028, and U+2029 (Unicode classes Cc, Zl, and Zp) "
/// See also: https://en.wikipedia.org/wiki/C0_and_C1_control_codes
/// @param codePoint    Code point for which name to return.
const char* controlCharacterName(uint32_t codePoint);

/// Describes return value of getFirstCodePoint() function.
struct CodePointInfo {
    bool valid;                          ///< True if valid code point was decoded. False otherwise.
    gsl::span<const char> consumedInput; ///< Bytes consumed from the input data. Length will be from 1 to 6.
    std::string info;                    ///< Arbitrary information about consumed bytes. If 'valid' is false will contain error information.
    uint32_t codePoint;                  ///< Decoded code point. Valid only if 'valid' is true.
};

/// Figures out what grapeheme is at the begining of data and returns it.
/// See: https://pl.wikipedia.org/wiki/UTF-8#Spos%C3%B3b_kodowania
/// @note All invalid byte sequences will be rendered as hex representations, with detailed explanation why it is invalid.
///       All control characters will be rendered as symbolic replacements.
///       Tab characters will be rendered as symbolic replacement.
/// @param data     Input string. It is assumed to be in UTF-8, but can contain invalid characters (which will be rendered as special Graphemes).
/// @param Returns number of bytes from input consumed.
CodePointInfo getFirstCodePoint(gsl::span<const char> data);

/// Parses a line of text into a list of CodePointInfos.
std::vector<CodePointInfo> parseLine(gsl::span<const char> inputData);

/// Analyzes given input data.
/// @param inputData     Input string. It is assumed to be in UTF-8, but can contain invalid characters (which will be annotated specially).
/// @return Valid UTF-8 string that describes what original string contains.
std::string analyzeData(gsl::span<const char> inputData);

/// Appends UTF-8 representation of given code point to a string.
/// @todo Move to some string utilities.
/// @param text         Text to append code point to.
/// @param codePoint    UTF-32 code point to append.
void appendCodePoint(std::string& text, uint32_t codePoint);

/// Appends UTF-8 representation of given code point to an output stream.
/// @todo Move to some string utilities.
/// @param os           Output stream to append code point to.
/// @param codePoint    UTF-32 code point to append.
void appendCodePoint(std::ostream& os, uint32_t codePoint);

} // namespace terminal_editor

max_stars_repo_name: Zbyl/terminal-editor

max_issues_repo_name: Zbyl/terminal-editor

max_forks_repo_name: Zbyl/terminal-editor


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>

int
main (void)
{
  const gsl_rng_type * T;
  gsl_rng * r;

  size_t i, k = 5, N = 100000;

  double * x = malloc (N * sizeof(double));
  double * small = malloc (k * sizeof(double));

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (i = 0; i < N; i++)
    {
      x[i] = gsl_rng_uniform(r);
    }

  gsl_sort_smallest (small, k, x, 1, N);

  printf ("%d smallest values from %d\n", k, N);

  for (i = 0; i < k; i++)
    {
      printf ("%d: %.18f\n", i, small[i]);
    }

  free (x);
  free (small);
  gsl_rng_free (r);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified