text:
/* specfunc/hyperg_2F0.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_hyperg.h>

#include "error.h"
#include "hyperg.h"

int
gsl_sf_hyperg_2F0_e(const double a, const double b, const double x, gsl_sf_result * result)
{
  if(x < 0.0) {
    /* Use "definition" 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x).
     */
    gsl_sf_result U;
    double pre = pow(-1.0/x, a);
    int stat_U = gsl_sf_hyperg_U_e(a, 1.0+a-b, -1.0/x, &U);
    result->val = pre * U.val;
    result->err = GSL_DBL_EPSILON * fabs(result->val) + pre * U.err;
    return stat_U;
  }
  else if(x == 0.0) {
    result->val = 1.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else {
    /* Use asymptotic series. ??
     */
    /* return hyperg_2F0_series(a, b, x, -1, result, &prec); */
    DOMAIN_ERROR(result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_hyperg_2F0(const double a, const double b, const double x)
{
  EVAL_RESULT(gsl_sf_hyperg_2F0_e(a, b, x, &result));
}


text:
// -*- C++ -*-  Copyright (c) Microsoft Corporation; see license.txt
#ifndef MESH_PROCESSING_LIBHH_MY_LAPACK_H_
#define MESH_PROCESSING_LIBHH_MY_LAPACK_H_

#include "Hh.h"

#if defined(HH_USE_LAPACK_INTEGER)      // Google

#include <lapack.h>
using lapack_int = lapack::integer;

#elif !(defined(_WIN32) || defined(__CYGWIN__)) // unix

#if defined(__GNUC__)
#pragma GCC diagnostic ignored "-Wredundant-decls"
#endif
#include <lapacke.h>            // in either /usr/include/ or /usr/include/lapacke/

#else  // cygwin or WIN32

#if defined(__CYGWIN__)

// cygwin: c:/cygwin/usr/src/lapack-3.5.0.tgz
//  lapacke_config.h
// #ifndef lapack_int
// #if defined(LAPACK_ILP64)
// #define lapack_int              long
// #else
// #define lapack_int              int
// #endif
// #endif
//
//  lapacke.h:  "#define lapack_int int"
// #define LAPACK_sgelss LAPACK_GLOBAL(sgelss, SGELSS)
// void LAPACK_sgelss( lapack_int* m, lapack_int* n, lapack_int* nrhs, float* a,
//                     lapack_int* lda, float* b, lapack_int* ldb, float* s,
//                     float* rcond, lapack_int* rank, float* work,
//                     lapack_int* lwork, lapack_int *info );

//  c:/cygwin/usr/src/lapack-3.5.0.tgz!lapack-3.5.0/lapacke/example/example_DGELS_colmajor.c
// #include <lapacke.h>
// {
//    double A[5][3] = {1, 2, 3, 4, 5, 1, 3, 5, 2, 4, 1, 4, 2, 5, 3};
//    double b[5][2] = {-10, 12, 14, 16, 18, -3, 14, 12, 16, 16};
//    lapack_int m = 5, n = 3, lda = 5, ldb = 5, nrhs = 2;
//    lapack_int info = LAPACKE_dgels(LAPACK_COL_MAJOR, 'N', m, n, nrhs, *A, lda, *b, ldb);
// }

using lapack_int = int;

#else  // _WIN32

using lapack_int = long;        // was defined as "integer" in f2c.h

#endif

// Continue: cygwin or WIN32

// clapack routines in liblapack
extern "C" {

// http://www.netlib.org/lapack/explore-html/d0/db8/group__real_g_esolve.html#ga206e3084597d088b31dc054a69aec93f
// SGELSS solves overdetermined or underdetermined systems for GE matrices, using SVD.
// ~/src/lib_src/CLAPACK/SRC/sgelss.c
    lapack_int sgelss_(lapack_int* m, lapack_int* n, lapack_int* nrhs, float* a, lapack_int* lda,
                       float* b, lapack_int* ldb, float* s, float* rcond, lapack_int* irank,
                       float* work, lapack_int* lwork, lapack_int* info);

// same for double-precision
    lapack_int dgelss_(lapack_int* m, lapack_int* n, lapack_int* nrhs, double* a, lapack_int* lda,
                       double* b, lapack_int* ldb, double* s, double* rcond, lapack_int* irank,
                       double* work, lapack_int* lwork, lapack_int* info);

// http://www.netlib.org/lapack/explore-html/d0/d2a/sgelsx_8f.html
// SGELSX solves overdetermined or underdetermined systems for GE matrices, using QR.
// This routine is deprecated and has been replaced by routine SGELSY.
// ~/src/lib_src/CLAPACK/SRC/sgelsx.c
// lapack_int sgelsx_(lapack_int* m, lapack_int* n, lapack_int* nrhs, float* a, lapack_int* lda,
//                   float* b, lapack_int* ldb, lapack_int* jpvt, float* rcond, lapack_int* irank,
//                   float* work, lapack_int* info);

// http://www.netlib.org/lapack/explore-html/dc/db6/sgelsy_8f.html
    lapack_int sgelsy_(lapack_int* m, lapack_int* n, lapack_int* nrhs, float* a, lapack_int* lda,
                       float* b, lapack_int* ldb, lapack_int* jpvt, float* rcond, lapack_int* irank,
                       float* work, lapack_int* lwork, lapack_int *info);

// http://www.netlib.org/lapack/explore-html/d4/dca/group__real_g_esing.html
// subroutine SGESVD (JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
// SGESVD computes the singular value decomposition (SVD) for GE matrices
// ~/src/lib_src/CLAPACK/SRC/sgesvd.c
    lapack_int sgesvd_(char* jobu, char* jobvt, lapack_int* m, lapack_int* n, float* a, lapack_int* lda,
                       float* s, float* u, lapack_int* ldu, float* vt, lapack_int* ldvt,
                       float* work, lapack_int* lwork, lapack_int* info);

// http://www.netlib.org/lapack/explore-html/d8/ddc/group__real_g_ecomputational.html#ga7cb54fa1727bf0166523036f4948bc56
// subroutine SGEQRF (M, N, A, LDA, TAU, WORK, LWORK, INFO)
// SGEQRF computes a QR factorization of a real M-by-N matrix A:
// A = Q * R.
// ~/src/lib_src/CLAPACK/SRC/sgeqrf.c
    lapack_int sgeqrf_(lapack_int* m, lapack_int* n, float* a, lapack_int* lda, float* tau,
                       float* work, lapack_int* lwork, lapack_int* info);

}

#if defined(_MSC_VER)
#if defined(_WIN64) && !defined(HH_NO_MKL) // for now, only use MKL on x64
#if !defined(HH_MKL_NOT_DLL)
// 2014
HH_REFERENCE_LIB("mkl_intel_lp64_dll.lib");
HH_REFERENCE_LIB("mkl_intel_thread_dll.lib");
HH_REFERENCE_LIB("mkl_core_dll.lib");
HH_REFERENCE_LIB("libiomp5md.lib");
// e.g.: depends MeshSimplify.exe
//  -> mkl_intel_thread.dll -> libiomp5md.dll mkl_core.dll ("core" dispatches among several mkl_*.dll based on CPU)
#else  // !defined(HH_MKL_NOT_DLL)
// 2012
#if defined(_WIN64)
HH_REFERENCE_LIB("mkl_intel_lp64.lib");
#else
HH_REFERENCE_LIB("mkl_intel_c.lib");
#endif
HH_REFERENCE_LIB("mkl_intel_thread.lib");
HH_REFERENCE_LIB("mkl_core.lib");
HH_REFERENCE_LIB("libiomp5md.lib");
#endif  // !defined(HH_MKL_NOT_DLL)
#else   // defined(_WIN64) && !defined(HH_NO_MKL)
HH_REFERENCE_LIB("liblapack.lib");
HH_REFERENCE_LIB("libBLAS.lib");
#if 1
HH_REFERENCE_LIB("libf2c.lib"); // CLAPACK
#else
HH_REFERENCE_LIB("libF77.lib"); // liblapack
HH_REFERENCE_LIB("libI77.lib");
#endif
#endif  // defined(_WIN64) && !defined(HH_NO_MKL)
#endif  // defined(_MSC_VER)

#endif  // !(defined(_WIN32) || defined(__CYGWIN__))

#if 0
// from mkl_lapack.h:
// It shows that the MKL lapack functions have almost the same signatures (except for "const char*").
using MKL_INT = int;            // 32-bit unless MKL_ILP64
extern "C" {
    void SGESVD(const char* jobu, const char* jobvt, const MKL_INT* m,
                const MKL_INT* n, float* a, const MKL_INT* lda, float* s,
                float* u, const MKL_INT* ldu, float* vt, const MKL_INT* ldvt,
                float* work, const MKL_INT* lwork, MKL_INT* info);
    void SGESVD_(const char* jobu, const char* jobvt, const MKL_INT* m,
                 const MKL_INT* n, float* a, const MKL_INT* lda, float* s,
                 float* u, const MKL_INT* ldu, float* vt, const MKL_INT* ldvt,
                 float* work, const MKL_INT* lwork, MKL_INT* info);
    void sgesvd(const char* jobu, const char* jobvt, const MKL_INT* m,
                const MKL_INT* n, float* a, const MKL_INT* lda, float* s,
                float* u, const MKL_INT* ldu, float* vt, const MKL_INT* ldvt,
                float* work, const MKL_INT* lwork, MKL_INT* info);
    void sgesvd_(const char* jobu, const char* jobvt, const MKL_INT* m,
                 const MKL_INT* n, float* a, const MKL_INT* lda, float* s,
                 float* u, const MKL_INT* ldu, float* vt, const MKL_INT* ldvt,
                 float* work, const MKL_INT* lwork, MKL_INT* info);
}

// from /usr/include/lapacke.h:

lapack_int LAPACKE_sgelss( int matrix_order, lapack_int m, lapack_int n,
                           lapack_int nrhs, float* a, lapack_int lda, float* b,
                           lapack_int ldb, float* s, float rcond,
                           lapack_int* rank );

void LAPACK_sgelss( lapack_int* m, lapack_int* n, lapack_int* nrhs, float* a,
                    lapack_int* lda, float* b, lapack_int* ldb, float* s,
                    float* rcond, lapack_int* rank, float* work,
                    lapack_int* lwork, lapack_int *info );

#endif

#endif // MESH_PROCESSING_LIBHH_MY_LAPACK_H_

ext: h

max_stars_repo_name: xianyuMeng/Mesh-processing-library

max_issues_repo_name: xianyuMeng/Mesh-processing-library

max_forks_repo_name: xianyuMeng/Mesh-processing-library


text:
//
// Created by pedram pakseresht on 2/18/21.
//

static char help[] = "settling in quiescent fluid";

#include <petsc.h>
#include "incompressibleFlow.h"
#include "mesh.h"
#include "particleInertial.h"
#include "particles.h"
#include "particleInitializer.h"
#include "petscviewer.h"

static PetscErrorCode uniform_u(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscInt d;
    for (d = 0; d < Dim; ++d)
        u[d] = 0.0;
    return 0;

}

static PetscErrorCode uniform_u_t(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *u, void *ctx) {
    PetscInt d;
    for (d = 0; d < Dim; ++d)
        u[d] = 0.0;
    return 0;
}

static PetscErrorCode uniform_p(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *p, void *ctx) {
    p[0] = 0.0;
    return 0;
}

static PetscErrorCode uniform_T(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *T, void *ctx) {
    T[0] = 0.0;
    return 0;
}
static PetscErrorCode uniform_T_t(PetscInt Dim, PetscReal time, const PetscReal *X, PetscInt Nf, PetscScalar *T, void *ctx) {
    T[0] = 0.0;
    return 0;
}

static PetscErrorCode SetInitialConditions(TS ts, Vec u) {

    PetscErrorCode (*initFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx) = {uniform_u,uniform_p,uniform_T};
    // u here is the solution vector including velocity, temperature and pressure fields.
    DM dm;
    PetscReal t;
    PetscErrorCode ierr;
    PetscFunctionBegin;

    ierr = TSGetDM(ts, &dm);
    CHKERRQ(ierr);
    ierr = TSGetTime(ts, &t);

    CHKERRQ(ierr);

    ierr = DMProjectFunction(dm, 0.0, initFuncs, NULL, INSERT_ALL_VALUES, u);
    CHKERRQ(ierr);
    // get the flow to apply the completeFlowInitialization method
    ierr = IncompressibleFlow_CompleteFlowInitialization(dm, u);
    CHKERRQ(ierr);

    PetscFunctionReturn(0);
}

static PetscErrorCode MonitorFlowAndParticleError(TS ts, PetscInt step, PetscReal crtime, Vec u, void *ctx) {
    // PetscErrorCode (*exactFuncs[3])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
    // void *ctxs[3];
    DM dm;
    PetscDS ds;
    Vec v;
    PetscReal ferrors[3];
    PetscInt f;
    PetscErrorCode ierr;
    PetscInt num;

    PetscFunctionBeginUser;
    ierr = TSGetDM(ts, &dm);
    CHKERRQ(ierr);
    ierr = DMGetDS(dm, &ds);
    CHKERRQ(ierr);

    // get the particle data from the context
    ParticleData particlesData = (ParticleData)ctx;
    PetscInt particleCount;
    ierr = DMSwarmGetSize(particlesData->dm, &particleCount);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // compute the average particle location
    const PetscReal *coords;
    PetscInt dims;
    PetscReal avg[3] = {0.0, 0.0, 0.0};
    ierr = DMSwarmGetField(particlesData->dm, DMSwarmPICField_coor, &dims, NULL, (void **)&coords);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    for (PetscInt p = 0; p < particleCount; p++) {
        for (PetscInt n = 0; n < dims; n++) {
            avg[n] += coords[p * dims + n] / particleCount;  // PetscReal
        }
    }
    ierr = DMSwarmRestoreField(particlesData->dm, DMSwarmPICField_coor, &dims, NULL, (void **)&coords);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD,
                       "Timestep: %04d time = %-8.4g \t L_2 Error: [%2.3g, %2.3g, %2.3g] ParticleCount: %d\n",
                       (int)step,
                       (double)crtime,
                       (double)ferrors[0],
                       (double)ferrors[1],
                       (double)ferrors[2],
                       particleCount,
                       (double)avg[0],
                       (double)avg[1]);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Avg Particle Location: [%2.3g, %2.3g, %2.3g]\n", (double)avg[0], (double)avg[1], (double)avg[2]);

    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = VecViewFromOptions(u, NULL, "-vec_view");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = DMSetOutputSequenceNumber(particlesData->dm, step, crtime);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    if (step == 0) {
        ierr = ParticleViewFromOptions(particlesData, "-particle_init");
        CHKERRABORT(PETSC_COMM_WORLD, ierr);
    } else {
        ierr = ParticleViewFromOptions(particlesData, "-particle_view");
        CHKERRABORT(PETSC_COMM_WORLD, ierr);
    }

    PetscFunctionReturn(0);
}

static PetscErrorCode ParticleInertialInitialize(ParticleData particles) {

    PetscFunctionBeginUser;
    PetscErrorCode  ierr;
    Vec vel,diam,dens;
    DM particleDm = particles->dm;

    // input parameters required for particles
    PetscScalar partVel = 0.0;
    PetscScalar partDiam = 0.22;
    PetscScalar partDens = 90.0;
    PetscScalar fluidDens = 1.0;
    PetscScalar fluidVisc = 1.0;
    PetscScalar gravity[3] = {0.0,1.0,0.0};

    ierr = DMSwarmCreateGlobalVectorFromField(particleDm,ParticleVelocity, &vel);CHKERRQ(ierr);
    ierr = DMSwarmCreateGlobalVectorFromField(particleDm,ParticleDiameter, &diam);CHKERRQ(ierr);
    ierr = DMSwarmCreateGlobalVectorFromField(particleDm,ParticleDensity, &dens);CHKERRQ(ierr);
    // set particle velocity, diameter and density
    ierr = VecSet(vel, partVel);CHKERRQ(ierr);
    ierr = VecSet(diam, partDiam);CHKERRQ(ierr);
    ierr = VecSet(dens, partDens);CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(particleDm,ParticleVelocity, &vel);CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(particleDm,ParticleDiameter, &diam);CHKERRQ(ierr);
    ierr = DMSwarmDestroyGlobalVectorFromField(particleDm,ParticleDensity, &dens);CHKERRQ(ierr);

    InertialParticleParameters *data;
    PetscNew(&data);
    particles->data =data;

    // set fluid parameters
    data->fluidDensity = fluidDens;
    data->fluidViscosity = fluidVisc;
    data->gravityField[0] = gravity[0];
    data->gravityField[1] = gravity[1];
    data->gravityField[2] = gravity[2];
    PetscFunctionReturn(0);

}

int main( int argc, char *argv[] )
{

    DM dm;     // domain definition
    TS ts;     // time-stepper
    PetscErrorCode ierr;

    PetscBag parameterBag; // constant flow parameters
    Vec flowField;         // flow solution vector
    FlowData flowData;

    PetscReal t;
    PetscInt Dim = 2;
    PetscReal dt = 0.05;       // dt for time stepper
    PetscInt max_steps = 10; // maximum time steps

    PetscReal Re = 0.1;       // Reynolds number
    PetscReal St = 1.0;       // Strouhal number
    PetscReal Pe = 1.0;       // Peclet number
    PetscReal Mu = 1.0;      // viscosity
    PetscReal K_input = 1.0;  // thermal conductivity
    PetscReal Cp = 1.0;    // heat capacity
    // PetscInt Np=10;

    // initialize Petsc ...
    ierr = PetscInitialize(&argc, &argv, NULL, "");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Settling particle in a quiescent fluid \n");CHKERRQ(ierr);

    // setup the ts
    ierr = TSCreate(PETSC_COMM_WORLD, &ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = CreateMesh(PETSC_COMM_WORLD, &dm, PETSC_TRUE, Dim);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetDM(ts, dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //output the mesh
    ierr = DMViewFromOptions(dm, NULL, "-dm_view");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup the flow data
    ierr = FlowCreate(&flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // setup problem
    ierr = IncompressibleFlow_SetupDiscretization(flowData, dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    IncompressibleFlowParameters flowParameters;

    // changing non-dimensional parameters manually here ...
    flowParameters.strouhal = St;
    flowParameters.reynolds = Re;
    flowParameters.peclet = Pe;
    flowParameters.mu = Mu;
    flowParameters.k = K_input;
    flowParameters.cp = Cp;

    // print out the parameters
    /*
    ierr = PetscPrintf(PETSC_COMM_WORLD, "*** non-dimensional parameters ***\n");CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "St=%g\n", flowParameters.strouhal);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Re=%g\n", flowParameters.reynolds);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Pe=%g\n", flowParameters.peclet);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Mu=%g\n", flowParameters.mu);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "K=%g\n", flowParameters.k);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Cp=%g\n", flowParameters.cp);CHKERRQ(ierr);
    */

    // Start the problem setup
    PetscScalar constants[TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS];
    ierr = IncompressibleFlow_PackParameters(&flowParameters, constants);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = IncompressibleFlow_StartProblemSetup(flowData, TOTAL_INCOMPRESSIBLE_FLOW_PARAMETERS, constants);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Override problem with source terms, boundary, and set the exact solution
    PetscDS prob;
    ierr = DMGetDS(dm, &prob);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup Boundary Conditions
    // Note: DM_BC_ESSENTIAL is a Dirichlet BC.
    PetscInt id;
    id = 3;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "top wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    id = 1;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "bottom wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 2;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "right wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 4;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "left wall velocity", "marker", VEL, 0, NULL, (void (*)(void))uniform_u, (void (*)(void))uniform_u_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 3;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "top wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 1;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "bottom wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 2;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "right wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    id = 4;
    ierr = PetscDSAddBoundary(
        prob, DM_BC_ESSENTIAL, "left wall temp", "marker", TEMP, 0, NULL, (void (*)(void))uniform_T, (void (*)(void))uniform_T_t, 1, &id, parameterBag);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = IncompressibleFlow_CompleteProblemSetup(flowData, ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Name the flow field
    ierr = PetscObjectSetName(((PetscObject)flowData->flowField), "Numerical Solution");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = SetInitialConditions(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = TSGetTime(ts, &t);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //  *** added by Pedram for dt and number of time steps ***
    ierr = TSSetTimeStep(ts,dt);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = TSSetMaxSteps(ts,max_steps);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    // ***********************************

    ierr = DMSetOutputSequenceNumber(dm, 0, t);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // checks for convergence ...
    ierr = DMTSCheckFromOptions(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ParticleData particles;

    // Setup the particle domain
    ierr = ParticleInertialCreate(&particles, Dim);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // link the flow to the particles
    ierr = ParticleInitializeFlow(particles, flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // name the particle domain
    ierr = PetscObjectSetOptionsPrefix((PetscObject)(particles->dm), "particles_");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscObjectSetName((PetscObject)particles->dm, "Particles");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // initialize the particles position
    ierr = ParticleInitialize(dm, particles->dm);CHKERRABORT(PETSC_COMM_WORLD, ierr);
    // initialize inertial particles velocity, diameter and density
    ierr = ParticleInertialInitialize(particles);CHKERRABORT(PETSC_COMM_WORLD, ierr);

    //ierr = ParticleInertialInitialize(particles->dm);CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // setup the flow monitor to also check particles
    ierr = TSMonitorSet(ts, MonitorFlowAndParticleError, particles, NULL);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSSetFromOptions(ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Setup particle position integrator
    TS particleTs;
    ierr = TSCreate(PETSC_COMM_WORLD, &particleTs);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscObjectSetOptionsPrefix((PetscObject)particleTs, "particle_");
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    ierr = ParticleInertialSetupIntegrator(particles, particleTs, flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    // ierr = VecView(particleVelocity,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

    // Solve the one way coupled system
    ierr = TSSolve(ts, flowData->flowField);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);

    // Cleanup
    ierr = DMDestroy(&dm);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSDestroy(&ts);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = TSDestroy(&particleTs);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = FlowDestroy(&flowData);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = ParticleInertialDestroy(&particles);
    CHKERRABORT(PETSC_COMM_WORLD, ierr);
    ierr = PetscFinalize();
    exit(ierr);

}

ext: c

max_stars_repo_name: pakserep/ablateClient

max_issues_repo_name: pakserep/ablateClient

max_forks_repo_name: pakserep/ablateClient


text:
/* interpolation/cspline.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include "integ_eval.h"
#include <gsl/gsl_interp.h>

typedef struct
{
  double * c;
  double * g;
  double * diag;
  double * offdiag;
} cspline_state_t;

/* common initialization */
static void *
cspline_alloc (size_t size)
{
  cspline_state_t * state = (cspline_state_t *) malloc (sizeof (cspline_state_t));

  if (state == NULL)
    {
      GSL_ERROR_NULL("failed to allocate space for state", GSL_ENOMEM);
    }

  state->c = (double *) malloc (size * sizeof (double));

  if (state->c == NULL)
    {
      free (state);
      GSL_ERROR_NULL("failed to allocate space for c", GSL_ENOMEM);
    }

  state->g = (double *) malloc (size * sizeof (double));

  if (state->g == NULL)
    {
      free (state->c);
      free (state);
      GSL_ERROR_NULL("failed to allocate space for g", GSL_ENOMEM);
    }

  state->diag = (double *) malloc (size * sizeof (double));

  if (state->diag == NULL)
    {
      free (state->g);
      free (state->c);
      free (state);
      GSL_ERROR_NULL("failed to allocate space for diag", GSL_ENOMEM);
    }

  state->offdiag = (double *) malloc (size * sizeof (double));

  if (state->offdiag == NULL)
    {
      free (state->diag);
      free (state->g);
      free (state->c);
      free (state);
      GSL_ERROR_NULL("failed to allocate space for offdiag", GSL_ENOMEM);
    }

  return state;
}

/* natural spline calculation
 * see [Engeln-Mullges + Uhlig, p. 254]
 */
static int
cspline_init (void * vstate, const double xa[], const double ya[],
              size_t size)
{
  cspline_state_t *state = (cspline_state_t *) vstate;

  size_t i;
  size_t num_points = size;
  size_t max_index = num_points - 1;  /* Engeln-Mullges + Uhlig "n" */
  size_t sys_size = max_index - 1;    /* linear system is sys_size x sys_size */

  state->c[0] = 0.0;
  state->c[max_index] = 0.0;

  for (i = 0; i < sys_size; i++)
    {
      const double h_i   = xa[i + 1] - xa[i];
      const double h_ip1 = xa[i + 2] - xa[i + 1];
      const double ydiff_i   = ya[i + 1] - ya[i];
      const double ydiff_ip1 = ya[i + 2] - ya[i + 1];
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
      state->offdiag[i] = h_ip1;
      state->diag[i] = 2.0 * (h_ip1 + h_i);
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 -  ydiff_i * g_i);
    }

  if (sys_size == 1)
    {
      state->c[1] = state->g[0] / state->diag[0];
      return GSL_SUCCESS;
    }
  else
    {
      gsl_vector_view g_vec = gsl_vector_view_array(state->g, sys_size);
      gsl_vector_view diag_vec = gsl_vector_view_array(state->diag, sys_size);
      gsl_vector_view offdiag_vec = gsl_vector_view_array(state->offdiag, sys_size - 1);
      gsl_vector_view solution_vec = gsl_vector_view_array ((state->c) + 1, sys_size);

      int status = gsl_linalg_solve_symm_tridiag(&diag_vec.vector,
                                                 &offdiag_vec.vector,
                                                 &g_vec.vector,
                                                 &solution_vec.vector);
      return status;
    }
}

/* periodic spline calculation
 * see [Engeln-Mullges + Uhlig, p. 256]
 */
static int
cspline_init_periodic (void * vstate, const double xa[], const double ya[],
                       size_t size)
{
  cspline_state_t *state = (cspline_state_t *) vstate;

  size_t i;
  size_t num_points = size;
  size_t max_index = num_points - 1;  /* Engeln-Mullges + Uhlig "n" */
  size_t sys_size = max_index;    /* linear system is sys_size x sys_size */

  if (sys_size == 2) {
    /* solve 2x2 system */

    const double h0 = xa[1] - xa[0];
    const double h1 = xa[2] - xa[1];

    const double A = 2.0*(h0 + h1);
    const double B = h0 + h1;
    double g[2];
    double det;

    g[0] = 3.0 * ((ya[2] - ya[1]) / h1 - (ya[1] - ya[0]) / h0);
    g[1] = 3.0 * ((ya[1] - ya[2]) / h0 - (ya[2] - ya[1]) / h1);

    det = 3.0 * (h0 + h1) * (h0 + h1);
    state->c[1] = ( A * g[0] - B * g[1])/det;
    state->c[2] = (-B * g[0] + A * g[1])/det;
    state->c[0] = state->c[2];

    return GSL_SUCCESS;
  } else {

    for (i = 0; i < sys_size-1; i++) {
      const double h_i       = xa[i + 1] - xa[i];
      const double h_ip1     = xa[i + 2] - xa[i + 1];
      const double ydiff_i   = ya[i + 1] - ya[i];
      const double ydiff_ip1 = ya[i + 2] - ya[i + 1];
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
      state->offdiag[i] = h_ip1;
      state->diag[i] = 2.0 * (h_ip1 + h_i);
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i);
    }

    i = sys_size - 1;

    {
      const double h_i       = xa[i + 1] - xa[i];
      const double h_ip1     = xa[1] - xa[0];
      const double ydiff_i   = ya[i + 1] - ya[i];
      const double ydiff_ip1 = ya[1] - ya[0];
      const double g_i = (h_i != 0.0) ? 1.0 / h_i : 0.0;
      const double g_ip1 = (h_ip1 != 0.0) ? 1.0 / h_ip1 : 0.0;
      state->offdiag[i] = h_ip1;
      state->diag[i] = 2.0 * (h_ip1 + h_i);
      state->g[i] = 3.0 * (ydiff_ip1 * g_ip1 - ydiff_i * g_i);
    }

    {
      gsl_vector_view g_vec = gsl_vector_view_array(state->g, sys_size);
      gsl_vector_view diag_vec = gsl_vector_view_array(state->diag, sys_size);
      gsl_vector_view offdiag_vec = gsl_vector_view_array(state->offdiag, sys_size);
      gsl_vector_view solution_vec = gsl_vector_view_array ((state->c) + 1, sys_size);

      int status = gsl_linalg_solve_symm_cyc_tridiag(&diag_vec.vector,
                                                     &offdiag_vec.vector,
                                                     &g_vec.vector,
                                                     &solution_vec.vector);
      state->c[0] = state->c[max_index];

      return status;
    }
  }
}

static
void
cspline_free (void * vstate)
{
  cspline_state_t *state = (cspline_state_t *) vstate;

  free (state->c);
  free (state->g);
  free (state->diag);
  free (state->offdiag);
  free (state);
}

/* function for common coefficient determination
 */
static inline void
coeff_calc (const double c_array[], double dy, double dx, size_t index,
            double * b, double * c, double * d)
{
  const double c_i = c_array[index];
  const double c_ip1 = c_array[index + 1];
  *b = (dy / dx) - dx * (c_ip1 + 2.0 * c_i) / 3.0;
  *c = c_i;
  *d = (c_ip1 - c_i) / (3.0 * dx);
}

static
int
cspline_eval (const void * vstate,
              const double x_array[], const double y_array[], size_t size,
              double x,
              gsl_interp_accel * a,
              double *y)
{
  const cspline_state_t *state = (const cspline_state_t *) vstate;

  double x_lo, x_hi;
  double dx;
  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  x_hi = x_array[index + 1];
  x_lo = x_array[index];
  dx = x_hi - x_lo;
  if (dx > 0.0)
    {
      const double y_lo = y_array[index];
      const double y_hi = y_array[index + 1];
      const double dy = y_hi - y_lo;
      double delx = x - x_lo;
      double b_i, c_i, d_i;
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
      *y = y_lo + delx * (b_i + delx * (c_i + delx * d_i));
      return GSL_SUCCESS;
    }
  else
    {
      *y = 0.0;
      return GSL_EINVAL;
    }
}

static
int
cspline_eval_deriv (const void * vstate,
                    const double x_array[], const double y_array[], size_t size,
                    double x,
                    gsl_interp_accel * a,
                    double *dydx)
{
  const cspline_state_t *state = (const cspline_state_t *) vstate;

  double x_lo, x_hi;
  double dx;
  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  x_hi = x_array[index + 1];
  x_lo = x_array[index];
  dx = x_hi - x_lo;
  if (dx > 0.0)
    {
      const double y_lo = y_array[index];
      const double y_hi = y_array[index + 1];
      const double dy = y_hi - y_lo;
      double delx = x - x_lo;
      double b_i, c_i, d_i;
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
      *dydx = b_i + delx * (2.0 * c_i + 3.0 * d_i * delx);
      return GSL_SUCCESS;
    }
  else
    {
      *dydx = 0.0;
      return GSL_EINVAL;
    }
}

static
int
cspline_eval_deriv2 (const void * vstate,
                     const double x_array[], const double y_array[], size_t size,
                     double x,
                     gsl_interp_accel * a,
                     double * y_pp)
{
  const cspline_state_t *state = (const cspline_state_t *) vstate;

  double x_lo, x_hi;
  double dx;
  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  x_hi = x_array[index + 1];
  x_lo = x_array[index];
  dx = x_hi - x_lo;
  if (dx > 0.0)
    {
      const double y_lo = y_array[index];
      const double y_hi = y_array[index + 1];
      const double dy = y_hi - y_lo;
      double delx = x - x_lo;
      double b_i, c_i, d_i;
      coeff_calc(state->c, dy, dx, index,  &b_i, &c_i, &d_i);
      *y_pp = 2.0 * c_i + 6.0 * d_i * delx;
      return GSL_SUCCESS;
    }
  else
    {
      *y_pp = 0.0;
      return GSL_EINVAL;
    }
}

static
int
cspline_eval_integ (const void * vstate,
                    const double x_array[], const double y_array[], size_t size,
                    gsl_interp_accel * acc,
                    double a, double b,
                    double * result)
{
  const cspline_state_t *state = (const cspline_state_t *) vstate;

  size_t i, index_a, index_b;

  if (acc != 0)
    {
      index_a = gsl_interp_accel_find (acc, x_array, size, a);
      index_b = gsl_interp_accel_find (acc, x_array, size, b);
    }
  else
    {
      index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
      index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
    }

  *result = 0.0;

  /* interior intervals */
  for(i=index_a; i<=index_b; i++) {
    const double x_hi = x_array[i + 1];
    const double x_lo = x_array[i];
    const double y_lo = y_array[i];
    const double y_hi = y_array[i + 1];
    const double dx = x_hi - x_lo;
    const double dy = y_hi - y_lo;
    if(dx != 0.0) {
      double b_i, c_i, d_i;
      coeff_calc(state->c, dy, dx, i,  &b_i, &c_i, &d_i);

      if (i == index_a || i == index_b)
        {
          double x1 = (i == index_a) ? a : x_lo;
          double x2 = (i == index_b) ? b : x_hi;
          *result += integ_eval(y_lo, b_i, c_i, d_i, x_lo, x1, x2);
        }
      else
        {
          *result += dx * (y_lo + dx*(0.5*b_i + dx*(c_i/3.0 + 0.25*d_i*dx)));
        }
    }
    else {
      *result = 0.0;
      return GSL_EINVAL;
    }
  }

  return GSL_SUCCESS;
}

static const gsl_interp_type cspline_type =
{
  "cspline",
  3,
  &cspline_alloc,
  &cspline_init,
  &cspline_eval,
  &cspline_eval_deriv,
  &cspline_eval_deriv2,
  &cspline_eval_integ,
  &cspline_free
};

const gsl_interp_type * gsl_interp_cspline = &cspline_type;

static const gsl_interp_type cspline_periodic_type =
{
  "cspline-periodic",
  2,
  &cspline_alloc,
  &cspline_init_periodic,
  &cspline_eval,
  &cspline_eval_deriv,
  &cspline_eval_deriv2,
  &cspline_eval_integ,
  &cspline_free
};

const gsl_interp_type * gsl_interp_cspline_periodic = &cspline_periodic_type;

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "cmopad.h"

/*
Copyright (C) 2010
Lars Krieger & Sebastian Heimann

Contact
lars.krieger@zmaw.de  &  sebastian.heimann@zmaw.de

Modifications - 2015
Converted original Python mopad.py to C - Ben Baker (ISTI)

Contact
benbaker@isti.com

#######################################################################

License:

GNU Lesser General Public License, Version 3

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
02110-1301, USA.
*/

#define rad2deg 180.0/M_PI
#define epsilon 1.e-13
static double numpy_mod(double a, double b);
static int __cmopad_argsort3(double *a, int *perm);
static int __cmopad_inv3(double Amat[3][3]);
static double __cmopad_determinant3x3(double A[3][3]);
static double __cmopad_trace3(double M[3][3]);
static int __cmopad_Eigs3x3(double a[3][3], int job, double eigs[3]);
static void __cmopad_swap8(double *a, double *b);
static void array_cross3(const double *__restrict__ u,
                         const double *__restrict__ v,
                         double *__restrict__ n);

//============================================================================//
/*!
 * @brief Generates basis for coordinate transform switch
 *
 * @param[in] in_system    input system: NED, USE, XYZ, NWU
 * @param[in] out_system   output system: NED, USE, XYZ, NWU
 *
 * @param[out] r           rotation matrix [3 x 3]
 *
 * @result 0 in indicates success
 *
 */
int cmopad_basis_switcher(enum cmopad_basis_enum in_system,
                          enum cmopad_basis_enum out_system, double r[3][3])
{
    const char *fcnm = "cmopad_basis_switcher\0";
    double From[9], To[9], r9[9];
    double alpha = 1.0;
    double beta  = 0.0;
    int n3 = 3;
    int i, j;
    //------------------------------------------------------------------------//
    //
    // null out to and from
    for (i=0; i<9; i++)
    {
        To[i]   = 0.0;
        From[i] = 0.0;
    }
    // Classify in
    if (in_system == NED)
    {
        // Inverse of identity: From[0] = 1.0; From[4] = 1.0; From[8] = 1.0 is:
        From[0] = 1.0; From[4] = 1.0; From[8] = 1.0;
    }
    else if (in_system == USE)
    {
        From[3] =-1.0; From[7] = 1.0; From[2] =-1.0;
    }
    else if (in_system == XYZ)
    {
        From[3] = 1.0; From[1] = 1.0; From[8] =-1.0;
    }
    else if (in_system == NWU)
    {
        From[0] = 1.0; From[4] =-1.0; From[8] =-1.0;
    }
    else
    {
        printf("%s: Invalid in coordinate system %d\n", fcnm, in_system);
        return -1;
    }
    // Classify out
    if (out_system == NED)
    {
        To[0] = 1.0; To[4] = 1.0; To[8] = 1.0;
    }
    else if (out_system == USE)
    {
        // Inverse of: To[3] =-1.0; To[7] = 1.0; To[2] =-1.0; is
        To[1] =-1.0; To[5] = 1.0; To[6] =-1.0;
    }
    else if (out_system == XYZ)
    {
        // Inverse of To[3] = 1.0; To[1] = 1.0; To[8] =-1.0; is:
        To[1] = 1.0; To[3] = 1.0; To[8] =-1.0;
    }
    else if (out_system == NWU)
    {
        //Inverse of: To[0] = 1.0; To[4] =-1.0; To[8] =-1.0; is
        To[0] = 1.0; To[4] =-1.0; To[8] =-1.0;
    }
    else
    {
        printf("%s: Invalid out coordinate system %d\n", fcnm, out_system);
        return -1;
    }
    // Multiply inv(in)*out
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                n3, n3, n3, alpha, From, n3, To, n3, beta,
                r9, n3);
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            r[i][j] = r9[3*j + i];
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Transform a matrix from in_sys to out_sys.  This is an interface
 *        routine to cmopad_basis__transformMatrixM33
 *        because many matrices are stored as a length 6 array.
 *
 * @param[in,out] m     on input length 6 matrix in in_sys.  the ordering is
 *                      \f$ {1,2,3,4,5,6} = {11,22,33,12,13,23} \f$.
 *                      on output length 6 matrix in out_sys.  the ordering is
 *                      \f$ {1,2,3,4,5,6} = {11,22,33,12,13,23} \f$.
 *
 * @param[in] in_sys    input system: NED, USE, XYZ, NWU
 * @param[in] out_sys   output system: NED, USE, XYZ, NWU
 *
 * @result 0 indicates success
 *
 */
int cmopad_basis_transformMatrixM6(double *m,
                                   enum cmopad_basis_enum in_sys,
                                   enum cmopad_basis_enum out_sys)
{
    const char *fcnm = "cmopad_basis_transformMatrixM6\0";
    double M33[3][3];
    int ierr;
    //------------------------------------------------------------------------//
    //
    // Convert array matrix to matrix matrix
    M33[0][0] = m[0];             //mxx; mrr
    M33[1][1] = m[1];             //myy; mtt
    M33[2][2] = m[2];             //mzz; mpp
    M33[0][1] = M33[1][0] = m[3]; //mxy; mrt
    M33[0][2] = M33[2][0] = m[4]; //mxz; mrp
    M33[1][2] = M33[2][1] = m[5]; //myz; mtp
    // Switch basis
    ierr = cmopad_basis_transformMatrixM33(M33, in_sys, out_sys);
    if (ierr != 0)
    {
        printf("%s: Error changing basis\n", fcnm);
        return -1;
    }
    // Convert matrix matrix back to array matrix
    m[0] = M33[0][0]; //mxx; mrr
    m[1] = M33[1][1]; //myy; mtt
    m[2] = M33[2][2]; //mzz; mpp
    m[3] = M33[0][1]; //mxy; mrt
    m[4] = M33[0][2]; //mxz; mrp
    m[5] = M33[1][2]; //myz; mtp
    return 0;
}
//============================================================================//
/*!
 * @brief Transform a symmetric 3x3 matrix from in_sys to out_sys.
 *
 * @param[in,out] m     on input [3x3] symmetric matrix in in_sys.
 *                      on output [3x3] symmetric matrix in out_sys
 *
 * @param[in] in_sys    input system: NED, USE, XYZ, NWU
 * @param[in] out_sys   output system: NED, USE, XYZ, NWU
 *
 * @result 0 indicates success
 *
 */
int cmopad_basis_transformMatrixM33(double m[3][3],
                                    enum cmopad_basis_enum in_sys,
                                    enum cmopad_basis_enum out_sys)
{
    const char *fcnm = "cmopad_basis_transformMatrixM33\0";
    double r[3][3], t9[9], invR[3][3], invR9[9], r9[9], m9[9];
    double alpha = 1.0;
    double beta = 0.0;
    int n3 = 3;
    int i, j, ierr;
    //------------------------------------------------------------------------//
    //
    // Compute change of basis matrix
    ierr = cmopad_basis_switcher(in_sys, out_sys, r);
    if (ierr != 0)
    {
        printf("%s: Error generating basis switching matrix\n",fcnm);
        return -1;
    }
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            r9[3*j + i] = r[i][j];
            m9[3*j + i] = m[i][j];
            //invR9[3*j + i] = r[i][j];
            invR[i][j] = r[i][j];
        }
    }
    // compute inv(R)
    ierr = __cmopad_inv3(invR); //double Amat[3][3])
    if (ierr != 0)
    {
        printf("%s: Unlikely error!\n", fcnm);
        return -1;
    }
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            invR9[3*j + i] = invR[i][j];
        }
    }
    // Compute R M inv(R)
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                n3, n3, n3, alpha, m9, n3, invR9, n3, beta,
                t9, n3);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                n3, n3, n3, alpha, r9, n3, t9, n3, beta,
                m9, n3);
    // Copy back
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            m[i][j] = m9[3*j + i];
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Transforms a vector v from basis in_sys to basis out_sys
 *
 * @param[in,out] v     on input length 3 vector in in_sys.
 *                      on output length 3 vector in out_sys.
 *
 * @param[in] in_sys    input system: NED, USE, XYZ, NWU
 * @param[in] out_sys   output system: NED, USE, XYZ, NWU
 *
 * @result 0 indicates success
 *
 */
int cmopad_basis_transformVector(double v[3],
                                 enum cmopad_basis_enum in_sys,
                                 enum cmopad_basis_enum out_sys)
{
    const char *fcnm = "cmopad_basis_transformVector\0";
    double r[3][3], r9[9], x[3];
    double alpha = 1.0;
    double beta = 0.0;
    int n3 = 3;
    int incx = 1;
    int incy = 1;
    int i,j,ierr;
    //------------------------------------------------------------------------//
    //
    // Compute change of basis matrix
    ierr = cmopad_basis_switcher(in_sys, out_sys, r);
    if (ierr != 0)
    {
        printf("%s: Error making basis switcher matrix\n", fcnm);
        return -1;
    }
    for (i=0; i<3; i++)
    {
        x[i] = v[i];
        for (j=0; j<3; j++)
        {
            r9[3*j + i] = r[i][j];
        }
    }
    // Compute v <- R*v
    cblas_dgemv(CblasColMajor, CblasNoTrans, n3, n3, alpha, r9,
                n3, x, incx, beta, v, incy);
    return 0;
}
//============================================================================//
/*!
 * @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]
 *
 */
static void array_cross3(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 Converts the unit eigenvector of length eig in coordinate system
 *        coord to length, azimuth, plunge.
 *        http://docs.obspy.org/_modules/obspy/imaging/beachball.html#PrincipalAxis
 *        This is an extension of MoPaD
 *
 * @param[in] coord    coordinate system (e.g., NED,USE,XYZ)
 * @param[in] eig      length of eigenvector
 * @param[in] ev       eigenvector in system coordinate system [3]
 *
 * @param[out] paxis   principal axis store azimuth, plunge, length [3]
 *
 * @result 0 indicates success
 *
 */
int cmopad_Eigenvector2PrincipalAxis(enum cmopad_basis_enum coord, double eig,
                                     double ev[3], double paxis[3])
{
    const char *fcnm = "cmopad_Eigenvector2PrincipalAxis\0";
    double v[3], az, plunge;
    double twopi = 2.0*M_PI;
    double pi180i = 180.0/M_PI;
    int ierr;
    //------------------------------------------------------------------------//
    //
    // Convert eigenvector ev to USE convention
    cblas_dcopy(3, ev, 1, v, 1);
    ierr = cmopad_basis_transformVector(v, coord, USE);
    if (ierr != 0)
    {
        printf("%s: Error switching basis\n", fcnm);
        return -1;
    }
    // Compute azimuth and plunge angles
    plunge = asin(-v[0]);       //arcsin(-z/r)
    az     = atan2(v[2],-v[1]); //atan(x/y)
    if (plunge <= 0.0)
    {
        plunge =-plunge;
        az = az + M_PI;
    }
    if (az < 0.0){az = az + twopi;}   //Shift to [0,360]
    if (az > twopi){az = az - twopi;} //Shift to [0,360]
    // Convert to degrees
    plunge = plunge*pi180i;
    az = az*pi180i;
    // Copy back result
    paxis[0] = az;     //First value is azimuth (degrees)
    paxis[1] = plunge; //Second value is plunge (degrees)
    paxis[2] = eig;    //Final value is eigenvalue (distance)
    return 0;
}
//============================================================================//
/*!
 * @brief Given coordinate system (x,y,z) and rotated system (xs,ys,zs)
 *        the line of nodes is the intersection between the x-y and the xs-ys
 *        planes.
 *
 * @param[in] alpha     is the angle between the z-axis and the zs-axis
 *                      and represents the dip angle (radians)
 * @param[in] beta      is the angle between the x-axis and the line of nodes
 *                      and represents the strike (radians)
 * @param[in] gamma     is the angle between the line of nodes and the xs-axis
 *                      and represents the rake angle (radians)
 *
 * @param[out] mat      rotation matrix built from euler angles [3 x 3]
 *
 * @note Original Python usage for moment tensors:
 *
 * m_unrot = numpy.matrix([[0,0,-1],[0,0,0],[-1,0,0]])
 * rotmat = cmopadEulerToMatrix(dip,strike,-rake)
 * m = rotmat.T * m_unrot * rotmat'''
 *
 */
void cmopad_eulerToMatrix(double alpha, double beta, double gamma,
                          double mat[3][3])
{
    double ca, cb, cg, sa, sb, sg;
    ca = cos(alpha);
    cb = cos(beta);
    cg = cos(gamma);
    sa = sin(alpha);
    sb = sin(beta);
    sg = sin(gamma);

    mat[0][0] = cb*cg-ca*sb*sg; mat[0][1] = sb*cg+ca*cb*sg; mat[0][2] = sa*sg;
    mat[1][0] =-cb*sg-ca*sb*cg; mat[1][1] =-sb*sg+ca*cb*cg; mat[1][2] = sa*cg;
    mat[2][0] = sa*sb;          mat[2][1] =-sa*cb;          mat[2][2] = ca;
    return;
}
//============================================================================//
/*!
 * @brief Computes the moment magnitude from a 6 vector
 *
 * @param[in] M6     moment tensor stored:
 *                    \f$ \{m_{11},m_{22},m_{33},m_{12},m_{13},m_{23} \} \f$
 *                   with the moment tensor terms of units Newton-meters
 *
 * @param[out] ierr  0 indicates success
 *
 * @result moment magnitue
 *
 * @author Ben Baker (ISTI)
 *
 */
double cmopad_momentMagnitudeM6(const double M6[6], int *ierr)
{
    const char *fcnm = "cmopad_momentMagnitudeM6\0";
    double mag;
    const double two_third = 2.0/3.0;
    *ierr = 0;
    mag = M6[0]*M6[0] + M6[1]*M6[1] + M6[2]*M6[2]
        + 2.0*(M6[3]*M6[3] + M6[4]*M6[4] + M6[5]*M6[5]);
    if (mag == 0.0)
    {
        *ierr = 1;
        printf("%s: Error magnitude is zero\n", fcnm);
    }
    else
    {
        mag = two_third*(log10(mag) - 9.1);
    }
    return mag;
}
//============================================================================//
/*!
 * @brief Computes the moment magnitude from a the symmetric moment tensor
 *
 * @param[in] M      symmetric moment tensor with terms given in Newton-meters
 *
 * @param[out] ierr  0 indicates success
 *
 * @result moment magnitue
 *
 * @author Ben Baker (ISTI)
 *
 */
double cmopad_momentMagnitudeM3x3(const double M[3][3], int *ierr)
{
    const char *fcnm = "cmopad_momentMagnitudeM33\0";
    double M6[6], mag;
    M6[0] = M[0][0];
    M6[1] = M[1][1];
    M6[2] = M[2][2];
    M6[3] = M[0][1];
    M6[4] = M[0][2];
    M6[5] = M[1][2];
    mag = cmopad_momentMagnitudeM6(M6, ierr);
    if (*ierr != 0)
    {
        printf("%s: Error computing moment magnitude\n", fcnm);
    }
    return mag;
}
//============================================================================//
/*!
 * @brief Sets the two angle-triples, describing the faultplanes of the
 *        Double Couple, defined by the eigenvectors P and T of the
 *        moment tensor object.
 *
 *        Define a reference Double Couple with strike = dip =
 *        slip-rake = 0, the moment tensor object's DC is transformed
 *        (rotated) w.r.t. this orientation. The respective rotation
 *        matrix yields the first fault plane angles as the Euler
 *        angles. After flipping the first reference plane by
 *        multiplying the appropriate flip-matrix, one gets the second fault
 *        plane's geometry.
 *
 *        All output angles are in degree
 *
 * @param[in] iverb     controls verbosity for debugging. 0 should be quiet.
 *
 * @param[in,out] src   on input holds plot_clr_order
 *                      on output holds the strike, dips, and rakes
 *                      describing fault planes fp1 and fp2
 *
 * @result 0 indicates success
 *
 */
int cmopad_findFaultPlanes(int iverb, struct cmopad_struct *src)
{
    const char *fcnm = "cmopad_findFaultPlanes\0";
    int n3 = 3;
    double rot_matrix_fp1[3][3], rot_matrix_fp2[3][3],
           refDC_evecs[9], flip_dc[9], work1[9], work2[9],
           pnt_sorted_EV_matrix[3][3],
           det;
    double sqrt22 = 0.7071067811865476; //sqrt(2)/2
    double alpha = 1.0, beta = 0.0;
    int i, j;
    //------------------------------------------------------------------------//
    //
    // reference Double Couple (in NED) basis w/ (strike,dip,slip) = (0,0,0)
    refDC_evecs[0] =-sqrt22; refDC_evecs[3] = 0.0; refDC_evecs[6] =-sqrt22;
    refDC_evecs[1] = 0.0;    refDC_evecs[4] = 1.0; refDC_evecs[7] = 0.0;
    refDC_evecs[2] =-sqrt22; refDC_evecs[5] = 0.0; refDC_evecs[8] = sqrt22;

    //refDC_evals[0] =-1.0;
    //refDC_evals[1] = 0.0;
    //refDC_evals[2] = 1.0;
    // Matrix which is turning from one fault plane to the other
    flip_dc[0] = 0.0; flip_dc[3] = 0.0; flip_dc[6] =-1.0;
    flip_dc[1] = 0.0; flip_dc[4] =-1.0; flip_dc[7] = 0.0;
    flip_dc[2] =-1.0; flip_dc[5] = 0.0; flip_dc[8] = 0.0;
    // Euler-tools need matrices of EV sorted in PNT (pressure,null,tension):
    for (i=0; i<3; i++)
    {
       for (j=0; j<3; j++)
       {
          pnt_sorted_EV_matrix[i][j] = src->rotation_matrix[i][j];
       }
    }
    // Re-sort only necessary if abs(p) <= abs(t)
    if (src->plot_clr_order < 0)
    {
        for (i=0; i<3; i++)
        {
            pnt_sorted_EV_matrix[i][0] = src->rotation_matrix[i][2];
            pnt_sorted_EV_matrix[i][2] = src->rotation_matrix[i][0];
        }
    }
    // Rotation matrix, describing the rotation of the eigenvector
    // system of the input moment tensor into the eigenvector
    // system of the reference Double Couple
    if (iverb > 1)
    {
        printf("%s: Pressure, null, tension eigenvectors:\n",fcnm);
    }
    for (j=0; j<3; j++)
    {
       for (i=0; i<3; i++)
       {
          work1[j*n3 + i] = pnt_sorted_EV_matrix[i][j];
       }
       if (iverb > 1)
       {
           printf(
           "%8.5f %8.5f %8.5f\n",pnt_sorted_EV_matrix[j][0],
                                 pnt_sorted_EV_matrix[j][1],
                                 pnt_sorted_EV_matrix[j][2]);
       }
    }
    if (iverb > 1)
    {
        printf("%s: Reference double couple eigenvectors:\n",fcnm);
        for (i=0; i<3; i++)
        {
            printf(
            "%8.5f %8.5f %8.5f\n",refDC_evecs[i],
                                  refDC_evecs[3+i],refDC_evecs[6+i]);
        }
    }
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n3, n3, n3,
                alpha, work1, n3, refDC_evecs, n3, beta,
                work2, n3);
    for (j=0; j<3; j++)
    {
        for (i=0; i<3; i++)
        {
            rot_matrix_fp1[j][i] = work2[j*n3 + i]; //Copy transpose
        }
    }
    // check if rotation has correct orientation
    det = __cmopad_determinant3x3(rot_matrix_fp1);
    if (det < 0.0)
    {
       for (i=0; i<3; i++)
          {
          for (j=0; j<3; j++)
          {
              rot_matrix_fp1[i][j] =-1.0*rot_matrix_fp1[i][j];
          }
       }
    }
    for (j=0; j<3; j++)
    {
        for (i=0; i<3; i++)
        {
            work1[j*3 + i] = rot_matrix_fp1[i][j];
        }
    }
    // adding a rotation into the (ambiguous) system of the 2nd fault plane
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n3, n3, n3,
                alpha, flip_dc, n3, work1, n3, beta,
                work2, n3);
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            rot_matrix_fp2[i][j] = work2[j*3 + i];
        }
    }
    // Calculate strike, dip, rake of fault planes
    cmopad_findStrikeDipRake(rot_matrix_fp1,
                             &src->fp1[1], &src->fp1[0], &src->fp1[2]);
    cmopad_findStrikeDipRake(rot_matrix_fp2,
                             &src->fp2[1], &src->fp2[0], &src->fp2[2]);
    // For convenience if choose the first fault plane to be the one with
    // the smaller strike
    if (src->fp2[0] < src->fp1[0])
    {
        for (i=0; i<3; i++)
        {
            __cmopad_swap8(&src->fp1[i], &src->fp2[i]);
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Returns dip, strike, and slip(rake) angles in degrees describing the
 *        fault plane.
 *
 * @param[in] rot_mat     matrix of eignevectors for fault plane [3 x 3]
 *
 * @param[out] alpha      dip angle (degrees)
 * @param[out] beta       strike angle (degrees)
 * @param[out] gamma      rake angle (degrees)
 *
 */
void cmopad_findStrikeDipRake(double rot_mat[3][3],
                              double *alpha, double *beta, double *gamma)
{
    double a = 0.0, b = 0.0, g = 0.0;
    cmopad_MatrixToEuler(rot_mat, &a, &b, &g);
    *alpha = a*rad2deg;
    *beta  = b*rad2deg;
    *gamma =-g*rad2deg;
    return;
}
//============================================================================//
/*!
 * @brief Returns three Euler angles alpha, beta, and gamma (radians) from a
 *        rotation matrix.
 *
 * @param[in] rotmat    rotation matrix of which to decompose into euler angles
 *                      [3 x 3]
 *
 * @param[out] alpha    euler angle representing dip (radians)
 * @param[out] beta     euler angle representing strike (radians)
 * @param[out] gamma    euler angle representing -rake (radians)
 *
 */
void cmopad_MatrixToEuler(double rotmat[3][3],
                          double *alpha, double *beta, double *gamma)
{
    const int n3 = 3;
    int inc = 1;
    double ex[3], ez[3], exs[3], ezs[3], enodes[3], enodess[3],
           xnorm, cos_alpha, a,b,g,a1,a2;
    double twopi = 2.0*M_PI;
    int i,j;
    //------------------------------------------------------------------------//
    //
    // transpose matrix vector multiply
    ex[0] = 1.0; ex[1] = 0.0; ex[2] = 0.0;
    ez[0] = 0.0; ez[1] = 0.0; ez[2] = 1.0;
    for (i=0; i<3; i++)
    {
        exs[i] = 0.0; ezs[i] = 0.0;
        for (j=0; j<3; j++)
        {
            exs[i] = exs[i] + rotmat[j][i]*ex[j];
            ezs[i] = ezs[i] + rotmat[j][i]*ez[j];
        }
    }
    // Cross product
    array_cross3(ez, ezs, enodes);
    xnorm = sqrt(pow(enodes[0],2) + pow(enodes[1],2) + pow(enodes[2],2));
    if (xnorm < 1.e-10)
    {
        for (i=0; i<3; i++)
        {
            enodes[i] = exs[i];
        }
    }
    for (i=0; i<3; i++)
    {
        enodess[i] = 0.0;
        for (j=0; j<3; j++)
        {
            enodess[i] = enodess[i] + rotmat[i][j]*enodes[j];
        }
    }
    // Multiply
    cos_alpha = cblas_ddot(n3, ez, inc, ezs, inc);
    cos_alpha = fmin( 1.0, cos_alpha);
    cos_alpha = fmax(-1.0, cos_alpha);
    a = acos(cos_alpha);
    a1 = atan2( enodes[1],  enodes[0]);
    a2 =-atan2(enodess[1], enodess[0]);
    b = numpy_mod(a1, twopi);
    g = numpy_mod(a2, twopi);
    // Compute unique Euler angles
    cmopad_uniqueEuler(&a, &b, &g);
    *alpha = a;
    *beta = b;
    *gamma = g;
    return;
}
//============================================================================//
/*!
 * @brief Clone of the numpy mod(a,b) function
 *
 * @result numpy.mod(a, b)
 *
 */
static double numpy_mod(double a, double b)
{
    if (b == 0.0){return a;}
    if (a >= 0.0)
    {
        if (b > 0.0)
        {
            return a - floor(a/b)*b;
        }
        else
        {
            return a - fabs(floor(a/b))*fabs(b);
        }
    }
    else
    {
        if (b > 0.0)
        {
            return a + fabs(floor(a/b))*fabs(b);
        }
        else
        {
            return a + floor(a/b)*fabs(b);
        }
    }
}
//============================================================================//
/*!
 * Read in Matrix M and set up eigenvalues (EW) and eigenvectors
 * (EV) for setting up the principal axis system.
 *
 * The internal convention is the 'HNS'-system: H is the
 * eigenvector for the smallest absolute eigenvalue, S is the
 * eigenvector for the largest absolute eigenvalue, N is the null
 * axis.
 *
 * Naming due to the geometry: a CLVD is
 * Symmetric to the S-axis,
 * Null-axis is common sense, and the third (auxiliary) axis
 * Helps to construct the RA3.
 *
 * Additionally builds matrix for basis transformation back to NED system.
 *
 * The eigensystem setup defines the colouring order for a later
 * plotting in the BeachBall class. This order is set by the
 * `_plot_clr_order' attribute.
 *
 * @param[in] iverb     controls verbosity.  0 is quiet.
 *
 * @param[in,out] src   on input holds the input moment tensor
 *                      on output has the eigenvalues and fault planes
 *
 */
int cmopad_MT2PrincipalAxisSystem(int iverb, struct cmopad_struct *src)
{
    const char *fcnm = "cmopad_MT2PrincipalAxisSystem\0";
    int ival = 0;
    double M[3][3], M_devi[3][3], EV_devi[3][3], EV[3][3], EW[3],
           EW_devi[3], EV1[3], EV2[3], EV3[3], EVh[3], EVs[3], EVn[3],
           EWs, EWh, EWn,
           EW1, EW2, EW3, EW1_devi, EW2_devi, EW3_devi,
           trace_M, xsum;
    int EW_order[3], symmetry_around_tension, clr, i, j, ierr;
    //------------------------------------------------------------------------//
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            M[i][j]      = src->M[i][j];
            M_devi[i][j] = src->M_devi[i][j];
            EV[i][j] = M[i][j];
            EV_devi[i][j] = M_devi[i][j];
        }
    }
    ierr = __cmopad_Eigs3x3(EV_devi, ival, EW_devi);
    if (ierr != 0)
    {
        printf("%s: Error computing eigenvalues (a)!", fcnm);
        return -1;
    }
    ierr = __cmopad_argsort3(EW_devi, EW_order);
    // Removed if statement, always true in python code
    trace_M = __cmopad_trace3(M); // M[0][0] + M[1][1] + M[2][2];
    if (trace_M < epsilon){trace_M = 0.0;}
    ierr = __cmopad_Eigs3x3(EV, ival, EW);
    if (ierr != 0)
    {
        printf("%s: Error computing eigenvalues (b)!", fcnm);
        return -1;
    }
    for (i=0; i<3; i++)
    {
        if (fabs(EW[i]) < epsilon){EW[i] = 0.0;}
    }
    //trace_M_devi = cmopadTrace(3,M_devi);

    // Save eigenvalues/eigenvectors for deviatoric and full moment tensor
    EW1_devi = EW_devi[EW_order[0]];
    EW2_devi = EW_devi[EW_order[1]];
    EW3_devi = EW_devi[EW_order[2]];
    if (fabs(EW2_devi) < epsilon){EW2_devi = 0.0;}

    EW1 = EW[EW_order[0]];
    EW2 = EW[EW_order[1]];
    EW3 = EW[EW_order[2]];
    for (i=0; i<3; i++)
    {
        //EV1_devi[i] = EV_devi[i][EW_order[0]];
        //EV2_devi[i] = EV_devi[i][EW_order[1]];
        //EV3_devi[i] = EV_devi[i][EW_order[2]];
        EV1[i] = EV[i][EW_order[0]];
        EV2[i] = EV[i][EW_order[1]];
        EV3[i] = EV[i][EW_order[2]];
    }

    // Classify tension axis symmetry
    symmetry_around_tension = 1;
    clr = 1;
    xsum = EW1 + EW2 + EW3;
    // implosion
    if (EW1 < 0.0 && EW2 < 0.0 && EW3 < 0.0)
    {
        if (iverb > 2){printf("%s: Implosion\n", fcnm);}
        symmetry_around_tension = 0;
        clr = 1;
    }
    // explosion
    else if (EW1 > 0.0 && EW2 > 0.0 && EW3 > 0.0)
    {
        if (iverb > 2){printf("%s: Explosion\n", fcnm);}
        symmetry_around_tension = 1;
        if (fabs(EW1_devi) > fabs(EW3_devi))
        {
            symmetry_around_tension = 0;
        }
        clr = -1;
    }
    // net-implosion
    else if (EW2 < 0.0 && xsum < 0.0)
    {
        if (iverb > 2){printf("%s: Net-implosion (1)\n", fcnm);}
        if (fabs(EW1_devi) < fabs(EW3_devi))
        {
            symmetry_around_tension = 1;
            clr = 1;
        }
        else
        {
            symmetry_around_tension = 1;
            clr = 1;
        }
    }
    // net-implosion
    else if (EW2_devi >= 0.0 && xsum < 0.0)
    {
        if (iverb > 2){printf("%s: Net-impolosion (2)\n", fcnm);}
        symmetry_around_tension = 0;
        clr = -1;
        if (fabs(EW1_devi) < fabs(EW3_devi))
        {
            symmetry_around_tension = 1;
            clr = 1;
        }
    }
    // net-explosion
    else if (EW2_devi < 0.0 && xsum > 0.0)
    {
        if (iverb > 2){printf("%s: Net-explosion (1)\n", fcnm);}
        symmetry_around_tension = 1;
        clr = 1;
        if (fabs(EW1_devi) > fabs(EW3_devi))
        {
            symmetry_around_tension = 0;
            clr = -1;
        }
    }
    // net-explosion
    else if ( EW2_devi >= 0.0 && xsum > 0.0)
    {
        if (iverb > 2){printf("%s: Net-explosion (2)\n", fcnm);}
        symmetry_around_tension = 0;
        clr = -1;
    }
    else
    {
        // If the trace is 0 to machine epsilon then purely deviatoric
        if (trace_M/fmax(fabs(EW1_devi),fabs(EW3_devi)) > 1.e-10){ //!= 0.0){
            printf("%s: Warning should not be here %f %f %f\n",
                      fcnm, EW2_devi, trace_M, xsum);
        }
    }
    if (iverb > 1)
    {
        printf("%s: sym/clr %d %d\n",
                  fcnm, symmetry_around_tension, clr);
    }
    // Deviatoric classification
    if (fabs(EW1_devi) < fabs(EW3_devi))
    {
        symmetry_around_tension = 1;
        clr = 1;
    }
    if (fabs(EW1_devi) >= fabs(EW3_devi))
    {
        symmetry_around_tension = 0;
        clr = -1;
    }
    //explosive src with all neg evals? - TODO: email lars.  an edge
    //case is (-1,-1,-1, 0,0,0) which has trace_M == 0.0 and is a pure
    //implosion.  in this case there is a breakdown of the code because
    //trace_M == 0.0
    if (EW3 < 0.0 && trace_M > 0.0)
    {
        printf("%s: Critical error!\n", fcnm);
        return -1;
    }

    // Classify pure deviatoric - pure deviatoric
    if (trace_M == 0.0)
    {
        if (iverb > 1){printf("%s: Pure-deviatoric\n", fcnm);}
        // pure shear
        if (EW2 == 0.0)
        {
            if (iverb > 1){printf("%s: Pure-shear\n", fcnm);}
            symmetry_around_tension = 1;
            clr = 1;
        }
        else if (2.0*fabs(EW2) == fabs(EW1) || 2.0*fabs(EW2) == fabs(EW3))
        {
            if (iverb > 1){printf("%s: Pure CLVD\n", fcnm);}
            // CLVD symmetry around tension
            if (fabs(EW1) < EW3)
            {
                if (iverb > 2)
                {
                    printf("%s: CLVD symmetry around tension\n", fcnm);
                }
                symmetry_around_tension = 1;
                clr = 1;
            }
            // CLVD symmetry around pressure
            else
            {
                if (iverb > 2)
                {
                    printf("%s: CLVD symmetry around pressure\n", fcnm);
                }
                symmetry_around_tension = 0;
                clr = -1;
            }
        }
        // mix of DC and CLVD
        else
        {
            if (iverb > 1)
            {
                printf("%s: Mix of DC and CLVD\n", fcnm);
            }
            // symmetry around tension
            if (fabs(EW1) < EW3)
            {
                if (iverb > 2)
                {
                    printf("%s: Symmetry around tension\n", fcnm);
                }
                symmetry_around_tension = 1;
                clr = 1;
            }
            // symmetry around pressure
            else
            {
                if (iverb > 2)
                {
                    printf("%s: Symmetry around pressure\n", fcnm);
                }
                symmetry_around_tension = 0;
                clr = -1;
            }
        }
    }
    if (iverb > 0)
    {
        printf("%s: Final symmetry %d %d\n",
                  fcnm, symmetry_around_tension, clr);
    }
    // Define order of eigenvectors and values according to symmetry axis
    if (symmetry_around_tension == 1)
    {
        EWs = EW3;
        EWh = EW1;
        EWn = EW2;
        for (i=0; i<3; i++)
        {
            EVs[i] = EV3[i];
            EVh[i] = EV1[i];
            EVn[i] = EV2[i];
        }
    }
    else
    {
        EWs = EW1;
        EWh = EW3;
        EWn = EW2;
        for (i=0; i<3; i++)
        {
            EVs[i] = EV1[i];
            EVh[i] = EV3[i];
            EVn[i] = EV2[i];
        }
    }
    // Order of eigenvector's basis: (H,N,S)
    for (i=0; i<3; i++)
    {
        // Matrix for basis transform
        src->rotation_matrix[i][0] = EVh[i];
        src->rotation_matrix[i][1] = EVn[i];
        src->rotation_matrix[i][2] = EVs[i];
        // Save eigenvectors
        src->eigenvectors[i][0] = EVh[i];
        src->eigenvectors[i][1] = EVn[i];
        src->eigenvectors[i][2] = EVs[i];
        // Principal axis
        src->p_axis[i] = EV1[i];    //Smallest (most negative) eval/evec pair
        src->null_axis[i] = EVn[i]; //Intermediate eval/evec pair
        src->t_axis[i] = EV3[i];    //Largest (most positive) eval/evec pair
        /* TODO the above is definitionally correct but the following works.
           What's going on here?  Did someone switch conventions on me? */
        //src->t_axis[i] = EV1[i];
        //src->null_axis[i] = EVn[i];
        //src->p_axis[i] = EV3[i];
    }
    src->eig_pnt[0] = EW1_devi; // Smallest - most negative eigenvalue
    src->eig_pnt[1] = EW2_devi;
    src->eig_pnt[2] = EW3_devi; // Largest eigenvalue
    // TODO - these are reversed to match above t_axis, null_axis, and p_axis
    //printf("%f %f\n", EW1_devi, EW3_devi);
    //src->eig_pnt[0] = EW3_devi; // Make pressure most positive eigenvalue?
    //src->eig_pnt[1] = EW2_devi;
    //src->eig_pnt[2] = EW1_devi; // Make tension most negative eigenvalue
    // Eigenvalues corresponding to eigenvectors (H,N,S) basis
    src->eigenvalues[0] = EWh;
    src->eigenvalues[1] = EWn;
    src->eigenvalues[2] = EWs;

    // Important for plot in BeachBall class (also used in FindFaultPlanes)
    src->plot_clr_order = clr;

    ierr = cmopad_findFaultPlanes(iverb, src);
    if (ierr != 0)
    {
        printf("%s: Error computing fault planes!\n", fcnm);
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Prints a 3x3 matrix in a Mopad-esque style.
 *
 * @param[in] lfact   if true then use a normalization factor
 * @param[in] m       the 3x3 matrix to write to standard out
 *
 */
void cmopad_printMatrix3x3(bool lfact, double m[3][3])
{
    const char *fcnm = "cmopad_printMatrix3x3\0";
    double norm_fact, m11,m22,m33,m12,m13,m23;
    int i,j;
    //------------------------------------------------------------------------//
    norm_fact = 1.0;
    if (lfact)
    {
        norm_fact = 0.0;
        for (i=0; i<3; i++)
        {
            for (j=0; j<3; j++)
            {
                //norm_fact = fmax(norm_fact,fabs(mathRound(m[i][j],5)));
                norm_fact = fabs(m[i][j]);
            }
        }
        if (norm_fact == 0.0)
        {
            printf("%s: Warning norm_fact = 0.0!\n", fcnm);
            return;
        }
    }
    m11 = m[0][0]/norm_fact;
    m22 = m[1][1]/norm_fact;
    m33 = m[2][2]/norm_fact;
    m12 = m[0][1]/norm_fact;
    m13 = m[0][2]/norm_fact;
    m23 = m[1][2]/norm_fact;
    if (lfact)
    {
        printf(" [ %5.2f %5.2f %5.2f ] \n",m11,m12,m13);
        printf(" [ %5.2f %5.2f %5.2f ] x %f\n",
                          m12,m22,m23,norm_fact);
        printf(" [ %5.2f %5.2f %5.2f ] \n",m13,m23,m33);
    }
    else
    {
        printf(" [ %5.2f %5.2f %5.2f ]\n",m11,m12,m13);
        printf(" [ %5.2f %5.2f %5.2f ]\n",m12,m22,m23);
        printf(" [ %5.2f %5.2f %5.2f ]\n",m13,m23,m33);
    }
    return;
}
//============================================================================//
/*!
 * @brief Brings the provided mechanism into symmetric 3x3 matrix form.
 *
 *        The source mechanism may be provided in different forms:
 *
 *        -- as 3x3 matrix - symmetry is checked - one basis
 *           system has to be chosen, or NED as default is taken
 *        -- as 3-element tuple or array - interpreted as
 *           strike, dip, slip-rake angles in degree
 *        -- as 4-element tuple or array - interpreted as strike, dip,
 *           slip-rake angles in degree + seismic scalar moment in Nm
 *        -- as 6-element tuple or array - interpreted as the 6 independent
 *           entries of the moment tensor
 *        -- as 7-element tuple or array - interpreted as the 6 independent
 *           entries of the moment tensor + seismic scalar moment in Nm
 *        -- as 9-element tuple or array - interpreted as the 9 entries
 *           of the moment tensor - checked for symmetry
 *        -- as a nesting of one of the upper types (e.g. a list of
 *           n-tuples) - first element of outer nesting is taken
 *
 */
int cmopad_SetupMT(int nmech, double *mech,
                   enum cmopad_basis_enum input_basisIn,
                   double Mech_out[3][3])
{
    const char *fcnm = "cmopad_SetupMT\0";
    enum cmopad_basis_enum input_basis;
    double rotmat[9], mtemp1[9], mtemp2[9], m_unrot[9], scalar_moment,
           strike, dip, rake;
    int i, j, ierr;
    double pi180 = M_PI/180.0;
    double alpha = 1.0;
    double beta = 0.0;
    int n3 = 3;
    //------------------------------------------------------------------------//
    //
    // All 9 elements are specified
    input_basis = input_basisIn;
    if (nmech == 9)
    {
        for (i=0; i<3; i++)
        {
            for (j=0; j<3; j++)
            {
                Mech_out[i][j] = mech[3*j + i];
            }
        }
    }
    // Mechanism given as 6 or 7 element array
    else if (nmech == 6 || nmech == 7)
    {
        scalar_moment = 1.0;
        if (nmech == 7){scalar_moment = mech[6];}
        Mech_out[0][0] = scalar_moment*mech[0];
        Mech_out[1][1] = scalar_moment*mech[1];
        Mech_out[2][2] = scalar_moment*mech[2];
        Mech_out[0][1] = Mech_out[1][0] = scalar_moment*mech[3];
        Mech_out[0][2] = Mech_out[2][0] = scalar_moment*mech[4];
        Mech_out[1][2] = Mech_out[2][1] = scalar_moment*mech[5];
    }
    // Strike dip rake; conventions from Jost and Herrmann; NED-basis
    else if (nmech == 3 || nmech == 4)
    {
        scalar_moment = 1.0;
        if (nmech == 4){scalar_moment = mech[3];}
        // Set matrix as function of Euler angles
        strike = mech[0]*pi180;
        dip    = mech[1]*pi180;
        rake   = mech[2]*pi180;
        cmopad_eulerToMatrix(dip, strike, -rake, Mech_out);
        for (i=0; i<3; i++)
        {
            for (j=0; j<3; j++)
            {
                m_unrot[3*j + i] = 0.0;
                rotmat[3*j + i] = Mech_out[i][j];
                mtemp1[3*j + i] = 0.0;
                mtemp2[3*j + i] = 0.0;
            }
        }
        m_unrot[6] =-1.0;
        m_unrot[2] =-1.0;
        // Multiply trans(Rotmat)*m_unrot*Rotmat
        cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, n3, n3, n3,
                    alpha, rotmat, n3, m_unrot, n3, beta,
                    mtemp1, n3);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n3, n3, n3,
                    alpha, rotmat, n3, mtemp1, n3, beta,
                    mtemp2, n3);
        // Copy back and include scalar moment
        for (i=0; i<3; i++)
        {
            for (j=0; j<3; j++)
            {
                Mech_out[i][j] = scalar_moment*mtemp2[3*j + i];
            }
        }
        // Assure right basis system
        input_basis = NED;
    }
    // Unknown mechanism
    else
    {
        printf("%s: Invalid mechanism!\n", fcnm);
        return -1;
    }
    // Rotate into NED basis for local processing
    ierr = cmopad_basis_transformMatrixM33(Mech_out, input_basis, NED);
    if (ierr != 0)
    {
        printf("%s: Error transforming matrix\n", fcnm);
        return -1;
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Sorts eigenvalues/eigenvectors of a matrix into ascending order.
 *        Optionally, can sort based on the absolute magnitude of eigenvalues
 *
 * @param[in] isabs     if true then sort on absolute value of eigenvalues
 *
 * @param[in,out] e     eigenvalues sorted in ascending order [3]
 * @param[in,out] ev    corresponding eigenvectors sorted along with eigenvalues
 *                      [3 x 3]
 *
 * @result 0 indicates success
 *
 */
int cmopad_sortEigenvAscending(bool isabs, double e[3], double ev[3][3])
{
    const char *fcnm = "cmopad_sortEigenvAscending\0";
    const int n = 3;
    double evwork[3][3], ework[3], esave[3];
    int iperm[3], i, ierr, j;
    ierr = 0;
    for (i=0; i<n; i++)
    {
        esave[i] = e[i];
        if (isabs)
        {
            ework[i] = fabs(e[i]);
        }
        else
        {
            ework[i] = e[i];
        }
        for (j=0; j<n; j++)
        {
            evwork[i][j] = ev[i][j];
        }
    }
    // Sort
    ierr = __cmopad_argsort3(ework, iperm);
    if (ierr != 0)
    {
        printf("%s: Error in permutation sort\n", fcnm);
        return -1;
    }
    // Copy back into ascending order
    for (i=0; i<n; i++)
    {
        e[i] = esave[iperm[i]];
        for (j=0; j<n; j++)
        {
            ev[i][j] = evwork[i][iperm[j]];
        }
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Decomposes the input moment tensor Min in the North, East, Down
 *        coordinates into an isotropic, compensated linear vector dipole,
 *        and pure double couple.
 *
 * @param[in] Min   3 x 3 moment tensor in NED coordinates to decompose
 *                  into a pure isotropic, CLVD, and DC
 *
 * @param[out] src  internal mopad structure.  this routine will compute:
 *                  M:      copy of the original moment tensor
 *                  M_iso:  diagonal isotropic moment tensor
 *                  M_devi: deviatoric moment tensor (M - M_iso)
 *                  M_CLVD: CLVD contribution to M_devi
 *                  M_DC:   DC contribution to M_devi
 *                  DC_percentage:   Percent (0,100) double couple
 *                  CLVD_percentage: Percent (0,100) CLVD
 *                  ISO_Percetnage:  Percent (0,100) isotropic
 *                  seismic_moment:  seismic moment (Nm) Bowers & Hudson
 *                                   1999
 *                  moment_magnitude: moment magnitude Hanks and Kanamori
 *                                    1979
 *
 * @result 0 indicates success
 *
 */
int cmopad_standardDecomposition(double Min[3][3], struct cmopad_struct *src)
{
    const char *fcnm = "cmopad_standardDecomposition\0";
    double M[3][3], M_devi[3][3], M_iso[3][3],
           M_DC[3][3], M_CLVD[3][3], eigenvtot[3][3], eigenvdevi[3][3],
           a3out, a2out, a1out, eigenwtot[3], eigenwdevi[3],
           a1[3], a2[3], a3[3];
    double M0, M0_iso, M0_devi, F, M_DC_percentage, M_iso_percentage, tracem;
    int lsub, i,j, ierr;
    int ival = 0; //Calculate eigenvalues only (1)
    // Copy
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            M[i][j] = Min[i][j];
        }
    }
    // Isotropic part
    tracem = __cmopad_trace3(M); //M[0][0] + M[1][1] + M[2][2];
    M_iso[0][0] = M_iso[0][1] = M_iso[0][2] = 0.0;
    M_iso[1][0] = M_iso[1][1] = M_iso[1][2] = 0.0;
    M_iso[2][0] = M_iso[2][1] = M_iso[2][2] = 0.0;
    M_iso[0][0] = M_iso[1][1] = M_iso[2][2] = (1.0/3.0)*tracem;
    M0_iso = fabs(1.0/3.0*tracem);

    // Deviatoric part
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            M_devi[i][j] = M[i][j] - M_iso[i][j];
        }
    }
    // Copy to structure
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            src->M[i][j] = M[i][j];
            src->M_iso[i][j]  = M_iso[i][j];
            src->M_devi[i][j] = M_devi[i][j];
        }
    }

    // Eigenvalues and eigenvectors
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            eigenvtot[i][j] = M_devi[i][j]; //I dont understand this
            eigenvdevi[i][j]= M_devi[i][j];
        }
    }
    ierr = __cmopad_Eigs3x3(eigenvtot, ival, eigenwtot);
    if (ierr != 0)
    {
        printf("%s: Error computing eigs (a)!\n",fcnm);
        return -1;
    }
    ierr = __cmopad_Eigs3x3(eigenvdevi, ival, eigenwdevi);
    if (ierr != 0)
    {
        printf("%s: Error computing eigs (b)!\n",fcnm);
        return -1;
    }

    // Eigenvalues in ascending order
    ierr = ierr + cmopad_sortEigenvAscending(true, eigenwtot , eigenvtot);
    ierr = ierr + cmopad_sortEigenvAscending(true, eigenwdevi, eigenvdevi);
    if (ierr != 0)
    {
        printf("%s: Error sorting eigenvalues/eigenvectors\n", fcnm);
        return -1;
    }

    /* TODO: Email Lars Krieger; confirm this is a mistake
           : Jost and Herrmann Eqn 19
    M0_devi =-1.0;
    for (i=0; i<3; i++){
        M0_devi = fmax(M0_devi, fabs(eigenwdevi[i]));
    }
    */
    // Compute mangnitude (Jost and Herrmann Eqn 19 noting the eigenvalues are
    // in ascending order)
    M0_devi = 0.5*(fabs(eigenwdevi[1]) + fabs(eigenwdevi[2]));

    // Named according to Jost and Herrmann
    for (i=0; i<3; i++)
    {
        a1[i] = eigenvtot[i][0];
        a2[i] = eigenvtot[i][1];
        a3[i] = eigenvtot[i][2];
    }

    // If only isotropic part exists
    if (M0_devi < epsilon)
    {
        F = 0.5;
    }
    else
    {
        F =-eigenwdevi[0]/eigenwdevi[2];
    }

    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            a3out = a3[i]*a3[j];
            a2out = a2[i]*a2[j];
            a1out = a1[i]*a1[j];
            // Second term of eqn 37 Jost and Herrmann
            M_DC[i][j] = eigenwdevi[2]*(1.0 - 2.0*F)*(a3out - a2out);
            // Remainder of deviatoric is CLVD (could do 3rd term of eqn 37)
            lsub = 1;
            if (lsub == 1)
            {
                M_CLVD[i][j] = M_devi[i][j] - M_DC[i][j];
            }
            else
            {
                M_CLVD[i][j] = eigenwdevi[2]*F
                              *(2.0*a3out - a2out - a1out);
            }
        }
    }

    //according to Bowers & Hudson:
    M0 = M0_iso + M0_devi;

    // Really shouldn't be rounding
    /*
    M_iso_percentage     = mathRound(M0_iso/M0*100.0,6);
    //cmopad.iso_percentage = M_iso_percentage;
    M_DC_percentage = ( mathRound( (1.0 - 2.0*fabs(F))
                                  *(1.0 - M_iso_percentage/100.0)*100,6));
    */
    M_iso_percentage = M0_iso/M0*100.0;
    M_DC_percentage = (1.0 - 2.0*fabs(F))*(1.0 - M_iso_percentage/100.0)*100.0;
    for (i=0; i<3; i++)
    {
        for (j=0; j<3; j++)
        {
            src->DC[i][j]   =  M_DC[i][j];
            src->CLVD[i][j] =  M_CLVD[i][j];
        }
    }
    src->DC_percentage   = M_DC_percentage;
    src->ISO_percentage  = M_iso_percentage;
    src->DEV_percentage  = 100.0 - src->ISO_percentage;
    src->CLVD_percentage = 100.0 - src->ISO_percentage - src->DC_percentage;
    ////src->seismic_moment   = sqrt(1./2.0*(pow(eigenw[0],2) + pow(eigenw[0],2) + pow(eigenw[0],2)));
    src->seismic_moment   = M0;
    src->moment_magnitude = log10(src->seismic_moment*1.0e7)/1.5 - 16.1/1.5;

    return 0;
}
//============================================================================//
/*!
 * @brief Return 6-tuple containing entries of M, calculated from fault plane
 *        angles (defined as in Jost&Herman), given in degrees.
 *        Basis for output is NED (= X,Y,Z)
 *
 * @param[in] strike    angle clockwise between north and plane (in [0,360])
 * @param[in] dip       angle between surface and dipping plane (in [0,90])
 *                      0 = horizontal, 90 = vertical
 * @param[in] rake      angle on the rupture plane between strike vector
 *                      and actual movement (defined mathematically
 *                      counter-clockwise rotation is positive)
 *
 * @param[out] moments  M = {M_nn, M_ee, M_dd, M_ne, M_nd, M_ed}
 *
 */
