diff --git a/Makefile b/Makefile index 3feabbd..ade3a5f 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,8 @@ MODULE_big = pg_sphere OBJS = sscan.o sparse.o sbuffer.o vector3d.o point.o \ euler.o circle.o line.o ellipse.o polygon.o \ path.o box.o output.o gq_cache.o gist.o key.o \ - gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o + gnomo.o healpix.o moc.o process_moc.o healpix_bare/healpix_bare.o \ + epochprop.o EXTENSION = pg_sphere RELEASE_SQL = $(EXTENSION)--$(PGSPHERE_VERSION).sql @@ -23,13 +24,13 @@ DATA_built = $(RELEASE_SQL) \ DOCS = README.pg_sphere COPYRIGHT.pg_sphere REGRESS = init tables points euler circle line ellipse poly path box index \ contains_ops contains_ops_compat bounding_box_gist gnomo healpix \ - moc mocautocast + moc mocautocast epochprop REGRESS_9_5 = index_9.5 # experimental for spoint3 TESTS = init_test tables points euler circle line ellipse poly path box index \ contains_ops contains_ops_compat bounding_box_gist gnomo healpix \ - moc mocautocast + moc mocautocast epochprop ifndef CXXFLAGS # no support for CXXFLAGS in PGXS before v11 @@ -39,7 +40,7 @@ endif EXTRA_CLEAN = $(PGS_SQL) pg_sphere.test.sql -CRUSH_TESTS = init_extended circle_extended +CRUSH_TESTS = init_extended circle_extended # order of sql files is important PGS_SQL = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \ @@ -47,7 +48,7 @@ PGS_SQL = pgs_types.sql pgs_point.sql pgs_euler.sql pgs_circle.sql \ pgs_box.sql pgs_contains_ops.sql pgs_contains_ops_compat.sql \ pgs_gist.sql gnomo.sql \ healpix.sql pgs_gist_spoint3.sql pgs_moc_type.sql pgs_moc_compat.sql pgs_moc_ops.sql \ - pgs_moc_geo_casts.sql + pgs_moc_geo_casts.sql pgs_epochprop.sql PGS_SQL_9_5 = pgs_9.5.sql # experimental for spoint3 USE_PGXS = 1 @@ -199,7 +200,7 @@ ifeq ($(has_parallel), n) sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL endif -pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in +pg_sphere--1.2.0--1.2.1.sql: pgs_moc_geo_casts.sql.in pgs_epochprop.sql.in cat $^ > $@ ifeq ($(has_parallel), n) sed -i -e '/PARALLEL/d' $@ # version $(pg_version) does not have support for PARALLEL diff --git a/doc/functions.sgm b/doc/functions.sgm index 1318141..30118bf 100644 --- a/doc/functions.sgm +++ b/doc/functions.sgm @@ -667,5 +667,133 @@ - + + + + Epoch propagation + + + + 6-Parameter Epoch Propagation + + + double precision[6] + epoch_prop + spoint pos + double precision parallax + double precision pm_long + double precision pm_lat + double precision radial_velocity + double precision delta_t + + + + Propagates a spherical phase vector in time (in particular, + applies proper motion to positions) + + + Following both pg_sphere and, where missing, astronomical + conventions makes units somewhat eclectic here; pm_long and pm_lat + need to be in rad/yr, whereas parallax is in mas, and + radial_velocity in km/s. The time difference must be in + (Julian) years. + + + + This function returns a 6-array of [long, lat, parallax, + pm_long, pm_lat, radial_velocity] of the corresponding values + delta_t years after the reference epoch for the original position. + As in the function arguments, long and lat are in rad, pm_lon and + pm_lat are in rad/yr, parallax is in mas, and radial_velocity is + in km/s. If you are only interested in the position, consider + the epoch_prop_pos functions below that have a somewhat less + contorted signature. + + + + It is an error to have either pos or delta_t NULL. For all + other arguments, NULLs are turned into 0s, except for parallax, + where some very small default is put in. In that case, + both parallax and radial_velocity will be NULL in the output + array. + + + + This uses the rigorous method derived in "The Hipparcos and Tycho + Catalogues", ESA Special Publication 1200 (1997), p 94f. It does + not take into account relativistic effects, and it also does not + account for secular aberration. + + + Propagating Barnard's star into the past + + + + + + Epoch Propagation of Positions Only + + + spoint + epoch_prop_pos + spoint pos + double precision parallax + double precision pm_long + double precision pm_lat + double precision radial_velocity + double precision delta_t + + + + + spoint + epoch_prop_pos + spoint pos + double precision pm_long + double precision pm_lat + double precision delta_t + + + + These are simplified versions of epoch_prop returning only spoints; + the propagated values for the other coordinates are discarded + (but still internallay computed; these functions do not + run any faster than epoch_prop itself). + + + + As with epoch_prop itself, missing values (except for pos and + delta_t) are substituted by 0 (or a very small value in the + case of parallax). + + + Barnard's star, position and proper motion + + + + diff --git a/epochprop.c b/epochprop.c new file mode 100644 index 0000000..4c7af21 --- /dev/null +++ b/epochprop.c @@ -0,0 +1,203 @@ +/* Handling of (conventional) proper motions. + +This code is largely based on a FORTRAN function written by Lennart Lindegren +(Lund Obs) in 1995 that implements the procedure described in The Hipparcos +and Tycho Catalogues (ESA SP-1200), Volume 1, Section 1.5.5, 'Epoch +Transformation: Rigorous Treatment'; cf. +. +*/ + +#include +#include + +#include "point.h" +#include "epochprop.h" +#include "vector3d.h" + +PG_FUNCTION_INFO_V1(epoch_prop); + +/* Astronomical unit in kilometers */ +#define AU 1.495978701e8 + +/* Julian year in seconds */ +#define J_YEAR (365.25*86400) + +/* A_nu as per ESA/SP-1200 */ +#define A_NU (AU/(J_YEAR)) + +/* Following SOFA, we use 1e-7 arcsec as minimal parallax + ("celestial sphere"); parallax=0 exactly means "infinite distance", which + leads to all sorts for problems; our parallaxes come in in mas, so: */ +#define PX_MIN 1e-7*1000 + +/* propagate an object at a phase vector over a time difference of delta_t, +stuffing an updated phase vector in result. + +This does not propagate errors. +*/ +static void propagate_phasevec( + const phasevec *pv, + const double delta_t, + phasevec *result) { + + double distance_factor, mu0abs, zeta0, parallax; + + Vector3D p0, r0, q0; + Vector3D mu0, pprime, qprime, mu, muprime, u, uprime; + + /* for very small or null parallaxes, our algorithm breaks; avoid that + and, if we did emergency measures, do not talk about parallax and + radial velocity in the output */ + if (pv->parallax_valid) { + parallax = pv->parallax; + } else { + parallax = PX_MIN; + } + result->parallax_valid = pv->parallax_valid; + + /* compute the normal triad as Vector3D-s, eq. (1.2.15)*/ + spoint_vector3d(&r0, &(pv->pos)); + + p0.x = -sin(pv->pos.lng); + p0.y = cos(pv->pos.lng); + p0.z = 0; + + q0.x = -sin(pv->pos.lat) * cos(pv->pos.lng); + q0.y = -sin(pv->pos.lat) * sin(pv->pos.lng); + q0.z = cos(pv->pos.lat); + + /* the original proper motion vector */ + mu0.x = mu0.y = mu0.z = 0; + vector3d_addwithscalar(&mu0, pv->pm[0], &p0); + vector3d_addwithscalar(&mu0, pv->pm[1], &q0); + mu0abs = vector3d_length(&mu0); + + /* radial velocity in mas/yr ("change of parallax per year"). eq. (1.5.24) + We're transforming this to rad/yr so the units work out below */ + zeta0 = (pv->rv * parallax / A_NU) / 3.6e6 / RADIANS; + /* distance factor eq. (1.5.25) */ + distance_factor = 1/sqrt(1 + + 2 * zeta0 * delta_t + + (mu0abs * mu0abs + zeta0 * zeta0) * delta_t * delta_t); + + /* the propagated proper motion vector, eq. (1.5.28) */ + muprime.x = muprime.y = muprime.z = 0; + vector3d_addwithscalar(&muprime, (1 + zeta0 * delta_t), &mu0); + vector3d_addwithscalar(&muprime, -mu0abs * mu0abs * delta_t, &r0); + mu.x = mu.y = mu.z = 0; + vector3d_addwithscalar(&mu, pow(distance_factor, 3), &muprime); + + /* parallax, eq. (1.5.27) */ + result->parallax = distance_factor*parallax; + /* zeta/rv, eq. (1.5.29); go back from rad to mas, too */ + result->rv = (zeta0 + (mu0abs * mu0abs + zeta0 * zeta0) * delta_t) + * distance_factor * distance_factor + * 3.6e6 * RADIANS + * A_NU / result->parallax; + + /* propagated position, eq. (1.5.26) */ + uprime.x = uprime.y = uprime.z = 0; + vector3d_addwithscalar(&uprime, (1 + zeta0 * delta_t), &r0); + vector3d_addwithscalar(&uprime, delta_t, &mu0); + u.x = u.y = u.z = 0; + vector3d_addwithscalar(&u, distance_factor, &uprime); + vector3d_spoint(&(result->pos), &u); + + /* compute a new triad for the propagated position, eq (1.5.15) */ + pprime.x = -sin(result->pos.lng); + pprime.y = cos(result->pos.lng); + pprime.z = 0; + qprime.x = -sin(result->pos.lat) * cos(result->pos.lng); + qprime.y = -sin(result->pos.lat) * sin(result->pos.lng); + qprime.z = cos(result->pos.lat); + + /* use it to compute the proper motions, eq. (1.5.32) */ + result->pm[0] = vector3d_scalar(&pprime, &mu); + result->pm[1] = vector3d_scalar(&qprime, &mu); +} + + +/* + Propagate a position with proper motions and optionally parallax + and radial velocity. + + Arguments: pos0 (spoint), pm_long, pm_lat (in rad/yr) + par (parallax, mas), rv (in km/s), delta_t (in years) + + This returns a 6-array of lat, long (in rad), parallax (in mas) + pmlat, pmlong (in rad/yr), rv (in km/s). +*/ +Datum +epoch_prop(PG_FUNCTION_ARGS) { + double delta_t; + phasevec input, output; + ArrayType *result; + Datum retvals[6]; + + if (PG_ARGISNULL(0)) { + ereport(ERROR, + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), + errmsg("NULL position not supported in epoch propagation"))); } + memcpy(&(input.pos), (void*)PG_GETARG_POINTER(0), sizeof(SPoint)); + if (PG_ARGISNULL(1)) { + input.parallax = 0; + } else { + input.parallax = PG_GETARG_FLOAT8(1); + } + input.parallax_valid = fabs(input.parallax) > PX_MIN; + + if (PG_ARGISNULL(2)) { + input.pm[0] = 0; + } else { + input.pm[0] = PG_GETARG_FLOAT8(2); + } + + if (PG_ARGISNULL(3)) { + input.pm[1] = 0; + } else { + input.pm[1] = PG_GETARG_FLOAT8(3); + } + + if (PG_ARGISNULL(4)) { + input.rv = 0; + } else { + input.rv = PG_GETARG_FLOAT8(4); + } + + if (PG_ARGISNULL(5)) { + ereport(ERROR, + (errcode(ERRCODE_NULL_VALUE_NOT_ALLOWED), + errmsg("NULL delta t not supported in epoch propagation"))); } + delta_t = PG_GETARG_FLOAT8(5); + + propagate_phasevec(&input, delta_t, &output); + + /* change to internal units: rad, rad/yr, mas, and km/s */ + retvals[0] = Float8GetDatum(output.pos.lng); + retvals[1] = Float8GetDatum(output.pos.lat); + retvals[2] = Float8GetDatum(output.parallax); + retvals[3] = Float8GetDatum(output.pm[0]); + retvals[4] = Float8GetDatum(output.pm[1]); + retvals[5] = Float8GetDatum(output.rv); + + { + bool isnull[6] = {0, 0, 0, 0, 0, 0}; + int lower_bounds[1] = {1}; + int dims[1] = {6}; +#ifdef USE_FLOAT8_BYVAL + bool embyval = true; +#else + bool embyval = false; +#endif + + if (! output.parallax_valid) { + /* invalidate parallax and rv */ + isnull[2] = 1; + isnull[5] = 1; + } + + result = construct_md_array(retvals, isnull, 1, dims, lower_bounds, + FLOAT8OID, sizeof(float8), embyval, 'd'); + } + PG_RETURN_ARRAYTYPE_P(result); +} diff --git a/epochprop.h b/epochprop.h new file mode 100644 index 0000000..9e1b87e --- /dev/null +++ b/epochprop.h @@ -0,0 +1,26 @@ +/* Definitions for handling proper motions */ + +#include + + +Datum epoch_prop(PG_FUNCTION_ARGS); + + +/* a cartesian point; this is like geo_decl's point, but you can't +have both geo_decls and pg_sphere right now (both define a type Point, +not to mention they have different ideas on EPSILON */ +typedef struct s_cpoint { + double x, y; +} CPoint; + +typedef struct s_phasevec { + SPoint pos; /* Position as an SPoint */ + double pm[2]; /* Proper motion long/lat in rad/year, + PM in longitude has cos(lat) applied */ + double parallax; /* in rad */ + double rv; /* radial velocity in km/s */ + int parallax_valid; /* 1 if the parallax really is a NULL */ +} phasevec; + + + diff --git a/expected/epochprop.out b/expected/epochprop.out new file mode 100644 index 0000000..8541601 --- /dev/null +++ b/expected/epochprop.out @@ -0,0 +1,107 @@ +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 546.9759, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+----------+------------+---------- + 269.4742714391 | 4.4072939987 | 543.624 | -791.442 | 10235.412 | -110.450 +(1 row) + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 0, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+---------+----------+------------+--------- + 269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | +(1 row) + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + NULL, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+---------+----------+------------+--------- + 269.4744079540 | 4.4055337210 | | -801.210 | 10361.762 | +(1 row) + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 23, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL, + 20) AS tp) AS q; + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+----------+------------+---------- + 269.4476085384 | 4.7509315989 | 23.000 | -801.617 | 10361.984 | 2.159 +(1 row) + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 23, + NULL, RADIANS(10362/3.6e6), -110, + 120) AS tp) AS q; + to_char | to_char | to_char | to_char | to_char | to_char +-----------------+-----------------+----------+----------+------------+---------- + 269.4520769500 | 5.0388680565 | 23.007 | -.000 | 10368.061 | -97.120 +(1 row) + +SELECT epoch_prop(NULL, + 23, + 0.01 , RADIANS(10362/3.6e6), -110, + 120); +ERROR: NULL position not supported in epoch propagation +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), + 23, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + 20) AS tp; + tp +----------------------------------------- + (4.70274792658313 , 0.0829194509345993) +(1 row) + +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), + 20) AS tp; + tp +----------------------------------------- + (4.70274793061952 , 0.0829193989380876) +(1 row) + diff --git a/pgs_epochprop.sql.in b/pgs_epochprop.sql.in new file mode 100644 index 0000000..754c51a --- /dev/null +++ b/pgs_epochprop.sql.in @@ -0,0 +1,41 @@ +-- *********************************** +-- +-- functions for propagating positions +-- +-- *********************************** + +CREATE FUNCTION epoch_prop( + pos spoint, + parallax DOUBLE PRECISION, + pmlng DOUBLE PRECISION, + pmlat DOUBLE PRECISION, + rv DOUBLE PRECISION, + delta_t DOUBLE PRECISION) + RETURNS double precision[6] + AS 'MODULE_PATHNAME', 'epoch_prop' + LANGUAGE 'c' + IMMUTABLE PARALLEL SAFE; + +CREATE FUNCTION epoch_prop_pos( + pos spoint, + parallax DOUBLE PRECISION, + pmlng DOUBLE PRECISION, + pmlat DOUBLE PRECISION, + rv DOUBLE PRECISION, + delta_t DOUBLE PRECISION) + RETURNS spoint + AS $body$ + SELECT spoint(pv[1], pv[2]) + FROM epoch_prop(pos, parallax, pmlng, pmlat, rv, delta_t) as pv + $body$ LANGUAGE SQL STABLE PARALLEL SAFE; + +CREATE FUNCTION epoch_prop_pos( + pos spoint, + pmlng DOUBLE PRECISION, + pmlat DOUBLE PRECISION, + delta_t DOUBLE PRECISION) + RETURNS spoint + AS $body$ + SELECT spoint(pv[1], pv[2]) + FROM epoch_prop(pos, 0, pmlng, pmlat, 0, delta_t) as pv + $body$ LANGUAGE SQL STABLE PARALLEL SAFE; diff --git a/sql/epochprop.sql b/sql/epochprop.sql new file mode 100644 index 0000000..fea9737 --- /dev/null +++ b/sql/epochprop.sql @@ -0,0 +1,78 @@ +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 546.9759, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 0, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + NULL, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + -100) AS tp) AS q; + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 23, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), NULL, + 20) AS tp) AS q; + +SELECT + to_char(DEGREES(tp[1]), '999D9999999999'), + to_char(DEGREES(tp[2]), '999D9999999999'), + to_char(tp[3], '999D999'), + to_char(DEGREES(tp[4])*3.6e6, '999D999'), + to_char(DEGREES(tp[5])*3.6e6, '99999D999'), + to_char(tp[6], '999D999') +FROM ( + SELECT epoch_prop(spoint(radians(269.45207695), radians(4.693364966)), + 23, + NULL, RADIANS(10362/3.6e6), -110, + 120) AS tp) AS q; + +SELECT epoch_prop(NULL, + 23, + 0.01 , RADIANS(10362/3.6e6), -110, + 120); + +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), + 23, + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), -110, + 20) AS tp; + +SELECT epoch_prop_pos(spoint(radians(269.45207695), radians(4.693364966)), + RADIANS(-801.551/3.6e6), RADIANS(10362/3.6e6), + 20) AS tp; diff --git a/vector3d.c b/vector3d.c index adaf5f7..27b7d8f 100644 --- a/vector3d.c +++ b/vector3d.c @@ -18,6 +18,14 @@ vector3d_cross(Vector3D *out, const Vector3D *v1, const Vector3D *v2) out->z = v1->x * v2->y - v1->y * v2->x; } +void +vector3d_addwithscalar(Vector3D *result, double scalar, const Vector3D *delta) +{ + result->x = result->x + delta->x * scalar; + result->y = result->y + delta->y * scalar; + result->z = result->z + delta->z * scalar; +} + float8 vector3d_scalar(Vector3D *v1, Vector3D *v2) { diff --git a/vector3d.h b/vector3d.h index 55d0f85..23ba9b9 100644 --- a/vector3d.h +++ b/vector3d.h @@ -37,4 +37,10 @@ float8 vector3d_scalar(Vector3D *v1, Vector3D *v2); */ float8 vector3d_length(const Vector3D *v); +/* + * Calculate result + scalar*delta + */ +void vector3d_addwithscalar( + Vector3D *result, double scalar, const Vector3D *delta); + #endif