function [thzeros, accur, wi, swi]=glgrid(n)
% ------------------------------------------------------------------------|
% use: [thzeros, accur, wi, swi]=glgrid(n)
% DESCRIPTION
% - "glgrid" determines the roots of Legendre polynomials for the 
%   Legendre-Gauss quadrature, i.e. points of Gaussian grid (L-G grid)
% - it is based on dividing the interval theta=<0,180> (what corresponds
%   with an abscissas=<-1,1>) by delta(theta)=180/(2n+1) bringing a very 
%   nice feature of Pn(cos(theta)) ~ sawtooth pattern, which is easy to
%   programme by iteration
% - to speed up an execution time, the function uses a symmetry with
%   respect to the equator
% - it is desired using a function for vector evaluation of the Legendre
%   polynomials, here "pl0" is desired and used
% - the fixed value of a number of interations is used, because the
%   convergence is quite stable
% ------------------------------------------------------------------------|
% INPUT:
% n - degree of Legendre polynomials
% ------------------------------------------------------------------------|
% OUTPUT:
% thzeros - vector(n) of polar angle of the roots, thzeros=(0,180), [deg]
% accur - vector(n) of the achieved accuracy given by Pn(thzeros)), what 
%         should be exactly a zero
% wi -  vector(n) of weights of the roots, gauss-legendre weights
% swi - scalar value for the sum of weights, it should hold sum(w_i)=2
% ------------------------------------------------------------------------|
% OTHERS:
% JS(11-2009), DGFI
% ------------------------------------------------------------------------|
% REFERENCES
% [1] - http://mathworld.wolfram.com/Legendre-GaussQuadrature.html
% ------------------------------------------------------------------------|

% tic;
ro=180/pi;
dth=180/(2*n+1);
th=(0:dth:90+2*dth);     % +2*dth is for symmetry algorithm

% ==============================================
% LEGENDRE POLYNOMIALS FOR ALL POLAR ANGLES

p=(pl0(n,th))';

% ==============================================
% MAIN ITERATION

if mod(n,2)==0 % n=even
 	nc=n;
elseif mod(n,2)==1 % n=odd
	nc=n+1;
end

thzeros=zeros(n,1);
wi=zeros(n,1);
accur=zeros(n,1);

for i=2:2:nc
		
	th1=th(i);
	th2=th(i+1);  % it holds th1 < th2

	% high precision is guaranteed even with the fixed ii
	for ii=1:7     % the first iteration step
		
		da=(th2-th1)/100; % the second iteration step
		th11=(th1:da:th2);
		pnew=pl0(n,th11);  
		
		if pnew(1)>0 && pnew(end)<0  % order of points with respect to sign of 2 consecutive points
			
			j1=find(pnew>0);
			j2=find(pnew<0);
			p1=pnew(j1(end));
			p2=pnew(j2(1));
			th1=th11(j1(end));
			th2=th11(j2(1));
			
		elseif pnew(1)<0 && pnew(end)>0
			
			j1=find(pnew<0);
			j2=find(pnew>0);
			p1=pnew(j1(end));
			p2=pnew(j2(1));
			th1=th11(j1(end));
			th2=th11(j2(1));
		else
			error('pnew has to contain (+) and (-) values, otherwise all is bad, CHECK dth!');
		end
	end
	
	if abs(p1)<=abs(p2)
		th00=th1;
		epsilon=abs(p1);
	else
		th00=th2;
		epsilon=abs(p2);
	end

	% Thetas, weights, accuracy
	pnp1=(pl0(n+1,th00))^2;
	thzeros(i/2)=th00;  
	wi(i/2)=2*((sin(th00/ro))^2)/(pnp1*(n+1)^2);
	accur(i/2)=epsilon;

	% applying the symmetry 
	thzeros(end-i/2+1)=180-th00;
	accur(end-i/2+1)=epsilon;
	wi(end-i/2+1)=wi(i/2);
		
end
% toc
swi=sum(wi);


% ==============================================
% ALL BELLOW CAN BE (UN)COMMENTED - for testing purpose
% ==============================================
% approximate estimate of theta for Pn=0

% % approximate formulae
% m=(1:n);
% thap=(4*m-1)./(4*n+2)*pi+1/(8*n*n)*cot((4*m-1)./(4*n+2));
% pn0ap=pl0(n,thap*ro);

%  PLOTTING
% figure;
% subplot(2,1,1);
% hold on;
% plot(th,p,'b.');
% plot(thap*ro,zeros(n,1),'ro');
% plot(thzeros(1:end),zeros(length(thap),1),'k.');
% plot(thzeros(1:end),wi,'k+');
% xlabel('Polar angle = acos(abscissas), [deg]');
% ylabel('Legendre polynomial');
% 
% legend('Pn','approximate points: Pn=0 - explicit formulae','precise points: Pn0=0 - iteration','Weights','Location','North')
% plot(th,p,'g');
% 
% subplot(2,1,2);
% hold off;
% semilogy(ro*thap,abs(pn0ap),'r');
% hold on;
% semilogy(thzeros(1:end),abs(accur),'k');
% legend('Approximate formulae','Iteration');
% xlabel('Polar angle = acos(abscissas), [deg]');
% ylabel('Accuracy: log10|Pn0(th)-0|');
% text(90,mean(accur),num2str(sum(wi),'sum(w_i)= %18.16f'));

