math.h File Reference

Header file for math.c. More...

#include "offload_acc_omp.h"
#include <math.h>
#include "ascot5.h"

Go to the source code of this file.

Macros
#define	math_degrad   0.0174532925199432957692
	One degree in radians.
#define	math_raddeg   57.295779513082320876798
	One radian in degrees.
#define	math_maxSimpsonDepth   20
	Maximum recursion depth for the simpson integral rule.
#define	math_copy(a, b)   do { a[0]=b[0];a[1]=b[1];a[2]=b[2]; } while(0)
	Copies elements of vector b to vector a.
#define	math_matcopy(a, b)
	Copies elements of matrix b to matrix a.
#define	math_dot(a, b)   (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
	Calculate dot product a[3] dot b[3].
#define	math_cross(a, b, c)
	Calculate cross product for 3D vectors c = a x b.
#define	math_scalar_triple_product(a, b, c)
	the mixed or triple product of three vectors (a cross b) dot c mathematica: a={a0,a1,a2};b={b0,b1,b2};c={c0,c1,c2};FullSimplify[Dot[Cross[a,b],c]] -(a2 b1 c0) + a1 b2 c0 + a2 b0 c1 - a0 b2 c1 - a1 b0 c2 + a0 b1 c2
#define	math_sumew(a, b)   do { a[0]=a[0]+b[0];a[1]=a[1]+b[1];a[2]=a[2]+b[2] } while(0)
	Elementwise vector sum a = a + b.
#define	math_prod(a, b)   do { a[0]=a[0]*b;a[1]=a[1]*b;a[2]=a[2]*b; } while(0)
	Multiply vector elements with a scalar a = b*a.
#define	math_matsumew(a, b)
	Elementwise matrix sum a = a + b.
#define	math_matprod(a, b)
	Multiply matrix a with scalar b: a = b*a.
#define	math_norm(a)   (sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]))
	Calculate norm of 3D vector a.
#define	math_normc(a1, a2, a3)   (sqrt(a1*a1+a2*a2+a3*a3))
	Calculate norm of 3D vector from its components a1, a2, a3.
#define	math_unit(a, b)
	Calculate unit vector b from a 3D vector a.
#define	math_xyz2rpz(xyz, rpz)
	Convert cartesian coordinates xyz to cylindrical coordinates rpz.
#define	math_rpz2xyz(rpz, xyz)
	Convert cylindrical coordinates rpz to cartesian coordinates xyz.
#define	math_vec_rpz2xyz(vrpz, vxyz, phi)
	Transform vector from cylindrical to cartesian basis: vrpz -> vxyz, phi is the toroidal angle in radians.
#define	math_vec_xyz2rpz(vxyz, vrpz, phi)
	Transform vector from cartesian to cylindrical basis: vxyz -> vrpz, phi is the toroidal angle in radians.
#define	math_determinant3x3(x1, x2, x3, y1, y2, y3, z1, z2, z3)
	Direct expansion of 3x3 matrix determinant.
#define	math_deg2rad(a)   (a * math_degrad)
	Convert degrees to radians.
#define	math_rad2deg(a)   (a * math_raddeg)
	Convert radians to degrees.

Functions
DECLARE_TARGET real	fmod (real x, real y)
	Compute the modulus of two real numbers.
DECLARE_TARGET_END DECLARE_TARGET_SIMD void	math_jac_rpz2xyz (real *rpz, real *xyz, real r, real phi)
	Convert a Jacobian from cylindrical to cartesian coordinates.
GPU_DECLARE_TARGET_SIMD void	math_jac_xyz2rpz (real *xyz, real *rpz, real r, real phi)
	Convert a Jacobian from cartesian to cylindrical coordinates.
DECLARE_TARGET_END GPU_DECLARE_TARGET_SIMD void	math_matmul (real *matA, real *matB, int d1, int d2, int d3, real *matC)
	Matrix multiplication.
DECLARE_TARGET_END DECLARE_TARGET_SIMD real	math_normal_rand ()
	Generate normally distributed random numbers.
DECLARE_TARGET_SIMD int	math_ipow (int a, int p)
	Calculate a^p where both a and p are integers (p >= 0)
