Modelica.Blocks.Continuous.Internal.Filter.Utilities

Utility functions for filter computations

Information

Extends from Modelica.Icons.Package (Icon for standard packages).

Package Content

Name Description

BesselBaseCoefficients Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB

toHighestPowerOne Transform filter to form with highest power of s equal 1

normalizationFactor Compute correction factor of low pass filter such that amplitude at cut-off frequency is -3db (=10^(-3/20) = 0.70794...)

bandPassAlpha Return alpha for band pass

Name	Description
BesselBaseCoefficients	Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB
toHighestPowerOne	Transform filter to form with highest power of s equal 1
normalizationFactor	Compute correction factor of low pass filter such that amplitude at cut-off frequency is -3db (=10^(-3/20) = 0.70794...)
bandPassAlpha	Return alpha for band pass

Modelica.Blocks.Continuous.Internal.Filter.Utilities.BesselBaseCoefficients

Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB

Information

The transfer function H(p) of a n 'th order Bessel filter is given by

         Bn(0)
 H(p) = -------
         Bn(p)

with the denominator polynomial

          n             n  (2n - k)!       p^k
 Bn(p) = sum c_k*p^k = sum ----------- * -------   (1)
         k=0           k=0 (n - k)!k!    2^(n-k)

and the numerator

               (2n)!     1
Bn(0) = c_0 = ------- * ---- .                     (2)
                n!      2^n

Although the coefficients c_k are integer numbers, it is not advisable to use the polynomials in an unfactorized form because the coefficients are fast growing with order n (c_0 is approximately 0.3e24 and 0.8e59 for order n=20 and order n=40 respectively).
Therefore, the polynomial Bn(p) is factorized to first and second order polynomials with real coefficients corresponding to zeros and poles representation that is used in this library.

The function returns the coefficients which resulted from factorization of the normalized transfer function

H'(p') = H(p),  p' = p/w0

as well as

alpha = 1/w0

the reciprocal of the cut of frequency w0 where the gain of the transfer function is decreased 3dB.

Both, coefficients and cut off frequency were calculated symbolically and were eventually evaluated with high precision calculation. The results were stored in this function as real numbers.

Calculation of normalized Bessel filter coefficients

Equation

   abs(H(j*w0)) = abs(Bn(0)/Bn(j*w0)) = 10^(-3/20)

which must be fulfilled for cut off frequency w = w0 leads to

   [Re(Bn(j*w0))]^2 + [Im(Bn(j*w0))]^2 - (Bn(0)^2)*10^(3/10) = 0

which has exactly one real solution w0 for each order n. This solutions of w0 are calculated symbolically first and evaluated by using high precise values of the coefficients c_k calculated by following (1) and (2).
With w0, the coefficients of the factorized polynomial can be computed by calculating the zeros of the denominator polynomial

         n
 Bn(p) = sum w0^k*c_k*(p/w0)^k
         k=0

of the normalized transfer function H'(p'). There exist n/2 of conjugate complex pairs of zeros (beta +-j*gamma) if n is even and one additional real zero (alpha) if n is odd. Finally, the coefficients a, b1_k, b2_k of the polynomials

 a*p + 1,  n is odd

and

 b2_k*p^2 + b1_k*p + 1,   k = 1,... div(n,2)

results from

 a = -1/alpha

and

 b2_k = 1/(beta_k^2 + gamma_k^2) b1_k = -2*beta_k/(beta_k^2 + gamma_k^2)

Inputs

Type Name Default Description

Integer order Order of filter in the range 1..41

Type	Name	Default	Description
Integer	order		Order of filter in the range 1..41

Outputs

Type Name Description

Real c1[mod(order, 2)] [p] coefficients of Bessel denominator polynomials (a*p + 1)

Real c2[integer(order/2), 2] [p^2, p] coefficients of Bessel denominator polynomials (b2*p^2 + b1*p + 1)

Real alpha Normalization factor

Type	Name	Description
Real	c1[mod(order, 2)]	[p] coefficients of Bessel denominator polynomials (a*p + 1)
Real	c2[integer(order/2), 2]	[p^2, p] coefficients of Bessel denominator polynomials (b2p^2 + b1p + 1)
Real	alpha	Normalization factor

Modelica definition

function BesselBaseCoefficients 
  "Return coefficients of normalized low pass Bessel filter (= gain at cut-off frequency 1 rad/s is decreased 3dB"

  import Modelica.Utilities.Streams;
  input Integer order "Order of filter in the range 1..41";
  output Real c1[mod(order, 2)] 
    "[p] coefficients of Bessel denominator polynomials (a*p + 1)";
  output Real c2[integer(order/2),2] 
    "[p^2, p] coefficients of Bessel denominator polynomials (b2*p^2 + b1*p + 1)";
  output Real alpha "Normalization factor";
