text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iniparser.h>
#include "parmt_utils.h"
#include "prepmt/prepmt_dataArchive.h"
#include "prepmt/prepmt_greens.h"
#include "prepmt/prepmt_hudson96.h"
#include "prepmt/prepmt_commands.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "ispl/process.h"
#include "iscl/array/array.h"
#include "iscl/fft/fft.h"
#include "iscl/memory/memory.h"
#include "iscl/os/os.h"
#include "iscl/signal/convolve.h"
#include "iscl/signal/filter.h"

static int getPrimaryArrival(const struct sacHeader_struct hdr,
                             double *time, char phaseName[8]);
static int getPrimaryArrivalNoPolarity(const struct sacHeader_struct hdr,
                                       double *time, char phaseName[8]);
static void shiftGreens(const int npgrns, const int lag,
                        double *__restrict__ G, double *__restrict__ work);

//============================================================================//
/*!
 * @brief Converts the fundamental faults Green's functions to Green's
 *        functions that can be used by parmt.
 *
 * @param[in] nobs      Number of observations.
 * @param[in] obs       Observed waveforms to contextualize.  This is an
 *                      array of length [nobs].
 * @param[in] ndepth    Number of depths.
 * @param[in] ntstar    Number of t*s'.
 * @param[in] ffGrns    Fundamental fault Green's functions for every t* and
 *                      depth in the grid search for each observation.  This
 *                      is an array of dimension [nobs x ndepth x ntstar x 10].
 *
 * @param[in,out] grns  Holds space for Green's functions.
 *                      Contains the Green's functions that can be applied to
 *                      a moment tensor to produce a synthetic for every t*
 *                      and depth in the grid search for each observation.
 *                      This is an array of length [nobs x ndepth x ntstar x 6].
 *
 * @result 0 indicates success.
 *
 * @author Ben Baker, ISTI
 *
 */
int prepmt_greens_ffGreensToGreens(const int nobs,
                                   const struct sacData_struct *obs,
                                   const int ndepth,
                                   const int ntstar,
                                   const struct sacData_struct *ffGrns,
                                   struct sacData_struct *grns)
{
    char knetwk[8], kstnm[8], kcmpnm[8], khole[8], phaseName[8],
         phaseNameGrns[8];
    double az, baz, cmpaz, cmpinc, cmpincSEED, epoch, epochNew,
           evla, evlo, o, pick, pickTime, pickTimeGrns, stel, stla, stlo;
    int indices[6], i, icomp, id, ierr, idx, iobs, it, kndx, l, npts;
    const char *kcmpnms[6] = {"GXX\0", "GYY\0", "GZZ\0",
                              "GXY\0", "GXZ\0", "GYZ\0"};
/*
    const double xmom = 1.0;     // no confusing `relative' magnitudes
    const double xcps = 1.e-20;  // convert dyne-cm mt to output cm
    const double cm2m = 1.e-2;   // cm to meters
    const double dcm2nm = 1.e+7; // magnitudes intended to be specified in
                                 // Dyne-cm but I work in N-m
    // Given a M0 in Newton-meters get a seismogram in meters
    const double xscal = xmom*xcps*cm2m*dcm2nm;
*/
    const int nTimeVars = 11;
    const enum sacHeader_enum pickVars[11]
       = {SAC_FLOAT_A,
          SAC_FLOAT_T0, SAC_FLOAT_T1, SAC_FLOAT_T2, SAC_FLOAT_T3,
          SAC_FLOAT_T4, SAC_FLOAT_T5, SAC_FLOAT_T6, SAC_FLOAT_T7,
          SAC_FLOAT_T8, SAC_FLOAT_T9};
    //memset(grns, 0, sizeof(struct tdSearchGreens_struct));
    for (iobs=0; iobs<nobs; iobs++)
    {
        ierr = 0;
        ierr += sacio_getFloatHeader(SAC_FLOAT_AZ,
                                     obs[iobs].header, &az);
        ierr += sacio_getFloatHeader(SAC_FLOAT_BAZ,
                                     obs[iobs].header, &baz);
        ierr += sacio_getFloatHeader(SAC_FLOAT_CMPINC,
                                     obs[iobs].header, &cmpinc);
        ierr += sacio_getFloatHeader(SAC_FLOAT_CMPAZ,
                                     obs[iobs].header, &cmpaz);
        ierr += sacio_getFloatHeader(SAC_FLOAT_EVLA,
                                     obs[iobs].header, &evla);
        ierr += sacio_getFloatHeader(SAC_FLOAT_EVLO,
                                     obs[iobs].header, &evlo);
        ierr += sacio_getFloatHeader(SAC_FLOAT_STLA,
                                     obs[iobs].header, &stla);
        ierr += sacio_getFloatHeader(SAC_FLOAT_STLO,
                                     obs[iobs].header, &stlo);
        ierr += sacio_getCharacterHeader(SAC_CHAR_KNETWK,
                                         obs[iobs].header, knetwk);
        ierr += sacio_getCharacterHeader(SAC_CHAR_KSTNM,
                                         obs[iobs].header, kstnm);
        ierr += sacio_getCharacterHeader(SAC_CHAR_KCMPNM,
                                         obs[iobs].header, kcmpnm);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error reading header variables\n", __func__);
            break;
        }
        // This one isn't critical but it would be nice to have
        sacio_getFloatHeader(SAC_FLOAT_STEL,
                             obs[iobs].header, &stel);
        // Station location code is not terribly important
        sacio_getCharacterHeader(SAC_CHAR_KHOLE, obs[iobs].header, khole);
        cmpincSEED = cmpinc - 90.0; // SAC to SEED convention
        // Get the primary arrival
        ierr = sacio_getEpochalStartTime(obs[iobs].header, &epoch);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting start time\n", __func__);
            break;
        }
        ierr += getPrimaryArrivalNoPolarity(obs[iobs].header, &pickTime,
                                            phaseName);
        //ierr += getPrimaryArrival(obs[iobs].header, &pickTime, phaseName);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting primary pick\n", __func__);
            break;
        }
        // Need to figure out the component
        ierr = parmt_utils_getComponent(kcmpnm, cmpinc, &icomp);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting component\n", __func__);
            break;
        }
/*
        icomp = 1;
        if (fabs(cmpinc - 0.0) < 1.e-4 || fabs(cmpinc - 180.0) < 1.e-4)
        {
            if (kcmpnm[2] == 'Z' || kcmpnm[2] == 'z' || kcmpnm[2] == '1')
            {
                icomp = 1;
            }
            else
            {
                fprintf(stderr, "%s: cmpinc is 0 but channel=%s isn't vertical",
                        __func__, kcmpnm);
            }
        }
        else if (fabs(cmpinc - 90.0) < 1.e-4)
        {
            icomp = 2;
            if (kcmpnm[2] == 'N' || kcmpnm[2] == 'n' ||
                kcmpnm[2] == '1' || kcmpnm[2] == '2')
            {
                icomp = 2;
            }
            else if (kcmpnm[2] == 'N' || kcmpnm[2] == 'n' ||
                     kcmpnm[2] == '2' || kcmpnm[2] == '3')
            {
                icomp = 3;
            }
            else
            {
                fprintf(stderr, "%s: cmpinc is 90 but channel=%s is weird",
                        __func__, kcmpnm);
                return -1;
            }
        }
*/
        /*
        else if (kcmpnm[2] == 'E' || kcmpnm[2] == 'e' || kcmpnm[2] == '3')
        {
            icomp = 3;
        }
        */
