0022550: Fixing data races
[occt.git] / src / Approx / Approx_CurvlinFunc.cxx
1 // File:        Approx_CurvlinFunc.cxx
2 // Created:     Tue May 12 14:03:09 1998
3 // Author:      Roman BORISOV
4 //              <rbv@sgi38>
5
6
7 #include <Approx_CurvlinFunc.ixx>
8 #include <Adaptor3d_CurveOnSurface.hxx>
9 #include <Adaptor3d_HCurveOnSurface.hxx>
10 #include <TColStd_SequenceOfReal.hxx>
11 #include <GeomLib.hxx>
12  #include <GCPnts_AbscissaPoint.hxx>
13 #include <Precision.hxx>
14
15 #ifdef __OCC_DEBUG_CHRONO
16 #include <OSD_Timer.hxx>
17 static OSD_Chronometer chr_uparam;
18 Standard_EXPORT Standard_Integer uparam_count;
19 Standard_EXPORT Standard_Real t_uparam;
20
21 //Standard_IMPORT extern void InitChron(OSD_Chronometer& ch);
22 Standard_IMPORT void InitChron(OSD_Chronometer& ch);
23 //Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
24 Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
25 #endif
26
27 static Standard_Real cubic(const Standard_Real X, const Standard_Real *Xi, const Standard_Real *Yi)
28 {
29   Standard_Real I1, I2, I3, I21, I22, I31, Result;
30
31   I1 = (Yi[0] - Yi[1])/(Xi[0] - Xi[1]);
32   I2 = (Yi[1] - Yi[2])/(Xi[1] - Xi[2]);
33   I3 = (Yi[2] - Yi[3])/(Xi[2] - Xi[3]);
34
35   I21 = (I1 - I2)/(Xi[0] - Xi[2]);
36   I22 = (I2 - I3)/(Xi[1] - Xi[3]);
37   
38   I31 = (I21 - I22)/(Xi[0] - Xi[3]);  
39
40   Result = Yi[0] + (X - Xi[0])*(I1 + (X - Xi[1])*(I21 + (X - Xi[2])*I31));
41
42   return Result;
43 }
44
45 //static void findfourpoints(const Standard_Real S, 
46 static void findfourpoints(const Standard_Real , 
47                            Standard_Integer NInterval, 
48                            const Handle(TColStd_HArray1OfReal)& Si, 
49                            Handle(TColStd_HArray1OfReal)& Ui, 
50                            const Standard_Real prevS, 
51                            const Standard_Real prevU, Standard_Real *Xi, 
52                            Standard_Real *Yi)
53 {
54   Standard_Integer i, j;
55   Standard_Integer NbInt = Si->Length() - 1;
56   if (NbInt < 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter");
57
58   if(NInterval < 1) NInterval = 1;
59   else if(NInterval > NbInt - 2) NInterval = NbInt - 2;
60
61   for(i = 0; i < 4; i++) {
62     Xi[i] = Si->Value(NInterval - 1 + i);
63     Yi[i] = Ui->Value(NInterval - 1 + i);
64   }
65 // try to insert (S, U)  
66   for(i = 0; i < 3; i++) {
67     if(Xi[i] < prevS && prevS < Xi[i+1]) {
68       for(j = 0; j < i; j++) {
69         Xi[j] = Xi[j+1];
70         Yi[j] = Yi[j+1];
71       }
72       Xi[i] = prevS;
73       Yi[i] = prevU;
74       break;
75     }
76   }
77 }
78
79 /*static Standard_Real curvature(const Standard_Real U, const Adaptor3d_Curve& C)
80 {
81   Standard_Real k, tau, mod1, mod2, OMEGA;
82   gp_Pnt P;
83   gp_Vec D1, D2, D3;
84   C.D3(U, P, D1, D2, D3);
85   mod1 = D1.Magnitude();
86   mod2 = D1.Crossed(D2).Magnitude();
87   k = mod2/(mod1*mod1*mod1);
88   tau = D1.Dot(D2.Crossed(D3));
89   tau /= mod2*mod2;
90   OMEGA = Sqrt(k*k + tau*tau);
91
92   return OMEGA;
93 }
94 */
95
96 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C), 
97                     myCase(1), 
98                     myFirstS(0), 
99                     myLastS(1), 
100                     myTolLen(Tol),
101                     myPrevS (0.0),
102                     myPrevU (0.0)
103 {
104   Init();
105 }
106
107 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D, const Handle(Adaptor3d_HSurface)& S, const Standard_Real Tol) : 
108                     myC2D1(C2D), 
109                     mySurf1(S), 
110                     myCase(2), 
111                     myFirstS(0), 
112                     myLastS(1), 
113                     myTolLen(Tol),
114                     myPrevS (0.0),
115                     myPrevU (0.0)
116 {  
117   Init();
118 }
119
120 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D1, const Handle(Adaptor2d_HCurve2d)& C2D2, const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_HSurface)& S2, const Standard_Real Tol) : 
121                     myC2D1(C2D1), 
122                     myC2D2(C2D2), 
123                     mySurf1(S1), 
124                     mySurf2(S2), 
125                     myCase(3), 
126                     myFirstS(0), 
127                     myLastS(1), 
128                     myTolLen(Tol),
129                     myPrevS (0.0),
130                     myPrevU (0.0)
131 {  
132   Init();
133 }
134
135 void Approx_CurvlinFunc::Init()
136 {
137   Adaptor3d_CurveOnSurface CurOnSur;
138   
139   switch(myCase) {
140   case 1:
141     Init(myC3D->GetCurve(), mySi_1, myUi_1);
142     myFirstU1 = myC3D->FirstParameter();
143     myLastU1 = myC3D->LastParameter();
144     myFirstU2 = myLastU2 = 0;
145     break;
146   case 2:
147     CurOnSur.Load(myC2D1);
148     CurOnSur.Load(mySurf1);
149     Init(CurOnSur, mySi_1, myUi_1);
150     myFirstU1 = CurOnSur.FirstParameter();
151     myLastU1 = CurOnSur.LastParameter();
152     myFirstU2 = myLastU2 = 0;
153     break;
154   case 3:
155     CurOnSur.Load(myC2D1);
156     CurOnSur.Load(mySurf1);
157     Init(CurOnSur, mySi_1, myUi_1);
158     myFirstU1 = CurOnSur.FirstParameter();
159     myLastU1 = CurOnSur.LastParameter();
160     CurOnSur.Load(myC2D2);
161     CurOnSur.Load(mySurf2);
162     Init(CurOnSur, mySi_2, myUi_2);
163     myFirstU2 = CurOnSur.FirstParameter();
164     myLastU2 = CurOnSur.LastParameter();
165   }
166
167   Length();
168 }
169
170
171 //=======================================================================
172 //function : Init
173 //purpose  : Init the values
174 //history  : 23/10/1998 PMN : Cut at curve's discontinuities
175 //=======================================================================
176 void Approx_CurvlinFunc::Init(Adaptor3d_Curve& C, Handle(TColStd_HArray1OfReal)& Si, 
177                               Handle(TColStd_HArray1OfReal)& Ui) const
178 {
179   Standard_Real Step, FirstU, LastU;
180   Standard_Integer i, j, k, NbInt, NbIntC3;
181   FirstU = C.FirstParameter();
182   LastU  = C.LastParameter();
183
184   NbInt = 10;
185   NbIntC3 = C.NbIntervals(GeomAbs_C3);
186   TColStd_Array1OfReal Disc(1, NbIntC3+1);
187
188   if (NbIntC3 >1) {
189      C.Intervals(Disc, GeomAbs_C3);
190   }
191   else {
192     Disc(1) = FirstU;
193     Disc(2) = LastU;
194   }
195
196   Ui = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
197   Si = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
198
199   Ui->SetValue(0, FirstU);
200   Si->SetValue(0, 0);
201
202   for(j = 1, i=1; j<=NbIntC3; j++) {
203     Step = (Disc(j+1) - Disc(j))/NbInt;
204     for(k = 1; k <= NbInt; k++, i++) {
205       Ui->ChangeValue(i) = Ui->Value(i-1) + Step;
206       Si->ChangeValue(i) = Si->Value(i-1) + Length(C, Ui->Value(i-1), Ui->Value(i));
207     }
208   }
209
210   Standard_Real Len = Si->Value(Si->Upper());
211   for(i = Si->Lower(); i<= Si->Upper(); i++)
212     Si->ChangeValue(i) /= Len;
213
214   // TODO - fields should be mutable
215   const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
216   const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
217 }
218
219 void  Approx_CurvlinFunc::SetTol(const Standard_Real Tol)
220 {
221   myTolLen = Tol;
222 }
223
224 Standard_Real Approx_CurvlinFunc::FirstParameter() const
225 {
226   return myFirstS;
227 }
228
229 Standard_Real Approx_CurvlinFunc::LastParameter() const
230 {
231   return myLastS;
232 }
233
234 Standard_Integer Approx_CurvlinFunc::NbIntervals(const GeomAbs_Shape S) const
235 {
236   Adaptor3d_CurveOnSurface CurOnSur;
237
238   switch(myCase) {
239   case 1:
240     return myC3D->NbIntervals(S);
241   case 2:
242     CurOnSur.Load(myC2D1);
243     CurOnSur.Load(mySurf1);
244     return CurOnSur.NbIntervals(S);
245   case 3:
246     Standard_Integer NbInt;
247     CurOnSur.Load(myC2D1);
248     CurOnSur.Load(mySurf1);
249     NbInt = CurOnSur.NbIntervals(S);
250     TColStd_Array1OfReal T1(1, NbInt+1);
251     CurOnSur.Intervals(T1, S);
252     CurOnSur.Load(myC2D2);
253     CurOnSur.Load(mySurf2);
254     NbInt = CurOnSur.NbIntervals(S);
255     TColStd_Array1OfReal T2(1, NbInt+1);
256     CurOnSur.Intervals(T2, S);
257
258     TColStd_SequenceOfReal Fusion;
259     GeomLib::FuseIntervals(T1, T2, Fusion);
260     return Fusion.Length() - 1;
261   }
262
263   //POP pour WNT
264   return 1;
265 }
266
267 void Approx_CurvlinFunc::Intervals(TColStd_Array1OfReal& T, const GeomAbs_Shape S) const
268 {
269   Adaptor3d_CurveOnSurface CurOnSur;
270   Standard_Integer i;
271     
272   switch(myCase) {
273   case 1:
274     myC3D->Intervals(T, S);
275     break;
276   case 2:
277     CurOnSur.Load(myC2D1);
278     CurOnSur.Load(mySurf1);
279     CurOnSur.Intervals(T, S);
280     break;
281   case 3:
282     Standard_Integer NbInt;
283     CurOnSur.Load(myC2D1);
284     CurOnSur.Load(mySurf1);
285     NbInt = CurOnSur.NbIntervals(S);
286     TColStd_Array1OfReal T1(1, NbInt+1);
287     CurOnSur.Intervals(T1, S);
288     CurOnSur.Load(myC2D2);
289     CurOnSur.Load(mySurf2);
290     NbInt = CurOnSur.NbIntervals(S);
291     TColStd_Array1OfReal T2(1, NbInt+1);
292     CurOnSur.Intervals(T2, S);
293
294     TColStd_SequenceOfReal Fusion;
295     GeomLib::FuseIntervals(T1, T2, Fusion);
296
297     for (i = 1; i <= Fusion.Length(); i++)
298       T.ChangeValue(i) = Fusion.Value(i);
299   }
300
301   for(i = 1; i <= T.Length(); i++)
302     T.ChangeValue(i) = GetSParameter(T.Value(i));
303 }
304
305 void Approx_CurvlinFunc::Trim(const Standard_Real First, const Standard_Real Last, const Standard_Real Tol)
306 {
307   if (First < 0 || Last >1) Standard_OutOfRange::Raise("Approx_CurvlinFunc::Trim");
308   if ((Last - First) < Tol) return;
309
310   Standard_Real FirstU, LastU;
311   Adaptor3d_CurveOnSurface CurOnSur;
312   Handle(Adaptor3d_HCurve) HCurOnSur;
313
314   switch(myCase) {
315   case 1:
316     myC3D = myC3D->Trim(myFirstU1, myLastU1, Tol);
317     FirstU = GetUParameter(myC3D->GetCurve(), First, 1);
318     LastU = GetUParameter(myC3D->GetCurve(), Last, 1);
319     myC3D = myC3D->Trim(FirstU, LastU, Tol);
320     break;
321   case 3:
322     CurOnSur.Load(myC2D2);
323     CurOnSur.Load(mySurf2);
324     HCurOnSur = CurOnSur.Trim(myFirstU2, myLastU2, Tol);
325     myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
326     mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
327     CurOnSur.Load(myC2D2);
328     CurOnSur.Load(mySurf2);
329
330     FirstU = GetUParameter(CurOnSur, First, 1);
331     LastU = GetUParameter(CurOnSur, Last, 1);
332     HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
333     myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
334     mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
335
336   case 2:
337     CurOnSur.Load(myC2D1);
338     CurOnSur.Load(mySurf1);
339     HCurOnSur = CurOnSur.Trim(myFirstU1, myLastU1, Tol);
340     myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
341     mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
342     CurOnSur.Load(myC2D1);
343     CurOnSur.Load(mySurf1);
344
345     FirstU = GetUParameter(CurOnSur, First, 1);
346     LastU = GetUParameter(CurOnSur, Last, 1);
347     HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
348     myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
349     mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
350   }
351   myFirstS = First;
352   myLastS = Last;  
353 }
354
355 void Approx_CurvlinFunc::Length()
356 {
357   Adaptor3d_CurveOnSurface CurOnSur;
358   Standard_Real FirstU, LastU;
359
360   switch(myCase){
361   case 1:
362     FirstU = myC3D->FirstParameter();
363     LastU = myC3D->LastParameter();
364     myLength = Length(myC3D->GetCurve(), FirstU, LastU);    
365     myLength1 = myLength2 = 0;    
366     break;
367   case 2:
368     CurOnSur.Load(myC2D1);
369     CurOnSur.Load(mySurf1);
370     FirstU = CurOnSur.FirstParameter();
371     LastU = CurOnSur.LastParameter();
372     myLength = Length(CurOnSur, FirstU, LastU);
373     myLength1 = myLength2 = 0;    
374     break;
375   case 3:
376     CurOnSur.Load(myC2D1);
377     CurOnSur.Load(mySurf1);
378     FirstU = CurOnSur.FirstParameter();
379     LastU = CurOnSur.LastParameter();
380     myLength1 = Length(CurOnSur, FirstU, LastU);
381     CurOnSur.Load(myC2D2);
382     CurOnSur.Load(mySurf2);
383     FirstU = CurOnSur.FirstParameter();
384     LastU = CurOnSur.LastParameter();
385     myLength2 = Length(CurOnSur, FirstU, LastU);
386     myLength = (myLength1 + myLength2)/2;
387   }
388 }
389
390
391 Standard_Real Approx_CurvlinFunc::Length(Adaptor3d_Curve& C, const Standard_Real FirstU, const Standard_Real LastU) const
392 {
393   Standard_Real Length;
394
395   Length = GCPnts_AbscissaPoint::Length(C, FirstU, LastU, myTolLen);
396   return Length;
397 }
398
399
400 Standard_Real Approx_CurvlinFunc::GetLength() const
401 {
402   return myLength;
403 }
404
405 Standard_Real Approx_CurvlinFunc::GetSParameter(const Standard_Real U) const
406 {
407   Standard_Real S=0, S1, S2;
408   Adaptor3d_CurveOnSurface CurOnSur;
409
410   switch (myCase) {
411   case 1:
412     S = GetSParameter(myC3D->GetCurve(), U, myLength);
413     break;
414   case 2:
415     CurOnSur.Load(myC2D1);
416     CurOnSur.Load(mySurf1);
417     S = GetSParameter(CurOnSur, U, myLength);
418     break;
419   case 3:
420     CurOnSur.Load(myC2D1);
421     CurOnSur.Load(mySurf1);
422     S1 = GetSParameter(CurOnSur, U, myLength1);    
423     CurOnSur.Load(myC2D2);
424     CurOnSur.Load(mySurf2);
425     S2 = GetSParameter(CurOnSur, U, myLength2);    
426     S = (S1 + S2)/2;
427   }
428   return S;
429 }
430
431
432
433 Standard_Real Approx_CurvlinFunc::GetUParameter(Adaptor3d_Curve& C, 
434                                                 const Standard_Real S, 
435                                                 const Standard_Integer NumberOfCurve) const
436 {
437   Standard_Real deltaS, base, U, Length;
438   Standard_Integer NbInt, NInterval, i;
439   Handle(TColStd_HArray1OfReal) InitUArray, InitSArray;
440 #ifdef __OCC_DEBUG_CHRONO
441   InitChron(chr_uparam);
442 #endif
443   if(S < 0 || S > 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter");
444
445   if(NumberOfCurve == 1) {
446     InitUArray = myUi_1;
447     InitSArray = mySi_1;
448     if(myCase == 3)
449       Length = myLength1;
450     else 
451       Length = myLength;
452   }
453   else {
454     InitUArray = myUi_2;
455     InitSArray = mySi_2;
456     Length = myLength2;
457   }
458
459   NbInt = InitUArray->Length() - 1;
460
461   if(S == 1) NInterval = NbInt - 1;
462   else {
463     for(i = 0; i < NbInt; i++) {
464       if((InitSArray->Value(i) <= S && S < InitSArray->Value(i+1)))
465         break;
466     }
467     NInterval = i;
468   }
469   if(S==InitSArray->Value(NInterval)) {
470     return InitUArray->Value(NInterval);
471   }
472   if(S==InitSArray->Value(NInterval+1)) {
473     return InitUArray->Value(NInterval+1);
474   } 
475
476   base = InitUArray->Value(NInterval);
477   deltaS = (S - InitSArray->Value(NInterval))*Length;
478
479 // to find an initial point
480   Standard_Real Xi[4], Yi[4], UGuess;
481   findfourpoints(S, NInterval, InitSArray, InitUArray, myPrevS, myPrevU, Xi, Yi);
482   UGuess = cubic(S , Xi, Yi);
483
484   U = GCPnts_AbscissaPoint(C, deltaS, base, UGuess, myTolLen).Parameter();
485
486   // TODO - fields should be mutable
487   const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
488   const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;
489
490 #ifdef __OCC_DEBUG_CHRONO
491   ResultChron(chr_uparam, t_uparam);
492   uparam_count++;
493 #endif
494
495   return U;
496 }
497
498 Standard_Real Approx_CurvlinFunc::GetSParameter(Adaptor3d_Curve& C, const Standard_Real U, const Standard_Real Len) const
499 {
500   Standard_Real S, Origin;
501
502   Origin = C.FirstParameter();
503   S = myFirstS + Length(C, Origin, U)/Len;    
504   return S;
505 }
506
507 Standard_Boolean Approx_CurvlinFunc::EvalCase1(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
508 {
509   if(myCase != 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase1");
510
511   gp_Pnt C;
512   gp_Vec dC_dU, dC_dS, d2C_dU2, d2C_dS2;
513   Standard_Real U, Mag, dU_dS, d2U_dS2;
514   
515   U = GetUParameter(myC3D->GetCurve(), S, 1);
516
517   switch(Order) {
518
519   case 0: 
520     myC3D->D0(U, C);
521
522     Result(0) = C.X();
523     Result(1) = C.Y();
524     Result(2) = C.Z();
525     break;
526     
527   case 1:
528     myC3D->D1(U, C, dC_dU);
529     Mag = dC_dU.Magnitude();
530     dU_dS = myLength/Mag;
531     dC_dS = dC_dU*dU_dS;
532
533     Result(0) = dC_dS.X();
534     Result(1) = dC_dS.Y();
535     Result(2) = dC_dS.Z();
536     break;
537
538   case 2:
539     myC3D->D2(U, C, dC_dU, d2C_dU2);
540     Mag = dC_dU.Magnitude();
541     dU_dS = myLength/Mag;
542     d2U_dS2 = -myLength*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
543     d2C_dS2 = d2C_dU2*dU_dS*dU_dS + dC_dU*d2U_dS2;
544
545     Result(0) = d2C_dS2.X();
546     Result(1) = d2C_dS2.Y();
547     Result(2) = d2C_dS2.Z();
548     break;
549
550   default: Result(0) = Result(1) = Result(2) = 0;
551     return Standard_False;
552   }
553   return Standard_True;
554 }
555
556 Standard_Boolean Approx_CurvlinFunc::EvalCase2(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
557 {
558   if(myCase != 2) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase2");
559
560   Standard_Boolean Done;
561
562   Done = EvalCurOnSur(S, Order, Result, 1);
563
564   return Done;
565 }
566
567 Standard_Boolean Approx_CurvlinFunc::EvalCase3(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result)
568 {
569   if(myCase != 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase3");
570   
571   TColStd_Array1OfReal tmpRes1(0, 4), tmpRes2(0, 4);
572   Standard_Boolean Done;
573
574   Done = EvalCurOnSur(S, Order, tmpRes1, 1);
575
576   Done = EvalCurOnSur(S, Order, tmpRes2, 2) && Done;
577
578   Result(0) = tmpRes1(0);
579   Result(1) = tmpRes1(1);
580   Result(2) = tmpRes2(0);
581   Result(3) = tmpRes2(1);
582   Result(4) = 0.5*(tmpRes1(2) + tmpRes2(2));
583   Result(5) = 0.5*(tmpRes1(3) + tmpRes2(3));
584   Result(6) = 0.5*(tmpRes1(4) + tmpRes2(4));
585
586   return Done;
587 }
588
589 Standard_Boolean Approx_CurvlinFunc::EvalCurOnSur(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result, const  Standard_Integer NumberOfCurve) const
590 {
591   Handle(Adaptor2d_HCurve2d) Cur2D;
592   Handle(Adaptor3d_HSurface) Surf;
593   Standard_Real U=0, Length=0;
594
595   if (NumberOfCurve == 1) {
596     Cur2D = myC2D1;
597     Surf = mySurf1;    
598     Adaptor3d_CurveOnSurface CurOnSur(myC2D1, mySurf1);
599     U = GetUParameter(CurOnSur, S, 1);
600     if(myCase == 3) Length = myLength1;
601     else Length = myLength;
602   }
603   else if (NumberOfCurve == 2) {
604     Cur2D = myC2D2;
605     Surf = mySurf2;    
606     Adaptor3d_CurveOnSurface CurOnSur(myC2D2, mySurf2);
607     U = GetUParameter(CurOnSur, S, 2);
608     Length = myLength2;
609   }
610   else 
611     Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCurOnSur");
612
613   Standard_Real Mag, dU_dS, d2U_dS2, dV_dU, dW_dU, dV_dS, dW_dS, d2V_dS2, d2W_dS2, d2V_dU2, d2W_dU2;
614   gp_Pnt2d C2D;
615   gp_Pnt C;
616   gp_Vec2d dC2D_dU, d2C2D_dU2;
617   gp_Vec dC_dU, d2C_dU2, dC_dS, d2C_dS2, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW;
618
619   switch(Order) {
620   case 0:
621     Cur2D->D0(U, C2D);
622     Surf->D0(C2D.X(), C2D.Y(), C);
623     
624     Result(0) = C2D.X();
625     Result(1) = C2D.Y();
626     Result(2) = C.X();
627     Result(3) = C.Y();
628     Result(4) = C.Z();
629     break;
630
631   case 1:
632     Cur2D->D1(U, C2D, dC2D_dU);
633     dV_dU = dC2D_dU.X();
634     dW_dU = dC2D_dU.Y();
635     Surf->D1(C2D.X(), C2D.Y(), C, dS_dV, dS_dW);
636     dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
637     Mag = dC_dU.Magnitude();
638     dU_dS = Length/Mag;
639
640     dV_dS = dV_dU*dU_dS;
641     dW_dS = dW_dU*dU_dS;
642     dC_dS = dC_dU*dU_dS;
643
644     Result(0) = dV_dS;
645     Result(1) = dW_dS;
646     Result(2) = dC_dS.X();
647     Result(3) = dC_dS.Y();
648     Result(4) = dC_dS.Z();    
649     break;
650
651   case 2:
652     Cur2D->D2(U, C2D, dC2D_dU, d2C2D_dU2);
653     dV_dU = dC2D_dU.X();
654     dW_dU = dC2D_dU.Y();
655     d2V_dU2 = d2C2D_dU2.X();
656     d2W_dU2 = d2C2D_dU2.Y();
657     Surf->D2(C2D.X(), C2D.Y(), C, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW);
658     dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
659     d2C_dU2 = (d2S_dV2*dV_dU + d2S_dVdW*dW_dU)*dV_dU + dS_dV*d2V_dU2 + 
660               (d2S_dVdW*dV_dU + d2S_dW2*dW_dU)*dW_dU + dS_dW*d2W_dU2;
661     Mag = dC_dU.Magnitude();
662     dU_dS = Length/Mag;
663     d2U_dS2 = -Length*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
664     
665     dV_dS = dV_dU * dU_dS;
666     dW_dS = dW_dU * dU_dS;
667     d2V_dS2 = d2V_dU2*dU_dS*dU_dS + dV_dU*d2U_dS2;
668     d2W_dS2 = d2W_dU2*dU_dS*dU_dS + dW_dU*d2U_dS2;
669     
670     d2U_dS2 = -dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag);
671     d2C_dS2 = (d2S_dV2 * dV_dS + d2S_dVdW * dW_dS) * dV_dS + dS_dV * d2V_dS2 +
672               (d2S_dW2 * dW_dS + d2S_dVdW * dV_dS) * dW_dS + dS_dW * d2W_dS2;
673     
674     Result(0) = d2V_dS2;
675     Result(1) = d2W_dS2;
676     Result(2) = d2C_dS2.X();
677     Result(3) = d2C_dS2.Y();
678     Result(4) = d2C_dS2.Z();      
679     break;
680
681   default: Result(0) = Result(1) = Result(2) = Result(3) = Result(4) = 0;
682     return Standard_False;
683   }
684   return Standard_True;
685 }