1 package org.lcsim.recon.tracking.trfcylplane;
2
3 import org.lcsim.recon.tracking.trfbase.Propagator;
4 import org.lcsim.recon.tracking.trfbase.PropDirected;
5 import org.lcsim.recon.tracking.trfbase.PropDir;
6 import org.lcsim.recon.tracking.trfbase.TrackVector;
7 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
8 import org.lcsim.recon.tracking.trfbase.PropStat;
9 import org.lcsim.recon.tracking.trfutil.TRFMath;
10 import org.lcsim.recon.tracking.trfutil.Assert;
11 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
12 import org.lcsim.recon.tracking.trfzp.SurfZPlane;
13 import org.lcsim.recon.tracking.trfbase.Surface;
14 import org.lcsim.recon.tracking.trfbase.VTrack;
15 import org.lcsim.recon.tracking.trfbase.ETrack;
16
17 /**
18 * Propagates tracks from an ZPlane to a Cylinder in a constant field.
19 *<p>
20 * Propagation will fail if either the origin is not an ZPlane
21 * or destination is not a Cylinder.
22 * Propagator works incorrectly for tracks with very small curvatures.
23 *<p>
24 *@author Norman A. Graf
25 *@version 1.0
26 *
27 */
28
29 public class PropZCyl extends PropDirected
30 {
31 // Assign track parameter indices.
32 private static final int IX = SurfZPlane.IX;
33 private static final int IY = SurfZPlane.IY;
34 private static final int IDXDZ = SurfZPlane.IDXDZ;
35 private static final int IDYDZ = SurfZPlane.IDYDZ;
36 private static final int IQP_Z = SurfZPlane.IQP;
37
38 private static final int IPHI = SurfCylinder.IPHI;
39 private static final int IZ = SurfCylinder.IZ;
40 private static final int IALF = SurfCylinder.IALF;
41 private static final int ITLM = SurfCylinder.ITLM;
42 private static final int IQPT = SurfCylinder.IQPT;
43
44 // attributes
45
46 private double _bfac;
47
48 // static methods
49
50 // Return the type name.
51 /**
52 *Return a String representation of the class' type name.
53 *Included for completeness with the C++ version.
54 *
55 * @return A String representation of the class' type name.
56 */
57 public static String typeName()
58 { return "PropZCyl";
59 }
60
61 // Return the type.
62 /**
63 *Return a String representation of the class' type name.
64 *Included for completeness with the C++ version.
65 *
66 * @return A String representation of the class' type name.
67 */
68 public static String staticType()
69 { return typeName();
70 }
71
72
73
74 /**
75 *Construct an instance from a constant solenoidal magnetic field in Tesla.
76 *
77 * @param bfield The magnetic field strength in Tesla.
78 */
79 public PropZCyl(double bfield)
80 {
81 _bfac = TRFMath.BFAC*bfield;
82 }
83
84 /**
85 *Clone an instance.
86 *
87 * @return A Clone of this instance.
88 */
89 public Propagator newPropagator( )
90 {
91 return new PropZCyl( bField() );
92 }
93
94 /**
95 *Propagate a track without error in the specified direction.
96 *
97 * The track parameters for a cylinder are:
98 * phi z alpha tan(lambda) curvature
99 *
100 * @param trv The VTrack to propagate.
101 * @param srf The Surface to which to propagate.
102 * @param dir The direction in which to propagate.
103 * @return The propagation status.
104 */
105 public PropStat vecDirProp( VTrack trv, Surface srf,
106 PropDir dir)
107 {
108 TrackDerivative deriv = null;
109 return vecDirProp(trv, srf, dir, deriv);
110
111 }
112
113 /**
114 *Propagate a track without error in the specified direction
115 *and return the derivative matrix in deriv.
116 *
117 * The track parameters for a cylinder are:
118 * phi z alpha tan(lambda) curvature
119 *
120 * @param trv The VTrack to propagate.
121 * @param srf The Surface to which to propagate.
122 * @param dir The direction in which to propagate.
123 * @param deriv The track derivatives to update at the surface srf.
124 * @return The propagation status.
125 */
126 public PropStat vecDirProp( VTrack trv, Surface srf,
127 PropDir dir,TrackDerivative deriv)
128 {
129 return vecPropagateZCyl( _bfac, trv, srf, dir, deriv );
130 }
131
132 /**
133 *Return the strength of the magnetic field in Tesla.
134 *
135 * @return The strength of the magnetic field in Tesla.
136 */
137 public double bField()
138 {
139 return _bfac/TRFMath.BFAC;
140 }
141
142 /**
143 *Return a String representation of the class' type name.
144 *Included for completeness with the C++ version.
145 *
146 * @return A String representation of the class' type name.
147 */
148 public String type()
149 { return staticType();
150 }
151
152 /**
153 *output stream
154 * @return A String representation of this instance.
155 */
156 public String toString()
157 {
158 return "ZPlane-Cylinder propagation with constant "
159 + bField() + " Tesla field";
160
161 }
162
163
164 //**********************************************************************
165
166 // Private function to propagate a track without error
167 //
168 // The track parameters for a cylinder are:
169 // phi z alpha tan(lambda) curvature
170 // (NOT [. . . lambda .] as in TRF and earlier version of TRF++.)
171 //
172 // If pderiv is nonzero, return the derivative matrix there.
173 // On Cylinder:
174 // r (cm) is fixed
175 // 0 - phi
176 // 1 - z (cm)
177 // 2 - alp
178 // 3 - tlam
179 // 4 - q/pt pt is transverse momentum of a track, q is its charge
180 // On ZPlane:
181 // 0 - x (cm)
182 // 1 - y (cm)
183 // 2 - dx/dz
184 // 3 - dy/dz
185 // 4 - q/p p is momentum of a track, q is its charge
186 // If pderiv is nonzero, return the derivative matrix there.
187
188 private PropStat
189 vecPropagateZCyl( double bfac, VTrack trv, Surface srf,
190 PropDir dir,
191 TrackDerivative deriv )
192 {
193
194 // construct return status
195 PropStat pstat = new PropStat();
196
197 // fetch the originating surface and vector
198 Surface srf1 = trv.surface();
199 TrackVector vec1 = trv.vector();
200
201 // Check origin is a Zplane.
202 Assert.assertTrue( srf1.pureType().equals(SurfZPlane.staticType()) );
203 if ( !srf1.pureType( ).equals(SurfZPlane.staticType()) )
204 return pstat;
205 SurfZPlane szp1 = ( SurfZPlane ) srf1;
206
207 // Check destination is a cylinder.
208 Assert.assertTrue( srf.pureType().equals(SurfCylinder.staticType()) );
209 if ( !srf.pureType( ).equals(SurfCylinder.staticType()) )
210 return pstat;
211 SurfCylinder scy2 = ( SurfCylinder ) srf;
212
213
214 // Fetch the z of the plane and the starting track vector.
215 int iz = SurfZPlane.ZPOS;
216 double z = szp1.parameter(iz);
217
218 TrackVector vec = trv.vector();
219 double x = vec.get(IX); // v
220 double y = vec.get(IY); // z
221 double b = vec.get(IDXDZ); // dv/du
222 double a = vec.get(IDYDZ); // dz/du
223 double e = vec.get(IQP_Z); // q/p
224
225 // Fetch the radii and the starting track vector.
226 int irad = SurfCylinder.RADIUS;
227 double r2 = scy2.parameter(irad);
228
229 int sign_dz = 0;
230 if(trv.isForward()) sign_dz = 1;
231 if(trv.isBackward()) sign_dz = -1;
232 if(sign_dz == 0)
233 {
234 System.out.println("PropZCyl._vec_propagate: Unknown direction of a track ");
235 System.exit(1);
236 }
237
238 // Calculate cylindrical coordinates
239
240 double cnst1=x>0?0.:Math.PI;
241 double cnst2=Math.PI;
242 if(a>0.&&b>0.&&sign_dz>0.) cnst2=0.;
243 if(a<0.&&b>0.&&sign_dz>0.) cnst2=0.;
244 if(a>0.&&b<0.&&sign_dz<0.) cnst2=0.;
245 if(a<0.&&b<0.&&sign_dz<0.) cnst2=0.;
246
247 double sign_y= y>0 ? 1.:-1;
248 double sign_a= a>0 ? 1.:-1;
249 if(y==0) sign_y=0.;
250 if(a==0) sign_a=0.;
251 double atnxy=(x!=0.?Math.atan(y/x):sign_y*Math.PI/2.);
252 double atnab=(b!=0.?Math.atan(a/b):sign_a*Math.PI/2.);
253
254 double phi1= TRFMath.fmod2(atnxy+cnst1,TRFMath.TWOPI);
255 double r1=Math.sqrt(x*x+y*y);
256 // 28jan01 dla
257 // Handle r1 = 0 to avaoid crash later on.
258 if ( r1 < 0.001 )
259 {
260 return pstat;
261 }
262 double z1=z;
263 double alp1= TRFMath.fmod2(atnab-phi1+cnst2,TRFMath.TWOPI);
264 double tlm1= sign_dz/Math.sqrt(a*a+b*b);
265 double qpt1=e*Math.sqrt((1+a*a+b*b)/(a*a+b*b));
266
267
268 // Check alpha range.
269 alp1 = TRFMath.fmod2( alp1, TRFMath.TWOPI );
270 Assert.assertTrue( Math.abs(alp1) <= Math.PI );
271
272 //if ( trv.is_forward() ) Assert.assertTrue( fabs(alp1) <= PI2 );
273 //else Assert.assertTrue( fabs(alp1) > PI2 );
274
275 // Calculate the cosine of lambda.
276 //double clam1 = 1.0/sqrt(1+tlm1*tlm1);
277
278 // Calculate curvature: C = _bfac*(q/p)/cos(lambda)
279 // and its derivatives
280 // Assert.assertTrue( clam1 != 0.0 );
281 // double dcrv1_dqp1 = bfac/clam1;
282 // double crv1 = dcrv1_dqp1*qp1;
283 // double dcrv1_dtlm1 = crv1*clam1*clam1*tlm1;
284
285 // Calculate the curvature = _bfac*(q/pt)
286 double dcrv1_dqpt1 = bfac;
287 double crv1 = dcrv1_dqpt1*qpt1;
288 //double dcrv1_dtlm1 = 0.0;
289
290 // Evaluate the new track vector.
291 // See dla log I-044
292
293 // lambda and curvature do not change
294 double tlm2 = tlm1;
295 double crv2 = crv1;
296 double qpt2 = qpt1;
297
298 // We can evaluate sin(alp2), leaving two possibilities for alp2
299 // 1st solution: alp21, phi21, phid21, tht21
300 // 2nd solution: alp22, phi22, phid22, tht22
301 // evaluate phi2 to choose
302 double salp1 = Math.sin( alp1 );
303 double calp1 = Math.cos( alp1 );
304 double salp2 = r1/r2*salp1 + 0.5*crv1/r2*(r2*r2-r1*r1);
305 // if salp2 > 1, track does not cross cylinder
306 if ( Math.abs(salp2) > 1.0 ) return pstat;
307 double alp21 = Math.asin( salp2 );
308 double alp22 = alp21>0 ? Math.PI-alp21 : -Math.PI-alp21;
309 double calp21 = Math.cos( alp21 );
310 double calp22 = Math.cos( alp22 );
311 double phi20 = phi1 + Math.atan2( salp1-r1*crv1, calp1 );
312 double phi21 = phi20 - Math.atan2( salp2-r2*crv2, calp21 ); // phi position
313 double phi22 = phi20 - Math.atan2( salp2-r2*crv2, calp22 );
314
315 // Construct an sT object for each solution.
316 STCalcZ sto1 = new STCalcZ(r1,phi1,alp1,crv1,r2,phi21,alp21);
317 STCalcZ sto2 = new STCalcZ(r1,phi1,alp1,crv1,r2,phi22,alp22);
318 // Check the two solutions are nonzero and have opposite sign
319 // or at least one is nonzero.
320
321 // Choose the correct solution
322 boolean use_first_solution = false;
323 if (dir.equals(PropDir.NEAREST))
324 {
325 use_first_solution = Math.abs(sto2.st()) > Math.abs(sto1.st());
326 }
327 else if (dir.equals(PropDir.FORWARD))
328 {
329 use_first_solution = sto1.st() > 0.0;
330 }
331 else if (dir.equals(PropDir.BACKWARD))
332 {
333 use_first_solution = sto1.st() < 0.0;
334 }
335 else
336 {
337 System.out.println( "PropCyl._vec_propagate: Unknown direction.");
338 System.exit(1);
339 }
340
341 // Assign phi2, alp2 and sto2 for the chosen solution.
342 double phi2, alp2;
343 STCalcZ sto;
344 double calp2;
345 if ( use_first_solution )
346 {
347 sto = sto1;
348 phi2 = phi21;
349 alp2 = alp21;
350 calp2 = calp21;
351 }
352 else
353 {
354 sto = sto2;
355 phi2 = phi22;
356 alp2 = alp22;
357 calp2 = calp22;
358 }
359
360 // fetch sT.
361 double st = sto.st();
362
363 // use sT to evaluate z2
364 double z2 = z1 + tlm1*st;
365
366 // Check alpha range.
367 Assert.assertTrue( Math.abs(alp2) <= Math.PI );
368
369 // put new values in vec
370 vec.set(IPHI , phi2);
371 vec.set(IZ , z2);
372 vec.set(IALF , alp2);
373 vec.set(ITLM , tlm2);
374 vec.set(IQPT , qpt2);
375
376 // Update trv
377 trv.setSurface(srf.newPureSurface());
378 trv.setVector(vec);
379 if ( Math.abs(alp2) <= TRFMath.PI2 ) trv.setForward();
380 else trv.setBackward();
381
382 // Evaluate s.
383 double s = st*Math.sqrt(1.0+tlm2*tlm2);
384
385 // Set the return status.
386 pstat.setPathDistance(s);
387 //st > 0 ? pstat.set_forward() : pstat.set_backward();
388
389 // exit now if user did not ask for error matrix.
390 if ( deriv == null ) return pstat;
391
392 // Calculate derivatives.
393 // dphi1
394
395 double dphi1_dx= -y/(x*x+y*y);
396 double dphi1_dy= x/(x*x+y*y);
397
398 // dz1
399
400 //double dz1_dz=0.;
401
402 // dalf1
403
404 double dalp1_da= b/(a*a+b*b);
405 double dalp1_db= -a/(a*a+b*b);
406 double dalp1_dy= -dphi1_dy;
407 double dalp1_dx= -dphi1_dx;
408
409 // dr1
410
411 double dr1_dx= x/Math.sqrt(x*x+y*y);
412 double dr1_dy= y/Math.sqrt(x*x+y*y);
413
414 // dtlm1
415
416 double dtlm1_da= -sign_dz*a/Math.sqrt(a*a+b*b)/(a*a+b*b);
417 double dtlm1_db= -sign_dz*b/(a*a+b*b)/Math.sqrt(a*a+b*b);
418
419 // dcrv1
420
421 double dqpt1_de= Math.sqrt((1+a*a+b*b)/(a*a+b*b));
422 double dqpt1_da= -a*e/Math.sqrt((a*a+b*b)*(1+a*a+b*b))/(a*a+b*b);
423 double dqpt1_db= -e*b/Math.sqrt(1+a*a+b*b)/Math.sqrt(a*a+b*b)/(a*a+b*b);
424
425 double dcrv1_de= dqpt1_de*bfac;
426 double dcrv1_da= dqpt1_da*bfac;
427 double dcrv1_db= dqpt1_db*bfac;
428
429 // alpha_2
430 double da2da1 = r1*calp1/r2/calp2;
431 double da2dc1 = (r2*r2-r1*r1)*0.5/r2/calp2;
432 double da2dr1 = (salp1-crv2*r1)/r2/calp2;
433
434 // phi2
435 double rcsal1 = r1*crv1*salp1;
436 double rcsal2 = r2*crv2*salp2;
437 double den1 = 1.0 + r1*r1*crv1*crv1 - 2.0*rcsal1;
438 double den2 = 1.0 + r2*r2*crv2*crv2 - 2.0*rcsal2;
439 double dp2dp1 = 1.0;
440 double dp2da1 = (1.0-rcsal1)/den1 - (1.0-rcsal2)/den2*da2da1;
441 double dp2dc1 = -r1*calp1/den1 + r2*calp2/den2
442 - (1.0-rcsal2)/den2*da2dc1;
443 double dp2dr1= -crv1*calp1/den1-(1.0-rcsal2)*da2dr1/den2;
444
445 // z2
446 //double dz2dz1 = 1.0;
447 double dz2dl1 = st;
448 double dz2da1 = tlm1*sto.d_st_dalp1(dp2da1, da2da1);
449 double dz2dc1 = tlm1*sto.d_st_dcrv1(dp2dc1, da2dc1);
450 double dz2dr1 = tlm1*sto.d_st_dr1( dp2dr1, da2dr1);
451
452
453 // final derivatives
454
455 // phi2
456 double dphi2_dx=dp2dp1*dphi1_dx+dp2da1*dalp1_dx+dp2dr1*dr1_dx;
457 double dphi2_dy=dp2dp1*dphi1_dy+dp2da1*dalp1_dy+dp2dr1*dr1_dy;
458 double dphi2_db=dp2da1*dalp1_db+dp2dc1*dcrv1_db;
459 double dphi2_da=dp2da1*dalp1_da+dp2dc1*dcrv1_da;
460 double dphi2_de=dp2dc1*dcrv1_de;
461
462 // alp2
463 double dalp2_dx= da2da1*dalp1_dx+da2dr1*dr1_dx;
464 double dalp2_dy= da2da1*dalp1_dy+da2dr1*dr1_dy;
465 double dalp2_db= da2da1*dalp1_db+da2dc1*dcrv1_db;
466 double dalp2_da= da2da1*dalp1_da+da2dc1*dcrv1_da;
467 double dalp2_de= da2dc1*dcrv1_de;
468
469 // crv2
470 double dqpt2_da=dqpt1_da;
471 double dqpt2_db=dqpt1_db;
472 double dqpt2_de=dqpt1_de;
473
474 // tlm2
475 double dtlm2_da= dtlm1_da;
476 double dtlm2_db= dtlm1_db;
477
478 // z2
479 double dz2_dx= dz2dr1*dr1_dx+dz2da1*dalp1_dx;
480 double dz2_dy= dz2dr1*dr1_dy+dz2da1*dalp1_dy;
481 double dz2_db= dz2da1*dalp1_db+dz2dl1*dtlm1_db+dz2dc1*dcrv1_db;
482 double dz2_da= dz2da1*dalp1_da+dz2dl1*dtlm1_da+dz2dc1*dcrv1_da;
483 double dz2_de= dz2dc1*dcrv1_de;
484
485
486 // Build derivative matrix.
487
488 deriv.set(IPHI,IX , dphi2_dx);
489 deriv.set(IPHI,IY , dphi2_dy);
490 deriv.set(IPHI,IDXDZ, dphi2_db);
491 deriv.set(IPHI,IDYDZ, dphi2_da);
492 deriv.set(IPHI,IQP_Z, dphi2_de);
493 deriv.set(IZ,IX , dz2_dx);
494 deriv.set(IZ,IY , dz2_dy);
495 deriv.set(IZ,IDXDZ , dz2_db);
496 deriv.set(IZ,IDYDZ , dz2_da);
497 deriv.set(IZ,IQP_Z, dz2_de);
498 deriv.set(IALF,IX , dalp2_dx);
499 deriv.set(IALF,IY , dalp2_dy);
500 deriv.set(IALF,IDXDZ , dalp2_db);
501 deriv.set(IALF,IDYDZ , dalp2_da);
502 deriv.set(IALF,IQP_Z , dalp2_de);
503 deriv.set(ITLM,IDXDZ , dtlm2_db);
504 deriv.set(ITLM,IDYDZ , dtlm2_da);
505 deriv.set(IQPT,IDXDZ , dqpt2_db);
506 deriv.set(IQPT,IDYDZ , dqpt2_da);
507 deriv.set(IQPT,IQP_Z , dqpt2_de);
508
509
510 return pstat;
511
512 }
513
514
515 //**********************************************************************
516 // helpers
517 //**********************************************************************
518
519 // Private class STCalcZ.
520 //
521 // An STCalcZ_ object calculates sT (the signed transverse path length)
522 // and its derivatives w.r.t. alf1 and crv1. It is constructed from
523 // the starting (r1, phi1, alf1, crv1) and final track parameters
524 // (r2, phi2, alf2) assuming these are consistent. Methods are
525 // provided to retrieve sT and the two derivatives.
526
527 class STCalcZ
528 {
529
530 private boolean _big_crv;
531 private double _st;
532 private double _dst_dphi21;
533 private double _dst_dcrv1;
534 private double _dst_dr1;
535 private double _cnst1,_cnst2;
536 public double _crv1;
537 // constructor
538 public STCalcZ()
539 {
540 }
541 public STCalcZ(double r1, double phi1, double alf1, double crv1,
542 double r2, double phi2, double alf2)
543 {
544
545 _crv1 = crv1;
546 Assert.assertTrue( r1 > 0.0 );
547 Assert.assertTrue( r2 > 0.0 );
548 double rmax = r1+r2;
549
550 // Calculate the change in xy direction
551 double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
552 Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
553
554 // Evaluate whether |C| is" big"
555 _big_crv = rmax*Math.abs(crv1) > 0.001 || Math.abs(r1-r2)<1.e-9;
556 if( Math.abs(crv1) < 1.e-10 ) _big_crv=false;
557
558 // If the curvature is big we can use
559 // sT = (phi_dir2 - phi_dir1)/crv1
560 if ( _big_crv )
561 {
562 Assert.assertTrue( crv1 != 0.0 );
563 _st = phi_dir_diff/crv1;
564 }
565
566 // Otherwise, we calculate the straight-line distance
567 // between the points and use an approximate correction
568 // for the (small) curvature.
569 else
570 {
571
572 // evaluate the distance
573 double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
574 double arg = 0.5*d*crv1;
575 double arg2 = arg*arg;
576 double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
577 _st = d + st_minus_d;
578
579 // evaluate the sign
580 // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
581 // Because sT*C = dphid and d = abs(sT):
582 // xsign = 0 for sT > 0
583 // xsign = 2 for sT < 0
584 // Numerical roundoff will smear these predictions.
585 double xsign = (crv1 == 0. ? 0.: Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1)) );
586 double sign = 0.0;
587 if ( crv1 != 0 )
588 {
589 if ( xsign < 0.5 ) sign = 1.0;
590 if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
591 }
592 // If the above is indeterminate, assume zero curvature.
593 // In this case abs(alpha) decreases monotonically
594 // with sT. Track passing through origin has alpha = 0 on one
595 // side and alpha = +/-pi on the other. If both points are on
596 // the same side, we use dr/ds > 0 for |alpha|<pi/2.
597 if ( sign == 0 )
598 {
599 sign = 1.0;
600 if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
601 if ( Math.abs(alf2) == Math.abs(alf1) )
602 {
603 if ( Math.abs(alf2) < TRFMath.PI2 )
604 {
605 if ( r2 < r1 ) sign = -1.0;
606 }
607 else
608 {
609 if ( r2 > r1 ) sign = -1.0;
610 }
611 }
612 }
613
614 if ( Math.abs(r2 - r1)<1e-12 && crv1 == 0. && Math.abs(alf1) < TRFMath.PI2 ) sign = -1.0;
615 if ( Math.abs(r2 - r1)<1e-12 && crv1 == 0. && Math.abs(alf1) > TRFMath.PI2 ) sign = 1.0;
616
617
618 // Correct _st using the above sign.
619 Assert.assertTrue( Math.abs(sign) == 1.0 );
620 _st = sign*_st;
621
622 // save derivatives
623 _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
624 double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
625 if( Math.abs(_st) > 1e-11 )
626 {
627 _dst_dphi21 = sign*(r1*r2*Math.sin(phi2-phi1))*root/d;
628 _dst_dr1= (1.+arg2/2.*(1+3./4.*arg2))/d*sign;
629 _cnst1=r1-r2*Math.cos(phi2-phi1);
630 _cnst2=r1*r2*Math.sin(phi2-phi1);
631 }
632 else
633 {
634 _dst_dphi21= r1;
635 _dst_dr1= sign;
636 _cnst1=0.;
637 _cnst2=0.;
638 }
639 }
640
641 }
642 public double st()
643 { return _st;
644 }
645 public double d_st_dalp1(double d_phi2_dalf1, double d_alf2_dalf1 )
646 {
647 if ( _big_crv ) return ( d_phi2_dalf1 + d_alf2_dalf1 - 1.0 ) / _crv1;
648 else return _dst_dphi21 * d_phi2_dalf1;
649 }
650 public double d_st_dcrv1(double d_phi2_dcrv1, double d_alf2_dcrv1 )
651 {
652 if ( _big_crv ) return ( d_phi2_dcrv1 + d_alf2_dcrv1 - _st ) / _crv1;
653 else return _dst_dcrv1 + _dst_dphi21*d_phi2_dcrv1;
654
655 }
656 public double d_st_dr1( double d_phi2_dr1, double d_alf2_dr1 )
657 {
658 if ( _big_crv ) return ( d_phi2_dr1 + d_alf2_dr1 ) / _crv1;
659 else
660 {
661 if( Math.abs(_cnst2)<1e-12 && Math.abs(_cnst1)<1e-12 )
662 {
663 return _dst_dr1*Math.sqrt(1. + _dst_dphi21*_dst_dphi21*d_phi2_dr1*d_phi2_dr1);
664 }
665 else
666 return _dst_dr1*(_cnst1+_cnst2*d_phi2_dr1);
667 }
668
669 }
670
671 }
672
673 }