diff --git a/src/conical.f90 b/src/conical.f90
index 81b2fee..fb117d2 100644
--- a/src/conical.f90
+++ b/src/conical.f90
@@ -1,17 +1,19 @@
 module conical
 
   use const_and_precisions, only : wp_
+
   implicit none
 
 contains
 
   function fconic(x,tau,m)
-!
-! this function subprogram computes the conical functions of the
-! first kind  P sub(-1/2 + i*tau) (x)  for m = 0 and m = 1.
-! Ref. in Kolbig, Comp. Phys. Commun. 23 (1981) 51
-!
+    ! Computes the conical functions of the first kind
+    ! P sub(-1/2 + i*tau) (x) for m = 0 and m = 1.
+    ! Ref. https://doi.org/10.1016/0010-4655(81)90129-6
+    use logger, only : log_error
+
     implicit none
+
     real(wp_), intent(in) :: x, tau
     integer, intent(in) :: m
     real(wp_) :: fconic
@@ -22,10 +24,10 @@ contains
     real(wp_), parameter :: rpi=1.7724538509055_wp_,pi2=0.63661977236758_wp_
     real(wp_), parameter :: eps=1.0e-14_wp_
     integer, parameter :: nmax=200
-!
+
     complex(wp_) a,b,c,ti,r,rr,q,u,u0,u1,u2,uu
     complex(wp_) v0,v1,v2,vv,w(19)
-!
+
     logical lm0,lm1,lta
 
     fconic=0.0_wp_
@@ -38,7 +40,7 @@ contains
     fm=m
     fconic=1.0_wp_-fm
     if(x == 1.0_wp_) return
-!
+
     fconic=0.0_wp_
     if(tau == 0.0_wp_ .and. abs(x-1.0_wp_) > 0.01_wp_) then
       if(x > 1.0_wp_) then
@@ -58,7 +60,7 @@ contains
       end if
     else
       ti=cmplx(0._wp_,tau,wp_)
-!
+
       if((-1._wp_ < x .and. x <= 0.0_wp_).or. &
          (0.0_wp_ < x .and. x <= 0.1_wp_ .and.tau<= 17.0_wp_).or. &
          (0.1_wp_ < x .and. x <= 0.2_wp_ .and.tau<=  5.0_wp_)) then
@@ -91,15 +93,15 @@ contains
           s=sqrt(x**2-1.0_wp_)
           t(1)=log(x+s)
           h(1)=tau*t(1)
-          b0=besj0l(h(1))
-          b1=besj1l(h(1))
+          b0=bessel_j0(h(1))
+          b1=bessel_j1(h(1))
           z=1.0_wp_
         else
           s=sqrt(1.0_wp_-x**2)
           t(1)=acos(x)
           h(1)=tau*t(1)
-          b0=besi0(h(1))
-          b1=besi1(h(1))
+          b0=bessel_i0(h(1))
+          b1=bessel_i1(h(1))
           z=-1.0_wp_
         end if
         h(1)=t(1)*x/s
@@ -111,7 +113,7 @@ contains
         do j = 2,11
           v(j)=v(j-1)*v(1)
         end do
-!
+
         if(lm1) then
           aa=-1.0_wp_
           a0=3.0_wp_*(1.0_wp_-h(1))/(8.0_wp_*t(1))
@@ -167,7 +169,7 @@ contains
         fconic=sqrt(t(1)/s)*(b0*s0+b1*s1)
         return
       end if
-!
+
       do
         if(lta) then
           y=-x1
@@ -189,14 +191,14 @@ contains
           w(14)=w(7)+6.0_wp_
           w(15)=2.0_wp_-w(8)
           w(16)=w(15)-2.0_wp_*w(7)
-!
+
           v0=1.0_wp_
           v1=1.0_wp_+(w(10)/(2.0_wp_*c))*y
           v2=w(12)+(w(10)*w(11)/(12.0_wp_*w(6)))*y2
           u0=1.0_wp_
           u1=v1-w(9)
           u2=v2-w(9)*w(12)+(w(8)*w(10)/(2.0_wp_*w(6)))*y2
