Skip to content

Working and correct CDM version #13

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
gracebato wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Working and correct CDM version #13

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion models/cdm/cdm/CDM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,7 +106,7 @@ def initialize(self, application):
theta = record.theta
phi = record.phi
# form the projection vectors and store them
self.los[obs, 0] = sin(theta) * cos(phi)
self.los[obs, 0] = -sin(theta) * cos(phi)
self.los[obs, 1] = sin(theta) * sin(phi)
self.los[obs, 2] = cos(theta)

Expand Down
51 changes: 44 additions & 7 deletions models/cdm/lib/libcdm/Source.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,6 +72,9 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
// clean up the resulting matrix
gsl_matrix_set_zero(predicted);

// allocate storage for the displacement vectors; we reuse this for all samples
gsl_matrix * disp = gsl_matrix_alloc(_locations->size1, 3);

// go through all the samples
for (auto sample=0; sample<nSamples; ++sample) {
// unpack the parameters
Expand All @@ -90,28 +93,62 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
auto omegaY = gsl_matrix_get(&samples->matrix, sample, _omegaYIdx);
auto omegaZ = gsl_matrix_get(&samples->matrix, sample, _omegaZIdx);

// compute the displacements
cdm(sample, _locations, _los,
// // compute the displacements
// cdm(sample, _locations, _los,
// xSrc, ySrc, dSrc,
// aX, aY, aZ,
// omegaX, omegaY, omegaZ,
// openingSrc,
// _nu,
// predicted);

// // apply the location specific projection to LOS vector and dataset shift
// for (auto loc=0; loc<_locations->size1; ++loc) {
// // get the current value
// auto u = gsl_matrix_get(predicted, sample, loc);
// // find the shift that corresponds to this observation
// auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]);
// // and apply it to the projected displacement
// u -= shift;
// // save
// gsl_matrix_set(predicted, sample, loc, u);

cdm(_locations,
xSrc, ySrc, dSrc,
openingSrc,
aX, aY, aZ,
omegaX, omegaY, omegaZ,
openingSrc,
_nu,
predicted);

disp);
// apply the location specific projection to LOS vector and dataset shift
for (auto loc=0; loc<_locations->size1; ++loc) {
// get the current value
auto u = gsl_matrix_get(predicted, sample, loc);
// compute the components of the unit LOS vector
auto nx = gsl_matrix_get(_los, loc, 0);
auto ny = gsl_matrix_get(_los, loc, 1);
auto nz = gsl_matrix_get(_los, loc, 2);

// get the three components of the predicted displacement for this location
auto ux = gsl_matrix_get(disp, loc, 0);
auto uy = gsl_matrix_get(disp, loc, 1);
auto ud = gsl_matrix_get(disp, loc, 2);

// project
auto u = ux*nx + uy*ny + ud*nz;
// find the shift that corresponds to this observation
auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]);
// and apply it to the projected displacement
u -= shift;

// save
gsl_matrix_set(predicted, sample, loc, u);

}
}

// clean up
gsl_matrix_free(disp);

// all done
return;
}
Expand Down
165 changes: 120 additions & 45 deletions models/cdm/lib/libcdm/cdm.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,12 @@ namespace altar {

// local helpers
static void
RDdispSurf(int sample,
const gsl_matrix * locations, const gsl_matrix * los,
// RDdispSurf(int sample,
// const gsl_matrix * locations, const gsl_matrix * los,
// const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4,
// double opening, double nu,
// gsl_matrix * results);
RDdispSurf(const gsl_matrix * locations,
const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4,
double opening, double nu,
gsl_matrix * results);
Expand Down Expand Up @@ -76,27 +80,34 @@ namespace altar {
// definitions
void
altar::models::cdm::
cdm(int sample,
const gsl_matrix * locations, const gsl_matrix * los,
cdm(const gsl_matrix * locations,
double x, double y, double depth,
double opening,
double aX, double aY, double aZ,
double omegaX, double omegaY, double omegaZ,
double opening,
double nu,
gsl_matrix * predicted)
// cdm(int sample,
// const gsl_matrix * locations, const gsl_matrix * los,
// double x, double y, double depth,
// double aX, double aY, double aZ,
// double omegaX, double omegaY, double omegaZ,
// double opening,
// double nu,
// gsl_matrix * predicted)
{
// convert semi-axes to axes
aX *= 2;
aY *= 2;
aZ *= 2;

// short circuit the trivial case
if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) {
// no displacements
gsl_matrix_set_zero(predicted);
// all done
return;
}
// if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) {
// // no displacements
// gsl_matrix_set_zero(predicted);
// // all done
// return;
// }

// the axis specific coordinate matrices
mat_t Rx = {1., 0., 0.,
Expand Down Expand Up @@ -141,34 +152,84 @@ cdm(int sample,
Q1[2] > 0 || Q2[2] > 0 || Q3[2] > 0 || Q4[2] > 0 ||
R1[2] > 0 || R2[2] > 0 || R3[2] > 0 || R4[2] > 0) {
// complain...
throw std::domain_error("the CDM must be below the surface");
// throw std::domain_error("the CDM must be below the surface");
printf("WARNING: The CDM must be below the surface \n");
}

// // dispatch the various cases
// if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) {
// RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted);
// } else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) {
// RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted);
// } else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) {
// RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted);
// } else {
// RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted);
// RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted);
// RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted);
// }

// dispatch the various cases
if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) {
RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted);
} else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) {
RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted);
} else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) {
RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted);
// short circuit the trivial case
if (std::abs(aX) == 0 && std::abs(aY) == 0 && std::abs(aZ) == 0) {
// no displacements
printf("I am if\n");
gsl_matrix_set_zero(predicted);
}
else if (std::abs(aX) == 0 && std::abs(aY) > 0 && std::abs(aZ) > 0) {
printf("I am elseif1 \n");
RDdispSurf(locations, P1, P2, P3, P4, opening, nu, predicted);
} else if (std::abs(aX) > 0 && std::abs(aY) == 0 && std::abs(aZ) > 0) {
printf("I am elseif2 \n");
RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, predicted);
} else if (std::abs(aX) > 0 && std::abs(aY) > 0 && std::abs(aZ) == 0) {
printf("I am elseif3 \n");
RDdispSurf(locations, R1, R2, R3, R4, opening, nu, predicted);
} else {
RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted);
RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted);
RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted);

