← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to routines/interf/intinterp.c

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2005-01-09 22:58:21 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20050109225821-473xr8vhgugxxx5j
Tags: 3.0-12
http://bugs.debian.org/255869
changed configure.in to build scilab's own malloc.o, closes: #255869

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/interf/intinterp.c
 
1
/*  
 
2
 *  PURPOSE
 
3
 *     interface code for some interpolation / approximation
 
4
 *     routines
 
5
 *
 
6
 *  AUTHOR
 
7
 *     Bruno Pincon
 
8
 */
 
9
 
 
10
#if WIN32
 
11
#include <string.h>
 
12
#endif
 
13
#include <math.h>
 
14
#include "../stack-c.h"
 
15
#include "interp.h"
 
16
 
 
17
 
 
18
 
 
19
#define min(a,b) ((a) < (b) ? (a) : (b))
 
20
#define max(a,b) ((a) < (b) ? (b) : (a))
 
21
 
 
22
 
 
23
#if WIN32
 
24
extern int C2F(dset)();
 
25
extern int C2F(bicubicinterpwithgrad)();
 
26
extern int C2F(bicubicinterpwithgradandhes)();
 
27
extern int C2F(db3ink)();
 
28
extern int C2F(driverdb3val)();
 
29
extern int C2F(driverdb3valwithgrad)();
 
30
#endif
 
31
 
 
32
 
 
33
 
 
34
enum {NOT_A_KNOT, NATURAL, CLAMPED, PERIODIC, FAST, FAST_PERIODIC, 
 
35
      MONOTONE, BY_ZERO, C0, LINEAR, BY_NAN, UNDEFINED};
 
36
 
 
37
typedef struct { char *str_type; int type; } TableType;
 
38
 
 
39
#define NB_SPLINE_TYPE 7
 
40
static TableType SplineTable[NB_SPLINE_TYPE] = { 
 
41
  { "not_a_knot"   , NOT_A_KNOT    },
 
42
  { "natural"      , NATURAL       },
 
43
  { "clamped"      , CLAMPED       },
 
44
  { "periodic"     , PERIODIC      },
 
45
  { "monotone"     , MONOTONE      },
 
46
  { "fast"         , FAST          },
 
47
  { "fast_periodic", FAST_PERIODIC }};
 
48
 
 
49
#define NB_OUTMODE 6
 
50
static TableType OutModeTable[NB_OUTMODE] = {
 
51
  { "C0"        , C0         },
 
52
  { "by_zero"   , BY_ZERO    },
 
53
  { "natural"   , NATURAL    },
 
54
  { "periodic"  , PERIODIC   },
 
55
  { "by_nan"    , BY_NAN     },
 
56
  { "linear"    , LINEAR     }};
 
57
 
 
58
static int good_order(double x[], int n)
 
59
{
 
60
  /*  test if x[i-1] < x[i] */
 
61
  static int first = 1;
 
62
  int i;
 
63
  static double inf;
 
64
 
 
65
  if ( first ) { inf = 1.0 / (double) (first - first) ; first = 0; }
 
66
 
 
67
  if (fabs(x[0]) == inf  ||  x[n-1] == inf)
 
68
    return ( 0 );
 
69
 
 
70
  for ( i = 1 ; i < n ; i++ ) 
 
71
    if ( ! (x[i-1] < x[i]) )   /* form of this test for detecting nan ... */ 
 
72
      return ( 0 );
 
73
 
 
74
  return ( 1 );
 
75
}
 
76
  
 
77
/*
 
78
 *  Routines sp�ciales : 
 
79
 *      r�cup�rer une hypermatrice r�elle
 
80
 *      r�cup�rer une string scilab sans convertion
 
81
 *      comparer une string scilab et une C string classique
 
82
 *  (elles ne modifient pas les variables correspondantes
 
83
 *   ce qui est fondamental pour que cette interface soit
 
84
 *   pass�e par r�f�rence).
 
85
 */
 
86
 
 
87
#define GetRhsScalarString(num,n,t) if (!get_rhs_scalar_string(num,n,t)) { return 0;}
 
88
#define GetRhsRealHMat(pos,H) if (!get_rhs_real_hmat(pos,H)) { return 0;}
 
89
 
 
90
typedef struct realhypermat {
 
91
  int dimsize;        /* number of dimensions of the hyper matrix */
 
92
  int size;           /* total number of elements : size = dims[0]x... x dims[dimsize-1] */
 
93
  int *dims;          /* number of elements in each dimension */
 
94
  double *R;          /* points to the elements  */
 
95
} RealHyperMat;
 
96
 
 
97
static int get_rhs_real_hmat(int num, RealHyperMat *H)
 
98
{
 
99
  int il, il1, il2, il3,/* it, */lw;
 
100
 
 
101
  lw = num + Top - Rhs;
 
102
  il = iadr(*lstk( lw )); 
 
103
  if ( *istk(il) < 0 )
 
104
    il = iadr(*istk(il+1));
 
105
          
 
106
  if ( *istk(il) != 17 )
 
107
    goto err;
 
108
  else if ( *istk(il+1) != 3 )  /* a hm mlist must have 3 fields */
 
109
    goto err;
 
110
 
 
111
  /*  get the pointers for the 3 fields */
 
112
  il1 = sadr(il+6);
 
113
  il2 = il1 + *istk(il+3) - 1;
 
114
  il3 = il1 + *istk(il+4) - 1;
 
115
  il1 = iadr(il1); il2 = iadr(il2); il3 = iadr(il3);
 
116
 
 
117
  /*  test if the first field is a matrix string with 3 components
 
118
   *  and that the first is "hm" (ie 17 22  in scilab char code)
 
119
   */
 
120
  if ( (*istk(il1) != 10)  |  ((*istk(il1+1))*(*istk(il1+2)) != 3)  )
 
121
    goto err;
 
122
  else if ( *istk(il1+5)-1 != 2 )  /* 1 str must have 2 chars */
 
123
    goto err;
 
124
  else if ( *istk(il1+8) != 17  || *istk(il1+9) != 22 )
 
125
    goto err;
 
126
 
 
127
  /*  get the 2d field */
 
128
  if ( (*istk(il2) != 8)  |  (*istk(il2+3) != 4) )
 
129
    goto err;
 
130
  H->dimsize = (*istk(il2+1))*(*istk(il2+2));
 
131
  H->dims = istk(il2+4);
 
132
 
 
133
  /*  get the 3d field */
 
134
  if ( !( *istk(il3) == 1 && *istk(il3+3)==0) )
 
135
    goto err;;
 
136
 
 
137
  H->size = (*istk(il3+1))*(*istk(il3+2));
 
138
  H->R = stk(sadr(il3+4));
 
139
 
 
140
  /* needed for Jpc stuff (putlhsvar) */
 
141
  Nbvars = Max(Nbvars,num);
 
142
  C2F(intersci).ntypes[num-1] = '$';
 
143
  C2F(intersci).iwhere[num-1] = *lstk(lw);
 
144
  C2F(intersci).lad[num-1] = 0;  /* a voir ? */
 
145
  return 1;
 
146
 
 
147
 err:
 
148
  Scierror(999,"Argument %d is not a real hypermatrix\r\n", num);
 
149
  return 0;
 
150
}
 
151
 
 
152
static int get_rhs_scalar_string(int num, int *length, int **tabchar)
 
153
{
 
154
  int il, lw;
 
155
 
 
156
  lw = num + Top - Rhs;
 
157
  il = iadr(*lstk( lw )); 
 
158
  if ( *istk(il) < 0 )
 
159
    il = iadr(*istk(il+1));
 
160
          
 
161
  if ( ! ( *istk(il) == 10  &&  (*istk(il+1))*(*istk(il+2)) == 1 ) ) 
 
162
    {
 
163
      /* we look for a scalar string */
 
164
      Scierror(999,"Argument %d is not a scalar string\r\n", num);
 
165
      return 0;
 
166
    }
 
167
  *length = *istk(il+5)-1;
 
168
  *tabchar = istk(il+6);
 
169
 
 
170
  Nbvars = Max(Nbvars,num);
 
171
  C2F(intersci).ntypes[num-1] = '$';
 
172
  C2F(intersci).iwhere[num-1] = *lstk(lw);
 
173
  C2F(intersci).lad[num-1] = 0;
 
174
  return 1;
 
175
}
 
176
 
 
177
static int equal_scistring_and_string(int length, int *scistr,  char *str)
 
