0001
0002 #include <complex.h>
0003 #include <math.h>
0004 #include "mex.h"
0005
0006 struct polynomial_piecewise_0D {
0007 double *data;
0008 };
0009
0010
0011 struct polynomial_piecewise_1D {
0012 mxChar *form_f;
0013 double *breaks_f;
0014 double *coefs_f;
0015 double *pieces_f;
0016 double *order_f;
0017 double *dim_f;
0018 mxChar *orient_f;
0019 mxChar *form_dfdrho;
0020 double *breaks_dfdrho;
0021 double *coefs_dfdrho;
0022 double *pieces_dfdrho;
0023 double *order_dfdrho;
0024 double *dim_dfdrho;
0025 mxChar *orient_dfdrho;
0026 mxChar *form_d2fdrho2;
0027 double *breaks_d2fdrho2;
0028 double *coefs_d2fdrho2;
0029 double *pieces_d2fdrho2;
0030 double *order_d2fdrho2;
0031 double *dim_d2fdrho2;
0032 mxChar *orient_d2fdrho2;
0033 };
0034
0035 struct polynomial_piecewise_2D {
0036 double *n;
0037 mxChar *form_a0;
0038 double *breaks_a0;
0039 double *coefs_a0;
0040 double *pieces_a0;
0041 double *order_a0;
0042 double *dim_a0;
0043 mxChar *orient_a0;
0044 mxChar *form_da0drho;
0045 double *breaks_da0drho;
0046 double *coefs_da0drho;
0047 double *pieces_da0drho;
0048 double *order_da0drho;
0049 double *dim_da0drho;
0050 mxChar *orient_da0drho;
0051 mxChar *form_d2a0drho2;
0052 double *breaks_d2a0drho2;
0053 double *coefs_d2a0drho2;
0054 double *pieces_d2a0drho2;
0055 double *order_d2a0drho2;
0056 double *dim_d2a0drho2;
0057 mxChar *orient_d2a0drho2;
0058 mxChar *form_an;
0059 double *breaks_an;
0060 double *coefs_an;
0061 double *pieces_an;
0062 double *order_an;
0063 double *dim_an;
0064 mxChar *orient_an;
0065 mxChar *form_dandrho;
0066 double *breaks_dandrho;
0067 double *coefs_dandrho;
0068 double *pieces_dandrho;
0069 double *order_dandrho;
0070 double *dim_dandrho;
0071 mxChar *orient_dandrho;
0072 mxChar *form_d2andrho2;
0073 double *breaks_d2andrho2;
0074 double *coefs_d2andrho2;
0075 double *pieces_d2andrho2;
0076 double *order_d2andrho2;
0077 double *dim_d2andrho2;
0078 mxChar *orient_d2andrho2;
0079 mxChar *form_bn;
0080 double *breaks_bn;
0081 double *coefs_bn;
0082 double *pieces_bn;
0083 double *order_bn;
0084 double *dim_bn;
0085 mxChar *orient_bn;
0086 mxChar *form_dbndrho;
0087 double *breaks_dbndrho;
0088 double *coefs_dbndrho;
0089 double *pieces_dbndrho;
0090 double *order_dbndrho;
0091 double *dim_dbndrho;
0092 mxChar *orient_dbndrho;
0093 mxChar *form_d2bndrho2;
0094 double *breaks_d2bndrho2;
0095 double *coefs_d2bndrho2;
0096 double *pieces_d2bndrho2;
0097 double *order_d2bndrho2;
0098 double *dim_d2bndrho2;
0099 mxChar *orient_d2bndrho2;
0100 };
0101
0102 typedef struct polynomial_piecewise_0D pp_0D_format;
0103 typedef struct polynomial_piecewise_1D pp_1D_format;
0104 typedef struct polynomial_piecewise_2D pp_2D_format;
0105
0106
0107 void ppvald_0D(pp_0D_format in,double y[])
0108 {
0109 y[0] = in.data[0];
0110 }
0111
0112
0113 void ppvald_1D(pp_1D_format in,double drho_lim,double x,double y[],double dydY[],double d2ydY2[])
0114 {
0115 register int j,k;
0116 int i0,i1,i2,i3;
0117
0118 if (x>1.0) {x = 1.0 - drho_lim/2.0;}
0119
0120 for (j = 0;j < (int)in.pieces_f[0];j++) {
0121 if ((x>=in.breaks_f[j]) && (x<=in.breaks_f[j+1])) {
0122 for (k = 1;k <= (int)in.dim_f[0]; k++) {
0123 if ((int)in.order_f[0] == 2) {
0124 i0 = k-1+j*(int)in.dim_f[0] + 0*(int)in.pieces_f[0]*(int)in.dim_f[0];
0125 i1 = k-1+j*(int)in.dim_f[0] + 1*(int)in.pieces_f[0]*(int)in.dim_f[0];
0126 y[k-1] = in.coefs_f[i0]*(x - in.breaks_f[j]) + in.coefs_f[i1];
0127 } else if ((int)in.order_f[0] == 3) {
0128 i0 = k-1+j*(int)in.dim_f[0] + 0*(int)in.pieces_f[0]*(int)in.dim_f[0];
0129 i1 = k-1+j*(int)in.dim_f[0] + 1*(int)in.pieces_f[0]*(int)in.dim_f[0];
0130 i2 = k-1+j*(int)in.dim_f[0] + 2*(int)in.pieces_f[0]*(int)in.dim_f[0];
0131 y[k-1] = in.coefs_f[i0]*pow((x - in.breaks_f[j]),2) + in.coefs_f[i1]*(x - in.breaks_f[j]) + in.coefs_f[i2];
0132 } else if ((int)in.order_f[0] == 4) {
0133 i0 = k-1+j*(int)in.dim_f[0] + 0*(int)in.pieces_f[0]*(int)in.dim_f[0];
0134 i1 = k-1+j*(int)in.dim_f[0] + 1*(int)in.pieces_f[0]*(int)in.dim_f[0];
0135 i2 = k-1+j*(int)in.dim_f[0] + 2*(int)in.pieces_f[0]*(int)in.dim_f[0];
0136 i3 = k-1+j*(int)in.dim_f[0] + 3*(int)in.pieces_f[0]*(int)in.dim_f[0];
0137 y[k-1] = in.coefs_f[i0]*pow((x - in.breaks_f[j]),3) + in.coefs_f[i1]*pow((x - in.breaks_f[j]),2) + in.coefs_f[i2]*(x - in.breaks_f[j]) + in.coefs_f[i3];
0138
0139 }
0140 }
0141 }
0142 }
0143
0144 for (j = 0;j < (int)in.pieces_dfdrho[0];j++) {
0145 if ((x>=in.breaks_dfdrho[j]) && (x<=in.breaks_dfdrho[j+1])) {
0146 for (k = 1;k <= (int)in.dim_dfdrho[0]; k++) {
0147 if ((int)in.order_dfdrho[0] == 2) {
0148 i0 = k-1+j*(int)in.dim_dfdrho[0] + 0*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0149 i1 = k-1+j*(int)in.dim_dfdrho[0] + 1*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0150 dydY[k-1] = in.coefs_dfdrho[i0]*(x - in.breaks_dfdrho[j]) + in.coefs_dfdrho[i1];
0151 } else if ((int)in.order_dfdrho[0] == 3) {
0152 i0 = k-1+j*(int)in.dim_dfdrho[0] + 0*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0153 i1 = k-1+j*(int)in.dim_dfdrho[0] + 1*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0154 i2 = k-1+j*(int)in.dim_dfdrho[0] + 2*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0155 dydY[k-1] = in.coefs_dfdrho[i0]*pow((x - in.breaks_dfdrho[j]),2) + in.coefs_dfdrho[i1]*(x - in.breaks_dfdrho[j]) + in.coefs_dfdrho[i2];
0156
0157 } else if ((int)in.order_dfdrho[0] == 4) {
0158 i0 = k-1+j*(int)in.dim_dfdrho[0] + 0*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0159 i1 = k-1+j*(int)in.dim_dfdrho[0] + 1*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0160 i2 = k-1+j*(int)in.dim_dfdrho[0] + 2*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0161 i3 = k-1+j*(int)in.dim_dfdrho[0] + 3*(int)in.pieces_dfdrho[0]*(int)in.dim_dfdrho[0];
0162 dydY[k-1] = in.coefs_dfdrho[i0]*pow((x - in.breaks_dfdrho[j]),3) + in.coefs_dfdrho[i1]*pow((x - in.breaks_dfdrho[j]),2) + in.coefs_dfdrho[i2]*(x - in.breaks_dfdrho[j]) + in.coefs_dfdrho[i3];
0163 }
0164 }
0165 }
0166 }
0167
0168 for (j = 0;j < (int)in.pieces_d2fdrho2[0];j++) {
0169 if ((x>=in.breaks_d2fdrho2[j]) && (x<=in.breaks_d2fdrho2[j+1])) {
0170 for (k = 1;k <= (int)in.dim_d2fdrho2[0]; k++) {
0171 if ((int)in.order_d2fdrho2[0] == 2) {
0172 i0 = k-1+j*(int)in.dim_d2fdrho2[0] + 0*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0173 i1 = k-1+j*(int)in.dim_d2fdrho2[0] + 1*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0174 d2ydY2[k-1] = in.coefs_d2fdrho2[i0]*(x - in.breaks_d2fdrho2[j]) + in.coefs_d2fdrho2[i1];
0175 } else if ((int)in.order_d2fdrho2[0] == 3) {
0176 i0 = k-1+j*(int)in.dim_d2fdrho2[0] + 0*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0177 i1 = k-1+j*(int)in.dim_d2fdrho2[0] + 1*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0178 i2 = k-1+j*(int)in.dim_d2fdrho2[0] + 2*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0179 d2ydY2[k-1] = in.coefs_d2fdrho2[i0]*pow((x - in.breaks_d2fdrho2[j]),2) + in.coefs_d2fdrho2[i1]*(x - in.breaks_d2fdrho2[j]) + in.coefs_d2fdrho2[i2];
0180
0181 } else if ((int)in.order_d2fdrho2[0] == 4) {
0182 i0 = k-1+j*(int)in.dim_d2fdrho2[0] + 0*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0183 i1 = k-1+j*(int)in.dim_d2fdrho2[0] + 1*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0184 i2 = k-1+j*(int)in.dim_d2fdrho2[0] + 2*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0185 i3 = k-1+j*(int)in.dim_d2fdrho2[0] + 3*(int)in.pieces_d2fdrho2[0]*(int)in.dim_d2fdrho2[0];
0186 d2ydY2[k-1] = in.coefs_d2fdrho2[i0]*pow((x - in.breaks_d2fdrho2[j]),3) + in.coefs_d2fdrho2[i1]*pow((x - in.breaks_d2fdrho2[j]),2) + in.coefs_d2fdrho2[i2]*(x - in.breaks_d2fdrho2[j]) + in.coefs_d2fdrho2[i3];
0187 }
0188 }
0189 }
0190 }
0191
0192 }
0193
0194
0195
0196
0197 void ppvald_2D(pp_2D_format in,double drho_lim, double x, double t, double* y, double* dydY, double* dydt, double* d2ydt2, double* d2ydtdY, double* d2ydY2)
0198 {
0199 register int j,k;
0200 int i0,i1,i2,i3;
0201 double *a0,*an,*bn,*da0dY,*dandY,*dbndY,*d2a0dY2,*d2andY2,*d2bndY2;
0202
0203 if (x>1.0) {x = 1.0 - drho_lim/2.0;}
0204
0205 a0 = mxCalloc((int)in.dim_a0[0],sizeof(double));
0206 da0dY = mxCalloc((int)in.dim_da0drho[0],sizeof(double));
0207 d2a0dY2 = mxCalloc((int)in.dim_d2a0drho2[0],sizeof(double));
0208 an = mxCalloc((int)in.dim_an[0],sizeof(double));
0209 dandY = mxCalloc((int)in.dim_dandrho[0],sizeof(double));
0210 d2andY2 = mxCalloc((int)in.dim_d2andrho2[0],sizeof(double));
0211 bn = mxCalloc((int)in.dim_bn[0],sizeof(double));
0212 dbndY = mxCalloc((int)in.dim_dbndrho[0],sizeof(double));
0213 d2bndY2 = mxCalloc((int)in.dim_d2bndrho2[0],sizeof(double));
0214
0215
0216 for (j = 0;j < (int)in.pieces_a0[0];j++) {
0217 if ((x>=in.breaks_a0[j]) && (x<=in.breaks_a0[j+1])) {
0218 for (k = 1;k <= (int)in.dim_a0[0]; k++) {
0219 if ((int)in.order_a0[0] == 2) {
0220 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0221 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0222 a0[k-1] = in.coefs_a0[i0]*(x - in.breaks_a0[j]) + in.coefs_a0[i1];
0223 } else if ((int)in.order_a0[0] == 3) {
0224 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0225 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0226 i2 = k-1+j*(int)in.dim_a0[0] + 2*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0227 a0[k-1] = in.coefs_a0[i0]*pow((x - in.breaks_a0[j]),2) + in.coefs_a0[i1]*(x - in.breaks_a0[j]) + in.coefs_a0[i2];
0228 } else if ((int)in.order_a0[0] == 4) {
0229 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0230 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0231 i2 = k-1+j*(int)in.dim_a0[0] + 2*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0232 i3 = k-1+j*(int)in.dim_a0[0] + 3*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0233 a0[k-1] = in.coefs_a0[i0]*pow((x - in.breaks_a0[j]),3) + in.coefs_a0[i1]*pow((x - in.breaks_a0[j]),2) + in.coefs_a0[i2]*(x - in.breaks_a0[j]) + in.coefs_a0[i3];
0234 }
0235 }
0236 }
0237 }
0238
0239 for (j = 0;j < (int)in.pieces_da0drho[0];j++) {
0240 if ((x>=in.breaks_da0drho[j]) && (x<=in.breaks_da0drho[j+1])) {
0241 for (k = 1;k <= (int)in.dim_da0drho[0]; k++) {
0242 if ((int)in.order_da0drho[0] == 2) {
0243 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0244 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0245 da0dY[k-1] = in.coefs_da0drho[i0]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i1];
0246 } else if ((int)in.order_da0drho[0] == 3) {
0247 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0248 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0249 i2 = k-1+j*(int)in.dim_da0drho[0] + 2*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0250 da0dY[k-1] = in.coefs_da0drho[i0]*pow((x - in.breaks_da0drho[j]),2) + in.coefs_da0drho[i1]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i2];
0251 } else if ((int)in.order_da0drho[0] == 4) {
0252 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0253 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0254 i2 = k-1+j*(int)in.dim_da0drho[0] + 2*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0255 i3 = k-1+j*(int)in.dim_da0drho[0] + 3*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0256 da0dY[k-1] = in.coefs_da0drho[i0]*pow((x - in.breaks_da0drho[j]),3) + in.coefs_da0drho[i1]*pow((x - in.breaks_da0drho[j]),2) + in.coefs_da0drho[i2]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i3];
0257 }
0258 }
0259 }
0260 }
0261
0262
0263 for (j = 0;j < (int)in.pieces_d2a0drho2[0];j++) {
0264 if ((x>=in.breaks_d2a0drho2[j]) && (x<=in.breaks_d2a0drho2[j+1])) {
0265 for (k = 1;k <= (int)in.dim_d2a0drho2[0]; k++) {
0266 if ((int)in.order_d2a0drho2[0] == 2) {
0267 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0268 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0269 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i1];
0270 } else if ((int)in.order_d2a0drho2[0] == 3) {
0271 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0272 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0273 i2 = k-1+j*(int)in.dim_d2a0drho2[0] + 2*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0274 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*pow((x - in.breaks_d2a0drho2[j]),2) + in.coefs_d2a0drho2[i1]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i2];
0275 } else if ((int)in.order_d2a0drho2[0] == 4) {
0276 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0277 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0278 i2 = k-1+j*(int)in.dim_d2a0drho2[0] + 2*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0279 i3 = k-1+j*(int)in.dim_d2a0drho2[0] + 3*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0280 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*pow((x - in.breaks_d2a0drho2[j]),3) + in.coefs_d2a0drho2[i1]*pow((x - in.breaks_d2a0drho2[j]),2) + in.coefs_d2a0drho2[i2]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i3];
0281 }
0282 }
0283 }
0284 }
0285
0286
0287 for (j = 0;j < (int)in.pieces_an[0];j++) {
0288 if ((x>=in.breaks_an[j]) && (x<=in.breaks_an[j+1])) {
0289 for (k = 1;k <= (int)in.dim_an[0]; k++) {
0290 if ((int)in.order_an[0] == 2) {
0291 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0292 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0293 an[k-1] = in.coefs_an[i0]*(x - in.breaks_an[j]) + in.coefs_an[i1];
0294 } else if ((int)in.order_an[0] == 3) {
0295 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0296 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0297 i2 = k-1+j*(int)in.dim_an[0] + 2*(int)in.pieces_an[0]*(int)in.dim_an[0];
0298 an[k-1] = in.coefs_an[i0]*pow((x - in.breaks_an[j]),2) + in.coefs_an[i1]*(x - in.breaks_an[j]) + in.coefs_an[i2];
0299 } else if ((int)in.order_an[0] == 4) {
0300 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0301 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0302 i2 = k-1+j*(int)in.dim_an[0] + 2*(int)in.pieces_an[0]*(int)in.dim_an[0];
0303 i3 = k-1+j*(int)in.dim_an[0] + 3*(int)in.pieces_an[0]*(int)in.dim_an[0];
0304 an[k-1] = in.coefs_an[i0]*pow((x - in.breaks_an[j]),3) + in.coefs_an[i1]*pow((x - in.breaks_an[j]),2) + in.coefs_an[i2]*(x - in.breaks_an[j]) + in.coefs_an[i3];
0305 }
0306 }
0307 }
0308 }
0309
0310 for (j = 0;j < (int)in.pieces_dandrho[0];j++) {
0311 if ((x>=in.breaks_dandrho[j]) && (x<=in.breaks_dandrho[j+1])) {
0312 for (k = 1;k <= (int)in.dim_dandrho[0]; k++) {
0313 if ((int)in.order_dandrho[0] == 2) {
0314 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0315 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0316 dandY[k-1] = in.coefs_dandrho[i0]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i1];
0317 } else if ((int)in.order_dandrho[0] == 3) {
0318 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0319 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0320 i2 = k-1+j*(int)in.dim_dandrho[0] + 2*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0321 dandY[k-1] = in.coefs_dandrho[i0]*pow((x - in.breaks_dandrho[j]),2) + in.coefs_dandrho[i1]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i2];
0322 } else if ((int)in.order_dandrho[0] == 4) {
0323 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0324 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0325 i2 = k-1+j*(int)in.dim_dandrho[0] + 2*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0326 i3 = k-1+j*(int)in.dim_dandrho[0] + 3*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0327 dandY[k-1] = in.coefs_dandrho[i0]*pow((x - in.breaks_dandrho[j]),3) + in.coefs_dandrho[i1]*pow((x - in.breaks_dandrho[j]),2) + in.coefs_dandrho[i2]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i3];
0328 }
0329 }
0330 }
0331 }
0332
0333
0334 for (j = 0;j < (int)in.pieces_d2andrho2[0];j++) {
0335 if ((x>=in.breaks_d2andrho2[j]) && (x<=in.breaks_d2andrho2[j+1])) {
0336 for (k = 1;k <= (int)in.dim_d2andrho2[0]; k++) {
0337 if ((int)in.order_d2andrho2[0] == 2) {
0338 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0339 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0340 d2andY2[k-1] = in.coefs_d2andrho2[i0]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i1];
0341 } else if ((int)in.order_d2andrho2[0] == 3) {
0342 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0343 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0344 i2 = k-1+j*(int)in.dim_d2andrho2[0] + 2*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0345 d2andY2[k-1] = in.coefs_d2andrho2[i0]*pow((x - in.breaks_d2andrho2[j]),2) + in.coefs_d2andrho2[i1]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i2];
0346 } else if ((int)in.order_d2andrho2[0] == 4) {
0347 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0348 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0349 i2 = k-1+j*(int)in.dim_d2andrho2[0] + 2*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0350 i3 = k-1+j*(int)in.dim_d2andrho2[0] + 3*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0351 d2andY2[k-1] = in.coefs_d2andrho2[i0]*pow((x - in.breaks_d2andrho2[j]),3) + in.coefs_d2andrho2[i1]*pow((x - in.breaks_d2andrho2[j]),2) + in.coefs_d2andrho2[i2]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i3];
0352 }
0353 }
0354 }
0355 }
0356
0357
0358 for (j = 0;j < (int)in.pieces_bn[0];j++) {
0359 if ((x>=in.breaks_bn[j]) && (x<=in.breaks_bn[j+1])) {
0360 for (k = 1;k <= (int)in.dim_bn[0]; k++) {
0361 if ((int)in.order_bn[0] == 2) {
0362 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0363 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0364 bn[k-1] = in.coefs_bn[i0]*(x - in.breaks_bn[j]) + in.coefs_bn[i1];
0365 } else if ((int)in.order_bn[0] == 3) {
0366 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0367 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0368 i2 = k-1+j*(int)in.dim_bn[0] + 2*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0369 bn[k-1] = in.coefs_bn[i0]*pow((x - in.breaks_bn[j]),2) + in.coefs_bn[i1]*(x - in.breaks_bn[j]) + in.coefs_bn[i2];
0370 } else if ((int)in.order_bn[0] == 4) {
0371 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0372 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0373 i2 = k-1+j*(int)in.dim_bn[0] + 2*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0374 i3 = k-1+j*(int)in.dim_bn[0] + 3*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0375 bn[k-1] = in.coefs_bn[i0]*pow((x - in.breaks_bn[j]),3) + in.coefs_bn[i1]*pow((x - in.breaks_bn[j]),2) + in.coefs_bn[i2]*(x - in.breaks_bn[j]) + in.coefs_bn[i3];
0376 }
0377 }
0378 }
0379 }
0380
0381 for (j = 0;j < (int)in.pieces_dbndrho[0];j++) {
0382 if ((x>=in.breaks_dbndrho[j]) && (x<=in.breaks_dbndrho[j+1])) {
0383 for (k = 1;k <= (int)in.dim_dbndrho[0]; k++) {
0384 if ((int)in.order_dbndrho[0] == 2) {
0385 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0386 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0387 dbndY[k-1] = in.coefs_dbndrho[i0]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i1];
0388 } else if ((int)in.order_dbndrho[0] == 3) {
0389 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0390 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0391 i2 = k-1+j*(int)in.dim_dbndrho[0] + 2*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0392 dbndY[k-1] = in.coefs_dbndrho[i0]*pow((x - in.breaks_dbndrho[j]),2) + in.coefs_dbndrho[i1]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i2];
0393 } else if ((int)in.order_dbndrho[0] == 4) {
0394 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0395 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0396 i2 = k-1+j*(int)in.dim_dbndrho[0] + 2*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0397 i3 = k-1+j*(int)in.dim_dbndrho[0] + 3*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0398 dbndY[k-1] = in.coefs_dbndrho[i0]*pow((x - in.breaks_dbndrho[j]),3) + in.coefs_dbndrho[i1]*pow((x - in.breaks_dbndrho[j]),2) + in.coefs_dbndrho[i2]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i3];
0399 }
0400 }
0401 }
0402 }
0403
0404 for (j = 0;j < (int)in.pieces_d2bndrho2[0];j++) {
0405 if ((x>=in.breaks_d2bndrho2[j]) && (x<=in.breaks_d2bndrho2[j+1])) {
0406 for (k = 1;k <= (int)in.dim_d2bndrho2[0]; k++) {
0407 if ((int)in.order_d2bndrho2[0] == 2) {
0408 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0409 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0410 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i1];
0411 } else if ((int)in.order_d2bndrho2[0] == 3) {
0412 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0413 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0414 i2 = k-1+j*(int)in.dim_d2bndrho2[0] + 2*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0415 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*pow((x - in.breaks_d2bndrho2[j]),2) + in.coefs_d2bndrho2[i1]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i2];
0416 } else if ((int)in.order_d2bndrho2[0] == 4) {
0417 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0418 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0419 i2 = k-1+j*(int)in.dim_d2bndrho2[0] + 2*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0420 i3 = k-1+j*(int)in.dim_d2bndrho2[0] + 3*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0421 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*pow((x - in.breaks_d2bndrho2[j]),3) + in.coefs_d2bndrho2[i1]*pow((x - in.breaks_d2bndrho2[j]),2) + in.coefs_d2bndrho2[i2]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i3];
0422 }
0423 }
0424 }
0425 }
0426
0427
0428
0429 y[0] = 0.0;
0430 dydY[0] = 0.0;
0431 dydt[0] = 0.0;
0432 d2ydt2[0] = 0.0;
0433 d2ydtdY[0] = 0.0;
0434 d2ydY2[0] = 0.0;
0435
0436 for (k = 1;k <= (int)in.dim_dbndrho[0]; k++) {
0437 y[0] = y[0] + an[k-1]*cos(k*t) + bn[k-1]*sin(k*t);
0438 dydY[0] = dydY[0] + dandY[k-1]*cos(k*t) + dbndY[k-1]*sin(k*t);
0439 dydt[0] = dydt[0] - k*an[k-1]*sin(k*t) + k*bn[k-1]*cos(k*t);
0440 d2ydt2[0] = d2ydt2[0] - k*k*an[k-1]*cos(k*t) - k*k*bn[k-1]*sin(k*t);
0441 d2ydtdY[0] = d2ydtdY[0] - k*dandY[k-1]*sin(k*t) + k*dbndY[k-1]*cos(k*t);
0442 d2ydY2[0] = d2ydY2[0] + d2andY2[k-1]*cos(k*t) + d2bndY2[k-1]*sin(k*t);
0443 }
0444
0445 y[0] = y[0] + a0[0];
0446 dydY[0] = dydY[0] + da0dY[0];
0447 d2ydY2[0] = d2ydY2[0] + d2a0dY2[0];
0448
0449 mxFree(a0);
0450 mxFree(da0dY);
0451 mxFree(d2a0dY2);
0452 mxFree(an);
0453 mxFree(dandY);
0454 mxFree(d2andY2);
0455 mxFree(bn);
0456 mxFree(dbndY);
0457 mxFree(d2bndY2);
0458
0459 }
0460
0461
0462 void ppvali_2D(pp_2D_format in,double drho_lim, double x, double t0, double t, double* iy, double* idydY, double* id2ydY2)
0463 {
0464 register int j,k;
0465 int i0,i1,i2,i3;
0466 double *a0,*an,*bn,*da0dY,*dandY,*dbndY,*d2a0dY2,*d2andY2,*d2bndY2;
0467
0468 if (x>1.0) {x = 1.0 - drho_lim/2.0;}
0469
0470 a0 = mxCalloc((int)in.dim_a0[0],sizeof(double));
0471 da0dY = mxCalloc((int)in.dim_da0drho[0],sizeof(double));
0472 d2a0dY2 = mxCalloc((int)in.dim_d2a0drho2[0],sizeof(double));
0473 an = mxCalloc((int)in.dim_an[0],sizeof(double));
0474 dandY = mxCalloc((int)in.dim_dandrho[0],sizeof(double));
0475 d2andY2 = mxCalloc((int)in.dim_d2andrho2[0],sizeof(double));
0476 bn = mxCalloc((int)in.dim_bn[0],sizeof(double));
0477 dbndY = mxCalloc((int)in.dim_dbndrho[0],sizeof(double));
0478 d2bndY2 = mxCalloc((int)in.dim_d2bndrho2[0],sizeof(double));
0479
0480 for (j = 0;j < (int)in.pieces_a0[0];j++) {
0481 if ((x>=in.breaks_a0[j]) && (x<=in.breaks_a0[j+1])) {
0482 for (k = 1;k <= (int)in.dim_a0[0]; k++) {
0483 if ((int)in.order_a0[0] == 2) {
0484 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0485 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0486 a0[k-1] = in.coefs_a0[i0]*(x - in.breaks_a0[j]) + in.coefs_a0[i1];
0487 } else if ((int)in.order_a0[0] == 3) {
0488 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0489 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0490 i2 = k-1+j*(int)in.dim_a0[0] + 2*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0491 a0[k-1] = in.coefs_a0[i0]*pow((x - in.breaks_a0[j]),2) + in.coefs_a0[i1]*(x - in.breaks_a0[j]) + in.coefs_a0[i2];
0492 } else if ((int)in.order_a0[0] == 4) {
0493 i0 = k-1+j*(int)in.dim_a0[0] + 0*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0494 i1 = k-1+j*(int)in.dim_a0[0] + 1*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0495 i2 = k-1+j*(int)in.dim_a0[0] + 2*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0496 i3 = k-1+j*(int)in.dim_a0[0] + 3*(int)in.pieces_a0[0]*(int)in.dim_a0[0];
0497 a0[k-1] = in.coefs_a0[i0]*pow((x - in.breaks_a0[j]),3) + in.coefs_a0[i1]*pow((x - in.breaks_a0[j]),2) + in.coefs_a0[i2]*(x - in.breaks_a0[j]) + in.coefs_a0[i3];
0498 }
0499 }
0500 }
0501 }
0502
0503 for (j = 0;j < (int)in.pieces_da0drho[0];j++) {
0504 if ((x>=in.breaks_da0drho[j]) && (x<=in.breaks_da0drho[j+1])) {
0505 for (k = 1;k <= (int)in.dim_da0drho[0]; k++) {
0506 if ((int)in.order_da0drho[0] == 2) {
0507 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0508 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0509 da0dY[k-1] = in.coefs_da0drho[i0]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i1];
0510 } else if ((int)in.order_da0drho[0] == 3) {
0511 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0512 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0513 i2 = k-1+j*(int)in.dim_da0drho[0] + 2*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0514 da0dY[k-1] = in.coefs_da0drho[i0]*pow((x - in.breaks_da0drho[j]),2) + in.coefs_da0drho[i1]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i2];
0515 } else if ((int)in.order_da0drho[0] == 4) {
0516 i0 = k-1+j*(int)in.dim_da0drho[0] + 0*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0517 i1 = k-1+j*(int)in.dim_da0drho[0] + 1*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0518 i2 = k-1+j*(int)in.dim_da0drho[0] + 2*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0519 i3 = k-1+j*(int)in.dim_da0drho[0] + 3*(int)in.pieces_da0drho[0]*(int)in.dim_da0drho[0];
0520 da0dY[k-1] = in.coefs_da0drho[i0]*pow((x - in.breaks_da0drho[j]),3) + in.coefs_da0drho[i1]*pow((x - in.breaks_da0drho[j]),2) + in.coefs_da0drho[i2]*(x - in.breaks_da0drho[j]) + in.coefs_da0drho[i3];
0521 }
0522 }
0523 }
0524 }
0525
0526
0527
0528 for (j = 0;j < (int)in.pieces_d2a0drho2[0];j++) {
0529 if ((x>=in.breaks_d2a0drho2[j]) && (x<=in.breaks_d2a0drho2[j+1])) {
0530 for (k = 1;k <= (int)in.dim_d2a0drho2[0]; k++) {
0531 if ((int)in.order_d2a0drho2[0] == 2) {
0532 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0533 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0534 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i1];
0535 } else if ((int)in.order_d2a0drho2[0] == 3) {
0536 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0537 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0538 i2 = k-1+j*(int)in.dim_d2a0drho2[0] + 2*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0539 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*pow((x - in.breaks_d2a0drho2[j]),2) + in.coefs_d2a0drho2[i1]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i2];
0540 } else if ((int)in.order_d2a0drho2[0] == 4) {
0541 i0 = k-1+j*(int)in.dim_d2a0drho2[0] + 0*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0542 i1 = k-1+j*(int)in.dim_d2a0drho2[0] + 1*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0543 i2 = k-1+j*(int)in.dim_d2a0drho2[0] + 2*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0544 i3 = k-1+j*(int)in.dim_d2a0drho2[0] + 3*(int)in.pieces_d2a0drho2[0]*(int)in.dim_d2a0drho2[0];
0545 d2a0dY2[k-1] = in.coefs_d2a0drho2[i0]*pow((x - in.breaks_d2a0drho2[j]),3) + in.coefs_d2a0drho2[i1]*pow((x - in.breaks_d2a0drho2[j]),2) + in.coefs_d2a0drho2[i2]*(x - in.breaks_d2a0drho2[j]) + in.coefs_d2a0drho2[i3];
0546 }
0547 }
0548 }
0549 }
0550
0551
0552
0553 for (j = 0;j < (int)in.pieces_an[0];j++) {
0554 if ((x>=in.breaks_an[j]) && (x<=in.breaks_an[j+1])) {
0555 for (k = 1;k <= (int)in.dim_an[0]; k++) {
0556 if ((int)in.order_an[0] == 2) {
0557 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0558 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0559 an[k-1] = in.coefs_an[i0]*(x - in.breaks_an[j]) + in.coefs_an[i1];
0560 } else if ((int)in.order_an[0] == 3) {
0561 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0562 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0563 i2 = k-1+j*(int)in.dim_an[0] + 2*(int)in.pieces_an[0]*(int)in.dim_an[0];
0564 an[k-1] = in.coefs_an[i0]*pow((x - in.breaks_an[j]),2) + in.coefs_an[i1]*(x - in.breaks_an[j]) + in.coefs_an[i2];
0565 } else if ((int)in.order_an[0] == 4) {
0566 i0 = k-1+j*(int)in.dim_an[0] + 0*(int)in.pieces_an[0]*(int)in.dim_an[0];
0567 i1 = k-1+j*(int)in.dim_an[0] + 1*(int)in.pieces_an[0]*(int)in.dim_an[0];
0568 i2 = k-1+j*(int)in.dim_an[0] + 2*(int)in.pieces_an[0]*(int)in.dim_an[0];
0569 i3 = k-1+j*(int)in.dim_an[0] + 3*(int)in.pieces_an[0]*(int)in.dim_an[0];
0570 an[k-1] = in.coefs_an[i0]*pow((x - in.breaks_an[j]),3) + in.coefs_an[i1]*pow((x - in.breaks_an[j]),2) + in.coefs_an[i2]*(x - in.breaks_an[j]) + in.coefs_an[i3];
0571 }
0572 }
0573 }
0574 }
0575
0576 for (j = 0;j < (int)in.pieces_dandrho[0];j++) {
0577 if ((x>=in.breaks_dandrho[j]) && (x<=in.breaks_dandrho[j+1])) {
0578 for (k = 1;k <= (int)in.dim_dandrho[0]; k++) {
0579 if ((int)in.order_dandrho[0] == 2) {
0580 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0581 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0582 dandY[k-1] = in.coefs_dandrho[i0]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i1];
0583 } else if ((int)in.order_dandrho[0] == 3) {
0584 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0585 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0586 i2 = k-1+j*(int)in.dim_dandrho[0] + 2*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0587 dandY[k-1] = in.coefs_dandrho[i0]*pow((x - in.breaks_dandrho[j]),2) + in.coefs_dandrho[i1]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i2];
0588 } else if ((int)in.order_dandrho[0] == 4) {
0589 i0 = k-1+j*(int)in.dim_dandrho[0] + 0*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0590 i1 = k-1+j*(int)in.dim_dandrho[0] + 1*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0591 i2 = k-1+j*(int)in.dim_dandrho[0] + 2*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0592 i3 = k-1+j*(int)in.dim_dandrho[0] + 3*(int)in.pieces_dandrho[0]*(int)in.dim_dandrho[0];
0593 dandY[k-1] = in.coefs_dandrho[i0]*pow((x - in.breaks_dandrho[j]),3) + in.coefs_dandrho[i1]*pow((x - in.breaks_dandrho[j]),2) + in.coefs_dandrho[i2]*(x - in.breaks_dandrho[j]) + in.coefs_dandrho[i3];
0594 }
0595 }
0596 }
0597 }
0598
0599
0600 for (j = 0;j < (int)in.pieces_d2andrho2[0];j++) {
0601 if ((x>=in.breaks_d2andrho2[j]) && (x<=in.breaks_d2andrho2[j+1])) {
0602 for (k = 1;k <= (int)in.dim_d2andrho2[0]; k++) {
0603 if ((int)in.order_d2andrho2[0] == 2) {
0604 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0605 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0606 d2andY2[k-1] = in.coefs_d2andrho2[i0]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i1];
0607 } else if ((int)in.order_d2andrho2[0] == 3) {
0608 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0609 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0610 i2 = k-1+j*(int)in.dim_d2andrho2[0] + 2*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0611 d2andY2[k-1] = in.coefs_d2andrho2[i0]*pow((x - in.breaks_d2andrho2[j]),2) + in.coefs_d2andrho2[i1]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i2];
0612 } else if ((int)in.order_d2andrho2[0] == 4) {
0613 i0 = k-1+j*(int)in.dim_d2andrho2[0] + 0*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0614 i1 = k-1+j*(int)in.dim_d2andrho2[0] + 1*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0615 i2 = k-1+j*(int)in.dim_d2andrho2[0] + 2*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0616 i3 = k-1+j*(int)in.dim_d2andrho2[0] + 3*(int)in.pieces_d2andrho2[0]*(int)in.dim_d2andrho2[0];
0617 d2andY2[k-1] = in.coefs_d2andrho2[i0]*pow((x - in.breaks_d2andrho2[j]),3) + in.coefs_d2andrho2[i1]*pow((x - in.breaks_d2andrho2[j]),2) + in.coefs_d2andrho2[i2]*(x - in.breaks_d2andrho2[j]) + in.coefs_d2andrho2[i3];
0618 }
0619 }
0620 }
0621 }
0622
0623
0624
0625 for (j = 0;j < (int)in.pieces_bn[0];j++) {
0626 if ((x>=in.breaks_bn[j]) && (x<=in.breaks_bn[j+1])) {
0627 for (k = 1;k <= (int)in.dim_bn[0]; k++) {
0628 if ((int)in.order_bn[0] == 2) {
0629 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0630 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0631 bn[k-1] = in.coefs_bn[i0]*(x - in.breaks_bn[j]) + in.coefs_bn[i1];
0632 } else if ((int)in.order_bn[0] == 3) {
0633 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0634 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0635 i2 = k-1+j*(int)in.dim_bn[0] + 2*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0636 bn[k-1] = in.coefs_bn[i0]*pow((x - in.breaks_bn[j]),2) + in.coefs_bn[i1]*(x - in.breaks_bn[j]) + in.coefs_bn[i2];
0637 } else if ((int)in.order_bn[0] == 4) {
0638 i0 = k-1+j*(int)in.dim_bn[0] + 0*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0639 i1 = k-1+j*(int)in.dim_bn[0] + 1*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0640 i2 = k-1+j*(int)in.dim_bn[0] + 2*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0641 i3 = k-1+j*(int)in.dim_bn[0] + 3*(int)in.pieces_bn[0]*(int)in.dim_bn[0];
0642 bn[k-1] = in.coefs_bn[i0]*pow((x - in.breaks_bn[j]),3) + in.coefs_bn[i1]*pow((x - in.breaks_bn[j]),2) + in.coefs_bn[i2]*(x - in.breaks_bn[j]) + in.coefs_bn[i3];
0643 }
0644 }
0645 }
0646 }
0647
0648 for (j = 0;j < (int)in.pieces_dbndrho[0];j++) {
0649 if ((x>=in.breaks_dbndrho[j]) && (x<=in.breaks_dbndrho[j+1])) {
0650 for (k = 1;k <= (int)in.dim_dbndrho[0]; k++) {
0651 if ((int)in.order_dbndrho[0] == 2) {
0652 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0653 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0654 dbndY[k-1] = in.coefs_dbndrho[i0]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i1];
0655 } else if ((int)in.order_dbndrho[0] == 3) {
0656 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0657 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0658 i2 = k-1+j*(int)in.dim_dbndrho[0] + 2*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0659 dbndY[k-1] = in.coefs_dbndrho[i0]*pow((x - in.breaks_dbndrho[j]),2) + in.coefs_dbndrho[i1]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i2];
0660 } else if ((int)in.order_dbndrho[0] == 4) {
0661 i0 = k-1+j*(int)in.dim_dbndrho[0] + 0*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0662 i1 = k-1+j*(int)in.dim_dbndrho[0] + 1*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0663 i2 = k-1+j*(int)in.dim_dbndrho[0] + 2*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0664 i3 = k-1+j*(int)in.dim_dbndrho[0] + 3*(int)in.pieces_dbndrho[0]*(int)in.dim_dbndrho[0];
0665 dbndY[k-1] = in.coefs_dbndrho[i0]*pow((x - in.breaks_dbndrho[j]),3) + in.coefs_dbndrho[i1]*pow((x - in.breaks_dbndrho[j]),2) + in.coefs_dbndrho[i2]*(x - in.breaks_dbndrho[j]) + in.coefs_dbndrho[i3];
0666 }
0667 }
0668 }
0669 }
0670
0671
0672 for (j = 0;j < (int)in.pieces_d2bndrho2[0];j++) {
0673 if ((x>=in.breaks_d2bndrho2[j]) && (x<=in.breaks_d2bndrho2[j+1])) {
0674 for (k = 1;k <= (int)in.dim_d2bndrho2[0]; k++) {
0675 if ((int)in.order_d2bndrho2[0] == 2) {
0676 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0677 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0678 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i1];
0679 } else if ((int)in.order_d2bndrho2[0] == 3) {
0680 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0681 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0682 i2 = k-1+j*(int)in.dim_d2bndrho2[0] + 2*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0683 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*pow((x - in.breaks_d2bndrho2[j]),2) + in.coefs_d2bndrho2[i1]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i2];
0684 } else if ((int)in.order_d2bndrho2[0] == 4) {
0685 i0 = k-1+j*(int)in.dim_d2bndrho2[0] + 0*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0686 i1 = k-1+j*(int)in.dim_d2bndrho2[0] + 1*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0687 i2 = k-1+j*(int)in.dim_d2bndrho2[0] + 2*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0688 i3 = k-1+j*(int)in.dim_d2bndrho2[0] + 3*(int)in.pieces_d2bndrho2[0]*(int)in.dim_d2bndrho2[0];
0689 d2bndY2[k-1] = in.coefs_d2bndrho2[i0]*pow((x - in.breaks_d2bndrho2[j]),3) + in.coefs_d2bndrho2[i1]*pow((x - in.breaks_d2bndrho2[j]),2) + in.coefs_d2bndrho2[i2]*(x - in.breaks_d2bndrho2[j]) + in.coefs_d2bndrho2[i3];
0690 }
0691 }
0692 }
0693 }
0694
0695
0696
0697
0698 iy[0] = 0.0;
0699 idydY[0] = 0.0;
0700 id2ydY2[0] = 0.0;
0701
0702 for (k = 1;k <= (int)in.dim_dbndrho[0]; k++) {
0703 iy[0] = iy[0] + an[k-1]*(sin(k*t) - sin(k*t0))/k + bn[k-1]*(cos(k*t0) - cos(k*t))/k;
0704 idydY[0] = idydY[0] + dandY[k-1]*(sin(k*t) - sin(k*t0))/k + dbndY[k-1]*(cos(k*t0) - cos(k*t))/k;
0705 id2ydY2[0] = id2ydY2[0] + d2andY2[k-1]*(sin(k*t) - sin(k*t0))/k + d2bndY2[k-1]*(cos(k*t0) - cos(k*t))/k;
0706 }
0707
0708 iy[0] = iy[0] + a0[0]*(t - t0);
0709 idydY[0] = idydY[0] + da0dY[0]*(t - t0);
0710 id2ydY2[0] = id2ydY2[0] + d2a0dY2[0]*(t - t0);
0711
0712 mxFree(a0);
0713 mxFree(da0dY);
0714 mxFree(d2a0dY2);
0715 mxFree(an);
0716 mxFree(dandY);
0717 mxFree(d2andY2);
0718 mxFree(bn);
0719 mxFree(dbndY);
0720 mxFree(d2bndY2);
0721 }
0722