algorithm 
  if order == 1 then
    alpha := 1.002377293007601;
    c1[1] := 0.9976283451109835;
  elseif order == 2 then
    alpha := 0.7356641785819585;
    c2[1, 1] := 0.6159132201783791;
    c2[1, 2] := 1.359315879600889;
  elseif order == 3 then
    alpha := 0.5704770156982642;
    c1[1] := 0.7548574865985343;
    c2[1, 1] := 0.4756958028827457;
    c2[1, 2] := 0.9980615136104388;
  elseif order == 4 then
    alpha := 0.4737978580281427;
    c2[1, 1] := 0.4873729247240677;
    c2[1, 2] := 1.337564170455762;
    c2[2, 1] := 0.3877724315741958;
    c2[2, 2] := 0.7730405590839861;
  elseif order == 5 then
    alpha := 0.4126226974763408;
    c1[1] := 0.6645723262620757;
    c2[1, 1] := 0.4115231900614016;
    c2[1, 2] := 1.138349926728708;
    c2[2, 1] := 0.3234938702877912;
    c2[2, 2] := 0.6205992985771313;
  elseif order == 6 then
    alpha := 0.3705098000736233;
    c2[1, 1] := 0.3874508649098960;
    c2[1, 2] := 1.219740879520741;
    c2[2, 1] := 0.3493298843155746;
    c2[2, 2] := 0.9670265529381365;
    c2[3, 1] := 0.2747419229514599;
    c2[3, 2] := 0.5122165075105700;
  elseif order == 7 then
    alpha := 0.3393452623586350;
    c1[1] := 0.5927147125821412;
    c2[1, 1] := 0.3383379423919174;
    c2[1, 2] := 1.092630816438030;
    c2[2, 1] := 0.3001025788696046;
    c2[2, 2] := 0.8289928256598656;
    c2[3, 1] := 0.2372867471539579;
    c2[3, 2] := 0.4325128641920154;
  elseif order == 8 then
    alpha := 0.3150267393795002;
    c2[1, 1] := 0.3151115975207653;
    c2[1, 2] := 1.109403015460190;
    c2[2, 1] := 0.2969344839572762;
    c2[2, 2] := 0.9737455812222699;
    c2[3, 1] := 0.2612545921889538;
    c2[3, 2] := 0.7190394712068573;
    c2[4, 1] := 0.2080523342974281;
    c2[4, 2] := 0.3721456473047434;
  elseif order == 9 then
    alpha := 0.2953310177184124;
    c1[1] := 0.5377196679501422;
    c2[1, 1] := 0.2824689124281034;
    c2[1, 2] := 1.022646191567475;
    c2[2, 1] := 0.2626824161383468;
    c2[2, 2] := 0.8695626454762596;
    c2[3, 1] := 0.2302781917677917;
    c2[3, 2] := 0.6309047553448520;
    c2[4, 1] := 0.1847991729757028;
    c2[4, 2] := 0.3251978031287202;
  elseif order == 10 then
    alpha := 0.2789426890619463;
    c2[1, 1] := 0.2640769908255582;
    c2[1, 2] := 1.019788132875305;
    c2[2, 1] := 0.2540802639216947;
    c2[2, 2] := 0.9377020417760623;
    c2[3, 1] := 0.2343577229427963;
    c2[3, 2] := 0.7802229808216112;
    c2[4, 1] := 0.2052193139338624;
    c2[4, 2] := 0.5594176813008133;
    c2[5, 1] := 0.1659546953748916;
    c2[5, 2] := 0.2878349616233292;
  elseif order == 11 then
    alpha := 0.2650227766037203;
    c1[1] := 0.4950265498954191;
    c2[1, 1] := 0.2411858478546218;
    c2[1, 2] := 0.9567800996387417;
    c2[2, 1] := 0.2296849355380925;
    c2[2, 2] := 0.8592523717113126;
    c2[3, 1] := 0.2107851705677406;
    c2[3, 2] := 0.7040216048898129;
    c2[4, 1] := 0.1846461385164021;
    c2[4, 2] := 0.5006729207276717;
    c2[5, 1] := 0.1504217970817433;
    c2[5, 2] := 0.2575070491320295;
  elseif order == 12 then
    alpha := 0.2530051198547209;
    c2[1, 1] := 0.2268294941204543;
    c2[1, 2] := 0.9473116570034053;
    c2[2, 1] := 0.2207657387793729;
    c2[2, 2] := 0.8933728946287606;
    c2[3, 1] := 0.2087600700376653;
    c2[3, 2] := 0.7886236252756229;
    c2[4, 1] := 0.1909959101492760;
    c2[4, 2] := 0.6389263649257017;
    c2[5, 1] := 0.1675208146048472;
    c2[5, 2] := 0.4517847275162215;
    c2[6, 1] := 0.1374257286372761;
    c2[6, 2] := 0.2324699157474680;
  elseif order == 13 then
    alpha := 0.2424910397561007;
    c1[1] := 0.4608848369928040;
    c2[1, 1] := 0.2099813050274780;
    c2[1, 2] := 0.8992478823790660;
    c2[2, 1] := 0.2027250423101359;
    c2[2, 2] := 0.8328117484224146;
    c2[3, 1] := 0.1907635894058731;
    c2[3, 2] := 0.7257379204691213;
    c2[4, 1] := 0.1742280397887686;
    c2[4, 2] := 0.5830640944868014;
    c2[5, 1] := 0.1530858190490478;
    c2[5, 2] := 0.4106192089751885;
    c2[6, 1] := 0.1264090712880446;
    c2[6, 2] := 0.2114980230156001;
  elseif order == 14 then
    alpha := 0.2331902368695848;
    c2[1, 1] := 0.1986162311411235;
    c2[1, 2] := 0.8876961808055535;
    c2[2, 1] := 0.1946683341271615;
    c2[2, 2] := 0.8500754229171967;
    c2[3, 1] := 0.1868331332895056;
    c2[3, 2] := 0.7764629313723603;
    c2[4, 1] := 0.1752118757862992;
    c2[4, 2] := 0.6699720402924552;
    c2[5, 1] := 0.1598906457908402;
    c2[5, 2] := 0.5348446712848934;
    c2[6, 1] := 0.1407810153019944;
    c2[6, 2] := 0.3755841316563539;
    c2[7, 1] := 0.1169627966707339;
    c2[7, 2] := 0.1937088226304455;
  elseif order == 15 then
    alpha := 0.2248854870552422;
    c1[1] := 0.4328492272335646;
    c2[1, 1] := 0.1857292591004588;
    c2[1, 2] := 0.8496337061962563;
    c2[2, 1] := 0.1808644178280136;
    c2[2, 2] := 0.8020517898136011;
    c2[3, 1] := 0.1728264404199081;
    c2[3, 2] := 0.7247449729331105;
    c2[4, 1] := 0.1616970125901954;
    c2[4, 2] := 0.6205369315943097;
    c2[5, 1] := 0.1475257264578426;
    c2[5, 2] := 0.4929612162355906;
    c2[6, 1] := 0.1301861023357119;
    c2[6, 2] := 0.3454770708040735;
    c2[7, 1] := 0.1087810777120188;
    c2[7, 2] := 0.1784526655428406;
  elseif order == 16 then
    alpha := 0.2174105053474761;
    c2[1, 1] := 0.1765637967473151;
    c2[1, 2] := 0.8377453068635511;
    c2[2, 1] := 0.1738525357503125;
    c2[2, 2] := 0.8102988957433199;
    c2[3, 1] := 0.1684627004613343;
    c2[3, 2] := 0.7563265923413258;
    c2[4, 1] := 0.1604519074815815;
    c2[4, 2] := 0.6776082294687619;
    c2[5, 1] := 0.1498828607802206;
    c2[5, 2] := 0.5766417034027680;
    c2[6, 1] := 0.1367764717792823;
    c2[6, 2] := 0.4563528264410489;
    c2[7, 1] := 0.1209810465419295;
    c2[7, 2] := 0.3193782657322374;
    c2[8, 1] := 0.1016312648007554;
    c2[8, 2] := 0.1652419227369036;
  elseif order == 17 then
    alpha := 0.2106355148193306;
    c1[1] := 0.4093223608497299;
    c2[1, 1] := 0.1664014345826274;
    c2[1, 2] := 0.8067173752345952;
    c2[2, 1] := 0.1629839591538256;
    c2[2, 2] := 0.7712924931447541;
    c2[3, 1] := 0.1573277802512491;
    c2[3, 2] := 0.7134213666303411;
    c2[4, 1] := 0.1494828185148637;
    c2[4, 2] := 0.6347841731714884;
    c2[5, 1] := 0.1394948812681826;
    c2[5, 2] := 0.5375594414619047;
    c2[6, 1] := 0.1273627583380806;
    c2[6, 2] := 0.4241608926375478;
    c2[7, 1] := 0.1129187258461290;
    c2[7, 2] := 0.2965752009703245;
    c2[8, 1] := 0.9533357359908857e-1;
    c2[8, 2] := 0.1537041700889585;
  elseif order == 18 then
    alpha := 0.2044575288651841;
    c2[1, 1] := 0.1588768571976356;
    c2[1, 2] := 0.7951914263212913;
    c2[2, 1] := 0.1569357024981854;
    c2[2, 2] := 0.7744529690772538;
    c2[3, 1] := 0.1530722206358810;
    c2[3, 2] := 0.7335304425992080;
    c2[4, 1] := 0.1473206710524167;
    c2[4, 2] := 0.6735038935387268;
    c2[5, 1] := 0.1397225420331520;
    c2[5, 2] := 0.5959151542621590;
    c2[6, 1] := 0.1303092459809849;
    c2[6, 2] := 0.5026483447894845;
    c2[7, 1] := 0.1190627367060072;
    c2[7, 2] := 0.3956893824587150;
    c2[8, 1] := 0.1058058030798994;
    c2[8, 2] := 0.2765091830730650;
    c2[9, 1] := 0.8974708108800873e-1;
    c2[9, 2] := 0.1435505288284833;
  elseif order == 19 then
    alpha := 0.1987936248083529;
    c1[1] := 0.3892259966869526;
    c2[1, 1] := 0.1506640012172225;
    c2[1, 2] := 0.7693121733774260;
    c2[2, 1] := 0.1481728062796673;
    c2[2, 2] := 0.7421133586741549;
    c2[3, 1] := 0.1440444668388838;
    c2[3, 2] := 0.6975075386214800;
    c2[4, 1] := 0.1383101628540374;
    c2[4, 2] := 0.6365464378910025;
    c2[5, 1] := 0.1310032283190998;
    c2[5, 2] := 0.5606211948462122;
    c2[6, 1] := 0.1221431166405330;
    c2[6, 2] := 0.4713530424221445;
    c2[7, 1] := 0.1116991161103884;
    c2[7, 2] := 0.3703717538617073;
    c2[8, 1] := 0.9948917351196349e-1;
    c2[8, 2] := 0.2587371155559744;
    c2[9, 1] := 0.8475989238107367e-1;
    c2[9, 2] := 0.1345537894555993;
  elseif order == 20 then
    alpha := 0.1935761760416219;
    c2[1, 1] := 0.1443871348337404;
    c2[1, 2] := 0.7584165598446141;
    c2[2, 1] := 0.1429501891353184;
    c2[2, 2] := 0.7423000962318863;
    c2[3, 1] := 0.1400877384920004;
    c2[3, 2] := 0.7104185332215555;
    c2[4, 1] := 0.1358210369491446;
    c2[4, 2] := 0.6634599783272630;
    c2[5, 1] := 0.1301773703034290;
    c2[5, 2] := 0.6024175491895959;
    c2[6, 1] := 0.1231826501439148;
    c2[6, 2] := 0.5285332736326852;
    c2[7, 1] := 0.1148465498575254;
    c2[7, 2] := 0.4431977385498628;
    c2[8, 1] := 0.1051289462376788;
    c2[8, 2] := 0.3477444062821162;
    c2[9, 1] := 0.9384622797485121e-1;
    c2[9, 2] := 0.2429038300327729;
    c2[10, 1] := 0.8028211612831444e-1;
    c2[10, 2] := 0.1265329974009533;
  elseif order == 21 then
    alpha := 0.1887494014766075;
    c1[1] := 0.3718070668941645;
    c2[1, 1] := 0.1376151928386445;
    c2[1, 2] := 0.7364290859445481;
    c2[2, 1] := 0.1357438914390695;
    c2[2, 2] := 0.7150167318935022;
    c2[3, 1] := 0.1326398453462415;
    c2[3, 2] := 0.6798001808470175;
    c2[4, 1] := 0.1283231214897678;
    c2[4, 2] := 0.6314663440439816;
    c2[5, 1] := 0.1228169159777534;
    c2[5, 2] := 0.5709353626166905;
    c2[6, 1] := 0.1161406100773184;
    c2[6, 2] := 0.4993087153571335;
    c2[7, 1] := 0.1082959649233524;
    c2[7, 2] := 0.4177766148584385;
    c2[8, 1] := 0.9923596957485723e-1;
    c2[8, 2] := 0.3274257287232124;
    c2[9, 1] := 0.8877776108724853e-1;
    c2[9, 2] := 0.2287218166767916;
    c2[10, 1] := 0.7624076527736326e-1;
    c2[10, 2] := 0.1193423971506988;
  elseif order == 22 then
    alpha := 0.1842668221199706;
    c2[1, 1] := 0.1323053462701543;
    c2[1, 2] := 0.7262446126765204;
    c2[2, 1] := 0.1312121721769772;
    c2[2, 2] := 0.7134286088450949;
    c2[3, 1] := 0.1290330911166814;
    c2[3, 2] := 0.6880287870435514;
    c2[4, 1] := 0.1257817990372067;
    c2[4, 2] := 0.6505015800059301;
    c2[5, 1] := 0.1214765261983008;
    c2[5, 2] := 0.6015107185211451;
    c2[6, 1] := 0.1161365140967959;
    c2[6, 2] := 0.5418983553698413;
    c2[7, 1] := 0.1097755171533100;
    c2[7, 2] := 0.4726370779831614;
    c2[8, 1] := 0.1023889478519956;
    c2[8, 2] := 0.3947439506537486;
    c2[9, 1] := 0.9392485861253800e-1;
    c2[9, 2] := 0.3090996703083202;
    c2[10, 1] := 0.8420273775456455e-1;
    c2[10, 2] := 0.2159561978556017;
    c2[11, 1] := 0.7257600023938262e-1;
    c2[11, 2] := 0.1128633732721116;
  elseif order == 23 then
    alpha := 0.1800893554453722;
    c1[1] := 0.3565232673929280;
    c2[1, 1] := 0.1266275171652706;
    c2[1, 2] := 0.7072778066734162;
    c2[2, 1] := 0.1251865227648538;
    c2[2, 2] := 0.6900676345785905;
    c2[3, 1] := 0.1227944815236645;
    c2[3, 2] := 0.6617011100576023;
    c2[4, 1] := 0.1194647013077667;
    c2[4, 2] := 0.6226432315773119;
    c2[5, 1] := 0.1152132989252356;
    c2[5, 2] := 0.5735222810625359;
    c2[6, 1] := 0.1100558598478487;
    c2[6, 2] := 0.5151027978024605;
    c2[7, 1] := 0.1040013558214886;
    c2[7, 2] := 0.4482410942032739;
    c2[8, 1] := 0.9704014176512626e-1;
    c2[8, 2] := 0.3738049984631116;
    c2[9, 1] := 0.8911683905758054e-1;
    c2[9, 2] := 0.2925028692588410;
    c2[10, 1] := 0.8005438265072295e-1;
    c2[10, 2] := 0.2044134600278901;
    c2[11, 1] := 0.6923832296800832e-1;
    c2[11, 2] := 0.1069984887283394;
  elseif order == 24 then
    alpha := 0.1761838665838427;
    c2[1, 1] := 0.1220804912720132;
    c2[1, 2] := 0.6978026874156063;
    c2[2, 1] := 0.1212296762358897;
    c2[2, 2] := 0.6874139794926736;
    c2[3, 1] := 0.1195328372961027;
    c2[3, 2] := 0.6667954259551859;
    c2[4, 1] := 0.1169990987333593;
    c2[4, 2] := 0.6362602049901176;
    c2[5, 1] := 0.1136409040480130;
    c2[5, 2] := 0.5962662188435553;
    c2[6, 1] := 0.1094722001757955;
    c2[6, 2] := 0.5474001634109253;
    c2[7, 1] := 0.1045052832229087;
    c2[7, 2] := 0.4903523180249535;
    c2[8, 1] := 0.9874509806025907e-1;
    c2[8, 2] := 0.4258751523524645;
    c2[9, 1] := 0.9217799943472177e-1;
    c2[9, 2] := 0.3547079765396403;
    c2[10, 1] := 0.8474633796250476e-1;
    c2[10, 2] := 0.2774145482392767;
    c2[11, 1] := 0.7627722381240495e-1;
    c2[11, 2] := 0.1939329108084139;
    c2[12, 1] := 0.6618645465422745e-1;
    c2[12, 2] := 0.1016670147947242;
  elseif order == 25 then
    alpha := 0.1725220521949266;
    c1[1] := 0.3429735385896000;
    c2[1, 1] := 0.1172525033170618;
    c2[1, 2] := 0.6812327932576614;
    c2[2, 1] := 0.1161194585333535;
    c2[2, 2] := 0.6671566071153211;
    c2[3, 1] := 0.1142375145794466;
    c2[3, 2] := 0.6439167855053158;
    c2[4, 1] := 0.1116157454252308;
    c2[4, 2] := 0.6118378416180135;
    c2[5, 1] := 0.1082654809459177;
    c2[5, 2] := 0.5713609763370088;
    c2[6, 1] := 0.1041985674230918;
    c2[6, 2] := 0.5230289949762722;
    c2[7, 1] := 0.9942439308123559e-1;
    c2[7, 2] := 0.4674627926041906;
    c2[8, 1] := 0.9394453593830893e-1;
    c2[8, 2] := 0.4053226688298811;
    c2[9, 1] := 0.8774221237222533e-1;
    c2[9, 2] := 0.3372372276379071;
    c2[10, 1] := 0.8075839512216483e-1;
    c2[10, 2] := 0.2636485508005428;
    c2[11, 1] := 0.7282483286646764e-1;
    c2[11, 2] := 0.1843801345273085;
    c2[12, 1] := 0.6338571166846652e-1;
    c2[12, 2] := 0.9680153764737715e-1;
  elseif order == 26 then
    alpha := 0.1690795702796737;
    c2[1, 1] := 0.1133168695796030;
    c2[1, 2] := 0.6724297955493932;
    c2[2, 1] := 0.1126417845769961;
    c2[2, 2] := 0.6638709519790540;
    c2[3, 1] := 0.1112948749545606;
    c2[3, 2] := 0.6468652038763624;
    c2[4, 1] := 0.1092823986944244;
    c2[4, 2] := 0.6216337070799265;
    c2[5, 1] := 0.1066130386697976;
    c2[5, 2] := 0.5885011413992190;
    c2[6, 1] := 0.1032969057045413;
    c2[6, 2] := 0.5478864278297548;
    c2[7, 1] := 0.9934388184210715e-1;
    c2[7, 2] := 0.5002885306054287;
    c2[8, 1] := 0.9476081523436283e-1;
    c2[8, 2] := 0.4462644847551711;
    c2[9, 1] := 0.8954648464575577e-1;
    c2[9, 2] := 0.3863930785049522;
    c2[10, 1] := 0.8368166847159917e-1;
    c2[10, 2] := 0.3212074592527143;
    c2[11, 1] := 0.7710664731701103e-1;
    c2[11, 2] := 0.2510470347119383;
    c2[12, 1] := 0.6965807988411425e-1;
    c2[12, 2] := 0.1756419294111342;
    c2[13, 1] := 0.6080674930548766e-1;
    c2[13, 2] := 0.9234535279274277e-1;
  elseif order == 27 then
    alpha := 0.1658353543067995;
    c1[1] := 0.3308543720638957;
    c2[1, 1] := 0.1091618578712746;
    c2[1, 2] := 0.6577977071169651;
    c2[2, 1] := 0.1082549561495043;
    c2[2, 2] := 0.6461121666520275;
    c2[3, 1] := 0.1067479247890451;
    c2[3, 2] := 0.6267937760991321;
    c2[4, 1] := 0.1046471079537577;
    c2[4, 2] := 0.6000750116745808;
    c2[5, 1] := 0.1019605976654259;
    c2[5, 2] := 0.5662734183049320;
    c2[6, 1] := 0.9869726954433709e-1;
    c2[6, 2] := 0.5257827234948534;
    c2[7, 1] := 0.9486520934132483e-1;
    c2[7, 2] := 0.4790595019077763;
    c2[8, 1] := 0.9046906518775348e-1;
    c2[8, 2] := 0.4266025862147336;
    c2[9, 1] := 0.8550529998276152e-1;
    c2[9, 2] := 0.3689188223512328;
    c2[10, 1] := 0.7995282239306020e-1;
    c2[10, 2] := 0.3064589322702932;
    c2[11, 1] := 0.7375174596252882e-1;
    c2[11, 2] := 0.2394754504667310;
    c2[12, 1] := 0.6674377263329041e-1;
    c2[12, 2] := 0.1676223546666024;
    c2[13, 1] := 0.5842458027529246e-1;
    c2[13, 2] := 0.8825044329219431e-1;
  elseif order == 28 then
    alpha := 0.1627710671942929;
    c2[1, 1] := 0.1057232656113488;
    c2[1, 2] := 0.6496161226860832;
    c2[2, 1] := 0.1051786825724864;
    c2[2, 2] := 0.6424661279909941;
    c2[3, 1] := 0.1040917964935006;
    c2[3, 2] := 0.6282470268918791;
    c2[4, 1] := 0.1024670101953951;
    c2[4, 2] := 0.6071189030701136;
    c2[5, 1] := 0.1003105109519892;
    c2[5, 2] := 0.5793175191747016;
    c2[6, 1] := 0.9762969425430802e-1;
    c2[6, 2] := 0.5451486608855443;
    c2[7, 1] := 0.9443223803058400e-1;
    c2[7, 2] := 0.5049796971628137;
    c2[8, 1] := 0.9072460982036488e-1;
    c2[8, 2] := 0.4592270546572523;
    c2[9, 1] := 0.8650956423253280e-1;
    c2[9, 2] := 0.4083368605952977;
    c2[10, 1] := 0.8178165740374893e-1;
    c2[10, 2] := 0.3527525188880655;
    c2[11, 1] := 0.7651838885868020e-1;
    c2[11, 2] := 0.2928534570013572;
    c2[12, 1] := 0.7066010532447490e-1;
    c2[12, 2] := 0.2288185204390681;
    c2[13, 1] := 0.6405358596145789e-1;
    c2[13, 2] := 0.1602396172588190;
    c2[14, 1] := 0.5621780070227172e-1;
    c2[14, 2] := 0.8447589564915071e-1;
  elseif order == 29 then
    alpha := 0.1598706626277596;
    c1[1] := 0.3199314513011623;
    c2[1, 1] := 0.1021101032532951;
    c2[1, 2] := 0.6365758882240111;
    c2[2, 1] := 0.1013729819392774;
    c2[2, 2] := 0.6267495975736321;
    c2[3, 1] := 0.1001476175660628;
    c2[3, 2] := 0.6104876178266819;
    c2[4, 1] := 0.9843854640428316e-1;
    c2[4, 2] := 0.5879603139195113;
    c2[5, 1] := 0.9625164534591696e-1;
    c2[5, 2] := 0.5594012291050210;
    c2[6, 1] := 0.9359356960417668e-1;
    c2[6, 2] := 0.5251016150410664;
    c2[7, 1] := 0.9047086748649986e-1;
    c2[7, 2] := 0.4854024475590397;
    c2[8, 1] := 0.8688856407189167e-1;
    c2[8, 2] := 0.4406826457109709;
    c2[9, 1] := 0.8284779224069856e-1;
    c2[9, 2] := 0.3913408089298914;
    c2[10, 1] := 0.7834154620997181e-1;
    c2[10, 2] := 0.3377643999400627;
    c2[11, 1] := 0.7334628941928766e-1;
    c2[11, 2] := 0.2802710651919946;
    c2[12, 1] := 0.6780290487362146e-1;
    c2[12, 2] := 0.2189770008083379;
    c2[13, 1] := 0.6156321231528423e-1;
    c2[13, 2] := 0.1534235999306070;
    c2[14, 1] := 0.5416797446761512e-1;
    c2[14, 2] := 0.8098664736760292e-1;
  elseif order == 30 then
    alpha := 0.1571200296252450;
    c2[1, 1] := 0.9908074847842124e-1;
    c2[1, 2] := 0.6289618807831557;
    c2[2, 1] := 0.9863509708328196e-1;
    c2[2, 2] := 0.6229164525571278;
    c2[3, 1] := 0.9774542692037148e-1;
    c2[3, 2] := 0.6108853364240036;
    c2[4, 1] := 0.9641490581986484e-1;
    c2[4, 2] := 0.5929869253412513;
    c2[5, 1] := 0.9464802912225441e-1;
    c2[5, 2] := 0.5693960175547550;
    c2[6, 1] := 0.9245027206218041e-1;
    c2[6, 2] := 0.5403402396359503;
    c2[7, 1] := 0.8982754584112941e-1;
    c2[7, 2] := 0.5060948065875106;
    c2[8, 1] := 0.8678535291732599e-1;
    c2[8, 2] := 0.4669749797983789;
    c2[9, 1] := 0.8332744242052199e-1;
    c2[9, 2] := 0.4233249626334694;
    c2[10, 1] := 0.7945356393775309e-1;
    c2[10, 2] := 0.3755006094498054;
    c2[11, 1] := 0.7515543969833788e-1;
    c2[11, 2] := 0.3238400339292700;
    c2[12, 1] := 0.7040879901685638e-1;
    c2[12, 2] := 0.2686072427439079;
    c2[13, 1] := 0.6515528854010540e-1;
    c2[13, 2] := 0.2098650589782619;
    c2[14, 1] := 0.5925168237177876e-1;
    c2[14, 2] := 0.1471138832654873;
    c2[15, 1] := 0.5225913954211672e-1;
    c2[15, 2] := 0.7775248839507864e-1;
  elseif order == 31 then
    alpha := 0.1545067022920929;
    c1[1] := 0.3100206996451866;
    c2[1, 1] := 0.9591020358831668e-1;
    c2[1, 2] := 0.6172474793293396;
    c2[2, 1] := 0.9530301275601203e-1;
    c2[2, 2] := 0.6088916323460413;
    c2[3, 1] := 0.9429332655402368e-1;
    c2[3, 2] := 0.5950511595503025;
    c2[4, 1] := 0.9288445429894548e-1;
    c2[4, 2] := 0.5758534119053522;
    c2[5, 1] := 0.9108073420087422e-1;
    c2[5, 2] := 0.5514734636081183;
    c2[6, 1] := 0.8888719137536870e-1;
    c2[6, 2] := 0.5221306199481831;
    c2[7, 1] := 0.8630901440239650e-1;
    c2[7, 2] := 0.4880834248148061;
    c2[8, 1] := 0.8335074993373294e-1;
    c2[8, 2] := 0.4496225358496770;
    c2[9, 1] := 0.8001502494376102e-1;
    c2[9, 2] := 0.4070602306679052;
    c2[10, 1] := 0.7630041338037624e-1;
    c2[10, 2] := 0.3607139804818122;
    c2[11, 1] := 0.7219760885744920e-1;
    c2[11, 2] := 0.3108783301229550;
    c2[12, 1] := 0.6768185077153345e-1;
    c2[12, 2] := 0.2577706252514497;
    c2[13, 1] := 0.6269571766328638e-1;
    c2[13, 2] := 0.2014081375889921;
    c2[14, 1] := 0.5710081766945065e-1;
    c2[14, 2] := 0.1412581515841926;
    c2[15, 1] := 0.5047740914807019e-1;
    c2[15, 2] := 0.7474725873250158e-1;
  elseif order == 32 then
    alpha := 0.1520196210848210;
    c2[1, 1] := 0.9322163554339406e-1;
    c2[1, 2] := 0.6101488690506050;
    c2[2, 1] := 0.9285233997694042e-1;
    c2[2, 2] := 0.6049832320721264;
    c2[3, 1] := 0.9211494244473163e-1;
    c2[3, 2] := 0.5946969295569034;
    c2[4, 1] := 0.9101176786042449e-1;
    c2[4, 2] := 0.5793791854364477;
    c2[5, 1] := 0.8954614071360517e-1;
    c2[5, 2] := 0.5591619969234026;
    c2[6, 1] := 0.8772216763680164e-1;
    c2[6, 2] := 0.5342177994699602;
    c2[7, 1] := 0.8554440426912734e-1;
    c2[7, 2] := 0.5047560942986598;
    c2[8, 1] := 0.8301735302045588e-1;
    c2[8, 2] := 0.4710187048140929;
    c2[9, 1] := 0.8014469519188161e-1;
    c2[9, 2] := 0.4332730387207936;
    c2[10, 1] := 0.7692807528893225e-1;
    c2[10, 2] := 0.3918021436411035;
    c2[11, 1] := 0.7336507157284898e-1;
    c2[11, 2] := 0.3468890521471250;
    c2[12, 1] := 0.6944555312763458e-1;
    c2[12, 2] := 0.2987898029050460;
    c2[13, 1] := 0.6514446669420571e-1;
    c2[13, 2] := 0.2476810747407199;
    c2[14, 1] := 0.6040544477732702e-1;
    c2[14, 2] := 0.1935412053397663;
    c2[15, 1] := 0.5509478650672775e-1;
    c2[15, 2] := 0.1358108994174911;
    c2[16, 1] := 0.4881064725720192e-1;
    c2[16, 2] := 0.7194819894416505e-1;
  elseif order == 33 then
    alpha := 0.1496489351138032;
    c1[1] := 0.3009752799176432;
    c2[1, 1] := 0.9041725460994505e-1;
    c2[1, 2] := 0.5995521047364046;
    c2[2, 1] := 0.8991117804113002e-1;
    c2[2, 2] := 0.5923764112099496;
    c2[3, 1] := 0.8906941547422532e-1;
    c2[3, 2] := 0.5804822013853129;
    c2[4, 1] := 0.8789442491445575e-1;
    c2[4, 2] := 0.5639663528946501;
    c2[5, 1] := 0.8638945831033775e-1;
    c2[5, 2] := 0.5429623519607796;
    c2[6, 1] := 0.8455834602616358e-1;
    c2[6, 2] := 0.5176379938389326;
    c2[7, 1] := 0.8240517431382334e-1;
    c2[7, 2] := 0.4881921474066189;
    c2[8, 1] := 0.7993380417355076e-1;
    c2[8, 2] := 0.4548502528082586;
    c2[9, 1] := 0.7714713890732801e-1;
    c2[9, 2] := 0.4178579388038483;
    c2[10, 1] := 0.7404596598181127e-1;
    c2[10, 2] := 0.3774715722484659;
    c2[11, 1] := 0.7062702339160462e-1;
    c2[11, 2] := 0.3339432938810453;
    c2[12, 1] := 0.6687952672391507e-1;
    c2[12, 2] := 0.2874950693388235;
    c2[13, 1] := 0.6277828912909767e-1;
    c2[13, 2] := 0.2382680702894708;
    c2[14, 1] := 0.5826808305383988e-1;
    c2[14, 2] := 0.1862073169968455;
    c2[15, 1] := 0.5321974125363517e-1;
    c2[15, 2] := 0.1307323751236313;
    c2[16, 1] := 0.4724820282032780e-1;
    c2[16, 2] := 0.6933542082177094e-1;
  elseif order == 34 then
    alpha := 0.1473858373968463;
    c2[1, 1] := 0.8801537152275983e-1;
    c2[1, 2] := 0.5929204288972172;
    c2[2, 1] := 0.8770594341007476e-1;
    c2[2, 2] := 0.5884653382247518;
    c2[3, 1] := 0.8708797598072095e-1;
    c2[3, 2] := 0.5795895850253119;
    c2[4, 1] := 0.8616320590689187e-1;
    c2[4, 2] := 0.5663615383647170;
    c2[5, 1] := 0.8493413175570858e-1;
    c2[5, 2] := 0.5488825092350877;
    c2[6, 1] := 0.8340387368687513e-1;
    c2[6, 2] := 0.5272851839324592;
    c2[7, 1] := 0.8157596213131521e-1;
    c2[7, 2] := 0.5017313864372913;
    c2[8, 1] := 0.7945402670834270e-1;
    c2[8, 2] := 0.4724089864574216;
    c2[9, 1] := 0.7704133559556429e-1;
    c2[9, 2] := 0.4395276256463053;
    c2[10, 1] := 0.7434009635219704e-1;
    c2[10, 2] := 0.4033126590648964;
    c2[11, 1] := 0.7135035113853376e-1;
    c2[11, 2] := 0.3639961488919042;
    c2[12, 1] := 0.6806813160738834e-1;
    c2[12, 2] := 0.3218025212900124;
    c2[13, 1] := 0.6448214312000864e-1;
    c2[13, 2] := 0.2769235521088158;
    c2[14, 1] := 0.6056719318430530e-1;
    c2[14, 2] := 0.2294693573271038;
    c2[15, 1] := 0.5626925196925040e-1;
    c2[15, 2] := 0.1793564218840015;
    c2[16, 1] := 0.5146352031547277e-1;
    c2[16, 2] := 0.1259877129326412;
    c2[17, 1] := 0.4578069074410591e-1;
    c2[17, 2] := 0.6689147319568768e-1;
  elseif order == 35 then
    alpha := 0.1452224267615486;
    c1[1] := 0.2926764667564367;
    c2[1, 1] := 0.8551731299267280e-1;
    c2[1, 2] := 0.5832758214629523;
    c2[2, 1] := 0.8509109732853060e-1;
    c2[2, 2] := 0.5770596582643844;
    c2[3, 1] := 0.8438201446671953e-1;
    c2[3, 2] := 0.5667497616665494;
    c2[4, 1] := 0.8339191981579831e-1;
    c2[4, 2] := 0.5524209816238369;
    c2[5, 1] := 0.8212328610083385e-1;
    c2[5, 2] := 0.5341766459916322;
    c2[6, 1] := 0.8057906332198853e-1;
    c2[6, 2] := 0.5121470053512750;
    c2[7, 1] := 0.7876247299954955e-1;
    c2[7, 2] := 0.4864870722254752;
    c2[8, 1] := 0.7667670879950268e-1;
    c2[8, 2] := 0.4573736721705665;
    c2[9, 1] := 0.7432449556218945e-1;
    c2[9, 2] := 0.4250013835198991;
    c2[10, 1] := 0.7170742126011575e-1;
    c2[10, 2] := 0.3895767735915445;
    c2[11, 1] := 0.6882488171701314e-1;
    c2[11, 2] := 0.3513097926737368;
    c2[12, 1] := 0.6567231746957568e-1;
    c2[12, 2] := 0.3103999917596611;
    c2[13, 1] := 0.6223804362223595e-1;
    c2[13, 2] := 0.2670123611280899;
    c2[14, 1] := 0.5849696460782910e-1;
    c2[14, 2] := 0.2212298104867592;
    c2[15, 1] := 0.5439628409499822e-1;
    c2[15, 2] := 0.1729443731341637;
    c2[16, 1] := 0.4981540179136920e-1;
    c2[16, 2] := 0.1215462157134930;
    c2[17, 1] := 0.4439981033536435e-1;
    c2[17, 2] := 0.6460098363520967e-1;
  elseif order == 36 then
    alpha := 0.1431515914458580;
    c2[1, 1] := 0.8335881847130301e-1;
    c2[1, 2] := 0.5770670512160201;
    c2[2, 1] := 0.8309698922852212e-1;
    c2[2, 2] := 0.5731929100172432;
    c2[3, 1] := 0.8257400347039723e-1;
    c2[3, 2] := 0.5654713811993058;
    c2[4, 1] := 0.8179117911600136e-1;
    c2[4, 2] := 0.5539556343603020;
    c2[5, 1] := 0.8075042173126963e-1;
    c2[5, 2] := 0.5387245649546684;
    c2[6, 1] := 0.7945413151258206e-1;
    c2[6, 2] := 0.5198817177723069;
    c2[7, 1] := 0.7790506514288866e-1;
    c2[7, 2] := 0.4975537629595409;
    c2[8, 1] := 0.7610613635339480e-1;
    c2[8, 2] := 0.4718884193866789;
    c2[9, 1] := 0.7406012816626425e-1;
    c2[9, 2] := 0.4430516443136726;
    c2[10, 1] := 0.7176927060205631e-1;
    c2[10, 2] := 0.4112237708115829;
    c2[11, 1] := 0.6923460172504251e-1;
    c2[11, 2] := 0.3765940116389730;
    c2[12, 1] := 0.6645495833489556e-1;
    c2[12, 2] := 0.3393522147815403;
    c2[13, 1] := 0.6342528888937094e-1;
    c2[13, 2] := 0.2996755899575573;
    c2[14, 1] := 0.6013361864949449e-1;
    c2[14, 2] := 0.2577053294053830;
    c2[15, 1] := 0.5655503081322404e-1;
    c2[15, 2] := 0.2135004731531631;
    c2[16, 1] := 0.5263798119559069e-1;
    c2[16, 2] := 0.1669320999865636;
    c2[17, 1] := 0.4826589873626196e-1;
    c2[17, 2] := 0.1173807590715484;
    c2[18, 1] := 0.4309819397289806e-1;
    c2[18, 2] := 0.6245036108880222e-1;
  elseif order == 37 then
    alpha := 0.1411669104782917;
    c1[1] := 0.2850271036215707;
    c2[1, 1] := 0.8111958235023328e-1;
    c2[1, 2] := 0.5682412610563970;
    c2[2, 1] := 0.8075727567979578e-1;
    c2[2, 2] := 0.5628142923227016;
    c2[3, 1] := 0.8015440554413301e-1;
    c2[3, 2] := 0.5538087696879930;
    c2[4, 1] := 0.7931239302677386e-1;
    c2[4, 2] := 0.5412833323304460;
    c2[5, 1] := 0.7823314328639347e-1;
    c2[5, 2] := 0.5253190555393968;
    c2[6, 1] := 0.7691895211595101e-1;
    c2[6, 2] := 0.5060183741977191;
    c2[7, 1] := 0.7537237072011853e-1;
    c2[7, 2] := 0.4835036020049034;
    c2[8, 1] := 0.7359601294804538e-1;
    c2[8, 2] := 0.4579149413954837;
    c2[9, 1] := 0.7159227884849299e-1;
    c2[9, 2] := 0.4294078049978829;
    c2[10, 1] := 0.6936295002846032e-1;
    c2[10, 2] := 0.3981491350382047;
    c2[11, 1] := 0.6690857785828917e-1;
    c2[11, 2] := 0.3643121502867948;
    c2[12, 1] := 0.6422751692085542e-1;
    c2[12, 2] := 0.3280684291406284;
    c2[13, 1] := 0.6131430866206096e-1;
    c2[13, 2] := 0.2895750997170303;
    c2[14, 1] := 0.5815677249570920e-1;
    c2[14, 2] := 0.2489521814805720;
    c2[15, 1] := 0.5473023527947980e-1;
    c2[15, 2] := 0.2062377435955363;
    c2[16, 1] := 0.5098441033167034e-1;
    c2[16, 2] := 0.1612849131645336;
    c2[17, 1] := 0.4680658811093562e-1;
    c2[17, 2] := 0.1134672937045305;
    c2[18, 1] := 0.4186928031694695e-1;
    c2[18, 2] := 0.6042754777339966e-1;
  elseif order == 38 then
    alpha := 0.1392625697140030;
    c2[1, 1] := 0.7916943373658329e-1;
    c2[1, 2] := 0.5624158631591745;
    c2[2, 1] := 0.7894592250257840e-1;
    c2[2, 2] := 0.5590219398777304;
    c2[3, 1] := 0.7849941672384930e-1;
    c2[3, 2] := 0.5522551628416841;
    c2[4, 1] := 0.7783093084875645e-1;
    c2[4, 2] := 0.5421574325808380;
    c2[5, 1] := 0.7694193770482690e-1;
    c2[5, 2] := 0.5287909941093643;
    c2[6, 1] := 0.7583430534712885e-1;
    c2[6, 2] := 0.5122376814029880;
    c2[7, 1] := 0.7451020436122948e-1;
    c2[7, 2] := 0.4925978555548549;
    c2[8, 1] := 0.7297197617673508e-1;
    c2[8, 2] := 0.4699889739625235;
    c2[9, 1] := 0.7122194706992953e-1;
    c2[9, 2] := 0.4445436860615774;
    c2[10, 1] := 0.6926216260386816e-1;
    c2[10, 2] := 0.4164072786327193;
    c2[11, 1] := 0.6709399961255503e-1;
    c2[11, 2] := 0.3857341621868851;
    c2[12, 1] := 0.6471757977022456e-1;
    c2[12, 2] := 0.3526828388476838;
    c2[13, 1] := 0.6213084287116965e-1;
    c2[13, 2] := 0.3174082831364342;
    c2[14, 1] := 0.5932799638550641e-1;
    c2[14, 2] := 0.2800495563550299;
    c2[15, 1] := 0.5629672408524944e-1;
    c2[15, 2] := 0.2407078154782509;
    c2[16, 1] := 0.5301264751544952e-1;
    c2[16, 2] := 0.1994026830553859;
    c2[17, 1] := 0.4942673259817896e-1;
    c2[17, 2] := 0.1559719194038917;
    c2[18, 1] := 0.4542996716979947e-1;
    c2[18, 2] := 0.1097844277878470;
    c2[19, 1] := 0.4070720755433961e-1;
    c2[19, 2] := 0.5852181110523043e-1;
  elseif order == 39 then
    alpha := 0.1374332900196804;
    c1[1] := 0.2779468246419593;
    c2[1, 1] := 0.7715084161825772e-1;
    c2[1, 2] := 0.5543001331300056;
    c2[2, 1] := 0.7684028301163326e-1;
    c2[2, 2] := 0.5495289890712267;
    c2[3, 1] := 0.7632343924866024e-1;
    c2[3, 2] := 0.5416083298429741;
    c2[4, 1] := 0.7560141319808483e-1;
    c2[4, 2] := 0.5305846713929198;
    c2[5, 1] := 0.7467569064745969e-1;
    c2[5, 2] := 0.5165224112570647;
    c2[6, 1] := 0.7354807648551346e-1;
    c2[6, 2] := 0.4995030679271456;
    c2[7, 1] := 0.7222060351121389e-1;
    c2[7, 2] := 0.4796242430956156;
    c2[8, 1] := 0.7069540462458585e-1;
    c2[8, 2] := 0.4569982440368368;
    c2[9, 1] := 0.6897453353492381e-1;
    c2[9, 2] := 0.4317502624832354;
    c2[10, 1] := 0.6705970959388781e-1;
    c2[10, 2] := 0.4040159353969854;
    c2[11, 1] := 0.6495194541066725e-1;
    c2[11, 2] := 0.3739379843169939;
    c2[12, 1] := 0.6265098412417610e-1;
    c2[12, 2] := 0.3416613843816217;
    c2[13, 1] := 0.6015440984955930e-1;
    c2[13, 2] := 0.3073260166338746;
    c2[14, 1] := 0.5745615876877304e-1;
    c2[14, 2] := 0.2710546723961181;
    c2[15, 1] := 0.5454383762391338e-1;
    c2[15, 2] := 0.2329316824061170;
    c2[16, 1] := 0.5139340231935751e-1;
    c2[16, 2] := 0.1929604256043231;
    c2[17, 1] := 0.4795705862458131e-1;
    c2[17, 2] := 0.1509655259246037;
    c2[18, 1] := 0.4412933231935506e-1;
    c2[18, 2] := 0.1063130748962878;
    c2[19, 1] := 0.3960672309405603e-1;
    c2[19, 2] := 0.5672356837211527e-1;
  elseif order == 40 then
    alpha := 0.1356742655825434;
    c2[1, 1] := 0.7538038374294594e-1;
    c2[1, 2] := 0.5488228264329617;
    c2[2, 1] := 0.7518806529402738e-1;
    c2[2, 2] := 0.5458297722483311;
    c2[3, 1] := 0.7480383050347119e-1;
    c2[3, 2] := 0.5398604576730540;
    c2[4, 1] := 0.7422847031965465e-1;
    c2[4, 2] := 0.5309482987446206;
    c2[5, 1] := 0.7346313704205006e-1;
    c2[5, 2] := 0.5191429845322307;
    c2[6, 1] := 0.7250930053201402e-1;
    c2[6, 2] := 0.5045099368431007;
    c2[7, 1] := 0.7136868456879621e-1;
    c2[7, 2] := 0.4871295553902607;
    c2[8, 1] := 0.7004317764946634e-1;
    c2[8, 2] := 0.4670962098860498;
    c2[9, 1] := 0.6853470921527828e-1;
    c2[9, 2] := 0.4445169164956202;
    c2[10, 1] := 0.6684507689945471e-1;
    c2[10, 2] := 0.4195095960479698;
    c2[11, 1] := 0.6497570123412630e-1;
    c2[11, 2] := 0.3922007419030645;
    c2[12, 1] := 0.6292726794917847e-1;
    c2[12, 2] := 0.3627221993494397;
    c2[13, 1] := 0.6069918741663154e-1;
    c2[13, 2] := 0.3312065181294388;
    c2[14, 1] := 0.5828873983769410e-1;
    c2[14, 2] := 0.2977798532686911;
    c2[15, 1] := 0.5568964389813015e-1;
    c2[15, 2] := 0.2625503293999835;
    c2[16, 1] := 0.5288947816690705e-1;
    c2[16, 2] := 0.2255872486520188;
    c2[17, 1] := 0.4986456327645859e-1;
    c2[17, 2] := 0.1868796731919594;
    c2[18, 1] := 0.4656832613054458e-1;
    c2[18, 2] := 0.1462410193532463;
    c2[19, 1] := 0.4289867647614935e-1;
    c2[19, 2] := 0.1030361558710747;
    c2[20, 1] := 0.3856310684054106e-1;
    c2[20, 2] := 0.5502423832293889e-1;
  elseif order == 41 then
    alpha := 0.1339811106984253;
    c1[1] := 0.2713685065531391;
    c2[1, 1] := 0.7355140275160984e-1;
    c2[1, 2] := 0.5413274778282860;
    c2[2, 1] := 0.7328319082267173e-1;
    c2[2, 2] := 0.5371064088294270;
    c2[3, 1] := 0.7283676160772547e-1;
    c2[3, 2] := 0.5300963437270770;
    c2[4, 1] := 0.7221298133014343e-1;
    c2[4, 2] := 0.5203345998371490;
    c2[5, 1] := 0.7141302173623395e-1;
    c2[5, 2] := 0.5078728971879841;
    c2[6, 1] := 0.7043831559982149e-1;
    c2[6, 2] := 0.4927768111819803;
    c2[7, 1] := 0.6929049381827268e-1;
    c2[7, 2] := 0.4751250308594139;
    c2[8, 1] := 0.6797129849758392e-1;
    c2[8, 2] := 0.4550083840638406;
    c2[9, 1] := 0.6648246325101609e-1;
    c2[9, 2] := 0.4325285673076087;
    c2[10, 1] := 0.6482554675958526e-1;
    c2[10, 2] := 0.4077964789091151;
    c2[11, 1] := 0.6300169683004558e-1;
    c2[11, 2] := 0.3809299858742483;
    c2[12, 1] := 0.6101130648543355e-1;
    c2[12, 2] := 0.3520508315700898;
    c2[13, 1] := 0.5885349417435808e-1;
    c2[13, 2] := 0.3212801560701271;
    c2[14, 1] := 0.5652528148656809e-1;
    c2[14, 2] := 0.2887316252774887;
    c2[15, 1] := 0.5402021575818373e-1;
    c2[15, 2] := 0.2545001287790888;
    c2[16, 1] := 0.5132588802608274e-1;
    c2[16, 2] := 0.2186415296842951;
    c2[17, 1] := 0.4841900639702602e-1;
    c2[17, 2] := 0.1811322622296060;
    c2[18, 1] := 0.4525419574485134e-1;
    c2[18, 2] := 0.1417762065404688;
    c2[19, 1] := 0.4173260173087802e-1;
    c2[19, 2] := 0.9993834530966510e-1;
    c2[20, 1] := 0.3757210572966463e-1;
    c2[20, 2] := 0.5341611499960143e-1;
  else
    Streams.error("Input argument order (= " + String(order) +
      ") of Bessel filter is not in the range 1..41");
  end if;

