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