function f=akin(n,x)
% Bickley function of order n >= 1 at x >= 0
% function f=akin(n,x)
% (c) 2008 Alain Hebert, Ecole Polytechnique de Montreal
   if n < 1 
      error('n must be > 0')
   elseif x < -1.0e-10
      error('x must be >= 0. x=%.5g.',x)
   end
   a=[10.584648, 6.545470, 1.976457, 0.333496, 0.035407, 0.002568, 0.000135, ...
      0.000005, 10.762666, 5.623335, 1.435437, 0.212504, 0.020365, 0.001360, ...
      0.000067, 0.000003, 0.0004364, -0.0003873, 0.0002740, -0.0001581, 0.0000763, ...
     -0.0000315, 0.0000113, -0.0000036, 0.0000010, -0.0000003, 0.0000001] ;
   b=[1.4534664,-0.2436620, 0.0258465, -0.0029653, 0.0005322, -0.0001499, ...
      0.0000560,-0.0000249, 0.0000125, -0.0000068, 0.0000040, -0.0000025, ...
      0.0000017,-0.0000012, 0.0000009, -0.0000008, 0.0000007,  0.3039967, ...
     -0.2136079, 0.0961280,-0.0324165,  0.0091054, -0.0023228, 0.0005813, ...
     -0.0001516, 0.0000426,-0.0000129,  0.0000042, -0.0000014, 0.0000005, ...
     -0.0000002, 0.0000001, 0.0031373, -0.0028564, 0.0021670, -0.0013874, ...
      0.0007615,-0.0003642, 0.0001541, -0.0000585, 0.0000202, -0.0000064, ...
      0.0000019,-0.0000005, 0.0000001] ;
   c=[0.4227843,  0.0534107, 0.08333333] ;
   d=[0.7853961, -0.9990226, 0.7266088, 0.7852024, -0.9912340, 0.6466375, ...
      0.7845986, -0.9791293, 0.5856605, 0.7834577, -0.9638914, 0.5346648, ...
      0.7817094, -0.9463843, 0.4907827, 0.7793031, -0.9271152, 0.4521752, ...
      0.7762107, -0.9064822, 0.4177388, 0.7724519, -0.8849865, 0.3869945, ...
      0.7679903, -0.8626685, 0.3590753, 0.7628988, -0.8400133, 0.3338676, ...
      0.7540982, -0.8054172, 0.2998569, 0.7401279, -0.7587821, 0.2609154, ...
      0.7239594, -0.7125290, 0.2278226, 0.7058777, -0.6672761, 0.1994999, ...
      0.6861762, -0.6234536, 0.1751248] ;
   e=[0.7247294, -0.7538355, 0.3203223, -5.337485E-2, 0.6663720, -0.6279752, ...
      0.2295280, -3.146833E-2, 0.5956163, -0.5094124, 0.1631667, -1.906198E-2, ...
      0.5191031, -0.4046007, 0.1152418, -1.174752E-2, 0.4425954, -0.3159648, ...
      8.097913E-2, -7.328415E-3, 0.3703178, -0.2434341,   5.669960E-2,-4.617254E-3, ...
      0.1684022,  -7.158569E-2, 7.923547E-3, 0.1278307,  -5.016344E-2, 5.095111E-3, ...
      9.611422E-2,-3.501524E-2, 3.286040E-3, 7.170491E-2,-2.437465E-2, 2.126242E-3, ...
      4.616317E-2,-1.425519E-2, 1.123687E-3, 2.475115E-2,-6.810124E-3, 4.762937E-4, ...
      1.302864E-2,-3.232035E-3, 2.031843E-4, 6.749972E-3,-1.524126E-3, 8.701440E-5, ...
      3.454768E-3,-7.157367E-4, 3.742673E-5, 1.749779E-3,-3.349194E-4, 1.615898E-5, ...
      7.522133E-4,-1.314173E-4, 5.777582E-6] ;
   inda=[3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 33, 36, 36, 39, 39, 42, 42, 45, 45] ;
   indb=[4, 8,12, 16, 20, 24, 27, 30, 33, 36, 39, 39, 42, 42, 45, 45, 48, 48, 51, ...
         51, 54, 54, 57, 57, 57] ;
   if n==1
     if x <= 0.0
       f=pi/2 ;
     elseif x < 0.04
       x2=x*x ; f=pi/2-x*((x2*c(2)+c(1))-log(0.5*x)*(x2*c(3)+1.0)) ;
     elseif x < 6.0
       u=0.11111111*x*x-2.0 ;
       i=8 ; imin=1 ;
       a0=a(i) ; a1=0.0 ; a2=0.0 ;
       while i > imin
         i=i-1 ;
         a2=a1 ; a1=a0 ; a0=u*a1-a2+a(i) ;
       end
       f1=0.5*(a0-a2) ;
       i=16 ; imin=9 ;
       a0=a(i) ; a1=0.0 ; a2=0.0 ;
       while i > imin
         i=i-1 ;
         a2=a1 ; a1=a0 ; a0=u*a1-a2+a(i) ;
       end
       kin=0.5*(a0-a2) ;
       f=pi/2-x*(f1-kin*log(0.5*x)) ;
     elseif x < 15.0
       u=0.44444444*x-4.6666667 ;
       i=27 ; imin=17 ;
       a0=a(i) ; a1=0.0 ; a2=0.0 ;
       while i > imin
         i=i-1 ;
         a2=a1 ; a1=a0 ; a0=u*a1-a2+a(i) ;
       end
       f=0.5*(a0-a2) ;
     else
       f=0 ;
     end
   elseif n==2
     if x <= 0.0
       f=1.0 ; return
     elseif x < 0.5
       u=8.0*x-2.0 ;
       i=17 ; imin=1 ;
     elseif x < 4.0
       u=1.1428571*x-2.5714286 ;
       i=32 ; imin=18 ;
     elseif x < 15.0
       u=0.36363636*x-3.4545455 ;
       i=45 ; imin=33 ;
     else
       f=0 ; return
     end
     b0=b(i) ; b1=0.0 ; b2=0.0 ;
     while i > imin
       i=i-1 ;
       b2=b1 ; b1=b0 ; b0=u*b1-b2+b(i) ;
     end
     f=0.5*(b0-b2) ;
   elseif n==3
     if x <= 0.0
       f=pi/4 ;
     elseif x < 1.0
       i=inda(fix(20.0*x)+1) ; f=x*(x*d(i)+d(i-1))+d(i-2) ;
     elseif x < 3.4
       i=indb(fix(2.5*(x-1.0))+1) ; f=x*(x*(x*e(i)+e(i-1))+e(i-2))+e(i-3) ;
     elseif x < 11.0
       i=indb(fix(2.5*(x-1.0))+1) ; f=x*(x*e(i)+e(i-1))+e(i-2) ;
     else
       f=0 ;
     end
   else
     ak1=akin(n-1,x) ; ak2=akin(n-2,x) ; ak3=akin(n-3,x) ;
     f=(ak3-ak1)*x/(n-1)+(n-2)*ak2/(n-1) ;
   end