178
{
 
179
  /* compare a scistring with a classic C string */
 
180
  int i, res;
 
181
  
 
182
  if ( strlen(str) != length )
 
183
    return 0;
 
184
 
 
185
  res = 1; i = 0;
 
186
  while (res && i < length)
 
187
    {
 
188
      res = (scistr[i] == (int)C2F(getfastcode)(str+i,1L));
 
189
      i++;
 
190
    }
 
191
  return (res);
 
192
}
 
193
 
 
194
static int get_type(TableType *Tab, int dim_table, int *scistr, int strlength)
 
195
{
 
196
  int i = 0, found = 0;
 
197
  while ( !found  &&  i < dim_table )
 
198
    {
 
199
      found =  equal_scistring_and_string(strlength, scistr, Tab[i].str_type);
 
200
      i++;
 
201
    } 
 
202
  if ( found )
 
203
    return ( Tab[i-1].type );
 
204
  else
 
205
    return ( UNDEFINED );
 
206
} 
 
207
 
 
208
 
 
209
static int intsplin(char * fname)
 
210
{
 
211
  int minrhs=2, maxrhs=4, minlhs=1, maxlhs=1;
 
212
 
 
213
  int mx, nx, lx, my, ny, ly, mc, nc, lc, n, spline_type;
 
214
  int *str_spline_type, ns;
 
215
  int /*i,*/ ld/*, flag*/;
 
216
  int mwk1, nwk1, lwk1, mwk2, nwk2, lwk2, mwk3, nwk3, lwk3, mwk4, nwk4, lwk4;
 
217
  double *x, *y, *d, *c;
 
218
 
 
219
  CheckRhs(minrhs,maxrhs);
 
220
  CheckLhs(minlhs,maxlhs);
 
221
 
 
222
  GetRhsVar(1,"d", &mx, &nx, &lx);
 
223
  GetRhsVar(2,"d", &my, &ny, &ly);
 
224
 
 
225
  if ( mx != my  ||  nx != ny  ||  mx != 1  &&  nx != 1 ) 
 
226
    { 
 
227
      Scierror(999,"%s: arg1 and arg2 must be 2 vectors with same size\r\n", fname);
 
228
      return 0;
 
229
    }
 
230
  n = mx*nx;    /* number of interpolation points */
 
231
  if ( n < 2 ) 
 
232
    { 
 
233
      Scierror(999,"%s: the number of interpolation points must be >= 2\r\n", fname);
 
234
      return 0;
 
235
    }
 
236
  x = stk(lx); y = stk(ly);
 
237
  if (! good_order(x, n))  /* verify strict increasing abscissae */
 
238
    {
 
239
      Scierror(999,"%s: elts of arg 1 not (strictly) increasing or +-inf detected\r\n", fname);
 
240
      return 0;
 
241
    }
 
242
 
 
243
  if ( Rhs >= 3 )   /* get the spline type */
 
244
    {
 
245
      GetRhsScalarString(3, &ns, &str_spline_type);
 
246
      spline_type =  get_type(SplineTable, NB_SPLINE_TYPE, str_spline_type, ns);
 
247
      if ( spline_type == UNDEFINED )
 
248
        {
 
249
          Scierror(999,"%s: unknown spline_type\n\r",fname);
 
250
          return 0;
 
251
        };
 
252
    }
 
253
  else
 
254
    spline_type = NOT_A_KNOT;
 
255
 
 
256
  if ( spline_type == CLAMPED ) /* get arg 4 which contains the end point slopes */
 
257
    {
 
258
      if ( Rhs != 4 )
 
259
        {
 
260
          Scierror(999,"%s: for a clamped spline you must give the endpoint slopes\n\r",fname);
 
261
          return 0;
 
262
        }
 
263
      GetRhsVar(4,"d", &mc, &nc, &lc);
 
264
      if ( mc*nc != 2 )
 
265
        {
 
266
          Scierror(999,"%s: bad dimension for the 4 arg (endpoint slopes)\n\r",fname);
 
267
          return 0;
 
268
        }
 
269
      c = stk(lc);
 
270
    }
 
271
  else if ( Rhs == 4 )
 
272
    {
 
273
      Scierror(999,"%s: 4 args are required only for a clamped spline\n\r",fname);
 
274
      return 0;
 
275
    }
 
276
    
 
277
  /*  verify y(1) = y(n) for periodic splines */
 
278
  if ( (spline_type == PERIODIC || spline_type == FAST_PERIODIC)  &&  y[0] != y[n-1] )
 
279
    {
 
280
      Scierror(999,"%s: for a periodic spline y(1) must be equal to y(n)\n\r",fname);
 
281
      return(0);
 
282
    };
 
283
 
 
284
  CreateVar(Rhs+1, "d", &mx,  &nx,   &ld); /* memory for d (only argument returned) */   
 
285
  d = stk(ld);
 
286
 
 
287
  switch(spline_type)
 
288
    {
 
289
    case(FAST) : case(FAST_PERIODIC) :
 
290
      nwk1 = 1;
 
291
      C2F(derivd) (x, y, d, &n, &nwk1, &spline_type);
 
292
      break;
 
293
 
 
294
    case(MONOTONE) :
 
295
      nwk1 = 1;
 
296
      C2F(dpchim) (&n, x, y, d, &nwk1);
 
297
      break;
 
298
 
 
299
    case(NOT_A_KNOT) : case(NATURAL) : case(CLAMPED) : case(PERIODIC) :
 
300
      /*  (the wk4 work array is used only in the periodic case) */
 
301
      mwk1 = n; nwk1 = 1; mwk2 = n-1; nwk2 = 1; mwk3 = n-1; nwk3 = 1; mwk4 = n-2; nwk4 = 1;
 
302
      CreateVar(Rhs+2, "d", &mwk1,  &nwk1,   &lwk1);
 
303
      CreateVar(Rhs+3, "d", &mwk2,  &nwk2,   &lwk2);
 
304
      CreateVar(Rhs+4, "d", &mwk3,  &nwk3,   &lwk3);
 
305
      lwk4 = lwk1;
 
306
      if (spline_type == CLAMPED) 
 
307
        { d[0] = c[0]; d[n-1] = c[1]; };
 
308
      if (spline_type == PERIODIC)
 
309
        CreateVar(Rhs+5, "d", &mwk4,  &nwk4,   &lwk4);
 
310
      C2F(splinecub) (x, y, d, &n, &spline_type, stk(lwk1), stk(lwk2), stk(lwk3), stk(lwk4));
 
311
      break;
 
312
    }
 
313
  LhsVar(1) = Rhs+1;
 
314
  PutLhsVar();
 
315
  return 0;
 
316
}
 
317
 
 
318
static int intlsq_splin(char * fname)
 
319
{
 
320
  /*   interface code for [y, d] = lsq_splin(xd, yd [, wd], x)  */
 
321
 
 
322
  int minrhs=3, maxrhs=4, minlhs=1, maxlhs=2;
 
323
 
 
324
  int mxd, nxd, lxd, myd, nyd, lyd, mx, nx, lx, mwd, nwd, lwd;
 
325
  int ly, ld, lwork, ndata, n, one=1, mwork, ierr;
 
326
  double un=1.0;
 
327
 
 
328
  CheckRhs(minrhs,maxrhs);
 
329
  CheckLhs(minlhs,maxlhs);
 
330
 
 
331
  GetRhsVar(1,"d", &mxd, &nxd, &lxd);
 
332
  GetRhsVar(2,"d", &myd, &nyd, &lyd);
 
333
  ndata = mxd*nxd;  /* number of data points */
 
334
  if ( ndata < 4  ||  mxd != myd  || nxd != nyd  ||  mxd != 1  &&  nxd != 1 ) 
 
335
    { 
 
336
      Scierror(999,"%s: arg 1 and 2 must be vectors of same size with at least %d elements\r\n",
 
337
               fname, 4);
 
338
      return 0;
 
339
    }
 
340
 
 
341
  if ( Rhs == 4 ) 
 
342
    {
 
343
      GetRhsVar(3,"d", &mwd, &nwd, &lwd);
 
344
      if ( mxd != mwd  ||  nxd != nwd )
 
345
        { 
 
346
          Scierror(999,"%s: bad input for arg 3\r\n", fname);
 
347
          return 0;
 
348
        }
 
349
    }
 
350
  GetRhsVar(Rhs,"d", &mx, &nx, &lx);
 
351
  n = mx*nx;
 
352
  if ( n < 2  ||  mx != 1  &&  nx != 1 )
 
353
    { 
 
354
      Scierror(999,"%s: bad input for x \r\n", fname);
 
355
      return 0;
 
356
    }
 
357
  
 
358
  if (! good_order(stk(lx), n))   /* verify strict increasing abscissae */
 
359
    {
 
360
      Scierror(999,"%s: elts of arg %d not (strictly) increasing or +-inf detected\r\n", fname, Rhs);
 
361
      return 0;
 
362
    }
 
363
 
 
364
  CreateVar(Rhs+1, "d", &mx,  &nx,   &ly);
 
365
  CreateVar(Rhs+2, "d", &mx,  &nx,   &ld);
 
366
  mwork = 7*n+18;
 
367
  CreateVar(Rhs+3, "d", &mwork, &one, &lwork); 
 
368
  if ( Rhs == 3 )
 
369
    {
 
370
      CreateVar(Rhs+4, "d", &mxd, &nxd, &lwd);
 
371
      C2F(dset)( &ndata, &un, stk(lwd), &one);  /* set all the weight = 1  */ 
 
372
    }
 
373
 
 
374
  C2F(spfit)(stk(lxd), stk(lyd), stk(lwd), &ndata, stk(lx), &n, stk(ly),
 
375
             stk(ld), stk(lwork), &ierr);
 
376
 
 
377
  if (ierr == -1)
 
378
    {
 
379
      Scierror(999,"%s: not enought points for the fit \r\n", fname);
 
380
      return 0;
 
381
    }
 
382
  else if (ierr == 1)
 
383
    sciprint("%s warning: rank deficiency of the least square matrix\r\n", fname);
 
384
 
 
385
  LhsVar(1) = Rhs+1;
 
386
  LhsVar(2) = Rhs+2;
 
387
  PutLhsVar();
 
388
  return 0;
 
389
}
 
