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