hkl/hkl-interval.c

   1 /* This file is part of the hkl library.
   2  *
   3  * The hkl library is free software: you can redistribute it and/or modify
   4  * it under the terms of the GNU General Public License as published by
   5  * the Free Software Foundation, either version 3 of the License, or
   6  * (at your option) any later version.
   7  *
   8  * The hkl library is distributed in the hope that it will be useful,
   9  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  10  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11  * GNU General Public License for more details.
  12  *
  13  * You should have received a copy of the GNU General Public License
  14  * along with the hkl library.  If not, see <http://www.gnu.org/licenses/>.
  15  *
  16  * Copyright (C) 2003-2010 Synchrotron SOLEIL
  17  *                         L'Orme des Merisiers Saint-Aubin
  18  *                         BP 48 91192 GIF-sur-YVETTE CEDEX
  19  *
  20  * Authors: Picca Frédéric-Emmanuel <picca@synchrotron-soleil.fr>
  21  */
  22 #define _GNU_SOURCE
  23 #include <math.h>
  24 #include <gsl/gsl_sf_trig.h>
  25
  26 #include "hkl/hkl-interval.h"
  27
  28 /** compare two intervals */
  29 int hkl_interval_cmp(HklInterval const *self, HklInterval const *interval)
  30 {
  31         return self->min == interval->min && self->max == interval->max;
  32 }
  33
  34 /** add two intervals */
  35 void hkl_interval_plus_interval(HklInterval *self, HklInterval const *interval)
  36 {
  37         self->min += interval->min;
  38         self->max += interval->max;
  39 }
  40
  41 /** add to an interval a double */
  42 void hkl_interval_plus_double(HklInterval *self, double const d)
  43 {
  44         self->min += d;
  45         self->max += d;
  46 }
  47
  48
  49 void hkl_interval_minus_interval(HklInterval *self, HklInterval const *interval)
  50 {
  51         self->min -= interval->max;
  52         self->max -= interval->min;
  53 }
  54
  55
  56 void hkl_interval_minus_double(HklInterval *self, double const d)
  57 {
  58         self->min -= d;
  59         self->max -= d;
  60 }
  61
  62 void hkl_interval_times_interval(HklInterval *self, HklInterval const *interval)
  63 {
  64         double min;
  65         double max;
  66         double m1 = self->min * interval->min;
  67         double m2 = self->min * interval->max;
  68         double m3 = self->max * interval->min;
  69         double m4 = self->max * interval->max;
  70
  71         min = m1;
  72         if (m2 < min)
  73                 min = m2;
  74         if (m3 < min)
  75                 min = m3;
  76         if (m4 < min)
  77                 min = m4;
  78
  79         max = m1;
  80         if (m2 > max)
  81                 max = m2;
  82         if (m3 > max)
  83                 max = m3;
  84         if (m4 > max)
  85                 max = m4;
  86
  87         self->min = min;
  88         self->max = max;
  89 }
  90
  91 void hkl_interval_times_double(HklInterval *self, double const d)
  92 {
  93         double min;
  94         double max;
  95         if (d < 0) {
  96                 min = self->max * d;
  97                 max = self->min * d;
  98         } else {
  99                 min = self->min * d;
 100                 max = self->max * d;
 101         }
 102         self->min = min;
 103         self->max = max;
 104 }
 105
 106 void hkl_interval_divides_double(HklInterval *self, double const d)
 107 {
 108         double min = self->min / d;
 109         double max = self->max / d;
 110         if (min > max){
 111                 double tmp = min;
 112                 min = max;
 113                 max = tmp;
 114         }
 115         self->min = min;
 116         self->max = max;
 117 }
 118
 119 int hkl_interval_contain_zero(HklInterval const *self)
 120 {
 121         if (self->min <= 0 && self->max >= 0)
 122                 return HKL_TRUE;
 123         else
 124                 return HKL_FALSE;
 125 }
 126
 127 void hkl_interval_cos(HklInterval *self)
 128 {
 129         double min = 0;
 130         double max = 0;
 131         double cmin;
 132         double cmax;
 133
 134         cmin = cos(self->min);
 135         cmax = cos(self->max);
 136
 137         if (self->max - self->min >= 2 * M_PI) {
 138                 min = -1;
 139                 max = 1;
 140         } else {
 141                 int quad_min;
 142                 int quad_max;
 143
 144                 quad_min = (int)floor(self->min / M_PI_2) % 4;
 145                 if (quad_min < 0)
 146                         quad_min += 4;
 147
 148                 quad_max = (int)floor(self->max / M_PI_2) % 4;
 149                 if (quad_max < 0)
 150                         quad_max += 4;
 151
 152                 switch (quad_max) {
 153                 case 0:
 154                         switch (quad_min) {
 155                         case 0:
 156                                 min = cmax;
 157                                 max = cmin;
 158                                 break;
 159                         case 1:
 160                                 min = -1;
 161                                 max = 1;
 162                                 break;
 163                         case 2:
 164                                 min = cmin;
 165                                 max = 1;
 166                                 break;
 167                         case 3:
 168                                 if (cmin < cmax) {
 169                                         min = cmin;
 170                                         max = 1;
 171                                 } else {
 172                                         min = cmax;
 173                                         max = 1;
 174                                 }
 175                                 break;
 176                         }
 177                         break;
 178                 case 1:
 179                         switch (quad_min) {
 180                         case 0:
 181                                 min = cmax;
 182                                 max = cmin;
 183                                 break;
 184                         case 1:
 185                                 min = -1;
 186                                 max = 1;
 187                                 break;
 188                         case 2:
 189                                 if (cmin < cmax) {
 190                                         min = cmin;
 191                                         max = 1;
 192                                 } else {
 193                                         min = cmax;
 194                                         max = 1;
 195                                 }
 196                                 break;
 197                         case 3:
 198                                 min = cmax;
 199                                 max = 1;
 200                                 break;
 201                         }
 202                         break;
 203                 case 2:
 204                         switch (quad_min) {
 205                         case 0:
 206                                 min = -1;
 207                                 max = cmin;
 208                                 break;
 209                         case 1:
 210                                 if (cmin < cmax) {
 211                                         min = -1;
 212                                         max = cmax;
 213                                 } else {
 214                                         min = -1;
 215                                         max = cmin;
 216                                 }
 217                                 break;
 218                         case 2:
 219                                 if (cmin < cmax) {
 220                                         min = cmin;
 221                                         max = cmax;
 222                                 } else {
 223                                         min = -1;
 224                                         max = 1;
 225                                 }
 226                                 break;
 227                         case 3:
 228                                 min = -1;
 229                                 max = 1;
 230                                 break;
 231                         }
 232                         break;
 233                 case 3:
 234                         switch (quad_min) {
 235                         case 0:
 236                                 if (cmin < cmax) {
 237                                         min = -1;
 238                                         max = cmax;
 239                                 } else {
 240                                         min = -1;
 241                                         max = cmin;
 242                                 }
 243                                 break;
 244                         case 1:
 245                                 min = -1;
 246                                 max = cmax;
 247                                 break;
 248                         case 2:
 249                                 min = cmin;
 250                                 max = cmax;
 251                                 break;
 252                         case 3:
 253                                 if (cmin < cmax) {
 254                                         min = cmin;
 255                                         max = cmax;
 256                                 } else {
 257                                         min = -1;
 258                                         max = 1;
 259                                 }
 260                                 break;
 261                         }
 262                         break;
 263                 }
 264         }
 265         self->min = min;
 266         self->max = max;
 267 }
 268
 269 void hkl_interval_acos(HklInterval *self)
 270 {
 271         double tmp;
 272
 273         tmp = self->min;
 274         self->min = acos(self->max);
 275         self->max = acos(tmp);
 276 }
 277
 278
 279 void hkl_interval_sin(HklInterval *self)
 280 {
 281         double min = 0;
 282         double max = 0;
 283         double smin;
 284         double smax;
 285
 286         smin = sin(self->min);
 287         smax = sin(self->max);
 288
 289         /* if there is at least one period in b, then a = [-1, 1] */
 290         if ( self->max - self->min >= 2 * M_PI) {
 291                 min = -1;
 292                 max = 1;
 293         } else {
 294                 int quad_min;
 295                 int quad_max;
 296
 297                 quad_min = (int)floor(self->min / M_PI_2) % 4;
 298                 if (quad_min < 0)
 299                         quad_min += 4;
 300
 301                 quad_max = (int)floor(self->max / M_PI_2) % 4;
 302                 if (quad_max < 0)
 303                         quad_max += 4;
 304
 305                 switch (quad_max) {
 306                 case 0:
 307                         switch (quad_min) {
 308                         case 0:
 309                                 if (smin < smax) {
 310                                         min = smin;
 311                                         max = smax;
 312                                 } else {
 313                                         min = -1;
 314                                         max = 1;
 315                                 }
 316                                 break;
 317                         case 3:
 318                                 min = smin;
 319                                 max = smax;
 320                                 break;
 321                         case 1:
 322                                 if (smin > smax) {
 323                                         min = -1;
 324                                         max = smin;
 325                                 } else {
 326                                         min = -1;
 327                                         max = smax;
 328                                 }
 329                                 break;
 330                         case 2:
 331                                 min = -1;
 332                                 max = smax;
 333                                 break;
 334                         }
 335                         break;
 336                 case 1:
 337                         switch (quad_min) {
 338                         case 0:
 339                                 if (smin < smax) {
 340                                         min = smin;
 341                                         max = 1;
 342                                 } else {
 343                                         min = smax;
 344                                         max = 1;
 345                                 }
 346                                 break;
 347                         case 1:
 348                                 if (smin < smax) {
 349                                         min = -1;
 350                                         max = 1;
 351                                 } else {
 352                                         min = smax;
 353                                         max = smin;
 354                                 }
 355                                 break;
 356                         case 2:
 357                                 min = -1;
 358                                 max = 1;
 359                                 break;
 360                         case 3:
 361                                 min = smin;
 362                                 max = 1;
 363                                 break;
 364                         }
 365                         break;
 366                 case 2:
 367                         switch (quad_min) {
 368                         case 0:
 369                                 min = smax;
 370                                 max = 1;
 371                                 break;
 372                         case 1:
 373                         case 2:
 374                                 if (smin < smax) {
 375                                         min = -1;
 376                                         max = 1;
 377                                 } else {
 378                                         min = smax;
 379                                         max = smin;
 380                                 }
 381                                 break;
 382                         case 3:
 383                                 if (smin < smax) {
 384                                         min = smin;
 385                                         max = 1;
 386                                 } else {
 387                                         min = smax;
 388                                         max = 1;
 389                                 }
 390                                 break;
 391                         }
 392                         break;
 393                 case 3:
 394                         switch (quad_min) {
 395                         case 0:
 396                                 min = -1;
 397                                 max = 1;
 398                                 break;
 399                         case 1:
 400                                 min = -1;
 401                                 max = smin;
 402                                 break;
 403                         case 2:
 404                                 if (smin < smax) {
 405                                         min = -1;
 406                                         max = smax;
 407                                 } else {
 408                                         min = -1;
 409                                         max = smin;
 410                                 }
 411                                 break;
 412                         case 3:
 413                                 if (smin < smax) {
 414                                         min = smin;
 415                                         max = smax;
 416                                 } else {
 417                                         min = -1;
 418                                         max = 1;
 419                                 }
 420                                 break;
 421                         }
 422                         break;
 423                 }
 424         }
 425         self->min = min;
 426         self->max = max;
 427 }
 428
 429 void hkl_interval_asin(HklInterval *self)
 430 {
 431         self->min = asin(self->min);
 432         self->max = asin(self->max);
 433 }
 434
 435 void hkl_interval_tan(HklInterval *self)
 436 {
 437         int quadrant_down = (int)floor(self->min / M_PI_2);
 438         int quadrant_up = (int)floor(self->max / M_PI_2);
 439
 440         /* if there is at least one period in b or if b contains a Pi/2 + k*Pi, */
 441         /* then a = ]-oo, +oo[ */
 442         if ( ((quadrant_up - quadrant_down) >= 2)
 443              || (!(quadrant_down % 2) && (quadrant_up % 2)) ) {
 444                 self->min = -INFINITY;
 445                 self->max = INFINITY;
 446         } else {
 447                 self->min = tan(self->min);
 448                 self->max = tan(self->max);
 449         }
 450 }
 451
 452 void hkl_interval_atan(HklInterval *self)
 453 {
 454         self->min = atan(self->min);
 455         self->max = atan(self->max);
 456 }
 457
 458 double hkl_interval_length(HklInterval const *self)
 459 {
 460         return self->max - self->min;
 461 }
 462
 463 void hkl_interval_angle_restrict_symm(HklInterval *self)
 464 {
 465         gsl_sf_angle_restrict_symm_e(&self->min);
 466         gsl_sf_angle_restrict_symm_e(&self->max);
 467 }