390
 
 
391
static int intinterp1(char * fname)
 
392
{
 
393
  int minrhs=4, maxrhs=5, minlhs=1, maxlhs=4;
 
394
 
 
395
  int mt, nt, lt, mx, nx, lx, my, ny, ly, md, nd, ld, ns, *str_outmode;
 
396
  int n, m, outmode, lst, ldst, lddst, ldddst;
 
397
/*  double *x;*/
 
398
 
 
399
  CheckRhs(minrhs,maxrhs);
 
400
  CheckLhs(minlhs,maxlhs);
 
401
 
 
402
  GetRhsVar(1,"d", &mt, &nt, &lt);
 
403
  GetRhsVar(2,"d", &mx, &nx, &lx);
 
404
  GetRhsVar(3,"d", &my, &ny, &ly);
 
405
  GetRhsVar(4,"d", &md, &nd, &ld);
 
406
 
 
407
  if ( mx != my  ||  nx != ny  ||  md != mx || nd != nx || mx != 1  &&  nx != 1 || mx*nx < 2) 
 
408
    { 
 
409
      Scierror(999,"%s: bad inputs \r\n", fname);
 
410
      return 0;
 
411
    }
 
412
  n = mx*nx;    /* number of interpolation points */
 
413
  m = mt*nt;    /* number of points to interpolate */
 
414
 
 
415
  if ( Rhs == 5 )   /* get the outmode */
 
416
    {
 
417
      GetRhsScalarString(5, &ns, &str_outmode);
 
418
      outmode =  get_type(OutModeTable, NB_OUTMODE, str_outmode, ns);
 
419
      if ( outmode == UNDEFINED )
 
420
        {
 
421
          Scierror(999,"%s: unknown outmode type\n\r",fname);
 
422
          return 0;
 
423
        };
 
424
    }
 
425
  else
 
426
    outmode = C0;  /* default outmode */
 
427
 
 
428
  /* memory for st, dst, ddst, dddst */
 
429
  CreateVar(Rhs+1, "d", &mt,  &nt, &lst);
 
430
  CreateVar(Rhs+2, "d", &mt,  &nt, &ldst);
 
431
  CreateVar(Rhs+3, "d", &mt,  &nt, &lddst);
 
432
  CreateVar(Rhs+4, "d", &mt,  &nt, &ldddst);
 
433
 
 
434
  /*      subroutine EvalPWHermite(t, st, dst, ddst, dddst, m, x, y, d, n, outmode)
 
435
   *      integer m, n, outmode
 
436
   *      double precision t(m), st(m), dst(m), ddst(m), dddst(m), x(n), y(n), d(n)
 
437
   */
 
438
  C2F(evalpwhermite) (stk(lt), stk(lst), stk(ldst), stk(lddst), stk(ldddst),
 
439
                      &m, stk(lx), stk(ly), stk(ld), &n, &outmode);
 
440
 
 
441
  LhsVar(1) = Rhs+1;
 
442
  LhsVar(2) = Rhs+2;
 
443
  LhsVar(3) = Rhs+3;
 
444
  LhsVar(4) = Rhs+4;
 
445
  PutLhsVar();
 
446
  return 0;
 
447
}
 
448
 
 
449
static int intlinear_interpn(char * fname)
 
450
{
 
451
/*  interpolation lineaire n-dimensionnelle
 
452
 *
 
453
 *   yp = linear_interpn(xp1, ..., xpn, x1, ..., xn, val, outmode)
 
454
 */
 
455
  int n, mxp, nxp, lxp, mxpn, nxpn, lxpn, mx, nx, lx, my, ny, ly, one=1;
 
456
  int ns, *str_outmode, np, *k, *ad, m, l, i, outmode;
 
457
  int *dim;
 
458
  double **xp, **x, *val, *u, *v, *yp;
 
459
  RealHyperMat U;
 
460
 
 
461
  n = (Rhs+1)/2 - 1;
 
462
  if ( n < 1 )
 
463
    { 
 
464
      Scierror(999,"%s: too few arg \r\n", fname);
 
465
      return 0;
 
466
    }
 
467
 
 
468
  /*  les points sur lesquels on evalue par interpolation */
 
469
/*   l = I_UINT32; CreateVar(Rhs+1, "I", &n, &one, &l); */
 
470
/*   xp = (double **) istk(l); */
 
471
  CreateVar(Rhs+1, "d", &n, &one, &l);  /* => lets store an array of pointers  */
 
472
  xp = (double **) stk(l);               /*   with size of 4 or 8 bytes */
 
473
 
 
474
  GetRhsVar(1,"d", &mxp, &nxp, &lxp);
 
475
  xp[0] = stk(lxp);
 
476
  np = mxp*nxp;
 
477
  for ( i = 2 ; i <= n ; i++ )
 
478
    {
 
479
      GetRhsVar(i,"d", &mxpn, &nxpn, &lxpn);
 
480
      if ( mxp != mxpn || nxp != nxpn )
 
481
        { 
 
482
          Scierror(999,"%s: bad inputs for xp1, xp2, ...., \r\n", fname);
 
483
          return 0;
 
484
        }
 
485
      xp[i-1] = stk(lxpn);
 
486
    }
 
487
 
 
488
  /* coordonn�es de la grille */
 
489
  l = I_INT32; CreateVar(Rhs+2, "I", &n, &one, &l);
 
490
  dim = istk(l);
 
491
/*   l = I_UINT32; CreateVar(Rhs+3, "I", &n, &one, &l); */
 
492
/*   x = (double **) istk(l); */
 
493
  CreateVar(Rhs+3, "d", &n, &one, &l);  /* => lets store an array of pointers  */
 
494
  x = (double **) stk(l);               /*    with size(void *) = 4 or 8 bytes */
 
495
 
 
496
  for ( i = 1 ; i <= n ; i++ )
 
497
    {
 
498
      GetRhsVar(n+i,"d", &mx, &nx, &lx);
 
499
      if ( (mx != 1 && nx != 1) && mx*nx < 2)
 
500
        { 
 
501
          Scierror(999,"%s: bad arg number %d \r\n", fname, n+i);
 
502
          return 0;
 
503
        }
 
504
      x[i-1] = stk(lx);
 
505
      dim[i-1] = mx*nx;
 
506
      /* verify strict increasing order  */
 
507
      if ( !good_order(x[i-1], mx*nx) )
 
508
        {
 
509
          Scierror(999,"%s: grid abscissae of dim %d not in strict increasing order \r\n", fname, n+i);
 
510
          return 0;
 
511
        }
 
512
    }
 
513
 
 
514
  /* les valeurs de la grille */
 
515
  if ( n >= 3 )
 
516
    {
 
517
      GetRhsRealHMat(2*n+1,&U);
 
518
      if ( U.dimsize != n )
 
519
        { 
 
520
          Scierror(999,"%s: U must be a real %d-dim hypermatrix  \n", fname, n);
 
521
          return 0;
 
522
        }
 
523
      for ( i = 0 ; i < n ; i++ )
 
524
        if ( U.dims[i] != dim[i] )
 
525
          { 
 
526
            Scierror(999,"%s: size incompatibility between grid points and grid values in dim %d \n\r", fname, i+1);
 
527
            return 0;
 
528
          }
 
529
      val = U.R;
 
530
    }
 
531
  else  /* n = 1 or 2 */
 
532
    {
 
533
      GetRhsVar(2*n+1,"d", &my, &ny, &ly);
 
534
      if ( n == 1  &&  my*ny != dim[0] )
 
535
        { 
 
536
          Scierror(999,"%s: size incompatibility between grid points and values in dim 1 \n\r", fname);
 
537
          return 0;
 
538
        }
 
539
      if ( n == 2  &&  (my != dim[0]  || ny != dim[1]) )
 
540
        { 
 
541
          Scierror(999,"%s: size incompatibility between grid points and values in dim 1 or 2 \n\r", fname);
 
542
          return 0;
 
543
        }
 
544
      val = stk(ly);
 
545
    }
 
546
 
 
547
  /* get the outmode */
 
548
  if ( Rhs == 2*n + 2 )
 
549
    {
 
550
      GetRhsScalarString(Rhs, &ns, &str_outmode);
 
551
      outmode =  get_type(OutModeTable, NB_OUTMODE, str_outmode, ns);
 
552
      if ( outmode == UNDEFINED || outmode == LINEAR )
 
553
        {
 
554
          Scierror(999,"%s: unsupported outmode type\n\r",fname);
 
555
          return 0;
 
556
        };
 
557
    }
 
558
  else
 
559
    outmode = C0;
 
560
 
 
561
  CreateVar(Rhs+4, "d", &n, &one, &l); u = stk(l);
 
562
  m = 1; for ( i = 1 ; i <= n ; i++) m = 2*m;
 
563
  CreateVar(Rhs+5, "d", &m, &one, &l); v = stk(l);
 
564
  l = 4; CreateVar(Rhs+6, "I", &m, &one, &l); ad = istk(l);
 
565
  l = 4; CreateVar(Rhs+7, "I", &n, &one, &l); k = istk(l);
 
566
  CreateVar(Rhs+8, "d", &mxp, &nxp, &l); yp = stk(l);
 
567
 
 
568
  nlinear_interp(x, val, dim, n, xp, yp, np, outmode, u, v, ad, k);
 
569
 
 
570
  LhsVar(1) = Rhs+8;
 
571
  PutLhsVar();
 
572
  /* correction Warning Allan CORNET */
 
573
  /* warning C4715: 'intlinear_interpn' : not all control paths return a value */
 
574
  return 0;
 
575
}
 
