C	ARTURO QUIRANTES SIERRA
C	Department of Applied Physics, Faculty of Sciences
C	University of Granada, 18071 Granada (SPAIN)
C	http://www.ugr.es/local/aquiran/codigos.htm
C	aquiran@ugr.es
C
C	Last update: 10 November 2.003	
C
C	Subroutine BURNS
C	to calculate weights and abscissas in Gauss integration.
C	Reference: "Light Scattering by Particles: Computational Methods",
C	P.W. Barber y S.C. Hill, pp. 172-173
C
C	You need to add a "CALL BURNS(NGAUSS,RX,RW)" in your main routine
C	Input:
C		NGAUSS (number of Gauss quadrature points)
C	Output:
C		RX(k),RW(k): Abscissas and weights; k(=1...NGAUSS)
	SUBROUTINE BURNS(NGAUSS,RX,RW)
	DOUBLE PRECISION RX(200),RW(200)
	DOUBLE PRECISION CN,CNN1,CON1,CON2,XI,X,PM1,PM2,RN,AUX,DER1P
	DOUBLE PRECISION DER2P,APPFCT,B,BISQ,BFROOT,RATIO,PROD
	DOUBLE PRECISION CONST,TOL,C1,C2,C3,C4,PI,P,PM1,PP
	INTEGER NDIV2,NP1,NM1,NM2,IN,NGAUSS,K
	DATA PI/3.141592653589793238462643D0/
	DATA CONST/0.148678816357D0/
	DATA TOL/1.0D-15/
	DATA C1,C2/0.125D0,-0.0807291666D0/
	DATA C3,C4/0.2460286458D0,-1.824438767D0/
C	Warning: if TOL is lower than your compuer's roundoff error,
C	your code might fall into an endless loop!
	IF(NGAUSS.EQ.1) THEN
		RX(1)=0.577350269189626D0
		RW(1)=1.0D0
		RETURN
	END IF
	CN=DFLOAT(NGAUSS)
	NDIV2=NGAUSS/2
	NP1=NGAUSS+1
	CNN1=CN*(CN+1.0D0)
	APPFCT=1.0D0/SQRT((CN+0.5D0)**2+CONST)
C	CON1, CON2 are the integration limits
C	For the usual Gauss quadrature, CON1=1, CON2=0
C	For zero-pi integration, CON1=CON2=pi/2
	CON1=0.5D0*PI
	CON2=0.5D0*PI
	DO 1030 K=1,NDIV2
		B=(DFLOAT(K)-0.25D0)*PI
		BISQ=1.0D0/(B*B)
C	Bisq is a first approximation to Rx
		BFROOT=B*(1.0D0+BISQ*(C1+BISQ*(C2+BISQ*(C3+C4*BISQ))))
	XI=COS(APPFCT*BFROOT)
1010		X=XI
		PM2=1.0D0
		PM1=X
		DO 1020 IN=2,NGAUSS
		    RN=DFLOAT(IN)
			P=((2.0D0*RN-1.0D0)*X*PM1-(RN-1.0D0)*PM2)/RN
			PM2=PM1
			PM1=P
1020		CONTINUE
		PM1=PM2
		AUX=1.0D0/(1.0D0-X*X)
		DER1P=CN*(PM1-X*P)*AUX
		DER2P=(2.0D0*X*DER1P-CNN1*P)*AUX
		RATIO=P/DER1P
		XI=X-RATIO*(1.0D0+RATIO*DER2P/(2.0D0*DER1P))
		IF(ABS(XI-X).GT.TOL) GO TO 1010
C	Now a bit of modification here to calculate Pn with higher accuracy
C	The above calculates PM1=Pn-1(X), not Pn-1(XI)
C	If X and XI differ, it might enlarge roundoff errors
C	So just to play in the safe side...
		X=XI
		PM2=1.0D0
		PM1=X
		PMP1=1.0D0
		DO IN=2,NGAUSS
			RN=DFLOAT(IN)
			P=((2.0D0*RN-1.0D0)*X*PM1-(RN-1.0D0)*PM2)/RN
			PM2=PM1
			PM1=P
			PP=X*PMP1+RN*PM2
			PMP1=PP
		END DO
C	End of the fine-tuning modification
C	The weights are usually calculated as wi=2*(1-xi^2)/(n*Pn-1)^2
C	However, for small xi values, the relative error of the weight is higher
C	Another way to calculate it is wi=2/[(1-xi^2)(Pn')^2]
C	This works better at the endpoints
		RX(K)=-XI
		RW(K)=2.0d0/((1.0D0-XI*XI)*PMP1*PMP1)
		RX(NP1-K)=-RX(K)
		RW(NP1-K)=RW(K)
1030		CONTINUE
		IF(MOD(NGAUSS,2).NE.0) THEN
		RX(NDIV2+1)=0.0D0
		NM1=NGAUSS-1
		NM2=NGAUSS-2
		PROD=CN
		DO 1040 K=1,NM2,2
			PROD=PROD*DFLOAT(NM1-K)/DFLOAT(N-K)
1040		CONTINUE
		RW(NDIV2+1)=2.0D0/(PROD*PROD)
		END IF
		DO 1050 K=1,NGAUSS
		RX(K)=CON1*RX(K)+CON2
		RW(K)=CON1*RW(K)
1050	CONTINUE
	RETURN
	END