summaryrefslogtreecommitdiff
path: root/src/c/specialFunctions/besselk
diff options
context:
space:
mode:
Diffstat (limited to 'src/c/specialFunctions/besselk')
-rw-r--r--src/c/specialFunctions/besselk/dbesselka.c141
-rw-r--r--src/c/specialFunctions/besselk/sbesselka.c141
2 files changed, 282 insertions, 0 deletions
diff --git a/src/c/specialFunctions/besselk/dbesselka.c b/src/c/specialFunctions/besselk/dbesselka.c
new file mode 100644
index 0000000..df6c070
--- /dev/null
+++ b/src/c/specialFunctions/besselk/dbesselka.c
@@ -0,0 +1,141 @@
+/* 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 "besselk.h"
+
+#define ACC 40.0
+#define BIGNO 1.0e10
+#define BIGNI 1.0e-10
+
+static double bessi0( double x )
+{
+ double ax,ans;
+ double y;
+
+
+ if ((ax=fabs(x)) < 3.75) {
+ y=x/3.75,y=y*y;
+ ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
+ +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
+ } else {
+ y=3.75/ax;
+ ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
+ +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
+ +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
+ +y*0.392377e-2))))))));
+ }
+ return ans;
+}
+
+
+
+
+static double bessi1( double x)
+{
+ double ax,ans;
+ double y;
+
+
+ if ((ax=fabs(x)) < 3.75) {
+ y=x/3.75,y=y*y;
+ ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
+ +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
+ } else {
+ y=3.75/ax;
+ ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
+ -y*0.420059e-2));
+ ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
+ +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
+ ans *= (exp(ax)/sqrt(ax));
+ }
+ return x < 0.0 ? -ans : ans;
+}
+
+
+static double bessk0( double x )
+{
+ double y,ans;
+
+ if (x <= 2.0) {
+ y=x*x/4.0;
+ ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
+ +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
+ +y*(0.10750e-3+y*0.74e-5))))));
+ } else {
+ y=2.0/x;
+ ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
+ +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
+ +y*(-0.251540e-2+y*0.53208e-3))))));
+ }
+ return ans;
+}
+
+
+
+
+static double bessk1( double x )
+{
+ double y,ans;
+
+ if (x <= 2.0) {
+ y=x*x/4.0;
+ ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
+ +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
+ +y*(-0.110404e-2+y*(-0.4686e-4)))))));
+ } else {
+ y=2.0/x;
+ ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
+ +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
+ +y*(0.325614e-2+y*(-0.68245e-3)))))));
+ }
+ return ans;
+}
+
+double bessk( int n, double x )
+{
+ int j;
+ double bk,bkm,bkp,tox;
+
+
+ if (n < 0 || x == 0.0)
+ {
+ double dblank;
+ return( dblank );
+ }
+ if (n == 0)
+ return( bessk0(x) );
+ if (n == 1)
+ return( bessk1(x) );
+
+ tox=2.0/x;
+ bkm=bessk0(x);
+ bk=bessk1(x);
+ for (j=1;j<n;j++) {
+ bkp=bkm+j*tox*bk;
+ bkm=bk;
+ bk=bkp;
+ }
+ return bk;
+}
+
+void dbesselka(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] = bessk(inp1[i],inp2[i]);
+ }
+}
+
diff --git a/src/c/specialFunctions/besselk/sbesselka.c b/src/c/specialFunctions/besselk/sbesselka.c
new file mode 100644
index 0000000..bfb35eb
--- /dev/null
+++ b/src/c/specialFunctions/besselk/sbesselka.c
@@ -0,0 +1,141 @@
+/* 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 "besselk.h"
+
+#define ACC 40.0
+#define BIGNO 1.0e10
+#define BIGNI 1.0e-10
+
+static float fbessi0( float x )
+{
+ float ax,ans;
+ float y;
+
+
+ if ((ax=fabs(x)) < 3.75) {
+ y=x/3.75,y=y*y;
+ ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
+ +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
+ } else {
+ y=3.75/ax;
+ ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
+ +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
+ +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
+ +y*0.392377e-2))))))));
+ }
+ return ans;
+}
+
+
+
+
+static float fbessi1( float x)
+{
+ float ax,ans;
+ float y;
+
+
+ if ((ax=fabs(x)) < 3.75) {
+ y=x/3.75,y=y*y;
+ ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
+ +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
+ } else {
+ y=3.75/ax;
+ ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
+ -y*0.420059e-2));
+ ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
+ +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
+ ans *= (exp(ax)/sqrt(ax));
+ }
+ return x < 0.0 ? -ans : ans;
+}
+
+
+static float fbessk0( float x )
+{
+ float y,ans;
+
+ if (x <= 2.0) {
+ y=x*x/4.0;
+ ans=(-log(x/2.0)*fbessi0(x))+(-0.57721566+y*(0.42278420
+ +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
+ +y*(0.10750e-3+y*0.74e-5))))));
+ } else {
+ y=2.0/x;
+ ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
+ +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
+ +y*(-0.251540e-2+y*0.53208e-3))))));
+ }
+ return ans;
+}
+
+
+
+
+static float fbessk1( float x )
+{
+ float y,ans;
+
+ if (x <= 2.0) {
+ y=x*x/4.0;
+ ans=(log(x/2.0)*fbessi1(x))+(1.0/x)*(1.0+y*(0.15443144
+ +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
+ +y*(-0.110404e-2+y*(-0.4686e-4)))))));
+ } else {
+ y=2.0/x;
+ ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
+ +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
+ +y*(0.325614e-2+y*(-0.68245e-3)))))));
+ }
+ return ans;
+}
+
+float fbessk( int n, float x )
+{
+ int j;
+ float bk,bkm,bkp,tox;
+
+
+ if (n < 0 || x == 0.0)
+ {
+ float dblank;
+ return( dblank );
+ }
+ if (n == 0)
+ return( fbessk0(x) );
+ if (n == 1)
+ return( fbessk1(x) );
+
+ tox=2.0/x;
+ bkm=fbessk0(x);
+ bk=fbessk1(x);
+ for (j=1;j<n;j++) {
+ bkp=bkm+j*tox*bk;
+ bkm=bk;
+ bk=bkp;
+ }
+ return bk;
+}
+
+void sbesselka(float* inp1,int size1, float* inp2,int size2, float* oup)
+{
+ int i;
+ if(size1 != size2)
+ printf("Error! arguments #1 and #2 have incompatible dimensions.");
+ for(i = 0; i<size1;i++)
+ {
+ oup[i] = fbessk(inp1[i],inp2[i]);
+ }
+}
+
generated by cgit (git 2.34.1) at 2024-12-19 19:16:17 +0000