/*
        else
        {
            fprintf(stderr, "%s: Can't classify component %s with cmpinc=%f\n",
                    __func__, kcmpnm, cmpinc);
        }
*/
        // Process all Green's functions in this block
        for (id=0; id<ndepth; id++)
        {
            for (it=0; it<ntstar; it++)
            {
                idx = prepmt_hudson96_observationDepthTstarToIndex(
                             ndepth, ntstar, iobs, id, it);
                kndx = 10*idx;
                sacio_getFloatHeader(SAC_FLOAT_O, ffGrns[kndx].header, &o);
                getPrimaryArrival(ffGrns[kndx].header,
                                  &pickTimeGrns, phaseNameGrns);
                if (strcasecmp(phaseNameGrns, phaseName) != 0)
                {
                    fprintf(stdout, "%s: Phase name mismatch %s %s\n",
                            __func__, phaseName, phaseNameGrns);
                }
                npts = ffGrns[kndx].npts;
                ierr = prepmt_greens_getHudson96GreensFunctionsIndices(
                                 nobs, ntstar, ndepth,
                                 iobs, it, id,
                                 indices);
                /*
                indx = prepmt_greens_getHudson96GreensFunctionIndex(G11_GRNS,
                                                            nobs, ntstar, ndepth,
                                                            iobs, it, id);
                */
                for (i=0; i<6; i++)
                {
                    sacio_copy(ffGrns[kndx], &grns[indices[i]]);
                    if (obs[iobs].pz.lhavePZ)
                    {
                        sacio_copyPolesAndZeros(obs[iobs].pz,
                                                &grns[indices[i]].pz);
                    }
                    sacio_setFloatHeader(SAC_FLOAT_AZ, az,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_BAZ, baz,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_CMPAZ, cmpaz,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_CMPINC, cmpinc,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_EVLA, evla,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_EVLO, evlo,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_STLA, stla,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_STLO, stlo,
                                         &grns[indices[i]].header);
                    sacio_setFloatHeader(SAC_FLOAT_STEL, stel,
                                         &grns[indices[i]].header);
                    sacio_setCharacterHeader(SAC_CHAR_KNETWK, knetwk,
                                             &grns[indices[i]].header);
                    sacio_setCharacterHeader(SAC_CHAR_KSTNM, kstnm,
                                             &grns[indices[i]].header);
                    sacio_setCharacterHeader(SAC_CHAR_KHOLE, khole,
                                             &grns[indices[i]].header);
                    sacio_setCharacterHeader(SAC_CHAR_KCMPNM, kcmpnms[i],
                                             &grns[indices[i]].header);
                    sacio_setCharacterHeader(SAC_CHAR_KEVNM, "SYNTHETIC\0",
                                             &grns[indices[i]].header);
                    // Set the start time by aligning on the arrival
                    epochNew = epoch + (pickTime - o) - pickTimeGrns;
                    sacio_setEpochalStartTime(epochNew,
                                              &grns[indices[i]].header);
                    // Update the pick times
                    for (l=0; l<nTimeVars; l++)
                    {
                        ierr = sacio_getFloatHeader(pickVars[l],
                                                    grns[indices[i]].header,
                                                    &pick);
                        if (ierr == 0)
                        {
                            pick = pick + o;
                            sacio_setFloatHeader(pickVars[l],
                                                 pick,
                                                 &grns[indices[i]].header);
                        }
                    }
                }
                //printf("%d %d\n", kndx, indices[0]);
                //printf("%e\n", array_max64f(npts, ffGrns[kndx+0].data));
                ierr = parmt_utils_ff2mtGreens64f(npts, icomp,
                                                  az, baz,
                                                  cmpaz, cmpincSEED,
                                                  ffGrns[kndx+0].data,
                                                  ffGrns[kndx+1].data,
                                                  ffGrns[kndx+2].data,
                                                  ffGrns[kndx+3].data,
                                                  ffGrns[kndx+4].data,
                                                  ffGrns[kndx+5].data,
                                                  ffGrns[kndx+6].data,
                                                  ffGrns[kndx+7].data,
                                                  ffGrns[kndx+8].data,
                                                  ffGrns[kndx+9].data,
                                                  grns[indices[0]].data,
                                                  grns[indices[1]].data,
                                                  grns[indices[2]].data,
                                                  grns[indices[3]].data,
                                                  grns[indices[4]].data,
                                                  grns[indices[5]].data);
                if (ierr != 0)
                {
                    fprintf(stderr, "%s: Failed to rotate Greens functions\n",
                            __func__);
                }
                // Fix the characteristic magnitude scaling in CPS
                // N.B. this is done early now
                /*
                for (i=0; i<6; i++)
                {
                    cblas_dscal(npts, xscal, grns[indices[i]].data, 1);
                    //printf("%d %d %e\n", kndx, indices[i],
                    //          array_max64f(npts,grns[indices[i]].data));
                }
                */
            }
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Convenience function which returns the index of the Green's
 *        function on the Green's function structure.
 *
 * @param[in] GMT_TERM    Name of the desired Green's function:
 *                        (G11_TERM, G22_TERM, ..., G23_TERM).
 * @param[in] nobs        Number of observations.
 * @param[in] ntstar      Number of t*'s.
 * @param[in] ndepth      Number of depths.
 * @param[in] iobs        Desired observation number (C numbering).
 * @param[in] itstar      Desired t* (C numbering).
 * @param[in] idepth      Desired depth (C numbering).
 *
 * @result Negative indicates failure.  Otherwise, this is the index in
 *         grns.grns corresponding to the desired
 *         (iobs, idepth, itstar, G??_GRNS) coordinate.
 *
 * @author Ben Baker, ISTI
 *
 */
int prepmt_greens_getHudson96GreensFunctionIndex(
    const enum prepmtGreens_enum GMT_TERM,
    const int nobs, const int ntstar, const int ndepth,
    const int iobs, const int itstar, const int idepth)
{
    int igx, indx, ngrns;
    ngrns = 6*nobs*ntstar*ndepth;
    indx =-1;
    igx = (int) GMT_TERM - 1;
    if (igx < 0 || igx > 5)
    {
        fprintf(stderr, "%s: Can't classify Green's functions index\n",
                __func__);
        return indx;
    }
    indx = iobs*(6*ntstar*ndepth)
         + idepth*(6*ntstar)
         + itstar*6
         + igx;
    if (indx < 0 || indx >= ngrns)
    {
        fprintf(stderr, "%s: indx out of bounds - segfault is coming\n",
                __func__);
        return -1;
    }
    return indx;
}
//============================================================================//
/*!
 * @brief Repicks the Green's functions onset time with an STA/LTA picker.
 *
 * @param[in] sta        Short term average window length (seconds).
 * @param[in] lta        Long term average window length (seconds).
 * @param[in] threshPct  Percentage of max STA/LTA after which an arrival
 *                       is declared.
 * @param[in] nobs       Number of observations.
 * @param[in] ntstar     Number of t*'s.
 * @param[in] ndepth     Number of depths.
 * @param[in] iobs       C numbered observation index.
 * @param[in] itstar     C numbered t* index.
 * @param[in] idepth     C numbered depth index.
 *
 * @param[in,out] grns   On input contains the Green's functions.
 *                       On output the first defined arrival time has been
 *                       modified with an STA/LTA picker.  The header variable
 *                       to be modified is most likely SAC_FLOAT_A.
 *                       This is an array of dimension [6].
 *
 * @brief 0 indicates success.
 *
 * @author Ben Baker, ISTI
 *
 */
int prepmt_greens_repickGreensWithSTALTA(
    const double sta, const double lta, const double threshPct,
    struct sacData_struct *grns)
{
    struct stalta_struct stalta;
    enum sacHeader_enum pickHeader;
    double *charFn, *g, *Gxx, *Gyy, *Gzz, *Gxy, *Gxz, *Gyz,
           *gxxPad, *gyyPad, *gzzPad, *gxyPad, *gxzPad, *gyzPad,
           apick, charMax, dt, tpick;
    int ierr, k, npts, nlta, npad, nsta, nwork, prePad;
    const int nTimeVars = 11;
    const enum sacHeader_enum pickVars[11]
       = {SAC_FLOAT_A,
          SAC_FLOAT_T0, SAC_FLOAT_T1, SAC_FLOAT_T2, SAC_FLOAT_T3,
          SAC_FLOAT_T4, SAC_FLOAT_T5, SAC_FLOAT_T6, SAC_FLOAT_T7,
          SAC_FLOAT_T8, SAC_FLOAT_T9};
    // Check STA/LTA
    ierr = 0;
    memset(&stalta, 0, sizeof(struct stalta_struct));
    if (lta < sta || sta < 0.0)
    {
        if (lta < sta){fprintf(stderr, "%s: Error lta < sta\n", __func__);}
        if (sta < 0.0){fprintf(stderr, "%s: Error sta is < 0\n", __func__);}
        return -1;
    }
    ierr = sacio_getIntegerHeader(SAC_INT_NPTS,
                                  grns[0].header, &npts);
    if (ierr != 0 || npts < 1)
    {
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting number of points from header\n",
                    __func__);
        }
        else
        {
            fprintf(stderr, "%s: Error no data points\n", __func__);
        }
        return -1;
    }
    ierr = sacio_getFloatHeader(SAC_FLOAT_DELTA,
                                grns[0].header, &dt);
    if (ierr != 0 || dt <= 0.0)
    {
        if (ierr != 0){fprintf(stderr, "%s: failed to get dt\n", __func__);}
        if (dt <= 0.0)
        {
            fprintf(stderr, "%s: invalid sampling period\n", __func__);
        }
        return -1;
    }
    // Define the windows
    nsta = (int) (sta/dt + 0.5);
    nlta = (int) (lta/dt + 0.5);
    prePad = MAX(64, fft_nextpow2(nlta, &ierr));
    npad = prePad + npts;
    // Set space
    gxxPad = memory_calloc64f(npad);
    gyyPad = memory_calloc64f(npad);
    gzzPad = memory_calloc64f(npad);
    gxyPad = memory_calloc64f(npad);
    gxzPad = memory_calloc64f(npad);
    gyzPad = memory_calloc64f(npad);
    charFn = memory_calloc64f(npad);
    // Reference pointers
    Gxx = grns[0].data;
    Gyy = grns[1].data;
    Gzz = grns[2].data;
    Gxy = grns[3].data;
    Gxz = grns[4].data;
    Gyz = grns[5].data;
    // Pre-pad signals
    array_set64f_work(prePad, Gxx[0], gxxPad);
    array_set64f_work(prePad, Gyy[0], gyyPad);
    array_set64f_work(prePad, Gzz[0], gzzPad);
    array_set64f_work(prePad, Gxy[0], gxyPad);
    array_set64f_work(prePad, Gxz[0], gxzPad);
    array_set64f_work(prePad, Gyz[0], gyzPad);
    // Copy rest of array
    array_copy64f_work(npts, Gxx, &gxxPad[prePad]);
    array_copy64f_work(npts, Gyy, &gyyPad[prePad]);
    array_copy64f_work(npts, Gzz, &gzzPad[prePad]);
    array_copy64f_work(npts, Gxy, &gxyPad[prePad]);
    array_copy64f_work(npts, Gxz, &gxzPad[prePad]);
    array_copy64f_work(npts, Gyz, &gyzPad[prePad]);
    // apply the sta/lta
    for (k=0; k<6; k++)
    {
        g = NULL;
        if (k == 0)
        {
            g = gxxPad;
        }
        else if (k == 1)
        {
            g = gyyPad;
        }
        else if (k == 2)
        {
            g = gzzPad;
        }
        else if (k == 3)
        {
            g = gxyPad;
        }
        else if (k == 4)
        {
            g = gxzPad;
        }
        else if (k == 5)
        {
            g = gyzPad;
        }
        ierr = stalta_setShortAndLongTermAverage(nsta, nlta, &stalta);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error setting STA/LTA\n", __func__);
            break;
        }
        ierr = stalta_setData64f(npad, g, &stalta);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error setting data\n", __func__);
            break;
        }
        ierr = stalta_applySTALTA(&stalta);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error applying STA/LTA\n", __func__);
            break;
        }
        ierr = stalta_getData64f(stalta, npad, &nwork, g);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting result\n", __func__);
            break;
        }
        cblas_daxpy(npad, 1.0, g, 1, charFn, 1);
        stalta_resetInitialConditions(&stalta);
        stalta_resetFinalConditions(&stalta);
        g = NULL;
    }
    // Compute the pick time
    charMax = array_max64f(npts, &charFn[prePad], &ierr);
    tpick =-1.0;
    for (k=prePad; k<npad; k++)
    {
        if (charFn[k] > 0.01*threshPct*charMax)
        {
            tpick = (double) (k - prePad)*dt;
            break;
        }
    }
    if (tpick ==-1.0)
    {
        tpick = (double) (array_argmax64f(npad, charFn, &ierr) - prePad)*dt;
    }
    // Overwrite the pick
    pickHeader = SAC_UNKNOWN_HDRVAR;
    for (k=0; k<nTimeVars; k++)
    {
        if (sacio_getFloatHeader(pickVars[k], grns[k].header, &apick) == 0)
        {
            pickHeader = pickVars[k];
            break;
         }
    }
    if (pickHeader == SAC_UNKNOWN_HDRVAR)
    {
        fprintf(stdout, "%s: Could not locate primary arrival - assume A\n",
                __func__);
        pickHeader = SAC_FLOAT_A;
    }
