1 package org.lcsim.recon.tracking.trfxyp; 2 3 import org.lcsim.recon.tracking.trfutil.TRFMath; 4 import org.lcsim.recon.tracking.trfutil.Assert; 5 6 import org.lcsim.recon.tracking.trfbase.Propagator; 7 import org.lcsim.recon.tracking.trfbase.PropDirected; 8 import org.lcsim.recon.tracking.trfbase.PropDir; 9 import org.lcsim.recon.tracking.trfbase.PropStat; 10 import org.lcsim.recon.tracking.trfbase.VTrack; 11 import org.lcsim.recon.tracking.trfbase.Surface; 12 import org.lcsim.recon.tracking.trfbase.TrackVector; 13 import org.lcsim.recon.tracking.trfbase.TrackDerivative; 14 15 /** 16 * Propagates tracks from one XYPlane to another in a constant field 17 * directed along the v axis (parallel to the surfaces of origin and 18 * destination and normal to the z axis). 19 *<p> 20 * Propagation will fail if either the origin or destination is 21 * not an XYPlane or if they are not parallel. 22 * 23 *@author Norman A. Graf 24 *@version 1.0 25 * 26 */ 27 public class PropXYXYBV extends PropDirected 28 { 29 30 // attributes 31 32 private double _bfac; 33 34 private static final int IV = SurfXYPlane.IV; 35 private static final int IZ = SurfXYPlane.IZ; 36 private static final int IDVDU = SurfXYPlane.IDVDU; 37 private static final int IDZDU = SurfXYPlane.IDZDU; 38 private static final int IQP = SurfXYPlane.IQP; 39 40 // static methods 41 42 /** 43 *Return a String representation of the class' type name. 44 *Included for completeness with the C++ version. 45 * 46 * @return A String representation of the class' type name. 47 */ 48 public static String typeName() 49 { return "PropXYXYBV"; } 50 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 staticType() 58 { return typeName(); } 59 60 /** 61 *Construct an instance from a constant solenoidal magnetic field in Tesla. 62 * 63 * @param bfield The magnetic field strength in Tesla. 64 */ 65 public PropXYXYBV(double bfield) 66 { 67 _bfac = TRFMath.BFAC*bfield; 68 } 69 70 /** 71 *Clone an instance. 72 * 73 * @return A Clone of this instance. 74 */ 75 public Propagator newPropagator( ) 76 { 77 return new PropXYXYBV( bField() ); 78 } 79 80 /** 81 *Propagate a track without error in the specified direction 82 *and return the derivative matrix in deriv. 83 * 84 * @param trv The VTrack to propagate. 85 * @param srf The Surface to which to propagate. 86 * @param dir The direction in which to propagate. 87 * @param deriv The track derivatives to update at the surface srf. 88 * @return The propagation status. 89 **/ 90 public PropStat vecDirProp( VTrack trv, Surface srf, 91 PropDir dir,TrackDerivative deriv ) 92 { 93 PropStat pstat = vec_propagatexyxy_bv_( _bfac, trv, srf, dir, deriv ); 94 return pstat; 95 } 96 97 /** 98 *Propagate a track without error in the specified direction. 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 *Return the strength of the magnetic field in Tesla. 115 * 116 * @return The strength of the magnetic field in Tesla. 117 */ 118 public double bField() 119 { 120 return _bfac/TRFMath.BFAC; 121 } 122 123 /** 124 *Return a String representation of the class' type name. 125 *Included for completeness with the C++ version. 126 * 127 * @return A String representation of the class' type name. 128 */ 129 public String type() 130 { return staticType(); } 131 132 133 /** 134 *output stream 135 * 136 * @return A String representation of this instance. 137 */ 138 public String toString() 139 { 140 return "XYPlane-XYPlane propagation with constant " + bField() 141 + " Tesla field normal to z and parallel to the planes"; 142 } 143 144 145 146 // Private function to determine dphi of the propagation 147 double directionBV(int flag_forward, int sign_dphi, 148 double du, 149 double norm, double cat, 150 double sinphi, double cosphi) 151 { 152 153 int sign_cat = 0; 154 if(du*flag_forward>0.) sign_cat = 1; 155 if(du*flag_forward<0.) sign_cat = -1; 156 if(du == 0.) 157 { 158 if(sinphi >=0. ) sign_cat = 1; 159 if(sinphi < 0. ) sign_cat = -1; 160 } 161 162 double sin_dphi = norm*sinphi + cat*cosphi*sign_cat; 163 double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat; 164 165 int sign_sindphi = 0; 166 if(sin_dphi> 0.) sign_sindphi = 1; 167 if(sin_dphi< 0.) sign_sindphi = -1; 168 if(sin_dphi == 0.) sign_sindphi = sign_dphi; 169 170 double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi); 171 return dphi; 172 173 } 174 //********************************************************************** 175 176 // Private function to propagate a track without error 177 // The corresponding track parameters are: 178 // u (cm) is fixed 179 // 0 - v (cm) 180 // 1 - z (cm) 181 // 2 - dv/du 182 // 3 - dz/du 183 // 4 - q/p p is momentum of a track, q is its charge 184 // If pderiv is nonzero, return the derivative matrix there. 185 186 // We use the code from PropXYXY written for the case of B^ || z^. 187 // First, we rotate the track to the c.f. where B^ || z^, propagate the 188 // track and rotate the propagated track back to the original frame. 189 // We require the two planes to be parallel to avoid the corner regions 190 // where the approximation B^ || v^ is not valid. 191 192 PropStat 193 vec_propagatexyxy_bv_( double B, VTrack trv, Surface srf, 194 PropDir dir, 195 TrackDerivative deriv ) 196 { 197 198 // construct return status 199 PropStat pstat = new PropStat(); 200 201 // fetch the originating surface and vector 202 Surface srf1 = trv.surface(); 203 // TrackVector vec1 = trv.get_vector(); 204 205 // Check origin is a XYPlane. 206 Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) ); 207 if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) ) 208 return pstat; 209 SurfXYPlane sxyp1 = ( SurfXYPlane ) srf1; 210 211 // Check destination is a XYPlane. 212 Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) ); 213 if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) ) 214 return pstat; 215 SurfXYPlane sxyp2 = ( SurfXYPlane ) srf; 216 217 // If surfaces are the same, we can return now. 218 if ( srf.pureEqual(srf1) ) 219 { 220 if(deriv != null) 221 { 222 deriv.setIdentity(); 223 } 224 pstat.setSame(); 225 return pstat; 226 } 227 228 // Fetch the u's and phi's of the planes and the starting track vector. 229 int iphi = SurfXYPlane.NORMPHI; 230 int idist = SurfXYPlane.DISTNORM; 231 double phi_0 = sxyp1.parameter(iphi); 232 double u_0 = sxyp1.parameter(idist); 233 double phi_n = sxyp2.parameter(iphi); 234 double u_n = sxyp2.parameter(idist); 235 double phi_u = phi_0 - phi_n; 236 237 // check here that first planes are parallel and second that 238 //planes have phi==0||PI , otherwise this propagator doesn't work. 239 240 if (Math.abs(phi_0)!=0. && Math.abs(phi_0)!=Math.PI ) return pstat; 241 if ( Math.abs(phi_u)!=0. && Math.abs(phi_u)!=Math.PI && Math.abs(phi_u)!=TRFMath.TWOPI ) return pstat; 242 243 //cout<<" PropXYXYBV "<<endl; 244 //cout<<trv<<endl; 245 TrackVector vec = trv.vector(); 246 247 // b_0 is rotated by 90 deg around the u axis wrt to the input track vector 248 // (v',z') = (-z,v) 249 double b1_0 = -vec.get(IZ); 250 double b2_0 = vec.get(IV); 251 double b3_0 = -vec.get(IDZDU); 252 double b4_0 = vec.get(IDVDU); 253 double b5_0 = vec.get(IQP); 254 255 256 double cosphi_u = Math.cos(phi_u); 257 double sinphi_u = Math.sin(phi_u); 258 259 // check if du == 0 ( that is track moves parallel to the destination plane ) 260 double du_du_0 = cosphi_u - b3_0*sinphi_u; 261 if(du_du_0==0.) return pstat; 262 263 double a_hat_u = 1./du_du_0; 264 double a_hat_u2 = a_hat_u*a_hat_u; 265 266 double u = u_0*cosphi_u - b1_0*sinphi_u; 267 double b1 = b1_0*cosphi_u + u_0*sinphi_u; 268 double b2 = b2_0; 269 double b3 = (b3_0*cosphi_u + sinphi_u)*a_hat_u; 270 double b4 = b4_0*a_hat_u; 271 double b5 = b5_0; 272 273 int sign_du0 = 0; 274 if(trv.isForward()) sign_du0 = 1; 275 if(trv.isBackward()) sign_du0 = -1; 276 if(sign_du0 == 0) 277 { 278 System.out.println("PropXYXYBV._vec_prop: Unknown direction of a track "); 279 System.exit(1); 280 } 281 int sign_du = 0; 282 if(du_du_0*sign_du0 > 0) sign_du = 1; 283 if(du_du_0*sign_du0 < 0) sign_du = -1; 284 285 // check that q/p != 0 286 Assert.assertTrue(b5 !=0. ); 287 288 // 1/curv of the track is r 289 double r = 1/(b5*B)*Math.sqrt(1 + b3_0*b3_0)/Math.sqrt(1 + b3_0*b3_0 + b4_0*b4_0); 290 double b3_hat = Math.sqrt(1 + b3*b3); 291 double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4); 292 double b3_hat2 = b3_hat*b3_hat; 293 double b34_hat2 = b34_hat*b34_hat; 294 double cosphi = -b3*sign_du/b3_hat; 295 double sinphi = sign_du/b3_hat; 296 double rsinphi = 1./(b5*B)*sign_du/b34_hat; 297 double rcosphi = -b3/(b5*B)*sign_du/b34_hat; 298 299 double du = u_n - u; 300 double norm = du/r - cosphi; 301 302 // xy-xy propagation failed : noway to the new plane 303 if(Math.abs(norm)>1.) return pstat; 304 305 double cat = Math.sqrt(1.-norm*norm); 306 int flag_forward = 0; 307 int sign_dphi = 0; 308 309 if ( dir.equals(PropDir.NEAREST) ) 310 { 311 { 312 // try forward propagation 313 flag_forward = 1; 314 if(b5>0.) sign_dphi = 1; 315 if(b5<0.) sign_dphi = -1; 316 double dphi1= 317 directionBV(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi); 318 // try backward propagation 319 flag_forward = -1; 320 if(b5>0.) sign_dphi = -1; 321 if(b5<0.) sign_dphi = 1; 322 double dphi2= 323 directionBV(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi); 324 if( Math.abs(dphi2)>Math.abs(dphi1) ) 325 { 326 flag_forward = -flag_forward; 327 sign_dphi = - sign_dphi; 328 } 329 } 330 } 331 else if ( dir.equals(PropDir.FORWARD) ) 332 { 333 flag_forward = 1; 334 if(b5>0.) sign_dphi = 1; 335 if(b5<0.) sign_dphi = -1; 336 } 337 else if ( dir.equals(PropDir.BACKWARD) ) 338 { 339 flag_forward = -1; 340 if(b5>0.) sign_dphi = -1; 341 if(b5<0.) sign_dphi = 1; 342 } 343 else 344 { 345 System.out.println( "PropXYXY._vec_prop: Unknown direction."); 346 System.exit(1); 347 } 348 349 int sign_cat = 0; 350 if(du*sign_dphi*b5>0.) sign_cat = 1; 351 if(du*sign_dphi*b5<0.) sign_cat = -1; 352 if(du == 0.) 353 { 354 if(sinphi >=0. ) sign_cat = 1; 355 if(sinphi < 0. ) sign_cat = -1; 356 } 357 358 double sin_dphi = norm*sinphi + cat*cosphi*sign_cat; 359 double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat; 360 if(cos_dphi>1.) cos_dphi=1.; 361 if(cos_dphi<-1.) cos_dphi= -1.; 362 363 int sign_sindphi = 0; 364 if(sin_dphi> 0.) sign_sindphi = 1; 365 if(sin_dphi< 0.) sign_sindphi = -1; 366 if(sin_dphi == 0.) sign_sindphi = sign_dphi; 367 368 double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi); 369 370 // check that I didn't make any mistakes 371 372 Assert.assertTrue(Math.abs(Math.sin(dphi) -sin_dphi)<1.e-5); 373 374 // check if dun == 0 ( that is track moves parallel to the destination plane) 375 double du_n_du = cos_dphi - b3*sin_dphi; 376 if(du_n_du==0.) return pstat; 377 378 double a_hat_dphi = 1./du_n_du; 379 double a_hat_dphi2 = a_hat_dphi*a_hat_dphi; 380 double c_hat_dphi = sin_dphi + b3*cos_dphi ; 381 382 double b1_n = b1 + rsinphi*(1-cos_dphi) - rcosphi*sin_dphi; 383 double b2_n = b2 + b4/(b5*B)*sign_du/b34_hat*dphi; 384 double b3_n = c_hat_dphi*a_hat_dphi; 385 double b4_n = b4*a_hat_dphi; 386 double b5_n = b5; 387 388 // double u_n_0 = u_n*cosphi_u + b1_n*sinphi_u; 389 390 // check if track crossed original plane during the propagation 391 // switch (dir) { 392 // if( dir.equals(PropDir.FORWARD: 393 // if((u_n_0 - u_0)*sign_du0<0) return pstat; 394 // break; 395 // if( dir.equals(PropDir.BACKWARD: 396 // if((u_n_0 - u_0)*sign_du0>0) return pstat; 397 // break; 398 // } 399 400 int sign_dun = 0; 401 if(du_n_du*sign_du > 0) sign_dun = 1; 402 if(du_n_du*sign_du < 0) sign_dun = -1; 403 404 // rotate the propagated track back by -90 deg around the u axis: 405 // (v,z) = (z',-v'). 406 vec.set(IV , b2_n); 407 vec.set(IZ , -b1_n); 408 vec.set(IDVDU , b4_n); 409 vec.set(IDZDU , -b3_n); 410 vec.set(IQP , b5_n); 411 412 // Update trv 413 trv.setSurface( srf.newPureSurface() ); 414 trv.setVector(vec); 415 416 // set new direction of the track 417 if(sign_dun == 1) trv.setForward(); 418 if(sign_dun == -1) trv.setBackward(); 419 //cout<<trv<<endl; 420 // Calculate sT. 421 double crv = B*b5; 422 double dv = b1_n - b1; 423 double dxy = Math.sqrt( du*du + dv*dv ); 424 double arg = 0.5*crv*dxy; 425 double dst = dxy*TRFMath.asinrat(arg); 426 427 // Calculate s. 428 double dz = b2_n - b2; 429 Assert.assertTrue( flag_forward==1 || flag_forward==-1 ); 430 double ds = ((double) flag_forward)*Math.sqrt(dst*dst+dz*dz); 431 432 // Set the return status. 433 pstat.setPathDistance(ds); 434 //(flag_forward==1)?pstat.set_forward():pstat.set_backward(); 435 436 // exit now if user did not ask for error matrix. 437 if ( deriv == null ) return pstat; 438 439 // du_d0 440 441 double du_db1_0 = - sinphi_u; 442 443 // db1_d0 444 445 double db1_db1_0 = cosphi_u; 446 447 // db3_d0 448 449 double db3_db3_0 = a_hat_u2; 450 451 // db4_d0 452 453 double db4_db3_0 = b4_0*sinphi_u*a_hat_u2; 454 double db4_db4_0 = a_hat_u; 455 456 // dr_d 457 458 double dr_db3 = r*b3*b4*b4/(b3_hat2*b34_hat2); 459 double dr_db4 = -r*b4/b34_hat2; 460 double dr_db5 = -r/b5; 461 462 // dcosphi_d 463 464 double dcosphi_db3 = - sign_du/b3_hat - cosphi*b3/b3_hat2; 465 466 // dsinphi_d 467 468 double dsinphi_db3 = - sinphi*b3/b3_hat2; 469 470 // dcat_d 471 472 double dcat_db3 = norm/cat*(du/(r*r)*dr_db3 + dcosphi_db3 ); 473 double dcat_db4 = norm/cat* du/(r*r)*dr_db4; 474 double dcat_db5 = norm/cat* du/(r*r)*dr_db5; 475 double dcat_du = norm/(cat*r); 476 477 // dnorm_d 478 479 double dnorm_db3 = - du/(r*r)*dr_db3 - dcosphi_db3; 480 double dnorm_db4 = - du/(r*r)*dr_db4; 481 double dnorm_db5 = - du/(r*r)*dr_db5; 482 double dnorm_du = - 1./r; 483 484 // dcos_dphi_d 485 486 double dcos_dphi_db3 = - cosphi*dnorm_db3 - norm*dcosphi_db3 + 487 sign_cat*(sinphi*dcat_db3 + cat*dsinphi_db3); 488 double dcos_dphi_db4 = - cosphi*dnorm_db4 + sign_cat*sinphi*dcat_db4; 489 double dcos_dphi_db5 = - cosphi*dnorm_db5 + sign_cat*sinphi*dcat_db5; 490 double dcos_dphi_du = - cosphi*dnorm_du + sign_cat*sinphi*dcat_du; 491 492 // dsin_dphi_d 493 494 double dsin_dphi_db3 = sinphi*dnorm_db3 + norm*dsinphi_db3 + 495 sign_cat*(cosphi*dcat_db3 + cat*dcosphi_db3); 496 double dsin_dphi_db4 = sinphi*dnorm_db4 + sign_cat*cosphi*dcat_db4; 497 double dsin_dphi_db5 = sinphi*dnorm_db5 + sign_cat*cosphi*dcat_db5; 498 double dsin_dphi_du = sinphi*dnorm_du + sign_cat*cosphi*dcat_du; 499 500 // ddphi_d 501 502 double ddphi_db3; 503 double ddphi_db4; 504 double ddphi_db5; 505 double ddphi_du; 506 if(Math.abs(sin_dphi)>0.5) 507 { 508 ddphi_db3 = - dcos_dphi_db3/sin_dphi; 509 ddphi_db4 = - dcos_dphi_db4/sin_dphi; 510 ddphi_db5 = - dcos_dphi_db5/sin_dphi; 511 ddphi_du = - dcos_dphi_du /sin_dphi; 512 } 513 else 514 { 515 ddphi_db3 = dsin_dphi_db3/cos_dphi; 516 ddphi_db4 = dsin_dphi_db4/cos_dphi; 517 ddphi_db5 = dsin_dphi_db5/cos_dphi; 518 ddphi_du = dsin_dphi_du /cos_dphi; 519 } 520 521 // da_hat_dphi_d 522 523 double da_hat_dphi_db3 = - a_hat_dphi2* 524 (dcos_dphi_db3 - sin_dphi - b3*dsin_dphi_db3); 525 double da_hat_dphi_db4 = - a_hat_dphi2*(dcos_dphi_db4 - b3*dsin_dphi_db4); 526 double da_hat_dphi_db5 = - a_hat_dphi2*(dcos_dphi_db5 - b3*dsin_dphi_db5); 527 double da_hat_dphi_du = - a_hat_dphi2*(dcos_dphi_du - b3*dsin_dphi_du ); 528 529 // dc_hat_dphi_d 530 531 double dc_hat_dphi_db3 = b3*dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi; 532 double dc_hat_dphi_db4 = b3*dcos_dphi_db4 + dsin_dphi_db4; 533 double dc_hat_dphi_db5 = b3*dcos_dphi_db5 + dsin_dphi_db5; 534 double dc_hat_dphi_du = b3*dcos_dphi_du + dsin_dphi_du ; 535 536 // db1_n_d 537 538 double db1_n_db1 = 1; 539 double db1_n_db3 = (dr_db3*sinphi+r*dsinphi_db3)*(1-cos_dphi) 540 - rsinphi*dcos_dphi_db3 541 - dr_db3*cosphi*sin_dphi - r*dcosphi_db3*sin_dphi 542 - rcosphi*dsin_dphi_db3; 543 double db1_n_db4 = dr_db4*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db4 544 - dr_db4*cosphi*sin_dphi - rcosphi*dsin_dphi_db4; 545 double db1_n_db5 = dr_db5*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db5 546 - dr_db5*cosphi*sin_dphi - rcosphi*dsin_dphi_db5; 547 double db1_n_du = - rsinphi*dcos_dphi_du - rcosphi*dsin_dphi_du; 548 549 // db2_n_d 550 551 double db2_n_db2 = 1.; 552 double db2_n_db3 = 1./(b5*B)*b4*sign_du/b34_hat* 553 ( - dphi*b3/b34_hat2 + ddphi_db3); 554 double db2_n_db4 = 1./(b5*B)*sign_du/b34_hat* 555 ( - dphi*b4*b4/b34_hat2 + b4*ddphi_db4 + dphi ); 556 double db2_n_db5 = 1./(b5*B)*b4*sign_du/b34_hat*( ddphi_db5 - dphi/b5); 557 double db2_n_du = 1./(b5*B)*b4*sign_du/b34_hat*ddphi_du; 558 559 // db3_n_d 560 561 double db3_n_db3 = a_hat_dphi*dc_hat_dphi_db3 + da_hat_dphi_db3*c_hat_dphi; 562 double db3_n_db4 = a_hat_dphi*dc_hat_dphi_db4 + da_hat_dphi_db4*c_hat_dphi; 563 double db3_n_db5 = a_hat_dphi*dc_hat_dphi_db5 + da_hat_dphi_db5*c_hat_dphi; 564 double db3_n_du = a_hat_dphi*dc_hat_dphi_du + da_hat_dphi_du *c_hat_dphi; 565 566 // db4_n_d 567 568 double db4_n_db3 = b4*da_hat_dphi_db3; 569 double db4_n_db4 = b4*da_hat_dphi_db4 + a_hat_dphi; 570 double db4_n_db5 = b4*da_hat_dphi_db5; 571 double db4_n_du = b4*da_hat_dphi_du; 572 573 // db5_n_d 574 575 // db1_n_d0 576 577 double db1_n_db1_0 = db1_n_du * du_db1_0 + db1_n_db1*db1_db1_0; 578 double db1_n_db2_0 = 0.; 579 double db1_n_db3_0 = db1_n_db3*db3_db3_0 + db1_n_db4*db4_db3_0; 580 double db1_n_db4_0 = db1_n_db4*db4_db4_0; 581 double db1_n_db5_0 = db1_n_db5; 582 583 // db2_n_d0 584 585 double db2_n_db1_0 = db2_n_du *du_db1_0; 586 double db2_n_db2_0 = db2_n_db2; 587 double db2_n_db3_0 = db2_n_db3*db3_db3_0 + db2_n_db4*db4_db3_0; 588 double db2_n_db4_0 = db2_n_db4*db4_db4_0; 589 double db2_n_db5_0 = db2_n_db5; 590 591 // db3_n_d0 592 593 double db3_n_db1_0 = db3_n_du*du_db1_0; 594 double db3_n_db2_0 = 0.; 595 double db3_n_db3_0 = db3_n_db3*db3_db3_0 + db3_n_db4*db4_db3_0; 596 double db3_n_db4_0 = db3_n_db4*db4_db4_0; 597 double db3_n_db5_0 = db3_n_db5; 598 599 // db4_n_d0 600 601 double db4_n_db1_0 = db4_n_du*du_db1_0; 602 double db4_n_db2_0 = 0.; 603 double db4_n_db3_0 = db4_n_db3*db3_db3_0 + db4_n_db4*db4_db3_0; 604 double db4_n_db4_0 = db4_n_db4*db4_db4_0; 605 double db4_n_db5_0 = db4_n_db5; 606 607 // db5_n_d0 608 609 double db5_n_db1_0 = 0.; 610 double db5_n_db2_0 = 0.; 611 double db5_n_db3_0 = 0.; 612 double db5_n_db4_0 = 0.; 613 double db5_n_db5_0 = 1.; 614 615 // apply the rotation to the matrix db to return to the original c.f. 616 deriv.set(IV,IV, db2_n_db2_0); 617 deriv.set(IV,IZ, -db2_n_db1_0); 618 deriv.set(IV,IDVDU, db2_n_db4_0); 619 deriv.set(IV,IDZDU, -db2_n_db3_0); 620 deriv.set(IV,IQP, db2_n_db5_0); 621 deriv.set(IZ,IV, -db1_n_db2_0); 622 deriv.set(IZ,IZ, db1_n_db1_0); 623 deriv.set(IZ,IDVDU, -db1_n_db4_0); 624 deriv.set(IZ,IDZDU, db1_n_db3_0); 625 deriv.set(IZ,IQP, -db1_n_db5_0); 626 deriv.set(IDVDU,IV, db4_n_db2_0); 627 deriv.set(IDVDU,IZ, -db4_n_db1_0); 628 deriv.set(IDVDU,IDVDU, db4_n_db4_0); 629 deriv.set(IDVDU,IDZDU, -db4_n_db3_0); 630 deriv.set(IDVDU,IQP, db4_n_db5_0); 631 deriv.set(IDZDU,IV, -db3_n_db2_0); 632 deriv.set(IDZDU,IZ, db3_n_db1_0); 633 deriv.set(IDZDU,IDVDU, -db3_n_db4_0); 634 deriv.set(IDZDU,IDZDU, db3_n_db3_0); 635 deriv.set(IDZDU,IQP, -db3_n_db5_0); 636 deriv.set(IQP,IV, db5_n_db2_0); 637 deriv.set(IQP,IZ, -db5_n_db1_0); 638 deriv.set(IQP,IDVDU, db5_n_db4_0); 639 deriv.set(IQP,IDZDU, -db5_n_db3_0); 640 deriv.set(IQP,IQP, db5_n_db5_0); 641 642 return pstat; 643 644 } 645 646 }