% function [r, w] = gauss(m)
% 
% Get the double Gauss weights and nodes (2m points) for the discrete 
% ordinates solution of the transport equation. As currently implemented
% the only legal values of m are
%
% m=8, 20, 40, 100 
%
function [r, w] = gauss(m)
%
% legal angular dimension?
%
if m ~=4 & m ~= 8 & m ~= 40 & m ~= 20 & m ~= 100
    disp(' dimension error ')
    return
end
%
% single 4 point Gauss on [-1,1]
%
if m==4
       ri(4)=.8611363115940526;
       ri(3)=.3399810435848563;
       wi(4)=.3478548451374539;
       wi(3)=.6521451548625461;
       ri(2)=-ri(3);
       ri(1)=-ri(4);
       wi(1)=wi(4);
       wi(2)=wi(3);
	   for i=1:4
	       r(i)=ri(i);
	       w(i)=wi(i);
	   end
end
%
% double Gauss on [-1,1]
%
if m == 8
       ri(4)=.8611363115940526;
       ri(3)=.3399810435848563;
       wi(4)=.3478548451374539;
       wi(3)=.6521451548625461;
       ri(2)=-ri(3);
       ri(1)=-ri(4);
       wi(1)=wi(4);
       wi(2)=wi(3);
       mm=m/2;
       for i=1:mm
       r(i+mm)=(1 + ri(i))*.5;
       w(i+mm)=wi(i)*.5;
       r(i)=(ri(i)-1)*.5;
       w(i)=wi(i)*.5;
       end
end
if m == 20
       ri(6) =               0.14887433898163121571;
       wi(6) =               0.29552422471475287002;
       ri(7) =               0.43339539412924721340;
       wi(7) =               0.26926671930999634963;
       ri(8) =               0.67940956829902443559;
       wi(8) =               0.21908636251598204159;
       ri(9) =               0.86506336668898453635;
       wi(9) =               0.14945134915058058689;
       ri(10) =              0.97390652851717174343;
       wi(10) =              0.066671344308688137992;
       for i=1:10,
          ri(i)=-ri(11-i);
          wi(i)=wi(11-i);
       end
       mm=m/2;
       for i=1:mm
       r(i+mm)=(1 + ri(i))*.5;
       w(i+mm)=wi(i)*.5;
       r(i)=(ri(i)-1)*.5;
       w(i)=wi(i)*.5;
       end
end
if m == 40
       ri(20)=.993128599185095;
       ri(19)=.963971927277914;
       ri(18)=.912234428251326;
       ri(17)=.839116971822218;
       ri(16)=.746331906460151;
       ri(15)=.636053680726515;
       ri(14)=.510867001950827;
       ri(13)=.373706088715420;
       ri(12)=.227785851141645;
       ri(11)=.076526521133497;
       wi(20)=.017614007139152;
       wi(19)=.040601429800387;
       wi(18)=.062672048334109;
       wi(17)=.083276741576705;
       wi(16)=.101930119817240;
       wi(15)=.118194531961518;
       wi(14)=.131688638449177;
       wi(13)=.142096109318382;
       wi(12)=.149172986472604;
       wi(11)=.152753387130726;
       for i=1:10,
          ri(i)=-ri(21-i);
          wi(i)=wi(21-i);
      end
      mm=m/2;
      for i=1:mm
       r(i+mm)=(1 + ri(i))*.5;
       w(i+mm)=wi(i)*.5;
       r(i)=(ri(i)-1)*.5;
       w(i)=wi(i)*.5;
       end
end
if m == 100
       ri(26) = .31098338327188876362D-1;
       wi(26) = .62176616655347259621D-1;
       ri(27) = .93174701560086142793D-1;
       wi(27) = .61936067420683242613D-1;
       ri(28) = .15489058999814589446;
       wi(28) = .61455899590316664938D-1;
       ri(29) = .21600723687604175827;
       wi(29) = .60737970841770218267D-1;
       ri(30) = .27628819377953200975;
       wi(30) = .59785058704265460483D-1;
       ri(31) = .33550024541943734846;
       wi(31) = .58600849813222444240D-1;
       ri(32) = .39341431189756514986;
       wi(32) = .57189925647728380442D-1;
       ri(33) = .44980633497403876842;
       wi(33) = .55557744806212519539D-1;
       ri(34) = .50445814490746421210;
       wi(34) = .53710621888996244733D-1;
       ri(35) = .55715830451465009343;
       wi(35) = .51655703069581136722D-1;
       ri(36) = .60770292718495022566;
       wi(36) = .49400938449466316738D-1;
       ri(37) = .65589646568543935601;
       wi(37) = .46955051303948433516D-1;
       ri(38) = .70155246870682219651;
       wi(38) = .44327504338803273842D-1;
       ri(39) = .74449430222606849394;
       wi(39) = .41528463090147696435D-1;
       ri(40) = .78455583290039931921;
       wi(40) = .38568756612587677679D-1;
       ri(41) = .82158207085933598890;
       wi(41) = .35459835615146151344D-1;
       ri(42) = .85542976942994608525;
       wi(42) = .32213728223578014243D-1;
       ri(43) = .88596797952361305839;
       wi(43) = .28842993580535197101D-1;
       ri(44) = .91307855665579185089;
       wi(44) = .25360673570012391698D-1;
       ri(45) = .93665661894487795003;
       wi(45) = .21780243170124793939D-1;
       ri(46) = .95661095524280792546;
       wi(46) = .18115560713489391870D-1;
       ri(47) = .97286438510669204227;
       wi(47) = .14380822761485573952D-1;
       ri(48) = .98535408404800584048;
       wi(48) = .10590548383650968609D-1;
       ri(49) = .99403196943209071179;
       wi(49) = .67597991957454011541D-2;
       ri(50) = .99886640442007101903;
       wi(50) = .29086225531551410800D-2;
       for i=1:25
           ri(i)=-ri(51-i);
           wi(i)=wi(51-i);
       end
       mm=m/2;
      for i=1:mm
       r(i+mm)=(1 + ri(i))*.5;
       w(i+mm)=wi(i)*.5;
       r(i)=(ri(i)-1)*.5;
       w(i)=wi(i)*.5;
       end
end
r=r';
w=w';
return
end
