Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F76684874
gradientgridCPU.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
Fri, Aug 9, 20:44
Size
3 KB
Mime Type
text/x-c
Expires
Sun, Aug 11, 20:44 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
19750993
Attached To
R1448 Lenstool-HPC
gradientgridCPU.cpp
View Options
#include "gradientgridCPU.hpp"
void
gradient_grid_CPU
(
double
*
grid_grad_x
,
double
*
grid_grad_y
,
const
struct
grid_param
*
frame
,
const
struct
Potential_SOA
*
lens
,
int
*
Nlens
){
gradient_grid_piemd_CPU
(
grid_grad_x
,
grid_grad_y
,
frame
,
Nlens
[
1
],
lens
[
1
].
position_x
,
lens
[
1
].
position_y
,
lens
[
1
].
b0
,
lens
[
1
].
ellipticity_angle
,
lens
[
1
].
ellipticity_potential
,
lens
[
1
].
rcore
,
lens
[
1
].
rcut
);
}
void
gradient_grid_piemd_CPU
(
double
*
grid_grad_x
,
double
*
grid_grad_y
,
const
struct
grid_param
*
frame
,
int
Nlens
,
double
*
lens_x
,
double
*
lens_y
,
double
*
lens_b0
,
double
*
lens_angle
,
double
*
lens_epot
,
double
*
lens_rcore
,
double
*
lens_rcut
){
int
bid
=
0
;
// index of the block (and of the set of images)
int
tid
=
0
;
// index of the thread within the block
double
dx
,
dy
,
x_pos
,
y_pos
;
//pixelsize
int
grid_dim
,
index
;
struct
point
true_coord
,
true_coord_rotation
;
complex
zis
;
dx
=
(
frame
->
xmax
-
frame
->
xmin
)
/
(
frame
->
nbgridcells
-
1
);
dy
=
(
frame
->
ymax
-
frame
->
ymin
)
/
(
frame
->
nbgridcells
-
1
);
grid_dim
=
(
frame
->
nbgridcells
);
index
=
bid
;
//while(index < grid_dim*grid_dim){
while
(
index
<
grid_dim
*
grid_dim
){
grid_grad_x
[
index
]
=
0.
;
grid_grad_y
[
index
]
=
0.
;
x_pos
=
frame
->
xmin
+
(
index
/
grid_dim
)
*
dx
;
y_pos
=
frame
->
ymin
+
(
index
%
grid_dim
)
*
dy
;
for
(
int
iterator
=
0
;
iterator
<
Nlens
;
iterator
++
){
true_coord
.
x
=
x_pos
-
lens_x
[
iterator
];
true_coord
.
y
=
y_pos
-
lens_y
[
iterator
];
// Change the origin of the coordinate system to the center of the clump
true_coord_rotation
=
rotateCoordinateSystem
(
true_coord
,
lens_angle
[
iterator
]);
double
x
=
true_coord_rotation
.
x
;
double
y
=
true_coord_rotation
.
y
;
double
eps
=
lens_epot
[
iterator
];
double
rc
=
lens_rcore
[
iterator
];
double
sqe
=
sqrt
(
eps
);
//
double
cx1
=
(
1.
-
eps
)
/
(
1.
+
eps
);
double
cxro
=
(
1.
+
eps
)
*
(
1.
+
eps
);
double
cyro
=
(
1.
-
eps
)
*
(
1.
-
eps
);
//
double
rem2
=
x
*
x
/
cxro
+
y
*
y
/
cyro
;
complex
zci
,
znum
,
zden
,
zres
;
double
norm
;
//
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) )
//
zis
.
re
=
zres
.
re
;
zis
.
im
=
zres
.
im
;
grid_grad_x
[
index
]
+=
lens_b0
[
iterator
]
*
zis
.
re
;
grid_grad_y
[
index
]
+=
lens_b0
[
iterator
]
*
zis
.
im
;
//printf("%d %f %f %f %f %f %f %f %f \n", index, x_pos, y_pos,true_coord.x,true_coord.y,grid_grad_x[index],grid_grad_y[index],lens_b0[iterator]*zis.re,lens_b0[iterator]*zis.im);
}
//printf("%d %d \n", index,grid_dim );
bid
+=
1
;
index
=
bid
*
1
+
tid
;
}
}
static
struct
point
rotateCoordinateSystem
(
struct
point
P
,
double
theta
)
{
struct
point
Q
;
Q
.
x
=
P
.
x
*
cos
(
theta
)
+
P
.
y
*
sin
(
theta
);
Q
.
y
=
P
.
y
*
cos
(
theta
)
-
P
.
x
*
sin
(
theta
);
return
(
Q
);
}
Event Timeline
Log In to Comment