/*
 * Academic License - for use in teaching, academic research, and meeting
 * course requirements at degree granting institutions only.  Not for
 * government, commercial, or other organizational use.
 *
 * moment_gen_function_cg.c
 *
 * Code generation for function 'moment_gen_function_cg'
 *
 */

/* Include files */
#include "moment_gen_function_cg.h"
#include "ODEsystem_2F_Gaussian_cg.h"
#include "digital_option_price_cg_data.h"
#include "digital_option_price_cg_emxutil.h"
#include "digital_option_price_cg_types.h"
#include "eml_int_forloop_overflow_check.h"
#include "rt_nonfinite.h"
#include "mwmathutil.h"

/* Variable Definitions */
static emlrtRSInfo s_emlrtRSI = {
    299,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo t_emlrtRSI = {
    300,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo u_emlrtRSI = {
    301,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo v_emlrtRSI = {
    302,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo kc_emlrtRSI = {
    318,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtDCInfo emlrtDCI = {
    276,                      /* lineNo */
    21,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    1 /* checkKind */
};

static emlrtDCInfo b_emlrtDCI = {
    280,                      /* lineNo */
    21,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    1 /* checkKind */
};

static emlrtRTEInfo emlrtRTEI = {
    298,                      /* lineNo */
    25,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtBCInfo g_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    295,                      /* lineNo */
    19,                       /* colNo */
    "t",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo h_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    299,                      /* lineNo */
    54,                       /* colNo */
    "t",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo i_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    296,                      /* lineNo */
    21,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo j_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    299,                      /* lineNo */
    66,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo k_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    300,                      /* lineNo */
    66,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo l_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    301,                      /* lineNo */
    66,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo m_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    302,                      /* lineNo */
    66,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo n_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    303,                      /* lineNo */
    35,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo o_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    303,                      /* lineNo */
    24,                       /* colNo */
    "w",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtECInfo yb_emlrtECI = {
    -1,                       /* nDims */
    310,                      /* lineNo */
    17,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtECInfo ac_emlrtECI = {
    -1,                       /* nDims */
    311,                      /* lineNo */
    17,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtECInfo bc_emlrtECI = {
    -1,                       /* nDims */
    312,                      /* lineNo */
    17,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtECInfo cc_emlrtECI = {
    -1,                       /* nDims */
    313,                      /* lineNo */
    17,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtDCInfo c_emlrtDCI = {
    47,                       /* lineNo */
    9,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    1 /* checkKind */
};

static emlrtDCInfo d_emlrtDCI = {
    47,                       /* lineNo */
    9,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    4 /* checkKind */
};

static emlrtBCInfo p_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    318,                      /* lineNo */
    24,                       /* colNo */
    "alpha",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo q_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    318,                      /* lineNo */
    39,                       /* colNo */
    "beta1",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo r_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    304,                      /* lineNo */
    22,                       /* colNo */
    "t",                      /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtDCInfo f_emlrtDCI = {
    47,                       /* lineNo */
    23,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    4 /* checkKind */
};

static emlrtDCInfo g_emlrtDCI = {
    47,                       /* lineNo */
    23,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    1 /* checkKind */
};

static emlrtBCInfo s_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    284,                      /* lineNo */
    41,                       /* colNo */
    "phi1",                   /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo t_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    285,                      /* lineNo */
    41,                       /* colNo */
    "phi2",                   /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo u_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    292,                      /* lineNo */
    44,                       /* colNo */
    "phi1",                   /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo v_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    292,                      /* lineNo */
    76,                       /* colNo */
    "phi2",                   /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo w_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    310,                      /* lineNo */
    25,                       /* colNo */
    "alpha",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo x_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    284,                      /* lineNo */
    26,                       /* colNo */
    "beta1",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo y_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    311,                      /* lineNo */
    25,                       /* colNo */
    "beta1",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtECInfo dc_emlrtECI = {
    2,                        /* nDims */
    318,                      /* lineNo */
    18,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtBCInfo ab_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    285,                      /* lineNo */
    26,                       /* colNo */
    "beta2",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo bb_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    312,                      /* lineNo */
    25,                       /* colNo */
    "beta2",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo cb_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    313,                      /* lineNo */
    25,                       /* colNo */
    "beta3",                  /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo db_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    299,                      /* lineNo */
    84,                       /* colNo */
    "uu",                     /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo eb_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    300,                      /* lineNo */
    84,                       /* colNo */
    "uu",                     /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo fb_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    301,                      /* lineNo */
    84,                       /* colNo */
    "uu",                     /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtBCInfo gb_emlrtBCI = {
    -1,                       /* iFirst */
    -1,                       /* iLast */
    302,                      /* lineNo */
    84,                       /* colNo */
    "uu",                     /* aName */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m", /* pName
                                                                         */
    0 /* checkKind */
};

static emlrtRTEInfo ac_emlrtRTEI = {
    47,                       /* lineNo */
    9,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo bc_emlrtRTEI = {
    48,                       /* lineNo */
    1,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo cc_emlrtRTEI = {
    49,                       /* lineNo */
    1,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo dc_emlrtRTEI = {
    276,                      /* lineNo */
    9,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo ec_emlrtRTEI = {
    280,                      /* lineNo */
    9,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo fc_emlrtRTEI = {
    310,                      /* lineNo */
    40,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo gc_emlrtRTEI = {
    310,                      /* lineNo */
    30,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo hc_emlrtRTEI = {
    311,                      /* lineNo */
    40,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo ic_emlrtRTEI = {
    311,                      /* lineNo */
    30,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo jc_emlrtRTEI = {
    312,                      /* lineNo */
    40,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo kc_emlrtRTEI = {
    312,                      /* lineNo */
    30,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo lc_emlrtRTEI = {
    1,                        /* lineNo */
    15,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo mc_emlrtRTEI = {
    1,                        /* lineNo */
    22,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo nc_emlrtRTEI = {
    1,                        /* lineNo */
    29,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo oc_emlrtRTEI = {
    310,                      /* lineNo */
    17,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo qc_emlrtRTEI = {
    50,                       /* lineNo */
    1,                        /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo rc_emlrtRTEI = {
    318,                      /* lineNo */
    33,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo sc_emlrtRTEI = {
    318,                      /* lineNo */
    18,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo tc_emlrtRTEI = {
    318,                      /* lineNo */
    54,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo uc_emlrtRTEI = {
    313,                      /* lineNo */
    30,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRTEInfo vc_emlrtRTEI = {
    313,                      /* lineNo */
    43,                       /* colNo */
    "moment_gen_function_cg", /* fName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pName
                                                                        */
};

static emlrtRSInfo dd_emlrtRSI = {
    313,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo ed_emlrtRSI = {
    310,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo fd_emlrtRSI = {
    311,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

static emlrtRSInfo gd_emlrtRSI = {
    312,                      /* lineNo */
    "moment_gen_function_cg", /* fcnName */
    "/Users/lazarus/Dropbox/Work/Research/Ned/VolatilityFinance/"
    "ReplicationCodeManual/Code/Temp/Toolbox/moment_gen_function_cg.m" /* pathName
                                                                        */
};

/* Function Declarations */
static void b_plus(const emlrtStack *sp, emxArray_creal_T *in1,
                   const emxArray_creal_T *in2);

static void binary_expand_op_12(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3);

static void binary_expand_op_13(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3);

static void binary_expand_op_14(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3);

static void binary_expand_op_21(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_creal_T *in2);

static void binary_expand_op_22(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2);

/* Function Definitions */
static void b_plus(const emlrtStack *sp, emxArray_creal_T *in1,
                   const emxArray_creal_T *in2)
{
  emxArray_creal_T *b_in1;
  const creal_T *in2_data;
  creal_T *b_in1_data;
  creal_T *in1_data;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_1;
  int32_T stride_1_1;
  in2_data = in2->data;
  in1_data = in1->data;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  emxInit_creal_T(sp, &b_in1, 2, &sc_emlrtRTEI);
  i = b_in1->size[0] * b_in1->size[1];
  b_in1->size[0] = 1;
  if (in2->size[1] == 1) {
    loop_ub = in1->size[1];
  } else {
    loop_ub = in2->size[1];
  }
  b_in1->size[1] = loop_ub;
  emxEnsureCapacity_creal_T(sp, b_in1, i, &sc_emlrtRTEI);
  b_in1_data = b_in1->data;
  stride_0_1 = (in1->size[1] != 1);
  stride_1_1 = (in2->size[1] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    int32_T i2;
    i1 = i * stride_0_1;
    i2 = i * stride_1_1;
    b_in1_data[i].re = in1_data[i1].re + in2_data[i2].re;
    b_in1_data[i].im = in1_data[i1].im + in2_data[i2].im;
  }
  i = in1->size[0] * in1->size[1];
  in1->size[0] = 1;
  in1->size[1] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &sc_emlrtRTEI);
  in1_data = in1->data;
  for (i = 0; i < loop_ub; i++) {
    in1_data[i] = b_in1_data[i];
  }
  emxFree_creal_T(sp, &b_in1);
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
}

static void binary_expand_op_12(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3)
{
  const creal_T *in3_data;
  creal_T *in1_data;
  const real_T *in2_data;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_0;
  int32_T stride_1_0;
  in3_data = in3->data;
  in2_data = in2->data;
  if (in3->size[0] == 1) {
    loop_ub = in2->size[1];
  } else {
    loop_ub = in3->size[0];
  }
  i = in1->size[0];
  in1->size[0] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &kc_emlrtRTEI);
  in1_data = in1->data;
  stride_0_0 = (in2->size[1] != 1);
  stride_1_0 = (in3->size[0] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    i1 = i * stride_1_0;
    in1_data[i].re = in2_data[6 * (i * stride_0_0) + 4] + in3_data[i1].re;
    in1_data[i].im = in3_data[i1].im;
  }
}

static void binary_expand_op_13(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3)
{
  const creal_T *in3_data;
  creal_T *in1_data;
  const real_T *in2_data;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_0;
  int32_T stride_1_0;
  in3_data = in3->data;
  in2_data = in2->data;
  if (in3->size[0] == 1) {
    loop_ub = in2->size[1];
  } else {
    loop_ub = in3->size[0];
  }
  i = in1->size[0];
  in1->size[0] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &ic_emlrtRTEI);
  in1_data = in1->data;
  stride_0_0 = (in2->size[1] != 1);
  stride_1_0 = (in3->size[0] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    i1 = i * stride_1_0;
    in1_data[i].re = in2_data[6 * (i * stride_0_0) + 2] + in3_data[i1].re;
    in1_data[i].im = in3_data[i1].im;
  }
}

static void binary_expand_op_14(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2,
                                const emxArray_creal_T *in3)
{
  const creal_T *in3_data;
  creal_T *in1_data;
  const real_T *in2_data;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_0;
  int32_T stride_1_0;
  in3_data = in3->data;
  in2_data = in2->data;
  if (in3->size[0] == 1) {
    loop_ub = in2->size[1];
  } else {
    loop_ub = in3->size[0];
  }
  i = in1->size[0];
  in1->size[0] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &gc_emlrtRTEI);
  in1_data = in1->data;
  stride_0_0 = (in2->size[1] != 1);
  stride_1_0 = (in3->size[0] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    i1 = i * stride_1_0;
    in1_data[i].re = in2_data[6 * (i * stride_0_0)] + in3_data[i1].re;
    in1_data[i].im = in3_data[i1].im;
  }
}

static void binary_expand_op_21(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_creal_T *in2)
{
  emxArray_creal_T *c_in2;
  const creal_T *in2_data;
  creal_T *b_in2_data;
  creal_T *in1_data;
  int32_T b_in2;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_1;
  int32_T stride_1_1;
  in2_data = in2->data;
  in1_data = in1->data;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  b_in2 = in2->size[0];
  emxInit_creal_T(sp, &c_in2, 2, &sc_emlrtRTEI);
  i = c_in2->size[0] * c_in2->size[1];
  c_in2->size[0] = 1;
  if (in1->size[1] == 1) {
    loop_ub = in2->size[1];
  } else {
    loop_ub = in1->size[1];
  }
  c_in2->size[1] = loop_ub;
  emxEnsureCapacity_creal_T(sp, c_in2, i, &sc_emlrtRTEI);
  b_in2_data = c_in2->data;
  stride_0_1 = (in2->size[1] != 1);
  stride_1_1 = (in1->size[1] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    int32_T i2;
    i1 = i * stride_0_1;
    i2 = i * stride_1_1;
    b_in2_data[i].re =
        in2_data[(b_in2 + in2->size[0] * i1) - 1].re + in1_data[i2].re;
    b_in2_data[i].im =
        in2_data[(b_in2 + in2->size[0] * i1) - 1].im + in1_data[i2].im;
  }
  i = in1->size[0] * in1->size[1];
  in1->size[0] = 1;
  in1->size[1] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &sc_emlrtRTEI);
  in1_data = in1->data;
  for (i = 0; i < loop_ub; i++) {
    in1_data[i] = b_in2_data[i];
  }
  emxFree_creal_T(sp, &c_in2);
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
}

static void binary_expand_op_22(const emlrtStack *sp, emxArray_creal_T *in1,
                                const emxArray_real_T *in2)
{
  emxArray_creal_T *b_in2;
  creal_T *b_in2_data;
  creal_T *in1_data;
  const real_T *in2_data;
  int32_T i;
  int32_T loop_ub;
  int32_T stride_0_0;
  int32_T stride_1_0;
  in2_data = in2->data;
  in1_data = in1->data;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  emxInit_creal_T(sp, &b_in2, 1, &uc_emlrtRTEI);
  if (in1->size[0] == 1) {
    loop_ub = in2->size[0];
  } else {
    loop_ub = in1->size[0];
  }
  i = b_in2->size[0];
  b_in2->size[0] = loop_ub;
  emxEnsureCapacity_creal_T(sp, b_in2, i, &uc_emlrtRTEI);
  b_in2_data = b_in2->data;
  stride_0_0 = (in2->size[0] != 1);
  stride_1_0 = (in1->size[0] != 1);
  for (i = 0; i < loop_ub; i++) {
    int32_T i1;
    i1 = i * stride_1_0;
    b_in2_data[i].re = in2_data[i * stride_0_0] + in1_data[i1].re;
    b_in2_data[i].im = in1_data[i1].im;
  }
  i = in1->size[0];
  in1->size[0] = loop_ub;
  emxEnsureCapacity_creal_T(sp, in1, i, &uc_emlrtRTEI);
  in1_data = in1->data;
  for (i = 0; i < loop_ub; i++) {
    in1_data[i] = b_in2_data[i];
  }
  emxFree_creal_T(sp, &b_in2);
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
}

creal_T b_moment_gen_function_cg(const emlrtStack *sp, real_T sv1, real_T sv2,
                                 real_T T, real_T rf, real_T d)
{
  emlrtStack st;
  emxArray_creal_T *alpha;
  emxArray_creal_T *b_r;
  emxArray_creal_T *beta1;
  emxArray_creal_T *beta2;
  emxArray_creal_T *r2;
  emxArray_real_T *t;
  emxArray_real_T *w;
  creal_T rr;
  creal_T *alpha_data;
  creal_T *beta1_data;
  creal_T *beta2_data;
  creal_T *r1;
  creal_T *r3;
  real_T h;
  real_T r;
  real_T *t_data;
  real_T *w_data;
  int32_T b_loop_ub;
  int32_T i;
  int32_T i1;
  int32_T loop_ub;
  int32_T loop_ub_tmp;
  boolean_T b;
  st.prev = sp;
  st.tls = sp->tls;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  /*  moment_gen_function.m computes the Joint Conditional MGF (or the */
  /*  Characteristic Function -CF-) of log(S_T/S_t), V1_T, V2_T, and V3_T */
  /*  corresponding to the model "distr" at points uu, phi1, phi2, and phi3.  */
  /*  */
  /*  Note: */
  /*  CF_t(x,u,phi1,phi2,phi3,T) = */
  /*  E_t[exp(uu*y_T+phi1*V_{1,T}+phi2*V_{2,T}+phi3*V_{3,T})] = */
  /*  exp(alpha(x,T) + beta1(x,T)V_{1,t} + beta1(x,T)V_{2,t} +
   * beta3(x,T)V_{3,t}) */
  /* ==========================================================================
   */
  /*  OUTPUT:  */
  /*  */
  /*            M(uu,phi1,phi2,t,T) = E_t[exp(uu*y_T+ alpha + phi1*V_{1,T} +
   * phi2*V_{2,T})]  */
  /*  INPUT: */
  /*        distr:  model distribution */
  /*        x:      model parameters */
  /*        uu:     evalutation points for y_T (can be a vector) */
  /*        phi1:   evalutation points for V_1 (can be a vector) */
  /*        phi2:   evalutation points for V_2 (can be a vector) */
  /*        phi3:   evalutation points for V_3 (can be a vector) */
  /*        sv1:    initial value of the first factor */
  /*        sv2:    initial value of the second factor */
  /*        sv3:    initial value of the third factor */
  /*        T:      maturity (in years, convention 365 days per year) */
  /*        rf:     interest rate */
  /*        d:      dividend yield */
  /*  */
  /* --------------------------------------------------------------------------
   */
  /*  Nicola Fusari (2015): */
  /*  You can distribute and use the code. Please cite the author if you use */
  /*  the code in in your research. */
  /*   */
  /*  The code used to Solve the Systme of ODEs is from Alain Kapitho */
  /*  (University of Pretoria), available at: */
  /*  http://m2matlabdb.ma.tum.de/rk4_systems.m?MP_ID=413 */
  /*  */
  /* --------------------------------------------------------------------------
   */
  r = muDoubleScalarRound(T * 365.0);
  emxInit_creal_T(sp, &alpha, 1, &lc_emlrtRTEI);
  if (!(r + 1.0 >= 0.0)) {
    emlrtNonNegativeCheckR2012b(r + 1.0, &d_emlrtDCI, (emlrtConstCTX)sp);
  }
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  i = alpha->size[0];
  alpha->size[0] = (int32_T)(r + 1.0);
  emxEnsureCapacity_creal_T(sp, alpha, i, &ac_emlrtRTEI);
  alpha_data = alpha->data;
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  loop_ub = (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub; i++) {
    alpha_data[i].re = 0.0;
    alpha_data[i].im = 0.0;
  }
  emxInit_creal_T(sp, &beta1, 1, &mc_emlrtRTEI);
  loop_ub = alpha->size[0];
  i = beta1->size[0];
  beta1->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta1, i, &bc_emlrtRTEI);
  beta1_data = beta1->data;
  emxInit_creal_T(sp, &beta2, 1, &nc_emlrtRTEI);
  i = beta2->size[0];
  beta2->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta2, i, &cc_emlrtRTEI);
  beta2_data = beta2->data;
  for (i = 0; i < loop_ub; i++) {
    beta1_data[i] = alpha_data[i];
    beta2_data[i] = alpha_data[i];
  }
  /*  The system of ODEs is solved between [a, b] */
  /*  Number of days -> Also, number of steps used to solve the ODEs system: */
  /*  Time steps */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &t, 2, &dc_emlrtRTEI);
  i = t->size[0] * t->size[1];
  t->size[0] = 1;
  b_loop_ub = (int32_T)(r + 1.0);
  t->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, t, i, &dc_emlrtRTEI);
  t_data = t->data;
  for (i = 0; i < b_loop_ub; i++) {
    t_data[i] = 0.0;
  }
  /*  Matrix that stores the entire solution of the ODEs system for a */
  /*  particular value of uu:  */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &b_emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &w, 2, &ec_emlrtRTEI);
  i = w->size[0] * w->size[1];
  w->size[0] = 6;
  w->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, w, i, &ec_emlrtRTEI);
  w_data = w->data;
  loop_ub_tmp = 6 * (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub_tmp; i++) {
    w_data[i] = 0.0;
  }
  if (r == 0.0) {
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i].re = 0.0;
      beta1_data[i].im = 0.0;
      beta2_data[i].re = 0.0;
      beta2_data[i].im = 0.0;
    }
  } else {
    /*  Vector of Initial Conditions */
    h = T / r;
    /* the step size */
    if ((int32_T)(r + 1.0) < 1) {
      emlrtDynamicBoundsCheckR2012b(1, 1, (int32_T)(r + 1.0), &g_emlrtBCI,
                                    (emlrtConstCTX)sp);
    }
    for (i = 0; i < 6; i++) {
      w_data[i] = 0.0;
    }
    /* initial conditions */
    i = (int32_T)r;
    emlrtForLoopVectorCheckR2021a(1.0, 1.0, r, mxDOUBLE_CLASS, (int32_T)r,
                                  &emlrtRTEI, (emlrtConstCTX)sp);
    for (loop_ub_tmp = 0; loop_ub_tmp < i; loop_ub_tmp++) {
      real_T b_b[6];
      real_T b_w[6];
      real_T k1[6];
      real_T k2[6];
      real_T k3[6];
      b = ((loop_ub_tmp + 1 < 1) || (loop_ub_tmp + 1 > t->size[1]));
      if (b) {
        emlrtDynamicBoundsCheckR2012b(loop_ub_tmp + 1, 1, t->size[1],
                                      &h_emlrtBCI, (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &j_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      ODEsystem_2F_Gaussian_cg(&w_data[6 * loop_ub_tmp], 0.0, rf, d, k1);
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &k_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k1[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.0, rf, d, k2);
      for (i1 = 0; i1 < 6; i1++) {
        k2[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &l_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k2[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.0, rf, d, k3);
      for (i1 = 0; i1 < 6; i1++) {
        k3[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &m_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + k3[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.0, rf, d, b_b);
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &n_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      (int32_T)(r + 1.0), &o_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] = w_data[i1 + 6 * loop_ub_tmp] +
                 (((k1[i1] + 2.0 * k2[i1]) + 2.0 * k3[i1]) + h * b_b[i1]) / 6.0;
      }
      for (i1 = 0; i1 < 6; i1++) {
        w_data[i1 + 6 * (loop_ub_tmp + 1)] = k1[i1];
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > t->size[1])) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      t->size[1], &r_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      t_data[loop_ub_tmp + 1] = ((real_T)loop_ub_tmp + 1.0) * h;
      if (*emlrtBreakCheckR2012bFlagVar != 0) {
        emlrtBreakCheckR2012b((emlrtConstCTX)sp);
      }
    }
    /*  Time x uu */
    /*  each column of alpha, beta1, beta2, and beta3 contain  */
    /*  the solution of the ODE system for a specific value of uu */
    emxInit_creal_T(sp, &b_r, 1, &gc_emlrtRTEI);
    i = b_r->size[0];
    b_r->size[0] = (int32_T)(r + 1.0);
    emxEnsureCapacity_creal_T(sp, b_r, i, &fc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 1;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    emxInit_creal_T(sp, &r2, 1, &oc_emlrtRTEI);
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = (int32_T)(r + 1.0);
      emxEnsureCapacity_creal_T(sp, r2, i, &gc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &ed_emlrtRSI;
      binary_expand_op_14(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &yb_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      alpha_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &hc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 3;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &ic_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 2] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &fd_emlrtRSI;
      binary_expand_op_13(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &ac_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &jc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 5;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &kc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 4] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &gd_emlrtRSI;
      binary_expand_op_12(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emxFree_creal_T(sp, &b_r);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &bc_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta2_data[i] = r3[i];
    }
    emxFree_creal_T(sp, &r2);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &w->size[1], 1,
                                  &cc_emlrtECI, (emlrtCTX)sp);
  }
  emxFree_real_T(sp, &w);
  emxFree_real_T(sp, &t);
  b = (alpha->size[0] < 1);
  if (b) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &p_emlrtBCI, (emlrtConstCTX)sp);
  }
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.re = (alpha_data[alpha->size[0] - 1].re +
           beta1_data[alpha->size[0] - 1].re * sv1) +
          beta2_data[alpha->size[0] - 1].re * sv2;
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.im = (alpha_data[alpha->size[0] - 1].im +
           beta1_data[alpha->size[0] - 1].im * sv1) +
          beta2_data[alpha->size[0] - 1].im * sv2;
  emxFree_creal_T(sp, &beta2);
  emxFree_creal_T(sp, &beta1);
  emxFree_creal_T(sp, &alpha);
  if (rr.re == 0.0) {
    h = rr.im;
    rr.re = muDoubleScalarCos(h);
    rr.im = muDoubleScalarSin(h);
  } else if (rr.im == 0.0) {
    h = rr.re;
    rr.re = muDoubleScalarExp(h);
    rr.im = 0.0;
  } else if (muDoubleScalarIsInf(rr.im) && muDoubleScalarIsInf(rr.re) &&
             (rr.re < 0.0)) {
    rr.re = 0.0;
    rr.im = 0.0;
  } else {
    r = muDoubleScalarExp(rr.re / 2.0);
    h = rr.im;
    rr.re = r * (r * muDoubleScalarCos(h));
    rr.im = r * (r * muDoubleScalarSin(h));
  }
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
  return rr;
}

creal_T c_moment_gen_function_cg(const emlrtStack *sp, real_T sv1, real_T sv2,
                                 real_T T, real_T rf, real_T d)
{
  emlrtStack st;
  emxArray_creal_T *alpha;
  emxArray_creal_T *b_r;
  emxArray_creal_T *beta1;
  emxArray_creal_T *beta2;
  emxArray_creal_T *r2;
  emxArray_real_T *t;
  emxArray_real_T *w;
  creal_T rr;
  creal_T *alpha_data;
  creal_T *beta1_data;
  creal_T *beta2_data;
  creal_T *r1;
  creal_T *r3;
  real_T h;
  real_T r;
  real_T *t_data;
  real_T *w_data;
  int32_T b_loop_ub;
  int32_T i;
  int32_T i1;
  int32_T loop_ub;
  int32_T loop_ub_tmp;
  boolean_T b;
  st.prev = sp;
  st.tls = sp->tls;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  /*  moment_gen_function.m computes the Joint Conditional MGF (or the */
  /*  Characteristic Function -CF-) of log(S_T/S_t), V1_T, V2_T, and V3_T */
  /*  corresponding to the model "distr" at points uu, phi1, phi2, and phi3.  */
  /*  */
  /*  Note: */
  /*  CF_t(x,u,phi1,phi2,phi3,T) = */
  /*  E_t[exp(uu*y_T+phi1*V_{1,T}+phi2*V_{2,T}+phi3*V_{3,T})] = */
  /*  exp(alpha(x,T) + beta1(x,T)V_{1,t} + beta1(x,T)V_{2,t} +
   * beta3(x,T)V_{3,t}) */
  /* ==========================================================================
   */
  /*  OUTPUT:  */
  /*  */
  /*            M(uu,phi1,phi2,t,T) = E_t[exp(uu*y_T+ alpha + phi1*V_{1,T} +
   * phi2*V_{2,T})]  */
  /*  INPUT: */
  /*        distr:  model distribution */
  /*        x:      model parameters */
  /*        uu:     evalutation points for y_T (can be a vector) */
  /*        phi1:   evalutation points for V_1 (can be a vector) */
  /*        phi2:   evalutation points for V_2 (can be a vector) */
  /*        phi3:   evalutation points for V_3 (can be a vector) */
  /*        sv1:    initial value of the first factor */
  /*        sv2:    initial value of the second factor */
  /*        sv3:    initial value of the third factor */
  /*        T:      maturity (in years, convention 365 days per year) */
  /*        rf:     interest rate */
  /*        d:      dividend yield */
  /*  */
  /* --------------------------------------------------------------------------
   */
  /*  Nicola Fusari (2015): */
  /*  You can distribute and use the code. Please cite the author if you use */
  /*  the code in in your research. */
  /*   */
  /*  The code used to Solve the Systme of ODEs is from Alain Kapitho */
  /*  (University of Pretoria), available at: */
  /*  http://m2matlabdb.ma.tum.de/rk4_systems.m?MP_ID=413 */
  /*  */
  /* --------------------------------------------------------------------------
   */
  r = muDoubleScalarRound(T * 365.0);
  emxInit_creal_T(sp, &alpha, 1, &lc_emlrtRTEI);
  if (!(r + 1.0 >= 0.0)) {
    emlrtNonNegativeCheckR2012b(r + 1.0, &d_emlrtDCI, (emlrtConstCTX)sp);
  }
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  i = alpha->size[0];
  alpha->size[0] = (int32_T)(r + 1.0);
  emxEnsureCapacity_creal_T(sp, alpha, i, &ac_emlrtRTEI);
  alpha_data = alpha->data;
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  loop_ub = (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub; i++) {
    alpha_data[i].re = 0.0;
    alpha_data[i].im = 0.0;
  }
  emxInit_creal_T(sp, &beta1, 1, &mc_emlrtRTEI);
  loop_ub = alpha->size[0];
  i = beta1->size[0];
  beta1->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta1, i, &bc_emlrtRTEI);
  beta1_data = beta1->data;
  emxInit_creal_T(sp, &beta2, 1, &nc_emlrtRTEI);
  i = beta2->size[0];
  beta2->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta2, i, &cc_emlrtRTEI);
  beta2_data = beta2->data;
  for (i = 0; i < loop_ub; i++) {
    beta1_data[i] = alpha_data[i];
    beta2_data[i] = alpha_data[i];
  }
  /*  The system of ODEs is solved between [a, b] */
  /*  Number of days -> Also, number of steps used to solve the ODEs system: */
  /*  Time steps */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &t, 2, &dc_emlrtRTEI);
  i = t->size[0] * t->size[1];
  t->size[0] = 1;
  b_loop_ub = (int32_T)(r + 1.0);
  t->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, t, i, &dc_emlrtRTEI);
  t_data = t->data;
  for (i = 0; i < b_loop_ub; i++) {
    t_data[i] = 0.0;
  }
  /*  Matrix that stores the entire solution of the ODEs system for a */
  /*  particular value of uu:  */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &b_emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &w, 2, &ec_emlrtRTEI);
  i = w->size[0] * w->size[1];
  w->size[0] = 6;
  w->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, w, i, &ec_emlrtRTEI);
  w_data = w->data;
  loop_ub_tmp = 6 * (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub_tmp; i++) {
    w_data[i] = 0.0;
  }
  if (r == 0.0) {
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i].re = 0.0;
      beta1_data[i].im = 0.0;
      beta2_data[i].re = 0.0;
      beta2_data[i].im = 0.0;
    }
  } else {
    /*  Vector of Initial Conditions */
    h = T / r;
    /* the step size */
    if ((int32_T)(r + 1.0) < 1) {
      emlrtDynamicBoundsCheckR2012b(1, 1, (int32_T)(r + 1.0), &g_emlrtBCI,
                                    (emlrtConstCTX)sp);
    }
    for (i = 0; i < 6; i++) {
      w_data[i] = 0.0;
    }
    /* initial conditions */
    i = (int32_T)r;
    emlrtForLoopVectorCheckR2021a(1.0, 1.0, r, mxDOUBLE_CLASS, (int32_T)r,
                                  &emlrtRTEI, (emlrtConstCTX)sp);
    for (loop_ub_tmp = 0; loop_ub_tmp < i; loop_ub_tmp++) {
      real_T b_b[6];
      real_T b_w[6];
      real_T k1[6];
      real_T k2[6];
      real_T k3[6];
      b = ((loop_ub_tmp + 1 < 1) || (loop_ub_tmp + 1 > t->size[1]));
      if (b) {
        emlrtDynamicBoundsCheckR2012b(loop_ub_tmp + 1, 1, t->size[1],
                                      &h_emlrtBCI, (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &j_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      ODEsystem_2F_Gaussian_cg(&w_data[6 * loop_ub_tmp], 0.005, rf, d, k1);
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &k_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k1[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.005, rf, d, k2);
      for (i1 = 0; i1 < 6; i1++) {
        k2[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &l_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k2[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.005, rf, d, k3);
      for (i1 = 0; i1 < 6; i1++) {
        k3[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &m_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + k3[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, 0.005, rf, d, b_b);
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &n_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      (int32_T)(r + 1.0), &o_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] = w_data[i1 + 6 * loop_ub_tmp] +
                 (((k1[i1] + 2.0 * k2[i1]) + 2.0 * k3[i1]) + h * b_b[i1]) / 6.0;
      }
      for (i1 = 0; i1 < 6; i1++) {
        w_data[i1 + 6 * (loop_ub_tmp + 1)] = k1[i1];
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > t->size[1])) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      t->size[1], &r_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      t_data[loop_ub_tmp + 1] = ((real_T)loop_ub_tmp + 1.0) * h;
      if (*emlrtBreakCheckR2012bFlagVar != 0) {
        emlrtBreakCheckR2012b((emlrtConstCTX)sp);
      }
    }
    /*  Time x uu */
    /*  each column of alpha, beta1, beta2, and beta3 contain  */
    /*  the solution of the ODE system for a specific value of uu */
    emxInit_creal_T(sp, &b_r, 1, &gc_emlrtRTEI);
    i = b_r->size[0];
    b_r->size[0] = (int32_T)(r + 1.0);
    emxEnsureCapacity_creal_T(sp, b_r, i, &fc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 1;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    emxInit_creal_T(sp, &r2, 1, &oc_emlrtRTEI);
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = (int32_T)(r + 1.0);
      emxEnsureCapacity_creal_T(sp, r2, i, &gc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &ed_emlrtRSI;
      binary_expand_op_14(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &yb_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      alpha_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &hc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 3;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &ic_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 2] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &fd_emlrtRSI;
      binary_expand_op_13(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &ac_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &jc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 5;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &kc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 4] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &gd_emlrtRSI;
      binary_expand_op_12(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emxFree_creal_T(sp, &b_r);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &bc_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta2_data[i] = r3[i];
    }
    emxFree_creal_T(sp, &r2);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &w->size[1], 1,
                                  &cc_emlrtECI, (emlrtCTX)sp);
  }
  emxFree_real_T(sp, &w);
  emxFree_real_T(sp, &t);
  b = (alpha->size[0] < 1);
  if (b) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &p_emlrtBCI, (emlrtConstCTX)sp);
  }
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.re = (alpha_data[alpha->size[0] - 1].re +
           beta1_data[alpha->size[0] - 1].re * sv1) +
          beta2_data[alpha->size[0] - 1].re * sv2;
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.im = (alpha_data[alpha->size[0] - 1].im +
           beta1_data[alpha->size[0] - 1].im * sv1) +
          beta2_data[alpha->size[0] - 1].im * sv2;
  emxFree_creal_T(sp, &beta2);
  emxFree_creal_T(sp, &beta1);
  emxFree_creal_T(sp, &alpha);
  if (rr.re == 0.0) {
    h = rr.im;
    rr.re = muDoubleScalarCos(h);
    rr.im = muDoubleScalarSin(h);
  } else if (rr.im == 0.0) {
    h = rr.re;
    rr.re = muDoubleScalarExp(h);
    rr.im = 0.0;
  } else if (muDoubleScalarIsInf(rr.im) && muDoubleScalarIsInf(rr.re) &&
             (rr.re < 0.0)) {
    rr.re = 0.0;
    rr.im = 0.0;
  } else {
    r = muDoubleScalarExp(rr.re / 2.0);
    h = rr.im;
    rr.re = r * (r * muDoubleScalarCos(h));
    rr.im = r * (r * muDoubleScalarSin(h));
  }
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
  return rr;
}

void d_moment_gen_function_cg(const emlrtStack *sp, const emxArray_creal_T *uu,
                              const emxArray_real_T *phi1,
                              const emxArray_real_T *phi2, real_T sv1,
                              real_T sv2, real_T T, real_T rf, real_T d,
                              emxArray_creal_T *rr, emxArray_creal_T *alpha,
                              emxArray_creal_T *beta1, emxArray_creal_T *beta2,
                              emxArray_creal_T *beta3)
{
  emlrtStack b_st;
  emlrtStack c_st;
  emlrtStack d_st;
  emlrtStack st;
  emxArray_creal_T *b_r;
  emxArray_creal_T *r2;
  emxArray_creal_T *r3;
  emxArray_real_T *r1;
  emxArray_real_T *t;
  emxArray_real_T *w;
  const creal_T *uu_data;
  creal_T b_uu;
  creal_T *alpha_data;
  creal_T *beta1_data;
  creal_T *beta2_data;
  creal_T *beta3_data;
  creal_T *r4;
  creal_T *r5;
  real_T h;
  real_T r;
  real_T *r6;
  real_T *t_data;
  real_T *w_data;
  int32_T b_i;
  int32_T b_loop_ub;
  int32_T g;
  int32_T i;
  int32_T i1;
  int32_T loop_ub;
  int32_T nx;
  st.prev = sp;
  st.tls = sp->tls;
  b_st.prev = &st;
  b_st.tls = st.tls;
  c_st.prev = &b_st;
  c_st.tls = b_st.tls;
  d_st.prev = &c_st;
  d_st.tls = c_st.tls;
  uu_data = uu->data;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  /*  moment_gen_function.m computes the Joint Conditional MGF (or the */
  /*  Characteristic Function -CF-) of log(S_T/S_t), V1_T, V2_T, and V3_T */
  /*  corresponding to the model "distr" at points uu, phi1, phi2, and phi3.  */
  /*  */
  /*  Note: */
  /*  CF_t(x,u,phi1,phi2,phi3,T) = */
  /*  E_t[exp(uu*y_T+phi1*V_{1,T}+phi2*V_{2,T}+phi3*V_{3,T})] = */
  /*  exp(alpha(x,T) + beta1(x,T)V_{1,t} + beta1(x,T)V_{2,t} +
   * beta3(x,T)V_{3,t}) */
  /* ==========================================================================
   */
  /*  OUTPUT:  */
  /*  */
  /*            M(uu,phi1,phi2,t,T) = E_t[exp(uu*y_T+ alpha + phi1*V_{1,T} +
   * phi2*V_{2,T})]  */
  /*  INPUT: */
  /*        distr:  model distribution */
  /*        x:      model parameters */
  /*        uu:     evalutation points for y_T (can be a vector) */
  /*        phi1:   evalutation points for V_1 (can be a vector) */
  /*        phi2:   evalutation points for V_2 (can be a vector) */
  /*        phi3:   evalutation points for V_3 (can be a vector) */
  /*        sv1:    initial value of the first factor */
  /*        sv2:    initial value of the second factor */
  /*        sv3:    initial value of the third factor */
  /*        T:      maturity (in years, convention 365 days per year) */
  /*        rf:     interest rate */
  /*        d:      dividend yield */
  /*  */
  /* --------------------------------------------------------------------------
   */
  /*  Nicola Fusari (2015): */
  /*  You can distribute and use the code. Please cite the author if you use */
  /*  the code in in your research. */
  /*   */
  /*  The code used to Solve the Systme of ODEs is from Alain Kapitho */
  /*  (University of Pretoria), available at: */
  /*  http://m2matlabdb.ma.tum.de/rk4_systems.m?MP_ID=413 */
  /*  */
  /* --------------------------------------------------------------------------
   */
  r = muDoubleScalarRound(T * 365.0);
  if (!(r + 1.0 >= 0.0)) {
    emlrtNonNegativeCheckR2012b(r + 1.0, &f_emlrtDCI, (emlrtConstCTX)sp);
  }
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &g_emlrtDCI, (emlrtConstCTX)sp);
  }
  loop_ub = (int32_T)(r + 1.0);
  i = alpha->size[0] * alpha->size[1];
  alpha->size[0] = (int32_T)(r + 1.0);
  b_loop_ub = uu->size[1];
  alpha->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, alpha, i, &ac_emlrtRTEI);
  alpha_data = alpha->data;
  nx = (int32_T)(r + 1.0) * uu->size[1];
  for (i = 0; i < nx; i++) {
    alpha_data[i].re = 0.0;
    alpha_data[i].im = 0.0;
  }
  i = beta1->size[0] * beta1->size[1];
  beta1->size[0] = (int32_T)(r + 1.0);
  beta1->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, beta1, i, &bc_emlrtRTEI);
  beta1_data = beta1->data;
  for (i = 0; i < nx; i++) {
    beta1_data[i] = alpha_data[i];
  }
  i = beta2->size[0] * beta2->size[1];
  beta2->size[0] = (int32_T)(r + 1.0);
  beta2->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, beta2, i, &cc_emlrtRTEI);
  beta2_data = beta2->data;
  for (i = 0; i < nx; i++) {
    beta2_data[i] = alpha_data[i];
  }
  i = beta3->size[0] * beta3->size[1];
  beta3->size[0] = (int32_T)(r + 1.0);
  beta3->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, beta3, i, &qc_emlrtRTEI);
  beta3_data = beta3->data;
  for (i = 0; i < nx; i++) {
    beta3_data[i] = alpha_data[i];
  }
  /*  The system of ODEs is solved between [a, b] */
  /*  Number of days -> Also, number of steps used to solve the ODEs system: */
  /*  Time steps */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &t, 2, &dc_emlrtRTEI);
  i = t->size[0] * t->size[1];
  t->size[0] = 1;
  t->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, t, i, &dc_emlrtRTEI);
  t_data = t->data;
  for (i = 0; i < loop_ub; i++) {
    t_data[i] = 0.0;
  }
  /*  Matrix that stores the entire solution of the ODEs system for a */
  /*  particular value of uu:  */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &b_emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &w, 2, &ec_emlrtRTEI);
  i = w->size[0] * w->size[1];
  w->size[0] = 6;
  w->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, w, i, &ec_emlrtRTEI);
  w_data = w->data;
  nx = 6 * (int32_T)(r + 1.0);
  for (i = 0; i < nx; i++) {
    w_data[i] = 0.0;
  }
  if (r == 0.0) {
    for (g = 0; g < b_loop_ub; g++) {
      if (g + 1 > phi1->size[1]) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, phi1->size[1], &s_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &x_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i = 0; i < loop_ub; i++) {
        beta1_data[i + beta1->size[0] * g].re = 0.0;
        beta1_data[i + beta1->size[0] * g].im = 0.0;
      }
      if (g + 1 > phi2->size[1]) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, phi2->size[1], &t_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &ab_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i = 0; i < loop_ub; i++) {
        beta2_data[i + beta2->size[0] * g].re = 0.0;
        beta2_data[i + beta2->size[0] * g].im = 0.0;
      }
      if (*emlrtBreakCheckR2012bFlagVar != 0) {
        emlrtBreakCheckR2012b((emlrtConstCTX)sp);
      }
    }
  } else {
    if (uu->size[1] - 1 >= 0) {
      h = T / r;
      i1 = (int32_T)r;
    }
    emxInit_creal_T(sp, &b_r, 1, &gc_emlrtRTEI);
    emxInit_real_T(sp, &r1, 1, &uc_emlrtRTEI);
    emxInit_creal_T(sp, &r2, 1, &oc_emlrtRTEI);
    for (g = 0; g < b_loop_ub; g++) {
      /*  Vector of Initial Conditions */
      if (g + 1 > phi1->size[1]) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, phi1->size[1], &u_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (g + 1 > phi2->size[1]) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, phi2->size[1], &v_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      /* the step size */
      if (t->size[1] < 1) {
        emlrtDynamicBoundsCheckR2012b(1, 1, t->size[1], &g_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      nx = w->size[1];
      if (w->size[1] < 1) {
        emlrtDynamicBoundsCheckR2012b(1, 1, w->size[1], &i_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i = 0; i < 6; i++) {
        w_data[i] = 0.0;
      }
      /* initial conditions */
      emlrtForLoopVectorCheckR2021a(1.0, 1.0, r, mxDOUBLE_CLASS, (int32_T)r,
                                    &emlrtRTEI, (emlrtConstCTX)sp);
      for (b_i = 0; b_i < i1; b_i++) {
        real_T b_b[6];
        real_T b_w[6];
        real_T k1[6];
        real_T k2[6];
        real_T k3[6];
        real_T b_uu_tmp;
        real_T uu_tmp;
        boolean_T b;
        b = ((b_i + 1 < 1) || (b_i + 1 > t->size[1]));
        if (b) {
          emlrtDynamicBoundsCheckR2012b(b_i + 1, 1, t->size[1], &h_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        if (((int32_T)((uint32_T)b_i + 1U) < 1) ||
            ((int32_T)((uint32_T)b_i + 1U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 1U), 1, nx,
                                        &j_emlrtBCI, (emlrtConstCTX)sp);
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &db_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        uu_tmp = uu_data[g].re;
        b_uu.re = uu_tmp;
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &db_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        b_uu_tmp = uu_data[g].im;
        b_uu.im = b_uu_tmp;
        st.site = &s_emlrtRSI;
        b_ODEsystem_2F_Gaussian_cg(&w_data[6 * b_i], b_uu, rf, d, k1);
        for (i = 0; i < 6; i++) {
          k1[i] *= h;
        }
        if (((int32_T)((uint32_T)b_i + 1U) < 1) ||
            ((int32_T)((uint32_T)b_i + 1U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 1U), 1, nx,
                                        &k_emlrtBCI, (emlrtConstCTX)sp);
        }
        for (i = 0; i < 6; i++) {
          b_w[i] = w_data[i + 6 * b_i] + 0.5 * k1[i];
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &eb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &eb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        st.site = &t_emlrtRSI;
        b_ODEsystem_2F_Gaussian_cg(b_w, b_uu, rf, d, k2);
        for (i = 0; i < 6; i++) {
          k2[i] *= h;
        }
        if (((int32_T)((uint32_T)b_i + 1U) < 1) ||
            ((int32_T)((uint32_T)b_i + 1U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 1U), 1, nx,
                                        &l_emlrtBCI, (emlrtConstCTX)sp);
        }
        for (i = 0; i < 6; i++) {
          b_w[i] = w_data[i + 6 * b_i] + 0.5 * k2[i];
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &fb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &fb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        st.site = &u_emlrtRSI;
        b_ODEsystem_2F_Gaussian_cg(b_w, b_uu, rf, d, k3);
        for (i = 0; i < 6; i++) {
          k3[i] *= h;
        }
        if (((int32_T)((uint32_T)b_i + 1U) < 1) ||
            ((int32_T)((uint32_T)b_i + 1U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 1U), 1, nx,
                                        &m_emlrtBCI, (emlrtConstCTX)sp);
        }
        for (i = 0; i < 6; i++) {
          b_w[i] = w_data[i + 6 * b_i] + k3[i];
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &gb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        if (g + 1 > b_loop_ub) {
          emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &gb_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        st.site = &v_emlrtRSI;
        b_ODEsystem_2F_Gaussian_cg(b_w, b_uu, rf, d, b_b);
        if (((int32_T)((uint32_T)b_i + 1U) < 1) ||
            ((int32_T)((uint32_T)b_i + 1U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 1U), 1, nx,
                                        &n_emlrtBCI, (emlrtConstCTX)sp);
        }
        if (((int32_T)((uint32_T)b_i + 2U) < 1) ||
            ((int32_T)((uint32_T)b_i + 2U) > nx)) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 2U), 1, nx,
                                        &o_emlrtBCI, (emlrtConstCTX)sp);
        }
        for (i = 0; i < 6; i++) {
          k1[i] = w_data[i + 6 * b_i] +
                  (((k1[i] + 2.0 * k2[i]) + 2.0 * k3[i]) + h * b_b[i]) / 6.0;
        }
        for (i = 0; i < 6; i++) {
          w_data[i + 6 * (b_i + 1)] = k1[i];
        }
        if (((int32_T)((uint32_T)b_i + 2U) < 1) ||
            ((int32_T)((uint32_T)b_i + 2U) > t->size[1])) {
          emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)b_i + 2U), 1,
                                        t->size[1], &r_emlrtBCI,
                                        (emlrtConstCTX)sp);
        }
        t_data[b_i + 1] = ((real_T)b_i + 1.0) * h;
        if (*emlrtBreakCheckR2012bFlagVar != 0) {
          emlrtBreakCheckR2012b((emlrtConstCTX)sp);
        }
      }
      /*  Time x uu */
      /*  each column of alpha, beta1, beta2, and beta3 contain  */
      /*  the solution of the ODE system for a specific value of uu */
      i = b_r->size[0];
      b_r->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, b_r, i, &fc_emlrtRTEI);
      r4 = b_r->data;
      for (i = 0; i < nx; i++) {
        b_i = 6 * i + 1;
        r4[i].re = w_data[b_i] * 0.0;
        r4[i].im = w_data[b_i];
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &w_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (b_r->size[0] == w->size[1]) {
        i = r2->size[0];
        r2->size[0] = w->size[1];
        emxEnsureCapacity_creal_T(sp, r2, i, &gc_emlrtRTEI);
        r5 = r2->data;
        for (i = 0; i < nx; i++) {
          r5[i].re = w_data[6 * i] + r4[i].re;
          r5[i].im = r4[i].im;
        }
      } else {
        st.site = &ed_emlrtRSI;
        binary_expand_op_14(&st, r2, w, b_r);
        r5 = r2->data;
      }
      emlrtSubAssignSizeCheckR2012b(&loop_ub, 1, &r2->size[0], 1, &yb_emlrtECI,
                                    (emlrtCTX)sp);
      for (i = 0; i < loop_ub; i++) {
        alpha_data[i + alpha->size[0] * g] = r5[i];
      }
      nx = w->size[1];
      i = b_r->size[0];
      b_r->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, b_r, i, &hc_emlrtRTEI);
      r4 = b_r->data;
      for (i = 0; i < nx; i++) {
        b_i = 6 * i + 3;
        r4[i].re = w_data[b_i] * 0.0;
        r4[i].im = w_data[b_i];
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &y_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (b_r->size[0] == w->size[1]) {
        i = r2->size[0];
        r2->size[0] = w->size[1];
        emxEnsureCapacity_creal_T(sp, r2, i, &ic_emlrtRTEI);
        r5 = r2->data;
        for (i = 0; i < nx; i++) {
          r5[i].re = w_data[6 * i + 2] + r4[i].re;
          r5[i].im = r4[i].im;
        }
      } else {
        st.site = &fd_emlrtRSI;
        binary_expand_op_13(&st, r2, w, b_r);
        r5 = r2->data;
      }
      emlrtSubAssignSizeCheckR2012b(&loop_ub, 1, &r2->size[0], 1, &ac_emlrtECI,
                                    (emlrtCTX)sp);
      for (i = 0; i < loop_ub; i++) {
        beta1_data[i + beta1->size[0] * g] = r5[i];
      }
      nx = w->size[1];
      i = b_r->size[0];
      b_r->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, b_r, i, &jc_emlrtRTEI);
      r4 = b_r->data;
      for (i = 0; i < nx; i++) {
        b_i = 6 * i + 5;
        r4[i].re = w_data[b_i] * 0.0;
        r4[i].im = w_data[b_i];
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &bb_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (b_r->size[0] == w->size[1]) {
        i = r2->size[0];
        r2->size[0] = w->size[1];
        emxEnsureCapacity_creal_T(sp, r2, i, &kc_emlrtRTEI);
        r5 = r2->data;
        for (i = 0; i < nx; i++) {
          r5[i].re = w_data[6 * i + 4] + r4[i].re;
          r5[i].im = r4[i].im;
        }
      } else {
        st.site = &gd_emlrtRSI;
        binary_expand_op_12(&st, r2, w, b_r);
        r5 = r2->data;
      }
      emlrtSubAssignSizeCheckR2012b(&loop_ub, 1, &r2->size[0], 1, &bc_emlrtECI,
                                    (emlrtCTX)sp);
      for (i = 0; i < loop_ub; i++) {
        beta2_data[i + beta2->size[0] * g] = r5[i];
      }
      nx = w->size[1];
      i = r1->size[0];
      r1->size[0] = w->size[1];
      emxEnsureCapacity_real_T(sp, r1, i, &uc_emlrtRTEI);
      r6 = r1->data;
      i = b_r->size[0];
      b_r->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, b_r, i, &vc_emlrtRTEI);
      r4 = b_r->data;
      for (i = 0; i < nx; i++) {
        r6[i] = 0.0 * w_data[6 * i];
        b_i = 6 * i + 1;
        r4[i].re = w_data[b_i] * 0.0;
        r4[i].im = w_data[b_i] * 0.0;
      }
      if (g + 1 > b_loop_ub) {
        emlrtDynamicBoundsCheckR2012b(g + 1, 1, b_loop_ub, &cb_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (r1->size[0] == b_r->size[0]) {
        i = b_r->size[0];
        b_r->size[0] = w->size[1];
        emxEnsureCapacity_creal_T(sp, b_r, i, &uc_emlrtRTEI);
        r4 = b_r->data;
        for (i = 0; i < nx; i++) {
          r4[i].re += r6[i];
        }
      } else {
        st.site = &dd_emlrtRSI;
        binary_expand_op_22(&st, b_r, r1);
        r4 = b_r->data;
      }
      emlrtSubAssignSizeCheckR2012b(&loop_ub, 1, &b_r->size[0], 1, &cc_emlrtECI,
                                    (emlrtCTX)sp);
      for (i = 0; i < loop_ub; i++) {
        beta3_data[i + beta3->size[0] * g] = r4[i];
      }
      if (*emlrtBreakCheckR2012bFlagVar != 0) {
        emlrtBreakCheckR2012b((emlrtConstCTX)sp);
      }
    }
    emxFree_creal_T(sp, &r2);
    emxFree_real_T(sp, &r1);
    emxFree_creal_T(sp, &b_r);
  }
  emxFree_real_T(sp, &w);
  emxFree_real_T(sp, &t);
  if ((int32_T)(r + 1.0) < 1) {
    emlrtDynamicBoundsCheckR2012b((int32_T)(r + 1.0), 1, (int32_T)(r + 1.0),
                                  &p_emlrtBCI, (emlrtConstCTX)sp);
  }
  i = rr->size[0] * rr->size[1];
  rr->size[0] = 1;
  rr->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, rr, i, &rc_emlrtRTEI);
  beta3_data = rr->data;
  for (i = 0; i < b_loop_ub; i++) {
    beta3_data[i].re =
        sv1 * beta1_data[(beta1->size[0] + beta1->size[0] * i) - 1].re;
    beta3_data[i].im =
        sv1 * beta1_data[(beta1->size[0] + beta1->size[0] * i) - 1].im;
  }
  if (alpha->size[1] == rr->size[1]) {
    i = rr->size[0] * rr->size[1];
    rr->size[0] = 1;
    rr->size[1] = uu->size[1];
    emxEnsureCapacity_creal_T(sp, rr, i, &sc_emlrtRTEI);
    beta3_data = rr->data;
    for (i = 0; i < b_loop_ub; i++) {
      beta3_data[i].re +=
          alpha_data[(alpha->size[0] + alpha->size[0] * i) - 1].re;
      beta3_data[i].im +=
          alpha_data[(alpha->size[0] + alpha->size[0] * i) - 1].im;
    }
  } else {
    st.site = &kc_emlrtRSI;
    binary_expand_op_21(&st, rr, alpha);
  }
  emxInit_creal_T(sp, &r3, 2, &sc_emlrtRTEI);
  i = r3->size[0] * r3->size[1];
  r3->size[0] = 1;
  r3->size[1] = uu->size[1];
  emxEnsureCapacity_creal_T(sp, r3, i, &tc_emlrtRTEI);
  r4 = r3->data;
  for (i = 0; i < b_loop_ub; i++) {
    r4[i].re = sv2 * beta2_data[(beta2->size[0] + beta2->size[0] * i) - 1].re;
    r4[i].im = sv2 * beta2_data[(beta2->size[0] + beta2->size[0] * i) - 1].im;
  }
  if ((rr->size[1] != uu->size[1]) &&
      ((rr->size[1] != 1) && (uu->size[1] != 1))) {
    emlrtDimSizeImpxCheckR2021b(rr->size[1], uu->size[1], &dc_emlrtECI,
                                (emlrtConstCTX)sp);
  }
  st.site = &kc_emlrtRSI;
  if (rr->size[1] == r3->size[1]) {
    loop_ub = rr->size[1] - 1;
    i = rr->size[0] * rr->size[1];
    rr->size[0] = 1;
    emxEnsureCapacity_creal_T(&st, rr, i, &sc_emlrtRTEI);
    beta3_data = rr->data;
    for (i = 0; i <= loop_ub; i++) {
      beta3_data[i].re += r4[i].re;
      beta3_data[i].im += r4[i].im;
    }
  } else {
    b_st.site = &kc_emlrtRSI;
    b_plus(&b_st, rr, r3);
    beta3_data = rr->data;
  }
  emxFree_creal_T(&st, &r3);
  b_st.site = &lc_emlrtRSI;
  nx = rr->size[1];
  c_st.site = &gc_emlrtRSI;
  if (rr->size[1] > 2147483646) {
    d_st.site = &hc_emlrtRSI;
    check_forloop_overflow_error(&d_st);
  }
  for (b_i = 0; b_i < nx; b_i++) {
    if (beta3_data[b_i].re == 0.0) {
      h = beta3_data[b_i].im;
      beta3_data[b_i].re = muDoubleScalarCos(h);
      beta3_data[b_i].im = muDoubleScalarSin(h);
    } else if (beta3_data[b_i].im == 0.0) {
      beta3_data[b_i].re = muDoubleScalarExp(beta3_data[b_i].re);
      beta3_data[b_i].im = 0.0;
    } else if (muDoubleScalarIsInf(beta3_data[b_i].im) &&
               muDoubleScalarIsInf(beta3_data[b_i].re) &&
               (beta3_data[b_i].re < 0.0)) {
      beta3_data[b_i].re = 0.0;
      beta3_data[b_i].im = 0.0;
    } else {
      r = muDoubleScalarExp(beta3_data[b_i].re / 2.0);
      h = beta3_data[b_i].im;
      beta3_data[b_i].re = r * (r * muDoubleScalarCos(h));
      beta3_data[b_i].im = r * (r * muDoubleScalarSin(h));
    }
  }
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
}

creal_T moment_gen_function_cg(const emlrtStack *sp, real_T sv1, real_T sv2,
                               real_T T, real_T rf, real_T d)
{
  emlrtStack st;
  emxArray_creal_T *alpha;
  emxArray_creal_T *b_r;
  emxArray_creal_T *beta1;
  emxArray_creal_T *beta2;
  emxArray_creal_T *r2;
  emxArray_real_T *t;
  emxArray_real_T *w;
  creal_T rr;
  creal_T *alpha_data;
  creal_T *beta1_data;
  creal_T *beta2_data;
  creal_T *r1;
  creal_T *r3;
  real_T h;
  real_T r;
  real_T *t_data;
  real_T *w_data;
  int32_T b_loop_ub;
  int32_T i;
  int32_T i1;
  int32_T loop_ub;
  int32_T loop_ub_tmp;
  boolean_T b;
  st.prev = sp;
  st.tls = sp->tls;
  emlrtHeapReferenceStackEnterFcnR2012b((emlrtConstCTX)sp);
  /*  moment_gen_function.m computes the Joint Conditional MGF (or the */
  /*  Characteristic Function -CF-) of log(S_T/S_t), V1_T, V2_T, and V3_T */
  /*  corresponding to the model "distr" at points uu, phi1, phi2, and phi3.  */
  /*  */
  /*  Note: */
  /*  CF_t(x,u,phi1,phi2,phi3,T) = */
  /*  E_t[exp(uu*y_T+phi1*V_{1,T}+phi2*V_{2,T}+phi3*V_{3,T})] = */
  /*  exp(alpha(x,T) + beta1(x,T)V_{1,t} + beta1(x,T)V_{2,t} +
   * beta3(x,T)V_{3,t}) */
  /* ==========================================================================
   */
  /*  OUTPUT:  */
  /*  */
  /*            M(uu,phi1,phi2,t,T) = E_t[exp(uu*y_T+ alpha + phi1*V_{1,T} +
   * phi2*V_{2,T})]  */
  /*  INPUT: */
  /*        distr:  model distribution */
  /*        x:      model parameters */
  /*        uu:     evalutation points for y_T (can be a vector) */
  /*        phi1:   evalutation points for V_1 (can be a vector) */
  /*        phi2:   evalutation points for V_2 (can be a vector) */
  /*        phi3:   evalutation points for V_3 (can be a vector) */
  /*        sv1:    initial value of the first factor */
  /*        sv2:    initial value of the second factor */
  /*        sv3:    initial value of the third factor */
  /*        T:      maturity (in years, convention 365 days per year) */
  /*        rf:     interest rate */
  /*        d:      dividend yield */
  /*  */
  /* --------------------------------------------------------------------------
   */
  /*  Nicola Fusari (2015): */
  /*  You can distribute and use the code. Please cite the author if you use */
  /*  the code in in your research. */
  /*   */
  /*  The code used to Solve the Systme of ODEs is from Alain Kapitho */
  /*  (University of Pretoria), available at: */
  /*  http://m2matlabdb.ma.tum.de/rk4_systems.m?MP_ID=413 */
  /*  */
  /* --------------------------------------------------------------------------
   */
  r = muDoubleScalarRound(T * 365.0);
  emxInit_creal_T(sp, &alpha, 1, &lc_emlrtRTEI);
  if (!(r + 1.0 >= 0.0)) {
    emlrtNonNegativeCheckR2012b(r + 1.0, &d_emlrtDCI, (emlrtConstCTX)sp);
  }
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  i = alpha->size[0];
  alpha->size[0] = (int32_T)(r + 1.0);
  emxEnsureCapacity_creal_T(sp, alpha, i, &ac_emlrtRTEI);
  alpha_data = alpha->data;
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &c_emlrtDCI, (emlrtConstCTX)sp);
  }
  loop_ub = (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub; i++) {
    alpha_data[i].re = 0.0;
    alpha_data[i].im = 0.0;
  }
  emxInit_creal_T(sp, &beta1, 1, &mc_emlrtRTEI);
  loop_ub = alpha->size[0];
  i = beta1->size[0];
  beta1->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta1, i, &bc_emlrtRTEI);
  beta1_data = beta1->data;
  emxInit_creal_T(sp, &beta2, 1, &nc_emlrtRTEI);
  i = beta2->size[0];
  beta2->size[0] = alpha->size[0];
  emxEnsureCapacity_creal_T(sp, beta2, i, &cc_emlrtRTEI);
  beta2_data = beta2->data;
  for (i = 0; i < loop_ub; i++) {
    beta1_data[i] = alpha_data[i];
    beta2_data[i] = alpha_data[i];
  }
  /*  The system of ODEs is solved between [a, b] */
  /*  Number of days -> Also, number of steps used to solve the ODEs system: */
  /*  Time steps */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &t, 2, &dc_emlrtRTEI);
  i = t->size[0] * t->size[1];
  t->size[0] = 1;
  b_loop_ub = (int32_T)(r + 1.0);
  t->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, t, i, &dc_emlrtRTEI);
  t_data = t->data;
  for (i = 0; i < b_loop_ub; i++) {
    t_data[i] = 0.0;
  }
  /*  Matrix that stores the entire solution of the ODEs system for a */
  /*  particular value of uu:  */
  if (r + 1.0 != (int32_T)(r + 1.0)) {
    emlrtIntegerCheckR2012b(r + 1.0, &b_emlrtDCI, (emlrtConstCTX)sp);
  }
  emxInit_real_T(sp, &w, 2, &ec_emlrtRTEI);
  i = w->size[0] * w->size[1];
  w->size[0] = 6;
  w->size[1] = (int32_T)(r + 1.0);
  emxEnsureCapacity_real_T(sp, w, i, &ec_emlrtRTEI);
  w_data = w->data;
  loop_ub_tmp = 6 * (int32_T)(r + 1.0);
  for (i = 0; i < loop_ub_tmp; i++) {
    w_data[i] = 0.0;
  }
  if (r == 0.0) {
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i].re = 0.0;
      beta1_data[i].im = 0.0;
      beta2_data[i].re = 0.0;
      beta2_data[i].im = 0.0;
    }
  } else {
    /*  Vector of Initial Conditions */
    h = T / r;
    /* the step size */
    if ((int32_T)(r + 1.0) < 1) {
      emlrtDynamicBoundsCheckR2012b(1, 1, (int32_T)(r + 1.0), &g_emlrtBCI,
                                    (emlrtConstCTX)sp);
    }
    for (i = 0; i < 6; i++) {
      w_data[i] = 0.0;
    }
    /* initial conditions */
    i = (int32_T)r;
    emlrtForLoopVectorCheckR2021a(1.0, 1.0, r, mxDOUBLE_CLASS, (int32_T)r,
                                  &emlrtRTEI, (emlrtConstCTX)sp);
    for (loop_ub_tmp = 0; loop_ub_tmp < i; loop_ub_tmp++) {
      real_T b_b[6];
      real_T b_w[6];
      real_T k1[6];
      real_T k2[6];
      real_T k3[6];
      b = ((loop_ub_tmp + 1 < 1) || (loop_ub_tmp + 1 > t->size[1]));
      if (b) {
        emlrtDynamicBoundsCheckR2012b(loop_ub_tmp + 1, 1, t->size[1],
                                      &h_emlrtBCI, (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &j_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      ODEsystem_2F_Gaussian_cg(&w_data[6 * loop_ub_tmp], -0.005, rf, d, k1);
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &k_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k1[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, -0.005, rf, d, k2);
      for (i1 = 0; i1 < 6; i1++) {
        k2[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &l_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + 0.5 * k2[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, -0.005, rf, d, k3);
      for (i1 = 0; i1 < 6; i1++) {
        k3[i1] *= h;
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &m_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        b_w[i1] = w_data[i1 + 6 * loop_ub_tmp] + k3[i1];
      }
      ODEsystem_2F_Gaussian_cg(b_w, -0.005, rf, d, b_b);
      if (((int32_T)((uint32_T)loop_ub_tmp + 1U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 1U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 1U), 1,
                                      (int32_T)(r + 1.0), &n_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > (int32_T)(r + 1.0))) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      (int32_T)(r + 1.0), &o_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      for (i1 = 0; i1 < 6; i1++) {
        k1[i1] = w_data[i1 + 6 * loop_ub_tmp] +
                 (((k1[i1] + 2.0 * k2[i1]) + 2.0 * k3[i1]) + h * b_b[i1]) / 6.0;
      }
      for (i1 = 0; i1 < 6; i1++) {
        w_data[i1 + 6 * (loop_ub_tmp + 1)] = k1[i1];
      }
      if (((int32_T)((uint32_T)loop_ub_tmp + 2U) < 1) ||
          ((int32_T)((uint32_T)loop_ub_tmp + 2U) > t->size[1])) {
        emlrtDynamicBoundsCheckR2012b((int32_T)((uint32_T)loop_ub_tmp + 2U), 1,
                                      t->size[1], &r_emlrtBCI,
                                      (emlrtConstCTX)sp);
      }
      t_data[loop_ub_tmp + 1] = ((real_T)loop_ub_tmp + 1.0) * h;
      if (*emlrtBreakCheckR2012bFlagVar != 0) {
        emlrtBreakCheckR2012b((emlrtConstCTX)sp);
      }
    }
    /*  Time x uu */
    /*  each column of alpha, beta1, beta2, and beta3 contain  */
    /*  the solution of the ODE system for a specific value of uu */
    emxInit_creal_T(sp, &b_r, 1, &gc_emlrtRTEI);
    i = b_r->size[0];
    b_r->size[0] = (int32_T)(r + 1.0);
    emxEnsureCapacity_creal_T(sp, b_r, i, &fc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 1;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    emxInit_creal_T(sp, &r2, 1, &oc_emlrtRTEI);
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = (int32_T)(r + 1.0);
      emxEnsureCapacity_creal_T(sp, r2, i, &gc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &ed_emlrtRSI;
      binary_expand_op_14(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &yb_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      alpha_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &hc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 3;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &ic_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 2] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &fd_emlrtRSI;
      binary_expand_op_13(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &ac_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta1_data[i] = r3[i];
    }
    b_loop_ub = w->size[1];
    i = b_r->size[0];
    b_r->size[0] = w->size[1];
    emxEnsureCapacity_creal_T(sp, b_r, i, &jc_emlrtRTEI);
    r1 = b_r->data;
    for (i = 0; i < b_loop_ub; i++) {
      i1 = 6 * i + 5;
      r1[i].re = w_data[i1] * 0.0;
      r1[i].im = w_data[i1];
    }
    if (b_r->size[0] == w->size[1]) {
      i = r2->size[0];
      r2->size[0] = w->size[1];
      emxEnsureCapacity_creal_T(sp, r2, i, &kc_emlrtRTEI);
      r3 = r2->data;
      for (i = 0; i < b_loop_ub; i++) {
        r3[i].re = w_data[6 * i + 4] + r1[i].re;
        r3[i].im = r1[i].im;
      }
    } else {
      st.site = &gd_emlrtRSI;
      binary_expand_op_12(&st, r2, w, b_r);
      r3 = r2->data;
    }
    emxFree_creal_T(sp, &b_r);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &r2->size[0], 1,
                                  &bc_emlrtECI, (emlrtCTX)sp);
    for (i = 0; i < loop_ub; i++) {
      beta2_data[i] = r3[i];
    }
    emxFree_creal_T(sp, &r2);
    emlrtSubAssignSizeCheckR2012b(&alpha->size[0], 1, &w->size[1], 1,
                                  &cc_emlrtECI, (emlrtCTX)sp);
  }
  emxFree_real_T(sp, &w);
  emxFree_real_T(sp, &t);
  b = (alpha->size[0] < 1);
  if (b) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &p_emlrtBCI, (emlrtConstCTX)sp);
  }
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.re = (alpha_data[alpha->size[0] - 1].re +
           beta1_data[alpha->size[0] - 1].re * sv1) +
          beta2_data[alpha->size[0] - 1].re * sv2;
  if (alpha->size[0] < 1) {
    emlrtDynamicBoundsCheckR2012b(alpha->size[0], 1, alpha->size[0],
                                  &q_emlrtBCI, (emlrtConstCTX)sp);
  }
  rr.im = (alpha_data[alpha->size[0] - 1].im +
           beta1_data[alpha->size[0] - 1].im * sv1) +
          beta2_data[alpha->size[0] - 1].im * sv2;
  emxFree_creal_T(sp, &beta2);
  emxFree_creal_T(sp, &beta1);
  emxFree_creal_T(sp, &alpha);
  if (rr.re == 0.0) {
    h = rr.im;
    rr.re = muDoubleScalarCos(h);
    rr.im = muDoubleScalarSin(h);
  } else if (rr.im == 0.0) {
    h = rr.re;
    rr.re = muDoubleScalarExp(h);
    rr.im = 0.0;
  } else if (muDoubleScalarIsInf(rr.im) && muDoubleScalarIsInf(rr.re) &&
             (rr.re < 0.0)) {
    rr.re = 0.0;
    rr.im = 0.0;
  } else {
    r = muDoubleScalarExp(rr.re / 2.0);
    h = rr.im;
    rr.re = r * (r * muDoubleScalarCos(h));
    rr.im = r * (r * muDoubleScalarSin(h));
  }
  emlrtHeapReferenceStackLeaveFcnR2012b((emlrtConstCTX)sp);
  return rr;
}

/* End of code generation (moment_gen_function_cg.c) */