-!
+
           r=1.0_wp_
           n=2
           do
@@ -228,7 +230,7 @@ contains
             w(10)=h(1)-c
             w(11)=w(1)*w(3)
             w(12)=w(5)*w(6)
-!
+
             w(17)=1.0_wp_+((h(9)+w(13)*fn+w(16))/(h(6)*w(5)))*y
             w(18)=-((w(11)*w(10)/h(6)+(h(9)-w(14)*fn+w(15))* &
                      w(11)*y/h(7))/w(12))*y
@@ -291,8 +293,12 @@ contains
 !
   end function fconic
 
+  !
+  ! Routines for conical function
+  !
+
   function clogam(z)
-!
+
     implicit none
     complex(wp_) :: clogam
     complex(wp_), intent(in) :: z
@@ -338,7 +344,7 @@ contains
      +r*(b(4)+r*(b(5)+r*(b(6)+r*(b(7)+r*(b(8)+r*(b(9)+r*b(10)))))))))) &
      /v-h
     if(x < 0.0_wp_) then
-!
+
       a=aint(x)-1.0_wp_
       c=pi*(x-a)
       d=pi*f
@@ -347,10 +353,10 @@ contains
       e=d+0.5_wp_*log(e*f**2+0.25_wp_*(1.0_wp_-e)**2)
       f=atan2(cos(c)*tanh(d),f)-a*pi
       clogam=1.1447298858494_wp_-cmplx(e,f,wp_)-clogam
-!
+
     end if
     if(t < 0.0_wp_) clogam=conjg(clogam)
-!
+
   end function clogam
 
   function ellick(xk)
@@ -380,7 +386,7 @@ contains
                3.3478943665762e-2_wp_, 2.6145014700314e-2_wp_, &
                1.6804023346363e-2_wp_, 6.4321465864383e-3_wp_, &
                9.8983328462254e-4_wp_, 3.1859195655502e-5_wp_/)
-!
+
     if(abs(xk) >= 1.0_wp_) then
       ellick=0.0_wp_
       return
@@ -398,9 +404,9 @@ contains
     pb=pb*eta
     ellick=1.3862943611199_wp_+pa-log(eta)*(0.5_wp_+pb)
     return
-!
+
     entry ellice(xk)
-!
+
     if (abs(xk) >= 1.0_wp_) then
       if (abs(xk) > 1.0_wp_) then
         ellick=0.0_wp_
@@ -423,210 +429,25 @@ contains
     ellick=1.0_wp_+pc-log(eta)*pd
   end function ellick
 
