# Functions for the ANSI version of Agena (this applies to the OS/2, Haiku, and DOS editions)

#####################################################################################################
#  BesselJY complex aux functions                                                                   #
#    taken from the FORTRAN mcjyna.for file written by Shanjie Zhang and Jianming Jin,              #
#    Computation of Special Functions, Copyright 1996 by John Wiley & Sons, Inc.,                   #
#    by kind permission of Jianming Jin                                                             #
#####################################################################################################

math.bessel_envj := << (n :: number, x :: number) -> 0.5*log10(6.28*n) - n*log10(1.36*x/n) >>;

# ===========================================================
# Function: msta1
# Purpose:  Determine the starting point for backward
#           recurrence such that the magnitude of
#           Jn(x) at that point is about 10^(-mp)
# Input:    x     - Argument of Jn(x)
#           mp    - Value of magnitude
# Output:   msta1 - Starting point
# ===========================================================

math.bessel_msta1 := proc(x :: number, mp :: number) is
   local a0, f, f0, f1, n0, n1, nn, it;
   a0 := abs(x);
   n0 := int(1.1*a0) + 1;
   f0 := math.bessel_envj(n0, a0) - mp;
   n1 := n0 + 5;
   f1 := math.bessel_envj(n1, a0) - mp;
   for it to 20 do
      nn := n1 - (n1-n0)/(1.0 - f0/f1);
      f := math.bessel_envj(nn,a0) - mp;
      if abs(nn-n1) < 1 then break fi;
      n0 := n1;
      f0 := f1;
      n1 := nn;
      f1 := f
   od;
   return int(nn)
end;

# ===========================================================
# Function: msta2
# Purpose: Determine the starting point for backward
#          recurrence such that all Jn(x) has mp
#          significant digits
# Input :  x     - Argument of Jn(x)
#          n     - Order of Jn(x)
#          mp    - Significant digit
# Output:  msta2 - Starting point
# ===========================================================

math.bessel_msta2 := proc(x :: number, n :: number, mp :: number) is
   local a0, hmp, ejn, obj, n0, n1, nn, f, f0, f1, it;
   a0 := abs(x);
   hmp := 0.5 * mp;
   ejn := math.bessel_envj(n, a0);
   if (ejn <= hmp) then
      obj := mp;
      n0 := int(1.1*a0)
   else
      obj := hmp + ejn;
      n0 := n
   fi;
   f0 := math.bessel_envj(n0, a0) - obj;
   n1 := n0 + 5;
   f1 := math.bessel_envj(n1, a0) - obj;
   for it to 20 do
      nn := n1 - (n1-n0)/(1.0-f0/f1);
      f := math.bessel_envj(nn, a0) - obj;
      if abs(nn-n1) < 1 then break fi;
      n0 := n1;
      f0 := f1;
      n1 := nn;
      f1 := f
   od;
   return int(nn) + 10
end;

# ===========================================================
# Purpose: Compute complex Bessel functions J0(z), J1(z)
#          Y0(z), Y1(z), and their derivatives
# Input :  z    - Complex argument
# Output:  cbj0 - J0(z)
#          cdj0 - J0'(z)
#          cbj1 - J1(z)
#          cdj1 - J1'(z)
#          cby0 - Y0(z)
#          cdy0 - Y0'(z)
#          cby1 - Y1(z)
#          cdy1 - Y1'(z)
# ===========================================================

environ.aux.cjy01_a := seq(
   -.703125e-01, .112152099609375e+00,
   -.5725014209747314e+00, .6074042001273483e+01,
   -.1100171402692467e+03, .3038090510922384e+04,
   -.1188384262567832e+06, .6252951493434797e+07,
   -.4259392165047669e+09, .3646840080706556e+11,
   -.3833534661393944e+13, .4854014686852901e+15);