576
 
 
577
 
 
578
static int intsplin2d(char * fname)
 
579
{
 
580
  /*    interface pour splin2d :
 
581
   *
 
582
   *    C = splin2d(x, y, z [, type])
 
583
   *
 
584
   */
 
585
 
 
586
  int minrhs=3, maxrhs=4, minlhs=1, maxlhs=1;
 
587
 
 
588
  int mx, nx, lx, my, ny, ly, mz, nz, lz, ns, mc, nc, lc, lp, lq, lr;
 
589
  int spline_type, *str_spline_type/*, i*/;
 
590
  int one = 1;
 
591
  double *x, *y, *C;
 
592
 
 
593
  CheckRhs(minrhs,maxrhs);
 
594
  CheckLhs(minlhs,maxlhs);
 
595
 
 
596
  GetRhsVar(1,"d", &mx, &nx, &lx);
 
597
  GetRhsVar(2,"d", &my, &ny, &ly);
 
598
  GetRhsVar(3,"d", &mz, &nz, &lz);
 
599
 
 
600
  if ( mx != 1 || my != 1 || mz != nx || nz != ny || nx < 2 || ny < 2)
 
601
    { 
 
602
      Scierror(999,"%s: bad inputs \r\n", fname);
 
603
      return 0;
 
604
    }
 
605
 
 
606
  /* verify strict increasing order for x and y */
 
607
  x = stk(lx);  y = stk(ly);
 
608
  if ( !good_order(x, nx) || !good_order(y, ny))
 
609
    {
 
610
      Scierror(999,"%s: x and/or y are not in strict increasing order (or +-inf detected) \r\n", fname);
 
611
      return 0;
 
612
    }
 
613
 
 
614
  /* get the spline type */
 
615
  if ( Rhs == 4 ) 
 
616
    {
 
617
      GetRhsScalarString(4, &ns, &str_spline_type);
 
618
      spline_type = get_type(SplineTable, NB_SPLINE_TYPE, str_spline_type, ns);
 
619
      if ( spline_type == UNDEFINED || spline_type == CLAMPED )
 
620
        {
 
621
          Scierror(999,"%s: unsupported spline type\n\r",fname);
 
622
          return 0;
 
623
        };
 
624
    }
 
625
  else
 
626
    spline_type = NOT_A_KNOT;
 
627
 
 
628
  /* memory for the big C array */
 
629
  mc = 16*(nx-1)*(ny-1); nc = 1;
 
630
  CreateVar( Rhs+1, "d", &mc,  &nc, &lc);
 
631
  C = stk(lc);
 
632
  /* memory for work arrays  */
 
633
  CreateVar( Rhs+2, "d", &nx, &ny, &lp);
 
634
  CreateVar( Rhs+3, "d", &nx, &ny, &lq);
 
635
  CreateVar( Rhs+4, "d", &nx, &ny, &lr);
 
636
 
 
637
  if (spline_type == MONOTONE || spline_type == FAST || spline_type == FAST_PERIODIC)
 
638
    {   
 
639
      C2F(bicubicsubspline)(x, y, stk(lz), &nx, &ny, stk(lc), 
 
640
                            stk(lp), stk(lq), stk(lr), &spline_type);
 
641
    }
 
642
 
 
643
  else   /*  spline */
 
644
    {
 
645
      int lA_d, lA_sd, ld, lqdu, lutemp, nxy, nxym1, nxym2, lll;
 
646
 
 
647
      nxy = max(nx,ny); nxym1 = nxy-1; nxym2 = nxy-2; 
 
648
 
 
649
      CreateVar( Rhs+5, "d", &nxy,   &one, &lA_d);
 
650
      CreateVar( Rhs+6, "d", &nxym1, &one, &lA_sd);
 
651
      CreateVar( Rhs+7, "d", &ny,    &one, &ld);
 
652
      CreateVar( Rhs+8, "d", &nxy,   &one, &lqdu);
 
653
      CreateVar( Rhs+9, "d", &ny,    &one, &lutemp);
 
654
 
 
655
      if (spline_type == PERIODIC)
 
656
        {
 
657
          CreateVar(Rhs+10, "d", &nxym2, &one, &lll) ;
 
658
        }
 
659
      else
 
660
        lll = lA_sd ;   /* bidon ... */
 
661
      C2F(bicubicspline)(x, y, stk(lz), &nx, &ny, stk(lc), stk(lp), stk(lq), stk(lr), 
 
662
                         stk(lA_d), stk(lA_sd), stk(ld), stk(lll), stk(lqdu), 
 
663
                         stk(lutemp), &spline_type);
 
664
    }
 
665
 
 
666
  LhsVar(1) = Rhs+1;
 
667
  PutLhsVar();
 
668
  return 0;
 
669
}
 
670
 
 
671
static int intinterp2d(char * fname)
 