double	math_simpson (double(*f)(double), double a, double b, double epsilon)
	Adaptive Simpsons rule for integral.
void	math_linspace (real *vec, real a, real b, int n)
	Generate linearly spaced vector.
DECLARE_TARGET_SIMD int	math_point_on_plane (real q[3], real t1[3], real t2[3], real t3[3])
	Find if a point is on a given plane.
DECLARE_TARGET_SIMD void	math_barycentric_coords_triangle (real AP[3], real AB[3], real AC[3], real n[3], real *s, real *t)
	Find barycentric coordinates for a given point.
void	math_uniquecount (int *in, int *unique, int *count, int n)
	Find unique numbers and their frequency in given array.
DECLARE_TARGET_SIMD real *	math_rsearch (const real key, const real *base, int num)
	Search for array element preceding a key value.
int	math_point_in_polygon (real r, real z, real *rv, real *zv, int n)
	Check if coordinates are within polygon.

Detailed Description

Header file for math.c.

Definition in file math.h.

Macro Definition Documentation

math_degrad

#define math_degrad   0.0174532925199432957692

One degree in radians.

Definition at line 13 of file math.h.

math_raddeg

#define math_raddeg   57.295779513082320876798

One radian in degrees.

Definition at line 15 of file math.h.

math_maxSimpsonDepth

#define math_maxSimpsonDepth   20

Maximum recursion depth for the simpson integral rule.

Definition at line 18 of file math.h.

math_copy

#define math_copy	(		a,
			b )   do { a[0]=b[0];a[1]=b[1];a[2]=b[2]; } while(0)

Copies elements of vector b to vector a.

Definition at line 21 of file math.h.

math_matcopy

#define math_matcopy	(		a,
			b )

Value:

        do { a[0]=b[0];a[1]=b[1];a[2]=b[2];a[3]=b[3]; \
        a[4]=b[4];a[5]=b[5];a[6]=b[6];a[7]=b[7];a[8]=b[8]; } while(0)

Copies elements of matrix b to matrix a.

Definition at line 24 of file math.h.

math_dot

#define math_dot	(		a,
			b )   (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])

Calculate dot product a[3] dot b[3].

Definition at line 28 of file math.h.

math_cross

#define math_cross	(		a,
			b,
			c )

Value:

        do { c[0]=a[1]*b[2]-a[2]*b[1];                \
        c[1]=a[2]*b[0]-a[0]*b[2];c[2]=a[0]*b[1]-a[1]*b[0]; } while(0)

Calculate cross product for 3D vectors c = a x b.

Definition at line 31 of file math.h.

math_scalar_triple_product

#define math_scalar_triple_product	(		a,
			b,
			c )

Value:

       ( - (a[2] * b[1] * c[0] ) \
         + (a[1] * b[2] * c[0] ) \
         + (a[2] * b[0] * c[1] ) \
         - (a[0] * b[2] * c[1] ) \
         - (a[1] * b[0] * c[2] ) \
         + (a[0] * b[1] * c[2] )     )

the mixed or triple product of three vectors (a cross b) dot c mathematica: a={a0,a1,a2};b={b0,b1,b2};c={c0,c1,c2};FullSimplify[Dot[Cross[a,b],c]] -(a2 b1 c0) + a1 b2 c0 + a2 b0 c1 - a0 b2 c1 - a1 b0 c2 + a0 b1 c2

Definition at line 39 of file math.h.

math_sumew

#define math_sumew	(		a,
			b )   do { a[0]=a[0]+b[0];a[1]=a[1]+b[1];a[2]=a[2]+b[2] } while(0)

Elementwise vector sum a = a + b.

Definition at line 48 of file math.h.

math_prod

#define math_prod	(		a,
			b )   do { a[0]=a[0]*b;a[1]=a[1]*b;a[2]=a[2]*b; } while(0)

Multiply vector elements with a scalar a = b*a.

Definition at line 51 of file math.h.

math_matsumew

#define math_matsumew	(		a,
			b )