-  function besjy(x)
-    implicit none
-    real(wp_), intent(in) :: x
-    real(wp_) :: besjy,besj0l,besj1l
-    real(wp_) :: besy0,besy1
-    logical :: l
-    real(wp_) :: v,f,a,b,p,q
-!
-    entry besj0l(x)
-!
-    l=.true.
-    v=abs(x)
-    if(v >= 8.0_wp_) go to 4
-      8 f=0.0625_wp_*x**2-2.0_wp_
-    a =           - 0.0000000000000008_wp_
-    b = f * a     + 0.0000000000000413_wp_
-    a = f * b - a - 0.0000000000019438_wp_
-    b = f * a - b + 0.0000000000784870_wp_
-    a = f * b - a - 0.0000000026792535_wp_
-    b = f * a - b + 0.0000000760816359_wp_
-    a = f * b - a - 0.0000017619469078_wp_
-    b = f * a - b + 0.0000324603288210_wp_
-    a = f * b - a - 0.0004606261662063_wp_
-    b = f * a - b + 0.0048191800694676_wp_
-    a = f * b - a - 0.0348937694114089_wp_
-    b = f * a - b + 0.1580671023320973_wp_
-    a = f * b - a - 0.3700949938726498_wp_
-    b = f * a - b + 0.2651786132033368_wp_
-    a = f * b - a - 0.0087234423528522_wp_
-    a = f * a - b + 0.3154559429497802_wp_
-    besjy=0.5_wp_*(a-b)
-    if(l) return
-!
-    a =           + 0.0000000000000016_wp_
-    b = f * a     - 0.0000000000000875_wp_
-    a = f * b - a + 0.0000000000040263_wp_
-    b = f * a - b - 0.0000000001583755_wp_
-    a = f * b - a + 0.0000000052487948_wp_
-    b = f * a - b - 0.0000001440723327_wp_
-    a = f * b - a + 0.0000032065325377_wp_
-    b = f * a - b - 0.0000563207914106_wp_
-    a = f * b - a + 0.0007531135932578_wp_
-    b = f * a - b - 0.0072879624795521_wp_
-    a = f * b - a + 0.0471966895957634_wp_
-    b = f * a - b - 0.1773020127811436_wp_
-    a = f * b - a + 0.2615673462550466_wp_
-    b = f * a - b + 0.1790343140771827_wp_
-    a = f * b - a - 0.2744743055297453_wp_
-    a = f * a - b - 0.0662922264065699_wp_
-    besjy=0.636619772367581_wp_*log(x)*besjy+0.5_wp_*(a-b)
-    return
-!
-      4 f=256.0_wp_/x**2-2.0_wp_
-    b =           + 0.0000000000000007_wp_
-    a = f * b     - 0.0000000000000051_wp_
-    b = f * a - b + 0.0000000000000433_wp_
-    a = f * b - a - 0.0000000000004305_wp_
-    b = f * a - b + 0.0000000000051683_wp_
-    a = f * b - a - 0.0000000000786409_wp_
-    b = f * a - b + 0.0000000016306465_wp_
-    a = f * b - a - 0.0000000517059454_wp_
-    b = f * a - b + 0.0000030751847875_wp_
-    a = f * b - a - 0.0005365220468132_wp_
-    a = f * a - b + 1.9989206986950373_wp_
-    p=a-b
-    b =           - 0.0000000000000006_wp_
-    a = f * b     + 0.0000000000000043_wp_
-    b = f * a - b - 0.0000000000000334_wp_
-    a = f * b - a + 0.0000000000003006_wp_
-    b = f * a - b - 0.0000000000032067_wp_
-    a = f * b - a + 0.0000000000422012_wp_
-    b = f * a - b - 0.0000000007271916_wp_
-    a = f * b - a + 0.0000000179724572_wp_
-    b = f * a - b - 0.0000007414498411_wp_
-    a = f * b - a + 0.0000683851994261_wp_
-    a = f * a - b - 0.0311117092106740_wp_
-    q=8.0_wp_*(a-b)/v
-    f=v-0.785398163397448_wp_
-    a=cos(f)
-    b=sin(f)
-    f=0.398942280401432_wp_/sqrt(v)
-    if(l) go to 6
-    besjy=f*(q*a+p*b)
-    return
-      6 besjy=f*(p*a-q*b)
-    return
-!
-    entry besj1l(x)
-!
-    l=.true.
-    v=abs(x)
-    if(v >= 8.0_wp_) go to 5
-      3 f=0.0625_wp_*x**2-2.0_wp_
-    b =           + 0.0000000000000114_wp_
-    a = f * b     - 0.0000000000005777_wp_
-    b = f * a - b + 0.0000000000252812_wp_
-    a = f * b - a - 0.0000000009424213_wp_
-    b = f * a - b + 0.0000000294970701_wp_
-    a = f * b - a - 0.0000007617587805_wp_
-    b = f * a - b + 0.0000158870192399_wp_
-    a = f * b - a - 0.0002604443893486_wp_
-    b = f * a - b + 0.0032402701826839_wp_
-    a = f * b - a - 0.0291755248061542_wp_
-    b = f * a - b + 0.1777091172397283_wp_
-    a = f * b - a - 0.6614439341345433_wp_
-    b = f * a - b + 1.2879940988576776_wp_
-    a = f * b - a - 1.1918011605412169_wp_
-    a = f * a - b + 1.2967175412105298_wp_
-    besjy=0.0625_wp_*(a-b)*x
-    if(l) return
-!
-    b =           - 0.0000000000000244_wp_
-    a = f * b     + 0.0000000000012114_wp_
-    b = f * a - b - 0.0000000000517212_wp_
-    a = f * b - a + 0.0000000018754703_wp_
-    b = f * a - b - 0.0000000568844004_wp_
-    a = f * b - a + 0.0000014166243645_wp_
-    b = f * a - b - 0.0000283046401495_wp_
-    a = f * b - a + 0.0004404786298671_wp_
-    b = f * a - b - 0.0051316411610611_wp_
-    a = f * b - a + 0.0423191803533369_wp_
-    b = f * a - b - 0.2266249915567549_wp_
-    a = f * b - a + 0.6756157807721877_wp_
-    b = f * a - b - 0.7672963628866459_wp_
-    a = f * b - a - 0.1286973843813500_wp_
-    a = f * a - b + 0.0406082117718685_wp_
-    besjy=0.636619772367581_wp_*log(x)*besjy-0.636619772367581_wp_/x &
-         +0.0625_wp_*(a-b)*x
-    return
-!
-      5 f=256.0_wp_/x**2-2.0_wp_
-    b =           - 0.0000000000000007_wp_
-    a = f * b     + 0.0000000000000055_wp_
-    b = f * a - b - 0.0000000000000468_wp_
-    a = f * b - a + 0.0000000000004699_wp_
-    b = f * a - b - 0.0000000000057049_wp_
-    a = f * b - a + 0.0000000000881690_wp_
-    b = f * a - b - 0.0000000018718907_wp_
-    a = f * b - a + 0.0000000617763396_wp_
-    b = f * a - b - 0.0000039872843005_wp_
-    a = f * b - a + 0.0008989898330859_wp_
-    a = f * a - b + 2.0018060817200274_wp_
-    p=a-b
-    b =           + 0.0000000000000007_wp_
-    a = f * b     - 0.0000000000000046_wp_
-    b = f * a - b + 0.0000000000000360_wp_
-    a = f * b - a - 0.0000000000003264_wp_
-    b = f * a - b + 0.0000000000035152_wp_
-    a = f * b - a - 0.0000000000468636_wp_
-    b = f * a - b + 0.0000000008229193_wp_
-    a = f * b - a - 0.0000000209597814_wp_
-    b = f * a - b + 0.0000009138615258_wp_
-    a = f * b - a - 0.0000962772354916_wp_
-    a = f * a - b + 0.0935555741390707_wp_
-    q=8.0_wp_*(a-b)/v
-    f=v-2.356194490192345_wp_
-    a=cos(f)
-    b=sin(f)
-    f=0.398942280401432_wp_/sqrt(v)
-    if(l) go to 7
-    besjy=f*(q*a+p*b)
-    return
-      7 besjy=f*(p*a-q*b)
-    if(x < 0.0_wp_) besjy=-besjy
-    return
-!
-    entry besy0(x)
-!
-    if(x <= 0.0_wp_) go to 9
-    l=.false.
-    v=x
-    if(v >= 8.0_wp_) go to 4
-    go to 8
-    entry besy1(x)
-!
-    if(x <= 0.0_wp_) go to 9
-    l=.false.
-    v=x
-    if(v >= 8.0_wp_) go to 5
-    go to 3
-!
-      9 besjy=0.0_wp_
-    write(*,"(1x,'besjy ... non-positive argument x = ',e15.4)") x
-  end function besjy
+  function bessel_i(x)
+    ! Computes the modified Bessel functions I₀(x), I₁(x):
+    !	  entry bessel_i0(X)= I0(x)
+    !	  entry bessel_i1(X)= I1(x)
 