environ.aux.cjy01_b := seq(
   .732421875e-01, -.2271080017089844e+00,
   .1727727502584457e+01, -.2438052969955606e+02,
   .5513358961220206e+03, -.1825775547429318e+05,
   .8328593040162893e+06, -.5006958953198893e+08,
   .3836255180230433e+10, -.3649010818849833e+12,
   .4218971570284096e+14, -.5827244631566907e+16);

environ.aux.cjy01_a1 := seq(
   .1171875e+00,-.144195556640625e+00,
   .6765925884246826e+00, -.6883914268109947e+01,
   .1215978918765359e+03, -.3302272294480852e+04,
   .1276412726461746e+06, -.6656367718817688e+07,
   .4502786003050393e+09, -.3833857520742790e+11,
   .4011838599133198e+13, -.5060568503314727e+15);

environ.aux.cjy01_b1 := seq(
   -.1025390625e+00,.2775764465332031e+00,
   -.1993531733751297e+01, .2724882731126854e+02,
   -.6038440767050702e+03, .1971837591223663e+05,
   -.8902978767070678e+06, .5310411010968522e+08,
   -.4043620325107754e+10, .3827011346598605e+12,
   -.4406481417852278e+14, .6065091351222699e+16);

cjy01 := proc(z) is
   local a0, k0, rp2, w0, w1, cr, ci, cp, cp0, cp1, cq0, cq1, cs, ct1, ct2, cu, z1, z2, k,
         cbj0, cbj1, cby0, cby1;
   rp2 := 2.0/Pi;
   ci := 0!1;
   a0 := abs(z);
   z2 := z*z;
   z1 := z;
   if a0 = 0.0 then
      cbj0 := 1!0;
      cbj1 := 0!0;
      cby0 := -infinity;
      cby1 := -infinity;
      return cbj0, cbj1, cby0, cby1
   fi;
   if real(z) < 0.0 then z1 := -z fi;
   if a0 <= 12.0 then
      cbj0 := 1!0;
      cr := 1!0;
      for k to 40 while abs(cr/cbj0) > 1.0e-15 do
         cr := -0.25*cr*z2/(k*k);
         cbj0 := cbj0+cr
      od;
      cbj1 := 1!0;
      cr := 1!0;
      for k to 40 while abs(cr/(cbj1)) > 1.0e-15 do
         cr := -0.25*cr*z2/(k*(k+1.0));
         cbj1 := cbj1+cr
      od;
      cbj1 := 0.5*z1*(cbj1);
      w0 := 0.0;
      cr := 1!0;
      cs := 0!0;
      for k to 40 do
        w0 := w0+1.0/k;
        cr := -0.25*cr/(k*k) * z2;
        cp := cr*w0;
        cs := cs+cp;
        if abs(cp/cs) < 1.0e-15 then break fi
      od;
      cby0 := rp2*(ln(z1/2.0)+EulerGamma)*cbj0 - rp2*cs;
      w1 := 0.0;
      cr := 1!0;
      cs := 1!0;
      for k to 40 do
         w1 := w1+1.0/k;
         cr := -0.25*cr/(k*(k+1)) * z2;
         cp := cr*(2.0*w1+1.0/(k+1.0));
         cs := cs+cp;
         if abs(cp/cs) < 1.0e-15 then break fi
      od;
      cby1 := rp2*((ln(z1/2.0)+EulerGamma)*cbj1 - 1.0/z1 - .25*z1*cs);
   else
      k0 := 12;
      if a0 >= 35.0 then k0 := 10 fi;
      if a0 >= 50.0 then k0 := 8 fi;
      ct1 := z1 - 0.25*Pi;
      cp0 := 1!0;
      for k to k0 do
         cp0 := cp0 + environ.aux.cjy01_a[k] * z1^(-2*k)
      od;
      cq0 := -0.125/z1;
      for k to k0 do
         cq0 := cq0 + environ.aux.cjy01_b[k] * z1^(-2*k-1)
      od;
      cu := sqrt(rp2/z1);
      cbj0 := cu * (cp0*cos(ct1) - cq0*sin(ct1));
      cby0 := cu * (cp0*sin(ct1) + cq0*cos(ct1));
      ct2 := z1-0.75*Pi;
      cp1 := 1!0;
      for k to k0 do
         cp1 := cp1 + environ.aux.cjy01_a1[k] * z1^(-2*k)
      od;
      cq1 := 0.375/z1;
      for k to k0 do
         cq1 := cq1 + environ.aux.cjy01_b1[k] * z1^(-2*k-1)
      od;
      cbj1 := cu * (cp1*cos(ct2) - cq1*sin(ct2));
      cby1 := cu * (cp1*sin(ct2) + cq1*cos(ct2));
   fi;
   if real(z) < 0.0 then
      if imag(z) < 0.0 then cby0 := cby0 - 2.0*ci*(cbj0) fi;
      if imag(z) > 0.0 then cby0 := cby0 + 2.0*ci*(cbj0) fi;
      if imag(z) < 0.0 then cby1 := -(cby1 - 2.0*ci*(cbj1)) fi;
      if imag(z) > 0.0 then cby1 := -(cby1 + 2.0*ci*(cbj1)) fi;
      cbj1 := -cbj1
   fi;
   return cbj0, cbj1, cby0, cby1;