void cmopad_strikeDipRake2MT6(double strike, double dip, double rake,
                              double moments[6])
{
    double M1,M2,M3,M4,M5,M6;
    double S_rad = strike/rad2deg;
    double D_rad = dip/rad2deg;
    double R_rad = rake/rad2deg;
    //-------------------------------------------------------------------------//
    if (fabs(S_rad) < epsilon){S_rad = 0.0;}
    if (fabs(D_rad) < epsilon){D_rad = 0.0;}
    if (fabs(R_rad) < epsilon){R_rad = 0.0;}

    M1 =-( sin(D_rad)*cos(R_rad)*sin(2.0*S_rad)
         + sin(2.0*D_rad)*sin(R_rad)*pow(sin(S_rad),2) );
    M2 = ( sin(D_rad)*cos(R_rad)*sin(2.0*S_rad)
         - sin(2.0*D_rad)*sin(R_rad)*pow(cos(S_rad),2) );
    M3 = sin(2.0*D_rad)*sin(R_rad);
    M4 = sin(D_rad)*cos(R_rad)*cos(2.0*S_rad)
       + 0.5*sin(2.0*D_rad)*sin(R_rad)*sin(2.0*S_rad);
    M5 =-( cos(D_rad)*cos(R_rad)*cos(S_rad)
         + cos(2.0*D_rad)*sin(R_rad)*sin(S_rad) );
    M6 =-( cos(D_rad)*cos(R_rad)*sin(S_rad)
         - cos(2.0*D_rad)*sin(R_rad)*cos(S_rad));

    moments[0] = M1;
    moments[1] = M2;
    moments[2] = M3;
    moments[3] = M4;
    moments[4] = M5;
    moments[5] = M6;
    return;
}
//============================================================================//
/*!
 * @brief Uniquify euler angle triplet.
 *
 * Puts euler angles into ranges compatible with (dip,strike,-rake) in
 * seismology:
 *
 * @param[in,out] alpha    dip angle to put into range [0, pi/2]
 * @param[in,out] beta     strike angle to put into range [0, 2*pi)
 * @param[in,out] gamma    -rake angle to put into range [-pi, pi)
 *
 * @note If alpha is near to zero, beta is replaced by beta+gamma and gamma is
 *       set to zero, to prevent that additional ambiguity.
 *
 *       If alpha is near to pi/2, beta is put into the range [0,pi).
 *
 */
