4 Bessil functions added

24 files changed, 1553 insertions, 0 deletions

diff --git a/2.3-1/demos/Brijesh_Demos/test_besseli.sci b/2.3-1/demos/Brijesh_Demos/test_besseli.sci
new file mode 100644
index 00000000..c08142fd
--- /dev/null
+++ b/2.3-1/demos/Brijesh_Demos/test_besseli.sci
@@ -0,0 +1,12 @@
+function test_besseli
+	disp('Datatype: Double');
+	i1 = [1 2 3 4 5];
+	i2 = [6 7 8 9 10];
+	o1 = besseli(i1,i2);
+	disp(o1);
+	disp('Datatype: float');
+	i3 = float([1 2 3 4 5]);
+	i4 = float([6 7 8 9 10]);
+	o2 = besseli(i3,i4);
+	disp(o2);
+endfunction
diff --git a/2.3-1/demos/Brijesh_Demos/test_besselj.sci b/2.3-1/demos/Brijesh_Demos/test_besselj.sci
new file mode 100644
index 00000000..a9b92ed1
--- /dev/null
+++ b/2.3-1/demos/Brijesh_Demos/test_besselj.sci
@@ -0,0 +1,12 @@
+function test_besselj
+	disp('Datatype: Double');
+	i1 = [1 2 3 4 5];
+	i2 = [6 7 8 9 10];
+	o1 = besselj(i1,i2);
+	disp(o1);
+	disp('Datatype: float');
+	i3 = float([1 2 3 4 5]);
+	i4 = float([6 7 8 9 10]);
+	o2 = besselj(i3,i4);
+	disp(o2);
+endfunction
diff --git a/2.3-1/demos/Brijesh_Demos/test_besselk.sci b/2.3-1/demos/Brijesh_Demos/test_besselk.sci
new file mode 100644
index 00000000..7beface2
--- /dev/null
+++ b/2.3-1/demos/Brijesh_Demos/test_besselk.sci
@@ -0,0 +1,12 @@
+function test_besselk
+	disp('Datatype: Double');
+	i1 = [1 2 3 4 5];
+	i2 = [6 7 8 9 10];
+	o1 = besselk(i1,i2);
+	disp(o1);
+	disp('Datatype: float');
+	i3 = float([1 2 3 4 5]);
+	i4 = float([6 7 8 9 10]);
+	o2 = besselk(i3,i4);
+	disp(o2);
+endfunction
diff --git a/2.3-1/demos/Brijesh_Demos/test_bessely.sci b/2.3-1/demos/Brijesh_Demos/test_bessely.sci
new file mode 100644
index 00000000..e4801690
--- /dev/null
+++ b/2.3-1/demos/Brijesh_Demos/test_bessely.sci
@@ -0,0 +1,12 @@
+function test_bessely
+	disp('Datatype: Double');
+	i1 = [1 2 3 4 5];
+	i2 = [6 7 8 9 10];
+	o1 = bessely(i1,i2);
+	disp(o1);
+	disp('Datatype: float');
+	i3 = float([1 2 3 4 5]);
+	i4 = float([6 7 8 9 10]);
+	o2 = bessely(i3,i4);
+	disp(o2);
+endfunction
diff --git a/2.3-1/demos/Brijesh_Demos/test_window.sci b/2.3-1/demos/Brijesh_Demos/test_window.sci
new file mode 100644
index 00000000..740d56b6
--- /dev/null
+++ b/2.3-1/demos/Brijesh_Demos/test_window.sci
@@ -0,0 +1,5 @@
+function test_window
+	disp('Datatype: Double');
+	o1 = window('hn', 10);
+	disp(o1);
+endfunction
diff --git a/2.3-1/src/c/signalProcessing/includes/window.h b/2.3-1/src/c/signalProcessing/includes/window.h
new file mode 100644
index 00000000..0a06f22e
--- /dev/null
+++ b/2.3-1/src/c/signalProcessing/includes/window.h
@@ -0,0 +1,33 @@
+ /* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __WINDOW_H__
+#define __WINDOW_H__
+#include "types.h"
+#include "floatComplex.h"
+#include "doubleComplex.h"
+#include "uint8.h"
+#include "uint16.h"
+#include "int16.h"
+#include "factorial.h"
+#include "gamma.h"
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+void dwindowa(char* inp1, int size, double n, double* out);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__WINDOW_H__*/
diff --git a/2.3-1/src/c/signalProcessing/interfaces/int_window.h b/2.3-1/src/c/signalProcessing/interfaces/int_window.h
new file mode 100644
index 00000000..77a90a2c
--- /dev/null
+++ b/2.3-1/src/c/signalProcessing/interfaces/int_window.h
@@ -0,0 +1,27 @@
+/* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_WINDOW_H__
+#define __INT_WINDOW_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+#define g2d0windowd2(in1, size, in2, out)  		dwindowa(in1, size[0]*size[1], in2,  out)
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_WINDOW_H__*/
diff --git a/2.3-1/src/c/signalProcessing/window/dwindowa.c b/2.3-1/src/c/signalProcessing/window/dwindowa.c
new file mode 100644
index 00000000..316f7138
--- /dev/null
+++ b/2.3-1/src/c/signalProcessing/window/dwindowa.c
@@ -0,0 +1,67 @@
+/* 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 <string.h>
+#include     "math.h"
+#include     "window.h"
+#include "ones.h"
+#include "abs.h"
+#define PI 3.1415927
+
+void dwindowa(char* inp1, int size, double n, double* out)
+{
+
+    double no2, un[(int)n], xt[(int)n];
+
+    no2 = (n-1)/2;
+    for(int i = 0; i < n; i++)
+	xt[i] = -no2 + i;
+
+    donesa(un, 1, n);
+
+
+    char tr[] = "tr", re[] = "re", hm[] = "hm", hn[] = "hn";
+    double flagtr = 0;
+
+    if(strncmp(re, inp1,2) == 0)
+    {
+	for(int i = 0; i < n; i++)
+	    out[i] = un[i];
+    }
+
+
+
+
+    if(strncmp(tr, inp1,2) == 0)
+    {
+	for(int i = 0; i < n; i++)
+	    out[i] = un[i] - (2 * dabss(xt[i])) / (n+1);
+    }
+
+
+
+    if(strncmp(hm, inp1,2) == 0)
+    {
+	for(int i = 0; i < n; i++)
+	    out[i] = 0.54 * un[i] + 0.46 * cos(2*PI*xt[i]/(n-1));
+    }
+
+
+    if(strncmp(hn, inp1,2) == 0)
+    {
+	for(int i = 0; i < n; i++)
+	    out[i] = 0.5 * un[i] + 0.5 * cos(2*PI*xt[i]/(n-1));
+    }
+
+
+
+}
diff --git a/2.3-1/src/c/specialFunctions/besseli/dbesselia.c b/2.3-1/src/c/specialFunctions/besseli/dbesselia.c
new file mode 100644
index 00000000..14ac1ed3
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/besseli/dbesselia.c
@@ -0,0 +1,111 @@
+/* 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     "besseli.h"
+
+#define ACC 40.0
+#define BIGNO 1.0e10
+#define BIGNI 1.0e-10
+
+
+double onebessi( double inp1, double inp2)
+{
+   int j;
+   double bi,bim,bip,tox,result,s,res,t;
+
+
+   if (inp1 < 0)
+   {
+      double   dblank;
+      return( dblank );
+   }
+   if (inp1 == 0)
+      {
+            if ((s=fabs(inp2)) < 3.75)
+            {
+                t=inp2/3.75,t=t*t;
+                res=1.0+t*(3.5156229+t*(3.0899424+t*(1.2067492+t*(0.2659732+t*(0.360768e-1+t*0.45813e-2)))));
+            }
+            else
+            {
+                t=3.75/s;
+                res=(exp(s)/sqrt(s))*(0.39894228+t*(0.1328592e-1+t*(0.225319e-2+t*(-0.157565e-2+t*(0.916281e-2+t*(-0.2057706e-1+t*(0.2635537e-1+t*(-0.1647633e-1+t*0.392377e-2))))))));
+            }
+            return res;
+      }
+   if (inp1 == 1)
+   {
+        if ((s=fabs(inp2)) < 3.75)
+        {
+            t=inp2/3.75,t=t*t;
+            res=s*(0.5+t*(0.87890594+t*(0.51498869+t*(0.15084934+t*(0.2658733e-1+t*(0.301532e-2+t*0.32411e-3))))));
+        }
+        else
+        {
+            t=3.75/s;
+            res=0.2282967e-1+t*(-0.2895312e-1+t*(0.1787654e-1-t*0.420059e-2));
+            res=0.39894228+t*(-0.3988024e-1+t*(-0.362018e-2+t*(0.163801e-2+t*(-0.1031555e-1+t*res))));
+            res *= (exp(s)/sqrt(s));
+        }
+        return inp2 < 0.0 ? -res : res;
+   }
+
+
+   if (inp2 == 0.0)
+      return 0.0;
+   else
+   {
+      tox=2.0/fabs(inp2);
+      bip=result=0.0;
+      bi=1.0;
+      for (j=2*(inp1+(int) sqrt(ACC*inp1));j>0;j--)
+      {
+         bim=bip+j*tox*bi;
+         bip=bi;
+         bi=bim;
+         if (fabs(bi) > BIGNO)
+         {
+            result *= BIGNI;
+            bi *= BIGNI;
+            bip *= BIGNI;
+         }
+         if (j == inp1) result=bip;
+      }
+
+      if ((s=fabs(inp2)) < 3.75)
+      {
+            t=inp2/3.75,t=t*t;
+            res=1.0+t*(3.5156229+t*(3.0899424+t*(1.2067492+t*(0.2659732+t*(0.360768e-1+t*0.45813e-2)))));
+      }
+      else
+      {
+            t=3.75/s;
+            res=(exp(s)/sqrt(s))*(0.39894228+t*(0.1328592e-1+t*(0.225319e-2+t*(-0.157565e-2+t*(0.916281e-2+t*(-0.2057706e-1+t*(0.2635537e-1+t*(-0.1647633e-1+t*0.392377e-2))))))));
+      }
+
+      result *= res/bi;
+      return  inp2 < 0.0 && (int)inp1%2 == 1 ? -result : result;
+   }
+}
+
+void dbesselia(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] = onebessi(inp1[i],inp2[i]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/besseli/sbesselia.c b/2.3-1/src/c/specialFunctions/besseli/sbesselia.c
new file mode 100644
index 00000000..6e185e0b
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/besseli/sbesselia.c
@@ -0,0 +1,111 @@
+/* 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     "besseli.h"
+
+#define ACC 40.0
+#define BIGNO 1.0e10
+#define BIGNI 1.0e-10
+
+
+float fonebessi( float inp1, float inp2)
+{
+   int j;
+   float bi,bim,bip,tox,result,s,res,t;
+
+
+   if (inp1 < 0)
+   {
+      float   dblank;
+      return( dblank );
+   }
+   if (inp1 == 0)
+      {
+            if ((s=fabs(inp2)) < 3.75)
+            {
+                t=inp2/3.75,t=t*t;
+                res=1.0+t*(3.5156229+t*(3.0899424+t*(1.2067492+t*(0.2659732+t*(0.360768e-1+t*0.45813e-2)))));
+            }
+            else
+            {
+                t=3.75/s;
+                res=(exp(s)/sqrt(s))*(0.39894228+t*(0.1328592e-1+t*(0.225319e-2+t*(-0.157565e-2+t*(0.916281e-2+t*(-0.2057706e-1+t*(0.2635537e-1+t*(-0.1647633e-1+t*0.392377e-2))))))));
+            }
+            return res;
+      }
+   if (inp1 == 1)
+   {
+        if ((s=fabs(inp2)) < 3.75)
+        {
+            t=inp2/3.75,t=t*t;
+            res=s*(0.5+t*(0.87890594+t*(0.51498869+t*(0.15084934+t*(0.2658733e-1+t*(0.301532e-2+t*0.32411e-3))))));
+        }
+        else
+        {
+            t=3.75/s;
+            res=0.2282967e-1+t*(-0.2895312e-1+t*(0.1787654e-1-t*0.420059e-2));
+            res=0.39894228+t*(-0.3988024e-1+t*(-0.362018e-2+t*(0.163801e-2+t*(-0.1031555e-1+t*res))));
+            res *= (exp(s)/sqrt(s));
+        }
+        return inp2 < 0.0 ? -res : res;
+   }
+
+
+   if (inp2 == 0.0)
+      return 0.0;
+   else
+   {
+      tox=2.0/fabs(inp2);
+      bip=result=0.0;
+      bi=1.0;
+      for (j=2*(inp1+(int) sqrt(ACC*inp1));j>0;j--)
+      {
+         bim=bip+j*tox*bi;
+         bip=bi;
+         bi=bim;
+         if (fabs(bi) > BIGNO)
+         {
+            result *= BIGNI;
+            bi *= BIGNI;
+            bip *= BIGNI;
+         }
+         if (j == inp1) result=bip;
+      }
+
+      if ((s=fabs(inp2)) < 3.75)
+      {
+            t=inp2/3.75,t=t*t;
+            res=1.0+t*(3.5156229+t*(3.0899424+t*(1.2067492+t*(0.2659732+t*(0.360768e-1+t*0.45813e-2)))));
+      }
+      else
+      {
+            t=3.75/s;
+            res=(exp(s)/sqrt(s))*(0.39894228+t*(0.1328592e-1+t*(0.225319e-2+t*(-0.157565e-2+t*(0.916281e-2+t*(-0.2057706e-1+t*(0.2635537e-1+t*(-0.1647633e-1+t*0.392377e-2))))))));
+      }
+
+      result *= res/bi;
+      return  inp2 < 0.0 && (int)inp1%2 == 1 ? -result : result;
+   }
+}
+
+void sbesselia(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] = fonebessi(inp1[i],inp2[i]);
+    }
+}
+
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]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/besselj/sbesselja.c b/2.3-1/src/c/specialFunctions/besselj/sbesselja.c
new file mode 100644
index 00000000..68d2ea4f
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/besselj/sbesselja.c
@@ -0,0 +1,142 @@
+/* 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 float fbessj0( float x )
+{
+   float ax,z;
+   float 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 float fbessj1( float x )
+{
+   float ax,z;
+   float 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;
+}
+
+
+float fbessj( int n, float x )
+{
+   int    j, jsum, m;
+   float ax, bj, bjm, bjp, sum, tox, ans;
+
+
+   if (n < 0)
+   {
+      float   dblank;
+      return( dblank );
+   }
+   ax=fabs(x);
+   if (n == 0)
+      return( fbessj0(ax) );
+   if (n == 1)
+      return( fbessj1(ax) );
+
+
+   if (ax == 0.0)
+      return 0.0;
+   else if (ax > (float) n) {
+      tox=2.0/ax;
+      bjm=fbessj0(ax);
+      bj=fbessj1(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 sbesselja(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] = fbessj(inp1[i],inp2[i]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/besselk/dbesselka.c b/2.3-1/src/c/specialFunctions/besselk/dbesselka.c
new file mode 100644
index 00000000..df6c0708
--- /dev/null
+++ b/2.3-1/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/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]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/bessely/dbesselya.c b/2.3-1/src/c/specialFunctions/bessely/dbesselya.c
new file mode 100644
index 00000000..76678629
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/bessely/dbesselya.c
@@ -0,0 +1,167 @@
+/* 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     "bessely.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;
+}
+
+
+static double bessy0( double x )
+{
+   double z;
+   double xx,y,ans,ans1,ans2;
+
+   if (x < 8.0) {
+      y=x*x;
+      ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
+         +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
+      ans2=40076544269.0+y*(745249964.8+y*(7189466.438
+         +y*(47447.26470+y*(226.1030244+y*1.0))));
+      ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
+   } else {
+      z=8.0/x;
+      y=z*z;
+      xx=x-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.934945152e-7))));
+      ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
+   }
+   return ans;
+}
+
+static double bessy1( double x )
+{
+   double z;
+   double xx,y,ans,ans1,ans2;
+
+   if (x < 8.0) {
+      y=x*x;
+      ans1=x*(-0.4900604943e13+y*(0.1275274390e13
+         +y*(-0.5153438139e11+y*(0.7349264551e9
+         +y*(-0.4237922726e7+y*0.8511937935e4)))));
+      ans2=0.2499580570e14+y*(0.4244419664e12
+         +y*(0.3733650367e10+y*(0.2245904002e8
+         +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
+      ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
+   } else {
+      z=8.0/x;
+      y=z*z;
+      xx=x-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/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
+   }
+   return ans;
+}
+
+double bessy( int n, double x )
+{
+   int j;
+   double by,bym,byp,tox;
+
+
+   if (n < 0 || x == 0.0)
+   {
+      double   dblank;
+      return( dblank );
+   }
+   if (n == 0)
+      return( bessy0(x) );
+   if (n == 1)
+      return( bessy1(x) );
+
+   tox=2.0/x;
+   by=bessy1(x);
+   bym=bessy0(x);
+   for (j=1;j<n;j++) {
+      byp=j*tox*by-bym;
+      bym=by;
+      by=byp;
+   }
+   return by;
+}
+
+void dbesselya(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] = bessy(inp1[i],inp2[i]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/bessely/sbesselya.c b/2.3-1/src/c/specialFunctions/bessely/sbesselya.c
new file mode 100644
index 00000000..1b42736c
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/bessely/sbesselya.c
@@ -0,0 +1,168 @@
+/* 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     "bessely.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;
+}
+
+
+
+static float fbessy0( float x )
+{
+   float z;
+   float xx,y,ans,ans1,ans2;
+
+   if (x < 8.0) {
+      y=x*x;
+      ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
+         +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
+      ans2=40076544269.0+y*(745249964.8+y*(7189466.438
+         +y*(47447.26470+y*(226.1030244+y*1.0))));
+      ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
+   } else {
+      z=8.0/x;
+      y=z*z;
+      xx=x-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.934945152e-7))));
+      ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
+   }
+   return ans;
+}
+
+static float fbessy1( float x )
+{
+   float z;
+   float xx,y,ans,ans1,ans2;
+
+   if (x < 8.0) {
+      y=x*x;
+      ans1=x*(-0.4900604943e13+y*(0.1275274390e13
+         +y*(-0.5153438139e11+y*(0.7349264551e9
+         +y*(-0.4237922726e7+y*0.8511937935e4)))));
+      ans2=0.2499580570e14+y*(0.4244419664e12
+         +y*(0.3733650367e10+y*(0.2245904002e8
+         +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
+      ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
+   } else {
+      z=8.0/x;
+      y=z*z;
+      xx=x-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/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
+   }
+   return ans;
+}
+
+float fbessy( int n, float x )
+{
+   int j;
+   float by,bym,byp,tox;
+
+
+   if (n < 0 || x == 0.0)
+   {
+      float   dblank;
+      return( dblank );
+   }
+   if (n == 0)
+      return( fbessy0(x) );
+   if (n == 1)
+      return( fbessy1(x) );
+
+   tox=2.0/x;
+   by=fbessy1(x);
+   bym=fbessy0(x);
+   for (j=1;j<n;j++) {
+      byp=j*tox*by-bym;
+      bym=by;
+      by=byp;
+   }
+   return by;
+}
+
+void sbesselya(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] = fbessy(inp1[i],inp2[i]);
+    }
+}
+
diff --git a/2.3-1/src/c/specialFunctions/includes/besseli.h b/2.3-1/src/c/specialFunctions/includes/besseli.h
new file mode 100644
index 00000000..6502b183
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/includes/besseli.h
@@ -0,0 +1,34 @@
+ /* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __BESSELI_H__
+#define __BESSELI_H__
+#include "types.h"
+#include "floatComplex.h"
+#include "doubleComplex.h"
+#include "uint8.h"
+#include "uint16.h"
+#include "int16.h"
+#include "factorial.h"
+#include "gamma.h"
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+void sbesselia(float* inp1,int size1, float* inp2,int size2, float* oup);
+void dbesselia(double* inp1,int size1, double* inp2,int size2, double* oup);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__BESSELI_H__*/
diff --git a/2.3-1/src/c/specialFunctions/includes/besselj.h b/2.3-1/src/c/specialFunctions/includes/besselj.h
new file mode 100644
index 00000000..de849a36
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/includes/besselj.h
@@ -0,0 +1,34 @@
+ /* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __BESSELJ_H__
+#define __BESSELJ_H__
+#include "types.h"
+#include "floatComplex.h"
+#include "doubleComplex.h"
+#include "uint8.h"
+#include "uint16.h"
+#include "int16.h"
+#include "factorial.h"
+#include "gamma.h"
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+void sbesselja(float* inp1,int size1, float* inp2,int size2, float* oup);
+void dbesselja(double* inp1,int size1, double* inp2,int size2, double* oup);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__BESSELJ_H__*/
diff --git a/2.3-1/src/c/specialFunctions/includes/besselk.h b/2.3-1/src/c/specialFunctions/includes/besselk.h
new file mode 100644
index 00000000..0be3c34e
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/includes/besselk.h
@@ -0,0 +1,34 @@
+ /* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __BESSELK_H__
+#define __BESSELK_H__
+#include "types.h"
+#include "floatComplex.h"
+#include "doubleComplex.h"
+#include "uint8.h"
+#include "uint16.h"
+#include "int16.h"
+#include "factorial.h"
+#include "gamma.h"
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+void sbesselka(float* inp1,int size1, float* inp2,int size2, float* oup);
+void dbesselka(double* inp1,int size1, double* inp2,int size2, double* oup);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__BESSELK_H__*/
diff --git a/2.3-1/src/c/specialFunctions/includes/bessely.h b/2.3-1/src/c/specialFunctions/includes/bessely.h
new file mode 100644
index 00000000..fb385efd
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/includes/bessely.h
@@ -0,0 +1,35 @@
+ /* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __BESSELY_H__
+#define __BESSELY_H__
+#include "types.h"
+#include "floatComplex.h"
+#include "doubleComplex.h"
+#include "uint8.h"
+#include "uint16.h"
+#include "int16.h"
+#include "factorial.h"
+#include "gamma.h"
+#include     "besselj.h"
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+void sbesselya(float* inp1,int size1, float* inp2,int size2, float* oup);
+void dbesselya(double* inp1,int size1, double* inp2,int size2, double* oup);
+
+#ifdef  __cplusplus 
+} /* extern "C" */
+#endif
+
+#endif /*__BESSELY_H__*/
diff --git a/2.3-1/src/c/specialFunctions/interfaces/int_besseli.h b/2.3-1/src/c/specialFunctions/interfaces/int_besseli.h
new file mode 100644
index 00000000..e11c48ca
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/interfaces/int_besseli.h
@@ -0,0 +1,28 @@
+/* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_BESSELI_H__
+#define __INT_BESSELI_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+#define d2d2besselid2(in1, size1, in2, size2, out)  		dbesselia(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+#define s2s2besselis2(in1, size1, in2, size2, out)  		sbesselia(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_BESSELI_H__*/
diff --git a/2.3-1/src/c/specialFunctions/interfaces/int_besselj.h b/2.3-1/src/c/specialFunctions/interfaces/int_besselj.h
new file mode 100644
index 00000000..cc0d85ee
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/interfaces/int_besselj.h
@@ -0,0 +1,28 @@
+/* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_BESSELJ_H__
+#define __INT_BESSELJ_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+#define d2d2besseljd2(in1, size1, in2, size2, out)  		dbesselja(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+#define s2s2besseljs2(in1, size1, in2, size2, out)  		sbesselja(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_BESSELJ_H__*/
diff --git a/2.3-1/src/c/specialFunctions/interfaces/int_besselk.h b/2.3-1/src/c/specialFunctions/interfaces/int_besselk.h
new file mode 100644
index 00000000..16b9c4b9
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/interfaces/int_besselk.h
@@ -0,0 +1,28 @@
+/* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_BESSELK_H__
+#define __INT_BESSELK_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+#define d2d2besselkd2(in1, size1, in2, size2, out)  		dbesselka(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+#define s2s2besselks2(in1, size1, in2, size2, out)  		sbesselka(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_BESSELK_H__*/
diff --git a/2.3-1/src/c/specialFunctions/interfaces/int_bessely.h b/2.3-1/src/c/specialFunctions/interfaces/int_bessely.h
new file mode 100644
index 00000000..749c8767
--- /dev/null
+++ b/2.3-1/src/c/specialFunctions/interfaces/int_bessely.h
@@ -0,0 +1,28 @@
+/* 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
+ Author: Brijesh Gupta C R
+ Organization: FOSSEE, IIT Bombay
+ Email: toolbox@scilab.in
+ */
+#ifndef __INT_BESSELY_H__
+#define __INT_BESSELY_H__
+
+#ifdef  __cplusplus
+extern "C" {
+#endif
+
+
+#define d2d2besselyd2(in1, size1, in2, size2, out)  		dbesselya(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+#define s2s2besselys2(in1, size1, in2, size2, out)  		sbesselya(in1,size1[0]*size1[1], in2, size2[0]*size2[1], out)
+
+
+#ifdef  __cplusplus
+} /* extern "C" */
+#endif
+
+#endif /*__INT_BESSELY_H__*/