Mathieu functions computational toolbox implemented in Matlab

Mathieu functions computational to olb o x implemen ted in Matlab E. Co jo caru Dep artment of The or etic al Physics, Horia Hulub ei National Institute of Physics and Nucle ar Engine ering, Magur ele-Buchar est P.O.Box MG-6, 077125 R omania ∗ The Mathieu functions are used to solv e analytically some problems in elliptical cylinder co ordinates. A computational to olbox w as implemented in Matlab. Since the notation and normalization for Mathieu functions v ary in the literature, w e ha ve included sufficient material to mak e this presentation self contained. Th us, all form ulas required to get the Mathieu functions are given explicitly . F ollo wing the outlines in this presentation, the Mathieu functions could b e readily implemen ted in other computer programs and used in differen t domains. T ables of n umerical v alues are provided. ∗ Electronic address: eco jo caru@theory .nipne.ro 2 I. INTR ODUCTION Some problems regarding the elliptical cylinders can b e solved by using an analytical approac h lik e that applied to circular cylinders: one separates the v ariables and the exact solution is giv en by expansions inv olving angular and radial Mathieu functions. These functions hav e b een introduced b y Emile Mathieu in 1868 by inv estigating the vibrating mo des in an elliptic mem brane [1]. Details (tables or relations) concerning the Mathieu functions can b e found for example in [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. F or circular cylinders the solutions in volv e readily av ailable trigonometric and Bessel functions, while for elliptical cylinders there are still contro versial and incomplete algorithms for computing the Mathieu functions. Largely applied computer programs provide only few or none routines refering to the Mathieu functions. A computational toolb ox for Mathieu functions w as implemented in Matlab [13]. Since not all p eople are familiarized with the Matlab program, in this presen tation the math- ematics is outlined. T ables of numerical v alues are provided. F ollo wing the outlines in this presen tation, it w ould b e a readily task to implement the Mathieu functions in other computer programs and use them in differen t domains. One reason for the lac k of algorithms for Mathieu functions was probably the complicated and v arious notation existen t in the literature. A main purpose for us w as to simplify as muc h as p ossible the notation. With a simplified and self-contained notation, the use of Mathieu functions should b e as simple as the use of Bessel functions. W e largely follow ed the notations used by Stratton [6] and Stamnes [11, 12], but w e in tro duced further simplifications. Since the notation and normalization for Mathieu functions v ary in the literature, w e hav e included sufficien t material to make this presentation self contained. Thus, all formulas required to get the Mathieu functions are given explicitly . T ables of numerical v alues are provided. Examples of Mathieu functions applied to plane w a ve scattering by elliptical cylinders are giv en in [14, 15]. 3 I I. FUND AMENT ALS A. Elliptical cylinder co ordinates Let consider an ellipse in the plane ( x, y ) defined by equation ( x/x 0 ) 2 + ( y /y 0 ) 2 = 1 with x 0 > y 0 . The semifo cal distance f is given by f 2 = x 2 0 − y 2 0 and the eccentricit y is e = f /x 0 < 1. The elliptic cylindrical co ordinates ( u, v , z ) are defined by relations x = f cosh u cosv , y = f sinh u sin v , z = z (1) with 0 ≤ u < ∞ and 0 ≤ v ≤ 2 π . In terms of ( ξ , η , z ), with ξ = cosh u and η = cos v , the elliptic cylindrical co ordinates are defined by relations x = f ξ η , y = f p ( ξ 2 − 1)(1 − η 2 ) , z = z . (2) The contours of constant u are confo cal ellipses (of semiaxes x 0 = f ξ , y 0 = f p ξ 2 − 1) and those of constan t v are confo cal h yp erb olas. The z axis coincides with the cylinder axis. The scale factors h j , with j = ξ , η , z , are defined lik e as for an y co ordinate transformation [6], h ξ = f p ξ 2 − η 2 p ξ 2 − 1 , h η = f p ξ 2 − η 2 p 1 − η 2 , h z = 1 . (3) B. W a v e equation in elliptic cylindrical co ordinates The scalar wa ve equation ( ∇ 2 + k 2 ) U ( r ) = 0 , where r is the p osition vector, k is the w av e n umber, k = 2 π √ /λ ,  is the p ermittivit y , and λ is the wa velength in v acuum, when expressed in elliptic cylindrical co ordinates b ecomes h 2 f 2 (cosh 2 u − cos 2 v )  ∂ 2 ∂ u 2 + ∂ 2 ∂ v 2  + ∂ 2 ∂ z 2 + k 2 i U ( u, v , z ) = 0 . (4) Using a solution of the form U = Z ( z ) S ( v ) R ( u ) giv es  d 2 d z 2 + k 2 z  Z ( z ) = 0 , (5) h d 2 d v 2 + ( a − 2 q cos 2 v ) i S ( v ) = 0 , (6) h d 2 d u 2 − ( a − 2 q cosh 2 u ) i R ( u ) = 0 , (7) 4 where k z is the w av e v ector comp onen t on z direction, q = k 2 τ f 2 / 4, with k 2 τ = k 2 − k 2 z , and a is separation constan t. Equation (5) has solution Z ( z ) = exp ( ik z z ). Equations (6) and (7) are kno wn as the angular and radial Mathieu equations, resp ectiv ely . I I I. ANGULAR MA THIEU FUNCTIONS In this presen tation, only the p erio dic solutions of p eriod π or 2 π are considered. F or a giv en order n , there are four categories of p eriodic solutions satisfying (6), 1 ev en-even: S ee ( v , q , n ) = ∞ X j =0 A (2 j ) ee ( q , n ) cos(2 j v ) , 2 ev en-o dd: S eo ( v , q , n ) = ∞ X j =0 A (2 j +1) eo ( q , n ) cos[(2 j + 1) v ] , (8) 3 o dd-ev en: S oe ( v , q , n ) = ∞ X j =1 A (2 j ) oe ( q , n ) sin(2 j v ) , 4 o dd-odd: S oo ( v , q , n ) = ∞ X j =0 A (2 j +1) oo ( q , n ) sin[(2 j + 1) v ] . A pm with p, m = e, o are expansion co efficien ts. In the follo wing, the angular Mathieu functions are denoted S pm ( v , q , n ), with p, m = e, o . Instead of t wo angular Mathieu func- tions, ev en S ep and o dd S op , with p = e, o [12], a single angular Mathieu function S pm , with p, m = e, o , is considered refering to all the four categories. F or a giv en v alue of q there exist four infinite sequences of c haracteristic v alues (eigenv alues) a , for either v alue of a corresp onding an infinite sequence (eigenv ector) of expansion co efficien ts. A. Characteristic v alues and co efficien ts By subsituting (8) in (6), the follo wing recurrence relations among the expansion co effi- cien ts result 1 ev en-even: aA (0) ee − q A (2) ee = 0 , ( a − 4) A (2) ee − q [2 A (0) ee + A (4) ee ] = 0 , [ a − (2 j ) 2 ] A (2 j ) ee − q [ A (2 j − 2) ee + A (2 j +2) ee ] = 0 , j = 2 , 3 , 4 · · · (9) 5 2 ev en-o dd: ( a − 1) A (1) eo − q [ A (1) eo + A (3) eo ] = 0 , [ a − (2 j + 1) 2 ] A (2 j +1) eo − q [ A (2 j − 1) eo + A (2 j +3) eo ] = 0 , j = 1 , 2 , 3 · · · (10) 3 o dd-ev en: ( a − 4) A (2) oe − q A (4) oe = 0 , [ a − (2 j ) 2 ] A (2 j ) oe − q [ A (2 j − 2) oe + A (2 j +2) oe ] = 0 , j = 2 , 3 , 4 · · · (11) 4 o dd-odd: ( a − 1) A (1) oo + q [ A (1) oo − A (3) oo ] = 0 , [ a − (2 j + 1) 2 ] A (2 j +1) oo − q [ A (2 j − 1) oo + A (2 j +3) oo ] = 0 , j = 1 , 2 , 3 · · · . (12) The recurrence relations can b e written in matrix form [11], 1 ev en-even:           − a q 0 0 0 0 · · · 2 q 2 2 − a q 0 0 0 · · · 0 q 4 2 − a q 0 0 · · · 0 0 q 6 2 − a q 0 · · · . . . . . . . . . . . . . . . . . . . . .                     A (0) ee A (2) ee A (4) ee A (6) ee . . .           = 0 , (13) 2 ev en-o dd:           1 + q − a q 0 0 0 0 · · · q 3 2 − a q 0 0 0 · · · 0 q 5 2 − a q 0 0 · · · 0 0 q 7 2 − a q 0 · · · . . . . . . . . . . . . . . . . . . . . .                     A (1) eo A (3) eo A (5) eo A (7) eo . . .           = 0 , (14) 3 o dd-ev en:           2 2 − a q 0 0 0 0 · · · q 4 2 − a q 0 0 0 · · · 0 q 6 2 − a q 0 0 · · · 0 0 q 8 2 − a q 0 · · · . . . . . . . . . . . . . . . . . . . . .                     A (2) oe A (4) oe A (6) oe A (8) oe . . .           = 0 , (15) 6 4 o dd-odd:           1 − q − a q 0 0 0 0 · · · q 3 2 − a q 0 0 0 · · · 0 q 5 2 − a q 0 0 · · · 0 0 q 7 2 − a q 0 · · · . . . . . . . . . . . . . . . . . . . . .                     A (1) oo A (3) oo A (5) oo A (7) oo . . .           = 0 . (16) The matrices are real, tridiagonal, and symmetric for all categories, with the exception of the “1 even-ev en” category where the matrix is slightly non-symmetric. The eigenv alue problem is accurately solv ed in Matlab. In other computer programs it could b e necessary to transform the sligh tly non-symmetric matrix in a symmetric one [11]. Both the eigen v alues a and the corresp onding eigen vectors ( A pm , with p, m = e, o ) are determined for either category at an y order n . The order n takes differen t v alues for eac h category of Mathieu functions. F or the purp ose of av oiding an y confusion, a distinction must b e done b et ween the n th order (in the succession of all orders) and the true v alue of that order. Thus, let denote n the order in the succession of all orders, and t the true v alue of order n . The v alues of n and t for the four categories of Mathieu functions are 1 ev en-even: n = 0 , 1 , 2 · · · t = 0 , 2 , 4 · · · , 2 ev en-o dd: n = 0 , 1 , 2 · · · t = 1 , 3 , 5 · · · , 3 o dd-ev en: n = 1 , 2 , 3 · · · t = 2 , 4 , 6 · · · , 4 o dd-odd: n = 0 , 1 , 2 · · · t = 1 , 3 , 5 · · · . Note that, if the notation is self-con tained b y all routines of Mathieu functions, there is no need to determine the sp ecific v alues of n and t for either category of Mathieu functions since it is done automatically . B. Normalization and orthogonality F ollo wing [6, 11], the angular Mathieu functions are normalized b y requiring that S ep (0 , q , n ) = 1 , h d S op ( v , q , n ) d v i v =0 = 1 , p = e, o. (17) 7 These requiremen ts imply that, 1 ev en-even: ∞ X j =0 A (2 j ) ee ( q , n ) = 1 , 2 ev en-o dd: ∞ X j =0 A (2 j +1) eo ( q , n ) = 1 , (18) 3 o dd-ev en: ∞ X j =1 2 j A (2 j ) oe ( q , n ) = 1 , 4 o dd-odd: ∞ X j =0 (2 j + 1) A (2 j +1) oo ( q , n ) = 1 . The orthogonalit y relation for the angular Mathieu functions is Z 2 π 0 S pm ( v , q , n ) S pm 0 ( v , q , n ) d v = N pm δ m m 0 , p, m, m 0 = e, o, (19) where N pm is normalization factor, δ m m 0 equals 1 if m = m 0 and equals 0 otherwise. Then, the follo wing relations for the normalization factor result, 1 ev en-even: N ee ( q , n ) = 2 π [ A (0) ee ( q , n )] 2 + π ∞ X j =1 [ A (2 j ) ee ( q , n )] 2 , 2 ev en-o dd: N eo ( q , n ) = π ∞ X j =0 [ A (2 j +1) eo ( q , n )] 2 , (20) 3 o dd-ev en: N oe ( q , n ) = π ∞ X j =1 [ A (2 j ) oe ( q , n )] 2 , 4 o dd-odd: N oo ( q , n ) = π ∞ X j =0 [ A (2 j +1) oo ( q , n )] 2 . Since different normalization sc hemes ha v e b een adopted in the literature, m uch attention should b e paid when numerical results provided by differen t authors are compared ones against the others. C. Correlation factors Let consider t wo regions of different p ermittivities,  and  0 . The parameter q b eing differen t in the tw o regions, q 6 = q 0 , the c haracteristic v alues and expansion co efficien ts are also differen t. Let S pm and S 0 pm b e the resp ectiv e angular Mathieu functions. The correlation 8 factors C pm ( q , q 0 , n ), with p, m = e, o , b et ween the angular Mathieu functions S pm and S 0 pm are defined by relation C pm ( q , q 0 , n ) = δ m m 0 Z 2 π 0 S pm 0 ( v , q , n ) S 0 pm ( v , q 0 , n ) d v , p, m, m 0 = e, o. (21) Using (8) gives 1 ev en-even: C ee ( q , q 0 , n ) = 2 π A (0) ee ( q , n ) A 0 (0) ee ( q 0 , n ) + π ∞ X j =1 A (2 j ) ee ( q , n ) A 0 (2 j ) ee ( q 0 , n ) , 2 ev en-o dd: C eo ( q , q 0 , n ) = π ∞ X j =0 A (2 j +1) eo ( q , n ) A 0 (2 j +1) eo ( q 0 , n ) , (22) 3 o dd-ev en: C oe ( q , q 0 , n ) = π ∞ X j =1 A (2 j ) oe ( q , n ) A 0 (2 j ) oe ( q 0 , n ) , 4 o dd-odd: C oo ( q , q 0 , n ) = π ∞ X j =0 A (2 j +1) oo ( q , n ) A 0 (2 j +1) oo ( q 0 , n ) . D. Deriv ativ es of angular Mathieu functions The deriv atives of the angular Mathieu functions follow readily from (8), 1 ev en-even: d S ee ( v , q , n ) d v = − ∞ X j =1 2 j A (2 j ) ee ( q , n ) sin(2 j v ) , 2 ev en-o dd: d S eo ( v , q , n ) d v = − ∞ X j =0 (2 j + 1) A (2 j +1) eo ( q , n ) sin[(2 j + 1) v ] , (23) 3 o dd-ev en: d S oe ( v , q , n ) d v = ∞ X j =1 2 j A (2 j ) oe ( q , n ) cos(2 j v ) , 4 o dd-odd: d S oo ( v , q , n ) d v = ∞ X j =0 (2 j + 1) A (2 j +1) oo ( q , n ) cos[(2 j + 1) v ] . IV. RADIAL MA THIEU FUNCTIONS Solutions of (7) can b e obtained from (8) by replacing v by iu . Instead of sin v and cos v , the terms of the series no w inv olve sinh u and cosh u . The conv ergence is lo w unless | u | is small. Better conv ergence of series results by expressing the solutions of (7) in terms of Bessel functions associated with the same expansion co efficien ts that are determined once for 9 b oth the angular and radial Mathieu functions. Either pair of angular and radial Mathieu functions are prop ortional to one another [6], S ep ( iu, q , n ) = √ 2 π g ep ( q , n ) J ep ( u, q , n ) , p = e, o, (24) where J ep are even radial Mathieu functions of the first kind and g ep are joining factors. When u = 0, S ep (0 , q , n ) = 1 , J ep (0 , q , n ) = 1 √ 2 π g ep ( q , n ) , p = e, o. (25) Th us, one obtains, 1 ev en-even: g ee ( q , n ) = ( − 1) r π A (0) ee ( q , n ) S ee ( π / 2 , q , n ) , r = t/ 2 , 2 ev en-o dd: g eo ( q , n ) = − ( − 1) r π √ q A (1) eo ( q , n ) h d S eo ( v , q , n ) d v i v = π/ 2 , r = ( t − 1) / 2 . (26) Similarly [6], − iS op ( iu, q , n ) = √ 2 π g op ( q , n ) J op ( u, q , n ) , p = e, o. (27) When u = 0, J op (0 , q , n ) = 0 , h d J op ( u, q , n ) d u i u =0 = 1 √ 2 π g op ( q , n ) , p = e, o. (28) Th us, one obtains, 3 o dd-ev en: g oe ( q , n ) = ( − 1) r π q A (2) oe ( q , n ) h d S oe ( v , q , n ) d v i v = π/ 2 , r = t/ 2 , 4 o dd-odd: g oo ( q , n ) = ( − 1) r π √ q A (1) oo ( q , n ) S oo ( π / 2 , q , n ) , r = ( t − 1) / 2 . (29) Remem b er that t is the true v alue of order n . A. Radial Mathieu functions of the first kind Since rapidly conv erging series are those expressed in terms of pro ducts of Bessel functions [10, 11], in the follo wing relations refer only to them. Similarly to the angular Mathieu functions, one ma y distinct four categories of radial Mathieu functions of the first kind 10 whic h are denoted J pm ( u, q , n ), with p, m = e, o , 1 ev en-even: J ee ( u, q , n ) = r π 2 ( − 1) r A (0) ee ( q , n ) ∞ X j =0 ( − 1) j A (2 j ) ee ( q , n ) J j ( v 1 ) J j ( v 2 ) , r = t/ 2 , 2 ev en-o dd: J eo ( u, q , n ) = r π 2 ( − 1) r A (1) eo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) eo ( q , n )[ J j ( v 1 ) J j +1 ( v 2 ) + J j ( v 2 ) J j +1 ( v 1 )] , r = ( t − 1) / 2 , (30) 3 o dd-ev en: J oe ( u, q , n ) = r π 2 ( − 1) r A (2) oe ( q , n ) ∞ X j =1 ( − 1) j A (2 j ) oe ( q , n )[ J j − 1 ( v 1 ) J j +1 ( v 2 ) − J j − 1 ( v 2 ) J j +1 ( v 1 )] , r = t/ 2 , 4 o dd-odd: J oo ( u, q , n ) = r π 2 ( − 1) r A (1) oo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) oo ( q , n )[ J j ( v 1 ) J j +1 ( v 2 ) − J j ( v 2 ) J j +1 ( v 1 )] , r = ( t − 1) / 2 , where v 1 = √ q exp ( − u ) and v 2 = √ q exp ( u ). The deriv atives of the radial Mathieu func- tions of the first kind are 1 ev en-even: r = t/ 2 , d J ee ( u, q , n ) d u = r π 2 ( − 1) r A (0) ee ( q , n ) ∞ X j =0 ( − 1) j A (2 j ) ee ( q , n )[ v 1 J j +1 ( v 1 ) J j ( v 2 ) − v 2 J j ( v 1 ) J j +1 ( v 2 )] , 2 ev en-o dd: r = ( t − 1) / 2 , d J eo ( u, q , n ) d u = r π 2 ( − 1) r A (1) eo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) eo ( q , n ) n ( v 2 − v 1 )[ J j ( v 1 ) J j ( v 2 ) − J j +1 ( v 1 ) J j +1 ( v 2 )] + (2 j + 1)[ J j +1 ( v 1 ) J j ( v 2 ) − J j ( v 1 ) J j +1 ( v 2 )] o , 3 o dd-ev en: r = t/ 2 , d J oe ( u, q , n ) d u = r π 2 ( − 1) r A (2) oe ( q , n ) ∞ X j =0 ( − 1) j +1 A (2 j +2) oe ( q , n )(4 j + 4) n J j ( v 1 ) J j ( v 2 ) (31) + cosh 2 uJ j +1 ( v 1 ) J j +1 ( v 2 ) − ( j + 1)[ 1 v 1 J j +1 ( v 1 ) J j ( v 2 ) + 1 v 2 J j ( v 1 ) J j +1 ( v 2 )] o , 11 4 o dd-odd: r = ( t − 1) / 2 , d J oo ( u, q , n ) d u = r π 2 ( − 1) r A (1) oo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) oo ( q , n ) n ( v 1 + v 2 )[ J j ( v 1 ) J j ( v 2 ) + J j +1 ( v 1 ) J j +1 ( v 2 )] − (2 j + 1)[ J j +1 ( v 1 ) J j ( v 2 ) + J j ( v 1 ) J j +1 ( v 2 )] o . B. Radial Mathieu functions of the second kind A second indep enden t solution of (7) is obtained by replacing the Bessel functions of the first kind J n ( v 2 ) in (30) by the Bessel functions of the second kind Y n ( v 2 ) [10, 11]. This solution is denoted Y pm ( u, q , n ), with p, m = e, o . 1 ev en-even: Y ee ( u, q , n ) = r π 2 ( − 1) r A (0) ee ( q , n ) ∞ X j =0 ( − 1) j A (2 j ) ee ( q , n ) J j ( v 1 ) Y j ( v 2 ) , r = t/ 2 , 2 ev en-o dd: Y eo ( u, q , n ) = r π 2 ( − 1) r A (1) eo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) eo ( q , n )[ J j ( v 1 ) Y j +1 ( v 2 ) + Y j ( v 2 ) J j +1 ( v 1 )] , r = ( t − 1) / 2 , (32) 3 o dd-ev en: Y oe ( u, q , n ) = r π 2 ( − 1) r A (2) oe ( q , n ) ∞ X j =1 ( − 1) j A (2 j ) oe ( q , n )[ J j − 1 ( v 1 ) Y j +1 ( v 2 ) − Y j − 1 ( v 2 ) J j +1 ( v 1 )] , r = t/ 2 , 4 o dd-odd: Y oo ( u, q , n ) = r π 2 ( − 1) r A (1) oo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) oo ( q , n )[ J j ( v 1 ) Y j +1 ( v 2 ) − Y j ( v 2 ) J j +1 ( v 1 )] , r = ( t − 1) / 2 , The deriv atives of the radial Mathieu functions of the second kind are 1 ev en-even: r = t/ 2 , d Y ee ( u, q , n ) d u = r π 2 ( − 1) r A (0) ee ( q , n ) ∞ X j =0 ( − 1) j A (2 j ) ee ( q , n )[ v 1 J j +1 ( v 1 ) Y j ( v 2 ) − v 2 J j ( v 1 ) Y j +1 ( v 2 )] , 12 2 ev en-o dd: r = ( t − 1) / 2 , d Y eo ( u, q , n ) d u = r π 2 ( − 1) r A (1) eo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) eo ( q , n ) n ( v 2 − v 1 )[ J j ( v 1 ) Y j ( v 2 ) − J j +1 ( v 1 ) Y j +1 ( v 2 )] + (2 j + 1)[ J j +1 ( v 1 ) Y j ( v 2 ) − J j ( v 1 ) Y j +1 ( v 2 )] o , 3 o dd-ev en: r = t/ 2 , d Y oe ( u, q , n ) d u = r π 2 ( − 1) r A (2) oe ( q , n ) ∞ X j =0 ( − 1) j +1 A (2 j +2) oe ( q , n )(4 j + 4) n J j ( v 1 ) Y j ( v 2 ) (33) + cosh 2 uJ j +1 ( v 1 ) Y j +1 ( v 2 ) − ( j + 1)[ 1 v 1 J j +1 ( v 1 ) Y j ( v 2 ) + 1 v 2 J j ( v 1 ) Y j +1 ( v 2 )] o , 4 o dd-odd: r = ( t − 1) / 2 , d Y oo ( u, q , n ) d u = r π 2 ( − 1) r A (1) oo ( q , n ) ∞ X j =0 ( − 1) j A (2 j +1) oo ( q , n ) n ( v 1 + v 2 )[ J j ( v 1 ) Y j ( v 2 ) + J j +1 ( v 1 ) Y j +1 ( v 2 )] − (2 j + 1)[ J j +1 ( v 1 ) Y j ( v 2 ) + J j ( v 1 ) Y j +1 ( v 2 )] o . C. Radial Mathieu functions of the third and the fourth kinds Radial Mathieu functions of the third kind, analogous to the Hank el functions of the first kind are defined as follows [6, 11] H pm 1 ( u, q , n ) = J pm ( u, q , n ) + iY pm ( u, q , n ) , p, m = e, o. (34) Similarly , radial Mathieu functions of the fourth kind, analogous to the Hank el functions of the second kind are defined as follo ws [6, 11] H pm 2 ( u, q , n ) = J pm ( u, q , n ) − iY pm ( u, q , n ) , p, m = e, o. (35) V. IMPLEMENT A TION OF MA THIEU FUNCTIONS IN MA TLAB F ollo wing the notation of the four categories of angular Mathieu functions, the implemen tation in Matlab or in any other computer program is readily done by in tro ducing a function co de K F . The first step in an y algorithm of Mathieu function computation is to find the c haracteristic v alues (eigenv alues) and the expansion co efficien ts (eigen vectors). In [13], this is done by routine “eig Spm” whic h has q as input parameter 13 (see T able I). Besides q , the function co de K F should b e sp ecified. Thus, if K F = 1, the routine “eig Spm” solves the eigen v alue problem for category “1 even-ev en” of Mathieu functions, if K F = 2 for category “2 ev en-o dd”, and so on. The num b er of expansion co efficien ts is the same, it is set equal to 25, for all categories of Mathieu functions. Concerning the outputs of routine “eig Spm”, v a is a line v ector representing the c haracteristic v alues a for all the 25 orders; mc is 25 × 25 matrix, where the columns represen t the eigenv ectors (that is, the expansion co efficien ts) for all orders; v t is a column v ector sp ecifying the true v alue t for all orders. Note that the eigenv ectors in mc were pro cessed to ob ey equation (18). F or the purp ose to sav e the time of computation, all the other routines hav e mc as input (see T able I), the routine “eig Spm” b eing called once, at the b eginning of the computation, for any v alues of co ordinates u and v that in tervene in that computation. Since in many cases the conv ergence is assured b y the first sev eral orders, all the other routines ha ve nmax ≤ 25 as input. It means that those routines tak e in to account only the first nmax orders, but for either order the length of the corresp onding eigenv ector is the same, equal to 25. The routine “extract one v alue” can b e used to get a single v alue, and the routine “extract one column” to get a single eigen vector, corresp onding to the order t . The deriv atives of S pm , with p, m = e, o , are computed b y routine “dSpm”. F or b oth “Spm” and “dSpm”, v is expressed in radians, with v alues in interv al (0 , 2 π ). The normalization, correlation, and joining factors are computed b y routines “Npm”, “Cpm”, and “gpm”, resp ectively . The four kinds of radial Mathieu functions, J pm , Y pm , H pm 1 , and H pm 2 , with p, m = e, o , are computed by routines “Jpm”,“Ypm”,“Hpm1”, and “Hpm2”, resp ectiv ely , and their deriv atives with resp ect to u b y routines “dJpm”,“dYpm”,“dHpm1”, and “dHpm2”, resp ectiv ely . Numerical v alues of the separation constan t a , of the angular Mathieu functions S pm and their deriv ativ es S 0 pm , with p, m = e, o , where the prime denotes differentia tion with resp ect to v , are given in T ables I I – IV. They can b e compared with data in [2]. With the purp ose to facilitate the comparison, since in [2] the normalization N pm = π is applied, the data of S pm and S 0 pm in T ables I I – IV are multiplied by p π / N pm . Concerning the radial Mathieu functions, n umerical v alues of S ep ( iu, q , n ) and − iS op ( iu, q , n ) are given for u = 0 . 5 in T ables V and VI. They are multiplied by p π / N pm and compared with data in [9]. Note that S ep is correlated to the radial Mathieu function of the first kind J ep b y Eq. (24), whereas S op is correlated to J op b y Eq. (27). W e found 14 that, for parameters in [9], the v alues of S ep ( iu, q , n ) and − iS op ( iu, q , n ) calculated with Eqs. (24) and (27) differ from those obtained with Eq. (8) b y less than 7 . 5 × 10 − 12 . [1] E. Mathieu “Le mouvemen t vibratoire d’une mem brane de forme elliptique,” Jour. de Math. Pures at Appliquees (Jour. de Liouville) 13 , 137–203 (1868). [2] M. Abramo witz and I. Stegun Handb o ok of Mathematic al F unctions (New Y ork, 1964). [3] I. S. Gradshteyn and I. M. Ryzhik T ables of Inte gr als, Series, and Pr o ducts (Academic Press, San Diego, 1994). [4] E. L. Ince Or dinary Differ ential Equations (New Y ork, 1967). [5] N. W. McLachlan The ory and Applic ation of Mathieu F unctions (Oxford Press, 1951). [6] J. A. Stratton Ele ctr omagnetic The ory (Mc-Gra w Hill New Y ork, 1941). [7] J. A. Stratton and P . M. Morse El liptic Cylinder and Spher oidal Wave F unctions Including T ables of Sep ar ation Constants and Co efficients (John Wiley & Sons, New Y ork, 1941). [8] E. T. Whittaker and G. N. W atson A Course of Mo dern A nalysis (Cam bridge Universit y Press, Cambridge, 1950). [9] E. T. Kirkpatrick, “T ables of v alues of the mo dified Mathieu functions,” Mathematics of Computation 14 118–129 (1960). [10] J. C. Gutierrez-V ega, F ormal analysis of the propagation of in v arian t optical fields in elliptic co ordinates, Ph. D. Thesis, INA OE, Mexico, 2000. [11] J. J. Stamnes and B. Sp jelk avik “New metho d for computing eigenfunctions (Mathieu functions) for scattering b y elliptical cylinders,” Pure Appl. Opt. 4 251–62 (1995). [12] J. J. Stamnes “Exact t wo-dimensional scattering b y p erfectly reflecting elliptical cylinders, strips and slits,” Pure Appl. Opt. 4 841–55 (1995). [13] E. Co jo caru, Matlab free a v ailable computer co de Mathieu F unctions T o olb o x v. 1.0; also free av ailable b y request at eco jo caru@theory .nipne.ro or co jo caru.e@gmail.com [14] E. Co jo caru, “Mathieu functions approac h to bidimensional scattering b y dielectric elliptical cylinders,” [15] E. Co jo caru, “Elliptical cylindrical in visibility cloak, a semianalytical approac h using Mathieu functions,” 15 T ABLE I: Routines comprised in the to olb o x [13]. Name of routine Routine call What the routine computes eig Spm [ v a, mc, v t ]=eig Spm( K F , q ) V ector of c haracteristic v alues v a , matrix of co efficien ts mc , and v ector of orders v t , at giv en function co de K F and elliptical parameter q ≥ 0. Spm y =Spm( K F , v , mc, nmax ) Angular Mathieu functions S pm , [Eq. (8)]. dSpm y =dSpm( K F , v , mc, nmax ) Deriv atives with resp ect to v of S pm , [Eq. (23)]. Npm y =Npm( K F , mc, nmax ) Normalizing factors of angular Mathieu functions S pm , [Eqs. (19) and (20)]. Cpm y =Cpm( K F , mc, mc 0 , nmax ) Correlation factors of S pm and S 0 pm , having matrices of co efficien ts mc and mc 0 , [Eqs. (21) and (22)]. Jpm y =Jpm( K F , u, q , mc, nmax ) Radial Mathieu functions of the first kind J pm , [Eq. (30)]. dJpm y =dJpm( K F , u, q , mc, nmax ) Deriv atives with resp ect to u of J pm , [Eq. (31)]. gpm y =gpm( K F , q , mc, nmax ) Joining factors for pairs of angular, S pm and radial, J pm Mathieu functions, [Eqs. (24)–(29)]. Ypm y =Ypm( K F , u, q , mc, nmax ) Radial Mathieu functions of the second kind Y pm , [Eq. (32)]. dYpm y =dYpm( K F , u, q , mc, nmax ) Deriv atives with resp ect to u of Y pm , [Eq. (33)]. Hpm1 y =Hpm1( K F , u, q , mc, nmax ) Radial Mathieu functions of the third kind H pm 1 , [Eq. (34)]. dHpm1 y =dHpm1( K F , u, q , mc, namax ) Deriv atives with resp ect to u of H pm 1 . Hpm2 y =Hpm2( K F , u, q , mc, nmax ) Radial Mathieu functions of the fourth kind H pm 2 , [Eq. (35)]. dHpm2 y =dHpm2( K F , u, q , mc, namax ) Deriv atives with resp ect to u of H pm 2 . extract one column y =extract one column( K F , t, mc ) Extracts one column from mc at given t . extract one v alue y =extract one v alue( K F , t, v ec ) Extracts one v alue from v ec at given t . 16 T ABLE I I: V alues of S ee m ultiplied b y γ ee = p π / N ee to b e compared with data in [2] t q a γ ee S ee (0 , q , n ) γ ee S ee ( π / 2 , q , n ) 0 0 0 0.7071067811865 0.7071067811865 5 -5.8000460208515 0.0448001816519 1.3348486746980 10 -13.9369799566589 0.0076265175709 1.4686604707129 15 -22.5130377608640 0.0019325083152 1.5501081466866 20 -31.3133900703364 0.0006037438292 1.6098908573959 25 -40.2567795465667 0.0002158630184 1.6575102983235 2 0 4.0000000000000 1.0000000000000 -1.0000000000000 5 7.4491097395292 0.7352943084007 -0.7244881519677 10 7.7173698497796 0.2458883492913 -0.9267592641263 15 5.0779831975435 0.0787928278464 -1.0199662260303 20 1.1542828852468 0.0286489431471 -1.0752932287797 25 -3.5221647271583 0.0115128663309 -1.1162789532953 10 0 100.0000000000000 1.0000000000000 -1.0000000000000 5 100.1263692161636 1.0259950270894 -0.9753474872360 10 100.5067700246816 1.0538159921009 -0.9516453181790 15 101.1452034473016 1.0841063118392 -0.9285480638845 20 102.0489160244372 1.1177886312594 -0.9057107845941 25 103.2302048044949 1.1562399186322 -0.8826919105637 17 T ABLE I I I: V alues of S eo and S 0 eo m ultiplied b y γ eo = p π / N eo to b e compared with data in [2] t q a γ eo S eo (0 , q , n ) γ eo S 0 eo ( π / 2 , q , n ) 1 0 1.0000000000000 1.0000000000000 -1.0000000000000 5 1.8581875415478 0.2565428793224 -3.4690420034057 10 -2.3991424000363 0.0535987477472 -4.8504383044964 15 -8.1011051316418 0.0150400664538 -5.7642064390510 20 -14.4913014251748 0.0050518137647 -6.4905657825800 25 -21.3148996906657 0.0019110515067 -7.1067412352901 5 0 25.0000000000000 1.0000000000000 -5.0000000000000 5 25.5499717499816 1.1248072506385 -5.3924861549882 10 27.7037687339393 1.2580199413083 -5.3212765411609 15 31.9578212521729 1.1934322304131 -5.1191498884064 20 36.6449897341328 0.9365755314226 -5.7786752500644 25 40.0501909858077 0.6106943100507 -7.0598842916553 15 0 225.0000000000000 1.0000000000000 15.0000000000000 5 225.0558124767096 1.0112937325296 15.1636574720602 10 225.2233569749644 1.0228782824382 15.3198803056623 15 225.5029562446541 1.0347936522369 15.4687435032830 20 225.8951534162079 1.0470843441629 15.6102785232380 25 226.4007200447481 1.0598004418139 15.7444725050679 18 T ABLE IV: V alues of S op and S 0 op m ultiplied b y p π / N op , p = e, o (see [2]) t q a p π / N oe S 0 oe (0 , q , n ) p π / N oe S 0 oe ( π / 2 , q , n ) 2 0 4.0000000000000 2.0000000000000 -2.0000000000000 5 2.0994604454867 0.7331661960372 -3.6405178524082 10 -2.3821582359570 0.2488228403985 -4.8634220691653 15 -8.0993467988959 0.0918197143696 -5.7655737717278 20 -14.4910632559807 0.0370277776852 -6.4907522240373 25 -21.3148606222498 0.0160562170491 -7.1067719073739 10 0 100.0000000000000 10.0000000000000 -10.0000000000000 5 100.1263692156019 9.7341731518695 -10.2396462566908 10 100.5067694628784 9.4404054347686 -10.4539475316485 15 101.1451722929092 9.1157513395126 -10.6428998776563 20 102.0483928609361 8.7555450801360 -10.8057241781325 25 103.2256800423735 8.3526783655914 -10.9413538308191 t q a p π / N oo S 0 oo (0 , q , n ) p π / N oo S oo ( π / 2 , q , n ) 1 0 1.0000000000000 1.0000000000000 1.0000000000000 5 -5.7900805986378 0.1746754006198 1.3374338870223 10 -13.9365524792501 0.0440225659111 1.4687556641029 15 -22.5130034974235 0.0139251347875 1.5501150743576 20 -31.3133861669129 0.0050778849001 1.6098915926038 25 -40.2567789846842 0.0020443593656 1.6575103983745 5 0 25.0000000000000 5.0000000000000 1.0000000000000 5 25.5108160463032 4.3395700104946 0.9060779302024 10 26.7664263604801 3.4072267604013 0.8460384335355 15 27.9678805967175 2.4116664728002 0.8379493400125 20 28.4682213251027 1.5688968684857 0.8635431218534 25 28.0627658994543 0.9640716219024 0.8992683245108 15 0 225.0000000000000 15.0000000000000 -1.0000000000000 5 225.0558124767096 14.8287889732852 -0.9889607027406 10 225.2233569749643 14.6498600449581 -0.9781423471832 15 225.5029562446537 14.4630006940372 -0.9675137031855 20 225.8951534161767 14.2679460909928 -0.9570452540613 25 226.4007200438825 14.0643732956172 -0.9467086958781 19 T ABLE V: V alues of S ep ( iu, q , n ) for u = 0 . 5 m ultiplied by p π / N ep , where p = e, o , compared with data in [9] t q V alues at p = e Data in [9] t q V alues at p = o Data in [9] 0 5 -0.019325304910071 -0.01932 1 5 0.021440743185527 0.02144 10 -0.007055239716193 -0.00705 10 -0.038634237458525 -0.03863 20 -0.000169411415735 -0.00016 20 -0.003373888309642 -0.00337 2 5 0.446937465741068 0.44693 3 5 1.205528267066838 1.2055 10 -0.063855921612085 -0.06385 10 0.235940782144547 0.23594 20 -0.024916657795101 -0.02491 20 -0.097385461808731 -0.09738 4 5 2.234088244534832 2.2341 5 5 3.864089377116713 3.8641 10 1.039103163573830 1.0391 10 2.285610444240526 2.2856 20 -0.143991090269732 -0.14399 20 0.274270780278172 0.27427 T ABLE VI: V alues of − iS op ( iu, q , n ) for u = 0 . 5 multiplied by p π / N op , with p = e, o , compared with data in [9] t q V alues at p = e Data in [9] t q V alues at p = o Data in [9] 2 5 0.238342768735937 0.23834 1 5 0.036613617783886 0.03661 10 0.028675814044625 0.02867 10 0.000750806874015 0.00075 20 -0.003176296415956 -0.00317 20 -0.000538258353937 -0.00053 4 5 1.883560277440876 1.8836 3 5 0.806555153528872 0.80655 10 0.769679129538722 0.76968 10 0.204495885546638 0.20449 20 0.040515136278697 0.04051 20 -0.005279473480675 -0.00527 6 5 6.6066602369876 6.6067 5 5 3.667530204538722 3.6675 10 4.1161420952367 4.1161 10 1.972361938552091 1.9724 20 1.1805904286267 1.1806 20 0.320398855944192 0.32040