// allocate room for the partial results
gsl_matrix * P = gsl_matrix_alloc(locations->size1, 3);
gsl_matrix * Q = gsl_matrix_alloc(locations->size1, 3);
gsl_matrix * R = gsl_matrix_alloc(locations->size1, 3);

// compute
RDdispSurf(locations, P1, P2, P3, P4, opening, nu, P);
RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, Q);
RDdispSurf(locations, R1, R2, R3, R4, opening, nu, R);

// assemble
for (auto loc=0; loc<locations->size1; ++loc) {
for (auto axis=0; axis<3; ++axis) {
// combine
auto result =
gsl_matrix_get(P, loc, axis) +
gsl_matrix_get(Q, loc, axis) +
gsl_matrix_get(R, loc, axis);

// assign
gsl_matrix_set(predicted, loc, axis, result);
}
}
}


// all done
return;
}

// implementations
static void
altar::models::cdm::
RDdispSurf(int sample,
const gsl_matrix * locations, const gsl_matrix * los,
// RDdispSurf(int sample,
// const gsl_matrix * locations, const gsl_matrix * los,
// const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4,
// double opening, double nu,
// gsl_matrix * results) {
RDdispSurf(const gsl_matrix * locations,
const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4,
double opening, double nu,
gsl_matrix * results) {

// cross
auto V = cross(P2-P1, P4-P1);
auto b = opening * V/norm(V);
Expand All @@ -184,23 +245,30 @@ RDdispSurf(int sample,
auto u3 = AngSetupFSC(x,y, b, P3, P4, nu);
auto u4 = AngSetupFSC(x,y, b, P4, P1, nu);

// // assemble
// auto u = u1 + u2 + u3 + u4;
// // compute the unit LOS vector
// vec_t n = { gsl_matrix_get(los, loc, 0),
// gsl_matrix_get(los, loc, 1),
// gsl_matrix_get(los, loc, 2) };

// // project the displacement to the LOS
// auto uLOS = dot(u, n);
// // save by accumulating my contribution to the slot
// // N.B.: note the "+=": the general case call this function three times
// // get the current value
// auto current = gsl_matrix_get(results, sample, loc);
// // update
// current += uLOS;
// // save
// gsl_matrix_set(results, sample, loc, current);
// }

// assemble
auto u = u1 + u2 + u3 + u4;
// compute the unit LOS vector
vec_t n = { gsl_matrix_get(los, loc, 0),
gsl_matrix_get(los, loc, 1),
gsl_matrix_get(los, loc, 2) };

// project the displacement to the LOS
auto uLOS = dot(u, n);
// save by accumulating my contribution to the slot
// N.B.: note the "+=": the general case call this function three times
// get the current value
auto current = gsl_matrix_get(results, sample, loc);
// update
current += uLOS;
// save
gsl_matrix_set(results, sample, loc, current);
for (auto axis=0; axis<3; ++axis) {
auto u = u1[axis] + u2[axis] + u3[axis] + u4[axis];
gsl_matrix_set(results, loc, axis, u);
}
}

// all done
Expand All @@ -215,7 +283,8 @@ AngSetupFSC(double x, double y,
double nu) {
vec_t SideVec = PB - PA;
vec_t eZ = {0, 0, 1};
auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec));
// auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec));
auto beta = std::acos(dot(-SideVec, eZ) / norm(SideVec));

if (std::abs(beta) < eps || std::abs(pi - beta) < eps) {
return { 0,0,0 };
Expand All @@ -239,14 +308,20 @@ AngSetupFSC(double x, double y,

vec_t vA, vB;
// distinguish the two configurations
if (beta*adcsA[0] > 0) {
// if (beta*adcsA[0] > 0) {
if (beta*adcsA[0] >= 0) {
// configuration I
vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]);
vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]);
// vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]);
// vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]);
vA = AngDisDispSurf(adcsA, -pi+beta, bADCS, nu, -PA[2]);
vB = AngDisDispSurf(adcsB, -pi+beta, bADCS, nu, -PB[2]);

} else {
// configuration II
vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]);
vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]);
// vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]);
// vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]);
vA = AngDisDispSurf(adcsA, beta, bADCS, nu, -PA[2]);
vB = AngDisDispSurf(adcsB, beta, bADCS, nu, -PB[2]);
}

vec_t v = xform(transpose(A), vB - vA);
Expand Down
14 changes: 11 additions & 3 deletions models/cdm/lib/libcdm/cdm.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,23 @@
namespace altar {
namespace models {
namespace cdm {
void cdm(int sample,
const gsl_matrix * locations, const gsl_matrix * los,
void cdm(const gsl_matrix * locations,
double X0, double Y0, double depth,
double opening,
double ax, double ay, double az,
double omegaX, double omegaY, double omegaZ,
double opening,
double nu,
gsl_matrix * predicted
);
// void cdm(int sample,
// const gsl_matrix * locations, const gsl_matrix * los,
// double X0, double Y0, double depth,
// double ax, double ay, double az,
// double omegaX, double omegaY, double omegaZ,
// double opening,
// double nu,
// gsl_matrix * predicted
// );
}
}
}
Expand Down