1 
2 /* move from Deconvolve.c into libmri.a         21 Jun 2010 [rickr] */
3 
legendre(double x,int m)4 double legendre( double x , int m )   /* Legendre polynomials over [-1,1] */
5 {
6    if( m < 0 ) return 1.0 ;    /* bad input */
7 
8    switch( m ){                /*** P_m(x) for m=0..20 ***/
9     case 0: return 1.0 ;
10     case 1: return x ;
11     case 2: return (3.0*x*x-1.0)/2.0 ;
12     case 3: return (5.0*x*x-3.0)*x/2.0 ;
13     case 4: return ((35.0*x*x-30.0)*x*x+3.0)/8.0 ;
14     case 5: return ((63.0*x*x-70.0)*x*x+15.0)*x/8.0 ;
15     case 6: return (((231.0*x*x-315.0)*x*x+105.0)*x*x-5.0)/16.0 ;
16     case 7: return (((429.0*x*x-693.0)*x*x+315.0)*x*x-35.0)*x/16.0 ;
17     case 8: return ((((6435.0*x*x-12012.0)*x*x+6930.0)*x*x-1260.0)*x*x+35.0)/128.0;
18 
19            /** 07 Feb 2005: this part generated by Maple, then hand massaged **/
20 
21     case 9:
22       return (0.24609375e1 + (-0.3609375e2 + (0.140765625e3 + (-0.20109375e3
23               + 0.949609375e2 * x * x) * x * x) * x * x) * x * x) * x;
24 
25     case 10:
26       return -0.24609375e0 + (0.1353515625e2 + (-0.1173046875e3 +
27               (0.3519140625e3 + (-0.42732421875e3 + 0.18042578125e3 * x * x)
28              * x * x) * x * x) * x * x) * x * x;
29 
30     case 11:
31       return (-0.270703125e1 + (0.5865234375e2 + (-0.3519140625e3 +
32              (0.8546484375e3 + (-0.90212890625e3 + 0.34444921875e3 * x * x)
33              * x * x) * x * x) * x * x) * x * x) * x;
34 
35     case 12:
36       return 0.2255859375e0 + (-0.17595703125e2 + (0.2199462890625e3 +
37              (-0.99708984375e3 + (0.20297900390625e4 + (-0.1894470703125e4
38              + 0.6601943359375e3 * x * x) * x * x) * x * x) * x * x) * x * x)
39              * x * x;
40 
41     case 13:
42       return (0.29326171875e1 + (-0.87978515625e2 + (0.7478173828125e3 +
43              (-0.270638671875e4 + (0.47361767578125e4 + (-0.3961166015625e4
44              + 0.12696044921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
45             * x * x) * x;
46 
47     case 14:
48       return -0.20947265625e0 + (0.2199462890625e2 + (-0.37390869140625e3 +
49              (0.236808837890625e4 + (-0.710426513671875e4 +
50              (0.1089320654296875e5 + (-0.825242919921875e4 +
51             0.244852294921875e4 * x * x) * x * x) * x * x) * x * x) * x * x)
52            * x * x) * x * x;
53 
54     case 15:
55       return (-0.314208984375e1 + (0.12463623046875e3 + (-0.142085302734375e4
56             + (0.710426513671875e4 + (-0.1815534423828125e5 +
57               (0.2475728759765625e5 + (-0.1713966064453125e5 +
58                0.473381103515625e4 * x * x) * x * x) * x * x) * x * x)
59              * x * x) * x * x) * x * x) * x;
60 
61     case 16:
62       return 0.196380615234375e0 + (-0.26707763671875e2 + (0.5920220947265625e3
63             + (-0.4972985595703125e4 + (0.2042476226806641e5 +
64               (-0.4538836059570312e5 + (0.5570389709472656e5 +
65                (-0.3550358276367188e5 + 0.9171758880615234e4 * x * x) * x * x)
66             * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
67 
68     case 17:
69       return (0.3338470458984375e1 + (-0.169149169921875e3 +
70              (0.2486492797851562e4 + (-0.1633980981445312e5 +
71              (0.5673545074462891e5 + (-0.1114077941894531e6 +
72              (0.1242625396728516e6 + (-0.7337407104492188e5 +
73               0.1780400253295898e5 * x * x) * x * x) * x * x) * x * x)
74            * x * x) * x * x) * x * x) * x * x) * x;
75 
76     case 18:
77       return -0.1854705810546875e0 + (0.3171546936035156e2 +
78             (-0.8880331420898438e3 + (0.9531555725097656e4 +
79             (-0.5106190567016602e5 + (0.153185717010498e6 +
80             (-0.2692355026245117e6 + (0.275152766418457e6 +
81             (-0.1513340215301514e6 + 0.3461889381408691e5 * x * x) * x * x)
82            * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x;
83 
84     case 19:
85       return (-0.3523941040039062e1 + (0.2220082855224609e3 +
86              (-0.4084952453613281e4 + (0.3404127044677734e5 +
87              (-0.153185717010498e6 + (0.4038532539367676e6 +
88              (-0.6420231216430664e6 + (0.6053360861206055e6 +
89              (-0.3115700443267822e6 + 0.6741574058532715e5 * x * x) * x * x)
90           * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x * x) * x;
91 
92     case 20:
93       return 0.1761970520019531e0 + (-0.3700138092041016e2 +
94             (0.127654764175415e4 + (-0.1702063522338867e5 +
95             (0.1148892877578735e6 + (-0.4442385793304443e6 +
96             (0.1043287572669983e7 + (-0.1513340215301514e7 +
97             (0.1324172688388824e7 + (-0.6404495355606079e6 +
98              0.1314606941413879e6 * x * x) * x * x) * x * x) * x * x) * x * x)
99             * x * x) * x * x) * x * x) * x * x) * x * x;
100    }
101 
102 #if 0
103    /* order out of range: return Chebyshev instead (it's easy) */
104 
105         if(  x >=  1.0 ) x = 0.0 ;
106    else if ( x <= -1.0 ) x = 3.14159265358979323846 ;
107    else                  x = acos(x) ;
108    return cos(m*x) ;
109 #else
110    /** if here, m > 20 ==> use recurrence relation **/
111 
112    { double pk=0, pkm1, pkm2 ; int k ;
113      pkm2 = legendre( x , 19 ) ;
114      pkm1 = legendre( x , 20 ) ;
115      for( k=21 ; k <= m ; k++ , pkm2=pkm1 , pkm1=pk )
116        pk = ((2.0*k-1.0)*x*pkm1 - (k-1.0)*pkm2)/k ;
117      return pk ;
118    }
119 #endif
120 }
121