672
{
 
673
  /*    interface pour interp2d :
 
674
   *
 
675
   *    [zp [, dzdxp, dzdyp [, d2zdx2, d2zdxy, d2zdy2]]] = interp2d(xp, yp, x, y, C[, outmode])
 
676
   */
 
677
 
 
678
  int minrhs=5, maxrhs=6, minlhs=1, maxlhs=6;
 
679
 
 
680
  int mxp, nxp, lxp, myp, nyp, lyp, mx, nx, lx, my, ny, ly;
 
681
  int mc, nc, lc, ns, *str_outmode, lzp, ldzdxp, ldzdyp, ld2zdx2p, ld2zdxyp, ld2zdy2p;
 
682
  int outmode, m;
 
683
 
 
684
  CheckRhs(minrhs,maxrhs);
 
685
  CheckLhs(minlhs,maxlhs);
 
686
 
 
687
  GetRhsVar(1,"d", &mxp, &nxp, &lxp);
 
688
  GetRhsVar(2,"d", &myp, &nyp, &lyp);
 
689
  GetRhsVar(3,"d", &mx, &nx, &lx);
 
690
  GetRhsVar(4,"d", &my, &ny, &ly);
 
691
  GetRhsVar(5,"d", &mc, &nc, &lc);
 
692
 
 
693
  if ( mxp != myp || nxp != nyp || mx != 1 || my != 1 || nc != 1 || nx < 2 || ny < 2
 
694
       || mc != 16*(nx-1)*(ny-1) )
 
695
    { 
 
696
      Scierror(999,"%s: bad inputs \r\n", fname);
 
697
      return 0;
 
698
    }
 
699
 
 
700
  /* get the outmode */
 
701
  if ( Rhs == 6 ) 
 
702
    {
 
703
      GetRhsScalarString(6, &ns, &str_outmode);
 
704
      outmode =  get_type(OutModeTable, NB_OUTMODE, str_outmode, ns);
 
705
      if ( outmode == UNDEFINED || outmode == LINEAR )
 
706
        {
 
707
          Scierror(999,"%s: unsupported outmode type\n\r",fname);
 
708
          return 0;
 
709
        };
 
710
    }
 
711
  else
 
712
    outmode = C0;
 
713
 
 
714
  /* memory for zp */
 
715
  CreateVar( Rhs+1, "d", &mxp,  &nxp, &lzp);
 
716
  m = mxp*nxp;
 
717
 
 
718
  if ( Lhs == 1 )
 
719
    {
 
720
      /*   subroutine BiCubicInterp(x, y, C, nx, ny, x_eval, y_eval, z_eval, m, outmode)
 
721
       *     integer nx, ny, m, outmode
 
722
       *     double precision x(nx), y(ny), C(4,4,nx-1,ny-1), x_eval(m), y_eval(m), z_eval(m)
 
723
       */
 
724
      C2F(bicubicinterp)(stk(lx), stk(ly), stk(lc), &nx, &ny, stk(lxp), stk(lyp), stk(lzp), 
 
725
                         &m, &outmode);
 
726
      LhsVar(1) = Rhs+1;
 
727
    }
 
728
  else   /* got also the derivatives */
 
729
    {
 
730
      /* memory for dzdxp and dzdyp */
 
731
      CreateVar( Rhs+2, "d", &mxp,  &nxp, &ldzdxp);
 
732
      CreateVar( Rhs+3, "d", &mxp,  &nxp, &ldzdyp);
 
733
 
 
734
      if ( Lhs <= 3 )
 
735
        {
 
736
          C2F(bicubicinterpwithgrad)(stk(lx), stk(ly), stk(lc), &nx, &ny, stk(lxp), 
 
737
                                     stk(lyp), stk(lzp), stk(ldzdxp), stk(ldzdyp), 
 
738
                                     &m, &outmode);
 
739
          LhsVar(1) = Rhs+1;
 
740
          LhsVar(2) = Rhs+2;
 
741
          LhsVar(3) = Rhs+3;
 
742
        }
 
743
      else /* got also 2d derivatives */
 
744
        {
 
745
          CreateVar( Rhs+4, "d", &mxp,  &nxp, &ld2zdx2p);
 
746
          CreateVar( Rhs+5, "d", &mxp,  &nxp, &ld2zdxyp);
 
747
          CreateVar( Rhs+6, "d", &mxp,  &nxp, &ld2zdy2p);
 
748
          C2F(bicubicinterpwithgradandhes)(stk(lx), stk(ly), stk(lc), &nx, &ny, stk(lxp), 
 
749
                                           stk(lyp), stk(lzp), stk(ldzdxp), stk(ldzdyp), 
 
750
                                           stk(ld2zdx2p), stk(ld2zdxyp), stk(ld2zdy2p), 
 
751
                                           &m, &outmode);
 
752
          LhsVar(1) = Rhs+1;
 
753
          LhsVar(2) = Rhs+2;
 
754
          LhsVar(3) = Rhs+3;
 
755
          LhsVar(4) = Rhs+4;
 
756
          LhsVar(5) = Rhs+5;
 
757
          LhsVar(6) = Rhs+6;
 
758
        }
 
759
    }
 
760
  PutLhsVar();
 
761
  return 0;
 
762
}
 
763
 
 
764
/*  
 
765
 *   interface sur le package cshep2d ....
 
766
 */
 
767
 
 
768
static int intcshep2d(char * fname)
 
769
{
 
770
  static char *Str[]={"cshep2d", "xyz", "lcell", "lnext", "grdim", "rmax", "rw", "a"};
 
771
  int minrhs=1, maxrhs=1, minlhs=1, maxlhs=1;
 
772
  int n, dim, nc, nw, nr, one=1, four=4, eight=8, nine=9, ier;
 
773
  int lxyz, lxyzn, lcell, lnext, lgrid, lrmax, lrw, la, ltlist, lar;
 
774
 
 
775
  double * xyz, * grid;
 
776
 
 
777
  CheckRhs(minrhs,maxrhs);
 
778
  CheckLhs(minlhs,maxlhs);
 
779
 
 
780
  GetRhsVar(1,"d", &n, &dim, &lxyz);
 
781
  if ( dim != 3  ||  n < 10 ) 
 
782
    { 
 
783
      Scierror(999,"%s: xyz must be a (n,3) real matrix with n >= 10 \r\n", fname);
 
784
      return 0;
 
785
    }
 
786
 
 
787
  /* choix pour nc (peut etre futur parametre optionnel) */
 
788
  nc = min( 17, n-1 );
 
789
  /* choix pour nw */
 
790
  nw = min( 30, n-1 );
 
791
  /* choix pour nr (grille nr x nr) */
 
792
  nr = (int) sqrt( n/3.0 ); /* comme n >= 10 nr >= 1 */
 
793
 
 
794
  /* all the information for the "interpolant" will be stored
 
795
   * in a tlist (which also contains the entry xyz)  
 
796
   */
 
797
  CreateVar(2,"t", &eight, &one, &ltlist);
 
798
  CreateListVarFromPtr(2, 1, "S", &one,  &eight, Str);
 
799
  CreateListVarFrom(2, 2, "d", &n ,   &dim,  &lxyzn, &lxyz);
 
800
  lcell = 4; lar = -1;
 
801
  CreateListVarFrom(2, 3, "I", &nr,   &nr,   &lcell, &lar); /* lcell */
 
802
  lnext = 4; lar = -1;
 
803
  CreateListVarFrom(2, 4, "I", &one,  &n,    &lnext, &lar); /* lnext */
 
804
  lar = -1;
 
805
  CreateListVarFrom(2, 5, "d", &one,  &four, &lgrid, &lar); /* xmin, ymin, dx, dy */
 
806
  lar = -1;
 
807
  CreateListVarFrom(2, 6, "d", &one,  &one,  &lrmax, &lar); /* rmax */
 
808
  lar = -1;
 
809
  CreateListVarFrom(2, 7, "d", &one,  &n,    &lrw,   &lar); /* rw */
 
810
  lar = -1;
 
811
  CreateListVarFrom(2, 8, "d", &nine, &n,    &la,    &lar); /* a */
 
812
  grid = stk(lgrid);
 
813
  xyz = stk(lxyz);
 
814
  
 
815
  /*      SUBROUTINE CSHEP2 (N,X,Y,F,NC,NW,NR, LCELL,LNEXT,XMIN,
 
816
   *                         YMIN,DX,DY,RMAX,RW,A,IER)
 
817
   */
 
818
  C2F(cshep2) ( &n, xyz, &xyz[n], &xyz[2*n], &nc, &nw, &nr, istk(lcell),
 
819
                istk(lnext), grid, &grid[1], &grid[2], &grid[3], stk(lrmax),
 
820
                stk(lrw), stk(la), &ier);
 
821
 
 
822
  if ( ier != 0 )
 
823
    {
 
824
      Scierror(999,"%s: duplicate nodes or all nodes colinears (ier = %d) \r\n", fname, ier);
 
825
      return 0;
 
826
    }
 
827
 
 
828
  LhsVar(1) = 2;
 
829
  PutLhsVar();
 
830
  return 0;
 
831
}
 
832
 
 
833
static int inteval_cshep2d(char * fname)
 
