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