! -*- F90 -*- subroutine ACFGPevolvep(xin,qin,p2in,ip2in,pdf) include 'parmsetup.inc' real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5 real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu real*4 par,par2,calc,celc common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32) character*16 name(nmxset) integer nmem(nmxset),ndef(nmxset),mmem common/NAME/name,nmem,ndef,mmem integer nset save if(imem.eq.1) then call ACFGP1(xin,qin,upv,dnv,usea,dsea,str,chm,glu) elseif(imem.eq.2) then call ACFGP2(xin,qin,upv,dnv,usea,dsea,str,chm,glu) elseif(imem.eq.3.or.imem.eq.0) then call SFAFG1(xin,qin,upv,dnv,usea,dsea,str,chm,glu) else CONTINUE endif pdf(-6)= 0.0d0 pdf(6)= 0.0d0 pdf(-5)= 0.0d0 pdf(5 )= 0.0d0 pdf(-4)= chm pdf(4 )= chm pdf(-3)= str pdf(3 )= str pdf(-2)= usea pdf(2 )= upv pdf(-1)= dsea pdf(1 )= dnv pdf(0 )= glu return !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc entry ACFGPread(nset) read(1,*)nmem(nset),ndef(nset) do iset = 1,3 read(1,*)(par(j,iset),j=1,4) read(1,*)(par(j,iset),j=5,8) read(1,*)(par(j,iset),j=9,12) read(1,*)(par(j,iset),j=13,16) read(1,*)(par(j,iset),j=17,20) read(1,*)(par(j,iset),j=21,24) read(1,*)(par(j,iset),j=25,28) read(1,*)(par(j,iset),j=29,30) do j=1,32 do k=1,20 read(1,*)(calc(i,k,j,iset),i=1,4) read(1,*)(calc(i,k,j,iset),i=5,8) enddo enddo enddo ! last one read(1,*)(par2(j),j=1,4) read(1,*)(par2(j),j=5,8) read(1,*)(par2(j),j=9,12) read(1,*)(par2(j),j=13,16) read(1,*)(par2(j),j=17,20) read(1,*)(par2(j),j=21,24) read(1,*)(par2(j),j=25,28) read(1,*)(par2(j),j=29,30) do j=1,32 do k=1,20 read(1,*)(celc(i,k,j),i=1,4) read(1,*)(celc(i,k,j),i=5,8) enddo enddo return ! !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc entry ACFGPalfa(alfas,qalfa) call getnset(iset) call GetOrderAsM(iset,iord) call Getlam4M(iset,imem,qcdl4) call Getlam5M(iset,imem,qcdl5) call aspdflib(alfas,Qalfa,iord,qcdl5) return ! !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc entry ACFGPinit(Eorder,Q2fit) return ! !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc entry ACFGPpdf(mem) imem = mem return ! 1000 format(5e13.5) END !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE ACFGP1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL) ! ! INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIV ! QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF ! ! THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MA ! CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABL ! ! THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS A ! BY A VDM "ANSATZ" AT Q2=.25 GEV**2. ! ! THE PROGRAM WORKS FOR 2. GEV**2 < Q2 <5.5E+5 AND .00137 < X < . ! ! THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME ! ! THE VALUE OF LAMBDA-MSB IS 200 MEV ! ! THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT': ! X*U=X*U(X,Q2) ! X*D= ... ! X*S= ... ! X*C= ... (MASSLESS CHARM WITH C( ! X*CM= ... (MASSIVE CHARM WITH MC=1.5 ! X*GLU=GLUON(X,Q2)*X ! ! ! F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM ! F2C= " " " WITH MASSIVE CHARM ! double precision & & DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL REAL Q2 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32) DIMENSION XPDF(7),CALCO(8,20,32) COMMON/W5051I7/CALCO EXTERNAL AFCPLU DATA ZERO/0.0/ !---------------------------------------------------------------------- DATA ISTART/0/ SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/ ! IF (ISTART.EQ.0) THEN ISTART=1 DO 10 K=1,32 DO 10 I=1,20 DO 10 M=1,8 10 CALCO(M,I,K) = CALC(M,I,K,1) OWLAM=PAR(1,1) OWLAM2=OWLAM**2 Q02=PAR(30,1) FLAV=PAR(25,1) DELTA=PAR(29,1) CALL WATE32 ENDIF ! X = DX Q = DQ Q2 = Q*Q IDQ2=2 SB=0. IF((Q2-Q02).LE.0) THEN GOTO 1 ELSE GOTO 2 ENDIF 2 IF((IDQ2-1).LE.0) THEN GOTO 1 ELSE GOTO 3 ENDIF 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2)) 1 CONTINUE CALL AURGAM(8,0,X,SB,XPDF(7)) CALL AURGAM(7,0,X,SB,SING) CALL AURGAM(4,0,X,SB,DPLUSNS) CALL AURGAM(3,0,X,SB,CPLUSNS) CALL AURGAM(5,0,X,SB,UPLUSNS) CALL AURGAM(6,0,X,SB,SPLUSNS) XPDF(3) = CPLUSNS XPDF(4) = DPLUSNS XPDF(5) = UPLUSNS XPDF(6) = SPLUSNS XPDF(1) = SING ! ADD = XPDF(1)/FLAV UPLUS=XPDF(5)+ADD DPLUS=-XPDF(4)+ADD SPLUS=-XPDF(6)+ADD CPLUS=-XPDF(3)+ADD UB=UPLUS*0.5 UV=UB DB=DPLUS*0.5 DV=DB SB=SPLUS*0.5 CB=CPLUS*0.5 SING=XPDF(1) GLU=XPDF(7) GL=GLU ! DUV=MAX(ZERO,UV) DDV=MAX(ZERO,DV) DUB=MAX(ZERO,UB) DDB=MAX(ZERO,DB) DSB=MAX(ZERO,SB) DCB=MAX(ZERO,CB) DGL=MAX(ZERO,GL) ! RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE ACFGP2(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL) ! ! INTERPOLATION PROGRAM WHICH INTERPOLATES THE GRID "DATAGA" AND GIV ! QUARK AND GLUON DISTRIBUTIONS IN THE REAL PHOTON, AS FUNCTIONS OF ! ! THE Q2-EVOLUTION IS PERFORMED WITH BLL AP-EQUATIONS AND NF=4. A MA ! CHARM DISTRIBUTION (BORROWED FROM GLUCK AND REYA) IS ALSO AVAILABL ! ! THE BOUNDARY CONDITIONS ARE SUCH THAT THE DISTRIBUTION FUNCTIONS A ! BY A VDM "ANSATZ" AT Q2=.25 GEV**2. ! ! THE PROGRAM WORKS FOR 2. GEV**2 < Q2 <5.5E+5 AND .00137 < X < . ! ! THE DISTRIBUTIONS ARE CALCULATED IN THE MSBAR FACTORIZATION SCHEME ! ! THE VALUE OF LAMBDA-MSB IS 200 MEV ! ! THE OUTPUT IS WRITTEN IN THE FILE 'FILEOUT': ! X*U=X*U(X,Q2) ! X*D= ... ! X*S= ... ! X*C= ... (MASSLESS CHARM WITH C( ! X*CM= ... (MASSIVE CHARM WITH MC=1.5 ! X*GLU=GLUON(X,Q2)*X ! ! ! F2 = PHOTON STRUCTURE FUNCTION WITHOUT CHARM ! F2C= " " " WITH MASSIVE CHARM ! double precision & & DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL REAL Q2 common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32) DIMENSION XPDF(7),CALCO(8,20,32) COMMON/W5051I7/CALCO EXTERNAL AFCPLU DATA ZERO/0.0/ !---------------------------------------------------------------------- DATA ISTART/0/ SAVE ISTART,OWLAM2,Q02,FLAV, /W5051I7/ ! IF (ISTART.EQ.0) THEN ISTART=1 DO 10 K=1,32 DO 10 I=1,20 DO 10 M=1,8 10 CALCO(M,I,K) = CALC(M,I,K,2) OWLAM=PAR(1,2) OWLAM2=OWLAM**2 Q02=PAR(30,2) FLAV=PAR(25,2) DELTA=PAR(29,2) CALL WATE32 ENDIF ! X = DX Q = DQ Q2 = Q*Q IDQ2=2 SB=0. IF((Q2-Q02).LE.0) THEN GOTO 1 ELSE GOTO 2 ENDIF 2 IF((IDQ2-1).LE.0) THEN GOTO 1 ELSE GOTO 3 ENDIF 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2)) 1 CONTINUE CALL AURGAM(8,0,X,SB,XPDF(7)) CALL AURGAM(7,0,X,SB,SING) CALL AURGAM(4,0,X,SB,DPLUSNS) CALL AURGAM(3,0,X,SB,CPLUSNS) CALL AURGAM(5,0,X,SB,UPLUSNS) CALL AURGAM(6,0,X,SB,SPLUSNS) XPDF(3) = CPLUSNS XPDF(4) = DPLUSNS XPDF(5) = UPLUSNS XPDF(6) = SPLUSNS XPDF(1) = SING ! ADD = XPDF(1)/FLAV UPLUS=XPDF(5)+ADD DPLUS=-XPDF(4)+ADD SPLUS=-XPDF(6)+ADD CPLUS=-XPDF(3)+ADD UB=UPLUS*0.5 UV=UB DB=DPLUS*0.5 DV=DB SB=SPLUS*0.5 CB=CPLUS*0.5 SING=XPDF(1) GLU=XPDF(7) GL=GLU !... get parton density with massive charm CPLUM=AFCPLU(X,Q2) CB=CPLUM*0.5 ! DUV=MAX(ZERO,UV) DDV=MAX(ZERO,DV) DUB=MAX(ZERO,UB) DDB=MAX(ZERO,DB) DSB=MAX(ZERO,SB) DCB=MAX(ZERO,CB) DGL=MAX(ZERO,GL) ! RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE SFAFG1(DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DGL) ! !*********************************************************************** ! ( 1st of February 1994) ! This is an interpolation program which reads the files GRPOL and ! GRVDM and gives the quark and gluon distributions in real photon ! as functions of x and Q**2. ! ! The Q**2 evolution is a BLL evolution (MSbar scheme) with Nf=4 ! and LAMBDA(MSbar)=.200 Gev. ! ! A massless charm distribution is generated for Q**2 > 2 Gev**2. ! ! The distributions are the sum of a pointlike part (PL) and of a ! Vdm part (VDM): ! dist=PL + KA*VDM ! KA is a factor which can be adjusted ( the default value is KA=1.0 ! The file GRPOL contains the pointlike part of the distributions. ! The file GRVDM contains the vdm part (A precise definition of this ! latter is given in the paper "PARTON DISTRIBUTIONS IN THE PHOTON", ! Preprint LPTHE Orsay 93-37, by P.Aurenche,M.Fontannaz and J.Ph.Gui ! ! The output of the program is written in the file GETOUT with the ! following conventions ! UPLUS=x(u+ubar) ! DPLUS=x(d+dbar) ! SPLUS=x(s+sbar) ! CPLUS=x(c+cbar) ! SING =UPLUS+DPLUS+SPLUS+CPLUS ! GLU =x*g ! ! The interpolation is valid for 2. < Q**2 < 5.5E+5 Gev**2, ! and for .0015< x < .99 ! ! The program also gives the structure function F2: ! F2 = q*Cq + g*Cg + Cgam ! Cq and Cg are the Wilson coeficients and Cgam is the direct term. ! ! Although the charm quark evolution is massless, the direct term ! Cgam includes the effects due to the charm quark mass. The charm ! quark threshold is therefore correctly described at the lowest ! ordre in alphastrong (Details are given in the preprint). ! ! !*********************************************************************** ! double precision & & DX,DQ,DUV,DDV,DUB,DDB,DSB,DCB,DBB,DGL REAL X, Q, UV, DV, UB, DB, SB, CB, BB, GL REAL Q2 DIMENSION XPDF(7) common/acfgp/PAR(30,3),CALC(8,20,32,3),PAR2(30),CELC(8,20,32) DIMENSION CALCO(8,20,32) DIMENSION CELCO(8,20,32) COMMON/W5051IA/CALCO COMMON/W5051IB/CELCO EXTERNAL AFCPLU DATA ZERO/0.0/ !---------------------------------------------------------------------- DATA ISTART/0/ SAVE ISTART,OWLAM2,Q02,FLAV,KA, /W5051IA/, /W5051IB/ ! IF (ISTART.EQ.0) THEN ISTART=1 DO 10 K=1,32 DO 10 I=1,20 DO 10 M=1,8 CALCO(M,I,K) = CALC(M,I,K,3) 10 CELCO(M,I,K) = CELC(M,I,K) OWLAM=PAR(1,3) OWLAM2=OWLAM**2 Q02=PAR(30,3) FLAV=PAR(25,3) DELTA=PAR(29,3) CALL WATE32 KA=1.0 ENDIF ! X = DX Q = DQ Q2 = Q*Q IDQ2=2 SB=0. IF((Q2-Q02).LE.0) THEN GOTO 1 ELSE GOTO 2 ENDIF 2 IF((IDQ2-1).LE.0) THEN GOTO 1 ELSE GOTO 3 ENDIF 3 SB= LOG( LOG( MAX(Q02,Q2)/OWLAM2)/ LOG(Q02/OWLAM2)) 1 CONTINUE CALL AFGINT(8,0,X,SB,XPDF(7)) CALL AFGINT(7,0,X,SB,SING) CALL AFGINT(4,0,X,SB,DPLUSNS) CALL AFGINT(3,0,X,SB,CPLUSNS) CALL AFGINT(5,0,X,SB,UPLUSNS) CALL AFGINT(6,0,X,SB,SPLUSNS) XPDF(3) = CPLUSNS XPDF(4) = DPLUSNS XPDF(5) = UPLUSNS XPDF(6) = SPLUSNS XPDF(1) = SING ! ADD = XPDF(1)/FLAV UPLUS= XPDF(5)+ADD DPLUS=-XPDF(4)+ADD SPLUS=-XPDF(6)+ADD CPLUS=-XPDF(3)+ADD SING=XPDF(1) GLU=XPDF(7) GL=GLU ! CALL AFGIN2(8,0,X,SB,XPDF(7)) CALL AFGIN2(7,0,X,SB,SING) CALL AFGIN2(4,0,X,SB,DPLUSNS) CALL AFGIN2(3,0,X,SB,CPLUSNS) CALL AFGIN2(5,0,X,SB,UPLUSNS) CALL AFGIN2(6,0,X,SB,SPLUSNS) XPDF(3) = CPLUSNS XPDF(4) = DPLUSNS XPDF(5) = UPLUSNS XPDF(6) = SPLUSNS XPDF(1) = SING ! ADD2 = XPDF(1)/FLAV UPLU2= XPDF(5)+ADD2 DPLU2=-XPDF(4)+ADD2 SPLU2=-XPDF(6)+ADD2 CPLU2=-XPDF(3)+ADD2 SING2=XPDF(1) GLU2=XPDF(7) UB=UPLUS+UPLU2*KA UB=UB/2.0 UV=UB DB=DPLUS+DPLU2*KA DB=DB/2.0 DV=DB SB=SPLUS+SPLU2*KA SB=SB/2.0 CB=CPLUS+CPLU2*KA CB=CB/2.0 SING=SING+SING2*KA GL=GLU+GLU2*KA ! DUV=MAX(ZERO,UV) DDV=MAX(ZERO,DV) DUB=MAX(ZERO,UB) DDB=MAX(ZERO,DB) DSB=MAX(ZERO,SB) DCB=MAX(ZERO,CB) DGL=MAX(ZERO,GL) ! RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE WATE32 ! 32 POINT GAUSSIAN QUADRATURE ROUTINE double precision & & X(16),W(16) double precision & & XI(32),WI(32),XX(33) COMMON/W5051I9/XI,WI,XX,NTERMS NTERMS=32 X(1)=0.048307665687738316235D0 X(2)=0.144471961582796493485D0 X(3)=0.239287362252137074545D0 X(4)=0.331868602282127649780D0 X(5)=0.421351276130635345364D0 X(6)=0.506899908932229390024D0 X(7)=0.587715757240762329041D0 X(8)=0.663044266930215200975D0 X(9)=0.732182118740289680387D0 X(10)=0.794483795967942406963D0 X(11)=0.849367613732569970134D0 X(12)=0.896321155766052123965D0 X(13)=0.934906075937739689171D0 X(14)=0.964762255587506430774D0 X(15)=0.985611511545268335400D0 X(16)=0.997263861849481563545D0 W(1)=0.096540088514727800567D0 W(2)=0.095638720079274859419D0 W(3)=0.093844399080804565639D0 W(4)=0.091173878695763884713D0 W(5)=0.087652093004403811143D0 W(6)=0.083311924226946755222D0 W(7)=0.078193895787070306472D0 W(8)=0.072345794108848506225D0 W(9)=0.065822222776361846838D0 W(10)=0.058684093478535547145D0 W(11)=0.050998059262376176196D0 W(12)=0.042835898022226680657D0 W(13)=0.034273862913021433103D0 W(14)=0.025392065309262059456D0 W(15)=0.016274394730905670605D0 W(16)=0.007018610009470096600D0 NTERMH = NTERMS/2 DO 1 I=1,NTERMH XI(I)=-X(17-I) WI(I)=W(17-I) XI(I+16)=X(I) WI(I+16)=W(I) 1 END DO DO 2 I=1,NTERMS 2 XX(I)=0.5D0*(XI(I)+1.0D0) XX(33)=1.0D0 RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE AURGAM(I,NDRV,X,S,ANS) DIMENSION F1(32),F2(32),F3(32) DIMENSION AF(3),AS(3) DIMENSION CALCO(8,20,32) COMMON/W5051I7/CALCO DATA DELTA/0.8000E-01/ ANS=0. IF(X.GT.0.9985) RETURN N=3 IS=S/DELTA+1 IF(IS.GE.17) IS=17 IS1=IS+1 IS2=IS1+1 DO 1 L=1,32 KL=L+32*NDRV F1(L)=CALCO(I,IS,KL) F2(L)=CALCO(I,IS1,KL) F3(L)=CALCO(I,IS2,KL) 1 END DO AF(1)=AFGETFV(X,F1) AF(2)=AFGETFV(X,F2) AF(3)=AFGETFV(X,F3) AS(1)=(IS-1)*DELTA AS(2)=AS(1)+DELTA AS(3)=AS(2)+DELTA CALL AFPOLIN(AS,AF,N,S,AANS,DY) ANS=AANS RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE AFGINT(I,NDRV,X,S,ANS) DIMENSION F1(32),F2(32),F3(32) DIMENSION AF(3),AS(3) DIMENSION CALCO(8,20,32) COMMON/W5051IA/CALCO DATA DELTA/0.8000E-01/ ANS=0. IF(X.GT.0.9985) RETURN N=3 IS=S/DELTA+1 ! IF(IS.GE.17) IS=17 IS1=IS+1 IS2=IS1+1 DO 1 L=1,32 KL=L+32*NDRV F1(L)=CALCO(I,IS,KL) F2(L)=CALCO(I,IS1,KL) F3(L)=CALCO(I,IS2,KL) 1 END DO AF(1)=AFGETFV(X,F1) AF(2)=AFGETFV(X,F2) AF(3)=AFGETFV(X,F3) AS(1)=(IS-1)*DELTA AS(2)=AS(1)+DELTA AS(3)=AS(2)+DELTA CALL AFPOLIN(AS,AF,N,S,AANS,DY) ANS=AANS RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE AFGIN2(I,NDRV,X,S,ANS) DIMENSION F1(32),F2(32),F3(32) DIMENSION AF(3),AS(3) DIMENSION CELCO(8,20,32) COMMON/W5051IB/CELCO DATA DELTA/0.8000E-01/ ANS=0. IF(X.GT.0.9985) RETURN N=3 IS=S/DELTA+1 ! IF(IS.GE.17) IS=17 IS1=IS+1 IS2=IS1+1 DO 1 L=1,32 KL=L+32*NDRV F1(L)=CELCO(I,IS,KL) F2(L)=CELCO(I,IS1,KL) F3(L)=CELCO(I,IS2,KL) 1 END DO AF(1)=AFGETFV(X,F1) AF(2)=AFGETFV(X,F2) AF(3)=AFGETFV(X,F3) AS(1)=(IS-1)*DELTA AS(2)=AS(1)+DELTA AS(3)=AS(2)+DELTA CALL AFPOLIN(AS,AF,N,S,AANS,DY) ANS=AANS RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION AFCPLU(X,Q2) CMS=1.5**2 BETS=1-4.*CMS*X/(1.-X)/Q2 IF(BETS.LE..0) THEN AFCPLU=.0 RETURN ENDIF BETA=SQRT(BETS) CPLU=(8.*X*(1.-X)-1.-4.*CMS*X*(1.-X)/Q2)*BETA CAU=X**2+(1.-X)**2+4.*CMS*X*(1.-3.*X)/Q2-8.*CMS**2*X**2/Q2**2 CPLU=CPLU+CAU* LOG((1.+BETA)/(1.-BETA)) AFCPLU=3.*(4./9.)*CPLU*X/(3.1415*137.) 1 RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION AFGETFV(X,FVL) ! NOUVEAU PROGRAMME D'INTERPOLATION UTILISANT UNE ROUTINE DE MATH. RECI DIMENSION FVL(32) double precision & & XI(32),WI(32),XX(33) COMMON/W5051I9/XI,WI,XX,NTERMS DIMENSION A(4),B(4) N=4 EPS=1.E-7 XAM=XX(1)-EPS XAP=XX(1)+EPS ! IF(X.LT.XAM) PRINT*,' X = ',X IF(X.GT.XAM.AND.X.LT.XAP) GOTO 50 GOTO 80 50 Y=FVL(1) GOTO 77 80 IF(X.LT.XX(2)) GOTO 51 IF(X.GT.XX(30)) GOTO 61 DO 1 I=3,30 IF(X.GT.XX(I)) GOTO 1 A(1)=XX(I-2) A(2)=XX(I-1) A(3)=XX(I) A(4)=XX(I+1) B(1)=FVL(I-2) B(2)=FVL(I-1) B(3)=FVL(I) B(4)=FVL(I+1) GOTO 70 1 CONTINUE 61 A(1)=XX(29) A(2)=XX(30) A(3)=XX(31) A(4)=XX(32) B(1)=FVL(29) B(2)=FVL(30) B(3)=FVL(31) B(4)=FVL(32) GOTO 70 51 A(1)=XX(1) A(2)=XX(2) A(3)=XX(3) A(4)=XX(4) B(1)=FVL(1) B(2)=FVL(2) B(3)=FVL(3) B(4)=FVL(4) ! 70 IF(X.GT..2.AND.X.LT..8) THEN ! CALL AFPOLIN(A,B,N,X,Y,DY) ! ELSE ! CALL AFRATIN(A,B,N,X,Y,DY) ! ENDIF 70 CONTINUE CALL AFPOLIN(A,B,N,X,Y,DY) 77 AFGETFV=Y RETURN END !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE AFPOLIN(XA,YA,N,X,Y,DY) PARAMETER (NMAX=10) DIMENSION XA(NMAX),YA(NMAX),C(NMAX),D(NMAX) Y=0. IF(N.GT.NMAX) RETURN NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 END DO Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP ! IF(DEN.EQ.0.)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 END DO RETURN END