//printf("%f %f\n", tpick, apick);
    //double apick;
    //sacio_getFloatHeader(SAC_FLOAT_A, grns[0].header, &apick);
    for (k=0; k<6; k++)
    {
        sacio_setFloatHeader(pickHeader, tpick,
                             &grns[k].header);
    }
    // Dereference pointers and free space
    Gxx = NULL;
    Gyy = NULL;
    Gzz = NULL;
    Gxy = NULL;
    Gxz = NULL;
    Gyz = NULL;
    memory_free64f(&gxxPad);
    memory_free64f(&gyyPad);
    memory_free64f(&gzzPad);
    memory_free64f(&gxyPad);
    memory_free64f(&gxzPad);
    memory_free64f(&gyzPad);
    memory_free64f(&charFn);
    stalta_free(&stalta);
    return ierr;
}
//============================================================================//
int prepmt_greens_processHudson96Greens(
    const int nobs, const int ntstar, const int ndepth,
    const struct prepmtCommands_struct cmds,
    struct sacData_struct *grns)
{
    struct serialCommands_struct commands;
    struct parallelCommands_struct parallelCommands;
    double *G, dt, dt0, epoch, epoch0, time;
    int *dataPtr, indices[6],
        i, i0, i1, i2, ierr, iobs, idep, it, kndx, l, npts, npts0, nq,
        nwork, ny, ns, nsuse;
    bool lnewDt, lnewStartTime;
    const int nTimeVars = 12;
    const enum sacHeader_enum timeVars[12]
       = {SAC_FLOAT_A, SAC_FLOAT_O,
          SAC_FLOAT_T0, SAC_FLOAT_T1, SAC_FLOAT_T2, SAC_FLOAT_T3,
          SAC_FLOAT_T4, SAC_FLOAT_T5, SAC_FLOAT_T6, SAC_FLOAT_T7,
          SAC_FLOAT_T8, SAC_FLOAT_T9};
    const int spaceInquiry =-1;
    // Loop on the observations
    ns = ndepth*ntstar*6;
    dataPtr = memory_calloc32i(ns);
    for (iobs=0; iobs<nobs; iobs++)
    {
        // Set the processing structure
        memset(&parallelCommands, 0,
               sizeof(struct parallelCommands_struct));
               memset(&commands, 0, sizeof(struct serialCommands_struct));
        kndx = prepmt_greens_getHudson96GreensFunctionIndex(G11_GRNS,
                                                      nobs, ntstar, ndepth,
                                                      iobs, 0, 0);
        // Parse the commands
        ierr = process_stringsToSerialCommandsOptions(
                                      cmds.cmds[iobs].ncmds,
                                      (const char **) cmds.cmds[iobs].cmds,
                                      &commands);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error setting serial command string\n",
                    __func__);
            goto ERROR;
        }
        // Determine some characteristics of the processing
        sacio_getEpochalStartTime(grns[kndx].header, &epoch0);
        sacio_getFloatHeader(SAC_FLOAT_DELTA, grns[kndx].header, &dt0);
        lnewDt = false;
        lnewStartTime = false;
        epoch = epoch0;
        dt = dt0;
        for (i=0; i<commands.ncmds; i++)
        {
            if (commands.commands[i].type == CUT_COMMAND)
            {
                i0 = commands.commands[i].cut.i0;
                epoch = epoch + (double) i0*dt;
                lnewStartTime = true;
            }
            if (commands.commands[i].type == DOWNSAMPLE_COMMAND)
            {
                nq = commands.commands[i].downsample.nq;
                dt = dt*(double) nq;
                lnewDt = true;
            }
            if (commands.commands[i].type == DECIMATE_COMMAND)
            {
                nq = commands.commands[i].decimate.nqAll;
                dt = dt*(double) nq;
                lnewDt = true;
            }
        }
        // Set the commands on the parallel processing structure
        ierr = process_setCommandOnAllParallelCommands(ns, commands,
                                                       &parallelCommands);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error setting the parallel commands\n",
                    __func__);
            goto ERROR;
        }
        // Get the data
        dataPtr[0] = 0;
        for (idep=0; idep<ndepth; idep++)
        {
            for (it=0; it<ntstar; it++)
            {
                kndx = prepmt_greens_getHudson96GreensFunctionIndex(G11_GRNS,
                                                          nobs, ntstar, ndepth,
                                                          iobs, idep, it);
                ierr = sacio_getIntegerHeader(SAC_INT_NPTS,
                                              grns[kndx].header, &npts);
                if (ierr != 0)
                {
                    fprintf(stderr, "%s: Error getting npts\n", __func__);
                    goto ERROR;
                }
                i1 = idep*6*ntstar + it*6 + 0;
                i2 = i1 + 6;
                for (i=i1; i<i2; i++)
                {
                    dataPtr[i+1] = dataPtr[i] + npts;
                }
            }
        }
        nwork = dataPtr[6*ndepth*ntstar];
        if (nwork < 1)
        {
            fprintf(stderr, "%s: Invalid workspace size: %d\n",
                    __func__, nwork);
            ierr = 1;
            goto ERROR;
        }
        G = memory_calloc64f(nwork);
        for (idep=0; idep<ndepth; idep++)
        {
            for (it=0; it<ntstar; it++)
            {
                ierr = prepmt_greens_getHudson96GreensFunctionsIndices(
                                        nobs, ntstar, ndepth,
                                        iobs, it, idep, indices);
                if (ierr != 0)
                {
                    fprintf(stderr, "%s: Error getting index\n", __func__);
                    goto ERROR;
                }
                for (i=0; i<6; i++)
                {
                    i1 = idep*6*ntstar + it*6 + i;
                    //if (i == 0){printf("fill: %d %d %d %e\n", i1, indices[i], dataPtr[i1], array_max64f(grns[indices[i]].npts, grns[indices[i]].data));}

                    cblas_dcopy(grns[indices[i]].npts,
                                grns[indices[i]].data, 1, &G[dataPtr[i1]], 1);
                }
            }
        }
        // Set the data
        ierr =  process_setParallelCommandsData64f(ns, dataPtr,
                                                   G, &parallelCommands);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error setting data\n", __func__);
            goto ERROR;
        }
        // Apply the commands
        ierr = process_applyParallelCommands(&parallelCommands);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error processing data\n", __func__);
            goto ERROR;
        }
        // Get the data
        ierr = process_getParallelCommandsData64f(parallelCommands,
                                                  spaceInquiry, spaceInquiry,
                                                  &ny, &nsuse,
                                                  dataPtr, G);
        if (ny < nwork)
        {
            memory_free64f(&G);
            G = memory_calloc64f(ny);
        }
        ny = nwork;
        ierr = process_getParallelCommandsData64f(parallelCommands,
                                                  nwork, ns,
                                                  &ny, &nsuse, dataPtr, G);
        //printf("%d\n", nsuse);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting data\n", __func__);
            goto ERROR;
        }
        // Unpack the data
        for (idep=0; idep<ndepth; idep++)
        {
            for (it=0; it<ntstar; it++)
            {
                ierr = prepmt_greens_getHudson96GreensFunctionsIndices(
                                        nobs, ntstar, ndepth,
                                        iobs, it, idep, indices);
                if (ierr != 0)
                {
                    fprintf(stderr, "%s: Error getting index\n", __func__);
                    goto ERROR;
                }
                for (i=0; i<6; i++)
                {
                    i1 = idep*6*ntstar + it*6 + i;
                    sacio_getIntegerHeader(SAC_INT_NPTS,
                                           grns[indices[i]].header, &npts0);
                    npts = dataPtr[i1+1] - dataPtr[i1];
                    // Resize event
                    if (npts != npts0)
                    {
                        sacio_freeData(&grns[indices[i]]);
                        grns[indices[i]].data = sacio_malloc64f(npts);
                        grns[indices[i]].npts = npts;
                        sacio_setIntegerHeader(SAC_INT_NPTS, npts,
                                               &grns[indices[i]].header);
                        ierr = array_copy64f_work(npts,
                                                  &G[dataPtr[i1]],
                                                  grns[indices[i]].data);
                    }
                    else
                    {
                        ierr = array_copy64f_work(npts,
                                                  &G[dataPtr[i1]],
                                                  grns[indices[i]].data);
                        //if (i == 0){printf("%d %d %e\n", iobs, indices[i], array_max64f(grns[indices[i]].npts, grns[indices[i]].data));}
                    }
                }
                // Update the times
                if (lnewStartTime)
                {
                    for (i=0; i<6; i++)
                    {
                        // Update the picks
                        for (l=0; l<nTimeVars; l++)
                        {
                            ierr = sacio_getFloatHeader(timeVars[l],
                                          grns[indices[i]].header, &time);
                            if (ierr == 0)
                            {
                                time = time + epoch0; // Turn to real time
                                time = time - epoch;  // Relative to new time
                                sacio_setFloatHeader(timeVars[l], time,
                                                    &grns[indices[i]].header);
                            }
                        } // Loop on picks
                        sacio_setEpochalStartTime(epoch,
                                                  &grns[indices[i]].header);
                    } // Loop on signals
                }
                // Update the sampling period
                if (lnewDt)
                {
                    for (i=0; i<6; i++)
                    {
                        sacio_setFloatHeader(SAC_FLOAT_DELTA, dt,
                                             &grns[indices[i]].header);
                    }
                }
            } // Loop on t*
        } // Loop on depths
        process_freeSerialCommands(&commands);
        process_freeParallelCommands(&parallelCommands);
        memory_free64f(&G);
    }
ERROR:;
    memory_free32i(&dataPtr);
    return 0;
}
//============================================================================//
/*!
 * @brief Writes a Green's functions archive that is appropriate for parmt.
 *
 * @param[in] archiveName   Name of archive file to write.
 * @param[in] nwaves        Number of waveforms.
 * @param[in] ndepths       Number of depths in grid-search.
 * @param[in] evla          Event latitude (degrees).
 * @param[in] evlo          Event longitude (degrees).
 * @param[in] depths        Depths (km) in grid search.  This is an array of
 *                          dimension [ndepths].
 * @param[in] sac           SAC observations to write.  This is an array of
 *                          dimension [nwaves].
 * @param[in] sacGrns       Corresponding Green's functions to write.  This
 *                          is an array of dimension [6 x ndepths x nwaves]
 *                          where the fastest dimension is 6 and the slowest
 *                          dimension is nwaves.  Each Green's function
 *                          group for the waveform, depth pair is packed:
 *                          \f$ \{ G_{xx}, G_{yy}, G_{zz},
 *                                 G_{xy}, G_{xz}, G_{yz} \} \f$.
 *
 * @result 0 indicates success.
 *
 * @author Ben Baker
 *
 */
int prepmt_greens_writeArchive(const char *archiveName, //const char *archiveDir, const char *projnm,
                               const int nwaves, const int ndepths,
                               const double evla, const double evlo,
                               const double *__restrict__ depths,
                               const struct sacData_struct *sac,
                               const struct sacData_struct *sacGrns)
{
    int id, ierr, indx, k;
    hid_t h5fl;
    // initialize the archive
    ierr = prepmt_dataArchive_createArchive(archiveName, //archiveDir, projnm,
                                            ndepths, evla, evlo, depths);
    if (ierr != 0)
    {
         fprintf(stderr, "%s: Error creating archive\n", __func__);
         return -1;
    }
    // open it for writing
    h5fl = prepmt_dataArchive_openArchive(archiveName, &ierr); //archiveDir, projnm, &ierr);
    // write each observation and greens fns
    for (k=0; k<nwaves; k++)
    {
        ierr = prepmt_dataArchive_addObservation(h5fl, sac[k]);
        for (id=0; id<ndepths; id++)
        {
            indx = k*ndepths*6 + id*6;
            ierr = prepmt_dataArchive_addGreensFunctions(h5fl, sac[k],
                                              sacGrns[indx+0], sacGrns[indx+1],
                                              sacGrns[indx+2], sacGrns[indx+3],
                                              sacGrns[indx+4], sacGrns[indx+5]);
            //printf("%f %e\n", sacGrns[indx+0].header.evdp, sacGrns[indx+0].data[12]);
            if (ierr != 0){return -1;}
        }
    }
    prepmt_dataArchive_closeArchive(h5fl);
    return ierr;
}
//============================================================================//
int prepmt_greens_cutHudson96FromData(const int nobs,
                                      const struct sacData_struct *data,
                                      const int ndepth, const int ntstar,
                                      struct sacData_struct *grns)
{
    char **newCmds;
    struct prepmtModifyCommands_struct options;
    double cut0, dt, epoch, epochData;
    int i, ierr, idep, indx, iobs, it, npts;
    const int ncmds = 2;
    //const char *cmds[2] = {"cut\0", "demean\0"};
    struct prepmtCommands_struct prepMTcmds;
    newCmds = (char **) calloc(2, sizeof(char *));
    for (i=0; i<2; i++)
    {
        newCmds[i] = (char *) calloc(MAX_CMD_LEN, sizeof(char));
    }
    strcpy(newCmds[1], "demean\0");

    prepMTcmds.cmds = (struct prepmtCommandsChars_struct *)
                      calloc(1, sizeof(struct prepmtCommandsChars_struct));
    prepMTcmds.cmds[0].ncmds = ncmds;
    memset(&options, 0, sizeof(struct prepmtModifyCommands_struct));
    for (iobs=0; iobs<nobs; iobs++)
    {
        // Get the primary pick
        ierr = sacio_getEpochalStartTime(data[iobs].header, &epoch);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error getting start time\n", __func__);
            break;
        }
        epochData = epoch;
        sacio_getFloatHeader(SAC_FLOAT_DELTA, data[iobs].header, &dt);
        sacio_getIntegerHeader(SAC_INT_NPTS, data[iobs].header, &npts);
        cut0 = epoch;
        //cut1 = cut0 + dt*(double) (npts - 1);
        for (idep=0; idep<ndepth; idep++)
        {
            for (it=0; it<ntstar; it++)
            {
                indx = prepmt_greens_getHudson96GreensFunctionIndex(G11_GRNS,
                                                           nobs, ntstar, ndepth,
                                                           iobs, it, idep);
//printf("%d %d %d %d\n", iobs, it, idep, indx);
                sacio_getEpochalStartTime(grns[indx].header, &epoch);
                if (cut0 - epoch < 0.0)
                {
                    fprintf(stdout, "%s: Cut may be funky\n", __func__);
                }
                options.cut0 = fmax(0.0, cut0 - epoch);
                options.cut1 = options.cut0 + dt*(double) (npts - 1);
/*
                newCmds = prepmt_commands_modifyCommands(ncmds, cmds,
                                                         options, grns[indx],
                                                         &ierr);
*/
                ierr = cut_cutEpochalTimesToString(grns[indx].header.delta,
                                                   0.0,
                                                   options.cut0,
                                                   options.cut1,
                                                   newCmds[0]);
                //strcpy(newCmds[1], "demean");
                //printf("%s %f\n", newCmds[0], grns[indx].header.delta);
                prepMTcmds.cmds[0].cmds = newCmds;
                ierr = prepmt_greens_processHudson96Greens(1, 1, 1,
                                                           prepMTcmds,
                                                           &grns[indx]);
                // Enforce the epochal time
                for (i=0; i<6; i++)
                {
                    sacio_setEpochalStartTime(epochData, &grns[indx+i].header);
                }
//printf("%d %e\n", indx, array_max64f(grns[indx].npts, grns[indx].data));
                for (i=0; i<ncmds; i++)
                {
               //     free(newCmds[i]);
                }
                //free(newCmds);
            }
        }
    }
    for (i=0; i<2; i++)
    {
        free(newCmds[i]);
    }
    free(newCmds[i]);
    free(prepMTcmds.cmds);
    return 0;
}
//============================================================================//
/*!
 * @brief Convenience function for extracting the:
 *        \$ \{ G_{xx}, G_{yy}, G_{zz}, G_{xy}, G_{xz}, G_{yz} \} \$
 *        Green's functions indices for the observation, t*, and depth.
 *
 * @param[in] nobs      Total number of observations.
 * @param[in] ntstar    Total number of t*.
 * @param[in] ndepth    Total number of depths.
 * @param[in] iobs      Observation number.  This is C indexed.
 * @param[in] itstar    t* index.  This is C indexed.
 * @param[in] idepth    Depth index.  This is C indexed.
 * @param[in] grns      Contains the Green's functions.
 * @param[out] indices  Contains the Green's functions indices defining
 *                      the indices that return the:
 *                       \$ \{ G_{xx}, G_{yy}, G_{zz},
 *                             G_{xy}, G_{xz}, G_{yz} \} \$
 *                      for this observation, t*, and depth.
 *
 * @result 0 indicates success.
 *
 * @author Ben Baker, ISTI
 *
 */