end BesselBaseCoefficients;

Modelica.Blocks.Continuous.Internal.Filter.Utilities.toHighestPowerOne

Transform filter to form with highest power of s equal 1

Inputs

Type Name Default Description

Real den1[:] [s] coefficients of polynomials (den1[i]*s + 1)

Real den2[:, 2] [s^2, s] coefficients of polynomials (den2[i,1]*s^2 + den2[i,2]*s + 1)

Type	Name	Default	Description
Real	den1[:]		[s] coefficients of polynomials (den1[i]*s + 1)
Real	den2[:, 2]		[s^2, s] coefficients of polynomials (den2[i,1]s^2 + den2[i,2]s + 1)

Outputs

Type Name Description

Real cr[size(den1, 1)] [s^0] coefficients of polynomials cr[i]*(s+1/cr[i])

Real c0[size(den2, 1)] [s^0] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))

Real c1[size(den2, 1)] [s^1] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))

Type	Name	Description
Real	cr[size(den1, 1)]	[s^0] coefficients of polynomials cr[i]*(s+1/cr[i])
Real	c0[size(den2, 1)]	[s^0] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))
Real	c1[size(den2, 1)]	[s^1] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))

Modelica definition

function toHighestPowerOne 
  "Transform filter to form with highest power of s equal 1"

  input Real den1[:] "[s] coefficients of polynomials (den1[i]*s + 1)";
  input Real den2[:,2] 
    "[s^2, s] coefficients of polynomials (den2[i,1]*s^2 + den2[i,2]*s + 1)";
  output Real cr[size(den1, 1)] 
    "[s^0] coefficients of polynomials cr[i]*(s+1/cr[i])";
  output Real c0[size(den2, 1)] 
    "[s^0] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))";
  output Real c1[size(den2, 1)] 
    "[s^1] coefficients of polynomials (s^2 + (den2[i,2]/den2[i,1])*s + (1/den2[i,1]))";
