242 {
243 HepVector cen1(lp1.
center());
244 HepVector cen2(lp2.
center());
245 double dx = cen1(1)-cen2(1);
246 double dy = cen1(2)-cen2(2);
247 double dc = sqrt(dx*dx+dy*dy);
249 double a1 = std::sqrt(lp1.alpha()) + std::sqrt(lp1.beta());
250 double a2 = std::sqrt(lp2.alpha()) + std::sqrt(lp2.beta());
251 double a3 = lp1.alpha()*lp2.alpha() + lp1.beta()*lp2.beta();
252 double det = lp1.alpha()*lp2.beta() - lp1.beta()*lp2.alpha();
253 if(fabs(det)>1e-12) {
254 double c1 = a2 * std::sqrt(lp1.
kappa()) + a1 * std::sqrt(lp2.
kappa()) -
257 double cinv = 1.0 /
c1;
258 double c2 = std::sqrt(a3) - 0.5 * (a1 + a2) - 2.0 * a3 *
259 (lp1.gamma() * lp2.
kappa() + lp2.gamma() * lp1.
kappa());
260 double c3 = a2 * std::sqrt(lp1.gamma()) + a1 * std::sqrt(lp2.gamma()) -
261 2.0 * a3 * lp1.gamma() * lp2.gamma();
262 double root = std::sqrt(c2) - 4.0 *
c1 * c3;
265 double rad2[2];
266 rad2[0] = 0.5 * cinv * (-c2 -
root);
267 rad2[1] = 0.5 * cinv * (-c2 +
root);
268 double ab1 = -(lp2.beta() * lp1.gamma() - lp1.beta() * lp2.gamma());
269 double ab2 = (lp2.alpha() * lp1.gamma() - lp1.alpha() * lp2.gamma());
270 double ac1 = -(lp2.beta() * lp1.
kappa() - lp1.beta() * lp2.
kappa());
271 double ac2 = (lp2.alpha() * lp1.
kappa() - lp1.alpha() * lp2.
kappa());
272 double dinv = 1.0 / det;
273 v1(1) = dinv * (ab1 + ac1 * rad2[0]);
274 v1(2) = dinv * (ab2 + ac2 * rad2[0]);
275 v1(3) = 0;
276 v2(1) = dinv * (ab1 + ac1 * rad2[1]);
277 v2(2) = dinv * (ab2 + ac2 * rad2[1]);
278 v2(3) = 0;
279 double d1 = lp1.
d(v1(1),v1(2));
280 double d2 = lp2.
d(v1(1),v1(2));
281 double d3 = lp1.
d(v2(1),v2(2));
282 double d4 = lp2.
d(v2(1),v2(2));
283 double r = sqrt(rad2[0]);
284 Lpar::Cpar cp1(lp1);
285 Lpar::Cpar cp2(lp2);
286 for(int j=0;j<2;j++) {
287 double s1,s2;
288 if(j==0) {
289 s1 = lp1.
s(v1(1),v1(2));
290 s2 = lp2.
s(v1(1),v1(2));
291 } else {
292 s1 = lp1.
s(v2(1),v2(2));
293 s2 = lp2.
s(v2(1),v2(2));
294 }
295 double phi1 = cp1.fi() + 2 * cp1.cu() * s1;
296 double phi2 = cp2.fi() + 2 * cp2.cu() * s2;
297 double f = (1 + 2 * cp1.cu() * cp1.da()) *
298 (1 + 2 * cp2.cu() * cp2.da()) *
cos(cp1.fi()-cp2.fi());
299 f -= 2 * (lp1.gamma() * lp2.
kappa() + lp2.gamma() * lp1.
kappa());
300 double cosphi12 = f;
301 }
302 return 2;
303 }
304 }
305 }
306 }
307 return 0;
308}
double cos(const BesAngle a)
double d(double x, double y) const
double s(double x, double y) const