int prepmt_greens_getHudson96GreensFunctionsIndices(
    const int nobs, const int ntstar, const int ndepth,
    const int iobs, const int itstar, const int idepth,
    int indices[6])
{
    int i, ierr;
    const enum prepmtGreens_enum mtTerm[6] =
       {G11_GRNS, G22_GRNS, G33_GRNS, G12_GRNS, G13_GRNS, G23_GRNS};
    ierr = 0;
    for (i=0; i<6; i++)
    {
        indices[i] = prepmt_greens_getHudson96GreensFunctionIndex(mtTerm[i],
                                                          nobs, ntstar, ndepth,
                                                          iobs, itstar, idepth);
        if (indices[i] < 0){ierr = ierr + 1;}
    }
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error getting indices (obs,t*,depth)=(%d,%d,%d)\n",
                __func__, iobs, itstar, idepth);
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Aligns the Green's functions to the data via cross-correlation.
 *
 * @param[in] data          Holds the observed data.
 * @param[in] luseEnvelope  If true then cross-correlate the envelopes of
 *                          the waveforms. \n
 *                          Otherwise, cross-correlate the waveforms
 *                          themselves.  In this case the absolute values
 *                          of the cross-correlation are used.
 * @param[in] lnorm         If true then use a normalized cross-correlation
 *                          so that all waveforms can contribute equally. \n
 *                          Otherwise, use the non-normalized cross-correlation.
 * @param[in] maxTimLag     Max time lag (seconds) allowed in cross-correlation.
 *
 * @param[in,out] grns      On input contains the Green's functions to
 *                          align to the data.
 *                          On output contains the Green's functions which
 *                          have been aligned to the data via cross-correlation.
 *
 * @result 0 indicates success.
 *
 * @author Ben Baker, ISTI
 *
 */
int prepmt_greens_xcAlignGreensToData(const struct sacData_struct data,
                                      const bool luseEnvelope, const bool lnorm,
                                      const double maxTimeLag,
                                      struct sacData_struct *grns)
{
    double *dataPad, *Gxx, *Gyy, *Gzz, *Gxy, *Gxz, *Gyz, *xc,
           dt, dtGrns, epochData, epochGrns;
    int ierr, insertData, lxc, maxShift, npadData, npts, npgrns;
    ierr = 0;
    sacio_getIntegerHeader(SAC_INT_NPTS, data.header, &npts);
    sacio_getIntegerHeader(SAC_INT_NPTS, grns[0].header, &npgrns);
    if (npts < 1 || npgrns < 1)
    {
        fprintf(stderr, "%s: Error no data\n", __func__);
        fprintf(stderr, "%s: Error no grns functions\n", __func__);
        return -1;
    }
    sacio_getFloatHeader(SAC_FLOAT_DELTA, data.header, &dt);
    sacio_getFloatHeader(SAC_FLOAT_DELTA, data.header, &dtGrns);
    if (fabs(dt - dtGrns) > 1.e-6 || dt <= 0.0)
    {
        if (dt <= 0.0){fprintf(stderr, "%s: dt is invalid\n", __func__);}
        if (fabs(dt - dtGrns) > 1.e-6)
        {
            fprintf(stderr, "%s: dt is inconsistent\n", __func__);
        }
        return -1;
    }
    // Need a common reference - choose the start time
    ierr = sacio_getEpochalStartTime(data.header,    &epochData);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error getting data start time\n", __func__);
        return -1;
    }
    ierr = sacio_getEpochalStartTime(grns[0].header, &epochGrns);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error getting grns start time\n", __func__);
        return -1;
    }
    if (epochData < epochGrns)
    {
        fprintf(stderr, "%s: epochData < epochGrns not programmed\n", __func__);
        return -1;
    }
    // Insert the data in a epochal time aligned array
    insertData = (int) ((epochData - epochGrns)/dt + 0.5);
    npadData = insertData + npts;
    //printf("%d %d\n", insertData, npts);
    dataPad = memory_calloc64f(npadData);
    array_copy64f_work(npts, data.data, &dataPad[insertData]);
    // Create pointers to Green's functions
    Gxx = grns[0].data;
    Gyy = grns[1].data;
    Gzz = grns[2].data;
    Gxy = grns[3].data;
    Gxz = grns[4].data;
    Gyz = grns[5].data;
    // Compute hte max lag
    maxShift =-1;
    if (maxTimeLag >= 0.0){maxShift = (int) (maxTimeLag/dt + 0.5);}
//printf("%d %f %f\n", maxShift, maxTimeLag/dt, maxTimeLag);
    // Align
    lxc = npadData + npgrns - 1; //npts + npgrns - 1;
    xc = memory_calloc64f(lxc);
    ierr = prepmt_greens_xcAlignGreensToData_work(npadData, dataPad, //data.data,
                                                  npgrns, luseEnvelope, lnorm,
                                                  maxShift,
                                                  Gxx, Gyy, Gzz,
                                                  Gxy, Gxz, Gyz,
                                                  lxc, xc);
    memory_free64f(&xc);
    memory_free64f(&dataPad);
    // Dereference pointers
    Gxx = NULL;
    Gyy = NULL;
    Gzz = NULL;
    Gxy = NULL;
    Gxz = NULL;
    Gyz = NULL;
    return ierr;
}
//============================================================================//
int prepmt_greens_xcAlignGreensToData_work(
    const int npts, double *__restrict__ data,
    const int npgrns,
    const bool luseEnvelope, const bool lnorm,
    const int maxShift,
    double *__restrict__ Gxx,
    double *__restrict__ Gyy,
    double *__restrict__ Gzz,
    double *__restrict__ Gxy,
    double *__restrict__ Gxz,
    double *__restrict__ Gyz,
    const int lxc, double *__restrict__ xc)
{
    double *dwork, *Gwork, *xcorrWork, egrns, esig, xdiv;
    int i, ierr, ierr1, kmt, l1, l2, lag, lcref, nlag;
    //------------------------------------------------------------------------//
    ierr = 0;
    lcref = npgrns + npts - 1;
    if (lcref > lxc){return -1;}
    xcorrWork = memory_calloc64f(lcref);
    array_zeros64f_work(lxc, xc);
    if (luseEnvelope)
    {
        Gwork = memory_calloc64f(npgrns);
        dwork = memory_calloc64f(npts);
        ierr = signal_filter_envelope64f_work(npts, data, dwork);
        esig = cblas_dnrm2(npts, dwork, 1);
    }
    else
    {
        Gwork = NULL;
        dwork = data;
        esig = cblas_dnrm2(npts, data, 1);
    }
    // Loop on the Green's functions
    for (kmt=0; kmt<6; kmt++)
    {
        if (kmt == 0)
        {
            if (!luseEnvelope)
            {
                Gwork = Gxx;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gxx, Gwork);
            }
        }
        else if (kmt == 1)
        {
            if (!luseEnvelope)
            {
                Gwork = Gyy;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gyy, Gwork);
            }
        }
        else if (kmt == 2)
        {
            if (!luseEnvelope)
            {
                Gwork = Gzz;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gzz, Gwork);
            }
        }
        else if (kmt == 3)
        {
            if (!luseEnvelope)
            {
                Gwork = Gxy;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gxy, Gwork);
            }
        }
        else if (kmt == 4)
        {
            if (!luseEnvelope)
            {
                Gwork = Gxz;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gxz, Gwork);
            }
        }
        else
        {
            if (!luseEnvelope)
            {
                Gwork = Gyz;
            }
            else
            {
                signal_filter_envelope64f_work (npgrns, Gyz, Gwork);
            }
        }
        egrns = cblas_dnrm2(npgrns, Gwork, 1);
        ierr1 = signal_convolve_correlate64f_work(npts, dwork,
                                                  npgrns,
                                                  Gwork,
                                                  CONVCOR_FULL,
                                                  lcref, xcorrWork);
        if (ierr1 != 0)
        {
            if (ierr1 != 0)
            {
                fprintf(stderr, "%s: Error correlating kmt %d\n",
                        __func__, kmt);
            }
            ierr = ierr + 1;
        }
        // sum in the result - absolute value handles polarity
        else
        {
            xdiv = 1.0;
            if (lnorm){xdiv = 1.0/(esig*egrns);}
            if (luseEnvelope)
            {
                #pragma omp simd
                for (i=0; i<lcref; i++)
                {
                    xc[i] = xc[i] + xcorrWork[i]*xdiv;
                }
            }
            else
            {
                #pragma omp simd
                for (i=0; i<lcref; i++)
                {
                    xc[i] = xc[i] + fabs(xcorrWork[i]*xdiv);
//if (kmt == 5){printf("%e\n", xc[i]);}
                }
            }
        }
    } // Loop on green's functions
    // Compute the lag where the times range from [-npgrns:npts].  Hence,
    // lag = npgrns is the unlagged cross-correlation.
    if (maxShift < 0)
    {
        lag = array_argmax64f(lcref, xc, &ierr);
    }
    else
    {
        l1 = MAX(0, npgrns - maxShift);
        l2 = MIN(lcref - 1, npgrns + maxShift);
        nlag = l2 - l1 + 1;
        lag = l1 + array_argmax64f(nlag, &xc[l1], &ierr);
    }
    lag =-npgrns + lag;
    //printf("%d %d %d\n", lag, npgrns, npts);
/*
for (i=0; i<npts; i++)
{
printf("%e\n", data[i]);
}
for (i=0; i<npgrns; i++)
{
printf("%e %e %e %e %e %e\n", Gxx[i], Gyy[i], Gzz[i], Gxy[i], Gxz[i], Gyz[i]);
}
for (i=0; i<lcref; i++)
{
printf("%d %e\n",-npgrns+i, xc[i]);
}
printf("%d\n", lag);
*/
//getchar();
    // Shift the greens' functions
    array_zeros64f_work(lcref, xcorrWork);
    shiftGreens(npgrns, lag, Gxx, xcorrWork);
    shiftGreens(npgrns, lag, Gyy, xcorrWork);
    shiftGreens(npgrns, lag, Gzz, xcorrWork);
    shiftGreens(npgrns, lag, Gxy, xcorrWork);
    shiftGreens(npgrns, lag, Gxz, xcorrWork);
    shiftGreens(npgrns, lag, Gyz, xcorrWork);
    if (luseEnvelope)
    {
        memory_free64f(&Gwork);
        memory_free64f(&dwork);
    }
    else
    {
        dwork = NULL;
        Gwork = NULL;
    }
    memory_free64f(&xcorrWork);
    return 0;
}
//============================================================================//
/*!
 * @brief Utility function for use after the Green's functions have been
 *        manually aligned to issue a measure of the alignment quality.
 *
 * @param[in] npts          Number of points in data and Green's functions.
 * @param[in] luseEnvelope  If true then cross-correlate the envelopes of
 *                          the waveforms. \n
 *                          Otherwise, cross-correlate the waveforms
 *                          themselves.  In this case the absolute values
 *                          of the cross-correlation are used.
 * @param[in] lnorm         If true then use a normalized cross-correlation
 *                          so that all waveforms can contribute equally. \n
 *                          Otherwise, use the non-normalized cross-correlation.
 * @param[in] data          Observed waveforms.  This is an array of dimension
 *                          [npts].
 * @param[in] Gxx           Aligned Gxx Green's function.  This is an array of
 *                          dimension [npts].
 * @param[in] Gyy           Aligned Gyy Green's function.  This is an array of
 *                          dimension [npts].
 * @param[in] Gzz           Aligned Gzz Green's function.  This is an array of
 *                          dimension [npts].
 * @param[in] Gxy           Aligned Gxy Green's function.  This is an array of
 *                          dimension [npts].
 * @param[in] Gxz           Aligned Gxz Green's function.  This is an array of
 *                          dimension [npts].
 * @param[in] Gyz           Aligned Gyz Green's function.  This is an array of
 *                          dimension [npts].
 *
 * @param[out] ierr         0 indicates success.
 *
 * @result The zero-lag cross-correlation score of the Green's function to
 *         data alignment.
 *
 * @author Ben Baker, ISTI
 *
 */