algorithm 
  for i in 1:size(den1, 1) loop
    cr[i] := 1/den1[i];
  end for;

  for i in 1:size(den2, 1) loop
    c1[i] := den2[i, 2]/den2[i, 1];
    c0[i] := 1/den2[i, 1];
  end for;
end toHighestPowerOne;

Modelica.Blocks.Continuous.Internal.Filter.Utilities.normalizationFactor

Compute correction factor of low pass filter such that amplitude at cut-off frequency is -3db (=10^(-3/20) = 0.70794...)

Inputs

Type Name Default Description

Real c1[:] [p] coefficients of denominator polynomials (c1[i}*p + 1)

Real c2[:, 2] [p^2, p] coefficients of denominator polynomials (c2[i,1]*p^2 + c2[i,2]*p + 1)

Type	Name	Default	Description
Real	c1[:]		[p] coefficients of denominator polynomials (c1[i}*p + 1)
Real	c2[:, 2]		[p^2, p] coefficients of denominator polynomials (c2[i,1]p^2 + c2[i,2]p + 1)

Outputs

Type Name Description

Real alpha Correction factor (replace p by alpha*p)

Type	Name	Description
Real	alpha	Correction factor (replace p by alpha*p)

Modelica definition

function normalizationFactor 
  "Compute correction factor of low pass filter such that amplitude at cut-off frequency is -3db (=10^(-3/20) = 0.70794...)"
  import Modelica;
  import Modelica.Utilities.Streams;

  input Real c1[:] "[p] coefficients of denominator polynomials (c1[i}*p + 1)";
  input Real c2[:,2] 
    "[p^2, p] coefficients of denominator polynomials (c2[i,1]*p^2 + c2[i,2]*p + 1)";
  output Real alpha "Correction factor (replace p by alpha*p)";
