Advertisement
Alehandreus

Langevin Preconditioning Matrix

Feb 26th, 2024 (edited)
811
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 1.17 KB | Science | 0 0
  1. MatrixXd precond_3dim_new(const VectorXd &p, const VectorXd &params) {
  2.     int n = p.size();
  3.  
  4.     double r1 = params(0);
  5.     double r2 = params(1);
  6.     double r3 = params(2);
  7.     double k  = params(3);
  8.  
  9.     double x = p(0), y = p(1), z = p(2);
  10.     double phi = atan2(y, x);
  11.     double theta = phi * k;
  12.    
  13.  
  14.     // ==== e1, e2, e3 ====
  15.     VectorXd e1(3);    
  16.     e1(0) = -x * sin(phi) / cos(phi) + k * z * cos(phi);
  17.     e1(1) = x * cos(phi) /  cos(phi) + k * z * sin(phi);
  18.     e1(2) = -k * x / cos(phi) + k * r1;
  19.     e1.normalize();
  20.  
  21.     VectorXd e2(3);
  22.     e2(0) = cos(theta) * cos(phi);
  23.     e2(1) = cos(theta) * sin(phi);
  24.     e2(2) = -sin(theta);
  25.  
  26.     VectorXd e3(3);
  27.     e3(0) = sin(theta) * cos(phi);
  28.     e3(1) = sin(theta) * sin(phi);
  29.     e3(2) = cos(theta);
  30.  
  31.    
  32.     // ==== e1, e2, e3 coefs ====
  33.     double de1 = sqrt(
  34.         2 * r3 * (r1 * r1 + k*k * r2*r2) / sqrt(r1*r1 + k*k*k*k * r2*r2)
  35.     ); // опасно r2 <-> r3
  36.     double de2 = r2;
  37.     double de3 = r3;
  38.  
  39.  
  40.     // ==== create preconditioning matrix ====
  41.     MatrixXd res = MatrixXd::Zero(3, 3);
  42.     res.col(0) = e1 * de1;
  43.     res.col(1) = e2 * de2;
  44.     res.col(2) = e3 * de3;
  45.  
  46.     return res;
  47. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement