/// The attractive part of the potential.
The attractive part of the potential.
The attractive part of the potential .
fn attractive_part(params: &Params, distsq: f64, r: V3) -> (f64, V3, M33) {
    debug_assert_eq!(distsq, r.sqnorm());
    // beta = r^2 / z0^2
    let z0sq = params.z0 * params.z0;
    let beta_d_r = 2.0 / z0sq * r;
    let beta_dd_r_r = 2.0 / z0sq; // times the identity matrix
    let beta_inv = z0sq / distsq;
    // -A * beta ** -3
    let value = -params.A * (beta_inv * beta_inv * beta_inv);
    let d_beta = (-3.0) * beta_inv * value;
    let dd_beta_beta = (-4.0) * beta_inv * d_beta;
    dbg!("VA: {:.9}", value);
    let d_r = d_beta * beta_d_r;
    let dd_r_r = d_beta * beta_dd_r_r * M33::eye() + outer(beta_d_r, beta_d_r) * dd_beta_beta;
    (value, d_r, dd_r_r)
}
fn attractive_part ( params : & Params , distsq : f64 , r : V3 ) -> ( f64 , V3 , M33 ) { debug_assert_eq ! ( distsq r sqnorm ( ) ) ; // beta = r^2 / z0^2 let z0sq = params . z0 * params . z0 ; let beta_d_r = 2.0 / z0sq * r ; let beta_dd_r_r = 2.0 / z0sq ; // times the identity matrix let beta_inv = z0sq / distsq ; // -A * beta ** -3 let value = - params . A * ( beta_inv * beta_inv * beta_inv ) ; let d_beta = ( - 3.0 ) * beta_inv * value ; let dd_beta_beta = ( - 4.0 ) * beta_inv * d_beta ; dbg ! ( " " value ) ; let d_r = d_beta * beta_d_r ; let dd_r_r = d_beta * beta_dd_r_r * M33 :: eye ( ) + outer ( beta_d_r , beta_d_r ) * dd_beta_beta ; ( value , d_r , dd_r_r ) }
The attractive part of the potential.
The attractive part of the potential .
params
&Params
distsq
f64
r
V3