protected 
  Real alpha_min;
  Real alpha_max;

  function normalizationResidue 
    "Residue of correction factor computation"
    input Real c1[:] 
      "[p] coefficients of denominator polynomials (c1[i]*p + 1)";
    input Real c2[:,2] 
      "[p^2, p] coefficients of denominator polynomials (c2[i,1]*p^2 + c2[i,2]*p + 1)";
    input Real alpha;
    output Real residue;
  protected 
    constant Real beta= 10^(-3/20) 
      "Amplitude of -3db required, i.e., -3db = 20*log(beta)";
    Real cc1;
    Real cc2;
    Real p;
    Real alpha2=alpha*alpha;
    Real alpha4=alpha2*alpha2;
    Real A2=1.0;
  algorithm 
    assert(size(c1,1) <= 1, "Internal error 2 (should not occur)");
    if size(c1, 1) == 1 then
      cc1 := c1[1]*c1[1];
      p := 1 + cc1*alpha2;
      A2 := A2*p;
    end if;
    for i in 1:size(c2, 1) loop
      cc1 := c2[i, 2]*c2[i, 2] - 2*c2[i, 1];
      cc2 := c2[i, 1]*c2[i, 1];
      p := 1 + cc1*alpha2 + cc2*alpha4;
      A2 := A2*p;
    end for;
    residue := 1/sqrt(A2) - beta;
  end normalizationResidue;

  function findInterval "Find interval for the root"
    input Real c1[:] "[p] coefficients of denominator polynomials (a*p + 1)";
    input Real c2[:,2] 
      "[p^2, p] coefficients of denominator polynomials (b*p^2 + a*p + 1)";
    output Real alpha_min;
    output Real alpha_max;
  protected 
    Real alpha = 1.0;
    Real residue;
  algorithm 
    alpha_min :=0;
    residue := normalizationResidue(c1, c2, alpha);
    if residue < 0 then
       alpha_max :=alpha;
    else
       while residue >= 0 loop
          alpha := 1.1*alpha;
          residue := normalizationResidue(c1, c2, alpha);
       end while;
       alpha_max :=alpha;
    end if;
  end findInterval;

