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-2017 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 <gsl/gsl_sf_trig.h>
  24 #include <math.h>                       // for floor, M_PI_2, acos, asin, etc
  25 #include <stdlib.h>                     // for free, NULL
  26 #include "hkl-interval-private.h"       // for HklInterval
  27 #include "hkl-macros-private.h"         // for HKL_MALLOC
  28 #include "hkl.h"                        // for FALSE, TRUE
  29
  30 /**
  31  * hkl_interval_dup: (skip)
  32  * @self:
  33  *
  34  * copy an #HklInterval
  35  *
  36  * Returns:
  37  **/
  38 HklInterval *hkl_interval_dup(const HklInterval *self)
  39 {
  40         if(!self)
  41                 return NULL;
  42
  43         HklInterval *dup = HKL_MALLOC(HklInterval);
  44
  45         *dup = *self;
  46
  47         return dup;
  48 }
  49
  50 /**
  51  * hkl_interval_free: (skip)
  52  * @self:
  53  *
  54  * delete an #HklInterval
  55  **/
  56 void hkl_interval_free(HklInterval *self)
  57 {
  58         if(self)
  59                 free(self);
  60 }
  61
  62 /**
  63  * hkl_interval_cmp: (skip)
  64  * @self:
  65  * @interval:
  66  *
  67  * compare two intervals
  68  *
  69  * Returns:
  70  **/
  71 int hkl_interval_cmp(const HklInterval *self, const HklInterval *interval)
  72 {
  73         return (fabs(self->min - interval->min) > HKL_EPSILON)
  74                 || (fabs(self->max - interval->max) > HKL_EPSILON);
  75 }
  76
  77 /**
  78  * hkl_interval_plus_interval: (skip)
  79  * @self:
  80  * @interval:
  81  *
  82  * add two ontervals
  83  **/
  84 void hkl_interval_plus_interval(HklInterval *self, const HklInterval *interval)
  85 {
  86         self->min += interval->min;
  87         self->max += interval->max;
  88 }
  89
  90 /**
  91  * hkl_interval_plus_double: (skip)
  92  * @self:
  93  * @d:
  94  *
  95  * add to an interval a double
  96  **/
  97 void hkl_interval_plus_double(HklInterval *self, double const d)
  98 {
  99         self->min += d;
 100         self->max += d;
 101 }
 102
 103
 104 /**
 105  * hkl_interval_minus_interval: (skip)
 106  * @self:
 107  * @interval:
 108  *
 109  * substract two #HklInterval
 110  **/
 111 void hkl_interval_minus_interval(HklInterval *self, const HklInterval *interval)
 112 {
 113         self->min -= interval->max;
 114         self->max -= interval->min;
 115 }
 116
 117
 118 /**
 119  * hkl_interval_minus_double: (skip)
 120  * @self:
 121  * @d:
 122  *
 123  * subst a double to an #HklInterval
 124  **/
 125 void hkl_interval_minus_double(HklInterval *self, double const d)
 126 {
 127         self->min -= d;
 128         self->max -= d;
 129 }
 130
 131 /**
 132  * hkl_interval_times_interval: (skip)
 133  * @self:
 134  * @interval:
 135  *
 136  * multiply two #HklInterval
 137  **/
 138 void hkl_interval_times_interval(HklInterval *self, const HklInterval *interval)
 139 {
 140         double min;
 141         double max;
 142         double m1 = self->min * interval->min;
 143         double m2 = self->min * interval->max;
 144         double m3 = self->max * interval->min;
 145         double m4 = self->max * interval->max;
 146
 147         min = m1;
 148         if (m2 < min)
 149                 min = m2;
 150         if (m3 < min)
 151                 min = m3;
 152         if (m4 < min)
 153                 min = m4;
 154
 155         max = m1;
 156         if (m2 > max)
 157                 max = m2;
 158         if (m3 > max)
 159                 max = m3;
 160         if (m4 > max)
 161                 max = m4;
 162
 163         self->min = min;
 164         self->max = max;
 165 }
 166
 167 /**
 168  * hkl_interval_times_double: (skip)
 169  * @self:
 170  * @d:
 171  *
 172  * multiply an #HklInterval by a double
 173  **/
 174 void hkl_interval_times_double(HklInterval *self, double const d)
 175 {
 176         double min;
 177         double max;
 178         if (d < 0) {
 179                 min = self->max * d;
 180                 max = self->min * d;
 181         } else {
 182                 min = self->min * d;
 183                 max = self->max * d;
 184         }
 185         self->min = min;
 186         self->max = max;
 187 }
 188
 189 /**
 190  * hkl_interval_divides_double: (skip)
 191  * @self:
 192  * @d:
 193  *
 194  * divide an #HklInterval by a double
 195  **/
 196 void hkl_interval_divides_double(HklInterval *self, double const d)
 197 {
 198         double min = self->min / d;
 199         double max = self->max / d;
 200         if (min > max){
 201                 double tmp = min;
 202                 min = max;
 203                 max = tmp;
 204         }
 205         self->min = min;
 206         self->max = max;
 207 }
 208
 209 /**
 210  * hkl_interval_contain_zero: (skip)
 211  * @self:
 212  *
 213  * check if an #HklInterval contain zero
 214  *
 215  * Returns:
 216  **/
 217 int hkl_interval_contain_zero(HklInterval const *self)
 218 {
 219         if (self->min <= 0 && self->max >= 0)
 220                 return TRUE;
 221         else
 222                 return FALSE;
 223 }
 224
 225 /**
 226  * hkl_interval_cos: (skip)
 227  * @self:
 228  *
 229  * compute the cosinus of an #HklInterval
 230  **/
 231 void hkl_interval_cos(HklInterval *self)
 232 {
 233         double min = 0;
 234         double max = 0;
 235         double cmin;
 236         double cmax;
 237
 238         cmin = cos(self->min);
 239         cmax = cos(self->max);
 240
 241         if (self->max - self->min >= 2 * M_PI) {
 242                 min = -1;
 243                 max = 1;
 244         } else {
 245                 int quad_min;
 246                 int quad_max;
 247
 248                 quad_min = (int)floor(self->min / M_PI_2) % 4;
 249                 if (quad_min < 0)
 250                         quad_min += 4;
 251
 252                 quad_max = (int)floor(self->max / M_PI_2) % 4;
 253                 if (quad_max < 0)
 254                         quad_max += 4;
 255
 256                 switch (quad_max) {
 257                 case 0:
 258                         switch (quad_min) {
 259                         case 0:
 260                                 min = cmax;
 261                                 max = cmin;
 262                                 break;
 263                         case 1:
 264                                 min = -1;
 265                                 max = 1;
 266                                 break;
 267                         case 2:
 268                                 min = cmin;
 269                                 max = 1;
 270                                 break;
 271                         case 3:
 272                                 if (cmin < cmax) {
 273                                         min = cmin;
 274                                         max = 1;
 275                                 } else {
 276                                         min = cmax;
 277                                         max = 1;
 278                                 }
 279                                 break;
 280                         }
 281                         break;
 282                 case 1:
 283                         switch (quad_min) {
 284                         case 0:
 285                                 min = cmax;
 286                                 max = cmin;
 287                                 break;
 288                         case 1:
 289                                 min = -1;
 290                                 max = 1;
 291                                 break;
 292                         case 2:
 293                                 if (cmin < cmax) {
 294                                         min = cmin;
 295                                         max = 1;
 296                                 } else {
 297                                         min = cmax;
 298                                         max = 1;
 299                                 }
 300                                 break;
 301                         case 3:
 302                                 min = cmax;
 303                                 max = 1;
 304                                 break;
 305                         }
 306                         break;
 307                 case 2:
 308                         switch (quad_min) {
 309                         case 0:
 310                                 min = -1;
 311                                 max = cmin;
 312                                 break;
 313                         case 1:
 314                                 if (cmin < cmax) {
 315                                         min = -1;
 316                                         max = cmax;
 317                                 } else {
 318                                         min = -1;
 319                                         max = cmin;
 320                                 }
 321                                 break;
 322                         case 2:
 323                                 if (cmin < cmax) {
 324                                         min = cmin;
 325                                         max = cmax;
 326                                 } else {
 327                                         min = -1;
 328                                         max = 1;
 329                                 }
 330                                 break;
 331                         case 3:
 332                                 min = -1;
 333                                 max = 1;
 334                                 break;
 335                         }
 336                         break;
 337                 case 3:
 338                         switch (quad_min) {
 339                         case 0:
 340                                 if (cmin < cmax) {
 341                                         min = -1;
 342                                         max = cmax;
 343                                 } else {
 344                                         min = -1;
 345                                         max = cmin;
 346                                 }
 347                                 break;
 348                         case 1:
 349                                 min = -1;
 350                                 max = cmax;
 351                                 break;
 352                         case 2:
 353                                 min = cmin;
 354                                 max = cmax;
 355                                 break;
 356                         case 3:
 357                                 if (cmin < cmax) {
 358                                         min = cmin;
 359                                         max = cmax;
 360                                 } else {
 361                                         min = -1;
 362                                         max = 1;
 363                                 }
 364                                 break;
 365                         }
 366                         break;
 367                 }
 368         }
 369         self->min = min;
 370         self->max = max;
 371 }
 372
 373 /**
 374  * hkl_interval_acos: (skip)
 375  * @self:
 376  *
 377  * compute the arc cosinus of an #HklInterval
 378  **/
 379 void hkl_interval_acos(HklInterval *self)
 380 {
 381         double tmp;
 382
 383         tmp = self->min;
 384         self->min = acos(self->max);
 385         self->max = acos(tmp);
 386 }
 387
 388
 389 /**
 390  * hkl_interval_sin: (skip)
 391  * @self:
 392  *
 393  * compute the sin of an #HklInterval
 394  **/
 395 void hkl_interval_sin(HklInterval *self)
 396 {
 397         double min = 0;
 398         double max = 0;
 399         double smin;
 400         double smax;
 401
 402         smin = sin(self->min);
 403         smax = sin(self->max);
 404
 405         /* if there is at least one period in b, then a = [-1, 1] */
 406         if ( self->max - self->min >= 2 * M_PI) {
 407                 min = -1;
 408                 max = 1;
 409         } else {
 410                 int quad_min;
 411                 int quad_max;
 412
 413                 quad_min = (int)floor(self->min / M_PI_2) % 4;
 414                 if (quad_min < 0)
 415                         quad_min += 4;
 416
 417                 quad_max = (int)floor(self->max / M_PI_2) % 4;
 418                 if (quad_max < 0)
 419                         quad_max += 4;
 420
 421                 switch (quad_max) {
 422                 case 0:
 423                         switch (quad_min) {
 424                         case 0:
 425                                 if (smin < smax) {
 426                                         min = smin;
 427                                         max = smax;
 428                                 } else {
 429                                         min = -1;
 430                                         max = 1;
 431                                 }
 432                                 break;
 433                         case 3:
 434                                 min = smin;
 435                                 max = smax;
 436                                 break;
 437                         case 1:
 438                                 if (smin > smax) {
 439                                         min = -1;
 440                                         max = smin;
 441                                 } else {
 442                                         min = -1;
 443                                         max = smax;
 444                                 }
 445                                 break;
 446                         case 2:
 447                                 min = -1;
 448                                 max = smax;
 449                                 break;
 450                         }
 451                         break;
 452                 case 1:
 453                         switch (quad_min) {
 454                         case 0:
 455                                 if (smin < smax) {
 456                                         min = smin;
 457                                         max = 1;
 458                                 } else {
 459                                         min = smax;
 460                                         max = 1;
 461                                 }
 462                                 break;
 463                         case 1:
 464                                 if (smin < smax) {
 465                                         min = -1;
 466                                         max = 1;
 467                                 } else {
 468                                         min = smax;
 469                                         max = smin;
 470                                 }
 471                                 break;
 472                         case 2:
 473                                 min = -1;
 474                                 max = 1;
 475                                 break;
 476                         case 3:
 477                                 min = smin;
 478                                 max = 1;
 479                                 break;
 480                         }
 481                         break;
 482                 case 2:
 483                         switch (quad_min) {
 484                         case 0:
 485                                 min = smax;
 486                                 max = 1;
 487                                 break;
 488                         case 1:
 489                         case 2:
 490                                 if (smin < smax) {
 491                                         min = -1;
 492                                         max = 1;
 493                                 } else {
 494                                         min = smax;
 495                                         max = smin;
 496                                 }
 497                                 break;
 498                         case 3:
 499                                 if (smin < smax) {
 500                                         min = smin;
 501                                         max = 1;
 502                                 } else {
 503                                         min = smax;
 504                                         max = 1;
 505                                 }
 506                                 break;
 507                         }
 508                         break;
 509                 case 3:
 510                         switch (quad_min) {
 511                         case 0:
 512                                 min = -1;
 513                                 max = 1;
 514                                 break;
 515                         case 1:
 516                                 min = -1;
 517                                 max = smin;
 518                                 break;
 519                         case 2:
 520                                 if (smin < smax) {
 521                                         min = -1;
 522                                         max = smax;
 523                                 } else {
 524                                         min = -1;
 525                                         max = smin;
 526                                 }
 527                                 break;
 528                         case 3:
 529                                 if (smin < smax) {
 530                                         min = smin;
 531                                         max = smax;
 532                                 } else {
 533                                         min = -1;
 534                                         max = 1;
 535                                 }
 536                                 break;
 537                         }
 538                         break;
 539                 }
 540         }
 541         self->min = min;
 542         self->max = max;
 543 }
 544
 545 /**
 546  * hkl_interval_asin: (skip)
 547  * @self:
 548  *
 549  * compute the arc sinus of an #HklInterval
 550  **/
 551 void hkl_interval_asin(HklInterval *self)
 552 {
 553         self->min = asin(self->min);
 554         self->max = asin(self->max);
 555 }
 556
 557 /**
 558  * hkl_interval_tan: (skip)
 559  * @self:
 560  *
 561  * compute the tangente of an #HklInterval
 562  **/
 563 void hkl_interval_tan(HklInterval *self)
 564 {
 565         int quadrant_down = (int)floor(self->min / M_PI_2);
 566         int quadrant_up = (int)floor(self->max / M_PI_2);
 567
 568         /* if there is at least one period in b or if b contains a Pi/2 + k*Pi, */
 569         /* then a = ]-oo, +oo[ */
 570         if ( ((quadrant_up - quadrant_down) >= 2)
 571              || (!(quadrant_down % 2) && (quadrant_up % 2)) ) {
 572                 self->min = -INFINITY;
 573                 self->max = INFINITY;
 574         } else {
 575                 self->min = tan(self->min);
 576                 self->max = tan(self->max);
 577         }
 578 }
 579
 580 /**
 581  * hkl_interval_atan: (skip)
 582  * @self:
 583  *
 584  * compute the arc tangente of an #HklInterval
 585  **/
 586 void hkl_interval_atan(HklInterval *self)
 587 {
 588         self->min = atan(self->min);
 589         self->max = atan(self->max);
 590 }
 591
 592 /**
 593  * hkl_interval_length: (skip)
 594  * @self:
 595  *
 596  * compute the length of an #HklInterval
 597  *
 598  * Returns:
 599  **/
 600 double hkl_interval_length(const HklInterval *self)
 601 {
 602         return self->max - self->min;
 603 }
 604
 605 /**
 606  * hkl_interval_angle_restrict_symm: (skip)
 607  * @self:
 608  *
 609  * restrict an #HklInterval into -pi, pi
 610  **/
 611 void hkl_interval_angle_restrict_symm(HklInterval *self)
 612 {
 613         gsl_sf_angle_restrict_symm_e(&self->min);
 614         gsl_sf_angle_restrict_symm_e(&self->max);
 615 }