Value:

        do { a[0]=a[0]+b[0];a[1]=a[1]+b[1];a[2]=a[2]+b[2]; \
        a[3]=a[3]+b[3];a[4]=a[4]+b[4];a[5]=a[5]+b[5];a[6]=a[6]+b[6];    \
        a[7]=a[7]+b[7];a[8]=a[8]+b[8]; } while(0)

Elementwise matrix sum a = a + b.

Definition at line 54 of file math.h.

math_matprod

#define math_matprod	(		a,
			b )

Value:

        do { a[0]=a[0]*b;a[1]=a[1]*b;a[2]=a[2]*b;  \
        a[3]=a[3]*b;a[4]=a[4]*b;a[5]=a[5]*b;a[6]=a[6]*b;a[7]=a[7]*b; \
        a[8]=a[8]*b } while(0)

Multiply matrix a with scalar b: a = b*a.

Definition at line 59 of file math.h.

math_norm

#define math_norm

(

a

)

(sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]))

Calculate norm of 3D vector a.

Definition at line 64 of file math.h.

math_normc

#define math_normc	(		a1,
			a2,
			a3 )   (sqrt(a1*a1+a2*a2+a3*a3))

Calculate norm of 3D vector from its components a1, a2, a3.

Definition at line 67 of file math.h.

math_unit

#define math_unit	(		a,
			b )

Value:

        do {real _n=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]); \
        b[0]=a[0]/_n;b[1]=a[1]/_n;b[2]=a[2]/_n; } while(0)

Calculate unit vector b from a 3D vector a.

Definition at line 70 of file math.h.

math_xyz2rpz

#define math_xyz2rpz	(		xyz,
			rpz )

Value:

        do { rpz[0]=sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]); \
        rpz[1]=atan2(xyz[1],xyz[0]);rpz[2]=xyz[2]; } while(0)

Convert cartesian coordinates xyz to cylindrical coordinates rpz.

Definition at line 74 of file math.h.

math_rpz2xyz

#define math_rpz2xyz	(		rpz,
			xyz )

Value:

        do { xyz[0]=rpz[0]*cos(rpz[1]); \
        xyz[1]=rpz[0]*sin(rpz[1]);xyz[2]=rpz[2]; } while(0)

Convert cylindrical coordinates rpz to cartesian coordinates xyz.

Definition at line 78 of file math.h.

math_vec_rpz2xyz

#define math_vec_rpz2xyz	(		vrpz,
			vxyz,
			phi )

Value:

        do {     \
        vxyz[0]=vrpz[0]*cos(phi)-vrpz[1]*sin(phi); \
        vxyz[1]=vrpz[0]*sin(phi)+vrpz[1]*cos(phi); \
        vxyz[2]=vrpz[2]; } while(0)

Transform vector from cylindrical to cartesian basis: vrpz -> vxyz, phi is the toroidal angle in radians.

Definition at line 83 of file math.h.

math_vec_xyz2rpz

#define math_vec_xyz2rpz	(		vxyz,
			vrpz,
			phi )

Value:

        do {      \
        vrpz[0]=vxyz[0]*cos(phi)+vxyz[1]*sin(phi);  \
        vrpz[1]=-vxyz[0]*sin(phi)+vxyz[1]*cos(phi); \
        vrpz[2]=vxyz[2]; } while(0)

Transform vector from cartesian to cylindrical basis: vxyz -> vrpz, phi is the toroidal angle in radians.

Definition at line 90 of file math.h.

math_determinant3x3

#define math_determinant3x3	(		x1,
			x2,
			x3,
			y1,
			y2,
			y3,
			z1,
			z2,
			z3 )

Value:

    (x1) * ( (y2)*(z3) - (y3)*(z2) ) + \
    (x2) * ( (y3)*(z1) - (y1)*(z3) ) + \
    (x3) * ( (y1)*(z2) - (y2)*(z1) )

Direct expansion of 3x3 matrix determinant.

Definition at line 97 of file math.h.

math_deg2rad

#define math_deg2rad

(

a

)

(a * math_degrad)

Convert degrees to radians.

Definition at line 107 of file math.h.

math_rad2deg

#define math_rad2deg

(

a

)

(a * math_raddeg)

Convert radians to degrees.

Definition at line 110 of file math.h.