-  function besik(x)
     implicit none
+
     real(wp_), intent(in) :: x
-    real(wp_) :: besik,ebesi0,besi0,ebesi1,besi1,ebesk0,besk0,ebesk1,besk1
+    real(wp_) :: bessel_i,bessel_i0,bessel_i1
     logical :: l,e
-    real(wp_) :: v,f,a,b,z
-!
-    entry ebesi0(x)
-!
-    e=.true.
-    go to 1
-!
-    entry besi0(x)
-!
+    real(wp_) :: v,f,a,b
+
+    entry bessel_i0(x)
+
     e=.false.
-      1 l=.true.
+    l=.true.
     v=abs(x)
     if(v >= 8.0_wp_) go to 4
-      8 f=0.0625_wp_*x**2-2.0_wp_
+      f=0.0625_wp_*x**2-2.0_wp_
     a = 0.000000000000002_wp_
     b = f * a     +   0.000000000000120_wp_
     a = f * b - a +   0.000000000006097_wp_
@@ -643,10 +464,10 @@ contains
     b = f * a - b +  82.489032744024100_wp_
     a = f * b - a + 190.494320172742844_wp_
     a = f * a - b + 255.466879624362167_wp_
-    besik=0.5_wp_*(a-b)
-    if(l .and. e) besik=exp(-v)*besik
+    bessel_i=0.5_wp_*(a-b)
+    if(l .and. e) bessel_i=exp(-v)*bessel_i
     if(l) return