function solveOneNonlinearEquation 
    "Solve f(u) = 0; f(u_min) and f(u_max) must have different signs"
    import Modelica.Utilities.Streams.error;

  input Real c1[:] "[p] coefficients of denominator polynomials (c1[i]*p + 1)";
  input Real c2[:,2] 
      "[p^2, p] coefficients of denominator polynomials (c2[i,1]*p^2 + c2[i,2]*p + 1)";
  input Real u_min "Lower bound of search intervall";
  input Real u_max "Upper bound of search intervall";
  input Real tolerance=100*Modelica.Constants.eps 
      "Relative tolerance of solution u";
  output Real u "Value of independent variable so that f(u) = 0";

  protected 
  constant Real eps=Modelica.Constants.eps "machine epsilon";
  Real a=u_min "Current best minimum interval value";
  Real b=u_max "Current best maximum interval value";
  Real c "Intermediate point a <= c <= b";
  Real d;
  Real e "b - a";
  Real m;
  Real s;
  Real p;
  Real q;
  Real r;
  Real tol;
  Real fa "= f(a)";
  Real fb "= f(b)";
  Real fc;
  Boolean found=false;
algorithm 
  // Check that f(u_min) and f(u_max) have different sign
  fa := normalizationResidue(c1,c2,u_min);
  fb := normalizationResidue(c1,c2,u_max);
  fc := fb;
  if fa > 0.0 and fb > 0.0 or fa < 0.0 and fb < 0.0 then
    error(
      "The arguments u_min and u_max to solveOneNonlinearEquation(..)\n" +
      "do not bracket the root of the single non-linear equation:\n" +
      "  u_min  = " + String(u_min) + "\n" + "  u_max  = " + String(u_max)
       + "\n" + "  fa = f(u_min) = " + String(fa) + "\n" +
      "  fb = f(u_max) = " + String(fb) + "\n" +
      "fa and fb must have opposite sign which is not the case");
  end if;

  // Initialize variables
  c := a;
  fc := fa;
  e := b - a;
  d := e;

  // Search loop
  while not found loop
    if abs(fc) < abs(fb) then
      a := b;
      b := c;
      c := a;
      fa := fb;
      fb := fc;
      fc := fa;
    end if;

    tol := 2*eps*abs(b) + tolerance;
    m := (c - b)/2;

    if abs(m) <= tol or fb == 0.0 then
      // root found (interval is small enough)
      found := true;
      u := b;
    else
      // Determine if a bisection is needed
      if abs(e) < tol or abs(fa) <= abs(fb) then
        e := m;
        d := e;
      else
        s := fb/fa;
        if a == c then
          // linear interpolation
          p := 2*m*s;
          q := 1 - s;
        else
          // inverse quadratic interpolation
          q := fa/fc;
          r := fb/fc;
          p := s*(2*m*q*(q - r) - (b - a)*(r - 1));
          q := (q - 1)*(r - 1)*(s - 1);
        end if;

        if p > 0 then
          q := -q;
        else
          p := -p;
        end if;

        s := e;
        e := d;
        if 2*p < 3*m*q - abs(tol*q) and p < abs(0.5*s*q) then
          // interpolation successful
          d := p/q;
        else
          // use bi-section
          e := m;
          d := e;
        end if;
      end if;

      // Best guess value is defined as "a"
      a := b;
      fa := fb;
      b := b + (if abs(d) > tol then d else if m > 0 then tol else -tol);
      fb := normalizationResidue(c1,c2,b);

      if fb > 0 and fc > 0 or fb < 0 and fc < 0 then
        // initialize variables
        c := a;
        fc := fa;
        e := b - a;
        d := e;
      end if;
    end if;
  end while;

