% Function for initializing the matrices 'tt' and 'w' for Gaussian
% Integration; with n = 5, 10, 20, 24, 32, 40 or 48
% SAMPLE CODE FOR IMPLEMENTING DOUBLE INTEGRATION
% gaussian_points = 10 ; % Define the number of points to be used
%
% a = 1; % Lower integration limit for y
% b = 2; % Upper integration limit for y
%
% c = 3; % Lower integration limit for x
% d = 4; % Upper integration limit for x
%
% sum = 0;
%
% for ii = 1:gaussian_points
%
% yy = ((d - c)*tt(ii) + d + c)/2;
%
% inner_sum = 0;
%
% for jj = 1:gaussian_points
% xx = ((b - a)*tt(jj) + b + a)/2;
% Fm1 = xx + yy*yy; % I N T E G R A N D
% inner_sum = inner_sum + 0.5*(b - a)*Fm1*w(jj);
% end
%
% sum = sum + 0.5*(d - c)*inner_sum*w(ii);
%
% end
% -------------------------------------------------------------------------
function [tt w] = func_gaussian_integration_weights (gaussian_points)
switch (gaussian_points)
case 5
tt = [-0.906179845938664
-0.538469310105683
0.0
0.538469310105683
0.906179845938664]; % Points where function will be evaluated
w = [0.236926885056189
0.478628670499366
0.568888888888889
0.478628670499366
0.236926885056189]; % Constants to be multiplied
case 10
tt = [-0.973906528517172
-0.865063366688985
-0.679409568299024
-0.433395394129247
-0.148874338981631
0.148874338981631
0.433395394129247
0.679409568299024
0.865063366688985
0.973906528517172]; % Points where function will be evaluated
w = [0.066671344308688
0.149451349150581
0.219086362515982
0.269266719309996
0.295524224714753
0.295524224714753
0.269266719309996
0.219086362515982
0.149451349150581
0.066671344308688]; % Constants to be multiplied
case 20
tt = [-0.993128599185094924786
-0.963971927277913791268
-0.912234428251325905868
-0.839116971822218823395
-0.746331906460150792614
-0.636053680726515025453
-0.510867001950827098004
-0.373706088715419560673
-0.227785851141645078080
-0.076526521133497333755
0.076526521133497333755
0.227785851141645078080
0.373706088715419560673
0.510867001950827098004
0.636053680726515025453
0.746331906460150792614
0.839116971822218823395
0.912234428251325905868
0.963971927277913791268
0.993128599185094924786]; % Points where function will be evaluated
w = [0.017614007139152118312
0.040601429800386941331
0.062672048334109063570
0.083276741576704748725
0.101930119817240435037
0.118194531961518417312
0.131688638449176626898
0.142096109318382051329
0.149172986472603746788
0.152753387130725850698
0.152753387130725850698
0.149172986472603746788
0.142096109318382051329
0.131688638449176626898
0.118194531961518417312
0.101930119817240435037
0.083276741576704748725
0.062672048334109063570
0.040601429800386941331
0.017614007139152118312]; % Constants to be multiplied
case 24
tt = zeros (24, 1);
tt(13:24) = [0.064056892862605626085
0.191118867473616309159
0.315042679696163374387
0.433793507626045138487
0.545421471388839535658
0.648093651936975569252
0.740124191578554364244
0.820001985973902921954
0.866415527004401034213
0.938274552002732758524
0.974728555971309498198
0.995187219997021360180];
tt(1:12) = -1*flipud (tt (13:24)); % Points where function will be evaluated
w = zeros (24, 1);
w(13:24) = [0.127938195346752156974
0.125837456346828296121
0.121670472927803391204
0.115505668053725601353
0.107444270115965634783
0.097618652104113888270
0.086190161531953275917
0.073346481411080305734
0.059298584915436780746
0.044277438817419806169
0.028531388628933663181
0.012341229799987199547];
w(1:12) = flipud (w (13:24));
case 32
tt = zeros (32, 1);
tt(17:32) = [0.048307665687738316235
0.144471961582796493485
0.239287362252137074545
0.331868602282127649780
0.421351276130635345364
0.506899908932229390024
0.587715757240762329041
0.663044266930215200975
0.732182118740289680387
0.794483795967942406963
0.849367613732569970134
0.896321155766052123965
0.934906075937739689171
0.964762255587506430774
0.985611511545268335400
0.997263861849481563545];
tt(1:16) = -1*flipud (tt (17:32)); % Points where function will be evaluated
w = zeros (32, 1);
w(17:32) = [0.096540088514727800567
0.095638720079274859419
0.093844399080804565639
0.091173878695763884713
0.087652093004403811143
0.083311924226946755222
0.078193895787070306472
0.072345794108848506225
0.065822222776361846838
0.058684093478535547145
0.050998059262376176196
0.042835898022226680657
0.034273862913021433103
0.025392065309262059456
0.016274394730905670605
0.007018610009470096600];
w(1:16) = flipud (w (17:32));
case 40
tt = zeros (40, 1);
tt(21:40) = [0.038772417506050821933
0.116084070675255208483
0.192697580701371099716
0.268152185007253681141
0.341994090825758473007
0.413779204371605001525
0.483075801686178712909
0.549467125095128202076
0.612553889667980237953
0.671956684614179548379
0.727318255189927103281
0.778305651426519387965
0.824612230833311663196
0.865959503212259503281
0.902098806968874296728
0.932812808278676533361
0.957916819213791655805
0.977259949983774262663
0.990726238699457006453
0.998237709710559200350];
tt(1:20) = -1*flipud (tt (21:40)); % Points where function will be evaluated
w = zeros (40, 1);
w(21:40) = [0.077505947978424811264
0.077039818164247965588
0.076110361900626242372
0.074723169057968264200
0.072886582395804059061
0.070611647391286779695
0.067912045815233903826
0.064804013456601038075
0.061306242492928939167
0.057439769099391551367
0.053227846983936824355
0.048695807635072232061
0.043870908185673271992
0.038782167974472017640
0.033460195282547847393
0.027937006980023401098
0.022245849194166957262
0.016421058381907888713
0.010498284531152813615
0.004521277098533191258];
w(1:20) = flipud (w (21:40));
case 48
tt = zeros (48, 1);
tt(25:48) = [0.032380170962869362033
0.097004699209462698930
0.161222356068891718056
0.224763790394689061225
0.287362487355455576736
0.348755886292160738160
0.408686481990716729916
0.466902904750958404545
0.523160974722233033678
0.577224726083972703818
0.628867396776513623995
0.677872379632663905212
0.724034130923814654674
0.767159032515740339254
0.807066204029442627083
0.843588261624393530711
0.876572020274247885906
0.905879136715569672822
0.931386690706554333114
0.952987703160430860723
0.970591592546247250461
0.984124583722826857745
0.993530172266350757548
0.998771007252426118601];
tt(1:24) = -1*flipud (tt (25:48)); % Points where function will be evaluated
w = zeros (48, 1);
w(25:48) = [0.064737696812683922503
0.064466164435950082207
0.063924238584648186624
0.063114192286254025657
0.062039423159892663904
0.06070443916589388053
0.059114839698395635746
0.057277292100403215705
0.055199503699984162868
0.052890189485193667096
0.050359035553854474958
0.047616658492490474826
0.044674560856694280419
0.041545082943464749214
0.038241351065830706317
0.034777222564770438893
0.031167227832798088902
0.027426509708356948200
0.023570760839324379141
0.019616160457355527814
0.015579315722943848728
0.011477234579234539490
0.007327553901276262102
0.003153346052305838633];
w(1:24) = flipud (w (25:48));
otherwise
error('Parameter "gaussian_points" can only take values 5, 10, 20, 24, 32, 40 or 48');
end
end