Function Documentation

fmod()

DECLARE_TARGET real fmod	(	real	x,
		real	y )

Compute the modulus of two real numbers.

Parameters

x	The dividend
y	The divisor

Returns

The modulus (remainder) of x and y

Definition at line 22 of file math.c.

math_jac_rpz2xyz()

DECLARE_TARGET_END DECLARE_TARGET_SIMD void math_jac_rpz2xyz	(	real *	rpz,
		real *	xyz,
		real	r,
		real	phi )

Convert a Jacobian from cylindrical to cartesian coordinates.

This function converts a Jacobian located at angle phi and radius r from cylindrical to cartesian coordinates. The input Jacobian is an array with

[Br dBr/dr dBr/dphi dBr/dz Bphi dBphi/dr dBphi/dphi dBphi/dz Bz dBz/dr dBz/dphi dBz/dz]

an the output is

[Bx dBx/dx dBx/dy dBx/dz By dBy/dx dBy/dy dBy/dz Bz dBz/dx dBz/dy dBz/dz].

Parameters

rpz	input rpz coordinates in a 3-length array
xyz	output xyz coordinates in a 3-length array
r	r coordinate of the gradient
phi	phi coordinate of the gradient

Definition at line 44 of file math.c.

math_jac_xyz2rpz()

GPU_DECLARE_TARGET_SIMD void math_jac_xyz2rpz	(	real *	xyz,
		real *	rpz,
		real	r,
		real	phi )

Convert a Jacobian from cartesian to cylindrical coordinates.

This function converts a Jacobian located at angle phi and radius r from cartesian to cylindrical coordinates. The input Jacobian is an array with

[Bx dBx/dx dBx/dy dBx/dz By dBy/dx dBy/dy dBy/dz Bz dBz/dx dBz/dy dBz/dz]

an the output is

[Br dBr/dr dBr/dphi dBr/dz Bphi dBphi/dr dBphi/dphi dBphi/dz Bz dBz/dr dBz/dphi dBz/dz].

Parameters

xyz	output xyz coordinates in a 3-length array
rpz	input rpz coordinates in a 3-length array
r	r coordinate of the gradient
phi	phi coordinate of the gradient

Definition at line 103 of file math.c.

math_matmul()

DECLARE_TARGET_END GPU_DECLARE_TARGET_SIMD void math_matmul	(	real *	matA,
		real *	matB,
		int	d1,
		int	d2,
		int	d3,
		real *	matC )

Matrix multiplication.

This function performs matrix multiplication. The matrices are given as arrays, the indexing being row-wise.

Parameters

matA	input array representing a d1 x d2 matrix
matB	input array representing a d2 x d3 matrix
d1	input scalar dimension (columns in matA and matC)
d2	input scalar dimension (rows in matA and columns in matB)
d3	input scalar dimension (rows in matB and matC)
matC	output array representing a d1 x d3 matrix

Definition at line 159 of file math.c.

math_normal_rand()

DECLARE_TARGET_END DECLARE_TARGET_SIMD real math_normal_rand

(

void

)

Generate normally distributed random numbers.

This function generates normally distributed random numbers with expected value of 0 and variance of 1 using the polar method [1]. The rand48 random number generator should be initialized with srand48() before calling this function.

[1] G. Marsaglia, T. A. Bray. A convenient method for generating normal variables. Siam Review 6.3 (1964): 260-264. http://epubs.siam.org/doi/abs/10.1137/1006063:w

Todo

Currently only one of the generated numbers is returned; a second independent variate is Y = v2 * sqrt(-2*log(s) / s)

Try other random number generators such as those in Intel MKL

Definition at line 188 of file math.c.

math_ipow()

DECLARE_TARGET_SIMD int math_ipow	(	int	a,
		int	p )

Calculate a^p where both a and p are integers (p >= 0)

Parameters

a	argument
p	power

Returns

a^b

Definition at line 210 of file math.c.

math_simpson()

double math_simpson	(	double(*)(double)	f,
		double	a,
		double	b,
		double	eps )

Adaptive Simpsons rule for integral.

This function uses recursive splitting of the interval until either maximum number of intervals is reached or given relative error tolerance is obtained

