anchor.diff:        r2b_d=two*(sigman2(2,imbb(i))-sigman2(1,imbb(i)))
         rb_d=sqrt(r2b_d)
         rkb_d=(strack(i)*pisqrt)/rb_d                            !hr03
  -
         xrb_d=abs(crkveb_d)/rb_d
         zrb_d=abs(cikveb_d)/rb_d
         call errf(zrb_d,xrb_d,crzb_d,crxb_d)
  
positive.code:       xbb_d=sigmanq(2,imbb(i))*xrb_d
        zbb_d=sigmanq(1,imbb(i))*zrb_d
        call errf(zbb_d,xbb_d,cbzb_d,cbxb_d)
      else if(ibtyp.eq.1) then
        r2b_d=two*(sigman2(2,imbb(i))-sigman2(1,imbb(i)))
        rb_d=sqrt(r2b_d)
  
negative.code: if(nbeaux(imbb(i)).eq.2) then
    ktrack(i)=42
    if(ibeco.eq.1) then
      r2b_d=two*(sigman2(1,imbb(i))-sigman2(2,imbb(i)))
      rb_d=sqrt(r2b_d)
      rkb_d=(strack(i)*pisqrt)/rb_d                            
      xrb_d=abs(crkveb_d)/rb_d
      zrb_d=abs(cikveb_d)/rb_d
      if(ibtyp.eq.0) then
  
  
        call errf(xrb_d,zrb_d,crxb_d,crzb_d)
  
  
        tkb_d=(crkveb_d**2/sigman2(1,imbb(i))+cikveb_d**2/sigman2(2,imbb(i)))*half                              
        xbb_d=sigmanq(2,imbb(i))*xrb_d
        zbb_d=sigmanq(1,imbb(i))*zrb_d
  
        call errf(xbb_d,zbb_d,cbxb_d,cbzb_d)
  
      else if(ibtyp.eq.1) then
  
  
        tkb_d=(crkveb_d**2/sigman2(1,imbb(i))+cikveb_d**2/sigman2(2,imbb(i)))*half                              
        xbb_d=sigmanq(2,imbb(i))*xrb_d
        zbb_d=sigmanq(1,imbb(i))*zrb_d
      endif
    endif
    beamoff(4,imbb(i))=(rkb_d*(crzb_d-exp_mb(-one*tkb_d)*cbzb_d))*sign(one,crkveb_d)
  