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