Skip to content

Commit

Permalink
Add 2 and 3 cmt matrices as eigen matrices
Browse files Browse the repository at this point in the history
  • Loading branch information
mattfidler committed Jul 19, 2024
1 parent 000ef9d commit 62abdb9
Showing 1 changed file with 145 additions and 2 deletions.
147 changes: 145 additions & 2 deletions src/compB.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ namespace stan {
Eigen::Matrix<T, Eigen::Dynamic, 2>
parTransB(const Eigen::Matrix<T, Eigen::Dynamic, 1>& p,
const int& ncmt,
const int& trans){
const int& trans) {
Eigen::Matrix<T, Eigen::Dynamic, 2> g(ncmt,3);
T btemp, ctemp, dtemp;
#define p1 p[0]
Expand Down Expand Up @@ -225,11 +225,154 @@ namespace stan {
#define k32 g(1, 1)
#define k24 g(2, 0)
#define k42 g(2, 1)

#define k10 g(0, 1)
#define k12 g(1, 0)
#define k21 g(1, 1)
#define k13 g(2, 0)
#define k31 g(2, 1)

template<class T>
int
solComp2Cpp(const Eigen::Matrix<T, Eigen::Dynamic, 2> g,
Eigen::Matrix<T, 2, 1>& L,
Eigen::Matrix<T, 2, 2>& C1,
Eigen::Matrix<T, 2, 2>& C2) {

T sum = k10 + k12 + k21;
T disc= sqrt(sum*sum - 4.0 * k10 * k21);

Eigen::Matrix<T, 2, 1> div;

L(0, 0) = 0.5 * (sum + disc);
L(1, 0) = 0.5 * (sum - disc);

div(0, 0) = L(1, 0) - L(0, 0);
div(1, 0) = L(0, 0) - L(1, 0);

if (div(0, 0)*div(1, 0) == 0) return 0;

C1(0, 0) = k21 - L(0, 0);
C1(0, 1) = k21 - L(1, 0);

C2(0, 0) = C2(0, 1) = k21;
C1(1, 0) = C1(1, 1) = k12;

T tmp = k10 + k12;

C2(1, 0) = tmp - L(0, 0);
C2(0, 1) = tmp - L(1, 0);

C1(0, 0) = C1(0, 0)/div(0, 0);
C1(1, 0) = C1(1, 0)/div(0, 0);

C2(0, 0) = C2(0, 0)/div(0, 0);
C2(1, 0) = C2(1, 0)/div(0, 0);

C1(0, 1) = C1(0, 1)/div(1, 0);
C1(1, 1) = C1(1, 1)/div(1, 0);

C2(0, 1) = C2(0, 1)/div(1, 0);
C2(1, 1) = C2(1, 1)/div(1, 0);

return 1;
}

template<class T>
int
solComp3Cpp(const Eigen::Matrix<T, Eigen::Dynamic, 2> g,
Eigen::Matrix<T, 3, 1>& L,
Eigen::Matrix<T, 3, 3>& C1,
Eigen::Matrix<T, 3, 3>& C2,
Eigen::Matrix<T, 3, 3>& C3) {

T A1 = k10 + k12 + k13 + k21 + k31;
T A2 = k10*k21 + k10*k31 + k12*k31 +
k13*k21 + k21*k31;
T A3 = k21*k31*k10;
T Q = (A1*A1 - 3.0*A2)/9.0;
T RQ = 2.0*sqrt(Q);
T R = (2.0*A1*A1*A1 - 9.0*A1*A2 + 27.0*A3)/54.0;
T M = Q*Q*Q - R*R;
if (M < 0) return 0;//stop("Error: Not real roots.")
T Th = acos(8.0*R/(RQ*RQ*RQ));
L(0, 0) = RQ*cos(Th/3.0) + A1/3.0;
L(1, 0) = RQ*cos((Th + M_2PI)/3.0) + A1/3.0;
L(2, 0) = RQ*cos((Th + 2*M_2PI)/3.0) + A1/3.0;

Eigen::Matrix<T, 3, 1> D;
D(0, 0) = (L(1, 0) - L(0, 0))*(L(2, 0) - L(0, 0));
D(1, 0) = (L(0, 0) - L(1, 0))*(L(2, 0) - L(1, 0));
D(2, 0) = (L(0, 0) - L(2, 0))*(L(1, 0) - L(2, 0));
if (D(0, 0)*D(1, 0)*D(2, 0) == 0.0) return 0;

C1(0, 0) = (k21 - L(0, 0))*(k31 - L(0, 0));
C1(0, 1) = (k21 - L(1, 0))*(k31 - L(1, 0));
C1(0, 2) = (k21 - L(2, 0))*(k31 - L(2, 0));

C2(0, 0) = (k21)*(k31 - L(0, 0));
C2(0, 1) = (k21)*(k31 - L(1, 0));
C2(0, 2) = (k21)*(k31 - L(2, 0));

C3(0, 0) = (k31)*(k21 - L(0, 0));
C3(0, 1) = (k31)*(k21 - L(1, 0));
C3(0, 2) = (k31)*(k21 - L(2, 0));

C1(1, 0) = (k12)*(k31 - L(0, 0));
C1(1, 1) = (k12)*(k31 - L(1, 0));
C1(1, 2) = (k12)*(k31 - L(2, 0));

C2(1, 0) = (k10 + k12 + k13 - L(0, 0))*(k31 - L(0, 0)) - (k31)*(k13);
C2(1, 1) = (k10 + k12 + k13 - L(1, 0))*(k31 - L(1, 0)) - (k31)*(k13);
C2(1, 2) = (k10 + k12 + k13 - L(2, 0))*(k31 - L(2, 0)) - (k31)*(k13);

C3(1, 0) = C3(1, 1) = C3(1, 2) = (k12)*(k31);

C1(2, 0) = (k13)*(k21 - L(0, 0));
C1(2, 1) = (k13)*(k21 - L(1, 0));
C1(2, 2) = (k13)*(k21 - L(2, 0));

C2(2, 0) = C2(2, 1) = C2(2, 2) = (k21)*(k13);

C3(2, 0) = (k10 + k12 + k13 - L(0, 0))*(k21 - L(0, 0)) - (k21)*(k12);
C3(2, 1) = (k10 + k12 + k13 - L(1, 0))*(k21 - L(1, 0)) - (k21)*(k12);
C3(2, 2) = (k10 + k12 + k13 - L(2, 0))*(k21 - L(2, 0)) - (k21)*(k12);

C1(0, 0) = C1(0, 0)/D(0, 0);
C1(1, 0) = C1(1, 0)/D(0, 0);
C1(2, 0) = C1(2, 0)/D(0, 0);

C2(0, 0) = C2(0, 0)/D(0, 0);
C2(1, 0) = C2(1, 0)/D(0, 0);
C2(2, 0) = C2(2, 0)/D(0, 0);

C3(0, 0) = C3(0, 0)/D(0, 0);
C3(1, 0) = C3(1, 0)/D(0, 0);
C3(2, 0) = C3(2, 0)/D(0, 0);

C1(0, 1) = C1(0, 1)/D(1, 0);
C1(1, 1) = C1(1, 1)/D(1, 0);
C1(2, 1) = C1(2, 1)/D(1, 0);

C2(0, 1) = C2(0, 1)/D(1, 0);
C2(1, 1) = C2(1, 1)/D(1, 0);
C2(2, 1) = C2(2, 1)/D(1, 0);

C3(0, 1) = C3(0, 1)/D(1, 0);
C3(1, 1) = C3(1, 1)/D(1, 0);
C3(2, 1) = C3(2, 1)/D(1, 0);

C1(0, 2) = C1(0, 2)/D(2, 0);
C1(1, 2) = C1(1, 2)/D(2, 0);
C1(2, 2) = C1(2, 2)/D(2, 0);

C2(0, 2) = C2(0, 2)/D(2, 0);
C2(1, 2) = C2(1, 2)/D(2, 0);
C2(2, 2) = C2(2, 2)/D(2, 0);

C3(0, 2) = C3(0, 2)/D(2, 0);
C3(1, 2) = C3(1, 2)/D(2, 0);
C3(2, 2) = C3(2, 2)/D(2, 0);

return 1;
}
}

0 comments on commit 62abdb9

Please sign in to comment.