end solveOneNonlinearEquation;

algorithm 
   // Find interval for alpha
   (alpha_min, alpha_max) :=findInterval(c1, c2);

   // Compute alpha, so that abs(G(p)) = -3db
   alpha :=solveOneNonlinearEquation(
    c1,
    c2,
    alpha_min,
    alpha_max);
end normalizationFactor;

Modelica.Blocks.Continuous.Internal.Filter.Utilities.bandPassAlpha

Return alpha for band pass

Information

A band pass with bandwidth "w" is determined from a low pass

  1/(p^2 + a*p + b)

with the transformation

  new(p) = (p + 1/p)/w

This results in the following derivation:

  1/(p^2 + a*p + b) -> 1/( (p+1/p)^2/w^2 + a*(p + 1/p)/w + b )
                     = 1 / ( p^2 + 1/p^2 + 2)/w^2 + (p + 1/p)*a/w + b )
                     = w^2*p^2 / (p^4 + 2*p^2 + 1 + (p^3 + p)a*w + b*w^2*p^2)
                     = w^2*p^2 / (p^4 + a*w*p^3 + (2+b*w^2)*p^2 + a*w*p + 1)

This 4th order transfer function shall be split in to two transfer functions of order 2 each for numerical reasons. With the following formulation, the fourth order polynomial can be represented (with the unknowns "c" and "alpha"):

  g(p) = w^2*p^2 / ( (p*alpha)^2 + c*(p*alpha) + 1) * (p/alpha)^2 + c*(p/alpha) + 1)
       = w^2*p^2 / ( p^4 + c*(alpha + 1/alpha)*p^3 + (alpha^2 + 1/alpha^2 + c^2)*p^2
                                                   + c*(alpha + 1/alpha)*p + 1 )

