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