SUBROUTINE SETUP(NGUIDE,GCON) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (MAXGUIDE=7,MAXPTS=65) DIMENSION GCON(MAXGUIDE) COMMON /PTS/ SPTS(MAXPTS),SVAL(MAXPTS),SINT(MAXPTS),SDER(MAXPTS) EXTERNAL POLY C C SET UP A BIASED AS RANDOM GUIDING FUNCTION AS EXPONENTIAL OF C A SPLINE -- USE BLI TO FIND WHERE TO EVALUATE IT -- C SPTS(I) CONTAINS THE X VALUES OF THE INTERPOLATING SPLINE C SVAL(I) CONTAINS THE Y VALUES OF THE INTERPOLATING SPLINE C SINT(I) CONTAINS THE INTEGRAL [0,X(I)] OF SVAL(I) C SDER(I) CONTAINS THE DERIVATIVE OF SVAL(I) C MAXPTS IS THE NUMBER OF POINTS IN THE INTERPOLATION C A=0 B=64 FTEST=BLI(A,B,POLY,SPTS,SVAL,SINT,BLIER,MAXPTS) WRITE(6,'('' IN SETUP A='',G10.3,'' B='',G10.3,'' NGUIDE='',I5/ 2 '' GCON=''/(4G20.6))') A,B,NGUIDE,(GCON(I),I=1,NGUIDE) C DO 20 I=1,MAXPTS C WRITE(6,'(2G20.6)') SPTS(I),SVAL(I) C20 CONTINUE WRITE(6,*)' FTEST=',FTEST,'+-',BLIER SINT(1)=0.0 DO 30 I=2,MAXPTS IM=I-1 SINT(I)=.5*(SPTS(I)-SPTS(IM))*(SVAL(I)+SVAL(IM))+SINT(IM) SDER(IM)=(SVAL(I)-SVAL(IM))/(SPTS(I)-SPTS(IM)) 30 CONTINUE WRITE(6,*)' ST LINE INTEGRAL ',SINT(MAXPTS-1) RETURN END C FUNCTION BLI(A,B,FUN,XI,FI,ERR,BLIER,NMAX) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION XI(NMAX),FI(NMAX),ERR(NMAX) C C *** USE HALF THE POINTS IN A MIDPOINT MESH C NDO=NMAX/4 H=(B-A)/(NDO-1) DO I=1,NDO XI(I)=A+H*(I-1) FI(I)=FUN(XI(I)) ENDDO AIC=0 DO I=2,NDO AIC=AIC+.5D0*H*(FI(I)+FI(I-1)) ENDDO AIN=AIC H3=H**3 C *** ERR(1) IS THE ERROR IN THE INTERVAL A TO X(1), FPP USES PTS 1,2,3 FPL=(FI(2)-FI(1))/(XI(2)-XI(1)) FPU=(FI(3)-FI(2))/(XI(3)-XI(2)) FPP=2*(FPU-FPL)/(XI(3)-XI(1)) C *** ERR(1) IS THE ERROR IN THE INTERVAL X(1) TO X(2), USES SAME FPP AS ERR(1)=ABS(8.333333D-2*H3*FPP) C *** MID REGION USES I+1 TO I+2 FOR FPU AND I-1 TO I FOR FPL FPTEMP=FPU DO I=2,NDO-2 FPU=(FI(I+2)-FI(I+1))/(XI(I+2)-XI(I+1)) FPP=2*(FPU-FPL)/(XI(I+2)+XI(I+1)-XI(I)-XI(I-1)) ERR(I)=ABS(8.333333D-2*H3*FPP) FPL=FPTEMP FPTEMP=FPU ENDDO C *** INTERVAL FROM XI(N-1) TO XI(N) USES FPU AND FPTEMP TO ESTIMATE FPP FPP=2*(FPU-FPL)/(XI(NDO)-XI(NDO-2)) ERR(NDO-1)=ABS(8.333333D-2*H3*FPP) ERR(NDO)=0 C *** INITIALIZING THE CURVE FIT PARTS SMA=0 SMB=0 SMD=0 B0=0 B1=0 SSQ=0 C *** FINDING THE LARGEST ERROR AND THE TOTAL SUM OF THE ERRORS 50 ERRM=ERR(1) ILE=1 ERRS=ERRM DO I=2,NDO-1 ERRS=ERRS+ERR(I) IF(ERR(I).GT.ERRM)THEN ERRM=ERR(I) ILE=I ENDIF ENDDO C *** UPDATING THE FITTING PARTS ERRS=DMAX1(ERRS,1.D-14*AIN) IF(ERRS.GT.0.D0)THEN ERRI=1/ERRS ERRI2=ERRI*ERRI SMA=SMA+ERRI2 SMB=SMB+ERRI SMD=SMD+1 AINMAC=AIN-AIC B0=B0+AINMAC*ERRI2 B1=B1+AINMAC*ERRI SSQ=SSQ+AINMAC*AINMAC*ERRI2 ENDIF IF(NDO.EQ.NMAX)THEN IF(SMD.EQ.0.D0)THEN BLI=0 BLIER=0 RETURN ENDIF DEN=SMA*SMD-SMB*SMB IF(DEN.LT.1.D-14*SMA*SMD)THEN BLI=B0/SMA+AIC BLIER=1.D-12*BLI RETURN ENDIF C0=(SMD*B0-SMB*B1)/DEN C1=(-SMB*B0+SMA*B1)/DEN CHIS=DMAX1(1.D-14*SSQ,(SSQ-C0*B0-C1*B1)) BLI=C0+AIC BLIER=DSQRT(CHIS*SMD/DEN) RETURN ENDIF C *** POINT TO BE INSERTED WILL BE ILE+1, THUS CURRENT XI(ILE+1) BECOMES C *** XI(ILE+2) DO I=NDO,ILE+1,-1 XI(I+1)=XI(I) FI(I+1)=FI(I) ERR(I+1)=ERR(I) ENDDO NDO=NDO+1 IADD=ILE+1 XI(IADD)=.5D0*(XI(ILE)+XI(ILE+2)) FI(IADD)=FUN(XI(IADD)) AIN=AIN-.5D0*(XI(ILE+2)-XI(ILE))*(FI(ILE)+FI(ILE+2)) 2 +.5D0*(XI(ILE+1)-XI(ILE))*(FI(ILE)+FI(ILE+1)) 3 +.5D0*(XI(ILE+2)-XI(ILE+1))*(FI(ILE+1)+FI(ILE+2)) IF(ILE.EQ.1)THEN FPL=(FI(2)-FI(1))/(XI(2)-XI(1)) FPU=(FI(3)-FI(2))/(XI(3)-XI(2)) FPP=2*(FPU-FPL)/(XI(3)-XI(1)) C *** ERR(1) IS THE ERROR IN THE INTERVAL X(1) TO X(2) ERR(1)=ABS(8.333333D-2*(XI(2)-XI(1))**3*FPP) C *** ERR3 IS BASED ON 12 FOR FPL AND 34 FOR FPU FTEMP=FPU FPU=(FI(4)-FI(3))/(XI(4)-XI(3)) FPP=2*(FPU-FPL)/(XI(4)+XI(3)-XI(2)-XI(1)) ERR(2)=ABS(8.333333D-2*(XI(3)-XI(2))**3*FPP) IF(NDO.GE.5)THEN FPL=FTEMP FPU=(FI(5)-FI(4))/(XI(5)-XI(4)) FPP=2*(FPU-FPL)/(XI(4)+XI(4)-XI(3)-XI(2)) ENDIF ERR(3)=ABS(8.333333D-2*(XI(4)-XI(3))**3*FPP) GOTO 50 ENDIF IF(ILE.EQ.NDO-1)THEN FPU=(FI(NDO)-FI(NDO-1))/(XI(NDO)-XI(NDO-1)) FPL=(FI(NDO-1)-FI(NDO-2))/(XI(NDO-1)-XI(NDO-2)) FTEMP=FPL FPP=2*(FPU-FPL)/(XI(NDO)-XI(NDO-2)) ERR(NDO-1)=ABS(8.333333D-2*FPP*(XI(NDO)-XI(NDO-1))**3) FPL=(FI(NDO-2)-FI(NDO-3))/(XI(NDO-2)-XI(NDO-3)) FPP=2*(FPU-FPL)/(XI(NDO)+XI(NDO-1)-XI(NDO-2)-XI(NDO-3)) ERR(NDO-2)=ABS(8.333333D-2*(XI(NDO-1)-XI(NDO-2))**3*FPP) IF(NDO-4.GT.0)THEN FPL=(FI(NDO-3)-FI(NDO-4))/(XI(NDO-3)-XI(NDO-4)) FPU=FTEMP FPP=2*(FPU-FPL)/(XI(NDO-1)+XI(NDO-2)-XI(NDO-3)-XI(NDO-4)) ENDIF ERR(NDO-3)=ABS(8.333333D-2*(XI(NDO-2)-XI(NDO-3))**3*FPP) GOTO 50 ENDIF C *** NEED TO RECALCULATE ERR(ILE-2) TO ERR(ILE+1) IF(ILE.GT.4)THEN FPL=(FI(ILE-2)-FI(ILE-3))/(XI(ILE-2)-XI(ILE-3)) FPU=(FI(ILE)-FI(ILE-1))/(XI(ILE)-XI(ILE-1)) FPP=2*(FPU-FPL)/(XI(ILE)+XI(ILE-1)-XI(ILE-2)-XI(ILE-3)) FPL=(FI(ILE-1)-FI(ILE-2))/(XI(ILE-1)-XI(ILE-2)) ELSE ILE=MAX0(3,ILE) FPL=(FI(ILE-1)-FI(ILE-2))/(XI(ILE-1)-XI(ILE-2)) FPU=(FI(ILE)-FI(ILE-1))/(XI(ILE)-XI(ILE-1)) FPP=2*(FPU-FPL)/(XI(ILE)-XI(ILE-2)) ENDIF ERR(ILE-2)=ABS(8.333333D-2*(XI(ILE-1)-XI(ILE-2))**3*FPP) FTEMP=FPU FPU=(FI(ILE+1)-FI(ILE))/(XI(ILE+1)-XI(ILE)) FPP=2*(FPU-FPL)/(XI(ILE+1)+XI(ILE)-XI(ILE-1)-XI(ILE-2)) ERR(ILE-1)=ABS(8.333333D-2*(XI(ILE)-XI(ILE-1))**3*FPP) FPL=FTEMP FTEMP=FPU FPU=(FI(ILE+2)-FI(ILE+1))/(XI(ILE+2)-XI(ILE+1)) FPP=2*(FPU-FPL)/(XI(ILE+2)+XI(ILE+1)-XI(ILE)-XI(ILE-1)) ERR(ILE)=ABS(8.333333D-2*(XI(ILE+1)-XI(ILE))**3*FPP) FPL=FTEMP IF(ILE+2.EQ.NDO)THEN FPP=2*(FPU-FPL)/(XI(ILE+2)-XI(ILE)) ELSE FPU=(FI(ILE+3)-FI(ILE+2))/(XI(ILE+3)-XI(ILE+2)) FPP=2*(FPU-FPL)/(XI(ILE+3)+XI(ILE+2)-XI(ILE+1)-XI(ILE)) ENDIF ERR(ILE+1)=ABS(8.333333D-2*(XI(ILE+2)-XI(ILE+1))**3*FPP) GOTO 50 END C FUNCTION POLY(X) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MAXGUIDE=7) COMMON /GUIDE/ GCON(MAXGUIDE),NGUIDE C *** The first coefficient is an exponential c(i)*x C *** The rest are in the form c(i)*exp((x-c(i+1))/c(i+2))**2 POLY=GCON(1)*X NEXP=(NGUIDE-1)/3 IOFF=1 DO I=1,NEXP POLY=POLY+GCON(IOFF+1)*EXP(-((X-GCON(IOFF+2))/GCON(IOFF+3))**2) IOFF=IOFF+3 ENDDO POLY=EXP(POLY) RETURN END C SUBROUTINE LOCATE(ARR,M,X,J) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION ARR(1) C C FIND THE LAST ELEMENT IN ARR(J) S.T. ARR(J)