void cmopad_uniqueEuler(double *alpha, double *beta, double *gamma)
{
    const char *fcnm = "cmopad_uniqueEuler\0";

    *alpha = numpy_mod(*alpha, 2.0*M_PI);

    if (0.5*M_PI < *alpha && *alpha <= M_PI)
    {
        *alpha = M_PI - *alpha;
        *beta  = *beta + M_PI;
        *gamma = 2.0*M_PI - *gamma;
    }
    else if(M_PI < *alpha && *alpha <= 1.5*M_PI)
    {
        *alpha = *alpha - M_PI;
        *gamma = M_PI - *gamma;
    }
    else if(1.5*M_PI < *alpha && *alpha <= 2.0*M_PI)
    {
        *alpha = 2.0*M_PI - *alpha;
        *beta  = *beta + M_PI;
        *gamma = M_PI + *gamma;
    }

    *alpha = numpy_mod(*alpha, 2.0*M_PI);
    *beta  = numpy_mod(*beta,  2.0*M_PI);
    *gamma = numpy_mod(*gamma+M_PI, 2.0*M_PI) - M_PI;

    // If dip is exactly 90 degrees, one is still
    // free to choose between looking at the plane from either side.
    // Choose to look at such that beta is in the range [0,180)

    // This should prevent some problems, when dip is close to 90 degrees:
    if (fabs(*alpha - 0.5*M_PI) < 1e-10){*alpha = 0.5*M_PI;}
    if (fabs(*beta - M_PI) < 1e-10){*beta = M_PI;}
    if (fabs(*beta - 2.*M_PI) < 1e-10){*beta = 0.;}
    if (fabs(*beta) < 1e-10){*beta = 0.;}

    if ((fabs(*alpha - 0.5*M_PI) < 1.e-14) && *beta >= M_PI)
    {
        *gamma =-*gamma;
        *beta  = numpy_mod(*beta-M_PI, 2.0*M_PI);
        *gamma = numpy_mod(*gamma+M_PI, 2.0*M_PI) - M_PI;
        if (!(0. <= *beta && *beta < M_PI))
        {
            printf("%s: Warning on bounds on beta!", fcnm);
        }
        if (!(-M_PI <= *gamma && *gamma < M_PI))
        {
            printf("%s: Warning on bounds on gamma!", fcnm);
        }
    }

    if (*alpha < 1e-7)
    {
        *beta = numpy_mod(*beta + *gamma, 2.0*M_PI);
        *gamma = 0.;
    }
    return;
}
//============================================================================//
/*!
 * @brief Ascending argument sort for three numbers.  For more information:
 *        http://stackoverflow.com/questions/4367745/simpler-way-of-sorting-three-numbers
 *
 * @param[in] x       array to sort into ascending[3]
 *
 * @param[out] iperm  permutation such that x3[iperm[:]] is in ascending
 *                    order
 *
 */
