anchor.diff:  
       else ! Crystal hit from above
   
  -      xp_tangent = asin_mb((rcurv*(one-cos_mb(cry_bend))-x)/sqrt(two*rcurv*(rcurv-x)*(one-cos_mb(cry_bend))+x**2))
  +      xp_tangent = asin_mb((c_rcurv*(one-cos_mb(cry_bend))-x)/sqrt(two*c_rcurv*(c_rcurv-x)*(one-cos_mb(cry_bend))+x**2))
         if(xp <= xp_tangent) then
           ! if it hits the crystal, calculate in which point and apply the
           ! transformation and drift to that point
           a_eq  = one + xp**2
  -        b_eq  = (two*xp)*(x - Rcurv)
  -        c_eq  = -(two*x)*Rcurv + x**2
  +        b_eq  = (two*xp)*(x - c_rcurv)
  +        c_eq  = -(two*x)*c_rcurv + x**2
           Delta = b_eq**2 - four*(a_eq*c_eq)
           S_int = (-b_eq-sqrt(Delta))/(two*a_eq)
   
  -        if(s_int < rcurv*sin_mb(cry_bend)) then
  +        if(s_int < c_rcurv*sin_mb(cry_bend)) then
             !  transform to a new ref system:shift and rotate
             x_int  = xp*s_int + x
             xp_int = xp
  
positive.code: 
        else
  
          s = rcurv*sin_mb(cry_length/rcurv)
          x = x+s*xp
          z = z+s*zp
  
  
negative.code:     y = y + yp*(s_length-s)
      return
  
    else if(x > c_xmax-c_alayer .and. x < c_xmax) then
      iProc = proc_AM
      call calc_ion_loss_cry(is,pc,s_length,dest)
      call move_am(is,nam,s_length,dest,dlyi(is),dlri(is),xp,yp,pc)
  