Parameters

f	pointer to the integrand function (of one variable)
a	lower limit for the interval
b	upper limit for the interval
eps	absolute tolerance

Bug

Probably does not work with SIMD

Definition at line 232 of file math.c.

math_linspace()

void math_linspace	(	real *	vec,
		real	a,
		real	b,
		int	n )

Generate linearly spaced vector.

Generates linearly space vector with n elements whose end points are a and b. If n = 1, return b.

Parameters

vec	n length array where result is stored
a	start point
b	end point
n	number of elements

Definition at line 251 of file math.c.

math_point_on_plane()

DECLARE_TARGET_SIMD int math_point_on_plane	(	real	q[3],
		real	t1[3],
		real	t2[3],
		real	t3[3] )

Find if a point is on a given plane.

Let t1=(x1, y1, z1), t2=(x2, y2, z2), and t3=(x3, y3, z3) be non-collinear points on the plane. Let q =(x, y, z) be an arbitrary point. It lies in the plane if:

det( |x -x1 y -y1 z -z1| |x2-x1 y2-y1 z2-z1| |x3-x1 y3-y1 z3-z1| ) == 0

det( |x -x1 y -y1 z -z1| |x -x2 y -y2 z -z2| |x -x3 y -y3 z -z3| ) == 0

Parameters

q	xyz coordinates of a query point
t1	xyz coordinates of the first point defining the plane
t2	xyz coordinates of the second point defining the plane
t3	xyz coordinates of the third point defining the plane

Definition at line 324 of file math.c.

math_barycentric_coords_triangle()

DECLARE_TARGET_SIMD void math_barycentric_coords_triangle	(	real	AP[3],
		real	AB[3],
		real	AC[3],
		real	n[3],
		real *	s,
		real *	t )

Find barycentric coordinates for a given point.

Let a,b,c be the vertices of a triangle and p an arbitrary point. The s,t are the barycentric coordinates of p if: p = (1-s-t)*a + s*b + t*c

Parameters

AP	point position vector [x,y,z]
AB	triangle edge AB vector
AC	triangle edge AC vector
n	normal vector [nx, ny, nz] of the triangle
s	evaluated s coordinate
t	evaluated t coordinate

Definition at line 360 of file math.c.

math_uniquecount()

void math_uniquecount	(	int *	in,
		int *	unique,
		int *	count,
		int	n )

Find unique numbers and their frequency in given array.

All argument arrays are assumed to have the same length as the array holding the input numbers. If not all numbers are unique, the remaining entries in unique and count arguments are zero.

Parameters

in	input array of length n
unique	array where unique numbers are stored, must be of length n
count	count[i] is how many times unique[i] appears in array in, must be of length n
n	length of the argument arrays

Definition at line 276 of file math.c.

math_rsearch()

DECLARE_TARGET_SIMD real * math_rsearch	(	const real	key,
		const real *	base,
		int	num )

Search for array element preceding a key value.

This function takes an array of real values and a key value, and finds the array element which precedes the key value. The array to be searched has to be in ascending sorted order (according to comparison function rcomp).

Parameters

key	real value that serves as key for the search
base	pointer to the first object of the array to be searched
num	number of elements in array

Returns

pointer to array element preceding the key, or NULL if search failed

Definition at line 405 of file math.c.

math_point_in_polygon()

int math_point_in_polygon	(	real	r,
		real	z,
		real *	rv,
		real *	zv,
		int	n )

Check if coordinates are within polygon.

This function checks if the given coordinates are within a 2D polygon using a modified axis crossing method [1]. Origin is moved to the coordinates and the number of wall segments crossing the positive r-axis are calculated. If this is odd, the point is inside the polygon.

[1] D.G. Alciatore, R. Miranda. A Winding Number and Point-in-Polygon Algorithm. Technical report, Colorado State University, 1995. http://www.engr.colostate.edu/~dga/dga/papers/point_in_polygon.pdf

Parameters

r	r coordinate
z	z coordinate
rv	array of r points for the polygon
zv	array of z points for the polygon
n	number of points in the polygon

Definition at line 427 of file math.c.