end;


# ==========================================================
# Purpose: Compute Bessel functions Jn(z), Yn(z) and
#          their derivatives for a complex argument
# Input :  z      - Complex argument of Jn(z) and Yn(z)
#          n      - Order of Jn(z) and Yn(z)
# Output:  cbj(n) - Jn(z)
#          cdj(n) - Jn'(z)
#          cby(n) - Yn(z)
#          cdy(n) - Yn'(z)
#          nm     - Highest order computed
# Routines called:
#          (1) CJY01 to calculate J0(z), J1(z), Y0(z), Y1(z)
#          (2) msta1 and msta2 to calculate the starting
#              point for backward recurrence
# ==========================================================

math.bessel_cjyna := proc(n :: number, z :: complex, mode :: number) is
   local cg0, cg1, cbj0, cbj1, cj0, cj1, cjk, cby0, cby1, cf, cf1, cf2, cs,
         cp11, cp12, cp21, cp22, cyk, cylk, cyl1, cyl2, ch0, ch1, ch2,
         a0, ya0, ya1, yak, m, k, lb, lb0, nm;
   create local table cbj(n+1), table cby(n+1);
   lb0 := 0;
   a0 := abs(z);
   nm := n;
   if a0 < 1.0e-100 then
      for k from 0 to n do
         cbj[k] := 0!0;
         cby[k] := -infinity
      od;
      cbj[0] := 1!0
   else
      cbj0, cbj1, cby0, cby1 := cjy01(z)
      cbj[0] := cbj0;
      cbj[1] := cbj1;
      cby[0] := cby0;
      cby[1] := cby1;
      if n <= 1 then
         return is mode = 0 then cbj[n] else cby[n] si
      fi;
      if n < int(0.25*a0) then
         cj0 := cbj0;
         cj1 := cbj1;
         for k from 2 to n do
            cjk := 2.0*(k-1.0)/z*cj1 - cj0;
            cbj[k] := cjk;
            cj0 := cj1;
            cj1 := cjk
         od
      else
         m := math.bessel_msta1(a0, 200);
         if m < n then
            nm := m
         else
            m := math.bessel_msta2(a0, n, 15)
         fi;
         cf2 := 0!0;
         cf1 := 1.0e-100!0;
         for k from m to 0 by -1 do
            cf := 2.0*(k+1.0)/z*cf1 - cf2;
            if k <= nm then cbj[k] := cf fi;
            cf2 := cf1;
            cf1 := cf
         od;
         if abs(cbj0) > abs(cbj1) then
            cs := cbj0/cf
         else
            cs := cbj1/cf2
         fi;
         for k from 0 to nm do
            cbj[k] := cs*cbj[k]
         od
      fi;
      ya0 := abs(cby0);
      lb := 0;
      cg0 := cby0;
      cg1 := cby1;
      for k from 2 to nm do
         cyk := 2.0*(k-1.0)/z*cg1 - cg0;
         if abs(cyk) > 1.0e+290 then skip fi;
         yak := abs(cyk);
         ya1 := abs(cg0);
         if yak < ya0 and yak < ya1 then lb := k fi;
         cby[k] := cyk;
         cg0 := cg1;
         cg1 := cyk
      od;
      if not(lb <= 4 or imag(z) = 0.0) then
         while lb <> lb0 do
            ch2 := 1.0+0.0*I;
            ch1 := 0.0+0.0*I;
            lb0 := lb;
            for k from lb to 1 by -1 do
               ch0 := 2.0*k/z*ch1-ch2;
               ch2 := ch1;
               ch1 := ch0
            od;
            cp12 := ch0;
            cp22 := ch2;
            ch2 := 0.0+0.0*I;
            ch1 := 1.0+0.0*I;
            for k from lb to 1 by -1 do
               ch0 := 2.0*k/z*ch1-ch2;
               ch2 := ch1;
               ch1 := ch0
            od;
            cp11 := ch0;
            cp21 := ch2;
            if lb = nm then
               cbj[lb+1] := 2.0*lb/z*cbj[lb]-cbj[lb-1]
            fi;
            if abs(cbj[0]) > abs(cbj[1]) then
               cby[lb+1] := (cbj[lb+1]*cby0-2.0*cp11/(Pi*z))/cbj[0];
               cby[lb] := (cbj[lb]*cby0+2.0*cp12/(Pi*z))/cbj[0]
            else
               cby[lb+1] := (cbj[lb+1]*cby1-2.0*cp21/(Pi*z))/cbj[1];
               cby[lb] := (cbj[lb]*cby1+2.0*cp22/(Pi*z))/cbj[1]
            fi;
            cyl2 := cby[lb+1];
            cyl1 := cby[lb];
            for k from lb-1 to 0 by -1 do
               cylk := 2.0*(k+1.0)/z*cyl1-cyl2;
               cby[k] := cylk;
               cyl2 := cyl1;
               cyl1 := cylk
            od;
            cyl1 := cby[lb];
            cyl2 := cby[lb+1];
            for k from lb+1 to nm-1 do
               cylk := 2.0*k/z*cyl2-cyl1;
               cby[k+1] := cylk;
               cyl1 := cyl2;
               cyl2 := cylk
            od;
            for k from 2 to nm do
               if abs(cby[k]) < abs(cby[k-1]) then lb := k fi;
            od
         od
      fi
   fi;
   return is mode = 0 then cbj[n] else cby[n] si;
end;


local oldbesselj := besselj;  # save C library function for rational values

besselj := proc(n :: number, z) is
   if n < 0 then
      error('Error in `besselj`: order must be nonnegative.')
   elif float(n)
      then n := int(n)
   fi;
   case type(z)
     of number then
        return oldbesselj(n, z)
     of complex then
        return math.bessel_cjyna(n, z, 0)
     else
        argerror(z, 'besselj', 'number or complex value expected')
   esac
end;


local oldbessely := bessely;  # save C library function for rational values

bessely := proc(n :: number, z) is
   if n < 0 then
      error('Error in `bessely`: order must be nonnegative.')
   elif float(n)
      then n := int(n)
   fi;
   case type(z)
     of number then
        return oldbessely(n, z)
     of complex then
        return math.bessel_cjyna(n, z, 1)
     else
        argerror(z, 'bessely', 'number or complex value expected')
   esac
end;


local old_gamma := gamma;

gamma := proc(x) is
   if x :: number then
     return old_gamma(x)
   elif x :: complex then
     return exp(lngamma(x))
   else
     argerror(x, 'gamma', 'number or complex value expected')
   fi
end;

