0022786: 64bit issue in AdvApp2Var
[occt.git] / src / AdvApp2Var / AdvApp2Var_MathBase.cxx
1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
2 //
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
7 //
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
10 //
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
17
18 // AdvApp2Var_MathBase.cxx
19 #include <math.h>
20 #include <AdvApp2Var_SysBase.hxx>
21 #include <AdvApp2Var_Data_f2c.hxx>
22 #include <AdvApp2Var_MathBase.hxx>
23 #include <AdvApp2Var_Data.hxx>
24
25 // statics 
26 static
27 int mmchole_(integer *mxcoef, 
28              integer *dimens, 
29              doublereal *amatri, 
30              integer *aposit, 
31              integer *posuiv, 
32              doublereal *chomat, 
33              integer *iercod);
34
35
36
37
38 static
39 int mmrslss_(integer *mxcoef, 
40              integer *dimens, 
41              doublereal *smatri, 
42              integer *sposit,
43              integer *posuiv, 
44              doublereal *mscnmbr,
45              doublereal *soluti, 
46              integer *iercod);
47
48 static
49 int mfac_(doublereal *f,
50           integer *n);
51
52 static
53 int mmaper0_(integer *ncofmx, 
54              integer *ndimen, 
55              integer *ncoeff, 
56              doublereal *crvlgd, 
57              integer *ncfnew, 
58              doublereal *ycvmax, 
59              doublereal *errmax);
60 static
61 int mmaper2_(integer *ncofmx,
62              integer *ndimen, 
63              integer *ncoeff, 
64              doublereal *crvjac, 
65              integer *ncfnew, 
66              doublereal *ycvmax, 
67              doublereal *errmax);
68
69 static
70 int mmaper4_(integer *ncofmx, 
71              integer *ndimen, 
72              integer *ncoeff, 
73              doublereal *crvjac, 
74              integer *ncfnew,
75              doublereal *ycvmax,
76              doublereal *errmax);
77
78 static
79 int mmaper6_(integer *ncofmx, 
80              integer *ndimen, 
81              integer *ncoeff, 
82              doublereal *crvjac, 
83              integer *ncfnew,
84              doublereal *ycvmax,
85              doublereal *errmax);
86
87 static
88 int mmarc41_(integer *ndimax, 
89              integer *ndimen, 
90              integer *ncoeff,
91              doublereal *crvold,
92              doublereal *upara0,
93              doublereal *upara1,
94              doublereal *crvnew,
95              integer *iercod);
96
97 static
98 int mmatvec_(integer *nligne, 
99              integer *ncolon,
100              integer *gposit,
101              integer *gnstoc, 
102              doublereal *gmatri,
103              doublereal *vecin, 
104              integer *deblig,
105              doublereal *vecout,
106              integer *iercod);
107
108 static
109 int mmcvstd_(integer *ncofmx, 
110              integer *ndimax, 
111              integer *ncoeff,
112              integer *ndimen, 
113              doublereal *crvcan, 
114              doublereal *courbe);
115
116 static
117 int mmdrvcb_(integer *ideriv,
118              integer *ndim, 
119              integer *ncoeff,
120              doublereal *courbe, 
121              doublereal *tparam,
122              doublereal *tabpnt, 
123              integer *iercod);
124
125 static
126 int mmexthi_(integer *ndegre, 
127              doublereal *hwgaus);
128
129 static
130 int mmextrl_(integer *ndegre,
131              doublereal *rootlg);
132
133
134
135 static
136 int mmherm0_(doublereal *debfin, 
137              integer *iercod);
138
139 static
140 int mmherm1_(doublereal *debfin, 
141              integer *ordrmx, 
142              integer *iordre, 
143              doublereal *hermit, 
144              integer *iercod);
145 static
146 int mmloncv_(integer *ndimax,
147              integer *ndimen,
148              integer *ncoeff,
149              doublereal *courbe, 
150              doublereal *tdebut, 
151              doublereal *tfinal, 
152              doublereal *xlongc, 
153              integer *iercod);
154 static
155 int mmpojac_(doublereal *tparam, 
156              integer *iordre, 
157              integer *ncoeff, 
158              integer *nderiv, 
159              doublereal *valjac, 
160              integer *iercod);
161
162 static
163 int mmrslw_(integer *normax, 
164             integer *nordre, 
165             integer *ndimen, 
166             doublereal *epspiv,
167             doublereal *abmatr,
168             doublereal *xmatri, 
169             integer *iercod);
170 static
171 int mmtmave_(integer *nligne, 
172              integer *ncolon, 
173              integer *gposit, 
174              integer *gnstoc, 
175              doublereal *gmatri,
176              doublereal *vecin, 
177              doublereal *vecout, 
178              integer *iercod);
179 static
180 int mmtrpj0_(integer *ncofmx,
181              integer *ndimen, 
182              integer *ncoeff, 
183              doublereal *epsi3d, 
184              doublereal *crvlgd, 
185              doublereal *ycvmax, 
186              doublereal *epstrc, 
187              integer *ncfnew);
188 static
189 int mmtrpj2_(integer *ncofmx,
190              integer *ndimen, 
191              integer *ncoeff, 
192              doublereal *epsi3d, 
193              doublereal *crvlgd, 
194              doublereal *ycvmax, 
195              doublereal *epstrc, 
196              integer *ncfnew);
197
198 static
199 int mmtrpj4_(integer *ncofmx,
200              integer *ndimen, 
201              integer *ncoeff, 
202              doublereal *epsi3d, 
203              doublereal *crvlgd, 
204              doublereal *ycvmax, 
205              doublereal *epstrc, 
206              integer *ncfnew);
207 static
208 int mmtrpj6_(integer *ncofmx,
209              integer *ndimen, 
210              integer *ncoeff, 
211              doublereal *epsi3d, 
212              doublereal *crvlgd, 
213              doublereal *ycvmax, 
214              doublereal *epstrc, 
215              integer *ncfnew);
216 static
217 integer  pow__ii(integer *x, 
218                  integer *n);
219
220 static
221 int mvcvin2_(integer *ncoeff, 
222              doublereal *crvold, 
223              doublereal *crvnew,
224              integer *iercod);
225
226 static
227 int mvcvinv_(integer *ncoeff,
228              doublereal *crvold, 
229              doublereal *crvnew, 
230              integer *iercod);
231
232 static
233 int mvgaus0_(integer *kindic, 
234              doublereal *urootl, 
235              doublereal *hiltab, 
236              integer *nbrval, 
237              integer *iercod);
238 static
239 int mvpscr2_(integer *ncoeff, 
240              doublereal *curve2, 
241              doublereal *tparam, 
242              doublereal *pntcrb);
243
244 static
245 int mvpscr3_(integer *ncoeff, 
246              doublereal *curve2, 
247              doublereal *tparam, 
248              doublereal *pntcrb);
249
250 static struct {
251     doublereal eps1, eps2, eps3, eps4;
252     integer niterm, niterr;
253 } mmprcsn_;
254
255 static struct {
256     doublereal tdebut, tfinal, verifi, cmherm[576];     
257 } mmcmher_;
258
259 //=======================================================================
260 //function : AdvApp2Var_MathBase::mdsptpt_
261 //purpose  : 
262 //=======================================================================
263 int AdvApp2Var_MathBase::mdsptpt_(integer *ndimen, 
264                                   doublereal *point1, 
265                                   doublereal *point2, 
266                                   doublereal *distan)
267
268 {
269   static integer c__8 = 8;
270   /* System generated locals */
271   integer i__1;
272   doublereal d__1;
273
274   /* Local variables */
275   static integer i__;
276   static doublereal differ[100];
277   static integer  ier;
278   intptr_t iofset, j;
279
280 /* ********************************************************************** 
281 */
282
283 /*     FUNCTION : */
284 /*     ---------- */
285 /*        CALCULATE DISTANCE BETWEEN TWO POINTS */
286
287 /*     KEYWORDS : */
288 /*     ----------- */
289 /*        DISTANCE,POINT. */
290
291 /*     INPUT ARGUMENTS : */
292 /*     ------------------ */
293 /*        NDIMEN: Space Dimension. */
294 /*        POINT1: Table of coordinates of the 1st point. */
295 /*        POINT2: Table of coordinates of the 2nd point. */
296
297 /*     OUTPUT ARGUMENTS : */
298 /*     ------------------- */
299 /*        DISTAN: Distance between 2 points. */
300
301 /*     COMMONS USED   : */
302 /*     ---------------- */
303
304 /*     REFERENCES CALLED   : */
305 /*     ----------------------- */
306
307 /*     DESCRIPTION/NOTES/LIMITATIONS : */
308 /*     ----------------------------------- */
309 /* > */
310 /* ********************************************************************** 
311 */
312
313
314 /* ***********************************************************************
315  */
316 /*                      INITIALISATION */
317 /* ***********************************************************************
318  */
319
320     /* Parameter adjustment */
321     --point2;
322     --point1;
323
324     /* Function Body */
325     iofset = 0;
326     ier = 0;
327
328 /* ***********************************************************************
329  */
330 /*                     TRAITEMENT */
331 /* ***********************************************************************
332  */
333
334     if (*ndimen > 100) {
335         AdvApp2Var_SysBase::mcrrqst_(&c__8, ndimen, differ, &iofset, &ier);
336     }
337
338 /* --- If allocation is refused, the trivial method is applied. */
339
340     if (ier > 0) {
341
342         *distan = 0.;
343         i__1 = *ndimen;
344         for (i__ = 1; i__ <= i__1; ++i__) {
345 /* Computing 2nd power */
346             d__1 = point1[i__] - point2[i__];
347             *distan += d__1 * d__1;
348         }
349         *distan = sqrt(*distan);
350
351 /* --- Otherwise MZSNORM is used to minimize the risks of overflow 
352 */
353
354     } else {
355         i__1 = *ndimen;
356         for (i__ = 1; i__ <= i__1; ++i__) {
357             j=iofset + i__ - 1;
358             differ[j] = point2[i__] - point1[i__];
359         }
360
361         *distan = AdvApp2Var_MathBase::mzsnorm_(ndimen, &differ[iofset]);
362
363     }
364
365 /* ***********************************************************************
366  */
367 /*                   RETURN CALLING PROGRAM */
368 /* ***********************************************************************
369  */
370
371 /* --- Dynamic Desallocation */
372
373     if (iofset != 0) {
374         AdvApp2Var_SysBase::mcrdelt_(&c__8, ndimen, differ, &iofset, &ier);
375     }
376
377  return 0 ;
378 } /* mdsptpt_ */
379
380 //=======================================================================
381 //function : mfac_
382 //purpose  : 
383 //=======================================================================
384 int mfac_(doublereal *f, 
385           integer *n)
386
387 {
388     /* System generated locals */
389     integer i__1;
390
391     /* Local variables */
392     static integer i__;
393
394 /*    FORTRAN CONFORME AU TEXT */
395 /*     CALCUL DE MFACTORIEL N */
396     /* Parameter adjustments */
397     --f;
398
399     /* Function Body */
400     f[1] = (float)1.;
401     i__1 = *n;
402     for (i__ = 2; i__ <= i__1; ++i__) {
403 /* L10: */
404         f[i__] = i__ * f[i__ - 1];
405     }
406     return 0;
407 } /* mfac_ */
408
409 //=======================================================================
410 //function : AdvApp2Var_MathBase::mmapcmp_
411 //purpose  : 
412 //=======================================================================
413 int AdvApp2Var_MathBase::mmapcmp_(integer *ndim, 
414                                   integer *ncofmx, 
415                                   integer *ncoeff, 
416                                   doublereal *crvold, 
417                                   doublereal *crvnew)
418
419 {
420   /* System generated locals */
421   integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1, 
422   i__2;
423
424   /* Local variables */
425   static integer ipair, nd, ndegre, impair, ibb, idg;
426   //extern  int  mgsomsg_();//mgenmsg_(),
427
428 /* ********************************************************************** 
429 */
430
431 /*     FUNCTION : */
432 /*     ---------- */
433 /*        Compression of curve CRVOLD in a table of  */
434 /*        coeff. of even : CRVNEW(*,0,*) */
435 /*        and uneven range : CRVNEW(*,1,*). */
436
437 /*     KEYWORDS : */
438 /*     ----------- */
439 /*        COMPRESSION,CURVE. */
440
441 /*     INPUT ARGUMENTS : */
442 /*     ------------------ */
443 /*     NDIM   : Space Dimension. */
444 /*     NCOFMX : Max nb of coeff. of the curve to compress. */
445 /*     NCOEFF : Max nb of coeff. of the compressed curve. */
446 /*     CRVOLD : The curve (0:NCOFMX-1,NDIM) to compress. */
447
448 /*     OUTPUT ARGUMENTS : */
449 /*     ------------------- */
450 /*     CRVNEW : Curve compacted in (0:(NCOEFF-1)/2,0,NDIM) (containing 
451 */
452 /*              even terms) and in (0:(NCOEFF-1)/2,1,NDIM) */
453 /*              (containing uneven terms). */
454
455 /*     COMMONS USED   : */
456 /*     ---------------- */
457
458 /*     REFERENCES CALLED   : */
459 /*     ----------------------- */
460
461 /*     DESCRIPTION/NOTES/LIMITATIONS : */
462 /*     ----------------------------------- */
463 /*     This routine is useful to prepare coefficients of a */
464 /*     curve in an orthogonal base (Legendre or Jacobi) before */
465 /*     calculating the coefficients in the canonical; base [-1,1] by */
466 /*     MMJACAN. */
467 /* ***********************************************************************
468  */
469
470 /*   Name of the routine */
471
472     /* Parameter adjustments */
473     crvold_dim1 = *ncofmx;
474     crvold_offset = crvold_dim1;
475     crvold -= crvold_offset;
476     crvnew_dim1 = (*ncoeff - 1) / 2 + 1;
477     crvnew_offset = crvnew_dim1 << 1;
478     crvnew -= crvnew_offset;
479
480     /* Function Body */
481     ibb = AdvApp2Var_SysBase::mnfndeb_();
482     if (ibb >= 3) {
483         AdvApp2Var_SysBase::mgenmsg_("MMAPCMP", 7L);
484     }
485
486     ndegre = *ncoeff - 1;
487     i__1 = *ndim;
488     for (nd = 1; nd <= i__1; ++nd) {
489         ipair = 0;
490         i__2 = ndegre / 2;
491         for (idg = 0; idg <= i__2; ++idg) {
492             crvnew[idg + (nd << 1) * crvnew_dim1] = crvold[ipair + nd * 
493                     crvold_dim1];
494             ipair += 2;
495 /* L200: */
496         }
497         if (ndegre < 1) {
498             goto L400;
499         }
500         impair = 1;
501         i__2 = (ndegre - 1) / 2;
502         for (idg = 0; idg <= i__2; ++idg) {
503             crvnew[idg + ((nd << 1) + 1) * crvnew_dim1] = crvold[impair + nd *
504                      crvold_dim1];
505             impair += 2;
506 /* L300: */
507         }
508
509 L400:
510 /* L100: */
511         ;
512     }
513
514 /* ---------------------------------- The end --------------------------- 
515 */
516
517     if (ibb >= 3) {
518         AdvApp2Var_SysBase::mgsomsg_("MMAPCMP", 7L);
519     }
520     return 0;
521 } /* mmapcmp_ */
522
523 //=======================================================================
524 //function : mmaper0_
525 //purpose  : 
526 //=======================================================================
527 int mmaper0_(integer *ncofmx, 
528              integer *ndimen, 
529              integer *ncoeff, 
530              doublereal *crvlgd, 
531              integer *ncfnew, 
532              doublereal *ycvmax, 
533              doublereal *errmax)
534
535 {
536   /* System generated locals */
537   integer crvlgd_dim1, crvlgd_offset, i__1, i__2;
538   doublereal d__1;
539
540   /* Local variables */
541   static integer ncut;
542   static doublereal bidon;
543   static integer ii, nd;
544
545 /* ***********************************************************************
546  */
547
548 /*     FUNCTION : */
549 /*     ---------- */
550 /*        Calculate the max error of approximation done when */
551 /*        only the first NCFNEW coefficients of a curve are preserved.  
552 */
553 /*        Degree NCOEFF-1 written in the base of Legendre (Jacobi */
554 /*        of  order 0). */
555
556 /*     KEYWORDS : */
557 /*     ----------- */
558 /*        LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
559
560 /*     INPUT ARGUMENTS : */
561 /*     ------------------ */
562 /*        NCOFMX : Max. degree of the curve. */
563 /*        NDIMEN : Space dimension. */
564 /*        NCOEFF : Degree +1 of the curve. */
565 /*        CRVLGD : Curve the degree which of should be lowered. */
566 /*        NCFNEW : Degree +1 of the resulting polynom. */
567
568 /*     OUTPUT ARGUMENTS : */
569 /*     ------------------- */
570 /*        YCVMAX : Auxiliary Table (max error on each dimension). 
571 */
572 /*        ERRMAX : Precision of the approximation. */
573
574 /*     COMMONS USED   : */
575 /*     ---------------- */
576
577 /*     REFERENCES CALLED   : */
578 /*     ----------------------- */
579
580 /*     DESCRIPTION/NOTES/LIMITATIONS : */
581 /*     ----------------------------------- */
582 /* ***********************************************************************
583  */
584
585
586 /* ------------------- Init to calculate an error ----------------------- 
587 */
588
589     /* Parameter adjustments */
590     --ycvmax;
591     crvlgd_dim1 = *ncofmx;
592     crvlgd_offset = crvlgd_dim1 + 1;
593     crvlgd -= crvlgd_offset;
594
595     /* Function Body */
596     i__1 = *ndimen;
597     for (ii = 1; ii <= i__1; ++ii) {
598         ycvmax[ii] = 0.;
599 /* L100: */
600     }
601
602 /* ------ Minimum that can be reached : Stop at 1 or NCFNEW ------ 
603 */
604
605     ncut = 1;
606     if (*ncfnew + 1 > ncut) {
607         ncut = *ncfnew + 1;
608     }
609
610 /* -------------- Elimination of high degree coefficients----------- 
611 */
612 /* ----------- Loop on the series of Legendre: NCUT --> NCOEFF -------- 
613 */
614
615     i__1 = *ncoeff;
616     for (ii = ncut; ii <= i__1; ++ii) {
617 /*   Factor of renormalization (Maximum of Li(t)). */
618         bidon = ((ii - 1) * 2. + 1.) / 2.;
619         bidon = sqrt(bidon);
620
621         i__2 = *ndimen;
622         for (nd = 1; nd <= i__2; ++nd) {
623             ycvmax[nd] += (d__1 = crvlgd[ii + nd * crvlgd_dim1], advapp_abs(d__1)) * 
624                     bidon;
625 /* L310: */
626         }
627 /* L300: */
628     }
629
630 /* -------------- The error is the norm of the vector error --------------- 
631 */
632
633     *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
634
635 /* --------------------------------- Fin -------------------------------- 
636 */
637
638     return 0;
639 } /* mmaper0_ */
640
641 //=======================================================================
642 //function : mmaper2_
643 //purpose  : 
644 //=======================================================================
645 int mmaper2_(integer *ncofmx,
646              integer *ndimen, 
647              integer *ncoeff, 
648              doublereal *crvjac, 
649              integer *ncfnew, 
650              doublereal *ycvmax, 
651              doublereal *errmax)
652
653 {
654   /* Initialized data */
655
656     static doublereal xmaxj[57] = { .9682458365518542212948163499456,
657             .986013297183269340427888048593603,
658             1.07810420343739860362585159028115,
659             1.17325804490920057010925920756025,
660             1.26476561266905634732910520370741,
661             1.35169950227289626684434056681946,
662             1.43424378958284137759129885012494,
663             1.51281316274895465689402798226634,
664             1.5878364329591908800533936587012,
665             1.65970112228228167018443636171226,
666             1.72874345388622461848433443013543,
667             1.7952515611463877544077632304216,
668             1.85947199025328260370244491818047,
669             1.92161634324190018916351663207101,
670             1.98186713586472025397859895825157,
671             2.04038269834980146276967984252188,
672             2.09730119173852573441223706382076,
673             2.15274387655763462685970799663412,
674             2.20681777186342079455059961912859,
675             2.25961782459354604684402726624239,
676             2.31122868752403808176824020121524,
677             2.36172618435386566570998793688131,
678             2.41117852396114589446497298177554,
679             2.45964731268663657873849811095449,
680             2.50718840313973523778244737914028,
681             2.55385260994795361951813645784034,
682             2.59968631659221867834697883938297,
683             2.64473199258285846332860663371298,
684             2.68902863641518586789566216064557,
685             2.73261215675199397407027673053895,
686             2.77551570192374483822124304745691,
687             2.8177699459714315371037628127545,
688             2.85940333797200948896046563785957,
689             2.90044232019793636101516293333324,
690             2.94091151970640874812265419871976,
691             2.98083391718088702956696303389061,
692             3.02023099621926980436221568258656,
693             3.05912287574998661724731962377847,
694             3.09752842783622025614245706196447,
695             3.13546538278134559341444834866301,
696             3.17295042316122606504398054547289,
697             3.2099992681699613513775259670214,
698             3.24662674946606137764916854570219,
699             3.28284687953866689817670991319787,
700             3.31867291347259485044591136879087,
701             3.35411740487202127264475726990106,
702             3.38919225660177218727305224515862,
703             3.42390876691942143189170489271753,
704             3.45827767149820230182596660024454,
705             3.49230918177808483937957161007792,
706             3.5260130200285724149540352829756,
707             3.55939845146044235497103883695448,
708             3.59247431368364585025958062194665,
709             3.62524904377393592090180712976368,
710             3.65773070318071087226169680450936,
711             3.68992700068237648299565823810245,
712             3.72184531357268220291630708234186 };
713
714     /* System generated locals */
715     integer crvjac_dim1, crvjac_offset, i__1, i__2;
716     doublereal d__1;
717
718     /* Local variables */
719     static integer idec, ncut;
720     static doublereal bidon;
721     static integer ii, nd;
722
723
724
725 /* ***********************************************************************
726  */
727
728 /*     FONCTION : */
729 /*     ---------- */
730 /*        Calculate max approximation error i faite lorsque l' on */
731 /*        ne conserve que les premiers NCFNEW coefficients d' une courbe 
732 */
733 /*        de degre NCOEFF-1 ecrite dans la base de Jacobi d' ordre 2. */
734
735 /*     KEYWORDS : */
736 /*     ----------- */
737 /*        JACOBI, POLYGON, APPROXIMATION, ERROR. */
738 /*
739 /*  INPUT ARGUMENTS : */
740 /*     ------------------ */
741 /*        NCOFMX : Max. degree of the curve. */
742 /*        NDIMEN : Space dimension. */
743 /*        NCOEFF : Degree +1 of the curve. */
744 /*        CRVLGD : Curve the degree which of should be lowered. */
745 /*        NCFNEW : Degree +1 of the resulting polynom. */
746
747 /*     OUTPUT ARGUMENTS : */
748 /*     ------------------- */
749 /*        YCVMAX : Auxiliary Table (max error on each dimension). 
750 */
751 /*        ERRMAX : Precision of the approximation. */
752
753 /*     COMMONS USED   : */
754 /*     ---------------- */
755
756 /*     REFERENCES CALLED   : */
757 /*     ----------------------- */
758 /*     DESCRIPTION/NOTES/LIMITATIONS : */
759 /*     ----------------------------------- */
760
761
762
763 /* ------------------ Table of maximums of (1-t2)*Ji(t) ---------------- 
764 */
765
766     /* Parameter adjustments */
767     --ycvmax;
768     crvjac_dim1 = *ncofmx;
769     crvjac_offset = crvjac_dim1 + 1;
770     crvjac -= crvjac_offset;
771
772     /* Function Body */
773
774
775
776 /* ------------------- Init for error  calculation ----------------------- 
777 */
778
779     i__1 = *ndimen;
780     for (ii = 1; ii <= i__1; ++ii) {
781         ycvmax[ii] = 0.;
782 /* L100: */
783     }
784
785 /* ------ Min. Degree that can be attained : Stop at 3 or NCFNEW ------ 
786 */
787
788     idec = 3;
789 /* Computing MAX */
790     i__1 = idec, i__2 = *ncfnew + 1;
791     ncut = advapp_max(i__1,i__2);
792
793 /* -------------- Removal of coefficients of high degree ----------- 
794 */
795 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ---------- 
796 */
797
798     i__1 = *ncoeff;
799     for (ii = ncut; ii <= i__1; ++ii) {
800 /*   Factor of renormalization. */
801         bidon = xmaxj[ii - idec];
802         i__2 = *ndimen;
803         for (nd = 1; nd <= i__2; ++nd) {
804             ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) * 
805                     bidon;
806 /* L310: */
807         }
808 /* L300: */
809     }
810
811 /* -------------- The error is the norm of the vector error --------------- 
812 */
813
814     *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
815
816 /* --------------------------------- Fin -------------------------------- 
817 */
818
819     return 0;
820 } /* mmaper2_ */
821
822 /* MAPER4.f -- translated by f2c (version 19960827).
823    You must link the resulting object file with the libraries:
824         -lf2c -lm   (in that order)
825 */
826
827 /* Subroutine */ 
828 //=======================================================================
829 //function : mmaper4_
830 //purpose  : 
831 //=======================================================================
832 int mmaper4_(integer *ncofmx, 
833              integer *ndimen, 
834              integer *ncoeff, 
835              doublereal *crvjac, 
836              integer *ncfnew,
837              doublereal *ycvmax,
838              doublereal *errmax)
839 {
840     /* Initialized data */
841
842     static doublereal xmaxj[55] = { 1.1092649593311780079813740546678,
843             1.05299572648705464724876659688996,
844             1.0949715351434178709281698645813,
845             1.15078388379719068145021100764647,
846             1.2094863084718701596278219811869,
847             1.26806623151369531323304177532868,
848             1.32549784426476978866302826176202,
849             1.38142537365039019558329304432581,
850             1.43575531950773585146867625840552,
851             1.48850442653629641402403231015299,
852             1.53973611681876234549146350844736,
853             1.58953193485272191557448229046492,
854             1.63797820416306624705258190017418,
855             1.68515974143594899185621942934906,
856             1.73115699602477936547107755854868,
857             1.77604489805513552087086912113251,
858             1.81989256661534438347398400420601,
859             1.86276344480103110090865609776681,
860             1.90471563564740808542244678597105,
861             1.94580231994751044968731427898046,
862             1.98607219357764450634552790950067,
863             2.02556989246317857340333585562678,
864             2.06433638992049685189059517340452,
865             2.10240936014742726236706004607473,
866             2.13982350649113222745523925190532,
867             2.17661085564771614285379929798896,
868             2.21280102016879766322589373557048,
869             2.2484214321456956597803794333791,
870             2.28349755104077956674135810027654,
871             2.31805304852593774867640120860446,
872             2.35210997297725685169643559615022,
873             2.38568889602346315560143377261814,
874             2.41880904328694215730192284109322,
875             2.45148841120796359750021227795539,
876             2.48374387161372199992570528025315,
877             2.5155912654873773953959098501893,
878             2.54704548720896557684101746505398,
879             2.57812056037881628390134077704127,
880             2.60882970619319538196517982945269,
881             2.63918540521920497868347679257107,
882             2.66919945330942891495458446613851,
883             2.69888301230439621709803756505788,
884             2.72824665609081486737132853370048,
885             2.75730041251405791603760003778285,
886             2.78605380158311346185098508516203,
887             2.81451587035387403267676338931454,
888             2.84269522483114290814009184272637,
889             2.87060005919012917988363332454033,
890             2.89823818258367657739520912946934,
891             2.92561704377132528239806135133273,
892             2.95274375377994262301217318010209,
893             2.97962510678256471794289060402033,
894             3.00626759936182712291041810228171,
895             3.03267744830655121818899164295959,
896             3.05886060707437081434964933864149 };
897
898     /* System generated locals */
899     integer crvjac_dim1, crvjac_offset, i__1, i__2;
900     doublereal d__1;
901
902     /* Local variables */
903     static integer idec, ncut;
904     static doublereal bidon;
905     static integer ii, nd;
906
907
908
909 /* ***********************************************************************
910  */
911
912 /*     FUNCTION : */
913 /*     ---------- */
914 /*        Calculate the max. error of approximation made when  */
915 /*        only first NCFNEW coefficients of a curve are preserved 
916 */
917 /*        degree NCOEFF-1 is written in the base of Jacobi of order 4. */
918 /*        KEYWORDS : */
919 /*     ----------- */
920 /*        LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
921
922 /*     INPUT ARGUMENTS : */
923 /*     ------------------ */
924 /*        NCOFMX : Max. degree of the curve. */
925 /*        NDIMEN : Space dimension. */
926 /*        NCOEFF : Degree +1 of the curve. */
927 /*        CRVJAC : Curve the degree which of should be lowered. */
928 /*        NCFNEW : Degree +1 of the resulting polynom. */
929
930 /*     OUTPUT ARGUMENTS : */
931 /*     ------------------- */
932 /*        YCVMAX : Auxiliary Table (max error on each dimension). 
933 */
934 /*        ERRMAX : Precision of the approximation. */
935
936 /*     COMMONS USED   : */
937 /*     ---------------- */
938
939 /*     REFERENCES CALLED   : */
940 /*     ----------------------- */
941
942 /*     DESCRIPTION/NOTES/LIMITATIONS : */
943
944
945 /* ***********************************************************************
946  */
947
948
949 /* ---------------- Table of maximums of ((1-t2)2)*Ji(t) --------------- 
950 */
951
952     /* Parameter adjustments */
953     --ycvmax;
954     crvjac_dim1 = *ncofmx;
955     crvjac_offset = crvjac_dim1 + 1;
956     crvjac -= crvjac_offset;
957
958     /* Function Body */
959
960
961
962 /* ------------------- Init for error calculation ----------------------- 
963 */
964
965     i__1 = *ndimen;
966     for (ii = 1; ii <= i__1; ++ii) {
967         ycvmax[ii] = 0.;
968 /* L100: */
969     }
970
971 /* ------ Min. Degree that can be attained : Stop at 5 or NCFNEW ------ 
972 */
973
974     idec = 5;
975 /* Computing MAX */
976     i__1 = idec, i__2 = *ncfnew + 1;
977     ncut = advapp_max(i__1,i__2);
978
979 /* -------------- Removal of high degree coefficients ----------- 
980 */
981 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ---------- 
982 */
983
984     i__1 = *ncoeff;
985     for (ii = ncut; ii <= i__1; ++ii) {
986 /*   Factor of renormalisation. */
987         bidon = xmaxj[ii - idec];
988         i__2 = *ndimen;
989         for (nd = 1; nd <= i__2; ++nd) {
990             ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) * 
991                     bidon;
992 /* L310: */
993         }
994 /* L300: */
995     }
996
997 /* -------------- The error is the norm of the error vector --------------- 
998 */
999
1000     *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1001
1002 /* --------------------------------- End -------------------------------- 
1003 */
1004
1005     return 0;
1006 } /* mmaper4_ */
1007
1008 //=======================================================================
1009 //function : mmaper6_
1010 //purpose  : 
1011 //=======================================================================
1012 int mmaper6_(integer *ncofmx, 
1013              integer *ndimen, 
1014              integer *ncoeff, 
1015              doublereal *crvjac, 
1016              integer *ncfnew,
1017              doublereal *ycvmax,
1018              doublereal *errmax)
1019
1020 {
1021     /* Initialized data */
1022
1023     static doublereal xmaxj[53] = { 1.21091229812484768570102219548814,
1024             1.11626917091567929907256116528817,
1025             1.1327140810290884106278510474203,
1026             1.1679452722668028753522098022171,
1027             1.20910611986279066645602153641334,
1028             1.25228283758701572089625983127043,
1029             1.29591971597287895911380446311508,
1030             1.3393138157481884258308028584917,
1031             1.3821288728999671920677617491385,
1032             1.42420414683357356104823573391816,
1033             1.46546895108549501306970087318319,
1034             1.50590085198398789708599726315869,
1035             1.54550385142820987194251585145013,
1036             1.58429644271680300005206185490937,
1037             1.62230484071440103826322971668038,
1038             1.65955905239130512405565733793667,
1039             1.69609056468292429853775667485212,
1040             1.73193098017228915881592458573809,
1041             1.7671112206990325429863426635397,
1042             1.80166107681586964987277458875667,
1043             1.83560897003644959204940535551721,
1044             1.86898184653271388435058371983316,
1045             1.90180515174518670797686768515502,
1046             1.93410285411785808749237200054739,
1047             1.96589749778987993293150856865539,
1048             1.99721027139062501070081653790635,
1049             2.02806108474738744005306947877164,
1050             2.05846864831762572089033752595401,
1051             2.08845055210580131460156962214748,
1052             2.11802334209486194329576724042253,
1053             2.14720259305166593214642386780469,
1054             2.17600297710595096918495785742803,
1055             2.20443832785205516555772788192013,
1056             2.2325216999457379530416998244706,
1057             2.2602654243075083168599953074345,
1058             2.28768115912702794202525264301585,
1059             2.3147799369092684021274946755348,
1060             2.34157220782483457076721300512406,
1061             2.36806787963276257263034969490066,
1062             2.39427635443992520016789041085844,
1063             2.42020656255081863955040620243062,
1064             2.44586699364757383088888037359254,
1065             2.47126572552427660024678584642791,
1066             2.49641045058324178349347438430311,
1067             2.52130850028451113942299097584818,
1068             2.54596686772399937214920135190177,
1069             2.5703922285006754089328998222275,
1070             2.59459096001908861492582631591134,
1071             2.61856915936049852435394597597773,
1072             2.64233265984385295286445444361827,
1073             2.66588704638685848486056711408168,
1074             2.68923766976735295746679957665724,
1075             2.71238965987606292679677228666411 };
1076
1077     /* System generated locals */
1078     integer crvjac_dim1, crvjac_offset, i__1, i__2;
1079     doublereal d__1;
1080
1081     /* Local variables */
1082     static integer idec, ncut;
1083     static doublereal bidon;
1084     static integer ii, nd;
1085
1086
1087
1088 /* ***********************************************************************
1089  */
1090 /*     FUNCTION : */
1091 /*     ---------- */
1092 /*        Calculate the max. error of approximation made when  */
1093 /*        only first NCFNEW coefficients of a curve are preserved 
1094 */
1095 /*        degree NCOEFF-1 is written in the base of Jacobi of order 6. */
1096 /*        KEYWORDS : */
1097 /*     ----------- */
1098 /*        JACOBI,POLYGON,APPROXIMATION,ERROR. */
1099
1100 /*     INPUT ARGUMENTS : */
1101 /*     ------------------ */
1102 /*        NCOFMX : Max. degree of the curve. */
1103 /*        NDIMEN : Space dimension. */
1104 /*        NCOEFF : Degree +1 of the curve. */
1105 /*        CRVJAC : Curve the degree which of should be lowered. */
1106 /*        NCFNEW : Degree +1 of the resulting polynom. */
1107
1108 /*     OUTPUT ARGUMENTS : */
1109 /*     ------------------- */
1110 /*        YCVMAX : Auxiliary Table (max error on each dimension). 
1111 */
1112 /*        ERRMAX : Precision of the approximation. */
1113
1114 /*     COMMONS USED   : */
1115 /*     ---------------- */
1116
1117 /*     REFERENCES CALLED   : */
1118 /*     ----------------------- */
1119
1120 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1121 /* > */
1122 /* ***********************************************************************
1123  */
1124
1125
1126 /* ---------------- Table of maximums of ((1-t2)3)*Ji(t) --------------- 
1127 */
1128
1129     /* Parameter adjustments */
1130     --ycvmax;
1131     crvjac_dim1 = *ncofmx;
1132     crvjac_offset = crvjac_dim1 + 1;
1133     crvjac -= crvjac_offset;
1134
1135     /* Function Body */
1136
1137
1138
1139 /* ------------------- Init for error calculation ----------------------- 
1140 */
1141
1142     i__1 = *ndimen;
1143     for (ii = 1; ii <= i__1; ++ii) {
1144         ycvmax[ii] = 0.;
1145 /* L100: */
1146     }
1147
1148 /* ------ Min Degree that can be attained : Stop at 3 or NCFNEW ------ 
1149 */
1150
1151     idec = 7;
1152 /* Computing MAX */
1153     i__1 = idec, i__2 = *ncfnew + 1;
1154     ncut = advapp_max(i__1,i__2);
1155
1156 /* -------------- Removal of high degree coefficients ----------- 
1157 */
1158 /* ----------- Loop on the series of Jacobi :NCUT --> NCOEFF ---------- 
1159 */
1160
1161     i__1 = *ncoeff;
1162     for (ii = ncut; ii <= i__1; ++ii) {
1163 /*   Factor of renormalization. */
1164         bidon = xmaxj[ii - idec];
1165         i__2 = *ndimen;
1166         for (nd = 1; nd <= i__2; ++nd) {
1167             ycvmax[nd] += (d__1 = crvjac[ii + nd * crvjac_dim1], advapp_abs(d__1)) * 
1168                     bidon;
1169 /* L310: */
1170         }
1171 /* L300: */
1172     }
1173
1174 /* -------------- The error is the norm of the vector error --------------- 
1175 */
1176
1177     *errmax = AdvApp2Var_MathBase::mzsnorm_(ndimen, &ycvmax[1]);
1178
1179 /* --------------------------------- END -------------------------------- 
1180 */
1181
1182     return 0;
1183 } /* mmaper6_ */
1184
1185 //=======================================================================
1186 //function : AdvApp2Var_MathBase::mmaperx_
1187 //purpose  : 
1188 //=======================================================================
1189 int AdvApp2Var_MathBase::mmaperx_(integer *ncofmx, 
1190                                   integer *ndimen, 
1191                                   integer *ncoeff, 
1192                                   integer *iordre, 
1193                                   doublereal *crvjac, 
1194                                   integer *ncfnew, 
1195                                   doublereal *ycvmax, 
1196                                   doublereal *errmax, 
1197                                   integer *iercod)
1198
1199 {
1200   /* System generated locals */
1201   integer crvjac_dim1, crvjac_offset;
1202
1203   /* Local variables */
1204   static integer jord;
1205
1206 /* ********************************************************************** 
1207 */
1208 /*     FUNCTION : */
1209 /*     ---------- */
1210 /*        Calculate the max. error of approximation made when  */
1211 /*        only first NCFNEW coefficients of a curve are preserved 
1212 */
1213 /*        degree NCOEFF-1 is written in the base of Jacobi of order IORDRE. */
1214 /*        KEYWORDS : */
1215 /*     ----------- */
1216 /*        JACOBI,LEGENDRE,POLYGON,APPROXIMATION,ERROR. */
1217
1218 /*     INPUT ARGUMENTS : */
1219 /*     ------------------ */
1220 /*        NCOFMX : Max. degree of the curve. */
1221 /*        NDIMEN : Space dimension. */
1222 /*        NCOEFF : Degree +1 of the curve. */ 
1223 /*        IORDRE : Order of continuity at the extremities. */
1224 /*        CRVJAC : Curve the degree which of should be lowered. */
1225 /*        NCFNEW : Degree +1 of the resulting polynom. */
1226
1227 /*     OUTPUT ARGUMENTS : */
1228 /*     ------------------- */
1229 /*        YCVMAX : Auxiliary Table (max error on each dimension). 
1230 */
1231 /*        ERRMAX : Precision of the approximation. */
1232 /*        IERCOD = 0, OK */
1233 /*               = 1, order of constraints (IORDRE) is not within the */
1234 /*                    autorized values. */
1235 /*     COMMONS USED   : */
1236 /*     ---------------- */
1237
1238 /*     REFERENCES CALLED   : */
1239 /*     ----------------------- */
1240
1241 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1242 /*     ----------------------------------- */
1243 /*     Canceled and replaced MMAPERR. */
1244 /* ***********************************************************************
1245  */
1246
1247
1248     /* Parameter adjustments */
1249     --ycvmax;
1250     crvjac_dim1 = *ncofmx;
1251     crvjac_offset = crvjac_dim1 + 1;
1252     crvjac -= crvjac_offset;
1253
1254     /* Function Body */
1255     *iercod = 0;
1256 /* --> Order of Jacobi polynoms */
1257     jord = ( *iordre + 1) << 1;
1258
1259     if (jord == 0) {
1260         mmaper0_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1261                 ycvmax[1], errmax);
1262     } else if (jord == 2) {
1263         mmaper2_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1264                 ycvmax[1], errmax);
1265     } else if (jord == 4) {
1266         mmaper4_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1267                 ycvmax[1], errmax);
1268     } else if (jord == 6) {
1269         mmaper6_(ncofmx, ndimen, ncoeff, &crvjac[crvjac_offset], ncfnew, &
1270                 ycvmax[1], errmax);
1271     } else {
1272         *iercod = 1;
1273     }
1274
1275 /* ----------------------------------- Fin ------------------------------ 
1276 */
1277
1278     return 0;
1279 } /* mmaperx_ */
1280
1281 //=======================================================================
1282 //function : mmarc41_
1283 //purpose  : 
1284 //=======================================================================
1285  int mmarc41_(integer *ndimax, 
1286               integer *ndimen, 
1287               integer *ncoeff,
1288               doublereal *crvold,
1289               doublereal *upara0,
1290               doublereal *upara1,
1291               doublereal *crvnew,
1292               integer *iercod)
1293
1294 {
1295   /* System generated locals */
1296     integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1, 
1297     i__2, i__3;
1298     
1299     /* Local variables */
1300     static integer nboct;
1301     static doublereal tbaux[61];
1302     static integer nd;
1303     static doublereal bid;
1304     static integer ncf, ncj;
1305
1306
1307 /*      IMPLICIT DOUBLE PRECISION(A-H,O-Z) */
1308 /*      IMPLICIT INTEGER (I-N) */
1309
1310 /* ***********************************************************************
1311  */
1312
1313 /*     FUNCTION : */
1314 /*     ---------- */
1315 /*     Creation of curve C2(v) defined on (0,1) identic to */
1316 /*     curve C1(u) defined on (U0,U1) (change of parameter */
1317 /*     of a curve). */
1318
1319 /*     KEYWORDS : */
1320 /*     ----------- */
1321 /*        LIMITATION, RESTRICTION, CURVE */
1322
1323 /*     INPUT ARGUMENTS : */
1324 /*     ------------------ */
1325 /*   NDIMAX : Space Dimensioning. */
1326 /*   NDIMEN : Curve Dimension. */
1327 /*   NCOEFF : Nb of coefficients of the curve. */
1328 /*   CRVOLD : Curve to be limited. */
1329 /*   UPARA0     : Min limit of the interval limiting the curve. 
1330 */
1331 /*   UPARA1     : Max limit of the interval limiting the curve. 
1332 */
1333
1334 /*     OUTPUT ARGUMENTS : */
1335 /*     ------------------- */
1336 /*   CRVNEW : Relimited curve, defined on (0,1) and equal to */
1337 /*            CRVOLD defined on (U0,U1). */
1338 /*   IERCOD : = 0, OK */
1339 /*            =10, Nb of coeff. <1 or > 61. */
1340
1341 /*     COMMONS USED   : */
1342 /*     ---------------- */
1343 /*     REFERENCES CALLED   : */
1344 /*     ---------------------- */
1345 /*     Type  Name */
1346 /*           MAERMSG              MCRFILL              MVCVIN2 */
1347 /*           MVCVINV */
1348
1349 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1350 /*     ----------------------------------- */
1351 /* ---> Algorithm used in this general case is based on the */
1352 /*     following principle  : */
1353 /*        Let S(t) = a0 + a1*t + a2*t**2 + ... of degree NCOEFF-1, and */
1354 /*               U(t) = b0 + b1*t, then the coeff. of */
1355 /*        S(U(t)) are calculated step by step with help of table TBAUX. */
1356 /*        At each step number N (N=2 to NCOEFF), TBAUX(n) contains */
1357 /*        the n-th coefficient of U(t)**N for n=1 to N. (RBD) */
1358 /* ---> Reference : KNUTH, 'The Art of Computer Programming', */
1359 /*                        Vol. 2/'Seminumerical Algorithms', */
1360 /*                        Ex. 11 p:451 et solution p:562. (RBD) */
1361
1362 /* ---> Removal of the input argument CRVOLD by CRVNEW is */
1363 /*     possible, which means that the call : */
1364 /*       CALL MMARC41(NDIMAX,NDIMEN,NCOEFF,CURVE,UPARA0,UPARA1 */
1365 /*                  ,CURVE,IERCOD) */
1366 /*     is absolutely LEGAL. (RBD) */
1367
1368 /* > */
1369 /* ********************************************************************** 
1370 */
1371
1372 /*   Name of the routine */
1373
1374 /*   Auxiliary table of coefficients of (UPARA1-UPARA0)T+UPARA0  */
1375 /*   with power N=1 to NCOEFF-1. */
1376
1377
1378     /* Parameter adjustments */
1379     crvnew_dim1 = *ndimax;
1380     crvnew_offset = crvnew_dim1 + 1;
1381     crvnew -= crvnew_offset;
1382     crvold_dim1 = *ndimax;
1383     crvold_offset = crvold_dim1 + 1;
1384     crvold -= crvold_offset;
1385
1386     /* Function Body */
1387     *iercod = 0;
1388 /* ********************************************************************** 
1389 */
1390 /*                CASE WHEN PROCESSING CAN'T BE DONE */
1391 /* ********************************************************************** 
1392 */
1393     if (*ncoeff > 61 || *ncoeff < 1) {
1394         *iercod = 10;
1395         goto L9999;
1396     }
1397 /* ********************************************************************** 
1398 */
1399 /*                         IF NO CHANGES */
1400 /* ********************************************************************** 
1401 */
1402     if (*ndimen == *ndimax && *upara0 == 0. && *upara1 == 1.) {
1403         nboct = (*ndimax << 3) * *ncoeff;
1404         AdvApp2Var_SysBase::mcrfill_(&nboct,
1405                  &crvold[crvold_offset], 
1406                  &crvnew[crvnew_offset]);
1407         goto L9999;
1408     }
1409 /* ********************************************************************** 
1410 */
1411 /*                    INVERSION 3D : FAST PROCESSING */
1412 /* ********************************************************************** 
1413 */
1414     if (*upara0 == 1. && *upara1 == 0.) {
1415         if (*ndimen == 3 && *ndimax == 3 && *ncoeff <= 21) {
1416             mvcvinv_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset], 
1417                     iercod);
1418             goto L9999;
1419         }
1420 /* ******************************************************************
1421 **** */
1422 /*                    INVERSION 2D : FAST PROCESSING */
1423 /* ******************************************************************
1424 **** */
1425         if (*ndimen == 2 && *ndimax == 2 && *ncoeff <= 21) {
1426             mvcvin2_(ncoeff, &crvold[crvold_offset], &crvnew[crvnew_offset], 
1427                     iercod);
1428             goto L9999;
1429         }
1430     }
1431 /* ********************************************************************** 
1432 */
1433 /*                          GENERAL PROCESSING */
1434 /* ********************************************************************** 
1435 */
1436 /* -------------------------- Initializations --------------------------- 
1437 */
1438
1439     i__1 = *ndimen;
1440     for (nd = 1; nd <= i__1; ++nd) {
1441         crvnew[nd + crvnew_dim1] = crvold[nd + crvold_dim1];
1442 /* L100: */
1443     }
1444     if (*ncoeff == 1) {
1445         goto L9999;
1446     }
1447     tbaux[0] = *upara0;
1448     tbaux[1] = *upara1 - *upara0;
1449
1450 /* ----------------------- Calculation of coeff. of CRVNEW ------------------ 
1451 */
1452
1453     i__1 = *ncoeff - 1;
1454     for (ncf = 2; ncf <= i__1; ++ncf) {
1455
1456 /* ------------ Take into account NCF-th coeff. of CRVOLD --------
1457 ---- */
1458
1459         i__2 = ncf - 1;
1460         for (ncj = 1; ncj <= i__2; ++ncj) {
1461             bid = tbaux[ncj - 1];
1462             i__3 = *ndimen;
1463             for (nd = 1; nd <= i__3; ++nd) {
1464                 crvnew[nd + ncj * crvnew_dim1] += crvold[nd + ncf * 
1465                         crvold_dim1] * bid;
1466 /* L400: */
1467             }
1468 /* L300: */
1469         }
1470
1471         bid = tbaux[ncf - 1];
1472         i__2 = *ndimen;
1473         for (nd = 1; nd <= i__2; ++nd) {
1474             crvnew[nd + ncf * crvnew_dim1] = crvold[nd + ncf * crvold_dim1] * 
1475                     bid;
1476 /* L500: */
1477         }
1478
1479 /* --------- Calculate (NCF+1) coeff. of ((U1-U0)*t + U0)**(NCF) ---
1480 ---- */
1481
1482         bid = *upara1 - *upara0;
1483         tbaux[ncf] = tbaux[ncf - 1] * bid;
1484         for (ncj = ncf; ncj >= 2; --ncj) {
1485             tbaux[ncj - 1] = tbaux[ncj - 1] * *upara0 + tbaux[ncj - 2] * bid;
1486 /* L600: */
1487         }
1488         tbaux[0] *= *upara0;
1489
1490 /* L200: */
1491     }
1492
1493 /* -------------- Take into account the last coeff. of CRVOLD ----------- 
1494 */
1495
1496     i__1 = *ncoeff - 1;
1497     for (ncj = 1; ncj <= i__1; ++ncj) {
1498         bid = tbaux[ncj - 1];
1499         i__2 = *ndimen;
1500         for (nd = 1; nd <= i__2; ++nd) {
1501             crvnew[nd + ncj * crvnew_dim1] += crvold[nd + *ncoeff * 
1502                     crvold_dim1] * bid;
1503 /* L800: */
1504         }
1505 /* L700: */
1506     }
1507     i__1 = *ndimen;
1508     for (nd = 1; nd <= i__1; ++nd) {
1509         crvnew[nd + *ncoeff * crvnew_dim1] = crvold[nd + *ncoeff * 
1510                 crvold_dim1] * tbaux[*ncoeff - 1];
1511 /* L900: */
1512     }
1513
1514 /* ---------------------------- The end --------------------------------- 
1515 */
1516
1517 L9999:
1518     if (*iercod != 0) {
1519         AdvApp2Var_SysBase::maermsg_("MMARC41", iercod, 7L);
1520     }
1521
1522  return 0 ;
1523 } /* mmarc41_ */
1524
1525 //=======================================================================
1526 //function : AdvApp2Var_MathBase::mmarcin_
1527 //purpose  : 
1528 //=======================================================================
1529 int AdvApp2Var_MathBase::mmarcin_(integer *ndimax, 
1530                                   integer *ndim, 
1531                                   integer *ncoeff, 
1532                                   doublereal *crvold, 
1533                                   doublereal *u0, 
1534                                   doublereal *u1, 
1535                                   doublereal *crvnew, 
1536                                   integer *iercod)
1537
1538 {
1539   /* System generated locals */
1540   integer crvold_dim1, crvold_offset, crvnew_dim1, crvnew_offset, i__1, 
1541   i__2, i__3;
1542   doublereal d__1;
1543   
1544   /* Local variables */
1545   static doublereal x0, x1;
1546   static integer nd;
1547   static doublereal tabaux[61];
1548   static integer ibb;
1549   static doublereal bid;
1550   static integer ncf;
1551   static integer ncj;
1552   static doublereal eps3;
1553   
1554
1555
1556 /* ********************************************************************** 
1557 *//*     FUNCTION : */
1558 /*     ---------- */
1559 /*     Creation of curve C2(v) defined on [U0,U1] identic to */
1560 /*     curve C1(u) defined on [-1,1] (change of parameter */
1561 /*     of a curve) with INVERSION of indices of the resulting table. */
1562
1563 /*     KEYWORDS : */
1564 /*     ----------- */
1565 /*        GENERALIZED LIMITATION, RESTRICTION, INVERSION, CURVE */
1566
1567 /*     INPUT ARGUMENTS : */
1568 /*     ------------------ */
1569 /*   NDIMAX : Maximum Space Dimensioning. */
1570 /*   NDIMEN : Curve Dimension. */
1571 /*   NCOEFF : Nb of coefficients of the curve. */
1572 /*   CRVOLD : Curve to be limited. */
1573 /*   U0     : Min limit of the interval limiting the curve. 
1574 */
1575 /*   U1     : Max limit of the interval limiting the curve. 
1576 */
1577
1578 /*     OUTPUT ARGUMENTS : */
1579 /*     ------------------- */
1580 /*   CRVNEW : Relimited curve, defined on  [U0,U1] and equal to */
1581 /*            CRVOLD defined on [-1,1]. */
1582 /*   IERCOD : = 0, OK */
1583 /*            =10, Nb of coeff. <1 or > 61. */
1584 /*            =13, the requested interval of variation is null. */
1585 /*     COMMONS USED   : */
1586 /*     ---------------- */
1587 /*     REFERENCES CALLED   : */
1588 /*     ---------------------- */
1589 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1590 /*     ----------------------------------- */
1591 /* > */
1592 /* ********************************************************************** 
1593 */
1594
1595 /*   Name of the routine */
1596
1597 /*   Auxiliary table of coefficients of X1*T+X0 */
1598 /*   with power N=1 to NCOEFF-1. */
1599
1600
1601     /* Parameter adjustments */
1602     crvnew_dim1 = *ndimax;
1603     crvnew_offset = crvnew_dim1 + 1;
1604     crvnew -= crvnew_offset;
1605     crvold_dim1 = *ncoeff;
1606     crvold_offset = crvold_dim1 + 1;
1607     crvold -= crvold_offset;
1608
1609     /* Function Body */
1610     ibb = AdvApp2Var_SysBase::mnfndeb_();
1611     if (ibb >= 2) {
1612         AdvApp2Var_SysBase::mgenmsg_("MMARCIN", 7L);
1613     }
1614
1615 /* At zero machine it is tested if the output interval is not null */
1616
1617     AdvApp2Var_MathBase::mmveps3_(&eps3);
1618     if ((d__1 = *u1 - *u0, advapp_abs(d__1)) < eps3) {
1619         *iercod = 13;
1620         goto L9999;
1621     }
1622     *iercod = 0;
1623
1624 /* ********************************************************************** 
1625 */
1626 /*                CASE WHEN THE PROCESSING IS IMPOSSIBLE */
1627 /* ********************************************************************** 
1628 */
1629     if (*ncoeff > 61 || *ncoeff < 1) {
1630         *iercod = 10;
1631         goto L9999;
1632     }
1633 /* ********************************************************************** 
1634 */
1635 /*          IF NO CHANGE OF THE INTERVAL OF DEFINITION */
1636 /*          (ONLY INVERSION OF INDICES OF TABLE CRVOLD) */
1637 /* ********************************************************************** 
1638 */
1639     if (*ndim == *ndimax && *u0 == -1. && *u1 == 1.) {
1640         AdvApp2Var_MathBase::mmcvinv_(ndim, ncoeff, ndim, &crvold[crvold_offset], &crvnew[
1641                 crvnew_offset]);
1642         goto L9999;
1643     }
1644 /* ********************************************************************** 
1645 */
1646 /*          CASE WHEN THE NEW INTERVAL OF DEFINITION IS [0,1] */
1647 /* ********************************************************************** 
1648 */
1649     if (*u0 == 0. && *u1 == 1.) {
1650         mmcvstd_(ncoeff, ndimax, ncoeff, ndim, &crvold[crvold_offset], &
1651                 crvnew[crvnew_offset]);
1652         goto L9999;
1653     }
1654 /* ********************************************************************** 
1655 */
1656 /*                          GENERAL PROCESSING */
1657 /* ********************************************************************** 
1658 */
1659 /* -------------------------- Initialization --------------------------- 
1660 */
1661
1662     x0 = -(*u1 + *u0) / (*u1 - *u0);
1663     x1 = 2. / (*u1 - *u0);
1664     i__1 = *ndim;
1665     for (nd = 1; nd <= i__1; ++nd) {
1666         crvnew[nd + crvnew_dim1] = crvold[nd * crvold_dim1 + 1];
1667 /* L100: */
1668     }
1669     if (*ncoeff == 1) {
1670         goto L9999;
1671     }
1672     tabaux[0] = x0;
1673     tabaux[1] = x1;
1674
1675 /* ----------------------- Calculation of coeff. of CRVNEW ------------------ 
1676 */
1677
1678     i__1 = *ncoeff - 1;
1679     for (ncf = 2; ncf <= i__1; ++ncf) {
1680
1681 /* ------------ Take into account the NCF-th coeff. of CRVOLD --------
1682 ---- */
1683
1684         i__2 = ncf - 1;
1685         for (ncj = 1; ncj <= i__2; ++ncj) {
1686             bid = tabaux[ncj - 1];
1687             i__3 = *ndim;
1688             for (nd = 1; nd <= i__3; ++nd) {
1689                 crvnew[nd + ncj * crvnew_dim1] += crvold[ncf + nd * 
1690                         crvold_dim1] * bid;
1691 /* L400: */
1692             }
1693 /* L300: */
1694         }
1695
1696         bid = tabaux[ncf - 1];
1697         i__2 = *ndim;
1698         for (nd = 1; nd <= i__2; ++nd) {
1699             crvnew[nd + ncf * crvnew_dim1] = crvold[ncf + nd * crvold_dim1] * 
1700                     bid;
1701 /* L500: */
1702         }
1703
1704 /* --------- Calculation of (NCF+1) coeff. of [X1*t + X0]**(NCF) --------
1705 ---- */
1706
1707         tabaux[ncf] = tabaux[ncf - 1] * x1;
1708         for (ncj = ncf; ncj >= 2; --ncj) {
1709             tabaux[ncj - 1] = tabaux[ncj - 1] * x0 + tabaux[ncj - 2] * x1;
1710 /* L600: */
1711         }
1712         tabaux[0] *= x0;
1713
1714 /* L200: */
1715     }
1716
1717 /* -------------- Take into account the last coeff. of CRVOLD ----------- 
1718 */
1719
1720     i__1 = *ncoeff - 1;
1721     for (ncj = 1; ncj <= i__1; ++ncj) {
1722         bid = tabaux[ncj - 1];
1723         i__2 = *ndim;
1724         for (nd = 1; nd <= i__2; ++nd) {
1725             crvnew[nd + ncj * crvnew_dim1] += crvold[*ncoeff + nd * 
1726                     crvold_dim1] * bid;
1727 /* L800: */
1728         }
1729 /* L700: */
1730     }
1731     i__1 = *ndim;
1732     for (nd = 1; nd <= i__1; ++nd) {
1733         crvnew[nd + *ncoeff * crvnew_dim1] = crvold[*ncoeff + nd * 
1734                 crvold_dim1] * tabaux[*ncoeff - 1];
1735 /* L900: */
1736     }
1737
1738 /* ---------------------------- The end --------------------------------- 
1739 */
1740
1741 L9999:
1742     if (*iercod > 0) {
1743         AdvApp2Var_SysBase::maermsg_("MMARCIN", iercod, 7L);
1744     }
1745     if (ibb >= 2) {
1746         AdvApp2Var_SysBase::mgsomsg_("MMARCIN", 7L);
1747     }
1748     return 0;
1749 } /* mmarcin_ */
1750
1751 //=======================================================================
1752 //function : mmatvec_
1753 //purpose  : 
1754 //=======================================================================
1755 int mmatvec_(integer *nligne, 
1756              integer *,//ncolon,
1757              integer *gposit,
1758              integer *,//gnstoc, 
1759              doublereal *gmatri,
1760              doublereal *vecin, 
1761              integer *deblig,
1762              doublereal *vecout,
1763              integer *iercod)
1764
1765 {
1766   /* System generated locals */
1767   integer i__1, i__2;
1768   
1769   /* Local variables */
1770     static logical ldbg;
1771   static integer jmin, jmax, i__, j, k;
1772   static doublereal somme;
1773   static integer aux;
1774
1775
1776 /* ***********************************************************************
1777  */
1778
1779 /*     FUNCTION : */
1780 /*     ---------- */
1781 /*      Produce vector matrix in form of profile */
1782
1783
1784 /*     MOTS CLES : */
1785 /*     ----------- */
1786 /*      RESERVE, MATRIX, PRODUCT, VECTOR, PROFILE */
1787
1788 /*     INPUT ARGUMENTS : */
1789 /*     -------------------- */
1790 /*       NLIGNE : Line number of the matrix of constraints */
1791 /*       NCOLON : Number of column of the matrix of constraints */
1792 /*       GNSTOC: Number of coefficients in the profile of matrix GMATRI */
1793
1794 /*       GPOSIT: Table of positioning of terms of storage */
1795 /*               GPOSIT(1,I) contains the number of terms-1 on the line I 
1796 /*               in the profile of the matrix. */
1797 /*              GPOSIT(2,I) contains the index of storage of diagonal term*/
1798 /*               of line I */
1799 /*               GPOSIT(3,I) contains the index of column of the first term of */
1800 /*                           profile of line I */
1801 /*       GNSTOC: Number of coefficients in the profile of matrix */
1802 /*               GMATRI */
1803 /*       GMATRI : Matrix of constraints in form of profile */
1804 /*       VECIN  : Input vector */
1805 /*       DEBLIG : Line indexusing which the vector matrix is calculated */
1806 /*               
1807 /*     OUTPUT ARGUMENTS */
1808 /*     --------------------- */
1809 /*       VECOUT : VECTOR PRODUCT */
1810
1811 /*       IERCOD : ERROR CODE */
1812
1813
1814 /*     COMMONS USED : */
1815 /*     ------------------ */
1816
1817
1818 /*     REFERENCES CALLED : */
1819 /*     --------------------- */
1820
1821
1822 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1823 /*     ----------------------------------- */
1824
1825 /* ***********************************************************************
1826  */
1827 /*                            DECLARATIONS */
1828 /* ***********************************************************************
1829  */
1830
1831
1832
1833 /* ***********************************************************************
1834  */
1835 /*                      INITIALISATIONS */
1836 /* ***********************************************************************
1837  */
1838
1839     /* Parameter adjustments */
1840     --vecout;
1841     gposit -= 4;
1842     --vecin;
1843     --gmatri;
1844
1845     /* Function Body */
1846     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1847     if (ldbg) {
1848         AdvApp2Var_SysBase::mgenmsg_("MMATVEC", 7L);
1849     }
1850     *iercod = 0;
1851
1852 /* ***********************************************************************
1853  */
1854 /*                    Processing */
1855 /* ***********************************************************************
1856  */
1857     AdvApp2Var_SysBase::mvriraz_(nligne, 
1858              &vecout[1]);
1859     i__1 = *nligne;
1860     for (i__ = *deblig; i__ <= i__1; ++i__) {
1861         somme = 0.;
1862         jmin = gposit[i__ * 3 + 3];
1863         jmax = gposit[i__ * 3 + 1] + gposit[i__ * 3 + 3] - 1;
1864         aux = gposit[i__ * 3 + 2] - gposit[i__ * 3 + 1] - jmin + 1;
1865         i__2 = jmax;
1866         for (j = jmin; j <= i__2; ++j) {
1867             k = j + aux;
1868             somme += gmatri[k] * vecin[j];
1869         }
1870         vecout[i__] = somme;
1871     }
1872
1873
1874
1875
1876
1877     goto L9999;
1878
1879 /* ***********************************************************************
1880  */
1881 /*                   ERROR PROCESSING */
1882 /* ***********************************************************************
1883  */
1884
1885
1886
1887
1888 /* ***********************************************************************
1889  */
1890 /*                   RETURN CALLING PROGRAM */
1891 /* ***********************************************************************
1892  */
1893
1894 L9999:
1895
1896 /* ___ DESALLOCATION, ... */
1897
1898     AdvApp2Var_SysBase::maermsg_("MMATVEC", iercod, 7L);
1899     if (ldbg) {
1900         AdvApp2Var_SysBase::mgsomsg_("MMATVEC", 7L);
1901     }
1902
1903  return 0 ;
1904 } /* mmatvec_ */
1905
1906 //=======================================================================
1907 //function : mmbulld_
1908 //purpose  : 
1909 //=======================================================================
1910 int AdvApp2Var_MathBase::mmbulld_(integer *nbcoln, 
1911                                   integer *nblign, 
1912                                   doublereal *dtabtr, 
1913                                   integer *numcle)
1914
1915 {
1916   /* System generated locals */
1917   integer dtabtr_dim1, dtabtr_offset, i__1, i__2;
1918   
1919   /* Local variables */
1920   static logical ldbg;
1921   static doublereal daux;
1922   static integer nite1, nite2, nchan, i1, i2;
1923   
1924 /* ***********************************************************************
1925  */
1926
1927 /*     FUNCTION : */
1928 /*     ---------- */
1929 /*        Parsing of columns of a table of integers in increasing order */
1930 /*     KEYWORDS : */
1931 /*     ----------- */
1932 /*     POINT-ENTRY, PARSING */
1933 /*     INPUT ARGUMENTS : */
1934 /*     -------------------- */
1935 /*       - NBCOLN : Number of columns in the table */
1936 /*       - NBLIGN : Number of lines in the table */
1937 /*       - DTABTR : Table of integers to be parsed */
1938 /*       - NUMCLE : Position of the key on the column */
1939
1940 /*     OUTPUT ARGUMENTS : */
1941 /*     --------------------- */
1942 /*       - DTABTR : Parsed table */
1943
1944 /*     COMMONS USED : */
1945 /*     ------------------ */
1946
1947
1948 /*     REFERENCES CALLED : */
1949 /*     --------------------- */
1950
1951
1952 /*     DESCRIPTION/NOTES/LIMITATIONS : */
1953 /*     ----------------------------------- */
1954 /*     Particularly performant if the table is almost parsed */
1955 /*     In the opposite case it is better to use MVSHELD */
1956 /* ***********************************************************************
1957  */
1958
1959     /* Parameter adjustments */
1960     dtabtr_dim1 = *nblign;
1961     dtabtr_offset = dtabtr_dim1 + 1;
1962     dtabtr -= dtabtr_offset;
1963
1964     /* Function Body */
1965     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 2;
1966     if (ldbg) {
1967         AdvApp2Var_SysBase::mgenmsg_("MMBULLD", 7L);
1968     }
1969     nchan = 1;
1970     nite1 = *nbcoln;
1971     nite2 = 2;
1972
1973 /* ***********************************************************************
1974  */
1975 /*                     PROCESSING */
1976 /* ***********************************************************************
1977  */
1978
1979 /* ---->ALGORITHM in N^2 / 2 additional iteration */
1980
1981     while(nchan != 0) {
1982
1983 /* ----> Parsing from left to the right */
1984
1985         nchan = 0;
1986         i__1 = nite1;
1987         for (i1 = nite2; i1 <= i__1; ++i1) {
1988             if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 - 1)
1989                      * dtabtr_dim1]) {
1990                 i__2 = *nblign;
1991                 for (i2 = 1; i2 <= i__2; ++i2) {
1992                     daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
1993                     dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 * 
1994                             dtabtr_dim1];
1995                     dtabtr[i2 + i1 * dtabtr_dim1] = daux;
1996                 }
1997                 if (nchan == 0) {
1998                     nchan = 1;
1999                 }
2000             }
2001         }
2002         --nite1;
2003
2004 /* ----> Parsing from right to the left */
2005
2006         if (nchan != 0) {
2007             nchan = 0;
2008             i__1 = nite2;
2009             for (i1 = nite1; i1 >= i__1; --i1) {
2010                 if (dtabtr[*numcle + i1 * dtabtr_dim1] < dtabtr[*numcle + (i1 
2011                         - 1) * dtabtr_dim1]) {
2012                     i__2 = *nblign;
2013                     for (i2 = 1; i2 <= i__2; ++i2) {
2014                         daux = dtabtr[i2 + (i1 - 1) * dtabtr_dim1];
2015                         dtabtr[i2 + (i1 - 1) * dtabtr_dim1] = dtabtr[i2 + i1 *
2016                                  dtabtr_dim1];
2017                         dtabtr[i2 + i1 * dtabtr_dim1] = daux;
2018                     }
2019                     if (nchan == 0) {
2020                         nchan = 1;
2021                     }
2022                 }
2023             }
2024             ++nite2;
2025         }
2026     }
2027
2028
2029     goto L9999;
2030
2031 /* ***********************************************************************
2032  */
2033 /*                   ERROR PROCESSING */
2034 /* ***********************************************************************
2035  */
2036
2037 /* ----> No errors at calling functions, only tests and loops. */
2038
2039 /* ***********************************************************************
2040  */
2041 /*                   RETURN CALLING PROGRAM */
2042 /* ***********************************************************************
2043  */
2044
2045 L9999:
2046
2047     if (ldbg) {
2048         AdvApp2Var_SysBase::mgsomsg_("MMBULLD", 7L);
2049     }
2050
2051  return 0 ;
2052 } /* mmbulld_ */
2053
2054
2055 //=======================================================================
2056 //function : AdvApp2Var_MathBase::mmcdriv_
2057 //purpose  : 
2058 //=======================================================================
2059 int AdvApp2Var_MathBase::mmcdriv_(integer *ndimen, 
2060                                   integer *ncoeff, 
2061                                   doublereal *courbe, 
2062                                   integer *ideriv, 
2063                                   integer *ncofdv, 
2064                                   doublereal *crvdrv)
2065
2066
2067 {
2068   /* System generated locals */
2069   integer courbe_dim1, courbe_offset, crvdrv_dim1, crvdrv_offset, i__1, 
2070   i__2;
2071   
2072   /* Local variables */
2073   static integer i__, j, k;
2074   static doublereal mfactk, bid;
2075   
2076
2077 /* ***********************************************************************
2078  */
2079
2080 /*     FUNCTION : */
2081 /*     ---------- */
2082 /*     Calculate matrix of a derivate curve of order IDERIV. */
2083 /*     with input parameters other than output parameters. */
2084
2085
2086 /*     KEYWORDS : */
2087 /*     ----------- */
2088 /*     COEFFICIENTS,CURVE,DERIVATE I-EME. */
2089
2090 /*     INPUT ARGUMENTS : */
2091 /*     ------------------ */
2092 /*   NDIMEN  : Space dimension (2 or 3 in general) */
2093 /*   NCOEFF  : Degree +1 of the curve. */
2094 /*   COURBE  : Table of coefficients of the curve. */
2095 /*   IDERIV  : Required order of derivation : 1=1st derivate, etc... */
2096
2097 /*     OUTPUT ARGUMENTS : */
2098 /*     ------------------- */
2099 /*   NCOFDV  : Degree +1 of the derivative of order IDERIV of the curve. */
2100 /*   CRVDRV  : Table of coefficients of the derivative of order IDERIV */
2101 /*            of the curve. */
2102
2103 /*     COMMONS USED   : */
2104 /*     ---------------- */
2105
2106 /*     REFERENCES CALLED   : */
2107 /*     ----------------------- */
2108
2109 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2110 /*     ----------------------------------- */
2111
2112 /* ---> It is possible to take as output argument the curve */
2113 /*     and the number of coeff passed at input by making : */
2114 /*        CALL MMCDRIV(NDIMEN,NCOEFF,COURBE,IDERIV,NCOEFF,COURBE). */
2115 /*     After this call, NCOEFF does the number of coeff of the derived */
2116 /*     curve the coefficients which of are stored in CURVE. */
2117 /*     Attention to the coefficients of CURVE of rank superior to */
2118 /*     NCOEFF : they are not set to zero. */
2119
2120 /* ---> Algorithm : */
2121 /*     The code below was written basing on the following algorithm: 
2122 */
2123
2124 /*     Let P(t) = a1 + a2*t + ... an*t**n. Derivate of order k of P */
2125 /*     (containing n-k coefficients) is calculated as follows : */
2126
2127 /*       Pk(t) = a(k+1)*CNP(k,k)*k! */
2128 /*             + a(k+2)*CNP(k+1,k)*k! * t */
2129 /*             . */
2130 /*             . */
2131 /*             . */
2132 /*             + a(n)*CNP(n-1,k)*k! * t**(n-k-1). */
2133 /* ***********************************************************************
2134  */
2135
2136
2137 /* -------------- Case when the order of derivative is  ------------------- 
2138 */
2139 /* ---------------- greater than the degree of the curve --------------------- 
2140 */
2141
2142 /* ********************************************************************** 
2143 */
2144
2145 /*     FUNCTION : */
2146 /*     ---------- */
2147 /*      Serves to provide the coefficients of binome (Pascal's triangle). */
2148
2149 /*     KEYWORDS : */
2150 /*     ----------- */
2151 /*      Binomial coeff from 0 to 60. read only . init par block data */
2152
2153 /*     DEMSCRIPTION/NOTES/LIMITATIONS : */
2154 /*     ----------------------------------- */
2155 /*     Binomial coefficients form a triangular matrix. */
2156 /*     This matrix is completed in table CNP by its transposition. */
2157 /*     So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
2158
2159 /*     Initialization is done by block-data MMLLL09.RES, */
2160 /*     created by program MQINICNP.FOR). */
2161 /* ********************************************************************** 
2162 */
2163
2164
2165
2166 /* ***********************************************************************
2167  */
2168
2169     /* Parameter adjustments */
2170     crvdrv_dim1 = *ndimen;
2171     crvdrv_offset = crvdrv_dim1 + 1;
2172     crvdrv -= crvdrv_offset;
2173     courbe_dim1 = *ndimen;
2174     courbe_offset = courbe_dim1 + 1;
2175     courbe -= courbe_offset;
2176
2177     /* Function Body */
2178     if (*ideriv >= *ncoeff) {
2179         i__1 = *ndimen;
2180         for (i__ = 1; i__ <= i__1; ++i__) {
2181             crvdrv[i__ + crvdrv_dim1] = 0.;
2182 /* L10: */
2183         }
2184         *ncofdv = 1;
2185         goto L9999;
2186     }
2187 /* ********************************************************************** 
2188 */
2189 /*                        General processing */
2190 /* ********************************************************************** 
2191 */
2192 /* --------------------- Calculation of Factorial(IDERIV) ------------------ 
2193 */
2194
2195     k = *ideriv;
2196     mfactk = 1.;
2197     i__1 = k;
2198     for (i__ = 2; i__ <= i__1; ++i__) {
2199         mfactk *= i__;
2200 /* L50: */
2201     }
2202
2203 /* ------------ Calculation of coeff of the derived of order IDERIV ---------- 
2204 */
2205 /* ---> Attention :  coefficient binomial C(n,m) is represented in */
2206 /*                 MCCNP by CNP(N+1,M+1). */
2207
2208     i__1 = *ncoeff;
2209     for (j = k + 1; j <= i__1; ++j) {
2210         bid = mmcmcnp_.cnp[j - 1 + k * 61] * mfactk;
2211         i__2 = *ndimen;
2212         for (i__ = 1; i__ <= i__2; ++i__) {
2213             crvdrv[i__ + (j - k) * crvdrv_dim1] = bid * courbe[i__ + j * 
2214                     courbe_dim1];
2215 /* L200: */
2216         }
2217 /* L100: */
2218     }
2219
2220     *ncofdv = *ncoeff - *ideriv;
2221
2222 /* -------------------------------- The end ----------------------------- 
2223 */
2224
2225 L9999:
2226     return 0;
2227 } /* mmcdriv_ */
2228
2229 //=======================================================================
2230 //function : AdvApp2Var_MathBase::mmcglc1_
2231 //purpose  : 
2232 //=======================================================================
2233 int AdvApp2Var_MathBase::mmcglc1_(integer *ndimax, 
2234                                   integer *ndimen, 
2235                                   integer *ncoeff, 
2236                                   doublereal *courbe, 
2237                                   doublereal *tdebut, 
2238                                   doublereal *tfinal, 
2239                                   doublereal *epsiln, 
2240                                   doublereal *xlongc, 
2241                                   doublereal *erreur, 
2242                                   integer *iercod)
2243
2244
2245 {
2246   /* System generated locals */
2247   integer courbe_dim1, courbe_offset, i__1;
2248   doublereal d__1;
2249   
2250   /* Local variables */
2251   static integer ndec;
2252   static doublereal tdeb, tfin;
2253   static integer iter;
2254   static doublereal oldso;
2255   static integer itmax;
2256   static doublereal sottc;
2257   static integer kk, ibb;
2258   static doublereal dif, pas;
2259   static doublereal som;
2260  
2261
2262 /* ***********************************************************************
2263  */
2264
2265 /*     FUNCTION : */
2266 /*     ---------- */
2267 /*      Allows calculating the length of an arc of curve POLYNOMIAL */
2268 /*      on an interval [A,B]. */
2269
2270 /*     KEYWORDS : */
2271 /*     ----------- */
2272 /*        LENGTH,CURVE,GAUSS,PRIVATE. */
2273
2274 /*     INPUT ARGUMENTS : */
2275 /*     ------------------ */
2276 /*      NDIMAX : Max. number of lines of tables */
2277 /*               (i.e. max. nb of polynoms). */
2278 /*      NDIMEN : Dimension of the space (nb of polynoms). */
2279 /*      NCOEFF : Nb of coefficients of the polynom. This is degree + 1. 
2280 */
2281 /*      COURBE(NDIMAX,NCOEFF) : Coefficients of the curve. */
2282 /*      TDEBUT : Lower limit of the interval of integration for  */
2283 /*               length calculation. */
2284 /*      TFINAL : Upper limit of the interval of integration for */
2285 /*               length calculation. */
2286 /*      EPSILN : REQIRED precision for length calculation. */
2287
2288 /*     OUTPUT ARGUMENTS : */
2289 /*     ------------------- */
2290 /*      XLONGC : Length of the arc of curve */
2291 /*      ERREUR : Precision OBTAINED for the length calculation. */
2292 /*      IERCOD : Error code, 0 OK, >0 Serious error. */
2293 /*               = 1 Too much iterations, the best calculated resultat */
2294 /*                   (is almost ERROR) */
2295 /*               = 2 Pb MMLONCV (no result) */
2296 /*               = 3 NDIM or NCOEFF invalid (no result) */
2297
2298 /*     COMMONS USED : */
2299 /*     ---------------- */
2300
2301 /*     REFERENCES CALLED : */
2302 /*     ----------------------- */
2303
2304 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2305 /*     ----------------------------------- */
2306 /*      The polynom is actually a set of polynoms with */
2307 /*      coefficients arranged in a table of 2 indices, */
2308 /*      each line relative to the polynom. */
2309 /*      The polynom is defined by these coefficients ordered */
2310 /*      by increasing power of the variable. */
2311 /*      All polynoms have the same number of coefficients (the */
2312 /*      same degree). */
2313
2314 /*      This program cancels and replaces LENGCV, MLONGC and MLENCV. */
2315
2316 /*      ATTENTION : if TDEBUT > TFINAL, the length is NEGATIVE. */
2317
2318 /* > */
2319 /* ***********************************************************************
2320  */
2321
2322 /*   Name of the routine */
2323
2324
2325 /* ------------------------ General Initialization --------------------- 
2326 */
2327
2328     /* Parameter adjustments */
2329     courbe_dim1 = *ndimax;
2330     courbe_offset = courbe_dim1 + 1;
2331     courbe -= courbe_offset;
2332
2333     /* Function Body */
2334     ibb = AdvApp2Var_SysBase::mnfndeb_();
2335     if (ibb >= 2) {
2336         AdvApp2Var_SysBase::mgenmsg_("MMCGLC1", 7L);
2337     }
2338
2339     *iercod = 0;
2340     *xlongc = 0.;
2341     *erreur = 0.;
2342
2343 /* ------ Test of equity of limits */
2344
2345     if (*tdebut == *tfinal) {
2346         *iercod = 0;
2347         goto L9999;
2348     }
2349
2350 /* ------ Test of the dimension and the number of coefficients */
2351
2352     if (*ndimen <= 0 || *ncoeff <= 0) {
2353         goto L9003;
2354     }
2355
2356 /* ----- Nb of current cutting, nb of iteration, */
2357 /*       max nb of iterations */
2358
2359     ndec = 1;
2360     iter = 1;
2361
2362     itmax = 13;
2363
2364 /* ------ Variation of the nb of intervals */
2365 /*       Multiplied by 2 at each iteration */
2366
2367 L5000:
2368     pas = (*tfinal - *tdebut) / ndec;
2369     sottc = 0.;
2370
2371 /* ------ Loop on all current NDEC intervals */
2372
2373     i__1 = ndec;
2374     for (kk = 1; kk <= i__1; ++kk) {
2375
2376 /* ------ Limits of the current integration interval */
2377
2378         tdeb = *tdebut + (kk - 1) * pas;
2379         tfin = tdeb + pas;
2380         mmloncv_(ndimax, ndimen, ncoeff, &courbe[courbe_offset], &tdeb, &tfin,
2381                  &som, iercod);
2382         if (*iercod > 0) {
2383             goto L9002;
2384         }
2385
2386         sottc += som;
2387
2388 /* L100: */
2389     }
2390
2391
2392 /* ----------------- Test of the maximum number of iterations ------------ 
2393 */
2394
2395 /*  Test if passes at least once ** */
2396
2397     if (iter == 1) {
2398         oldso = sottc;
2399         ndec <<= 1;
2400         ++iter;
2401         goto L5000;
2402     } else {
2403
2404 /* ------ Take into account DIF - Test of convergence */
2405
2406         ++iter;
2407         dif = (d__1 = sottc - oldso, advapp_abs(d__1));
2408
2409 /* ------ If DIF is OK, leave..., otherwise: */
2410
2411         if (dif > *epsiln) {
2412
2413 /* ------ If nb iteration exceeded, leave */
2414
2415             if (iter > itmax) {
2416                 *iercod = 1;
2417                 goto L9000;
2418             } else {
2419
2420 /* ------ Otherwise continue by cutting the initial interval.
2421  */
2422
2423                 oldso = sottc;
2424                 ndec <<= 1;
2425                 goto L5000;
2426             }
2427         }
2428     }
2429
2430 /* ------------------------------ THE END ------------------------------- 
2431 */
2432
2433 L9000:
2434     *xlongc = sottc;
2435     *erreur = dif;
2436     goto L9999;
2437
2438 /* ---> PB in MMLONCV */
2439
2440 L9002:
2441     *iercod = 2;
2442     goto L9999;
2443
2444 /* ---> NCOEFF or NDIM invalid. */
2445
2446 L9003:
2447     *iercod = 3;
2448     goto L9999;
2449
2450 L9999:
2451     if (*iercod > 0) {
2452         AdvApp2Var_SysBase::maermsg_("MMCGLC1", iercod, 7L);
2453     }
2454     if (ibb >= 2) {
2455         AdvApp2Var_SysBase::mgsomsg_("MMCGLC1", 7L);
2456     }
2457     return 0;
2458 } /* mmcglc1_ */
2459
2460 //=======================================================================
2461 //function : mmchole_
2462 //purpose  : 
2463 //=======================================================================
2464 int mmchole_(integer *,//mxcoef, 
2465              integer *dimens, 
2466              doublereal *amatri, 
2467              integer *aposit, 
2468              integer *posuiv, 
2469              doublereal *chomat, 
2470              integer *iercod)
2471
2472 {
2473   /* System generated locals */
2474   integer i__1, i__2, i__3;
2475   doublereal d__1;
2476   
2477   /* Builtin functions */
2478   //double sqrt();
2479   
2480     /* Local variables */
2481   static logical ldbg;
2482   static integer kmin, i__, j, k;
2483   static doublereal somme;
2484   static integer ptini, ptcou;
2485
2486
2487 /* ***********************************************************************
2488  */
2489
2490 /*     FUNCTION : */
2491 /*     ----------                                                  T */
2492 /*     Produce decomposition of choleski of matrix A in S.S */
2493 /*     Calculate inferior triangular matrix S. */
2494
2495 /*     KEYWORDS : */
2496 /*     ----------- */
2497 /*     RESOLUTION, MFACTORISATION, MATRIX_PROFILE, CHOLESKI */
2498
2499 /*     INPUT ARGUMENTS : */
2500 /*     -------------------- */
2501 /*     MXCOEF : Max number of terms in the hessian profile */
2502 /*     DIMENS : Dimension of the problem */
2503 /*     AMATRI(MXCOEF) : Coefficients of the matrix profile */
2504 /*        APOSIT(1,*) : Distance diagonal-left extremity of the line 
2505 */
2506 /*        APOSIT(2,*) : Position of diagonal terms in HESSIE */
2507 /*     POSUIV(MXCOEF) :  first line inferior not out of profile */
2508
2509 /*     OUTPUT ARGUMENTS : */
2510 /*     --------------------- */
2511 /*      CHOMAT(MXCOEF) : Inferior triangular matrix preserving the */
2512 /*                       profile of AMATRI. */
2513 /*      IERCOD : error code */
2514 /*               = 0 : ok */
2515 /*               = 1 : non-defined positive matrix */
2516
2517 /*     COMMONS USED : */
2518 /*     ------------------ */
2519
2520 /*      .Neant. */
2521
2522 /*     REFERENCES CALLED   : */
2523 /*     ---------------------- */
2524
2525 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2526 /*     ----------------------------------- */
2527 /*     DEBUG LEVEL = 4 */
2528 /* ***********************************************************************
2529  */
2530 /*                            DECLARATIONS */
2531 /* ***********************************************************************
2532  */
2533
2534
2535
2536 /* ***********************************************************************
2537  */
2538 /*                      INITIALISATIONS */
2539 /* ***********************************************************************
2540  */
2541
2542     /* Parameter adjustments */
2543     --chomat;
2544     --posuiv;
2545     --amatri;
2546     aposit -= 3;
2547
2548     /* Function Body */
2549     ldbg = AdvApp2Var_SysBase::mnfndeb_() >= 4;
2550     if (ldbg) {
2551         AdvApp2Var_SysBase::mgenmsg_("MMCHOLE", 7L);
2552     }
2553     *iercod = 0;
2554
2555 /* ***********************************************************************
2556  */
2557 /*                    PROCESSING */
2558 /* ***********************************************************************
2559  */
2560
2561     i__1 = *dimens;
2562     for (j = 1; j <= i__1; ++j) {
2563
2564         ptini = aposit[(j << 1) + 2];
2565
2566         somme = 0.;
2567         i__2 = ptini - 1;
2568         for (k = ptini - aposit[(j << 1) + 1]; k <= i__2; ++k) {
2569 /* Computing 2nd power */
2570             d__1 = chomat[k];
2571             somme += d__1 * d__1;
2572         }
2573
2574         if (amatri[ptini] - somme < 1e-32) {
2575             goto L9101;
2576         }
2577         chomat[ptini] = sqrt(amatri[ptini] - somme);
2578
2579         ptcou = ptini;
2580
2581         while(posuiv[ptcou] > 0) {
2582
2583             i__ = posuiv[ptcou];
2584             ptcou = aposit[(i__ << 1) + 2] - (i__ - j);
2585
2586 /*           Calculate the sum of S  .S   for k =1 a j-1 */
2587 /*                               ik  jk */
2588             somme = 0.;
2589 /* Computing MAX */
2590             i__2 = i__ - aposit[(i__ << 1) + 1], i__3 = j - aposit[(j << 1) + 
2591                     1];
2592             kmin = advapp_max(i__2,i__3);
2593             i__2 = j - 1;
2594             for (k = kmin; k <= i__2; ++k) {
2595                 somme += chomat[aposit[(i__ << 1) + 2] - (i__ - k)] * chomat[
2596                         aposit[(j << 1) + 2] - (j - k)];
2597             }
2598
2599             chomat[ptcou] = (amatri[ptcou] - somme) / chomat[ptini];
2600         }
2601     }
2602
2603     goto L9999;
2604
2605 /* ***********************************************************************
2606  */
2607 /*                   ERROR PROCESSING */
2608 /* ***********************************************************************
2609  */
2610
2611 L9101:
2612     *iercod = 1;
2613     goto L9999;
2614
2615 /* ***********************************************************************
2616  */
2617 /*                  RETURN CALLING PROGRAM */
2618 /* ***********************************************************************
2619  */
2620
2621 L9999:
2622
2623     AdvApp2Var_SysBase::maermsg_("MMCHOLE", iercod, 7L);
2624     if (ldbg) {
2625         AdvApp2Var_SysBase::mgsomsg_("MMCHOLE", 7L);
2626     }
2627
2628  return 0 ;
2629 } /* mmchole_ */
2630
2631 //=======================================================================
2632 //function : AdvApp2Var_MathBase::mmcvctx_
2633 //purpose  : 
2634 //=======================================================================
2635 int AdvApp2Var_MathBase::mmcvctx_(integer *ndimen, 
2636                                   integer *ncofmx, 
2637                                   integer *nderiv, 
2638                                   doublereal *ctrtes, 
2639                                   doublereal *crvres, 
2640                                   doublereal *tabaux, 
2641                                   doublereal *xmatri, 
2642                                   integer *iercod)
2643
2644 {
2645   /* System generated locals */
2646   integer ctrtes_dim1, ctrtes_offset, crvres_dim1, crvres_offset, 
2647   xmatri_dim1, xmatri_offset, tabaux_dim1, tabaux_offset, i__1, 
2648   i__2;
2649   
2650   /* Local variables */
2651   static integer moup1, nordr;
2652   static integer nd;
2653   static integer ibb, ncf, ndv;
2654   static doublereal eps1;
2655
2656
2657 /* ***********************************************************************
2658  */
2659
2660 /*     FUNCTION : */
2661 /*     ---------- */
2662 /*        Calculate a polynomial curve checking the  */
2663 /*        passage constraints (interpolation) */
2664 /*        from first derivatives, etc... to extremities. */
2665 /*        Parameters at the extremities are supposed to be -1 and 1. */
2666
2667 /*     KEYWORDS : */
2668 /*     ----------- */
2669 /*     ALL, AB_SPECIFI::CONSTRAINTS&,INTERPOLATION,&CURVE */
2670
2671 /*     INPUT ARGUMENTS : */
2672 /*     ------------------ */
2673 /*     NDIMEN : Space Dimension. */
2674 /*     NCOFMX : Nb of coeff. of curve CRVRES on each */
2675 /*              dimension. */
2676 /*     NDERIV : Order of constraint with derivatives : */
2677 /*              0 --> interpolation simple. */
2678 /*              1 --> interpolation+constraints with 1st. */
2679 /*              2 --> cas (0)+ (1) +   "         "   2nd derivatives. */
2680 /*                 etc... */
2681 /*     CTRTES : Table of constraints. */
2682 /*              CTRTES(*,1,*) = contraints at -1. */
2683 /*              CTRTES(*,2,*) = contraints at  1. */
2684
2685 /*     OUTPUT ARGUMENTS : */
2686 /*     ------------------- */
2687 /*     CRVRES : Resulting curve defined on (-1,1). */
2688 /*     TABAUX : Auxilliary matrix. */
2689 /*     XMATRI : Auxilliary matrix. */
2690
2691 /*     COMMONS UTILISES   : */
2692 /*     ---------------- */
2693
2694 /*      .Neant. */
2695
2696 /*     REFERENCES CALLED   : */
2697 /*     ---------------------- */
2698 /*     Type  Name */
2699 /*           MAERMSG         R*8  DFLOAT              MGENMSG */
2700 /*           MGSOMSG              MMEPS1               MMRSLW */
2701 /*      I*4  MNFNDEB */
2702
2703 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2704 /*     ----------------------------------- */
2705 /*        The polynom (or the curve) is calculated by solving a */
2706 /*        system of linear equations. If the imposed degree is great */
2707 /*        it is preferable to call a routine based on */
2708 /*        Lagrange or Hermite interpolation depending on the case. */
2709 /*        (for a high degree the matrix of the system can be badly */
2710 /*        conditionned). */
2711 /*        This routine returns a curve defined in (-1,1). */
2712 /*        In general case, it is necessary to use MCVCTG. */
2713 /* > */
2714 /* ***********************************************************************
2715  */
2716
2717 /*   Name of the routine */
2718
2719
2720     /* Parameter adjustments */
2721     crvres_dim1 = *ncofmx;
2722     crvres_offset = crvres_dim1 + 1;
2723     crvres -= crvres_offset;
2724     xmatri_dim1 = *nderiv + 1;
2725     xmatri_offset = xmatri_dim1 + 1;
2726     xmatri -= xmatri_offset;
2727     tabaux_dim1 = *nderiv + 1 + *ndimen;
2728     tabaux_offset = tabaux_dim1 + 1;
2729     tabaux -= tabaux_offset;
2730     ctrtes_dim1 = *ndimen;
2731     ctrtes_offset = ctrtes_dim1 * 3 + 1;
2732     ctrtes -= ctrtes_offset;
2733
2734     /* Function Body */
2735     ibb = AdvApp2Var_SysBase::mnfndeb_();
2736     if (ibb >= 3) {
2737         AdvApp2Var_SysBase::mgenmsg_("MMCVCTX", 7L);
2738     }
2739 /*   Precision. */
2740     AdvApp2Var_MathBase::mmeps1_(&eps1);
2741
2742 /* ****************** CALCULATION OF EVEN COEFFICIENTS ********************* 
2743 */
2744 /* ------------------------- Initialization ----------------------------- 
2745 */
2746
2747     nordr = *nderiv + 1;
2748     i__1 = nordr;
2749     for (ncf = 1; ncf <= i__1; ++ncf) {
2750         tabaux[ncf + tabaux_dim1] = 1.;
2751 /* L100: */
2752     }
2753
2754 /* ---------------- Calculation of terms corresponding to derivatives ------- 
2755 */
2756
2757     i__1 = nordr;
2758     for (ndv = 2; ndv <= i__1; ++ndv) {
2759         i__2 = nordr;
2760         for (ncf = 1; ncf <= i__2; ++ncf) {
2761             tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) * 
2762                     tabaux_dim1] * (doublereal) ((ncf << 1) - ndv);
2763 /* L300: */
2764         }
2765 /* L200: */
2766     }
2767
2768 /* ------------------ Writing the second member ----------------------- 
2769 */
2770
2771     moup1 = 1;
2772     i__1 = nordr;
2773     for (ndv = 1; ndv <= i__1; ++ndv) {
2774         i__2 = *ndimen;
2775         for (nd = 1; nd <= i__2; ++nd) {
2776             tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1) 
2777                     + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2778                      * ctrtes_dim1]) / 2.;
2779 /* L500: */
2780         }
2781         moup1 = -moup1;
2782 /* L400: */
2783     }
2784
2785 /* -------------------- Resolution of the system --------------------------- 
2786 */
2787
2788     mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2789             xmatri_offset], iercod);
2790     if (*iercod > 0) {
2791         goto L9999;
2792     }
2793     i__1 = *ndimen;
2794     for (nd = 1; nd <= i__1; ++nd) {
2795         i__2 = nordr;
2796         for (ncf = 1; ncf <= i__2; ++ncf) {
2797             crvres[(ncf << 1) - 1 + nd * crvres_dim1] = xmatri[ncf + nd * 
2798                     xmatri_dim1];
2799 /* L700: */
2800         }
2801 /* L600: */
2802     }
2803
2804 /* ***************** CALCULATION OF UNEVEN COEFFICIENTS ******************** 
2805 */
2806 /* ------------------------- Initialization ----------------------------- 
2807 */
2808
2809
2810     i__1 = nordr;
2811     for (ncf = 1; ncf <= i__1; ++ncf) {
2812         tabaux[ncf + tabaux_dim1] = 1.;
2813 /* L1100: */
2814     }
2815
2816 /* ---------------- Calculation of terms corresponding to derivatives ------- 
2817 */
2818
2819     i__1 = nordr;
2820     for (ndv = 2; ndv <= i__1; ++ndv) {
2821         i__2 = nordr;
2822         for (ncf = 1; ncf <= i__2; ++ncf) {
2823             tabaux[ncf + ndv * tabaux_dim1] = tabaux[ncf + (ndv - 1) * 
2824                     tabaux_dim1] * (doublereal) ((ncf << 1) - ndv + 1);
2825 /* L1300: */
2826         }
2827 /* L1200: */
2828     }
2829
2830 /* ------------------ Writing of the second member ----------------------- 
2831 */
2832
2833     moup1 = -1;
2834     i__1 = nordr;
2835     for (ndv = 1; ndv <= i__1; ++ndv) {
2836         i__2 = *ndimen;
2837         for (nd = 1; nd <= i__2; ++nd) {
2838             tabaux[nordr + nd + ndv * tabaux_dim1] = (ctrtes[nd + ((ndv << 1) 
2839                     + 2) * ctrtes_dim1] + moup1 * ctrtes[nd + ((ndv << 1) + 1)
2840                      * ctrtes_dim1]) / 2.;
2841 /* L1500: */
2842         }
2843         moup1 = -moup1;
2844 /* L1400: */
2845     }
2846
2847 /* -------------------- Solution of the system --------------------------- 
2848 */
2849
2850     mmrslw_(&nordr, &nordr, ndimen, &eps1, &tabaux[tabaux_offset], &xmatri[
2851             xmatri_offset], iercod);
2852     if (*iercod > 0) {
2853         goto L9999;
2854     }
2855     i__1 = *ndimen;
2856     for (nd = 1; nd <= i__1; ++nd) {
2857         i__2 = nordr;
2858         for (ncf = 1; ncf <= i__2; ++ncf) {
2859             crvres[(ncf << 1) + nd * crvres_dim1] = xmatri[ncf + nd * 
2860                     xmatri_dim1];
2861 /* L1700: */
2862         }
2863 /* L1600: */
2864     }
2865
2866 /* --------------------------- The end ---------------------------------- 
2867 */
2868
2869 L9999:
2870     if (*iercod != 0) {
2871         AdvApp2Var_SysBase::maermsg_("MMCVCTX", iercod, 7L);
2872     }
2873     if (ibb >= 3) {
2874         AdvApp2Var_SysBase::mgsomsg_("MMCVCTX", 7L);
2875     }
2876
2877  return 0 ;
2878 } /* mmcvctx_ */
2879
2880 //=======================================================================
2881 //function : AdvApp2Var_MathBase::mmcvinv_
2882 //purpose  : 
2883 //=======================================================================
2884  int AdvApp2Var_MathBase::mmcvinv_(integer *ndimax, 
2885                             integer *ncoef,
2886                             integer *ndim, 
2887                             doublereal *curveo, 
2888                             doublereal *curve)
2889
2890 {
2891   /* Initialized data */
2892   
2893   static char nomprg[8+1] = "MMCVINV ";
2894   
2895   /* System generated locals */
2896   integer curve_dim1, curve_offset, curveo_dim1, curveo_offset, i__1, i__2;
2897   
2898   /* Local variables */
2899   static integer i__, nd, ibb;
2900   
2901
2902 /* ***********************************************************************
2903  */
2904
2905 /*     FUNCTION : */
2906 /*     ---------- */
2907 /*        Inversion of arguments of the final curve. */
2908
2909 /*     KEYWORDS : */
2910 /*     ----------- */
2911 /*        SMOOTHING,CURVE */
2912
2913
2914 /*     INPUT ARGUMENTS : */
2915 /*     ------------------ */
2916
2917 /*        NDIM: Space Dimension. */
2918 /*        NCOEF: Degree of the polynom. */
2919 /*        CURVEO: The curve before inversion. */
2920
2921 /*     OUTPUT ARGUMENTS : */
2922 /*     ------------------- */
2923 /*        CURVE: The curve after inversion. */
2924
2925 /*     COMMONS USED : */
2926 /*     ---------------- */
2927 /*     REFERENCES APPELEES   : */
2928 /*     ----------------------- */
2929 /*     DESCRIPTION/NOTES/LIMITATIONS : */
2930 /*     ----------------------------------- */
2931 /* ***********************************************************************
2932  */
2933
2934 /*   The name of the routine */
2935     /* Parameter adjustments */
2936     curve_dim1 = *ndimax;
2937     curve_offset = curve_dim1 + 1;
2938     curve -= curve_offset;
2939     curveo_dim1 = *ncoef;
2940     curveo_offset = curveo_dim1 + 1;
2941     curveo -= curveo_offset;
2942
2943     /* Function Body */
2944
2945     ibb = AdvApp2Var_SysBase::mnfndeb_();
2946     if (ibb >= 2) {
2947         AdvApp2Var_SysBase::mgenmsg_(nomprg, 6L);
2948     }
2949
2950     i__1 = *ncoef;
2951     for (i__ = 1; i__ <= i__1; ++i__) {
2952         i__2 = *ndim;
2953         for (nd = 1; nd <= i__2; ++nd) {
2954             curve[nd + i__ * curve_dim1] = curveo[i__ + nd * curveo_dim1];
2955 /* L300: */
2956         }
2957     }
2958
2959 /* L9999: */
2960     return 0;
2961 } /* mmcvinv_ */
2962
2963 //=======================================================================
2964 //function : mmcvstd_
2965 //purpose  : 
2966 //=======================================================================
2967 int mmcvstd_(integer *ncofmx, 
2968              integer *ndimax, 
2969              integer *ncoeff,
2970              integer *ndimen, 
2971              doublereal *crvcan, 
2972              doublereal *courbe)
2973
2974 {
2975   /* System generated locals */
2976   integer courbe_dim1, crvcan_dim1, crvcan_offset, i__1, i__2, i__3;
2977   
2978   /* Local variables */
2979   static integer ndeg, i__, j, j1, nd, ibb;
2980   static doublereal bid;
2981   
2982
2983 /* ***********************************************************************
2984  */
2985
2986 /*     FUNCTION : */
2987 /*     ---------- */
2988 /*        Transform curve defined between [-1,1] into [0,1]. */
2989
2990 /*     KEYWORDS : */
2991 /*     ----------- */
2992 /*        LIMITATION,RESTRICTION,CURVE */
2993
2994 /*     INPUT ARGUMENTS : */
2995 /*     ------------------ */
2996 /*        NDIMAX : Dimension of the space. */
2997 /*        NDIMEN : Dimension of the curve. */
2998 /*        NCOEFF : Degree of the curve. */
2999 /*        CRVCAN(NCOFMX,NDIMEN): The curve is defined at the interval [-1,1]. */
3000
3001 /*     OUTPUT ARGUMENTS : */
3002 /*     ------------------- */
3003 /*        CURVE(NDIMAX,NCOEFF): Curve defined at the interval [0,1]. */
3004
3005 /*     COMMONS USED   : */
3006 /*     ---------------- */
3007
3008 /*     REFERENCES CALLED   : */
3009 /*     ----------------------- */
3010
3011 /*     DESCRIPTION/NOTES/LIMITATIONS : */
3012 /*     ----------------------------------- */
3013 /* > */
3014 /* ***********************************************************************
3015  */
3016
3017 /*   Name of the program. */
3018
3019
3020 /* ********************************************************************** 
3021 */
3022
3023 /*     FUNCTION : */
3024 /*     ---------- */
3025 /*      Provides binomial coefficients (Pascal triangle). */
3026
3027 /*     KEYWORDS : */
3028 /*     ----------- */
3029 /*      Binomial coefficient from 0 to 60. read only . init by block data */
3030
3031 /*     DEMSCRIPTION/NOTES/LIMITATIONS : */
3032 /*     ----------------------------------- */
3033 /*     Binomial coefficients form a triangular matrix. */
3034 /*     This matrix is completed in table CNP by its transposition. */
3035 /*     So: CNP(I,J) = CNP(J,I) for I and J = 0, ..., 60. */
3036
3037 /*     Initialization is done with block-data MMLLL09.RES, */
3038 /*     created by the program MQINICNP.FOR. */
3039 /* > */
3040 /* ********************************************************************** 
3041 */
3042
3043
3044
3045 /* ***********************************************************************
3046  */
3047
3048     /* Parameter adjustments */
3049     courbe_dim1 = *ndimax;
3050     --courbe;
3051     crvcan_dim1 = *ncofmx;
3052     crvcan_offset = crvcan_dim1;
3053     crvcan -= crvcan_offset;
3054
3055     /* Function Body */
3056     ibb = AdvApp2Var_SysBase::mnfndeb_();
3057     if (ibb >= 3) {
3058         AdvApp2Var_SysBase::mgenmsg_("MMCVSTD", 7L);
3059     }
3060     ndeg = *ncoeff - 1;
3061
3062 /* ------------------ Construction of the resulting curve ---------------- 
3063 */
3064
3065     i__1 = *ndimen;
3066     for (nd = 1; nd <= i__1; ++nd) {
3067         i__2 = ndeg;
3068         for (j = 0; j <= i__2; ++j) {
3069             bid = 0.;
3070             i__3 = ndeg;
3071             for (i__ = j; i__ <= i__3; i__ += 2) {
3072                 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j 
3073                         * 61];
3074 /* L410: */
3075             }
3076             courbe[nd + j * courbe_dim1] = bid;
3077
3078             bid = 0.;
3079             j1 = j + 1;
3080             i__3 = ndeg;
3081             for (i__ = j1; i__ <= i__3; i__ += 2) {
3082                 bid += crvcan[i__ + nd * crvcan_dim1] * mmcmcnp_.cnp[i__ + j 
3083                         * 61];
3084 /* L420: */
3085             }
3086             courbe[nd + j * courbe_dim1] -= bid;
3087 /* L400: */
3088         }
3089 /* L300: */
3090     }
3091
3092 /* ------------------- Renormalization of the CURVE -------------------------
3093  */
3094
3095     bid = 1.;
3096     i__1 = ndeg;
3097     for (i__ = 0; i__ <= i__1; ++i__) {
3098         i__2 = *ndimen;
3099         for (nd = 1; nd <= i__2; ++nd) {
3100             courbe[nd + i__ * courbe_dim1] *= bid;
3101 /* L510: */
3102         }
3103         bid *= 2.;
3104 /* L500: */
3105     }
3106
3107 /* ----------------------------- The end -------------------------------- 
3108 */
3109
3110     if (ibb >= 3) {
3111         AdvApp2Var_SysBase::mgsomsg_("MMCVSTD", 7L);
3112     }
3113     return 0;
3114 } /* mmcvstd_ */
3115
3116 //=======================================================================
3117 //function : AdvApp2Var_MathBase::mmdrc11_
3118 //purpose  : 
3119 //=======================================================================
3120 int AdvApp2Var_MathBase::mmdrc11_(integer *iordre, 
3121                                   integer *ndimen, 
3122                                   integer *ncoeff, 
3123                                   doublereal *courbe, 
3124                                   doublereal *points, 
3125                                   doublereal *mfactab)
3126
3127 {
3128   /* System generated locals */
3129   integer courbe_dim1, courbe_offset, points_dim2, points_offset, i__1, 
3130   i__2;
3131   
3132   /* Local variables */
3133   
3134   static integer ndeg, i__, j, ndgcb, nd, ibb;
3135   
3136
3137 /* ********************************************************************** 
3138 */
3139
3140 /*     FUNCTION : */
3141 /*     ---------- */
3142 /*        Calculation of successive derivatives of equation CURVE with */
3143 /*        parameters -1, 1 from order 0 to order IORDRE */
3144 /*        included. The calculation is produced without knowing the coefficients of */
3145 /*        derivatives of the curve. */
3146
3147 /*     KEYWORDS : */
3148 /*     ----------- */
3149 /*        POSITIONING,EXTREMITIES,CURVE,DERIVATIVE. */
3150
3151 /*     INPUT ARGUMENTS : */
3152 /*     ------------------ */
3153 /*        IORDRE  : Maximum order of calculation of derivatives. */
3154 /*        NDIMEN  : Dimension of the space. */
3155 /*        NCOEFF  : Number of coefficients of the curve (degree+1). */
3156 /*        COURBE  : Table of coefficients of the curve. */
3157
3158 /*     OUTPUT ARGUMENTS : */
3159 /*     ------------------- */
3160 /*        POINTS    : Table of values of consecutive derivatives */
3161 /*                 of parameters -1.D0 and 1.D0. */
3162 /*        MFACTAB : Auxiliary table for calculation of factorial(I). 
3163 */
3164
3165 /*     COMMONS USED   : */
3166 /*     ---------------- */
3167 /*        None. */
3168
3169 /*     REFERENCES CALLED   : */
3170 /*     ----------------------- */
3171
3172 /*     DESCRIPTION/NOTES/LIMITATIONS : */
3173 /*     ----------------------------------- */
3174
3175 /* ---> ATTENTION, the coefficients of the curve are  */
3176 /*     in a reverse order. */
3177
3178 /* ---> The algorithm of calculation of derivatives is based on */
3179 /*     generalization of Horner scheme : */
3180 /*                          k             2 */
3181 /*          Let C(t) = uk.t  + ... + u2.t  + u1.t + u0 . */
3182
3183
3184 /*      a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3185
3186 /*          aj = a(j-1).x + u(k-j) */
3187 /*          bj = b(j-1).x + a(j-1) */
3188 /*          cj = c(j-1).x + b(j-1) */
3189
3190 /*     So : C(x) = ak, C'(x) = bk, C"(x) = 2.ck  . */
3191
3192 /*     The algorithm is generalized easily for calculation of */
3193
3194 /*               (n) */
3195 /*              C  (x)   . */
3196 /*             --------- */
3197 /*                n! */
3198
3199 /*      Reference : D. KNUTH, "The Art of Computer Programming" */
3200 /*      ---------              Vol. 2/Seminumerical Algorithms */
3201 /*                             Addison-Wesley Pub. Co. (1969) */
3202 /*                             pages 423-425. */
3203 /* > */
3204 /* ********************************************************************** 
3205 */
3206
3207 /*   Name of the routine */
3208
3209     /* Parameter adjustments */
3210     points_dim2 = *iordre + 1;
3211     points_offset = (points_dim2 << 1) + 1;
3212     points -= points_offset;
3213     courbe_dim1 = *ncoeff;
3214     courbe_offset = courbe_dim1;
3215     courbe -= courbe_offset;
3216
3217     /* Function Body */
3218     ibb = AdvApp2Var_SysBase::mnfndeb_();
3219     if (ibb >= 2) {
3220         AdvApp2Var_SysBase::mgenmsg_("MMDRC11", 7L);
3221     }
3222
3223     if (*iordre < 0 || *ncoeff < 1) {
3224         goto L9999;
3225     }
3226
3227 /* ------------------- Initialization of table POINTS ----------------- 
3228 */
3229
3230     ndgcb = *ncoeff - 1;
3231     i__1 = *ndimen;
3232     for (nd = 1; nd <= i__1; ++nd) {
3233         points[(nd * points_dim2 << 1) + 1] = courbe[ndgcb + nd * courbe_dim1]
3234                 ;
3235         points[(nd * points_dim2 << 1) + 2] = courbe[ndgcb + nd * courbe_dim1]
3236                 ;
3237 /* L100: */
3238     }
3239
3240     i__1 = *ndimen;
3241     for (nd = 1; nd <= i__1; ++nd) {
3242         i__2 = *iordre;
3243         for (j = 1; j <= i__2; ++j) {
3244             points[((j + nd * points_dim2) << 1) + 1] = 0.;
3245             points[((j + nd * points_dim2) << 1) + 2] = 0.;
3246 /* L400: */
3247         }
3248 /* L300: */
3249     }
3250
3251 /*    Calculation with parameter -1 and 1 */
3252
3253     i__1 = *ndimen;
3254     for (nd = 1; nd <= i__1; ++nd) {
3255         i__2 = ndgcb;
3256         for (ndeg = 1; ndeg <= i__2; ++ndeg) {
3257             for (i__ = *iordre; i__ >= 1; --i__) {
3258                 points[((i__ + nd * points_dim2) << 1) + 1] = -points[((i__ + nd 
3259                         * points_dim2) << 1) + 1] + points[((i__ - 1 + nd * 
3260                         points_dim2) << 1) + 1];
3261                 points[((i__ + nd * points_dim2) << 1) + 2] += points[((i__ - 1 
3262                         + nd * points_dim2) << 1) + 2];
3263 /* L800: */
3264             }
3265             points[(nd * points_dim2 << 1) + 1] = -points[(nd * points_dim2 <<
3266                      1) + 1] + courbe[ndgcb - ndeg + nd * courbe_dim1];
3267             points[(nd * points_dim2 << 1) + 2] += courbe[ndgcb - ndeg + nd * 
3268                     courbe_dim1];
3269 /* L700: */
3270         }
3271 /* L600: */
3272     }
3273
3274 /* --------------------- Multiplication by factorial(I) -------------- 
3275 */
3276
3277     if (*iordre > 1) {
3278         mfac_(&mfactab[1], iordre);
3279
3280         i__1 = *ndimen;
3281         for (nd = 1; nd <= i__1; ++nd) {
3282             i__2 = *iordre;
3283             for (i__ = 2; i__ <= i__2; ++i__) {
3284                 points[((i__ + nd * points_dim2) << 1) + 1] = mfactab[i__] * 
3285                         points[((i__ + nd * points_dim2) << 1) + 1];
3286                 points[((i__ + nd * points_dim2) << 1) + 2] = mfactab[i__] * 
3287                         points[((i__ + nd * points_dim2) << 1) + 2];
3288 /* L1000: */
3289             }
3290 /* L900: */
3291         }
3292     }
3293
3294 /* ---------------------------- End ------------------------------------- 
3295 */
3296
3297 L9999:
3298     if (ibb >= 2) {
3299         AdvApp2Var_SysBase::mgsomsg_("MMDRC11", 7L);
3300     }
3301     return 0;
3302 } /* mmdrc11_ */
3303
3304 //=======================================================================
3305 //function : mmdrvcb_
3306 //purpose  : 
3307 //=======================================================================
3308 int mmdrvcb_(integer *ideriv,
3309              integer *ndim, 
3310              integer *ncoeff,
3311              doublereal *courbe, 
3312              doublereal *tparam,
3313              doublereal *tabpnt, 
3314              integer *iercod)
3315
3316 {
3317   /* System generated locals */
3318   integer courbe_dim1, tabpnt_dim1, i__1, i__2, i__3;
3319   
3320   /* Local variables */
3321   static integer ndeg, i__, j, nd, ndgcrb, iptpnt, ibb;
3322   
3323
3324 /* ***********************************************************************
3325 /*     FUNCTION : */
3326 /*     ---------- */
3327
3328 /*        Calculation of successive derivatives of equation CURVE with */
3329 /*        parameter TPARAM from order 0 to order IDERIV included. */
3330 /*        The calculation is produced without knowing the coefficients of */
3331 /*        derivatives of the CURVE. */
3332
3333 /*     KEYWORDS : */
3334 /*     ----------- */
3335 /*        POSITIONING,PARAMETER,CURVE,DERIVATIVE. */
3336
3337 /*     INPUT ARGUMENTS : */
3338 /*     ------------------ */
3339 /*        IORDRE  : Maximum order of calculation of derivatives. */
3340 /*        NDIMEN  : Dimension of the space. */
3341 /*        NCOEFF  : Number of coefficients of the curve (degree+1). */
3342 /*        COURBE  : Table of coefficients of the curve. */
3343 /*        TPARAM  : Value of the parameter where the curve should be evaluated. */
3344
3345 /*     OUTPUT ARGUMENTS : */
3346 /*     ------------------- */
3347 /*        TABPNT  : Table of values of consecutive derivatives */
3348 /*                  of parameter TPARAM. */
3349   /*        IERCOD  : 0 = OK, */
3350 /*                    1 = incoherent input. */
3351
3352 /*     COMMONS USED  : */
3353 /*     ---------------- */
3354 /*        None. */
3355
3356 /*     REFERENCES CALLED   : */
3357 /*     ----------------------- */
3358
3359 /*     DESCRIPTION/NOTES/LIMITATIONS : */
3360 /*     ----------------------------------- */
3361
3362 /*     The algorithm of  calculation of derivatives is based on */
3363 /*     generalization of the Horner scheme : */
3364 /*                          k             2 */
3365 /*          Let C(t) = uk.t  + ... + u2.t  + u1.t + u0 . */
3366
3367
3368 /*      a0 = uk, b0 = 0, c0 = 0 and for 1<=j<=k, it is calculated : */
3369
3370 /*          aj = a(j-1).x + u(k-j) */
3371 /*          bj = b(j-1).x + a(j-1) */
3372 /*          cj = c(j-1).x + b(j-1) */
3373
3374 /*     So, it is obtained : C(x) = ak, C'(x) = bk, C"(x) = 2.ck  . */
3375
3376 /*     The algorithm can be easily generalized for the calculation of */
3377
3378 /*               (n) */
3379 /*              C  (x)   . */
3380 /*             --------- */
3381 /*                n! */
3382
3383 /*      Reference : D. KNUTH, "The Art of Computer Programming" */
3384 /*      ---------              Vol. 2/Seminumerical Algorithms */
3385 /*                             Addison-Wesley Pub. Co. (1969) */
3386 /*                             pages 423-425. */
3387
3388 /* ---> To evaluare derivatives at 0 and 1, it is preferable */
3389 /*      to use routine MDRV01.FOR . */
3390 /* > */
3391 /* ********************************************************************** 
3392 */
3393
3394 /*   Name of the routine */
3395
3396     /* Parameter adjustments */
3397     tabpnt_dim1 = *ndim;
3398     --tabpnt;
3399     courbe_dim1 = *ndim;
3400     --courbe;
3401
3402     /* Function Body */
3403     ibb = AdvApp2Var_SysBase::mnfndeb_();
3404     if (ibb >= 2) {
3405         AdvApp2Var_SysBase::mgenmsg_("MMDRVCB", 7L);
3406     }
3407
3408     if (*ideriv < 0 || *ncoeff < 1) {
3409         *iercod = 1;
3410         goto L9999;
3411     }
3412     *iercod = 0;
3413
3414 /* ------------------- Initialization of table TABPNT ----------------- 
3415 */
3416
3417     ndgcrb = *ncoeff - 1;
3418     i__1 = *ndim;
3419     for (nd = 1; nd <= i__1; ++nd) {
3420         tabpnt[nd] = courbe[nd + ndgcrb * courbe_dim1];
3421 /* L100: */
3422     }
3423
3424     if (*ideriv < 1) {
3425         goto L200;
3426     }
3427     iptpnt = *ndim * *ideriv;
3428     AdvApp2Var_SysBase::mvriraz_(&iptpnt, 
3429              &tabpnt[tabpnt_dim1 + 1]);
3430 L200:
3431
3432 /* ------------------------ Calculation of parameter TPARAM ------------------ 
3433 */
3434
3435     i__1 = ndgcrb;
3436     for (ndeg = 1; ndeg <= i__1; ++ndeg) {
3437         i__2 = *ndim;
3438         for (nd = 1; nd <= i__2; ++nd) {
3439             for (i__ = *ideriv; i__ >= 1; --i__) {
3440                 tabpnt[nd + i__ * tabpnt_dim1] = tabpnt[nd + i__ * 
3441                         tabpnt_dim1] * *tparam + tabpnt[nd + (i__ - 1) * 
3442                         tabpnt_dim1];
3443 /* L700: */
3444             }
3445             tabpnt[nd] = tabpnt[nd] * *tparam + courbe[nd + (ndgcrb - ndeg) * 
3446                     courbe_dim1];
3447 /* L600: */
3448         }
3449 /* L500: */
3450     }
3451
3452 /* --------------------- Multiplication by factorial(I) ------------- 
3453 */
3454
3455     i__1 = *ideriv;
3456     for (i__ = 2; i__ <= i__1; ++i__) {
3457         i__2 = i__;
3458         for (j = 2; j <= i__2; ++j) {
3459             i__3 = *ndim;
3460             for (nd = 1; nd <= i__3; ++nd) {
3461                 tabpnt[nd + i__ * tabpnt_dim1] = (doublereal) j * tabpnt[nd + 
3462                         i__ * tabpnt_dim1];
3463 /* L1200: */
3464             }
3465 /* L1100: */
3466         }
3467 /* L1000: */
3468     }
3469
3470 /* --------------------------- The end --------------------------------- 
3471 */
3472
3473 L9999:
3474     if (*iercod > 0) {
3475         AdvApp2Var_SysBase::maermsg_("MMDRVCB", iercod, 7L);
3476     }
3477     return 0;
3478 } /* mmdrvcb_ */
3479
3480 //=======================================================================
3481 //function : AdvApp2Var_MathBase::mmdrvck_
3482 //purpose  : 
3483 //=======================================================================
3484 int AdvApp2Var_MathBase::mmdrvck_(integer *ncoeff, 
3485                                   integer *ndimen, 
3486                                   doublereal *courbe, 
3487                                   integer *ideriv, 
3488                                   doublereal *tparam, 
3489                                   doublereal *pntcrb)
3490
3491 {
3492   /* Initialized data */
3493   
3494   static doublereal mmfack[21] = { 1.,2.,6.,24.,120.,720.,5040.,40320.,
3495             362880.,3628800.,39916800.,479001600.,6227020800.,87178291200.,
3496             1.307674368e12,2.0922789888e13,3.55687428096e14,6.402373705728e15,
3497             1.21645100408832e17,2.43290200817664e18,5.109094217170944e19 };
3498   
3499   /* System generated locals */
3500   integer courbe_dim1, courbe_offset, i__1, i__2;
3501   
3502   /* Local variables */
3503   static integer i__, j, k, nd;
3504   static doublereal mfactk, bid;
3505