diff options
Diffstat (limited to '2.3-1/src/c/specialFunctions/besselj/dbesselja.c')
| -rw-r--r-- | 2.3-1/src/c/specialFunctions/besselj/dbesselja.c | 143 |
1 files changed, 143 insertions, 0 deletions
diff --git a/2.3-1/src/c/specialFunctions/besselj/dbesselja.c b/2.3-1/src/c/specialFunctions/besselj/dbesselja.c new file mode 100644 index 00000000..23a355f5 --- /dev/null +++ b/2.3-1/src/c/specialFunctions/besselj/dbesselja.c @@ -0,0 +1,143 @@ +/* Copyright (C) 2017 - IIT Bombay - FOSSEE + + This file must be used under the terms of the CeCILL. + This source file is licensed as described in the file COPYING, which + you should have received as part of this distribution. The terms + are also available at + http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt + Organization: FOSSEE, IIT Bombay + Author: Brijesh Gupta C R + Email: toolbox@scilab.in +*/ +#include <stdio.h> +#include "math.h" +#include "besselj.h" + +#define ACC 40.0 +#define BIGNO 1.0e10 +#define BIGNI 1.0e-10 + + +static double bessj0( double x ) +{ + double ax,z; + double xx,y,ans,ans1,ans2; + + if ((ax=fabs(x)) < 8.0) { + y=x*x; + ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7 + +y*(-11214424.18+y*(77392.33017+y*(-184.9052456))))); + ans2=57568490411.0+y*(1029532985.0+y*(9494680.718 + +y*(59272.64853+y*(267.8532712+y*1.0)))); + ans=ans1/ans2; + } else { + z=8.0/ax; + y=z*z; + xx=ax-0.785398164; + ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4 + +y*(-0.2073370639e-5+y*0.2093887211e-6))); + ans2 = -0.1562499995e-1+y*(0.1430488765e-3 + +y*(-0.6911147651e-5+y*(0.7621095161e-6 + -y*0.934935152e-7))); + ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); + } + return ans; +} + + +static double bessj1( double x ) +{ + double ax,z; + double xx,y,ans,ans1,ans2; + + if ((ax=fabs(x)) < 8.0) { + y=x*x; + ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1 + +y*(-2972611.439+y*(15704.48260+y*(-30.16036606)))))); + ans2=144725228442.0+y*(2300535178.0+y*(18583304.74 + +y*(99447.43394+y*(376.9991397+y*1.0)))); + ans=ans1/ans2; + } else { + z=8.0/ax; + y=z*z; + xx=ax-2.356194491; + ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 + +y*(0.2457520174e-5+y*(-0.240337019e-6)))); + ans2=0.04687499995+y*(-0.2002690873e-3 + +y*(0.8449199096e-5+y*(-0.88228987e-6 + +y*0.105787412e-6))); + ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2); + if (x < 0.0) ans = -ans; + } + return ans; +} + + +double bessj( int n, double x ) +{ + int j, jsum, m; + double ax, bj, bjm, bjp, sum, tox, ans; + + + if (n < 0) + { + double dblank; + //setdblack_c( &dblank ); + return( dblank ); + } + ax=fabs(x); + if (n == 0) + return( bessj0(ax) ); + if (n == 1) + return( bessj1(ax) ); + + + if (ax == 0.0) + return 0.0; + else if (ax > (double) n) { + tox=2.0/ax; + bjm=bessj0(ax); + bj=bessj1(ax); + for (j=1;j<n;j++) { + bjp=j*tox*bj-bjm; + bjm=bj; + bj=bjp; + } + ans=bj; + } else { + tox=2.0/ax; + m=2*((n+(int) sqrt(ACC*n))/2); + jsum=0; + bjp=ans=sum=0.0; + bj=1.0; + for (j=m;j>0;j--) { + bjm=j*tox*bj-bjp; + bjp=bj; + bj=bjm; + if (fabs(bj) > BIGNO) { + bj *= BIGNI; + bjp *= BIGNI; + ans *= BIGNI; + sum *= BIGNI; + } + if (jsum) sum += bj; + jsum=!jsum; + if (j == n) ans=bjp; + } + sum=2.0*sum-bj; + ans /= sum; + } + return x < 0.0 && n%2 == 1 ? -ans : ans; +} + +void dbesselja(double* inp1,int size1, double* inp2,int size2, double* oup) +{ + int i; + if(size1 != size2) + printf("Error! arguments #1 and #2 have incompatible dimensions."); + for(i = 0; i<size1;i++) + { + oup[i] = bessj(inp1[i],inp2[i]); + } +} + |