834
{
 
835
  /*
 
836
   *   [f [,dfdx, dfdy [, dffdxx, dffdxy, dffdyy]]] = eval_cshep2d(xp, yp, tlcoef)  
 
837
   */
 
838
 
 
839
  int minrhs=3, maxrhs=3, minlhs=1, maxlhs=6;
 
840
  int mx, nx, lx, my, ny, ly, mt, nt, lt;
 
841
  char **Str;
 
842
  int m1, n1, m2, n2, m3, n3, m4, n4, m5, n5, m6, n6, m7, n7, m8, n8;
 
843
  int lxyz, lgrid, lrmax, lrw, la;
 
844
  double *xp, *yp, *xyz, *grid, *f, *dfdx, *dfdy, *dffdxx, *dffdyy, *dffdxy;
 
845
  int i, ier, n, np, nr, lf, ldfdx, ldfdy, ldffdxx, ldffdyy, ldffdxy;
 
846
  SciIntMat Cell, Next;
 
847
  int *cell, *next;
 
848
 
 
849
  CheckRhs(minrhs,maxrhs);
 
850
  CheckLhs(minlhs,maxlhs);
 
851
 
 
852
  GetRhsVar(1,"d", &mx, &nx, &lx);
 
853
  GetRhsVar(2,"d", &my, &ny, &ly);
 
854
  if ( mx != my  ||  nx != ny ) 
 
855
    { 
 
856
      Scierror(999,"%s: xp and yp must have the same dimension \r\n", fname);
 
857
      return 0;
 
858
    }
 
859
 
 
860
  GetRhsVar(3,"t",&mt, &nt, &lt);
 
861
  GetListRhsVar(3, 1, "S", &m1,  &n1, &Str);    /* m1 = 1, n1 = 8 ? a verifier */
 
862
  if ( strcmp(Str[0],"cshep2d") != 0) 
 
863
    {
 
864
      FreeRhsSVar(Str);
 
865
      Scierror(999,"%s: Argument 2 is not an cshep2d tlist\r\n", fname);
 
866
      return 0;
 
867
    }
 
868
  FreeRhsSVar(Str);  
 
869
  GetListRhsVar(3, 2, "d", &m2, &n2,  &lxyz);   /* m2 = n , n2 = 3  */
 
870
  GetListRhsVar(3, 3, "I", &m3, &n3,  (int *)&Cell);  /* m3 = nr, n3 = nr */
 
871
  GetListRhsVar(3, 4, "I", &m4, &n4,  (int *)&Next);  /* m4 = 1 , n4 = n  */
 
872
  GetListRhsVar(3, 5, "d", &m5, &n5,  &lgrid);  /* m5 = 1 , n5 = 4  */
 
873
  GetListRhsVar(3, 6, "d", &m6, &n6,  &lrmax);  /* m6 = 1 , n6 = 1  */
 
874
  GetListRhsVar(3, 7, "d", &m7, &n7,  &lrw);    /* m7 = 1 , n7 = n  */
 
875
  GetListRhsVar(3, 8, "d", &m8, &n8,  &la);     /* m8 = 9 , n8 = n  */
 
876
 
 
877
  cell = (int *)Cell.D; next = (int *)Next.D; 
 
878
  xp = stk(lx); yp = stk(ly); np = mx*nx; 
 
879
  n = m2; nr = m3;
 
880
  xyz = stk(lxyz); grid = stk(lgrid);
 
881
 
 
882
  CreateVar(4, "d", &mx, &nx, &lf); f = stk(lf);
 
883
  if ( Lhs > 1 ) 
 
884
    { 
 
885
      CreateVar(5, "d", &mx, &nx, &ldfdx); dfdx = stk(ldfdx);
 
886
      CreateVar(6, "d", &mx, &nx, &ldfdy); dfdy = stk(ldfdy);
 
887
    }
 
888
  if ( Lhs > 3 ) 
 
889
    {
 
890
      CreateVar(7, "d", &mx, &nx, &ldffdxx); dffdxx = stk(ldffdxx);
 
891
      CreateVar(8, "d", &mx, &nx, &ldffdxy); dffdyy = stk(ldffdxy);
 
892
      CreateVar(9, "d", &mx, &nx, &ldffdyy); dffdxy = stk(ldffdyy);
 
893
    }
 
894
 
 
895
  switch( Lhs )
 
896
    {
 
897
    case ( 1 ) :
 
898
      for ( i = 0 ; i < np ; i++ )
 
899
        /*            DOUBLE PRECISION FUNCTION CS2VAL (PX,PY,N,X,Y,F,NR,
 
900
         *                          LCELL,LNEXT,XMIN,YMIN,DX,DY,RMAX,RW,A)
 
901
         */
 
902
        f[i] = C2F(cs2val)(&xp[i], &yp[i], &n, xyz, &xyz[n], &xyz[2*n], &nr,
 
903
                           cell, next, grid, &grid[1], &grid[2], &grid[3],
 
904
                           stk(lrmax), stk(lrw), stk(la));
 
905
      LhsVar(1) = 4;
 
906
      break;
 
907
 
 
908
    case ( 2 ) :
 
909
    case ( 3 ) :
 
910
      for ( i = 0 ; i < np ; i++ )
 
911
        /*      SUBROUTINE CS2GRD (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
 
912
         *.                   YMIN,DX,DY,RMAX,RW,A, C,CX,CY,IER)
 
913
         */
 
914
        C2F(cs2grd) (&xp[i], &yp[i], &n, xyz, &xyz[n], &xyz[2*n], &nr,
 
915
                     cell, next, grid, &grid[1], &grid[2], &grid[3],
 
916
                     stk(lrmax), stk(lrw), stk(la), &f[i], &dfdx[i], &dfdy[i], &ier);
 
917
      LhsVar(1) = 4;
 
918
      LhsVar(2) = 5;
 
919
      LhsVar(3) = 6;
 
920
      break;
 
921
 
 
922
    case ( 4 ) :
 
923
    case ( 5 ) :
 
924
    case ( 6 ) :
 
925
      for ( i = 0 ; i < np ; i++ )
 
926
        {
 
927
        /*   SUBROUTINE CS2HES (PX,PY,N,X,Y,F,NR,LCELL,LNEXT,XMIN,
 
928
         *.                     YMIN,DX,DY,RMAX,RW,A, C,CX,CY,CXX,CXY,CYY,IER)
 
929
         */
 
930
        C2F(cs2hes) (&xp[i], &yp[i], &n, xyz, &xyz[n], &xyz[2*n], &nr,
 
931
                     cell, next, grid, &grid[1], &grid[2], &grid[3],
 
932
                     stk(lrmax), stk(lrw), stk(la), &f[i], &dfdx[i], &dfdy[i], 
 
933
                     &dffdxx[i], &dffdxy[i], &dffdyy[i], &ier);
 
934
        }
 
935
      LhsVar(1) = 4;
 
936
      LhsVar(2) = 5;
 
937
      LhsVar(3) = 6;
 
938
      LhsVar(4) = 7;
 
939
      LhsVar(5) = 8;
 
940
      LhsVar(6) = 9;
 
941
      break;
 
942
    } 
 
943
  PutLhsVar();
 
944
  return 0;
 
945
}
 
946
 
 
947
static int intsplin3d(char * fname)
 
