////////////////////////////////////////////////////////////////
//
//     FloatMathieuSe.q
//
////////////////////////////////////////////////////////////////

(Def (_FloatMathieuSe nu x q)
  (Local answer)
  (Setq x (Float x))
  (Setq q (Float q))
  (Cond ((And (IsInteger nu) (IsFloat x) (IsFloat q))
         (Setq answer (_FloatMathieuSeCalc nu x q)))
        (TRUE
         (Setq answer `(MathieuSe ,nu ,x ,q))))
  (Return answer)
)

(Def (_FloatMathieuSeCalc nu x q)
  (Local b Bs answer)
  (WithMorePrecision (GetPrecision)
    (Setq answer 0.)
    (Setq b (MathieuB nu q))
    (Cond ((IsEven nu)
           (Setq Bs (_FloatMathieuOddPeriodPiSolution b q))
           (Loop (To k 0 (- (Length Bs) 1))
             (Do (Unless (IsOdd k)
                   (Setq answer (+ answer (* (Elt Bs k) (Sin (* k x)))))))))
          (TRUE
           (Setq Bs (_FloatMathieuOddPeriodTwoPiSolution b q))
           (Loop (To k 0 (- (Length Bs) 1))
             (Do (Unless (IsEven k)
                   (Setq answer (+ answer (* (Elt Bs k) (Sin (* k x)))))))))))
  (Setq answer (Float answer))
  (Return answer)
)

(Def (_FloatMathieuOddPeriodPiSolution b q)
  ////////////////////////////////////////////////////////////////
  // Odd solutions of period _Pi .  (A&S 20.2.4 with p=0)
  // Input:
  //     b = Real characteristic value (?) of Mathieu function, e.g. b = .1
  //     q = Real parameter of Mathieu function, e.g. q = 5.0
  // Output:
  //     answer = [B_0 B_1 B_2 ... B_j] mentioned in A&S 20.2
  // "Determination of Characteristic Values" .  In particular,
  // 20.2.9-10 "Recurrence Relations Among the Coefficients"
  ////////////////////////////////////////////////////////////////
  (Local j B0 B2 B4 m Bm Bm-2 Bm+2 s s0 answer)
  (Setq answer [])
  (Setq j (Max 6 (ISqrt (Ceiling (Abs b)))))
  (Loop (To i 0)
    (Do (Cond ((IsOdd i)
               (Push answer 0))
              ((= i 0)
               // A&S 20.2.1 "where B0 can be taken as zero"
               (Setq B0 0)
               (Push answer B0))
              ((= i 2)
               // Arbitrarily set B2 == 1
               (Setq B2 1)
               (Push answer B2))
              ((= i 4)
               // A&S 20.2.9
               (Setq B4 (/ (* (- b 4) B2) q))
               (Push answer B4))
              (TRUE
               // A&S 20.2.10
               (Setq m (- i 2))
               (Setq Bm (Elt answer m))
               (Setq Bm-2 (Elt answer (- m 2)))
               (Setq Bm+2 (- (/ (* (- b (^ m 2)) Bm) q) Bm-2))
               (When (And (>= i j) (>= (Abs Bm+2) (Abs Bm)))
                 (Break))
               (Push answer Bm+2)))))
  // Normalize to agree with A&S p750 Table 20.2
  (Setq s (Apply '+ (^ answer 2)))
  (Setq s0 (Sqrt (* 1. (^ s -1))))
  (Setq answer (* s0 answer))
  (Return answer)
)

(Def (_FloatMathieuOddPeriodTwoPiSolution b q)
  ////////////////////////////////////////////////////////////////
  // Odd solutions of period 2*_Pi .  (A&S 20.2.4 with p=1)
  // Input:
  //     b = Real characteristic value (?) of Mathieu function, e.g. b = .1
  //     q = Real parameter of Mathieu function, e.g. q = 5.0
  // Output:
  //     answer = [B_0 B_1 B_2 ... B_j] mentioned in A&S 20.2
  // "Determination of Characteristic Values" .  In particular,
  // 20.2.10-11 "Recurrence Relations Among the Coefficients"
  ////////////////////////////////////////////////////////////////
  (Local j B1 B3 m Bm Bm-2 Bm+2 s s0 answer)
  (Setq answer [])
  (Setq j (Max 6 (ISqrt (Ceiling (Abs b)))))
  (Loop (To i 0)
    (Do (Cond ((IsEven i)
               (Push answer 0))
              ((= i 1)
               // Arbitrarily set B1 == 1
               (Setq B1 1)
               (Push answer B1))
              ((= i 3)
               // A&S 20.2.11
               (Setq B3 (/ (* B1 (+ q b -1)) q))
               (Push answer B3))
              (TRUE
               // A&S 20.2.10
               (Setq m (- i 2))
               (Setq Bm (Elt answer m))
               (Setq Bm-2 (Elt answer (- m 2)))
               (Setq Bm+2 (- (/ (* (- b (^ m 2)) Bm) q) Bm-2))
               (When (And (>= i j) (>= (Abs Bm+2) (Abs Bm)))
                 (Break))
               (Push answer Bm+2)))))
  // Normalize to agree with A&S p750 Table 20.2
  (Setq s (Apply '+ (^ answer 2)))
  (Setq s0 (Sqrt (* 1. (^ s -1))))
  (Setq answer (* s0 answer))
  (Return answer)
)