-!
+
     a =           +   0.000000000000003_wp_
     b = f * a     +   0.000000000000159_wp_
     a = f * b - a +   0.000000000007658_wp_
@@ -663,10 +484,10 @@ contains
     b = f * a - b +   8.005368868700334_wp_
     a = f * b - a -   4.563433586448395_wp_
     a = f * a - b -  21.057660177402440_wp_
-    besik=-log(0.125_wp_*x)*besik+0.5_wp_*(a-b)
-    if(e) besik=exp(x)*besik
+    bessel_i=-log(0.125_wp_*x)*bessel_i+0.5_wp_*(a-b)
+    if(e) bessel_i=exp(x)*bessel_i
     return
-!
+
       4 f=32.0_wp_/v-2.0_wp_
     b =           - 0.000000000000001_wp_
     a = f * b     - 0.000000000000001_wp_
@@ -689,23 +510,18 @@ contains
     b = f * a - b + 0.000172700630778_wp_
     a = f * b - a + 0.008445122624921_wp_
     a = f * a - b + 2.016558410917480_wp_
-    besik=0.199471140200717_wp_*(a-b)/sqrt(v)
+    bessel_i=0.199471140200717_wp_*(a-b)/sqrt(v)
     if(e) return
-    besik=exp(v)*besik
+    bessel_i=exp(v)*bessel_i
     return
-!
-    entry ebesi1(x)
-!
-    e=.true.
-    go to 2
-!
-    entry besi1(x)
-!
+
+    entry bessel_i1(x)
+
     e=.false.
-      2 l=.true.
+      l=.true.
     v=abs(x)
     if(v >= 8.0_wp_) go to 3
-      7 f=0.0625_wp_*x**2-2.0_wp_
+      f=0.0625_wp_*x**2-2.0_wp_
     a =           +   0.000000000000001_wp_
     b = f * a     +   0.000000000000031_wp_
     a = f * b - a +   0.000000000001679_wp_
@@ -722,10 +538,10 @@ contains
     b = f * a - b +  69.395917633734448_wp_
     a = f * b - a + 181.312616040570265_wp_
     a = f * a - b + 259.890237806477292_wp_
-    besik=0.0625_wp_*(a-b)*x
-    if(l .and. e) besik=exp(-v)*besik
+    bessel_i=0.0625_wp_*(a-b)*x
+    if(l .and. e) bessel_i=exp(-v)*bessel_i
     if(l) return
-!
+
     a =           +   0.000000000000001_wp_
     b = f * a     +   0.000000000000042_wp_
     a = f * b - a +   0.000000000002163_wp_
@@ -742,10 +558,10 @@ contains
     b = f * a - b +   9.361617831395389_wp_
     a = f * b - a -   1.839239224286199_wp_
     a = f * a - b -  26.688095480862668_wp_
