Commit 43f1d3be authored by Lilli Frison's avatar Lilli Frison

update readme, add source files

parent fcd72d31
CONTENTS
Implementation of gPC (uniform, normal and truncated normal distribution) for
solving optimal control problems with uncertain parameters and initial conditions.
This repository only contains the gPC extension files that must be used on top
of a numerical optimal control solver (at this time only MUSCOD-II possible).
In the future, we hope make an open source substitute publicly available.
src-pce
Contains necessary extension files for gPC projection.
REQUIREMENTS
A copy of the source tree of the MUSCOD-II optimal control problem solver is
required to use the uocp extension.
As MUSCOD-II is not in the public domain, we are unfortunately bound by a
closed-source license at this time and its source code cannot be provided.
In the future, we hope make an open source substitute publicly available that
is currently under active development, and we hope to use ICMS (International
Conference on Mathematical Software) as one possible outlet.
examples
A collection of OCPs with uncertain parameters and initial conditions. NISP (non-intrusive spectral projection) is implemented on top of problem equations.
- van der Pol problem (Bock/Plitt, 1984)
- container crane (Goh/Teo, 1988)
- fermentation process
/*
*
* MUSCOD-II/MC2_TEST/SRC/ccrane.c
* (c) Daniel B. Leineweber, 1995
*
* container crane (Goh/Teo, 1988)
*
* $Id: ccrane.cpp 725 2011-08-10 23:04:16Z ckirches $
*
*/
#include <cmath>
#include "def_usrmod.hpp"
#define NMOS 1
#define NP 0
#define NRC 0
#define NRCE 0
#define NXD 7
#define NXA 0
#define NU 2
#define NPR 0
static void lfcn(double *t, double *xd, double *xa, double *u,
double *p, double *lval, double *rwh, long *iwh, InfoPtr *info)
{
*lval = 0.5*(xd[2]*xd[2] + xd[5]*xd[5]);
// *lval = 0.5*(xd[2]*xd[2] + xd[5]*xd[5]) + 0.01*u[0]*u[0] + 0.1*u[1]*u[1];
}
static void mfcn(double *ts, double *sd, double *sa, double *p, double *pr, double *mval, long *dpnd, InfoPtr *info)
{
if (*dpnd) { *dpnd = MFCN_DPND(0, *sd, 0, 0, 0); return; }
*mval = sd[NXD-1];
}
static void ffcn(double *t, double *xd, double *xa, double *u,double *p, double *rhs, double *rwh, long *iwh, InfoPtr *info)
{
rhs[0] = xd[3];
rhs[1] = xd[4];
rhs[2] = xd[5];
rhs[3] = u[0] + 1.76*9.81*xd[2];
rhs[4] = u[1];
rhs[5] = -(u[0] + (1.0+1.76)*9.81*xd[2] + 2*xd[4]*xd[5])/xd[1];
rhs[NXD-1] = 0.5*(xd[2]*xd[2] + xd[5]*xd[5] + 1e-4*u[1]*u[1]);
}
static void rdfcnRe(double *ts, double *sd, double *sa, double *u, double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) { *dpnd = RFCN_DPND(0, *sd, 0, 0, 0, *pr); return; }
res[0] = sd[0] - 9.9;
res[1] = sd[1] - 13.9;
res[2] = sd[2] + 0.001;
res[3] = sd[3] - 2.4;
res[4] = sd[4] + 0.001;
res[5] = sd[5] + 0.001;
res[6] = -sd[0] + 11.0;
res[7] = -sd[1] + 14.1;
res[8] = -sd[2] + 0.001;
res[9] = -sd[3] + 2.6;
res[10] = -sd[4] + 0.001;
res[11] = -sd[5] + 0.001;
}
/* -------------------------------------------- */
extern "C" void def_model(void);
void def_model(void)
{
def_mdims(NMOS, NP, NRC, NRCE);
def_mstage(0,NXD, NXA, NU, mfcn, NULL, 0, 0, 0, NULL, ffcn, NULL,NULL, NULL);
def_mpc(0, "End Point", NPR, 12,0, rdfcnRe, NULL);
//def_mpc(0, "Inner Point", NPR, 2,0, rdfcnR, NULL);
}
/* old */
static void rdfcn_s(double *ts, double *sd, double *sa, double *u,
double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) {
*dpnd = RFCN_DPND(0, *sd, 0, 0, 0, 0);
return;
}
res[0] = sd[0];
res[1] = sd[1] - 22.0;
res[2] = sd[2];
res[3] = sd[3];
res[4] = sd[4] + 1.0;
res[5] = sd[5];
}
static void rdfcn_e(double *ts, double *sd, double *sa, double *u,double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) {
*dpnd = RFCN_DPND(0, *sd, 0, 0, 0, 0);
return;
}
res[0] = sd[0] - 10.0;
res[1] = sd[1] - 14.0;
res[2] = sd[2];
res[3] = sd[3] - 2.5;
res[4] = sd[4];
res[5] = sd[5];
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
/* model: Model is adapted for simulation study by IWR
* author: Dr. Hergen Schultze, GVC/S - B009, BASF AG
* date: 2006-02-17
*/
#include <math.h>
#include <stdio.h>
#include "def_usrmod.hpp"
#define NMOS 1
#define NXD 9
#define NXA 0
#define NU 3
#define NP 7
#define NPR 0
#define NRC 0
#define NRCE 0
#define NRD_E 1
#define NRDE_E 1
static void ffcn(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, InfoPtr *info)
{
/*
c * x(1) <--> Pr ... product concentration
c * x(2) <--> S1 ... (fermentable) substrate1 concentration
c * x(3) <--> S2 ... substrate2 concentration
c * x(4) <--> En ... (living) biomass concentration
c * x(5) <--> Vo ... (absolute) volume
c * x(6) <--> Go ... goodies (out of batch) concentration
c #P
c #S1
c #S2
c
c * u(1) <--> FS1 ... feed rate substrate1
c * u(2) <--> FS2 ... feed rate substrate2
c * u(3) <--> Zw ... early harvest rate
c
c * p(1) <--> mu_X10 ... growth rate biomass
c * p(2) <--> mu_P10 ... production rate
c * p(3) <--> Gamma_XS1 ... consumption of substrate1 by growth
c * p(4) <--> Gamma_PS1 ... consumption of substrate1 by production
c * p(5) <--> Gamma_PS2 ... consumption of substrate2 by production
c * p(6) <--> Gamma_XS2 ... consumption of substrate2 by growth
c * p(7) <--> Gamma_XG ... consumption of goodies by growth
*/
double pr, s1, s2, en, vo, go;
double mu_x10, mu_p10, gamma_xs1, gamma_ps1, gamma_ps2, gamma_xs2, gamma_xg;
double fs1, fs2, zw, s11, s22;
/* differential states in right hand side */
pr = xd[0];
s1 = xd[1];
s2 = xd[2];
en = xd[3];
vo = xd[4];
go = xd[5];
/* parameters */
mu_x10 = p[0] * 1e4;
mu_p10 = p[1] * 1e4; // 1e3;
gamma_xs1= p[2] * 1e4;
gamma_ps1= p[3] * 1e4; // 1e3;
gamma_ps2= p[4] * 1e4; // 1e3;
gamma_xs2= p[5] * 1e4;
gamma_xg = p[6] * 1e4;
/* process controls */
fs1= u[0] / 25.0;
fs2= u[1] / 25.0;
zw = u[2] / 25.0;
/* intern parameters */
s11 = 42.0 / 1e+2;
s22 = 33.3 / 1e+2;
rhs[0] = mu_p10 * en * s1 * s2 - pr * (fs1 + fs2) / vo; // product concentration
rhs[1] = - gamma_xs1 * en * s1 * s2 * go - gamma_ps1 * en * s1 * s2 + s11 * fs1 / vo - s1 * (fs1 + fs2) / vo; // substrate1 fermentable
rhs[2] = - gamma_ps2 * en * s1 * s2 - gamma_xs2 * en * s1 * s2 * go + s22 * fs2 / vo - s2 * (fs1 + fs2) / vo; // substrate2
rhs[3] = mu_x10 * en * s1 * s2 * go - en * (fs1 + fs2) / vo; // cell volume
rhs[4] = fs1 + fs2 - zw; // fermentation volume
rhs[5] = - gamma_xg * en * s1 * s2 * go - go * (fs1 + fs2) / vo; // goodies
rhs[6] = zw * pr + (fs1 + fs2 - zw) * pr + vo * (mu_p10 * en * s1 * s2 - pr * (fs1 + fs2)/vo); // product accumulated
rhs[7] = fs1 * s11; // substrate1 accumulated
rhs[8] = fs2 * s22; // substrate2 accumulated
// product accumulated:
// #P = int (z * p) + v * p
// dP = z * p + dv * p + v * dp
return;
}
static void iffcn(double *t, double *xd, double *xdd, double *xa, double *u,double *p, double *rhs, double *rwh, long *iwh, InfoPtr *info)
{
double pr, s1, s2, en, vo, go;
double mu_x10, mu_p10, gamma_xs1, gamma_ps1, gamma_ps2, gamma_xs2, gamma_xg;
double fs1, fs2, zw, s11, s22;
/* differential states in right hand side */
pr = xd[0];
s1 = xd[1];
s2 = xd[2];
en = xd[3];
vo = xd[4];
go = xd[5];
/* parameters */
mu_x10 = p[0] * 1e4;
mu_p10 = p[1] * 1e4; // 1e3;
gamma_xs1= p[2] * 1e4;
gamma_ps1= p[3] * 1e4; // 1e3;
gamma_ps2= p[4] * 1e4; // 1e3;
gamma_xs2= p[5] * 1e4;
gamma_xg = p[6] * 1e4;
/* process controls */
fs1= u[0] / 25.0;
fs2= u[1] / 25.0;
zw = u[2] / 25.0;
/* intern parameters */
s11 = 42.0 / 1e+2;
s22 = 33.3 / 1e+2;
rhs[0] = mu_p10 * en * s1 * s2 - pr * (fs1 + fs2) / vo - xdd[0]; // product concentration
rhs[1] = - gamma_xs1 * en * s1 * s2 * go - gamma_ps1 * en * s1 * s2 + s11 * fs1 / vo - s1 * (fs1 + fs2) / vo - xdd[1]; // substrate1 fermentable
rhs[2] = - gamma_ps2 * en * s1 * s2 - gamma_xs2 * en * s1 * s2 * go + s22 * fs2 / vo - s2 * (fs1 + fs2) / vo - xdd[2]; // substrate2
rhs[3] = mu_x10 * en * s1 * s2 * go - en * (fs1 + fs2) / vo - xdd[3]; // cell volume
rhs[4] = fs1 + fs2 - zw - xdd[4]; // fermentation volume
rhs[5] = - gamma_xg * en * s1 * s2 * go - go * (fs1 + fs2) / vo - xdd[5]; // goodies
rhs[6] = zw * pr + (fs1 + fs2 - zw) * pr + vo * (mu_p10 * en * s1 * s2 - pr * (fs1 + fs2)/vo) - xdd[6]; // product accumulated
rhs[7] = fs1 * s11 - xdd[7]; // substrate1 accumulated
rhs[8] = fs2 * s22 - xdd[8]; // substrate2 accumulated
//rhs[6] = vo * mu_p10 * en * s1 * s2 - xdd[6];
}
static void mfcn(double *ts, double *sd, double *sa,
double *p, double *pr, double *mval, long *dpnd, InfoPtr *info)
{
if (*dpnd) {*dpnd = MFCN_DPND(0, *sd, 0, 0, 0); return;}
/*
double p_sub1, p_sub2; // price
double a_pro, a_sub1, a_sub2; // amount
p_sub1 = 1.0;
p_sub2 = 0.5;
a_pro = sd[6];
a_sub1 = sd[7];
a_sub2 = sd[8];
// variable costs:
*mval = (p_sub1 * a_sub1 + p_sub2 * a_sub2) / a_pro;*/
*mval = (sd[7] + 0.5*sd[8])/sd[6];
}
static void rdfcnR(double *ts, double *sd, double *sa, double *u, double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) {*dpnd = RFCN_DPND(0, *sd, 0, 0, 0, 0); return;}
res[0] = 0.04 - sd[1];
}
/* main function */
extern "C" void def_model(void){
def_mdims(NMOS, NP, NRC, NRCE);
def_mstage(0, NXD, NXA, NU, mfcn, NULL, 0, 0, 0, NULL, ffcn, NULL, NULL, NULL );
//def_imstage( 0, NXD, NXA, NU, 0, 0, 0, 0, mfcn, 2, iffcn, NULL,NULL, NULL, NULL, NULL);
//def_mpc( 0, "I", NPR, 1, 0, rdfcnR, NULL );
//def_mpc( 0, "E", NPR, 1, 0, rdfcnR, NULL );
}
/* unsued */
static void rdfcn_i(double *ts, double *sd, double *sa, double *u, double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) {*dpnd = RFCN_DPND(0, *sd, 0, 0, 0, 0); return;}
res[0] = 0.1 - sd[0]; // lb critical
res[1] = 0.06 - sd[1]; // ub critical (0.04 orig)
res[2] = 0.03 - sd[2];
res[3] = 0.1 - sd[3]; // Cell Volume
res[4] = 0.45 - sd[4]; // deterministic anyway
res[5] = 0.1 - sd[5];
res[6] = 0.05 - sd[6];
res[7] = 0.2 - sd[7];
res[8] = 0.025 - sd[8];
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
/*
*
* MUSCOD-II/MC2_TEST/SRC/vdpol.c
* (c) Daniel B. Leineweber, 1995
*
* van der Pol problem (Bock/Plitt 1984)
*
* $Id: vdpol.cpp 725 2011-08-10 23:04:16Z ckirches $
*
*/
#include <cmath>
#include "def_usrmod.hpp"
#define NMOS 1
#define NP 0
#define NRC 0
#define NRCE 0
#define NXD 2
#define NXA 0
#define NU 1
#define NPR 0
#define NRD_E 1
#define NRDE_E 1
static void lfcn(double *t, double *xd, double *xa, double *u,
double *p, double *lval, double *rwh, long *iwh, InfoPtr *info)
{
*lval = 0.5*(xd[0]*xd[0] + xd[1]*xd[1] + u[0]*u[0]);
}
static void ffcn(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, InfoPtr *info)
{
rhs[0] = xd[1];
rhs[1] = -xd[0] + (1.0 - xd[0]*xd[0])*xd[1] + u[0];
}
const double EPS = 0.0;
static void rdfcn_e(double *ts, double *sd, double *sa, double *u,
double *p, double *pr, double *res, long *dpnd, InfoPtr *info)
{
if (*dpnd) {
*dpnd = RFCN_DPND(0, *sd, 0, 0, 0, 0);
return;
}
res[0] = sd[1] - sd[0] - 1.0;
// res[0] = sd[1] - sd[0] - 1.0 + EPS;
// res[1] = -(sd[1] - sd[0] - 1.0) + EPS;
}
extern "C" void def_model(void)
{
def_mdims(NMOS, NP, NRC, NRCE);
def_mstage(
0,
NXD, NXA, NU,
NULL, lfcn,
0, 0, 0, NULL, ffcn, NULL,
NULL, NULL
);
def_mpc(0, "End Point", NPR, 1,0, rdfcn_e, NULL);
}
*
*
* MUSCOD-II/MDAT/vdpol.dat
* (c) Daniel B. Leineweber, 1995
*
* van der Pol problem (Bock/Plitt, 1984)
*
* $Id: vdpol.dat 658 2009-06-24 14:16:16Z potschka $
*
*
* # of multiple shooting intervals on each model stage
nshoot
0: 20
* model stage duration start values, scale factors, and bounds
h
0: 5.0
h_sca
0: 5.0
h_min
0: 5.0
h_max
0: 5.0
* model stage duration fixed value flags
h_fix
0: 1
* differential state start values, scale factors, and bounds
sd(*,*)
0: 1.0
1: 0.0
sd_fix(0,0)
ALL: 1
sd_sca(*,*)
ALL: 1
sd_min(*,*)
0: -10.0
1: -10.0
sd_max(*,*)
0: 10.0
1: 10.0
* control parameterization types
u_type(*)
0: 0
* control start values, scale factors, and bounds
u(*,*)
0: 0.5
u_sca(*,*)
0: 1.0
u_min(*,*)
0: -2.0
u_max(*,*)
0: 2.0
* control slope start values, scale factors, and bounds
udot(*,*)
0: 0.0
udot_sca(*,*)
0: 1.0
udot_min(*,*)
0: -10.0
udot_max(*,*)
0: 10.0
* decoupled i.p.c. scale factors
rd_sca(0,S)
0: 1.0
1: 1.0
rd_sca(0,E)
0: 1.0
1: 1.0
* objective scale and expected range; # of values in history plot
of_sca
1.0
of_min
0.0
of_max
10.0
nhist
30
**********************
* Choosing libraries *
**********************
libmodel
SRC/libvdpol0
libhessian
hess_update
libsolve
solve_tbox
libcond
cond_std
libtchk
tchk
libmssqp
mssqp_standard
libeval
eval_ind
libind
0: ind_rkfswt
libqps
qps_qpopt
libplot
plot_noplot
**********************************
* Setting algorithmic parameters *
**********************************
options_acc
1e-6
options_ftol
-1.0
options_itol
-1.0
options_rfac
0.0
options_levmar
0.0
options_qp_featol
1.0e-8
options_qp_relax
1.1
options_nhtopy
0
options_frstart
0
options_frmax
0
options_itmax
100
options_plevel_screen
0
options_plevel_file
1
options_plevel_matlab
0
options_bflag
-1
options_qp_itmax
10000
options_qp_expand
99999999
options_sflag
0
options_options_wflag
0
options_cflag
0
options_output_ps
0
options_output_gif
0
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
extern const double C;
extern const long T[];
extern double W[],X[];
void fnisp(int k, int NFUN, int M, double *W, double ** evalpoly, double **rhs, double *res);
void dfnisp(int k, int d, int NFUN, int M, double * W, double ** evalpoly, double **rhs,double *res);
/* non-normalized version */
void fnisp0(int k, int NFUN, int M, double *W, double ** evalpoly, double **rhs, double *res);
void dfnisp0(int k, int d, int NFUN, int M, double * W, double ** evalpoly, double **rhs,double *res);
double evalpce(int ND, int k, double *x, double *xi, int NORM);
double evaldpce(int ND, int k, double *x, double *xi, int NORM);
double evald2pce(int ND, int k, double *x, double *xi, int NORM);
double evalhermite(int d, int k, double *xi, int NORM);
/* ND=1 */
//void hermitequad(int k, int m, double x[], double w[]);
double evalpce_ND1(int k, double *x, double *xi);
double evaldpce_ND1(int k, double *x, double *xi);
double evald2pce_ND1(int k, double *x, double *xi);
double evalpcenorm_ND1(int k, double *x, double *xi);
double evaldpcenorm_ND1(int k, double *x, double *xi);
double evald2pcenorm_ND1(int k, double *x, double *xi);
double evalhermite_ND1(int k, double *xi);
double evalhermitenorm_ND1(int k, double *xi);
double evaltrunchermite_ND1(int k, double *xi);
double evaltruncpce_ND1(int k, double *x, double *xi);
/* ND=2 */
double evalpce_ND2(int k, double *x, double *xi);
double evalpcenorm_ND2(int k, double *x, double *xi);
double evaldpce_ND2(int k, double *x, double *xi);
double evaldpcenorm_ND2(int k, double *x, double *xi);
double evalhermite_ND2(int k, double *xi);
double evalhermitenorm_ND2(int k, double *xi);
double evaltrunchermite_ND2(int k, double *xi);
double evaltruncpce_ND2(int k, double *x, double *xi);
/* ND=3 */
//double evalpce_ND3(int k, double *x, double *xi);
double evalpcenorm_ND3(int k, double *x, double *xi);
//double evaldpce_ND3(int k, double *x, double *xi);
//double evalhermite_ND3(int k, double *xi);
double evalhermitenorm_ND3(int k, double *xi);
//double evaldpcenorm_ND3(int k, double *x, double *xi);
This diff is collapsed.
extern double W[],X[],T[];
void fnisp(int k, int NFUN, int M, double *W, double ** evalpoly, double **rhs, double *res);
void dfnisp(int k, int d, int NFUN, int M, double * W, double ** evalpoly, double **rhs,double *res);
double evalpce(int ND, int k, double *x, double *xi, int NORM);
double evaldpce(int ND, int k, double *x, double *xi, int NORM);
//double evald2pce(int ND, int k, double *x, double *xi, int NORM)
double evallegendre(int d, int k, double *xi, int NORM);
/******* ND=1 *******/
//double X[1] = {0.0};
//double X[5] = {-2.85697001387280558049, -1.35562617997426571037, -0.00000000000000003599, 1.35562617997426615446, 2.85697001387280558049};
//double X[7] = {-3.75043971772574113999, -2.36675941073454065844, -1.15440539473996817144, -0.00000000000000008676, 1.15440539473996817144, 2.36675941073454021435, 3.75043971772574113999};
//double X[12] = {-5.50090170446774529012, -4.27182584793228325992, -3.22370982877009648604, -2.25946445100079928991, -1.34037519715161645983, -0.44440300194413895341, 0.44440300194413884238, 1.34037519715161823619, 2.25946445100079840174, 3.22370982877009959466, 4.27182584793227970721, 5.50090170446774529012};
//double X[18] = {-7.13946484914648, -6.007745911359597, -5.05407268544274, -4.188020231629403, -3.374736535778091, -2.595833688911241, -1.839779921508647, -1.098395518091501, -0.3652457555076976, 0.3652457555076979, 1.098395518091501, 1.839779921508645, 2.59583368891124, 3.37473653577809, 4.188020231629404, 5.054072685442738, 6.007745911359601, 7.139464849146481};
double X[20] = {-7.61904854167975464918, -6.51059015701365240147, -5.57873880589320503276, -4.73458133404605518990, -3.94396735065731762759, -3.18901481655339003041, -2.45866361117236653655, -1.74524732081412636830, -1.04294534880275069355, -0.34696415708135580624, 0.34696415708135608380, 1.04294534880275113764, 1.74524732081412614626, 2.45866361117236831291, 3.18901481655339091859, 3.94396735065731629533, 4.73458133404605430172, 5.57873880589320414458, 6.51059015701365151330, 7.61904854167975997825};
/* N=30 */
//double X[28] = {-8.680837722732214, -7.825051744352809, -7.055396866960289, -6.339997686869598, -5.662381850082875, -5.012600596486519, -4.384020365898052, -3.771894423159235, -3.172634639420403, -2.583402100229275, -2.001858612956431, -1.426005658374114, -0.8540733517109721, -0.2844387607362086, 0.2844387607362088, 0.8540733517109731, 1.426005658374115, 2.001858612956432, 2.583402100229273, 3.172634639420402, 3.771894423159237, 4.384020365898054, 5.012600596486519, 5.662381850082874, 6.339997686869599, 7.055396866960293, 7.825051744352809, 8.680837722732212};
//double X[30] = {-9.706235997359524, -8.680837722732214, -7.825051744352809, -7.055396866960289, -6.339997686869598, -5.662381850082875, -5.012600596486519, -4.384020365898052, -3.771894423159235, -3.172634639420403, -2.583402100229275, -2.001858612956431, -1.426005658374114, -0.8540733517109721, -0.2844387607362086, 0.2844387607362088, 0.8540733517109731, 1.426005658374115, 2.001858612956432, 2.583402100229273, 3.172634639420402, 3.771894423159237, 4.384020365898054, 5.012600596486519, 5.662381850082874, 6.339997686869599, 7.055396866960293, 7.825051744352809, 8.680837722732212, 9.706235997359524};
/* N=40 */
//double X[34] = {-8.94950454385555538295, -8.27894062365947291937, -7.64616376454145907360, -7.04173840645383286585, -6.45942337758376794454, -5.89480567537201594064, -5.34460544572008444675, -4.80628719209387256228, -4.27782615636274776705, -3.75755977616898428906, -3.24408873299987066119, -2.73620834046543182083, -2.23285921863487191175, -1.73309059063172132831, -1.23603200479915797949, -0.74087072528593000964, -0.24683289602272429075, 0.24683289602272487362, 0.74087072528593178600, 1.23603200479916019994, 1.73309059063172243853, 2.23285921863487235584, 2.73620834046543048856, 3.24408873299987066119, 3.75755977616898562133, 4.27782615636274954340, 4.80628719209387256228, 5.34460544572008622310, 5.89480567537201860517, 6.45942337758376794454, 7.04173840645382842496, 7.64616376454145996178, 8.27894062365947469573, 8.94950454385555538295};
/*normalized */
const double C = 1.0; // scaling factor for weights
const long int T[32] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};
//double W[1] = {1.0};
//double W[5] = {0.01125741132772069969, 0.22207592200561251938, 0.53333333333333321491, 0.22207592200561246387, 0.01125741132772068581};
//double W[7] = {0.00054826885597221984, 0.03075712396758649436, 0.24012317860501264377, 0.45714285714285690654, 0.24012317860501275479, 0.03075712396758662967, 0.00054826885597221865};
//double W[12] = {0.00000014999271676372, 0.00004837184922590639, 0.00220338068753319748, 0.02911668791236416212, 0.14696704804533000654, 0.32166436151283034350, 0.32166436151283034350, 0.14696704804532972899, 0.02911668791236421069, 0.00220338068753319532, 0.00004837184922590633, 0.00000014999271676372};
//double W[18] = {0.00000000000441658877, 0.00000000590548847884, 0.00000102155239763698, 0.00005179896144116208, 0.00106548479629165068, 0.01051651775194133338, 0.05489663248022256387, 0.16068530389351259879, 0.27278323465428777617, 0.27278323465428755412, 0.16068530389351307064, 0.05489663248022253611, 0.01051651775194134032, 0.00106548479629165480, 0.00005179896144116197, 0.00000102155239763699, 0.00000000590548847884, 0.00000000000441658877};
double W[20] = {0.00000000000012578007, 0.00000000024820623623, 0.00000006127490259983, 0.00000440212109023081, 0.00012882627996192966, 0.00183010313108048677, 0.01399783744710104451, 0.06150637206397697315, 0.16173933398400003325, 0.26079306344955488495, 0.26079306344955477392, 0.16173933398399958916, 0.06150637206397704254, 0.01399783744710100808, 0.00183010313108048677, 0.00012882627996192923, 0.00000440212109023086, 0.00000006127490259983, 0.00000000024820623623, 0.00000000000012578007};
//double W[28] = {0.00000000000000001586, 0.00000000000001624080, 0.00000000000457342587, 0.00000000051784594672, 0.00000002882175154048, 0.00000089090888686211, 0.00001657998163067345, 0.00019651294398482469, 0.00154470733986608860, 0.00829574755772326082, 0.03111177018350123891, 0.08278683671562220292, 0.15804695320902073519, 0.21799997181557764780, 0.21799997181557731474, 0.15804695320902087396, 0.08278683671562188373, 0.03111177018350117993, 0.00829574755772333715, 0.00154470733986608903, 0.00019651294398482407, 0.00001657998163067352, 0.00000089090888686212, 0.00000002882175154048, 0.00000000051784594672, 0.00000000000457342587, 0.00000000000001624080, 0.00000000000000001586};
//double W[30] = {0.00000000000000000000, 0.00000000000000001586, 0.00000000000001624080, 0.00000000000457342587, 0.00000000051784594672, 0.00000002882175154048, 0.00000089090888686211, 0.00001657998163067345, 0.00019651294398482469, 0.00154470733986608860, 0.00829574755772326082, 0.03111177018350123891, 0.08278683671562220292, 0.15804695320902073519, 0.21799997181557764780, 0.21799997181557731474, 0.15804695320902087396, 0.08278683671562188373, 0.03111177018350117993, 0.00829574755772333715, 0.00154470733986608903, 0.00019651294398482407, 0.00001657998163067352, 0.00000089090888686212, 0.00000002882175154048, 0.00000000051784594672, 0.00000000000457342587, 0.00000000000001624080, 0.00000000000000001586, 0.00000000000000000000};
//N=40 points