double prepmt_greens_scoreXCAlignment(const int npts,
                                      const bool luseEnvelope,
                                      const bool lnorm,
                                      const double *__restrict__ data,
                                      const double *__restrict__ Gxx,
                                      const double *__restrict__ Gyy,
                                      const double *__restrict__ Gzz,
                                      const double *__restrict__ Gxy,
                                      const double *__restrict__ Gxz,
                                      const double *__restrict__ Gyz,
                                      int *ierr)
{
    double *dwork, *Gwork, egrns, esig, xcStack, xdiv;
    const double *G;
    int i;
    enum isclError_enum isclError;
    *ierr = 0;
    xcStack = 0.0;
    dwork = NULL;
    G = NULL;
    Gwork = memory_calloc64f(npts);
    if (luseEnvelope)
    {
        dwork = signal_filter_envelope64f(npts, data, &isclError);
        if (isclError != ISCL_SUCCESS)
        {
            fprintf(stderr, "%s: Error computing data envelope\n", __func__);
            *ierr = 1;
            return xcStack;
        }
        esig = cblas_dnrm2(npts, dwork, 1);
    }
    else
    {
        dwork = array_copy64f(npts, data, &isclError);
        if (isclError != ISCL_SUCCESS)
        {
            fprintf(stderr, "%s: Error copying data\n", __func__);
            *ierr = 1;
            return xcStack;
        }
        esig = cblas_dnrm2(npts, data, 1);
    }
    for (i=0; i<6; i++)
    {
        if (i == 0)
        {
            G = Gxx;
        }
        else if (i == 1)
        {
            G = Gyy;
        }
        else if (i == 2)
        {
            G = Gzz;
        }
        else if (i == 3)
        {
            G = Gxy;
        }
        else if (i == 4)
        {
            G = Gxz;
        }
        else if (i == 5)
        {
            G = Gyz;
        }
        if (luseEnvelope)
        {
            isclError = signal_filter_envelope64f_work(npts, G, Gwork);
            if (isclError != ISCL_SUCCESS)
            {
                fprintf(stderr, "%s: Error computing Grns envelope\n",
                        __func__);
                *ierr = 1;
                return xcStack;
            }
        }
        else
        {
            isclError = array_copy64f_work(npts, G, Gwork);
            if (isclError != ISCL_SUCCESS)
            {
                fprintf(stderr, "%s: Error copying G to Gwork\n", __func__);
                *ierr = 1;
                return xcStack;
            }
        }
        egrns = cblas_dnrm2(npts, Gwork, 1);
        xdiv = 1.0;
        if (lnorm){xdiv = 1.0/(esig*egrns);}
        xdiv = xdiv*1.0/6.0; // Normalization term for sum
        if (luseEnvelope)
        {
            xcStack = xcStack + xdiv*cblas_ddot(npts, dwork, 1, Gwork, 1);
        }
        else
        {
            xcStack = xcStack + xdiv*fabs(cblas_ddot(npts, dwork, 1, Gwork, 1));
        }
        G = NULL;
    }
    memory_free64f(&dwork);
    memory_free64f(&Gwork);
    return xcStack;
}
//============================================================================//
static int getPrimaryArrival(const struct sacHeader_struct hdr,
                             double *time, char phaseName[8])
{
    const enum sacHeader_enum timeVars[11]
       = {SAC_FLOAT_A,
          SAC_FLOAT_T0, SAC_FLOAT_T1, SAC_FLOAT_T2, SAC_FLOAT_T3,
          SAC_FLOAT_T4, SAC_FLOAT_T5, SAC_FLOAT_T6, SAC_FLOAT_T7,
          SAC_FLOAT_T8, SAC_FLOAT_T9};
    const enum sacHeader_enum timeVarNames[11]
       = {SAC_CHAR_KA,
          SAC_CHAR_KT0, SAC_CHAR_KT1, SAC_CHAR_KT2, SAC_CHAR_KT3,
          SAC_CHAR_KT4, SAC_CHAR_KT5, SAC_CHAR_KT6, SAC_CHAR_KT7,
          SAC_CHAR_KT8, SAC_CHAR_KT9};
    int i, ifound1, ifound2;
    memset(phaseName, 0, 8*sizeof(char));
    for (i=0; i<11; i++)
    {
        ifound1 = sacio_getFloatHeader(timeVars[i], hdr, time);
        ifound2 = sacio_getCharacterHeader(timeVarNames[i], hdr, phaseName);
        if (ifound1 == 0 && ifound2 == 0){return 0;}
    }
    fprintf(stderr, "%s: Failed to get primary pick\n", __func__);
    *time =-12345.0;
    memset(phaseName, 0, 8*sizeof(char));
    strcpy(phaseName, "-12345");
    return -1;
}
//============================================================================//
static int getPrimaryArrivalNoPolarity(const struct sacHeader_struct hdr,
                                       double *time, char phaseName[8])
{
    int ierr;
    size_t lenos;
    ierr = getPrimaryArrival(hdr, time, phaseName);
    if (ierr == 0)
    {
        lenos = strlen(phaseName);
        if (lenos > 0)
        {
            if (phaseName[lenos-1] == '+' || phaseName[lenos-1] == '-')
            {
                phaseName[lenos-1] = '\0';
            }
        }
    }
    return ierr;
}
//============================================================================//
static void shiftGreens(const int npgrns, const int lag,
                        double *__restrict__ G, double *__restrict__ work)
{
    int ncopy;
    if (lag == 0){return;}
    array_copy64f_work(npgrns, G, work);
    array_zeros64f_work(npgrns, G);
    if (lag > 0)
    {
        ncopy = npgrns - lag;
        array_copy64f_work(ncopy, work, &G[lag]);
    }
    else
    {
        ncopy = npgrns + lag;
        array_copy64f_work(ncopy, &work[-lag], G);
    }
    return;
}


text:
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "cconfigspace_internal.h"
#include "distribution_internal.h"

struct _ccs_distribution_normal_data_s {
	_ccs_distribution_common_data_t common_data;
	ccs_float_t                     mu;
	ccs_float_t                     sigma;
	ccs_scale_type_t                scale_type;
	ccs_numeric_t                   quantization;
	int                             quantize;
};
typedef struct _ccs_distribution_normal_data_s _ccs_distribution_normal_data_t;