Comparison of coefficients:

  c*(alpha + 1/alpha) = a*w           -> c = a*w / (alpha + 1/alpha)
  alpha^2 + 1/alpha^2 + c^2 = 2+b*w^2 -> equation to determine alpha

  alpha^4 + 1 + a^2*w^2*alpha^4/(1+alpha^2)^2 = (2+b*w^2)*alpha^2
    or z = alpha^2
  z^2 + a^2*w^2*z^2/(1+z)^2 - (2+b*w^2)*z + 1 = 0

Therefore the last equation has to be solved for "z" (basically, this means to compute a real zero of a fourth order polynomal):

   solve: 0 = f(z)  = z^2 + a^2*w^2*z^2/(1+z)^2 - (2+b*w^2)*z + 1  for "z"
              f(0)  = 1  > 0
              f(1)  = 1 + a^2*w^2/4 - (2+b*w^2) + 1
                    = (a^2/4 - b)*w^2  < 0
                    // since b - a^2/4 > 0 requirement for complex conjugate poles
   -> 0 < z < 1

This function computes the solution of this equation and returns "alpha = sqrt(z)";

Inputs

Type Name Default Description

Real a Coefficient of s^1

Real b Coefficient of s^0

AngularVelocity w Bandwidth angular frequency [rad/s]

Type	Name	Description
Real	a	Coefficient of s^1
Real	b	Coefficient of s^0
AngularVelocity	w	Bandwidth angular frequency [rad/s]

Outputs

Type Name Description

Real alpha Alpha factor to build up band pass

Type	Name	Description
Real	alpha	Alpha factor to build up band pass

Modelica definition

encapsulated function bandPassAlpha "Return alpha for band pass"
  import Modelica;
   input Real a "Coefficient of s^1";
   input Real b "Coefficient of s^0";
   input Modelica.SIunits.AngularVelocity w "Bandwidth angular frequency";
   output Real alpha "Alpha factor to build up band pass";

protected 
  Real alpha_min;
  Real alpha_max;
  Real z_min;
  Real z_max;
  Real z;

  function residue "Residue of non-linear equation"
    input Real a;
    input Real b;
    input Real w;
    input Real z;
    output Real res;
  algorithm 
    res := z^2 + (a*w*z/(1+z))^2 - (2+b*w^2)*z + 1;
  end residue;

function solveOneNonlinearEquation 
    "Solve f(u) = 0; f(u_min) and f(u_max) must have different signs"
    import Modelica.Utilities.Streams.error;

  input Real aa;
  input Real bb;
  input Real ww;
  input Real u_min "Lower bound of search intervall";
  input Real u_max "Upper bound of search intervall";
  input Real tolerance=100*Modelica.Constants.eps 
      "Relative tolerance of solution u";
  output Real u "Value of independent variable so that f(u) = 0";

  protected 
  constant Real eps=Modelica.Constants.eps "machine epsilon";
  Real a=u_min "Current best minimum interval value";
  Real b=u_max "Current best maximum interval value";
  Real c "Intermediate point a <= c <= b";
  Real d;
  Real e "b - a";
  Real m;
  Real s;
  Real p;
  Real q;
  Real r;
  Real tol;
  Real fa "= f(a)";
  Real fb "= f(b)";
  Real fc;
  Boolean found=false;
algorithm 
  // Check that f(u_min) and f(u_max) have different sign
  fa := residue(aa,bb,ww,u_min);
  fb := residue(aa,bb,ww,u_max);
  fc := fb;
  if fa > 0.0 and fb > 0.0 or fa < 0.0 and fb < 0.0 then
    error(
      "The arguments u_min and u_max to solveOneNonlinearEquation(..)\n" +
      "do not bracket the root of the single non-linear equation:\n" +
      "  u_min  = " + String(u_min) + "\n" + "  u_max  = " + String(u_max)
       + "\n" + "  fa = f(u_min) = " + String(fa) + "\n" +
      "  fb = f(u_max) = " + String(fb) + "\n" +
      "fa and fb must have opposite sign which is not the case");
  end if;

  // Initialize variables
  c := a;
  fc := fa;
  e := b - a;
  d := e;

  // Search loop
  while not found loop
    if abs(fc) < abs(fb) then
      a := b;
      b := c;
      c := a;
      fa := fb;
      fb := fc;
      fc := fa;
    end if;

    tol := 2*eps*abs(b) + tolerance;
    m := (c - b)/2;

    if abs(m) <= tol or fb == 0.0 then
      // root found (interval is small enough)
      found := true;
      u := b;
    else
      // Determine if a bisection is needed
      if abs(e) < tol or abs(fa) <= abs(fb) then
        e := m;
        d := e;
      else
        s := fb/fa;
        if a == c then
          // linear interpolation
          p := 2*m*s;
          q := 1 - s;
        else
          // inverse quadratic interpolation
          q := fa/fc;
          r := fb/fc;
          p := s*(2*m*q*(q - r) - (b - a)*(r - 1));
          q := (q - 1)*(r - 1)*(s - 1);
        end if;

        if p > 0 then
          q := -q;
        else
          p := -p;
        end if;

        s := e;
        e := d;
        if 2*p < 3*m*q - abs(tol*q) and p < abs(0.5*s*q) then
          // interpolation successful
          d := p/q;
        else
          // use bi-section
          e := m;
          d := e;
        end if;
      end if;

      // Best guess value is defined as "a"
      a := b;
      fa := fb;
      b := b + (if abs(d) > tol then d else if m > 0 then tol else -tol);
      fb := residue(aa,bb,ww,b);

      if fb > 0 and fc > 0 or fb < 0 and fc < 0 then
        // initialize variables
        c := a;
        fc := fa;
        e := b - a;
        d := e;
      end if;
    end if;
  end while;

end solveOneNonlinearEquation;

algorithm 
  assert( a^2/4 - b <= 0,  "Band pass transformation cannot be computed");
  z :=solveOneNonlinearEquation(a, b, w, 0, 1);
  alpha := sqrt(z);

end bandPassAlpha;

Automatically generated Fri Nov 12 16:27:36 2010.