948
{
 
949
  /*
 
950
   *   [tlist] = splin3d(x, y, z, v [,orderxyz])
 
951
   */
 
952
  static char *Str[]={"tensbs3d", "tx", "ty", "tz", "order", "bcoef", "xyzminmax"};
 
953
 
 
954
  int minrhs=4, maxrhs=5, minlhs=1, maxlhs=1;
 
955
 
 
956
  int mx, nx, lx, my, ny, ly, mz, nz, lz, mo, no, lo, kx, ky, kz;
 
957
  int ntx, nty, ntz, ltx, lty, ltz, lbcoef, lxyzminmax, mwk, mwkx, mwky, mwkz;
 
958
  int flag, one=1, three=3, six=6, seven=7,ltlist, nxyz, lwork, lar, lorder, *order;
 
959
  double *x, *y, *z, *xyzminmax;
 
960
  RealHyperMat V;
 
961
 
 
962
  CheckRhs(minrhs,maxrhs);
 
963
  CheckLhs(minlhs,maxlhs);
 
964
 
 
965
  GetRhsVar(1,"d", &mx, &nx, &lx);
 
966
  CheckVector(1, mx, nx); x = stk(lx);
 
967
  GetRhsVar(2,"d", &my, &ny, &ly);
 
968
  CheckVector(2, my, ny); y = stk(ly);
 
969
  GetRhsVar(3,"d", &mz, &nz, &lz);
 
970
  CheckVector(2, mz, nz); z = stk(lz);
 
971
 
 
972
  nx = mx*nx; ny = my*ny; nz = mz*nz;
 
973
 
 
974
  if ( nx < 3  ||  ny < 3  ||  nz < 3 ) 
 
975
    { 
 
976
      Scierror(999,"%s: the x, y and z grids must have at least 3 points \r\n", fname);
 
977
      return 0;
 
978
    }
 
979
 
 
980
  GetRhsRealHMat(4, &V);
 
981
  if ( V.dimsize != 3 )
 
982
    { 
 
983
      Scierror(999,"%s: 4 th argument must be a real 3-dim hypermatrix  \n", fname);
 
984
      return 0;
 
985
    }
 
986
  if ( V.dims[0] != nx  ||  V.dims[1] != ny  ||  V.dims[2] != nz  )
 
987
    { 
 
988
      Scierror(999,"%s: size incompatibility between grid points and grid values\n\r", fname);
 
989
      return 0;
 
990
    }
 
991
 
 
992
  if ( Rhs == 5 )
 
993
    {
 
994
      GetRhsVar(5,"d", &mo, &no, &lo);
 
995
      if ( (mo != 1 && no != 1)  ||  mo*no != 3 )
 
996
        { 
 
997
          Scierror(999,"%s: the 4 th arg must be a vector with 3 components \r\n", fname);
 
998
          return 0;
 
999
        }
 
1000
      kx = (int)*stk(lo); ky = (int)*stk(lo+1); kz = (int)*stk(lo+2);
 
1001
      if ( kx < 2  ||  kx >= nx  ||  ky < 2  ||  ky >= ny  ||  kz < 2  ||  kz >= nz )
 
1002
        { 
 
1003
          Scierror(999,"%s: bad 5 th arg [kx ky kz]\n\r", fname);
 
1004
          return 0;
 
1005
        }
 
1006
    }
 
1007
  else
 
1008
    {
 
1009
      kx = 4; ky = 4; kz = 4;
 
1010
    }
 
1011
 
 
1012
  ntx = nx + kx;
 
1013
  nty = ny + ky;
 
1014
  ntz = nz + kz;
 
1015
  mwkx = kx*(nx+1); mwky = ky*(ny+1); mwkz = kz*(nz+1); 
 
1016
  mwkx = max(mwkx, mwky);
 
1017
  mwk = nx*ny*nz + 2*(max(mwkx, mwkz));
 
1018
  nxyz = nx*ny*nz;
 
1019
 
 
1020
  CreateVar(Rhs+1,"t", &seven, &one, &ltlist);
 
1021
  CreateListVarFromPtr(Rhs+1, 1, "S", &one,  &seven, Str);
 
1022
  lar = -1; CreateListVarFrom(Rhs+1, 2, "d", &ntx, &one, &ltx, &lar);
 
1023
  lar = -1; CreateListVarFrom(Rhs+1, 3, "d", &nty, &one, &lty, &lar);
 
1024
  lar = -1; CreateListVarFrom(Rhs+1, 4, "d", &ntz, &one, &ltz, &lar);
 
1025
  lorder = 4; 
 
1026
  lar = -1; CreateListVarFrom(Rhs+1, 5, "I", &three, &one, &lorder, &lar);
 
1027
  order = istk(lorder); order[0] = kx; order[1] = ky; order[2] = kz;
 
1028
  lar = -1; CreateListVarFrom(Rhs+1, 6, "d", &nxyz,  &one, &lbcoef, &lar);
 
1029
  lar = -1; CreateListVarFrom(Rhs+1, 7, "d", &six,  &one, &lxyzminmax, &lar); 
 
1030
  xyzminmax = stk(lxyzminmax); 
 
1031
  xyzminmax[0] = x[0]; xyzminmax[1] = x[nx-1];  
 
1032
  xyzminmax[2] = y[0]; xyzminmax[3] = y[ny-1];  
 
1033
  xyzminmax[4] = z[0]; xyzminmax[5] = z[nz-1];  
 
1034
  CreateVar(Rhs+2, "d", &mwk, &one, &lwork);    /* work */
 
1035
 
 
1036
  flag = 0;
 
1037
  C2F(db3ink) ( stk(lx), &nx, stk(ly), &ny, stk(lz), &nz, V.R,
 
1038
                &nx, &ny, &kx, &ky, &kz, stk(ltx), stk(lty), stk(ltz), 
 
1039
                stk(lbcoef), stk(lwork), &flag);
 
1040
 
 
1041
  if ( flag != 1 )
 
1042
    {
 
1043
      Scierror(999,"%s: problem : flag = %d \r\n", fname, flag);
 
1044
      return 0;
 
1045
    }
 
1046
 
 
1047
  /*  Return only the tlist  */
 
1048
  LhsVar(1) = Rhs+1;
 
1049
  PutLhsVar();
 
1050
  return 0;
 
1051
}
 
1052
 
 
1053
static int intinterp3d(char * fname)  /* a suivre */
 
1054
{
 
1055
  /*
 
1056
   *   [f [, dfdx, dfdy, dfdz]] = interp3d(xp, yp, zp, tlcoef [,outmode])  
 
1057
   */
 
1058
 
 
1059
  int minrhs=4, maxrhs=5, minlhs=1, maxlhs=4;
 
1060
 
 
1061
  int mxp, nxp, lxp, myp, nyp, lyp, mzp, nzp, lzp, mt, nt, lt, np;
 
1062
  int zero=0, one=1, kx, ky, kz;
 
1063
  int nx, ny, nz, nxyz, mtx, mty, mtz, m, n, ltx, lty, ltz, lbcoef, mwork, lwork, lfp;
 
1064
  int lxyzminmax, nsix, outmode, ns, *str_outmode;
 
1065
  int /*i,*/ m1, n1, ldfpdx, ldfpdy, ldfpdz;
 
1066
  double *fp, *xp, *yp, *zp, *dfpdx, *dfpdy, *dfpdz;
 
1067
  double *xyzminmax, xmin, xmax, ymin, ymax, zmin, zmax;
 
1068
  SciIntMat Order; int *order;
 
1069
  char **Str;
 
1070
 
 
1071
  CheckRhs(minrhs,maxrhs);
 
1072
  CheckLhs(minlhs,maxlhs);
 
1073
 
 
1074
  GetRhsVar(1,"d", &mxp, &nxp, &lxp); xp = stk(lxp);
 
1075
  GetRhsVar(2,"d", &myp, &nyp, &lyp); yp = stk(lyp);
 
1076
  GetRhsVar(3,"d", &mzp, &nzp, &lzp); zp = stk(lzp);
 
1077
  if ( mxp != myp  ||  nxp != nyp || mxp != mzp  ||  nxp != nzp) 
 
1078
    { 
 
1079
      Scierror(999,"%s: xp, yp and zp must have the same dimensions \r\n", fname);
 
1080
      return 0;
 
1081
    }
 
1082
  np = mxp * nxp;
 
1083
 
 
1084
  GetRhsVar(4,"t",&mt, &nt, &lt);
 
1085
  GetListRhsVar(4, 1, "S", &m1,  &n1, &Str);
 
1086
  if ( strcmp(Str[0],"tensbs3d") != 0) 
 
1087
    {
 
1088
      FreeRhsSVar(Str);
 
1089
      Scierror(999,"%s: 4 th argument is not an tensbs3d tlist \r\n", fname);
 
1090
      return 0;
 
1091
    }
 
1092
  FreeRhsSVar(Str);
 
1093
  GetListRhsVar(4, 2, "d", &mtx, &n,  &ltx);
 
1094
  GetListRhsVar(4, 3, "d", &mty, &n,  &lty);
 
1095
  GetListRhsVar(4, 4, "d", &mtz, &n,  &ltz);
 
1096
  GetListRhsVar(4, 5, "I", &m  , &n,  (int *)&Order);
 
1097
  GetListRhsVar(4, 6, "d", &nxyz,&n,  &lbcoef);
 
1098
  GetListRhsVar(4, 7, "d", &nsix,&n,  &lxyzminmax);
 
1099
  xyzminmax = stk(lxyzminmax);
 
1100
  xmin = xyzminmax[0];  xmax = xyzminmax[1]; 
 
1101
  ymin = xyzminmax[2];  ymax = xyzminmax[3]; 
 
1102
  zmin = xyzminmax[4];  zmax = xyzminmax[5]; 
 
1103
 
 
1104
 
 
1105
  /* get the outmode */
 
1106
  if ( Rhs == 5 ) 
 
1107
    {
 
1108
      GetRhsScalarString(5, &ns, &str_outmode);
 
1109
      outmode =  get_type(OutModeTable, NB_OUTMODE, str_outmode, ns);
 
1110
      if ( outmode == UNDEFINED || outmode == LINEAR || outmode == NATURAL )
 
1111
        {
 
1112
          Scierror(999,"%s: unsupported outmode type\n\r",fname);
 
1113
          return 0;
 
1114
        };
 
1115
    }
 
1116
  else
 
1117
    outmode = C0;
 
1118
 
 
1119
  CreateVar(Rhs+1, "d", &mxp, &nxp, &lfp); fp = stk(lfp);
 
1120
 
 
1121
  order = (int *)Order.D;
 
1122
  kx = order[0]; ky = order[1]; kz = order[2];
 
1123
  nx = mtx - kx; ny = mty - ky; nz = mtz - kz; 
 
1124
 
 
1125
  mwork = ky*kz + 3*max(kx,max(ky,kz)) + kz;
 
1126
  CreateVar(Rhs+2, "d", &mwork, &one, &lwork);
 
1127
 
 
1128
  if ( Lhs == 1 )
 
1129
    {
 
1130
      C2F(driverdb3val)(xp,yp,zp,fp,&np,stk(ltx), stk(lty), stk(ltz),
 
1131
                        &nx, &ny, &nz, &kx, &ky, &kz, stk(lbcoef), stk(lwork),
 
1132
                        &xmin, &xmax, &ymin, &ymax, &zmin, &zmax, &outmode);
 
1133
      LhsVar(1) = Rhs+1;
 
1134
    }
 
1135
  else
 
1136
    {
 
1137
      CreateVar(Rhs+3, "d", &mxp, &nxp, &ldfpdx); dfpdx = stk(ldfpdx);
 
1138
      CreateVar(Rhs+4, "d", &mxp, &nxp, &ldfpdy); dfpdy = stk(ldfpdy);
 
1139
      CreateVar(Rhs+5, "d", &mxp, &nxp, &ldfpdz); dfpdz = stk(ldfpdz);
 
1140
      C2F(driverdb3valwithgrad)(xp,yp,zp,fp,dfpdx, dfpdy, dfpdz, &np,
 
1141
                                stk(ltx), stk(lty), stk(ltz),
 
1142
                                &nx, &ny, &nz, &kx, &ky, &kz, stk(lbcoef), stk(lwork),
 
1143
                                &xmin, &xmax, &ymin, &ymax, &zmin, &zmax, &outmode);
 
1144
      LhsVar(1) = Rhs+1;
 
1145
      LhsVar(2) = Rhs+3;
 
1146
      LhsVar(3) = Rhs+4;
 
1147
      LhsVar(4) = Rhs+5;
 
1148
    }
 
1149
 
 
1150
  PutLhsVar();
 
1151
  return 0;
 
1152
}
 
