src/svector.c

   1 #include "svecmat.h"
   2
   3 void hkl_svector_fprintf(FILE * file, struct hkl_svector const * v)
   4 {
   5   fprintf(file, "|%f, %f, %f|", v->data[0], v->data[1], v->data[2]);
   6 }
   7
   8 int hkl_svector_cmp(struct hkl_svector const * v, struct hkl_svector const * v1)
   9 {
  10   unsigned int i;
  11
  12   for (i=0; i<3; i++)
  13     if ( fabs(v->data[i] - v1->data[i]) > HKL_EPSILON )
  14       return HKL_FALSE;
  15   return HKL_TRUE;
  16 }
  17
  18 void hkl_svector_minus_svector(struct hkl_svector * v, struct hkl_svector const * v1)
  19 {
  20   unsigned int i;
  21   for (i=0;i<3;i++)
  22     v->data[i] -= v1->data[i];
  23 }
  24
  25 void hkl_svector_div_double(struct hkl_svector * v, double const d)
  26 {
  27   unsigned int i;
  28   for (i=0;i<3;i++)
  29     v->data[i] /= d;
  30 }
  31
  32 void hkl_svector_times_double(struct hkl_svector * v, double const d)
  33 {
  34   unsigned int i;
  35   for (i=0;i<3;i++)
  36     v->data[i] *= d;
  37 }
  38
  39 void hkl_svector_times_svector(struct hkl_svector * v, struct hkl_svector const * v1)
  40 {
  41   unsigned int i;
  42   for (i=0;i<3;i++)
  43     v->data[i] *= v1->data[i];
  44 }
  45
  46 void hkl_svector_times_smatrix(struct hkl_svector * v, struct hkl_smatrix const * m)
  47 {
  48   struct hkl_svector tmp;
  49   tmp = *v;
  50
  51   v->data[0] = tmp.data[0] * m->data[0][0] + tmp.data[1] * m->data[1][0] + tmp.data[2] * m->data[2][0];
  52   v->data[1] = tmp.data[0] * m->data[0][1] + tmp.data[1] * m->data[1][1] + tmp.data[2] * m->data[2][1];
  53   v->data[2] = tmp.data[0] * m->data[0][2] + tmp.data[1] * m->data[1][2] + tmp.data[2] * m->data[2][2];
  54 }
  55
  56 double hkl_svector_sum(struct hkl_svector const * v)
  57 {
  58   return v->data[0] + v->data[1] + v->data[2];
  59 }
  60
  61 double hkl_svector_scalar_product(struct hkl_svector const * v, struct hkl_svector const * v1)
  62 {
  63   unsigned int i;
  64   double scalar = 0;
  65
  66   for (i=0;i<3;i++)
  67     scalar += v->data[i] * v1->data[i];
  68   return scalar;
  69 }
  70
  71 void hkl_svector_vectorial_product(struct hkl_svector * v, struct hkl_svector const * v1)
  72 {
  73   struct hkl_svector tmp;
  74
  75   tmp = *v;
  76   v->data[0] = tmp.data[1] * v1->data[2] - tmp.data[2] * v1->data[1];
  77   v->data[1] = tmp.data[2] * v1->data[0] - tmp.data[0] * v1->data[2];
  78   v->data[2] = tmp.data[0] * v1->data[1] - tmp.data[1] * v1->data[0];
  79 }
  80
  81
  82 double hkl_svector_angle(struct hkl_svector const * v, struct hkl_svector const * v1)
  83 {
  84   double angle;
  85   double norm = hkl_svector_norm2(v) * hkl_svector_norm2(v1);
  86
  87   double cosinus_angle = hkl_svector_scalar_product(v, v1) / norm;
  88
  89   // problem with round
  90   if (cosinus_angle >= 1 )
  91     angle = 0;
  92   else
  93     if (cosinus_angle <= -1 )
  94       angle = M_PI;
  95     else
  96       angle = acos(cosinus_angle);
  97
  98   return angle;
  99 }
 100
 101 double hkl_svector_norm2(struct hkl_svector const * v)
 102 {
 103   return sqrt(v->data[0] * v->data[0] + v->data[1] * v->data[1] + v->data[2] * v->data[2]);
 104 }
 105
 106 /**
 107  * @brief normalize a hkl_svector
 108  * @return true if the hkl_svector can be normalized, false otherwise
 109  * @todo check the status
 110  */
 111 int hkl_svector_normalize(struct hkl_svector * v)
 112 {
 113   int status = HKL_FAIL;
 114
 115   double norm = hkl_svector_norm2(v);
 116   if ( norm > HKL_EPSILON )
 117     {
 118       hkl_svector_div_double(v, norm);
 119       status = HKL_SUCCESS;
 120     }
 121
 122   return status;
 123 }
 124
 125 int hkl_svector_is_colinear(struct hkl_svector const * v, struct hkl_svector const * v1)
 126 {
 127   int is_colinear = 0;
 128   struct hkl_svector tmp = *v;
 129
 130   hkl_svector_vectorial_product(&tmp, v1);
 131   if (hkl_svector_norm2(&tmp) < HKL_EPSILON)
 132     is_colinear = 1;
 133
 134   return is_colinear;
 135 }
 136
 137
 138 void hkl_svector_randomize(struct hkl_svector * v)
 139 {
 140   v->data[0] = -1 + 2 * rand()/(RAND_MAX+1.0);
 141   v->data[1] = -1 + 2 * rand()/(RAND_MAX+1.0);
 142   v->data[2] = -1 + 2 * rand()/(RAND_MAX+1.0);
 143 }
 144
 145 void hkl_svector_randomize_svector(struct hkl_svector * v, struct hkl_svector const * v1)
 146 {
 147   do
 148     hkl_svector_randomize(v);
 149   while (hkl_svector_cmp(v, v1) == HKL_TRUE);
 150 }
 151
 152 void hkl_svector_randomize_svector_svector(struct hkl_svector * v, struct hkl_svector const * v1, struct hkl_svector const * v2)
 153 {
 154   do
 155     hkl_svector_randomize(v);
 156   while (hkl_svector_cmp(v, v1) == HKL_TRUE || hkl_svector_cmp(v, v2) == HKL_TRUE);
 157 }
 158
 159
 160 void hkl_svector_rotated_around_vector(struct hkl_svector * v, struct hkl_svector const * axe, double angle)
 161 {
 162   double c = cos(angle);
 163   double s = sin(angle);
 164   struct hkl_svector axe_n;
 165   struct hkl_svector tmp;
 166
 167   axe_n = *axe;
 168   hkl_svector_normalize(&axe_n);
 169
 170   tmp = *v;
 171
 172   v->data[0] = (c + (1 - c) * axe_n.data[0] * axe_n.data[0]) * tmp.data[0];
 173   v->data[0] += ((1 - c) * axe_n.data[0] * axe_n.data[1] - axe_n.data[2] * s) * tmp.data[1];
 174   v->data[0] += ((1 - c) * axe_n.data[0] * axe_n.data[2] + axe_n.data[1] * s) * tmp.data[2];
 175
 176   v->data[1] = ((1 - c) * axe_n.data[0] * axe_n.data[1] + axe_n.data[2] * s) * tmp.data[0];
 177   v->data[1] += (c + (1 - c) * axe_n.data[1] * axe_n.data[1]) * tmp.data[1];
 178   v->data[1] += ((1 - c) * axe_n.data[1] * axe_n.data[2] - axe_n.data[0] * s) * tmp.data[2];
 179
 180   v->data[2] = ((1 - c) * axe_n.data[0] * axe_n.data[2] - axe_n.data[1] * s) * tmp.data[0];
 181   v->data[2] += ((1 - c) * axe_n.data[1] * axe_n.data[2] + axe_n.data[0] * s) * tmp.data[1];
 182   v->data[2] += (c + (1 - c) * axe_n.data[2] * axe_n.data[2]) * tmp.data[2];
 183 }
 184
 185 /**
 186  * @brief check if the hkl_svector is null
 187  * @return true if all |elements| are below HKL_EPSILON, false otherwise
 188  * @todo test
 189  */
 190 int hkl_svector_is_null(struct hkl_svector const * v)
 191 {
 192   unsigned int i;
 193   for (i=0; i<3; i++)
 194     if ( fabs(v->data[i]) > HKL_EPSILON )
 195       return HKL_FALSE;
 196   return HKL_TRUE;
 197 }