static ccs_result_t
_ccs_distribution_del(ccs_object_t o) {
	(void)o;
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_normal_get_bounds(_ccs_distribution_data_t *data,
                                    ccs_interval_t           *interval_ret);

static ccs_result_t
_ccs_distribution_normal_samples(_ccs_distribution_data_t *data,
                                 ccs_rng_t                 rng,
                                 size_t                    num_values,
                                 ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_normal_strided_samples(_ccs_distribution_data_t *data,
                                         ccs_rng_t                 rng,
                                         size_t                    num_values,
                                         size_t                    stride,
                                         ccs_numeric_t            *values);

static ccs_result_t
_ccs_distribution_normal_soa_samples(_ccs_distribution_data_t  *data,
                                     ccs_rng_t                  rng,
                                     size_t                     num_values,
                                     ccs_numeric_t            **values);

static _ccs_distribution_ops_t _ccs_distribution_normal_ops = {
	{ &_ccs_distribution_del },
	&_ccs_distribution_normal_samples,
	&_ccs_distribution_normal_get_bounds,
	&_ccs_distribution_normal_strided_samples,
	&_ccs_distribution_normal_soa_samples
};

static ccs_result_t
_ccs_distribution_normal_get_bounds(_ccs_distribution_data_t *data,
                                    ccs_interval_t           *interval_ret) {
	_ccs_distribution_normal_data_t *d = (_ccs_distribution_normal_data_t *)data;
	const ccs_numeric_type_t  data_type   = d->common_data.data_types[0];
	const ccs_scale_type_t scale_type     = d->scale_type;
	const ccs_numeric_t    quantization   = d->quantization;
	const int              quantize       = d->quantize;
	ccs_numeric_t          l;
	ccs_bool_t             li;
	ccs_numeric_t          u;
	ccs_bool_t             ui;

	if (scale_type == CCS_LOGARITHMIC) {
		if (data_type == CCS_NUM_FLOAT) {
			if (quantize) {
				l.f = quantization.f;
				li = CCS_TRUE;
			} else {
				l.f = 0.0;
				li = CCS_FALSE;
			}
			u.f = CCS_INFINITY;
			ui = CCS_FALSE;
		} else {
			if (quantize) {
				l.i = quantization.i;
				u.i = (CCS_INT_MAX/quantization.i)*quantization.i;
			} else {
				l.i = 1;
				u.i = CCS_INT_MAX;
			}
			li = CCS_TRUE;
			ui = CCS_TRUE;
		}
	} else {
		if (data_type == CCS_NUM_FLOAT) {
			l.f = -CCS_INFINITY;
			li = CCS_FALSE;
			u.f = CCS_INFINITY;
			ui = CCS_FALSE;
		} else {
			if (quantize) {
				l.i = (CCS_INT_MIN/quantization.i)*quantization.i;
				u.i = (CCS_INT_MAX/quantization.i)*quantization.i;
			} else {
				l.i = CCS_INT_MIN;
				u.i = CCS_INT_MAX;
			}
			li = CCS_TRUE;
			ui = CCS_TRUE;
		}
	}
	interval_ret->type = data_type;
	interval_ret->lower = l;
	interval_ret->upper = u;
	interval_ret->lower_included = li;
	interval_ret->upper_included = ui;
	return CCS_SUCCESS;
}

static inline ccs_result_t
_ccs_distribution_normal_samples_float(gsl_rng                *grng,
                                       const ccs_scale_type_t  scale_type,
                                       const ccs_float_t       quantization,
                                       const ccs_float_t       mu,
                                       const ccs_float_t       sigma,
                                       const int               quantize,
                                       size_t                  num_values,
                                       ccs_float_t            *values) {
	size_t i;
	if (scale_type == CCS_LOGARITHMIC && quantize) {
		ccs_float_t lq = log(quantization*0.5);
		if (mu - lq >= 0.0)
			//at least 50% chance to get a valid value
			for (i = 0; i < num_values; i++)
				do {
					values[i] = gsl_ran_gaussian(grng, sigma) + mu;
				} while (values[i] < lq);
		else
			//use tail distribution
			for (i = 0; i < num_values; i++)
				values[i] = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
	} else
		for (i = 0; i < num_values; i++)
			values[i] = gsl_ran_gaussian(grng, sigma) + mu;
	if (scale_type == CCS_LOGARITHMIC)
		for (i = 0; i < num_values; i++)
			values[i] = exp(values[i]);
	if (quantize) {
			ccs_float_t rquantization = 1.0 / quantization;
			for (i = 0; i < num_values; i++)
				values[i] = round(values[i] * rquantization) * quantization;
	}
	return CCS_SUCCESS;
}

static inline ccs_result_t
_ccs_distribution_normal_samples_int(gsl_rng                *grng,
                                     const ccs_scale_type_t  scale_type,
                                     const ccs_int_t         quantization,
                                     const ccs_float_t       mu,
                                     const ccs_float_t       sigma,
                                     const int               quantize,
                                     size_t                  num_values,
                                     ccs_numeric_t          *values) {
	size_t i;
	ccs_float_t q;
	if (quantize)
		q = quantization*0.5;
	else
		q = 0.5;
	if (scale_type == CCS_LOGARITHMIC) {
		ccs_float_t lq = log(q);
		if (mu - lq >= 0.0)
			for (i = 0; i < num_values; i++)
				do {
					do {
						values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
					} while (values[i].f < lq);
					values[i].f = exp(values[i].f);
				} while (CCS_UNLIKELY(values[i].f - q > (ccs_float_t)CCS_INT_MAX));
		else
			for (i = 0; i < num_values; i++)
				do {
					values[i].f = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
					values[i].f = exp(values[i].f);
				} while (CCS_UNLIKELY(values[i].f - q > (ccs_float_t)CCS_INT_MAX));
	}
	else
		for (i = 0; i < num_values; i++)
			do {
				values[i].f = gsl_ran_gaussian(grng, sigma) + mu;
			} while (CCS_UNLIKELY(values[i].f - q > (ccs_float_t)CCS_INT_MAX || values[i].f + q < (ccs_float_t)CCS_INT_MIN));
	if (quantize) {
		ccs_float_t rquantization = 1.0 / quantization;
		for (i = 0; i < num_values; i++)
			values[i].i = (ccs_int_t)round(values[i].f * rquantization) * quantization;
	} else
		for (i = 0; i < num_values; i++)
			values[i].i = round(values[i].f);
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_normal_samples(_ccs_distribution_data_t *data,
                                  ccs_rng_t                 rng,
                                  size_t                    num_values,
                                  ccs_numeric_t            *values) {
	_ccs_distribution_normal_data_t *d = (_ccs_distribution_normal_data_t *)data;
	const ccs_numeric_type_t  data_type   = d->common_data.data_types[0];
	const ccs_scale_type_t scale_type     = d->scale_type;
	const ccs_numeric_t    quantization   = d->quantization;
	const ccs_float_t      mu             = d->mu;
	const ccs_float_t      sigma          = d->sigma;
	const int              quantize       = d->quantize;
	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));
	if (data_type == CCS_NUM_FLOAT)
		return _ccs_distribution_normal_samples_float(grng, scale_type,
                                                              quantization.f, mu,
                                                              sigma, quantize,
                                                              num_values,
                                                              (ccs_float_t*) values);
	else
		return _ccs_distribution_normal_samples_int(grng, scale_type,
                                                            quantization.i, mu,
                                                            sigma, quantize,
                                                            num_values, values);
}

static inline ccs_result_t
_ccs_distribution_normal_strided_samples_float(gsl_rng                *grng,
                                               const ccs_scale_type_t  scale_type,
                                               const ccs_float_t       quantization,
                                               const ccs_float_t       mu,
                                               const ccs_float_t       sigma,
                                               const int               quantize,
                                               size_t                  num_values,
                                               size_t                  stride,
                                               ccs_float_t            *values) {
	size_t i;
	if (scale_type == CCS_LOGARITHMIC && quantize) {
		ccs_float_t lq = log(quantization*0.5);
		if (mu - lq >= 0.0)
			//at least 50% chance to get a valid value
			for (i = 0; i < num_values; i++)
				do {
					values[i*stride] = gsl_ran_gaussian(grng, sigma) + mu;
				} while (values[i*stride] < lq);
		else
			//use tail distribution
			for (i = 0; i < num_values; i++)
				values[i*stride] = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
	} else
		for (i = 0; i < num_values; i++)
			values[i*stride] = gsl_ran_gaussian(grng, sigma) + mu;
	if (scale_type == CCS_LOGARITHMIC)
		for (i = 0; i < num_values; i++)
			values[i*stride] = exp(values[i*stride]);
	if (quantize) {
			ccs_float_t rquantization = 1.0 / quantization;
			for (i = 0; i < num_values; i++)
				values[i*stride] = round(values[i*stride] * rquantization) * quantization;
	}
	return CCS_SUCCESS;
}

static inline ccs_result_t
_ccs_distribution_normal_strided_samples_int(gsl_rng                *grng,
                                             const ccs_scale_type_t  scale_type,
                                             const ccs_int_t         quantization,
                                             const ccs_float_t       mu,
                                             const ccs_float_t       sigma,
                                             const int               quantize,
                                             size_t                  num_values,
                                             size_t                  stride,
                                             ccs_numeric_t          *values) {
	size_t i;
	ccs_float_t q;
	if (quantize)
		q = quantization*0.5;
	else
		q = 0.5;
	if (scale_type == CCS_LOGARITHMIC) {
		ccs_float_t lq = log(q);
		if (mu - lq >= 0.0)
			for (i = 0; i < num_values; i++)
				do {
					do {
						values[i*stride].f = gsl_ran_gaussian(grng, sigma) + mu;
					} while (values[i*stride].f < lq);
					values[i*stride].f = exp(values[i*stride].f);
				} while (CCS_UNLIKELY(values[i*stride].f - q > (ccs_float_t)CCS_INT_MAX));
		else
			for (i = 0; i < num_values; i++)
				do {
					values[i*stride].f = gsl_ran_gaussian_tail(grng, lq - mu, sigma) + mu;
					values[i*stride].f = exp(values[i*stride].f);
				} while (CCS_UNLIKELY(values[i*stride].f - q > (ccs_float_t)CCS_INT_MAX));
	}
	else
		for (i = 0; i < num_values; i++)
			do {
				values[i*stride].f = gsl_ran_gaussian(grng, sigma) + mu;
			} while (CCS_UNLIKELY(values[i*stride].f - q > (ccs_float_t)CCS_INT_MAX || values[i*stride].f + q < (ccs_float_t)CCS_INT_MIN));
	if (quantize) {
		ccs_float_t rquantization = 1.0 / quantization;
		for (i = 0; i < num_values; i++)
			values[i*stride].i = (ccs_int_t)round(values[i*stride].f * rquantization) * quantization;
	} else
		for (i = 0; i < num_values; i++)
			values[i*stride].i = round(values[i*stride].f);
	return CCS_SUCCESS;
}

static ccs_result_t
_ccs_distribution_normal_strided_samples(_ccs_distribution_data_t *data,
                                         ccs_rng_t                 rng,
                                         size_t                    num_values,
                                         size_t                    stride,
                                         ccs_numeric_t            *values) {
	_ccs_distribution_normal_data_t *d = (_ccs_distribution_normal_data_t *)data;
	const ccs_numeric_type_t  data_type   = d->common_data.data_types[0];
	const ccs_scale_type_t scale_type     = d->scale_type;
	const ccs_numeric_t    quantization   = d->quantization;
	const ccs_float_t      mu             = d->mu;
	const ccs_float_t      sigma          = d->sigma;
	const int              quantize       = d->quantize;
	gsl_rng *grng;
	CCS_VALIDATE(ccs_rng_get_gsl_rng(rng, &grng));
	if (data_type == CCS_NUM_FLOAT)
		return _ccs_distribution_normal_strided_samples_float(grng, scale_type,
		                                                      quantization.f, mu,
		                                                      sigma, quantize,
		                                                      num_values, stride,
		                                                      (ccs_float_t*) values);
	else
		return _ccs_distribution_normal_strided_samples_int(grng, scale_type,
		                                                    quantization.i, mu,
		                                                    sigma, quantize,
		                                                    num_values, stride, values);
}

static ccs_result_t
_ccs_distribution_normal_soa_samples(_ccs_distribution_data_t  *data,
                                     ccs_rng_t                  rng,
                                     size_t                     num_values,
                                     ccs_numeric_t            **values) {
	if (*values)
		return _ccs_distribution_normal_samples(data, rng, num_values, *values);
	return CCS_SUCCESS;
}

extern ccs_result_t
ccs_create_normal_distribution(ccs_numeric_type_t  data_type,
                               ccs_float_t         mu,
                               ccs_float_t         sigma,
                               ccs_scale_type_t    scale_type,
                               ccs_numeric_t       quantization,
                               ccs_distribution_t *distribution_ret) {
	CCS_CHECK_PTR(distribution_ret);
	if (data_type != CCS_NUM_FLOAT && data_type != CCS_NUM_INTEGER)
		return -CCS_INVALID_TYPE;
	if (scale_type != CCS_LINEAR && scale_type != CCS_LOGARITHMIC)
		return -CCS_INVALID_SCALE;
	if (data_type == CCS_NUM_INTEGER && quantization.i < 0 )
		return -CCS_INVALID_VALUE;
	if (data_type == CCS_NUM_FLOAT && quantization.f < 0.0 )
		return -CCS_INVALID_VALUE;
	uintptr_t mem = (uintptr_t)calloc(1, sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_normal_data_t) + sizeof(ccs_numeric_type_t));

	if (!mem)
		return -CCS_OUT_OF_MEMORY;
	ccs_distribution_t distrib = (ccs_distribution_t)mem;
	_ccs_object_init(&(distrib->obj), CCS_DISTRIBUTION, (_ccs_object_ops_t *)&_ccs_distribution_normal_ops);
        _ccs_distribution_normal_data_t * distrib_data = (_ccs_distribution_normal_data_t *)(mem + sizeof(struct _ccs_distribution_s));
	distrib_data->common_data.data_types    = (ccs_numeric_type_t *)(mem + sizeof(struct _ccs_distribution_s) + sizeof(_ccs_distribution_normal_data_t));
	distrib_data->common_data.type          = CCS_NORMAL;
	distrib_data->common_data.dimension     = 1;
	distrib_data->common_data.data_types[0] = data_type;
	distrib_data->scale_type                = scale_type;
	distrib_data->quantization              = quantization;
	distrib_data->mu                        = mu;
	distrib_data->sigma                     = sigma;
	if (data_type == CCS_NUM_FLOAT) {
		if (quantization.f != 0.0)
			distrib_data->quantize = 1;
	} else {
		if (quantization.i != 0)
			distrib_data->quantize = 1;
	}
	distrib->data = (_ccs_distribution_data_t *)distrib_data;
	*distribution_ret = distrib;
	return CCS_SUCCESS;
}

extern ccs_result_t
ccs_normal_distribution_get_parameters(ccs_distribution_t  distribution,
                                       ccs_float_t        *mu_ret,
                                       ccs_float_t        *sigma_ret,
                                       ccs_scale_type_t   *scale_type_ret,
                                       ccs_numeric_t      *quantization_ret) {
	CCS_CHECK_DISTRIBUTION(distribution, CCS_NORMAL);
	if (!mu_ret && !sigma_ret && !scale_type_ret && !quantization_ret)
		return -CCS_INVALID_VALUE;
	_ccs_distribution_normal_data_t * data = (_ccs_distribution_normal_data_t *)distribution->data;

	if (mu_ret)
		*mu_ret = data->mu;
	if (sigma_ret)
		*sigma_ret = data->sigma;
	if (scale_type_ret)
		*scale_type_ret = data->scale_type;
	if (quantization_ret)
		*quantization_ret = data->quantization;
	return CCS_SUCCESS;
}

ext: c

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

#define qh_QHimport
#include <libqhull_r/qhull_ra.h>

void print_summary(qhT *qh);
double findDelaunay(qhT *qh, int dim, coordT *TargetCoord);
double findNearestVertices(qhT *qh, int dim, coordT *TargetCoord);
double findNearestCentrum(qhT *qh, int dim, coordT *TargetCoord);