1153
 
 
1154
static int intbsplin3val(char * fname)
 
1155
{
 
1156
  /*
 
1157
   *   [fp] = bsplin3val(xp, yp, zp, tlcoef, der)  
 
1158
   */
 
1159
 
 
1160
  int minrhs=5, maxrhs=5, minlhs=1, maxlhs=1;
 
1161
 
 
1162
  int mxp, nxp, lxp, myp, nyp, lyp, mzp, nzp, lzp, mt, nt, lt, m1, n1, np;
 
1163
  int zero=0, one=1, kx, ky, kz;
 
1164
  int nx, ny, nz, nxyz, mtx, mty, mtz, m, n, ltx, lty, ltz, lbcoef, mwork, lwork, lfp;
 
1165
  int lxyzminmax, nsix;
 
1166
  int i, mder,nder,lder, ox, oy, oz;
 
1167
  double *fp, *xp, *yp, *zp, *der;
 
1168
  double *xyzminmax, xmin, xmax, ymin, ymax, zmin, zmax;
 
1169
  SciIntMat Order; int *order;
 
1170
  char **Str;
 
1171
 
 
1172
  CheckRhs(minrhs,maxrhs);
 
1173
  CheckLhs(minlhs,maxlhs);
 
1174
 
 
1175
  GetRhsVar(1,"d", &mxp, &nxp, &lxp); xp = stk(lxp);
 
1176
  GetRhsVar(2,"d", &myp, &nyp, &lyp); yp = stk(lyp);
 
1177
  GetRhsVar(3,"d", &mzp, &nzp, &lzp); zp = stk(lzp);
 
1178
  if ( mxp != myp  ||  nxp != nyp || mxp != mzp  ||  nxp != nzp) 
 
1179
    { 
 
1180
      Scierror(999,"%s: xp, yp and zp must have the same dimensions \r\n", fname);
 
1181
      return 0;
 
1182
    }
 
1183
  np = mxp * nxp;
 
1184
 
 
1185
  GetRhsVar(4,"t",&mt, &nt, &lt);
 
1186
  GetListRhsVar(4, 1, "S", &m1,  &n1, &Str);
 
1187
  if ( strcmp(Str[0],"tensbs3d") != 0) 
 
1188
    {
 
1189
      FreeRhsSVar(Str);
 
1190
      Scierror(999,"%s: 4 th argument is not an tensbs3d tlist \r\n", fname);
 
1191
      return 0;
 
1192
    }
 
1193
  FreeRhsSVar(Str);
 
1194
  GetListRhsVar(4, 2, "d", &mtx, &n,  &ltx);
 
1195
  GetListRhsVar(4, 3, "d", &mty, &n,  &lty);
 
1196
  GetListRhsVar(4, 4, "d", &mtz, &n,  &ltz);
 
1197
  GetListRhsVar(4, 5, "I", &m  , &n,  (int *)&Order);
 
1198
  GetListRhsVar(4, 6, "d", &nxyz,&n,  &lbcoef);
 
1199
  GetListRhsVar(4, 7, "d", &nsix,&n,  &lxyzminmax);
 
1200
  xyzminmax = stk(lxyzminmax);
 
1201
  xmin = xyzminmax[0];  xmax = xyzminmax[1]; 
 
1202
  ymin = xyzminmax[2];  ymax = xyzminmax[3]; 
 
1203
  zmin = xyzminmax[4];  zmax = xyzminmax[5]; 
 
1204
 
 
1205
  GetRhsVar(5,"d", &mder, &nder, &lder);
 
1206
  der = stk(lder);
 
1207
  if (   mder*nder != 3
 
1208
      || der[0] != floor(der[0]) || der[0] < 0.0 
 
1209
      || der[1] != floor(der[1]) || der[1] < 0.0 
 
1210
      || der[2] != floor(der[2]) || der[2] < 0.0 )
 
1211
    {
 
1212
      Scierror(999,"%s: bad 5 th argument \r\n", fname);
 
1213
      return 0;
 
1214
    }
 
1215
  ox = (int) der[0];  oy = (int) der[1];  oz = (int) der[2];
 
1216
 
 
1217
 
 
1218
  CreateVar(Rhs+1, "d", &mxp, &nxp, &lfp); fp = stk(lfp);
 
1219
 
 
1220
  order = (int *)Order.D;
 
1221
  kx = order[0]; ky = order[1]; kz = order[2];
 
1222
  nx = mtx - kx; ny = mty - ky; nz = mtz - kz; 
 
1223
 
 
1224
  mwork = ky*kz + 3*max(kx,max(ky,kz)) + kz;
 
1225
  CreateVar(Rhs+2, "d", &mwork, &one, &lwork);
 
1226
 
 
1227
  for ( i=0; i<np ; i++ )
 
1228
    {
 
1229
      fp[i] = C2F(db3val)(&(xp[i]), &(yp[i]), &(zp[i]), &ox, &oy, &oz, 
 
1230
                          stk(ltx), stk(lty), stk(lty), &nx, &ny, &nz,
 
1231
                          &kx, &ky, &kz, stk(lbcoef), stk(lwork));
 
1232
    }
 
1233
 
 
1234
  LhsVar(1) = Rhs+1;
 
1235
  PutLhsVar();
 
1236
  return 0;
 
1237
}
 
1238
 
 
1239
static TabF Tab[]={ 
 
1240
  {intsplin,           "splin"},
 
1241
  {intlsq_splin,       "lsq_splin"},
 
1242
  {intinterp1,          "interp"},
 
1243
  {intlinear_interpn,  "linear_interpn"},
 
1244
  {intsplin2d,         "splin2d"},
 
1245
  {intinterp2d,        "interp2d"},
 
1246
  {intcshep2d,         "cshep2d"},
 
1247
  {inteval_cshep2d,    "eval_cshep2d" },
 
1248
  {intsplin3d,         "splin3d"},
 
1249
  {intinterp3d,        "interp3d"},
 
1250
  {intbsplin3val,      "bsplin3val"}
 
1251
};
 
1252
 
 
1253
int C2F(intinterp)(void)
 
1254
{
 
1255
  Rhs = Max(0, Rhs);
 
1256
  (*(Tab[Fin-1].f))(Tab[Fin-1].name);
 
1257
  return 0;
 
1258
}