-    besik=log(0.125_wp_*x)*besik+1.0_wp_/x-0.0625_wp_*(a-b)*x
-    if(e) besik=exp(x)*besik
+    bessel_i=log(0.125_wp_*x)*bessel_i+1.0_wp_/x-0.0625_wp_*(a-b)*x
+    if(e) bessel_i=exp(x)*bessel_i
     return
-!
+
       3 f=32.0_wp_/v-2.0_wp_
     b =           + 0.000000000000001_wp_
     a = f * b     + 0.000000000000001_wp_
@@ -768,83 +584,11 @@ contains
     b = f * a - b - 0.000277205360764_wp_
     a = f * b - a - 0.024467442963276_wp_
     a = f * a - b + 1.951601204652572_wp_
-    besik=0.199471140200717_wp_*(a-b)/sqrt(v)
-    if(x < 0.0_wp_) besik=-besik
+    bessel_i=0.199471140200717_wp_*(a-b)/sqrt(v)
+    if(x < 0.0_wp_) bessel_i=-bessel_i
     if(e) return
-    besik=exp(v)*besik
+    bessel_i=exp(v)*bessel_i
     return
-!
-    entry ebesk0 (x)
-!
-    e=.true.
-    go to 11
-!
-    entry besk0(x)
-!
-    e=.false.
-     11 if(x <= 0.0_wp_) go to 9
-    l=.false.
-    v=x
-    if(x < 5.0_wp_) go to 8
-    f=20.0_wp_/x-2.0_wp_
-    a =           - 0.000000000000002_wp_
-    b = f * a     + 0.000000000000011_wp_
-    a = f * b - a - 0.000000000000079_wp_
-    b = f * a - b + 0.000000000000581_wp_
-    a = f * b - a - 0.000000000004580_wp_
-    b = f * a - b + 0.000000000039044_wp_
-    a = f * b - a - 0.000000000364547_wp_
-    b = f * a - b + 0.000000003792996_wp_
-    a = f * b - a - 0.000000045047338_wp_
-    b = f * a - b + 0.000000632575109_wp_
-    a = f * b - a - 0.000011106685197_wp_
-    b = f * a - b + 0.000269532612763_wp_
-    a = f * b - a - 0.011310504646928_wp_
-    a = f * a - b + 1.976816348461652_wp_
-    besik=0.626657068657750_wp_*(a-b)/sqrt(x)
-    if(e) return
-    z=besik
-    besik=0.0_wp_
-    if(x < 180.0_wp_) besik=exp(-x)*z
-    return
-!
-    entry ebesk1(x)
-!
-    e=.true.
-    go to 12
-!
-    entry besk1(x)
-!
-    e=.false.
-     12 if(x <= 0.0_wp_) go to 9
-    l=.false.
-    v=x
-    if(x < 5.0_wp_) go to 7
-    f=20.0_wp_/x-2.0_wp_
-    a =           + 0.000000000000002_wp_
-    b = f * a     - 0.000000000000013_wp_
-    a = f * b - a + 0.000000000000089_wp_
-    b = f * a - b - 0.000000000000663_wp_
-    a = f * b - a + 0.000000000005288_wp_
-    b = f * a - b - 0.000000000045757_wp_
-    a = f * b - a + 0.000000000435417_wp_
-    b = f * a - b - 0.000000004645555_wp_
-    a = f * b - a + 0.000000057132218_wp_
-    b = f * a - b - 0.000000845172048_wp_
-    a = f * b - a + 0.000016185063810_wp_
-    b = f * a - b - 0.000468475028167_wp_
-    a = f * b - a + 0.035465291243331_wp_
-    a = f * a - b + 2.071901717544716_wp_
-    besik=0.626657068657750_wp_*(a-b)/sqrt(x)
-    if(e) return
-    z=besik
-    besik=0.0_wp_
-    if(x < 180.0_wp_) besik=exp(-x)*z
-    return
-      9 besik=0.0_wp_
-    write(*,"(1x,'besik ... non-positive argument x = ',e15.4)") x
-  end function besik
-!
-! routines for conical function: end
-!
-end module conical
\ No newline at end of file
+  end function bessel_i
+
+end module conical