src/mathstuff.c

   1 /*
   2  * Grace - GRaphing, Advanced Computation and Exploration of data
   3  *
   4  * Home page: http://plasma-gate.weizmann.ac.il/Grace/
   5  *
   6  * Copyright (c) 1991-1995 Paul J Turner, Portland, OR
   7  * Copyright (c) 1996-2001 Grace Development Team
   8  *
   9  * Maintained by Evgeny Stambulchik
  10  *
  11  *
  12  *                           All Rights Reserved
  13  *
  14  *    This program is free software; you can redistribute it and/or modify
  15  *    it under the terms of the GNU General Public License as published by
  16  *    the Free Software Foundation; either version 2 of the License, or
  17  *    (at your option) any later version.
  18  *
  19  *    This program is distributed in the hope that it will be useful,
  20  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  21  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  22  *    GNU General Public License for more details.
  23  *
  24  *    You should have received a copy of the GNU General Public License
  25  *    along with this program; if not, write to the Free Software
  26  *    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  27  */
  28
  29 /* Wrappers for some math functions */
  30
  31 #include <stdlib.h>
  32
  33 #include "protos.h"
  34 #include "defines.h"
  35 #include "utils.h"
  36 #include "cephes.h"
  37
  38 #include "mathstuff.h"
  39
  40 double ai_wrap(double x)
  41 {
  42     double retval, dummy1, dummy2, dummy3;
  43     (void) airy(x, &retval, &dummy1, &dummy2, &dummy3);
  44     return retval;
  45 }
  46
  47 double bi_wrap(double x)
  48 {
  49     double retval, dummy1, dummy2, dummy3;
  50     (void) airy(x, &dummy1, &dummy2, &retval, &dummy3);
  51     return retval;
  52 }
  53
  54 double ci_wrap(double x)
  55 {
  56     double retval, dummy1;
  57     (void) sici(x, &dummy1, &retval);
  58     return retval;
  59 }
  60
  61 double si_wrap(double x)
  62 {
  63     double retval, dummy1;
  64     (void) sici(x, &retval, &dummy1);
  65     return retval;
  66 }
  67
  68 double chi_wrap(double x)
  69 {
  70     double retval, dummy1;
  71     (void) shichi(x, &dummy1, &retval);
  72     return retval;
  73 }
  74
  75 double shi_wrap(double x)
  76 {
  77     double retval, dummy1;
  78     (void) shichi(x, &retval, &dummy1);
  79     return retval;
  80 }
  81
  82 double fresnlc_wrap(double x)
  83 {
  84     double retval, dummy1;
  85     (void) fresnl(x, &dummy1, &retval);
  86     return retval;
  87 }
  88
  89 double fresnls_wrap(double x)
  90 {
  91     double retval, dummy1;
  92     (void) fresnl(x, &retval, &dummy1);
  93     return retval;
  94 }
  95
  96 double iv_wrap(double v, double x)
  97 {
  98     double retval;
  99     if (v == 0) {
 100         retval = i0(x);
 101     } else if (v == 1) {
 102         retval = i1(x);
 103     } else {
 104         retval = iv(v, x);
 105     }
 106     return retval;
 107 }
 108
 109 double jv_wrap(double v, double x)
 110 {
 111     double retval;
 112     if (v == rint(v)) {
 113         retval = jn((int) v, x);
 114     } else {
 115         retval = jv(v, x);
 116     }
 117     return retval;
 118 }
 119
 120 double kn_wrap(int n, double x)
 121 {
 122     double retval;
 123     if (n == 0) {
 124         retval = k0(x);
 125     } else if (n == 1) {
 126         retval = k1(x);
 127     } else {
 128         retval = kn(n, x);
 129     }
 130     return retval;
 131 }
 132
 133 double yv_wrap(double v, double x)
 134 {
 135     double retval;
 136     if (v == rint(v)) {
 137         retval = yn((int) v, x);
 138     } else {
 139         retval = yv(v, x);
 140     }
 141     return retval;
 142 }
 143
 144 double sqr_wrap(double x)
 145 {
 146     return x*x;
 147 }
 148
 149 double max_wrap(double x, double y)
 150 {
 151     return MAX2(x, y);
 152 }
 153
 154 double min_wrap(double x, double y)
 155 {
 156     return MIN2(x, y);
 157 }
 158
 159 double irand_wrap(int x)
 160 {
 161     return (double) (lrand48() % x);
 162 }
 163
 164 double pi_const(void)
 165 {
 166     return M_PI;
 167 }
 168
 169 double deg_uconst(void)
 170 {
 171     return M_PI / 180.0;
 172 }
 173
 174 double rad_uconst(void)
 175 {
 176     return 1.0;
 177 }
 178
 179 #define C1 0.1978977093962766
 180 #define C2 0.1352915131768107
 181
 182 double rnorm(double mean, double sdev)
 183 {
 184     double u = drand48();
 185
 186     return mean + sdev * (pow(u, C2) - pow(1.0 - u, C2)) / C1;
 187 }
 188
 189 double fx(double x)
 190 {
 191     return 1.0 / sqrt(2.0 * M_PI) * exp(-x * x * 0.5);
 192 }
 193
 194
 195 /* Wrapper around the cephes ellpe() call which definition
 196    is different to the standard definition. See cephes/ellpe.c */
 197 double ellpe_wrap ( double x ) {
 198     double y;
 199
 200     y = ellpe(1.0 - x);
 201     return y;
 202 }
 203
 204
 205 /* Wrapper around the cephes ellpk() call which definition
 206    is different to the standard definition. See cephes/ellpk.c */
 207 double ellpk_wrap ( double x ) {
 208     double y;
 209
 210     y = ellpk(1.0 - x);
 211     return y;
 212 }