void obtain_convex_h_qhull(int n_data, double *ef2cvh, int *x, double *y, int NCatDim, int *nCat_min, int *nCat_max, double *E0, char *qhullflags) {

    int i, j, k, l, nid, Ndim;
    double tmp;

    Ndim = NCatDim - 1;     // because nV is dependent on other nCats
    printf("Convex hull is created using qhull.\n");
    printf("Input data\n");
    printf("---------------------------\n");
    printf("  nL   nN   nM   nC   nV    Ef\n");
    for (i=0;i<n_data;i++) {
        for (j=0;j<NCatDim;j++)
            printf("%4d ",x[NCatDim*i+j]);
        printf("  % .4f\n",y[i]);
    }
    printf("---------------------------\n");

    // Convex Hull by calling qhull
    int DIM;
    int SIZEcube, SIZEdiamond, TOTpoints;
    int dim;            /* dimension of points */
    int numpoints;      /* number of points */
    DIM = Ndim + 1;     /* because of ef as one additional dimension */
    dim = DIM;
    TOTpoints = n_data;
    numpoints= n_data;

//    coordT points[(DIM+1)*TOTpoints]; /* array of coordinates for each point */
    coordT points[DIM*TOTpoints]; /* array of coordinates for each point */
    coordT *rows[TOTpoints];
    boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
    char flags[250];          /* option flags for qhull, see qh-quick.htm */
    FILE *outfile= stdout;    /* output from qh_produce_output()
                               use NULL to skip qh_produce_output() */
    FILE *errfile= stderr;    /* error messages from qhull code */
    int exitcode;             /* 0 if no error from qhull */
    facetT *facet;            /* set by FORALLfacets */
    int curlong, totlong;     /* memory remaining after qh_memfreeshort */

    qhT qh_qh;                /* Qhull's data structure.  First argument of most calls */
    qhT *qh= &qh_qh;

    QHULL_LIB_CHECK

    qh_zero(qh, errfile);

    printf( "\ncompute %d-d convexhull, options:%s\n", dim, qhullflags);
//    sprintf(flags, "qhull o Qx ");     // Error generated
    sprintf(flags, "qhull o %s ",qhullflags);
//    sprintf(flags, "qhull o QJ ");
//    sprintf(flags, "qhull o C-0 ");     // Error generated
//    sprintf(flags, "qhull o C-n ");     // Error generated
//    sprintf(flags, "qhull o Cn ");     // Error generated
//    sprintf(flags, "qhull o Qt ");     // Error generated
//    sprintf(flags, "qhull o Qx C-n ");     // Error generated
    for (i=0;i<n_data;i++) {    //size of points is n_data by (Ndim+1)
        for (j=0;j<Ndim;j++)    //size of x is n_data by Ndim
            points[dim*i+j] = x[NCatDim*i+j];
        points[dim*i+Ndim] = y[i];
    }
    for (i=0;i<numpoints;i++)   //size of rows is n_data x (Ndim+1)
        rows[i]= points+dim*i;
//    for (i=numpoints; i--; )
//        rows[i]= points+dim*i;
    qh_printmatrix(qh, outfile, "input", rows, numpoints, dim);
    exitcode= qh_new_qhull(qh, dim, numpoints, points, ismalloc,
                           flags, outfile, errfile);
    if (!exitcode) {                  /* if no error */
        print_summary(qh);
        printf( "\nThe relative Ef w.r.t. the convex hull:\n", dim);
        exitcode= setjmp(qh->errexit);
        if (!exitcode) {
            /* Trap Qhull errors in findDelaunay().  Without the setjmp(), Qhull
             will exit() after reporting an error */
            qh->NOerrexit= False;
            for (i=0;i<numpoints;i++) {
//                ef2cvh[i] = findDelaunay(qh, DIM, rows[i]);
//                ef2cvh[i] = findNearestVertices(qh, DIM, rows[i]);
                ef2cvh[i] = findNearestCentrum(qh, DIM, rows[i]);
                printf("ef2cvh[%d] = % .4e\n", i, ef2cvh[i]);
            }
        }
        qh->NOerrexit= True;
    }
    qh_freeqhull(qh, !qh_ALL);                 /* free long memory */
    qh_memfreeshort(qh, &curlong, &totlong);  /* free short memory and memory allocator */
    if (curlong || totlong)
        fprintf(errfile, "qhull internal warning (user_eg, #2): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
}

void print_summary(qhT *qh) {
    facetT *facet;
    vertexT *vertex;
    int k;

    printf("\n%d facets with normals:\n", qh->num_facets);
    FORALLfacets {
        for (k=0; k < qh->hull_dim; k++) {
            printf("%6.2g ", facet->normal[k]);
//            printf("", facet->vertices);
        }
        printf("\n");
    }
    printf("\n%d vertices:\n",qh->num_vertices);
    FORALLvertices {
        qh_printvertex(qh, qh->fout, vertex);
    }
}

double findDelaunay(qhT *qh, int dim, coordT *TargetCoord) {
    int j, k;
    coordT tmppoint[ 100];
    boolT isoutside;
    realT bestdist;
    facetT *facet;
    vertexT *vertex, **vertexp;

    for (k= 0; k < dim; k++)
        tmppoint[k]= TargetCoord[k];
    qh_setdelaunay(qh, dim+1, 1, tmppoint);
    facet= qh_findbestfacet(qh, tmppoint, qh_ALL, &bestdist, &isoutside);
    if (facet->tricoplanar) {
        fprintf(stderr, "findDelaunay: not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).\n",
                facet->id);
        qh_errexit(qh, qh_ERRqhull, facet, NULL);
    }

    for (k= 0; k < dim-1; k++)
        printf("% .3f ",tmppoint[k]);
    printf("% .3f, ",tmppoint[dim-1]);
    printf("facet id: %d\n",facet->id);
/*    FOREACHvertex_(facet->vertices) {
        for (k=0; k < dim; k++)
            printf("%5.3f ", vertex->point[k]);
        printf("vertex id: %d, point id: %d\n", vertex->id, qh_pointid(qh, vertex->point));
    }*/
    FOREACHvertex_(facet->vertices) {
        printf("vertex id: %d, point id: %d\n", vertex->id, qh_pointid(qh, vertex->point));
    }

    gsl_matrix * A = gsl_matrix_alloc (dim,dim);
    gsl_matrix * V = gsl_matrix_alloc (dim,dim);
    gsl_vector * S = gsl_vector_alloc (dim);
    gsl_vector * work = gsl_vector_alloc (dim);
    gsl_vector * b = gsl_vector_alloc (dim);
    gsl_vector * x = gsl_vector_alloc (dim);

    j= 0;
    FOREACHvertex_(facet->vertices) {
        for (k= 0; k < dim; k++)
            gsl_matrix_set (A,k,j, vertex->point[k]);
        gsl_vector_set(b,j,tmppoint[j]);
        j++;
    }
    double ef[dim];
    // dismiss the last row corresponding to Ef
    // & apply the constriction on the coefficients sum_x_i = 1
    for (k= 0; k < dim; k++) {
        ef[k] = gsl_matrix_get (A,dim-1,k);
        gsl_matrix_set (A,dim-1,k,1.0);
    }
    gsl_vector_set(b,dim-1,2.0);

    gsl_linalg_SV_decomp(A, V, S, work);
    gsl_linalg_SV_solve (A, V, S, b, x);

    double ef_cvh = 0.0;
    for (k= 0; k < dim; k++)
        ef_cvh += gsl_vector_get(x,k)*ef[k];

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(x);

    return(TargetCoord[dim-1]-ef_cvh);
}

double findNearestVertices(qhT *qh, int dim, coordT *TargetCoord) {
    int i, j, k;
    int nid;
    double ef_cvh;
    double tmp;
    double tmppoint[dim];
    boolT isoutside;
    realT bestdist;
    facetT *facet;
    vertexT *vertex, **vertexp;

    nid = qh->num_vertices;
    for (k= 0; k < dim; k++)
        tmppoint[k]= TargetCoord[k];

    double verticespoint[nid][dim];
    j = 0;
    FORALLvertices {
        for (k= 0; k < dim; k++)
            verticespoint[j][k] = vertex->point[k];
        j++;
    }

    double dij2[nid];
    int ids[nid];
    for (j= 0; j< nid; j++) {
        tmp = 0.0;
        for (k= 0; k < dim-1; k++) {    //ef is excluded for configurational coordinate
            tmp += pow(verticespoint[j][k] - tmppoint[k],2);
        }
        dij2[j] = tmp;
        ids[j] = j;
    }
    sort_array_ids(dij2,ids,nid);

    gsl_matrix * A = gsl_matrix_alloc (dim,dim);
    gsl_matrix * V = gsl_matrix_alloc (dim,dim);
    gsl_vector * S = gsl_vector_alloc (dim);
    gsl_vector * work = gsl_vector_alloc (dim);
    gsl_vector * b = gsl_vector_alloc (dim);
    gsl_vector * x = gsl_vector_alloc (dim);

    if (dij2[0]==0.0)
        ef_cvh = 0.0;
    else {
        for (j= 0; j< dim; j++) {
            for (k= 0; k < dim; k++)
                gsl_matrix_set (A,k,j,verticespoint[ids[j]][k]);
            gsl_vector_set(b,j,tmppoint[j]);
        }
        for (k= 0; k < dim; k++) {      //dismiss the last row corresponding to ef
            gsl_matrix_set (A,dim-1,k,1.0);
        }
        gsl_vector_set(b,dim-1,2.0);    //dismiss the last row corresponding to ef
        gsl_linalg_SV_decomp(A, V, S, work);
        gsl_linalg_SV_solve (A, V, S, b, x);
        tmp = 0.0;
        for (j= 0; j< dim; j++)
            tmp += gsl_vector_get(x,j)*verticespoint[ids[j]][dim-1];
        ef_cvh = tmppoint[dim-1]-tmp;
    }

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(x);

    return(ef_cvh);
}

double findNearestCentrum(qhT *qh, int dim, coordT *TargetCoord) {
    int i, j, k;
    int nf, nv;
    double ef_cvh;
    double tmp;
    double tmppoint[dim];
    boolT isoutside;
    realT bestdist;
    facetT *facet;
    vertexT *vertex, **vertexp;
    pointT *centrum;

    nf = qh->num_facets;
    nv = qh->num_vertices;

    for (k= 0; k < dim; k++)
        tmppoint[k]= TargetCoord[k];

    double facet_normal[nf][dim];
    double facet_center[nf][dim];
    int facet_id[nf];
    int num_lowerhullfacets;

    j = 0;
    FORALLfacets {
        if (facet->normal[dim-1] < 0.0) {//is the slope in ef direction negative?
            for (k= 0; k < dim; k++) {
                facet_normal[j][k] = facet->normal[k];
                if (qh->CENTERtype == qh_AScentrum)
                    centrum= facet->center;
                else
                    centrum= qh_getcentrum(qh, facet);
                facet_center[j][k] = centrum[k];
                facet_id[j] = facet->id;
            }
            j++;
        }
    }
    num_lowerhullfacets = j;

    double dij2[num_lowerhullfacets];
    int ids[num_lowerhullfacets];
    for (j= 0; j< num_lowerhullfacets; j++) {
        tmp = 0.0;
        for (k= 0; k < dim-1; k++) {    //ef is excluded for configurational coordinate
            tmp += pow(facet_center[j][k] - tmppoint[k],2);
        }
        dij2[j] = tmp;
        ids[j] = j;
    }
    sort_array_ids(dij2,ids,num_lowerhullfacets);

    int projectedfacet = facet_id[ids[0]];
    double vertices_point[dim][dim];
    FORALLfacets {
        if (facet->id == projectedfacet) {
            j = 0;
            FOREACHvertex_(facet->vertices) {
                for (k=0; k < dim; k++)
                    vertices_point[j][k] = vertex->point[k];
                j++;
            }
        }
    }

    gsl_matrix * A = gsl_matrix_alloc (dim,dim);
    gsl_matrix * V = gsl_matrix_alloc (dim,dim);
    gsl_vector * S = gsl_vector_alloc (dim);
    gsl_vector * work = gsl_vector_alloc (dim);
    gsl_vector * b = gsl_vector_alloc (dim);
    gsl_vector * x = gsl_vector_alloc (dim);

    for (j= 0; j< dim; j++) {
        for (k= 0; k < dim; k++)
            gsl_matrix_set (A,k,j,vertices_point[j][k]);
        gsl_vector_set(b,j,tmppoint[j]);
    }
    for (k= 0; k < dim; k++) {      //dismiss the last row corresponding to ef
        gsl_matrix_set (A,dim-1,k,1.0);
    }
    gsl_vector_set(b,dim-1,1.0);    //dismiss the last row corresponding to ef

    printf("A:\n");
    for (j= 0; j< dim; j++) {
        for (k= 0; k < dim; k++)
            printf("%g ", gsl_matrix_get (A,j,k));
        printf("\n");
    }
    printf("b:\n");
    for (j= 0; j< dim; j++)
        printf("%g\n", gsl_vector_get(b,j));

    gsl_linalg_SV_decomp(A, V, S, work);
    gsl_linalg_SV_solve (A, V, S, b, x);
    printf("solution:\n");
    for (j= 0; j< dim; j++)
        printf("%g\n", gsl_vector_get(x,j));
    // need to check if the components of x is in [0,2]
    for (j= 0; j< dim; j++) {
        tmp = gsl_vector_get(x,j);
        if (tmp < -0.1 || tmp > 1.0) {
            printf("\n%d-component of x is %f, out of range\n", j, tmp);
            exit(-1);
        }
    }
    tmp = 0.0;
    for (j= 0; j< dim; j++)
        tmp += gsl_vector_get(x,j)*vertices_point[j][dim-1];
    ef_cvh = tmppoint[dim-1]-tmp;

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(x);

    return(ef_cvh);
}


text:
#ifndef _PROBABILITY_H_
#define _PROBABILITY_H_

#include <gsl/gsl_rng.h>
// Non-uniform distributions:
//There are a HUGE variety of distributions to choose from, which are
//  documented at the gsl website.
#include <gsl/gsl_randist.h>
// Seeding:
#include <sys/time.h>

//When using ANY gsl code, must add -lgsl -lgslcblas [-lm ???] to compile lines
//USE_LINK_FLAG -lgsl -lgslcblas

//RandGen is a virtual class.  Use it in function definitions,
//  but when DECLARING a random number generator, use one of its derived
//  classes, such as gsl_RandGen.
#define UNUSED_ARGUMENT(x) (void)x

class RandGen
{
 public:
  //Create a random number generator with Seed generated by system clock:
  RandGen() {}
  //Create a random number generator, and specify seed
  RandGen(const long & seed) {UNUSED_ARGUMENT(seed);}
  virtual ~RandGen() {}

  //uniform double on [0, 1):
  virtual double uniform(void) = 0;
  //uniform double on (0, 1):
  virtual double openUniform(void) = 0;
  //unsigned long int on [0, 1, 2, ..., n - 1]:
  virtual unsigned long unsignedInt(const unsigned long & n) = 0;
  // bool, true with probability p, false with probability 1-p:
  virtual bool bernoulli(const double & p) = 0;
  //exponentially distributed, p(x) ~ exp(-x/mu)
  virtual double exponential(const double & mu) = 0;
  //gaussian double with zero mean and standard deviation sigma:
  virtual double gaussian(const double & sigma) = 0;
  //poisson unsigned int with mean Mu:
  virtual unsigned int poisson(const double & mu) = 0;
  //binomial unsigned int for n trials with probability P of success
  virtual unsigned int binomial(const double & p, unsigned intnN) = 0;
};

class GslRandGen : public RandGen
{
 public:
  GslRandGen();
  GslRandGen(const long & seed);
  virtual ~GslRandGen();

  //uniform double on [0, 1):
  virtual double uniform(void);
  //uniform double on (0, 1):
  virtual double openUniform(void);
  //unsigned long int on [0, 1, 2, ..., n - 1]:
  virtual unsigned long unsignedInt(const unsigned long & n);

  // bool, true with probability p, false with probability 1-p:
  bool bernoulli(const double & p);

  //exponentially distributed, p(x) ~ exp(-x/mu)
  double exponential(const double & mu);
  //gaussian double with zero mean and standard deviation sigma:
  double gaussian(const double & sigma);
  //poisson unsigned int with mean Mu:
  unsigned int poisson(const double & mu);
  //binomial unsigned int for n trials with probability P of success
  unsigned int binomial(const double & p, unsigned intnN);

 private:
  gsl_rng* r;
  static bool firstInit;
};

#endif

ext: h


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 *  \file allvars.h
 *  \brief header file of allvars.c.
 */

#ifndef _ALLVARS_H
#define _ALLVARS_H

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <mpi.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_interp.h>

/*! \def GRAVITY
 *  \brief Gravitational constant.

 *  \def SOLAR_MASS
 *  \brief Solar mass.

 *  \def C
 *  \brief speed of light.

 *  \def SEC_PER_YEAR
 *  \brief seconds per year.

 *  \def CM_PER_LD
 *  \brief light day in centimeter.

 *  \def PI
 *  \brief Pi.

 *  \def BRAINS_MAX_STR_LENGTH
 *  \brief maximum string length.
 */
#define GRAVITY       6.672e-8
#define SOLAR_MASS    1.989e33
#define C             2.9979e10
#define SEC_PER_YEAR  3.155e7
#define CM_PER_LD    (C*8.64e4)
#define CM_PER_PC    (3.08567758e+18)
#define CM_PER_MPC   (3.08567758e+24)

#define PI            M_PI
#define BRAINS_MAX_STR_LENGTH  (256)

#define EPS (DBL_MIN)  /* epsilon of the machine as in Matlab */

/* for PARDICT */
#define MAXTAGS 300
#define DOUBLE 1
#define STRING 2
#define INT 3

enum PRIOR_TYPE {GAUSSIAN=1, UNIFORM=2};

/* variables for MPICH */
extern int thistask, totaltask, namelen;
extern int roottask;
extern char proc_name[MPI_MAX_PROCESSOR_NAME];

/*! \struct PARSET
 *  \brief the configuration parameters.
 */
typedef struct
{
  char param_file[BRAINS_MAX_STR_LENGTH];
  char continuum_file[BRAINS_MAX_STR_LENGTH],
       line_file[BRAINS_MAX_STR_LENGTH],
       line2d_file[BRAINS_MAX_STR_LENGTH],
       pcon_out_file[BRAINS_MAX_STR_LENGTH],
       pline_out_file[BRAINS_MAX_STR_LENGTH],
       pline2d_out_file[BRAINS_MAX_STR_LENGTH],
       pline2d_data_out_file[BRAINS_MAX_STR_LENGTH],
       cloud_out_file[BRAINS_MAX_STR_LENGTH],
       tran_out_file[BRAINS_MAX_STR_LENGTH],
       tran2d_out_file[BRAINS_MAX_STR_LENGTH],
       tran2d_data_out_file[BRAINS_MAX_STR_LENGTH];

  char file_dir[BRAINS_MAX_STR_LENGTH];

  int flag_dim;

  int flag_trend, flag_trend_diff;

  int n_con_recon, n_line_recon, n_vel_recon;

  int n_tau;
  double rcloud_max, time_back;

  int flag_blrmodel;

  int n_cloud_per_task, n_vel_per_cloud;

  int flag_save_clouds;

  int flag_InstRes;
  double InstRes, InstRes_err;
  char file_instres[BRAINS_MAX_STR_LENGTH];

  int num_particles;

  char str_par_fix[BRAINS_MAX_STR_LENGTH], str_par_fix_val[BRAINS_MAX_STR_LENGTH];

  int flag_narrowline;
  double flux_narrowline, width_narrowline, shift_narrowline;
  double flux_narrowline_err, width_narrowline_err, shift_narrowline_err;

  int flag_linecenter;
  double linecenter_err;

  int flag_postprc;
  int flag_restart;
  int flag_sample_info;
  int flag_temp;
  int flag_exam_prior;
  int flag_rng_seed, rng_seed;
  int flag_fixvar;
  int flag_nonlinear;

  int flag_force_update;
  int flag_force_run;

  int flag_con_sys_err, flag_line_sys_err;

  double temperature;

  int flag_help, flag_end;
  int flag_para_name;

  double redshift;
  double linecenter;
#ifdef SA
  char sa_file[BRAINS_MAX_STR_LENGTH];
  int flag_sa_blrmodel, flag_sa_par_mutual;
  double sa_linecenter;
  double sa_InstRes;
  char sa_str_par_fix[BRAINS_MAX_STR_LENGTH], sa_str_par_fix_val[BRAINS_MAX_STR_LENGTH];
  int n_sa_vel_recon, n_sa_base_recon;
#endif
}PARSET;
extern PARSET parset;

typedef struct
{
  int id;
  void *addr;
  char tag[50];
  int isset;
}PARDICT;
extern PARDICT *pardict;
extern int num_pardict;

extern double VelUnit, C_Unit;

extern int n_con_data, n_line_data, n_vel_data, n_vel_data_incr, n_vel_data_ext, n_con_max;
extern double *Tcon_data, *Fcon_data,  *Fcerrs_data;
extern double *Tline_data, *Fline_data, *Flerrs_data;
extern double *Vline_data, *Fline2d_data, *Flerrs2d_data, *Wline_data;
extern double *Vline_data_ext, *Wline_data_ext;
extern double con_scale, line_scale;

extern double con_error_mean, line_error_mean, line_error_mean_sq;

extern char dnest_options_file[BRAINS_MAX_STR_LENGTH];

extern int which_parameter_update, which_particle_update;  // which parameter and particle to be updated
extern int which_level_update;
extern double *limits;  // external from dnest

/* continuum reconstruction */
extern int nq;
extern double *Tcon, *Fcon, *Fcerrs, *Fcon_rm;
extern double Tcon_min, Tcon_max;
extern double *PSmat, *USmat, *PSmat_data;
extern double *PCmat_data, *IPCmat_data, *PQmat, *PEmat1, *PEmat2;
extern double *workspace, *workspace_uv;
extern double *var_param, *var_param_std;
extern double logz_con;
extern double *Larr_data, *Larr_rec;
extern double *pow_Tcon_data;

extern double Tspan_data, Tspan_data_con, Tcad_data, Tmed_data;

/* line reconstruction */
extern double *Fline_at_data;
extern double *Tline, *Fline, *Flerrs;
extern double Tline_min, Tline_max;
extern double logz_line;

/* line 2d reconstruction */
extern double *Fline2d_at_data;
extern double *Fline2d, *Flerrs2d;
extern double logz_line2d;

extern int BLRmodel_size;
extern int *par_fix, npar_fix;
extern double *par_fix_val;

extern int num_params_blr_tot;
extern int num_params, num_params_blr, num_params_blr_model, num_params_var, num_params_difftrend, num_params_nlr, num_params_res;
extern int num_params_drw, num_params_trend, num_params_resp;
extern int num_params_linecenter;

extern double **blr_range_model, **par_range_model;
extern int *par_prior_model;  /* prior type for model parameters */
extern double **par_prior_gaussian; /* center and std of Gaussian priors */
extern double var_range_model[15][2];
extern double nlr_range_model[3][2];
extern int nlr_prior_model[3];
extern double mass_range[2];
extern double sys_err_line_range[2]; /* range for systematic error of line */
extern double resp_range[2][2];
extern int idx_resp, idx_difftrend, idx_linecenter;

extern double *instres_epoch, *instres_err_epoch;

/* transfer function / velocity-delay map */
extern double *TransTau, *TransV, *TransW, *Trans1D, *Trans2D_at_veldata, *Trans2D;
extern double rcloud_min_set, rcloud_max_set, time_back_set;

extern double **Fcon_rm_particles, **Fcon_rm_particles_perturb;
extern double *prob_con_particles, *prob_con_particles_perturb;

extern int force_update;
extern double **TransTau_particles, **TransTau_particles_perturb;
extern double **Trans1D_particles, **Trans1D_particles_perturb;
extern double **Trans2D_at_veldata_particles, **Trans2D_at_veldata_particles_perturb;
extern double **Fline_at_data_particles, **Fline_at_data_particles_perturb;

extern double *clouds_tau, *clouds_weight, *clouds_vel;

extern FILE *fcloud_out;
extern int icr_cloud_save;

/* GSL */
extern const gsl_rng_type * gsl_T;
extern gsl_rng * gsl_r;

extern gsl_interp_accel *gsl_acc;
extern gsl_interp  *gsl_linear;

#ifdef SA
extern double PhaseFactor;

extern int num_params_rm;
extern int num_params_sa, num_params_sa_blr_model, num_params_sa_extpar, num_params_sa_blr;
extern int n_epoch_sa_data, n_vel_sa_data, n_base_sa_data;
extern double *vel_sa_data, *base_sa_data, *Fline_sa_data, *Flerrs_sa_data, *phase_sa_data, *pherrs_sa_data;
extern double *wave_sa_data;

extern double sa_flux_norm;

extern double *phase_sa, *Fline_sa, *base_sa, *vel_sa, *wave_sa;

extern double *clouds_alpha, *clouds_beta;

extern double **sa_extpar_range, **sa_blr_range_model;

extern int SABLRmodel_size;

extern int *idx_sa_par_mutual, *idx_rm_par_mutual;

extern double *prob_sa_particles, *prob_sa_particles_perturb;

extern double sa_phase_error_mean, sa_line_error_mean;

extern double *workspace_phase;
#endif
#endif

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS