cal/cal/gauss.cal

{
        Gaussian Integration

                10/20/87

        order                   - order of integral (# of sample points)
        xmin                    - interval begin
        xmax                    - interval end
        f(x)                    - integrand
        integ                   - integral value
}
                        { Defaults }
order = 24 ;
xmin = -1 ;
xmax = 1 ;

integ = q1(ceil(order/2-.25)) * (xmax-xmin)/2 ;

q1(i) = if( i-1.5, wght(i)*(f1(absc(i))+f1(-absc(i))) + q1(i-1),
                wght(1)*(f1(absc(1))+if(even,f1(-absc(1)),0)) ) ;

even = floor(order/2+.25) - order/2 + .25 ;

f1(u) = f((xmax*(1+u)+xmin*(1-u))/2) ;

                        { Abscissas }
absc(i) =       select( order,
                        0,
                        .577350269189626,               { n = 2 }
                        select( i,                      { n = 3 }
                                0,
                                .774596669241483
                        ),
                        0, 0,
                        select( i,                      { n = 6 }
                                .238619186083197,
                                .661209386466265,
                                .932469514203152
                        ),
                        0, 0, 0, 0, 0,
                        select( i,                      { n = 12 }
                                .125233408511469,
                                .367831498998180,
                                .587317954286617,
                                .769902674194305,
                                .904117256370475,
                                .981560634246719
                        ),
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        select( i,                      { n = 24 }
                                .064056892862606,
                                .191118867473616,
                                .315042679696163,
                                .433793507626045,
                                .545421471388840,
                                .648093651936976,
                                .740124191578554,
                                .820001985973903,
                                .886415527004401,
                                .938274552002733,
                                .974728555971309,
                                .995187219997021
                        )
                ) ;

                        { Weights }
wght(i) =       select( order,
                        0,
                        1,                              { n = 2 }
                        select( i,                      { n = 3 }
                                .888888888888888,
                                .555555555555555
                        ),
                        0, 0,
                        select( i,                      { n = 6 }
                                .467913934572691,
                                .360761573048139,
                                .171324492379170
                        ),
                        0, 0, 0, 0, 0,
                        select( i,                      { n = 12 }
                                .249147045813403,
                                .233492536538355,
                                .203167426723066,
                                .160078328543346,
                                .106939325995318,
                                .047175336386512
                        ),
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        select( i,                      { n = 24 }
                                .127938195346752,
                                .125837456346828,
                                .121670472927803,
                                .115505668053726,
                                .107444270115966,
                                .097618652104114,
                                .086190161531953,
                                .073346481411080,
                                .059298584915437,
                                .044277438817420,
                                .028531388628934,
                                .012341229799987
                        )
                ) ;
Revision:	1.1
Committed:	Sat Feb 22 02:07:21 2003 UTC (21 years, 2 months ago) by greg
Branch:	MAIN
CVS Tags:	rad5R2, rad4R2P2, rad5R0, rad5R1, rad3R7P2, rad3R7P1, rad4R2, rad4R1, rad4R0, rad3R5, rad3R6, rad3R6P1, rad3R8, rad3R9, rad4R2P1
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	1.1	{
2			Gaussian Integration
3
4			10/20/87
5
6			order - order of integral (# of sample points)
7			xmin - interval begin
8			xmax - interval end
9			f(x) - integrand
10			integ - integral value
11			}
12			{ Defaults }
13			order = 24 ;
14			xmin = -1 ;
15			xmax = 1 ;
16
17			integ = q1(ceil(order/2-.25)) * (xmax-xmin)/2 ;
18
19			q1(i) = if( i-1.5, wght(i)*(f1(absc(i))+f1(-absc(i))) + q1(i-1),
20			wght(1)*(f1(absc(1))+if(even,f1(-absc(1)),0)) ) ;
21
22			even = floor(order/2+.25) - order/2 + .25 ;
23
24			f1(u) = f((xmax(1+u)+xmin(1-u))/2) ;
25
26			{ Abscissas }
27			absc(i) = select( order,
28			0,
29			.577350269189626, { n = 2 }
30			select( i, { n = 3 }
31			0,
32			.774596669241483
33			),
34			0, 0,
35			select( i, { n = 6 }
36			.238619186083197,
37			.661209386466265,
38			.932469514203152
39			),
40			0, 0, 0, 0, 0,
41			select( i, { n = 12 }
42			.125233408511469,
43			.367831498998180,
44			.587317954286617,
45			.769902674194305,
46			.904117256370475,
47			.981560634246719
48			),
49			0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
50			select( i, { n = 24 }
51			.064056892862606,
52			.191118867473616,
53			.315042679696163,
54			.433793507626045,
55			.545421471388840,
56			.648093651936976,
57			.740124191578554,
58			.820001985973903,
59			.886415527004401,
60			.938274552002733,
61			.974728555971309,
62			.995187219997021
63			)
64			) ;
65
66			{ Weights }
67			wght(i) = select( order,
68			0,
69			1, { n = 2 }
70			select( i, { n = 3 }
71			.888888888888888,
72			.555555555555555
73			),
74			0, 0,
75			select( i, { n = 6 }
76			.467913934572691,
77			.360761573048139,
78			.171324492379170
79			),
80			0, 0, 0, 0, 0,
81			select( i, { n = 12 }
82			.249147045813403,
83			.233492536538355,
84			.203167426723066,
85			.160078328543346,
86			.106939325995318,
87			.047175336386512
88			),
89			0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
90			select( i, { n = 24 }
91			.127938195346752,
92			.125837456346828,
93			.121670472927803,
94			.115505668053726,
95			.107444270115966,
96			.097618652104114,
97			.086190161531953,
98			.073346481411080,
99			.059298584915437,
100			.044277438817420,
101			.028531388628934,
102			.012341229799987
103			)
104			) ;