1 files changed, 141 insertions, 0 deletions
diff --git a/2.3-1/src/c/specialFunctions/besselk/sbesselka.c b/2.3-1/src/c/specialFunctions/besselk/sbesselka.c
new file mode 100644
index 00000000..bfb35eb1
--- /dev/null
+++ b/2.3-1/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]);
+    }
+}
+