static int __cmopad_argsort3(double x[3], int iperm[3])
{
    const char *fcnm = "__cmopad_argsort3\0";
    int i, temp;
    const int a = 0;
    const int b = 1;
    const int c = 2;
    // Copy
    iperm[a] = a;
    iperm[b] = b;
    iperm[c] = c;
    if (x[iperm[a]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[a];
        iperm[a] = temp;
    }
    if (x[iperm[a]] > x[iperm[b]])
    {
        temp = iperm[b];
        iperm[b] = iperm[a];
        iperm[a] = temp;
    }
    //Now the smallest element is the first one. Just check the 2-nd and 3-rd
    if (x[iperm[b]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[b];
        iperm[b] = temp;
    }
    // Verify
    for (i=1; i<3; i++)
    {
        if (x[iperm[i-1]] > x[iperm[i]])
        {
            printf("%s: Failed to sort numbers in ascending order\n", fcnm);
            return -1;
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Determinant of a 3 x 3 matrix
 *
 * @result determinant of a 3 x 3 matrix
 *
 */
static double __cmopad_determinant3x3(double A[3][3])
{
    double determinant;
    determinant = A[0][0]*( A[1][1]*A[2][2] - A[1][2]*A[2][1])
                - A[0][1]*( A[1][0]*A[2][2] - A[1][2]*A[2][0])
                + A[0][2]*( A[1][0]*A[2][1] - A[1][1]*A[2][0]);
    return determinant;
}
//============================================================================//
/*!
 * @brief Computes the inverse of a 3 x 3 matrix
 *
 * @param[in,out] Amat   on input the 3 x 3 matrix to invert.
 *                       on output the inverted matrix
 *
 * @result 0 indicates success
 *
 */
static int __cmopad_inv3(double Amat[3][3])
{
    const char *fcnm = "__cmopad_inv3\0";
    double a11, a12, a13, a21, a22, a23, a31, a32, a33, detA, detAi;
    // Determinant of a
    detA = __cmopad_determinant3x3(Amat);
    if (detA == 0.0)
    {
        printf("%s: Error matrix is singular!\n", fcnm);
        return -1;
    }
    detAi = 1.0/detA;
    a11 = Amat[0][0]; a12 = Amat[0][1]; a13 = Amat[0][2];
    a21 = Amat[1][0]; a22 = Amat[1][1]; a23 = Amat[1][2];
    a31 = Amat[2][0]; a32 = Amat[2][1]; a33 = Amat[2][2];
    // Row 1
    Amat[0][0] = detAi*(a22*a33 - a32*a23);
    Amat[0][1] = detAi*(a13*a32 - a12*a33);
    Amat[0][2] = detAi*(a12*a23 - a13*a22);
    // Row 2
    Amat[1][0] = detAi*(a23*a31 - a21*a33);
    Amat[1][1] = detAi*(a11*a33 - a13*a31);
    Amat[1][2] = detAi*(a13*a21 - a11*a23);
    // Row 3
    Amat[2][0] = detAi*(a21*a32 - a22*a31);
    Amat[2][1] = detAi*(a12*a31 - a11*a32);
    Amat[2][2] = detAi*(a11*a22 - a12*a21);
    return 0;
}
//============================================================================//
/*!
 * @brief Computes the trace (sum of diagonal) of a 3x3 matrix
 *
 * @param[in] M     matrix of which to compute trace
 *
 * @result trace of matrix
 *
 */
static double __cmopad_trace3(double M[3][3])
{
    double trace;
    trace = M[0][0] + M[1][1] + M[2][2];
    return trace;
}
//===========================================================================//
/*!
 * @brief Calculates all eigenvalues and eigenvectors (job = 0) of a
 *        symmetric 3 x 3 matrix using LAPACK.
 *
 * @param[in,out] a  on input matrix of which to calculate eigenvalues
 *                   and optionally eigenvectors.
 *                   on output if the eigenvectors are desired they are
 *                   stored on a
 *
 * @param[in] job    = 0 calculate eigenvectors, otherwise, just eigenvalues
 *
 * @param[out] eigs  eigenvalues of matrix a [3]
 *
 * @result 0 indicates success
 *
 * @author B. Baker (ISTI)
 * @date April 2016
 */
static int __cmopad_Eigs3x3(double a[3][3], int job, double eigs[3])
{
    const char *fcnm = "__cmopad_Eigs3x3\0";
    double avec[9];
    int n = 3;
    int i, j, indx, ierr;
    int lda = n;
    for (i = 0; i < n; i++)
    {
        for (j = i; j < n; j++)
        {
            // fill lower half
            if (i != j)
            {
                indx = i * lda + j;
                avec[indx] = a[j][i];
            }
            indx = j * lda + i;
            avec[indx] = a[i][j];
        }
    }
    // Calculate eigenvalues/eigenvectors
    if (job != 0)
    {
        ierr = LAPACKE_dsyev(LAPACK_COL_MAJOR, 'N', 'L', n, avec, lda, eigs);
    }
    else
    {
        ierr = LAPACKE_dsyev(LAPACK_COL_MAJOR, 'V', 'L', n, avec, lda, eigs);
    }
    if (ierr != 0)
    {
        printf("%s: Error computing eigenvalues\n", fcnm);
        ierr = 1;
    }
    // Copy eigenvectors back onto A
    if (job == 0 && ierr == 0)
    {
        indx = 0;
        for (j = 0; j < n; j++)
        {
            for (i = 0; i < n; i++)
            {
                indx = j * lda + i;
                a[i][j] = avec[indx];
            }
        }
    }
    return 0;
}
/*!
 * @brief Function for swapping to two numbers a, b
 *
 * @param[in,out] a    on input a = a, on output a = b
 * @param[in,out] b    on input b = b, on output b = a
 *
 */
static void __cmopad_swap8(double *a, double *b)
{
    double temp;
    temp = *a;
    *a = *b;
    *b = temp;
    return;
}

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach


text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main(int argc, char *argv[]) {
    unsigned long seed = 1234;
    int n = 10;
    if (argc > 1)
        n = atoi(argv[1]);
    gsl_rng_env_setup();
    gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937);
    gsl_rng_set(rng, seed);
    for (int i = 0; i < n; i++) {
        double x = gsl_rng_uniform(rng);
        printf("%.12f\n", x);
    }
    gsl_rng_free(rng);
    return EXIT_SUCCESS;
}

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material