Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F92462534
GradientTest.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Nov 20, 12:42
Size
9 KB
Mime Type
text/x-c
Expires
Fri, Nov 22, 12:42 (2 d)
Engine
blob
Format
Raw Data
Handle
22441421
Attached To
R1448 Lenstool-HPC
GradientTest.cpp
View Options
#include <iostream>
#include <string.h>
#include "structure.h"
#include <math.h>
#include <sys/time.h>
#include <fstream>
#include <../GradientBenchmark/Grad.h>
complex piemd_1derivatives(double x, double y, double eps, double rc);
void rotatecoordinatestest();
void piemdtest();
void sistest();
int main()
{
//Rotation of coordinates
rotatecoordinatestest();
//Tests for PIEMD potential
piemdtest();
//Tests for SIS potential
sistest();
}
void rotatecoordinatestest(){
//Init
point P;
P.x = P.y = 2;
double theta = 0.3; //rad
point coord = rotateCoordinateSystem(P,theta);
if (coord.x != P.x*cos(theta) + P.y*sin(theta)){
std::cerr << "Rotatecoordinatestest X: Theoretically we should have "
<< P.x*cos(theta) + P.y*sin(theta) << " and "
<< coord.x << std::endl;
exit(0);
}
if (coord.y != P.x*cos(theta) - P.y*sin(theta)){
std::cerr << "Rotatecoordinatestest Y: Theoretically we should have "
<< P.x*cos(theta) - P.y*sin(theta) << " and "
<< coord.y << std::endl;
exit(0);
}
}
void piemdtest() {
//////Setting the problem
int big(1);
runmode_param runmode;
point image;
runmode.nhalos = big;
image.x = image.y = 2;
/////// Test for 1 PIEMD mass distribution:
//1: module_readParameters_calculatePotentialparameter
//2: piemd_1derivatives_ci05
//3: gradient
//1: module_readParameters_calculatePotentialparameter
///////theoretical b0, ellipticity calculation
//Init
double b0(0), ellipticity_potential(0);
double vdisp(1.), ellipticity(0.11);
//Computation
ellipticity_potential = (1. - sqrt(1 - ellipticity * ellipticity))
/ ellipticity;
b0 = 6. * pi_c2 * vdisp * vdisp;
///////numerical b0, ellipticity_potential calculation
//Init
Potential* ilens;
Potential lens[big];
for (int i = 0; i < big; ++i) {
ilens = &lens[i];
ilens->position.x = ilens->position.y = -1.;
ilens->type = 8;
ilens->ellipticity = 0.11;
ilens->ellipticity_potential = 0.;
ilens->ellipticity_angle = 0.4;
ilens->vdisp = 1.;
ilens->rcut = 5.;
ilens->rcore = 1;
ilens->weight = 0;
ilens->rscale = 0;
ilens->exponent = 0;
ilens->alpha = 0.;
ilens->einasto_kappacritic = 0;
ilens->z = 0.4;
//Computation
module_readParameters_calculatePotentialparameter(ilens);
}
if (b0 != lens[0].b0) {
std::cerr << "PIEMD b0: Theoretically we should have " << b0 << " and "
<< lens[0].b0 << std::endl;
exit(0);
}
if (ellipticity_potential != lens[0].ellipticity_potential) {
std::cerr << "PIEMD ellipticity_potential: Theoretically we should have "
<< ellipticity_potential << " and "
<< lens[0].ellipticity_potential << std::endl;
exit(0);
}
//2: piemd_1derivatives_ci05
//Init
complex ztheo, zcomput;
///////theoretical piemd
ztheo = piemd_1derivatives(image.x, image.y, ellipticity_potential, 1);
///////numerical piemd
zcomput = piemd_1derivatives_ci05(image.x, image.y, ellipticity_potential,1);
if (ztheo.re != zcomput.re) {
std::cerr << "PIEMD z_piemd.re: Theoretically we should have " << ztheo.re
<< " and " << zcomput.re << std::endl;
exit(0);
}
if (ztheo.im != zcomput.im) {
std::cerr << "PIEMD z_piemd.im: Theoretically we should have " << ztheo.im
<< " and " << zcomput.im << std::endl;
exit(0);
}
//3: gradient
///////theoretical
//Init
point gradtheo, grad, true_coord, true_coord_rot;
grad.x = grad.y = 0;
//Computation
true_coord.x = image.x + 1; // Change the origin of the coordinate system to the center of the clump
true_coord.y = image.y + 1;
true_coord_rot = rotateCoordinateSystem(true_coord,0.4);
ztheo = piemd_1derivatives(true_coord_rot.x, true_coord_rot.y, ellipticity_potential, 1);
gradtheo.x = b0 * ztheo.re;
gradtheo.y = b0 * ztheo.im;
///////Computation
//Init SoA
PotentialSet lenses;
lenses.type = new int[big];
lenses.x = new double[big];
lenses.y = new double[big];
lenses.b0 = new double[big];
lenses.ellipticity_angle = new double[big];
lenses.ellipticity = new double[big];
lenses.ellipticity_potential = new double[big];
lenses.rcore = new double[big];
lenses.rcut = new double[big];
lenses.z = new double[big];
for (int i = 0; i < big; ++i) {
lenses.type[i] = lens[i].type;
lenses.x[i] = lens[i].position.x;
lenses.y[i] = lens[i].position.y;
lenses.b0[i] = lens[i].b0;
lenses.ellipticity_angle[i] = lens[i].ellipticity_angle;
lenses.ellipticity[i] = lens[i].ellipticity;
lenses.ellipticity_potential[i] = lens[i].ellipticity_potential;
lenses.rcore[i] = lens[i].rcore;
lenses.rcut[i] = lens[i].rcut;
lenses.z[i] = lens[i].z;
}
// Computation
grad = module_potentialDerivatives_totalGradient(&runmode, &image, &lenses);
if (gradtheo.x != grad.x) {
std::cerr << "PIEMD grad.x: Theoretically we should have " << gradtheo.x
<< " and " << grad.x << std::endl;
exit(0);
}
if (gradtheo.y != grad.y) {
std::cerr << "PIEMD grad.y: Theoretically we should have " << gradtheo.y
<< " and " << grad.y << std::endl;
exit(0);
}
std::cout << "PIEMD Everythings fine, No test failed " << std::endl;
}
void sistest() {
//////Setting the problem
int big(1);
runmode_param runmode;
point image;
runmode.nhalos = big;
image.x = image.y = 2;
/////// Test for 1 SIS mass distribution:
//1: module_readParameters_calculatePotentialparameter
//2: gradient
//1: module_readParameters_calculatePotentialparameter
///////theoretical b0, ellipticity calculation
//Init
double b0(0), ellipticity_potential(0);
double vdisp(1.), ellipticity(0.11);
//Computation
ellipticity_potential = ellipticity/3;
b0 = 4. * pi_c2 * vdisp * vdisp;
///////numerical b0, ellipticity_potential calculation
//Init
Potential* ilens;
Potential lens[big];
for (int i = 0; i < big; ++i) {
ilens = &lens[i];
ilens->position.x = ilens->position.y = -1.;
ilens->type = 5;
ilens->ellipticity = 0.11;
ilens->ellipticity_potential = 0.;
ilens->ellipticity_angle = 0.4;
ilens->vdisp = 1.;
ilens->rcut = 5.;
ilens->rcore = 1;
ilens->weight = 0;
ilens->rscale = 0;
ilens->exponent = 0;
ilens->alpha = 0.;
ilens->einasto_kappacritic = 0;
ilens->z = 0.4;
//Computation
module_readParameters_calculatePotentialparameter(ilens);
}
if (b0 != lens[0].b0) {
std::cerr << "SIS b0: Theoretically we should have " << b0 << " and "
<< lens[0].b0 << std::endl;
exit(0);
}
if (ellipticity_potential != lens[0].ellipticity_potential) {
std::cerr << "SIS ellipticity_potential: Theoretically we should have "
<< ellipticity_potential << " and "
<< lens[0].ellipticity_potential << std::endl;
exit(0);
}
//2: gradient
///////theoretical
//Init
double R;
point gradtheo, grad, true_coord, true_coord_rot;
grad.x = grad.y = 0;
//Computation
true_coord.x = image.x + 1; // Change the origin of the coordinate system to the center of the clump
true_coord.y = image.y + 1;
true_coord_rot = rotateCoordinateSystem(true_coord,0.4);
R=sqrt(true_coord_rot.x*true_coord_rot.x*(1-ellipticity/3.)+true_coord_rot.y*true_coord_rot.y*(1+ellipticity/3.));
gradtheo.x=(1-ellipticity/3.)*b0*true_coord_rot.x/(R);
gradtheo.y=(1+ellipticity/3.)*b0*true_coord_rot.y/(R);
///////Computation
//Init SoA
PotentialSet lenses;
lenses.type = new int[big];
lenses.x = new double[big];
lenses.y = new double[big];
lenses.b0 = new double[big];
lenses.ellipticity_angle = new double[big];
lenses.ellipticity = new double[big];
lenses.ellipticity_potential = new double[big];
lenses.rcore = new double[big];
lenses.rcut = new double[big];
lenses.z = new double[big];
for (int i = 0; i < big; ++i) {
lenses.type[i] = lens[i].type;
lenses.x[i] = lens[i].position.x;
lenses.y[i] = lens[i].position.y;
lenses.b0[i] = lens[i].b0;
lenses.ellipticity_angle[i] = lens[i].ellipticity_angle;
lenses.ellipticity[i] = lens[i].ellipticity;
lenses.ellipticity_potential[i] = lens[i].ellipticity_potential;
lenses.rcore[i] = lens[i].rcore;
lenses.rcut[i] = lens[i].rcut;
lenses.z[i] = lens[i].z;
}
// Computation
grad = module_potentialDerivatives_totalGradient(&runmode, &image, &lenses);
if (gradtheo.x != grad.x) {
std::cerr << "SIS grad.x: Theoretically we should have " << gradtheo.x
<< " and " << grad.x << std::endl;
exit(0);
}
if (gradtheo.y != grad.y) {
std::cerr << "SIS grad.y: Theoretically we should have " << gradtheo.y
<< " and " << grad.y << std::endl;
exit(0);
}
std::cout << "SIS Everythings fine, No test failed " << std::endl;
}
complex piemd_1derivatives(double x, double y, double eps, double rc)
{
double sqe, cx1, cxro, cyro, rem2;
complex zci, znum, zden, zis, zres;
double norm;
sqe = sqrt(eps);
cx1 = (1. - eps) / (1. + eps);
cxro = (1. + eps) * (1. + eps);
cyro = (1. - eps) * (1. - eps);
rem2 = x * x / cxro + y * y / cyro;
/*zci=cpx(0.,-0.5*(1.-eps*eps)/sqe);
znum=cpx(cx1*x,(2.*sqe*sqrt(rc*rc+rem2)-y/cx1));
zden=cpx(x,(2.*rc*sqe-y));
zis=pcpx(zci,lncpx(dcpx(znum,zden)));
zres=pcpxflt(zis,b0);*/
// --> optimized code
zci.re = 0;
zci.im = -0.5 * (1. - eps * eps) / sqe;
znum.re = cx1 * x;
znum.im = 2.*sqe * sqrt(rc * rc + rem2) - y / cx1;
zden.re = x;
zden.im = 2.*rc * sqe - y;
norm = zden.re * zden.re + zden.im * zden.im; // zis = znum/zden
zis.re = (znum.re * zden.re + znum.im * zden.im) / norm;
zis.im = (znum.im * zden.re - znum.re * zden.im) / norm;
norm = zis.re;
zis.re = log(sqrt(norm * norm + zis.im * zis.im)); // ln(zis) = ln(|zis|)+i.Arg(zis)
zis.im = atan2(zis.im, norm);
// norm = zis.re;
zres.re = zci.re * zis.re - zci.im * zis.im; // Re( zci*ln(zis) )
zres.im = zci.im * zis.re + zis.im * zci.re; // Im( zci*ln(zis) )
//zres.re = zis.re*b0;
//zres.im = zis.im*b0;
return(zres);
}
Event Timeline
Log In to Comment