* *-----INIT------------------------------------------------------------------ * BLOCK DATA INIT IMPLICIT REAL*8 (A-H,O,P,R-Z) IMPLICIT REAL*16(Q) REAL*8 MM,NM * COMMON/FN/NF COMMON/PARAM/PI,PIS,DELTA COMMON/QVARIA/QALPHA,QSLLC COMMON/CORRM/AVEC1(25),AVEC2(36) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/PAR/STH2,STH4,CTH2,ZM2,WM2,HM2 COMMON/QNUM/BQL,BQN,BQUQ,BQDQ,ZIU,ZID COMMON/COUPLINGS/ALPHA,GF,ALS,RS0,CONV COMMON/FMASSES/EM,MM,TM,NM,UQM,DQM,CQM,SQM,BQM COMMON/QFMASSES/QEM,QMM,QTM,QNM,QUQM,QDQM,QCQM,QSQM,QBQM COMMON/EXP/ZME,WTE,RLE,SHADE,ALE,GBGHE,GCGHE,ABE,ACE, # ST2EE,ATAUE,WME,SW2E,ALRBE,ALRCE,ST2BE COMMON/ERR/EEZM,EWT,ERL,ESHAD,EAL,EGBGH,EGCGH,EAB,EAC, # EST2E,EATAU,EWM,ESW2,EALRB,EALRC,EST2B * DATA QDELTA/9.0258247Q0/,QEPS/1.Q-90/, # QPI/3.141592653589793238462643Q0/, # QALPHA/7.2973530796Q-3/ * DATA PI/3.141592653589793238462643D0/ * *-----GF IS THE FERMI COUPLING CONSTANT (BEWARE OF THE DEFINITION) * RS0 IS THE REFERENCE ENERGY FOR ALPHA_S * NF IS THE NUMBER OF ACTIVE FLAVORS * DATA ALPHA/7.2973530796D-3/,GF/8.2476172696D-6/,RS0/91.1884D0/, # NF/5/,DELTA/9.0258247D0/ * *-----FERMION MASSES * DATA EM/0.51099906D-3/,MM/0.10565839D0/,NM/1.D-10/, # TM/1.7841D0/,UQM/0.041D0/,DQM/0.041D0/, # CQM/1.5D0/,SQM/0.15D0/,BQM/0.0D0/ * DATA QEM/0.51099906Q-3/,QMM/0.10565839Q0/,QNM/1.Q-10/, # QTM/1.7841Q0/,QUQM/0.041Q0/,QDQM/0.041Q0/, # QCQM/1.5Q0/,QSQM/0.15Q0/,QBQM/0.0Q0/ * *-----FERMION QUANTUM NUMBERS * DATA BQL/-1.D0/,BQN/0.D0/,ZIU/0.5D0/,ZID/-0.5D0/, # BQDQ/-0.333333333333333333333333D0/, # BQUQ/0.666666666666666666666666D0/ * *----EXPERIMENTAL VALUES&ERRORS * *-----'96 AVERAGE * DATA WTE/2.4946D0/,SHADE/41.508D0/,RLE/20.778D0/,ALE/0.0174D0/ DATA EWT/0.0027D0/,ESHAD/0.056D0/,ERL/0.029D0/,EAL/0.0010D0/ DATA ZME/91.1863D0/,EEZM/0.0020D0/ DATA GBGHE/0.2178D0/,ABE/0.09790D0/ DATA EGBGH/0.001144D0/,EAB/0.002314D0/ DATA GCGHE/0.1715D0/,EGCGH/0.005594D0/ DATA ST2EE/0.23061D0/,EST2E/0.00047D0/ DATA WME/80.356D0/,EWM/0.125D0/ DATA SW2E/0.2244D0/,ESW2/0.0042D0/ DATA ACE/0.07351D0/,EAC/0.004839D0/ DATA ATAUE/0.14008D0/,EATAU/0.00666D0/ DATA ALRBE/0.8628D0/,EALRB/0.04867D0/ DATA ALRCE/0.6249D0/,EALRC/0.08446D0/ DATA ST2BE/0.2320D0/,EST2B/0.0010D0/ * *-----TH * * DATA WTE/2.4963D0/,SHADE/41.466D0/,RLE/20.754D0/,ALE/0.0161D0/ * DATA EWT/0.0027D0/,ESHAD/0.056D0/,ERL/0.029D0/,EAL/0.0010D0/ * DATA ZME/91.1862D0/,EEZM/0.0020D0/ * DATA GBGHE/0.2158D0/,ABE/0.10263D0/ * DATA EGBGH/0.001144D0/,EAB/0.002314D0/ * DATA GCGHE/0.1723D0/,EGCGH/0.005594D0/ * DATA ST2EE/0.23159D0/,EST2E/0.00047D0/ * DATA WME/80.350D0/,EWM/0.125D0/ * DATA SW2E/0.2236D0/,ESW2/0.0042D0/ * DATA ACE/0.07334D0/,EAC/0.004839D0/ * DATA ATAUE/0.14645D0/,EATAU/0.00666D0/ * DATA ALRBE/0.93467D0/,EALRB/0.04867D0/ * DATA ALRCE/0.66762D0/,EALRC/0.08446D0/ * DATA ST2BE/0.2328D0/,EST2B/0.0010D0/ * *-----'95 CORRELATION MATRICES * DATA AVEC1/1.00D0,0.089D0,-0.012D0,-0.013D0,0.076D0, # 0.089D0,1.0D0,-0.143D0,-0.007D0,0.00D0, # -0.012D0,-0.143D0,1.0D0,0.152D0,0.006D0, # -0.013D0,-0.007D0,0.152D0,1.0D0,0.005D0, # 0.076D0,0.000D0,0.006D0,0.005D0,1.00D0/ DATA AVEC2/1.00D0,-.23D0,.00D0,.00D0,-.03D0,.01D0, # -.23D0,1.00D0,.04D0,-.06D0,.05D0,-.07D0, # .00D0,.04D0,1.00D0,.10D0,.04D0,.02D0, # .00D0,-.06D0,.10D0,1.00D0,.01D0,.10D0, # -.03D0,.05D0,.04D0,.01D0,1.00D0,.12D0, # .01D0,-.07D0,.02D0,.10D0,.12D0,1.00D0/ END * *------MAIN--------------------------------------------------------------- * IMPLICIT REAL*8(A-H,O-Z) CHARACTER*1 OMON,OAAS,OPOP,OPE CHARACTER*50 FILEOUT CHARACTER*1 OU1,OU2,OU3,OU4,OU5,OU6,OU7 * PARAMETER(MMDEC=20,MMDEC0=18,MNDEC=6,LIW=1,LV=MNDEC,LJ=MMDEC, # LW=6*MNDEC+MMDEC*MNDEC+2*MMDEC+MNDEC*(MNDEC-1)/2, # NIN=20,NOUT=21,NOUTC=80,NC=19) PARAMETER(NMAX1=5,NMAX12=NMAX1*NMAX1,IA1=NMAX1,IV1=NMAX1) PARAMETER(NMAX2=6,NMAX22=NMAX2*NMAX2,IA2=NMAX2,IV2=NMAX2) PARAMETER(LWORK1=25*NMAX1,LWORK2=36*NMAX2) * COMMON/AAS/SC COMMON/OE/OPE COMMON/ECM/RS COMMON/TC/TCHI2 COMMON/HIGGSM/HM COMMON/AASO/OAAS COMMON/FIXED/AX(3) COMMON/MONITOR/OMON COMMON/EXC/IUO,IUA(16) COMMON/PARAM/PI,PIS,DELTA COMMON/CORRM/AVEC1(25),AVEC2(36) COMMON/RES1/VM1(IV1,NMAX1),R1(NMAX1) COMMON/RES2/VM2(IV2,NMAX2),R2(NMAX2) COMMON/TYPE/ITYPE,IPENT,IPENA,IPENAL COMMON/SCALEP/AEMIL0,DAEMIL,BQM0,DBQM COMMON/THU/OU1,OU2,OU3,OU4,OU5,OU6,OU7 COMMON/SCALE/TQM0,DTQM,ALS0,DALS,ZM0,DZM,HM0,DHM COMMON/TH/ZMT,WTT,RLT,SHADT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ATAUT,WMT,SW2T,ALRBT,ALRCT,ST2BT COMMON/EXP/ZME,WTE,RLE,SHADE,ALE,GBGHE,GCGHE,ABE,ACE, # ST2EE,ATAUE,WME,SW2E,ALRBE,ALRCE,ST2BE COMMON/ERR/EEZM,EWT,ERL,ESHAD,EAL,EGBGH,EGCGH,EAB,EAC, # EST2E,EATAU,EWM,ESW2,EALRB,EALRC,EST2B * DIMENSION XC(MNDEC) DIMENSION CHITAB(32) DIMENSION OTAB(16,32) DIMENSION OTABMN(16),OTABMX(16) DIMENSION FJAC(MMDEC,MNDEC),FVEC(MMDEC),G(MNDEC),S(MNDEC), # V(LV,MNDEC),W(LW),X(MNDEC),IW(LIW),CJ(MNDEC),WORK(LW) DIMENSION HCNTRL(MNDEC),HESIAN(MNDEC,MNDEC),HFORW(MNDEC), # OBJGRD(MNDEC),USER(1),WORKP(MNDEC*MNDEC+MNDEC), # INFO(MNDEC),IUSER(1) DIMENSION A1(IA1,NMAX1),E1(NMAX1) DIMENSION WORKI1(LWORK1),IPIV1(NMAX1) DIMENSION A2(IA2,NMAX2),E2(NMAX2) DIMENSION WORKI2(LWORK2),IPIV2(NMAX2) * *-----OUTPUT FILE NAME * PRINT*,' CHOOSE OUTPUT FILENAME ' READ '(A)',FILEOUT OPEN(NOUT,FILE=FILEOUT,STATUS='NEW') * *-----OBSERVABLES TO BE FITTED * THE FOLLOWING OBSERVABLES ARE FITTED: * M_Z, G_Z, SIGMA_H, R, A_FB(L), * R_B, R_C, A_FB(B), A_FB(C), * ALR(B), ALR(C), * SIN^2(E), A^0_E, * M_W, 1-M_W^2/M_Z^2 * SIN^2(B) * DO IU=1,16 IUA(IU)= 1 ENDDO WRITE(NOUT,FMT=*) ' ENTER NUMBER OF UNWANTED OBSERVABLES ' PRINT*,' ENTER NUMBER OF UNWANTED OBSERVABLES ' READ*,IUO WRITE(NOUT,FMT=*) ' ENTER [0/1] FOR EACH OBSERVABLE ' PRINT*,' ENTER [0/1] FOR EACH OBSERVABLE ' READ*,(IUA(IU),IU=1,16) * *-----INPUT/OUTPUT FACILITIES * WRITE(NOUT,FMT=*) ' MONITORING [Y/N] ' PRINT*,' MONITORING [Y/N] ' READ '(A)',OMON WRITE(NOUT,FMT=*) ' PROPAGATION OF ERRORS [Y/N] ' PRINT*,' PROPAGATION OF ERRORS [Y/N] ' READ '(A)',OPE * *-----CHOICE OF FIT * PRINT*,' ENTER TYPE OF FIT: ' PRINT*,' M_Z, M_T, M_H AND ALPHA_S [6] ' PRINT*,' M_Z, M_T AND M_H [5] ' PRINT*,' M_Z, M_T AND ALPHA_S [4] ' PRINT*,' M_Z, ALPHA_S AT M_T FIXED [2] ' PRINT*,' M_Z, M_T AT ALPHA_S FIXED [3] ' READ*,ITYPE * *-----INVERSION AND DIAGONALIZATION OF THE CORRELATION MATRICES * ND1= NMAX1 DO I=1,NMAX1 DO J=1,NMAX1 K= NMAX1*(I-1)+J A1(I,J)= AVEC1(K) ENDDO ENDDO * CALL F07ADF(NMAX1,NMAX1,A1,IA1,IPIV1,INFOP1) IF(INFOP1.EQ.0) THEN CALL F07AJF(NMAX1,A1,IA1,IPIV1,WORKI1,LWORK1,INFO1) ELSE WRITE(NOUT,FMT=*) 'THE FACTOR U IS SINGULAR' STOP ENDIF * IFAIL= 1 CALL F02ABF(A1,IA1,ND1,R1,VM1,IV1,E1,IFAIL) * IF(IFAIL.NE.0) THEN WRITE(NOUT,FMT=*) ' ERROR DETECTED BY F02ABF, IFAIL = ' # ,IFAIL ENDIF * ND2= NMAX2 DO I=1,NMAX2 DO J=1,NMAX2 K= NMAX2*(I-1)+J A2(I,J)= AVEC2(K) ENDDO ENDDO * CALL F07ADF(NMAX2,NMAX2,A2,IA2,IPIV2,INFOP2) IF(INFOP2.EQ.0) THEN CALL F07AJF(NMAX2,A2,IA2,IPIV2,WORKI2,LWORK2,INFO2) ELSE WRITE(NOUT,FMT=*) 'THE FACTOR U IS SINGULAR' STOP ENDIF * IFAIL= 1 CALL F02ABF(A2,IA2,ND2,R2,VM2,IV2,E2,IFAIL) * IF(IFAIL.NE.0) THEN WRITE(NOUT,FMT=*) ' ERROR DETECTED BY F02ABF, IFAIL = ' # ,IFAIL ENDIF * TQM0= 170.D0 DTQM= 80.D0 ALS0= 0.125D0 DALS= 0.011D0 ZM0= 91.190D0 DZM= 0.15D0 BQM0= 4.8D0 DBQM= 0.05D0 AEMIL0= 128.87D0 DAEMIL= 0.1D0 HM0= 150.D0 DHM= 100.D0 * IF(ITYPE.EQ.2.OR.ITYPE.EQ.3.OR.ITYPE.EQ.4) THEN WRITE(NOUT,FMT=*) ' ENTER HIGGS MASS (GEV) ' PRINT*,' ENTER HIGGS MASS (GEV) ' READ*,HM ENDIF * IF(ITYPE.EQ.6) THEN NDEC= 6 ELSE IF(ITYPE.EQ.4.OR.ITYPE.EQ.5) THEN NDEC= 5 ELSE NDEC= 4 ENDIF * *-----INITIALIZES PARAMETERS (FIXED AND TO BE FITTED) * IF(ITYPE.EQ.2) THEN PRINT*,' ENTER M_T (GEV) AND ' PRINT*,' INITIAL GUESS FOR M_Z AND ALPHA_S ' READ*,TQM,GZM,GALS XC(1)= (GZM-ZM0)/DZM XC(2)= TQM XC(3)= HM XC(4)= (GALS-ALS0)/DALS ELSE IF(ITYPE.EQ.3) THEN PRINT*,' ENTER ALPHA_S AND ' PRINT*,' INITIAL GUESS FOR M_Z AND M_T (GEV) ' READ*,ALPHAS,GZM,GTQM XC(1)= (GZM-ZM0)/DZM XC(2)= (GTQM-TQM0)/DTQM XC(3)= HM XC(4)= ALPHAS ELSE IF(ITYPE.EQ.4) THEN PRINT*,' ENTER INITIAL GUESS FOR ' PRINT*,' M_Z, M_T (GEV) AND ALPHA_S ' READ*,GZM,GTQM,GALS XC(1)= (GZM-ZM0)/DZM XC(2)= (GTQM-TQM0)/DTQM XC(3)= HM XC(4)= (GALS-ALS0)/DALS ELSE IF(ITYPE.EQ.5) THEN PRINT*,' ENTER ALPHA_S ' READ*,ALS PRINT*,' ENTER INITIAL GUESS FOR ' PRINT*,' M_Z, M_T AND M_H (GEV) ' READ*,GZM,GTQM,GHM XC(1)= (GZM-ZM0)/DZM XC(2)= (GTQM-TQM0)/DTQM XC(4)= (GHM-HM0)/DHM XC(3)= ALS ELSE IF(ITYPE.EQ.6) THEN PRINT*,' ENTER INITIAL GUESS FOR ' PRINT*,' M_Z, M_T, M_H (GEV) AND ALPHA_S ' READ*,GZM,GTQM,GHM,GALS XC(1)= (GZM-ZM0)/DZM XC(2)= (GTQM-TQM0)/DTQM XC(3)= (GHM-HM0)/DHM XC(4)= (GALS-ALS0)/DALS ENDIF PRINT*,' ENTER INITIAL GUESS FOR ' PRINT*,' 1/ALPHA_EM, M_B (GEV) ' READ*,GAEMIL,GBQM XC(5)= (GAEMIL-AEMIL0)/DAEMIL XC(6)= (GBQM-BQM0)/DBQM * MDEC= 20 * IPENT= 0 IPENA= 0 * WRITE(NOUT,FMT=*) ' ENTER PENALTY OPTION [Y=1,N=0]' WRITE(NOUT,FMT=*) ' PENALTY FOR M_T ' PRINT*,' ENTER PENALTY OPTION [Y=1,N=0]' PRINT*,' PENALTY FOR M_T ' READ*,IPENT WRITE(NOUT,FMT=*) ' PENALTY FOR ALPHA_S ' PRINT*,' PENALTY FOR ALPHA_S ' READ*,IPENA * IPENAL= IPENT+IPENA * * *-----OPTION FOR THE SCALE IN ALPHA_S IN THE ALPHA*ALPHA_S * CORRECTIONS: IF OAAS=N THEN SCALE=M_T OTHERWISE ENTER * SC SUCH THAT SCALE=SC*M_T * PRINT*,' ENTER OPTION FOR THE SCALE OF ALPHA_S [Y/N] ' READ '(A)',OAAS IF(OAAS.EQ.'N') THEN SC= 1.D0 ELSE IF(OAAS.EQ.'Y') THEN PRINT*,' ENTER SCALE ' READ*,SC ENDIF * *-----LOOPS ON OPTIONS START * PRINT*,' ENTER OPTION FOR OPTIONS [Y/N] ' READ '(A)',OPOP * DO I1=1,2 DO I4=1,2 DO I5=1,2 DO I6=1,2 IF(I1.EQ.1) THEN OU1= 'N' ELSE IF(I1.EQ.2) THEN OU1= 'Y' ENDIF IF(I4.EQ.1) THEN OU4= 'N' ELSE IF(I4.EQ.2) THEN OU4= 'Y' ENDIF IF(I5.EQ.1) THEN OU5= 'N' ELSE IF(I5.EQ.2) THEN OU5= 'Y' ENDIF IF(I6.EQ.1) THEN OU6= 'N' ELSE IF(I6.EQ.2) THEN OU6= 'Y' ENDIF DO I7=1,2 IF(I7.EQ.1) THEN OU7= 'N' ELSE IF(I7.EQ.2) THEN OU7= 'Y' ENDIF K= 2*(2*(2*(2*(I1-1)+I4-1)+I5-1)+ # I6-1)+I7 OU2= 'N' OU3= 'Y' IF(OPOP.EQ.'N'.AND.K.NE.NC) THEN CHITAB(K)= 0.D0 GO TO 1000 ELSE CALL FITTER(MDEC,NDEC,XC,FJAC, # FVEC,G,S,V,W,X,IW, # CJ,WORK,HCNTRL,HESIAN, # HFORW,OBJGRD,USER,WORKP, # INFO,IUSER) CHITAB(K)= TCHI2 OTAB(1,K)= ZMT OTAB(2,K)= WTT OTAB(3,K)= SHADT OTAB(4,K)= RLT OTAB(5,K)= ALT OTAB(6,K)= GBGHT OTAB(7,K)= GCGHT OTAB(8,K)= ABT OTAB(9,K)= ACT OTAB(10,K)= ST2ET OTAB(11,K)= ST2BT OTAB(12,K)= ATAUT OTAB(13,K)= WMT OTAB(14,K)= SW2T OTAB(15,K)= ALRBT OTAB(16,K)= ALRCT ENDIF 1000 CONTINUE ENDDO ENDDO ENDDO ENDDO ENDDO * CHI2MN= CHITAB(1) CHI2MX= CHITAB(1) DO I=1,32 IF(CHITAB(I).LT.CHI2MN) THEN CHI2MN= CHITAB(I) ENDIF IF(CHITAB(I).GT.CHI2MX) THEN CHI2MX= CHITAB(I) ENDIF ENDDO * DO I=1,16 OTABMN(I)= OTAB(I,1) OTABMX(I)= OTAB(I,1) DO J=1,32 IF(OTAB(I,J).LT.OTABMN(I)) THEN OTABMN(I)= OTAB(I,J) ENDIF IF(OTAB(I,J).GT.OTABMX(I)) THEN OTABMX(I)= OTAB(I,J) ENDIF ENDDO ENDDO * DO I=1,16 UER= OTABMX(I)-OTAB(I,NC) DER= OTAB(I,NC)-OTABMN(I) PRINT 102,DER,UER ENDDO 102 FORMAT(1X,2E20.5) * WRITE(NOUT,FMT=100) CHITAB(NC),CHI2MN,CHI2MX PRINT 100,CHITAB(NC),CHI2MN,CHI2MX 100 FORMAT(///' CHI^2 = ',E20.5,/ # ' CHI^2 MIN = ',E20.5,/ # ' CHI^2 MAX = ',E20.5,/) * DO I=1,16 WRITE(NOUT,FMT=101) OTAB(I,NC),OTABMN(I),OTABMX(I) PRINT 101,OTAB(I,NC),OTABMN(I),OTABMX(I) PRINT*,'----------------------------------' ENDDO 101 FORMAT(///' OBS = ',E20.7,/ # ' OBS MIN = ',E20.7,/ # ' OBS MAX = ',E20.7,/) * STOP END * SUBROUTINE FITTER(MDEC,NDEC,XC,FJAC,FVEC,G,S,V,W,X,IW,CJ,WORK, # HCNTRL,HESIAN,HFORW,OBJGRD,USER,WORKP, # INFO,IUSER) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*1 OMON,OPE CHARACTER*4 CLABSX(6) * PARAMETER(LIW=1,NOUT=21) PARAMETER(NMAX1=6,IA1=NMAX1,IV1=NMAX1,LDCM=6) PARAMETER(NMAX2=6,IA2=NMAX2,IV2=NMAX2) PARAMETER(KMAX=6) * COMMON/OE/OPE COMMON/ECM/RS COMMON/TC/TCHI2 COMMON/FIXED/AX(3) COMMON/MONITOR/OMON COMMON/CGAMMA/IGAMMA COMMON/EXC/IUO,IUA(16) COMMON/PARAM/PI,PIS,DELTA COMMON/RES1/VM1(IV1,NMAX1),R1(NMAX1) COMMON/RES2/VM2(IV2,NMAX2),R2(NMAX2) COMMON/TYPE/ITYPE,IPENT,IPENA,IPENAL COMMON/SCALEP/AEMIL0,DAEMIL,BQM0,DBQM COMMON/URES/FVECCL(NMAX1),FVECCQ(NMAX2) COMMON/SCALE/TQM0,DTQM,ALS0,DALS,ZM0,DZM,HM0,DHM COMMON/TH/ZMT,WTT,RLT,SHADT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ATAUT,WMT,SW2T,ALRBT,ALRCT,ST2BT COMMON/EXP/ZME,WTE,RLE,SHADE,ALE,GBGHE,GCGHE,ABE,ACE, # ST2EE,ATAUE,WME,SW2E,ALRBE,ALRCE,ST2BE COMMON/ERR/EEZM,EWT,ERL,ESHAD,EAL,EGBGH,EGCGH,EAB,EAC, # EST2E,EATAU,EWM,ESW2,EALRB,EALRC,EST2B * DIMENSION XC(6) DIMENSION FJAC(MDEC,NDEC),FVEC(MDEC),G(NDEC),S(NDEC), # V(NDEC,NDEC),W(6*NDEC+MDEC*NDEC+2*MDEC+NDEC*(NDEC-1)/2), # X(NDEC),IW(LIW),CJ(NDEC),WORK(6*NDEC+MDEC*NDEC+2*MDEC+ # NDEC*(NDEC-1)/2),VOR(KMAX,KMAX),DELV(KMAX), # VORC(KMAX,KMAX) DIMENSION HCNTRL(NDEC),HESIAN(NDEC,NDEC),HFORW(NDEC),OBJGRD(NDEC), # USER(1),WORKP(NDEC*NDEC+NDEC),INFO(NDEC),IUSER(1) * EXTERNAL E04FCF,LSQFUN,LSQGRD,LSQMON,X02AJF,E04YCF EXTERNAL E04XAF,OBJFUN EXTERNAL X04CBF * DATA CLABSX/'X_1', ' X_2', ' X_3', ' X_4', 'X_5', 'X_6'/ * LV= NDEC LJ= MDEC LW= 6*NDEC+MDEC*NDEC+2*MDEC+NDEC*(NDEC-1)/2 N= NDEC IPENALD= 2-IPENAL M= MDEC-IPENALD MEFF= M-IUO NDOF= MEFF-N * IF(OMON.EQ.'Y') THEN IPRINT= 1 ELSE IF(OMON.EQ.'N') THEN IPRINT= -1 ENDIF * MAXCAL= 400*N ETA= 0.5D0 XTOL= 1.D-1 STEPMX= SQRT(4.D0*N) * IF(ITYPE.EQ.2) THEN X(1)= XC(1) X(2)= XC(4) X(3)= XC(5) X(4)= XC(6) AX(1)= XC(2) AX(2)= XC(3) ELSE IF(ITYPE.EQ.3) THEN X(1)= XC(1) X(2)= XC(2) X(3)= XC(5) X(4)= XC(6) AX(1)= XC(3) AX(2)= XC(4) ELSE IF(ITYPE.EQ.4) THEN X(1)= XC(1) X(2)= XC(2) X(3)= XC(4) X(4)= XC(5) X(5)= XC(6) AX(1)= XC(3) ELSE IF(ITYPE.EQ.5) THEN X(1)= XC(1) X(2)= XC(2) X(3)= XC(4) X(4)= XC(5) X(5)= XC(6) AX(1)= XC(3) ELSE IF(ITYPE.EQ.6) THEN X(1)= XC(1) X(2)= XC(2) X(3)= XC(4) X(4)= XC(5) X(5)= XC(6) X(6)= XC(3) ENDIF * IFAIL= 1 CALL E04FCF(M,N,LSQFUN,LSQMON,IPRINT,MAXCAL,ETA,XTOL,STEPMX,X, # FSUMSQ,FVEC,FJAC,LJ,S,V,LV,NITER,NF,IW,LIW,W,LW,IFAIL) IF(IFAIL.NE.0) THEN WRITE(NOUT,FMT=2) IFAIL PRINT 2,IFAIL ENDIF IF(IFAIL.NE.1) THEN TCHI2= FSUMSQ*NDOF WRITE (NOUT,FMT=3) FSUMSQ,TCHI2 PRINT 3,FSUMSQ,TCHI2 IF(OMON.EQ.'Y') THEN WRITE (NOUT,FMT=4) (X(J),J=1,N) WRITE (NOUT,FMT=5) PRINT 4,(X(J),J=1,N) PRINT 5 ENDIF ENDIF * IF(ITYPE.EQ.2) THEN ZMB= ZM0+DZM*X(1) ALSB= ALS0+DALS*X(2) TQMB= AX(1) AEMILB= AEMIL0+DAEMIL*X(3) BQMB= BQM0+DBQM*X(4) ELSE IF(ITYPE.EQ.3) THEN ZMB= ZM0+DZM*X(1) TQMB= TQM0+DTQM*X(2) ALSB= AX(2) AEMILB= AEMIL0+DAEMIL*X(3) BQMB= BQM0+DBQM*X(4) ELSE IF(ITYPE.EQ.4) THEN ZMB= ZM0+DZM*X(1) TQMB= TQM0+DTQM*X(2) ALSB= ALS0+DALS*X(3) AEMILB= AEMIL0+DAEMIL*X(4) BQMB= BQM0+DBQM*X(5) ELSE IF(ITYPE.EQ.5) THEN ZMB= ZM0+DZM*X(1) TQMB= TQM0+DTQM*X(2) HMB= HM0+DHM*X(3) AEMILB= AEMIL0+DAEMIL*X(4) BQMB= BQM0+DBQM*X(5) ELSE IF(ITYPE.EQ.6) THEN ZMB= ZM0+DZM*X(1) TQMB= TQM0+DTQM*X(2) ALSB= ALS0+DALS*X(3) AEMILB= AEMIL0+DAEMIL*X(4) BQMB= BQM0+DBQM*X(5) HMB= HM0+DHM*X(6) ENDIF IF(OMON.EQ.'Y') THEN IF(ITYPE.LE.4) THEN WRITE(NOUT,FMT=66) TQMB,ALSB,ZMB PRINT 66,TQMB,ALSB,ZMB ELSE IF(ITYPE.EQ.5) THEN WRITE(NOUT,FMT=66) TQMB,ZMB,HMB PRINT 67,TQMB,ZMB,HMB ELSE IF(ITYPE.EQ.6) THEN WRITE(NOUT,FMT=66) TQMB,ALSB,ZMB,HMB PRINT 68,TQMB,ALSB,ZMB,HMB ENDIF WRITE(NOUT,FMT=60) AEMILB,BQMB PRINT 60,AEMILB,BQMB ENDIF * CALL LSQGRD(M,N,FVEC,FJAC,LJ,G) IF(OMON.EQ.'Y') THEN * WRITE (NOUT,FMT=23) DO J=1,5 IF(J.EQ.1) THEN WRITE(NOUT,FMT=8) ELSE IF(J.EQ.2) THEN WRITE(NOUT,FMT=9) ELSE IF(J.EQ.3) THEN WRITE(NOUT,FMT=10) ELSE IF(J.EQ.4) THEN WRITE(NOUT,FMT=11) ELSE IF(J.EQ.5) THEN WRITE(NOUT,FMT=12) ENDIF WRITE (NOUT,FMT=22) FVECCL(J) ENDDO DO J=6,11 IF(J.EQ.6) THEN WRITE(NOUT,FMT=13) ELSE IF(J.EQ.7) THEN WRITE(NOUT,FMT=14) ELSE IF(J.EQ.8) THEN WRITE(NOUT,FMT=15) ELSE IF(J.EQ.9) THEN WRITE(NOUT,FMT=16) ELSE IF(J.EQ.10) THEN WRITE(NOUT,FMT=160) ELSE IF(J.EQ.11) THEN WRITE(NOUT,FMT=161) ENDIF J0= J-5 WRITE (NOUT,FMT=22) FVECCQ(J0) ENDDO DO J=12,16 IF(J.EQ.12) THEN WRITE(NOUT,FMT=17) ELSE IF(J.EQ.13) THEN WRITE(NOUT,FMT=18) ELSE IF(J.EQ.14) THEN WRITE(NOUT,FMT=19) ELSE IF(J.EQ.15) THEN WRITE(NOUT,FMT=20) ELSE IF(J.EQ.16) THEN WRITE(NOUT,FMT=21) ENDIF WRITE (NOUT,FMT=22) FVEC(J) ENDDO PRINT 23 DO J=1,5 IF(J.EQ.1) THEN PRINT 8 ELSE IF(J.EQ.2) THEN PRINT 9 ELSE IF(J.EQ.3) THEN PRINT 10 ELSE IF(J.EQ.4) THEN PRINT 11 ELSE IF(J.EQ.5) THEN PRINT 12 ENDIF PRINT 22,FVECCL(J) ENDDO DO J=6,11 IF(J.EQ.6) THEN PRINT 13 ELSE IF(J.EQ.7) THEN PRINT 14 ELSE IF(J.EQ.8) THEN PRINT 15 ELSE IF(J.EQ.9) THEN PRINT 16 ELSE IF(J.EQ.10) THEN PRINT 160 ELSE IF(J.EQ.11) THEN PRINT 161 ENDIF J0= J-5 PRINT 22,FVECCQ(J0) ENDDO DO J=12,16 IF(J.EQ.12) THEN PRINT 17 ELSE IF(J.EQ.13) THEN PRINT 18 ELSE IF(J.EQ.14) THEN PRINT 19 ELSE IF(J.EQ.15) THEN PRINT 20 ELSE IF(J.EQ.16) THEN PRINT 21 ENDIF PRINT 22,FVEC(J) ENDDO IF(IPENT.NE.0.AND.IPENA.EQ.0) THEN WRITE(NOUT,FMT=70) FVEC(19) PRINT 70,FVEC(19) ENDIF WRITE (NOUT,FMT=24) PRINT 24 * WRITE (NOUT,FMT=25) ZMT,WTT,SHADT,RLT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ST2BT,ATAUT,WMT,SW2T,ALRBT,ALRCT PRINT 25,ZMT,WTT,SHADT,RLT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ST2BT,ATAUT,WMT,SW2T,ALRBT,ALRCT * ENDIF * IFAIL= 1 JOB= -1 CALL E04YCF(JOB,M,N,FSUMSQ,S,V,LV,CJ,WORK,IFAIL) IF(IFAIL.NE.0) THEN WRITE(NOUT,FMT=200) IFAIL PRINT 200,IFAIL ENDIF IF(IFAIL.EQ.1.OR.IFAIL.EQ.2) STOP IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=120) PRINT 120 * *-----THE COVARIANCE MATRIX IN X-SPACE IS PRINTED * NCOLS= 80 INDENT= 0 IFAIL= 0 * CALL X04CBF('U','N',N,N,V,N,'E15.3',' ','C',CLABSX,'C', # CLABSX,NCOLS,INDENT,IFAIL) CALL X04ABF(1,NOUT) CALL X04CBF('U','N',N,N,V,N,'E15.3',' ','C',CLABSX,'C', # CLABSX,NCOLS,INDENT,IFAIL) CALL X04ABF(1,6) ENDIF * *-----COVARIANCE MATRIX FOR PHYSICAL VARIABLES * IF(ITYPE.EQ.2) THEN DELV(1)= DZM DELV(2)= DALS DELV(3)= DAEMIL DELV(4)= DBQM ELSE IF(ITYPE.EQ.3) THEN DELV(1)= DZM DELV(2)= DTQM DELV(3)= DAEMIL DELV(4)= DBQM ELSE IF(ITYPE.EQ.4) THEN DELV(1)= DZM DELV(2)= DTQM DELV(3)= DALS DELV(4)= DAEMIL DELV(5)= DBQM ELSE IF(ITYPE.EQ.5) THEN DELV(1)= DZM DELV(2)= DTQM DELV(3)= DHM DELV(4)= DAEMIL DELV(5)= DBQM ELSE IF(ITYPE.EQ.6) THEN DELV(1)= DZM DELV(2)= DTQM DELV(3)= DALS DELV(4)= DAEMIL DELV(5)= DBQM DELV(6)= DHM ENDIF * DO LV=1,N DO JV=1,N VOR(LV,JV)= DELV(LV)*DELV(JV)*V(LV,JV) ENDDO ENDDO DO LV=1,N DO JV=1,N VORC(LV,JV)= VOR(LV,JV)/SQRT(VOR(LV,LV)*VOR(JV,JV)) ENDDO ENDDO * *-----THE COVARIANCE MATRIX FOR PHYSICAL VARIABLES IS PRINTED * IF(OMON.EQ.'Y') THEN * WRITE(NOUT,FMT=220) PRINT 220 NCOLS= 80 INDENT= 0 IFAIL= 0 * CALL X04CBF('U','N',N,N,VORC,LDCM,'E15.3',' ','C',CLABSX,'C', # CLABSX,NCOLS,INDENT,IFAIL) CALL X04ABF(1,NOUT) CALL X04CBF('U','N',N,N,VORC,LDCM,'E15.3',' ','C',CLABSX,'C', # CLABSX,NCOLS,INDENT,IFAIL) CALL X04ABF(1,6) ENDIF * OMBETA= 31.67D0/100.D0 PB= OMBETA/2.D0 IFAIL= 1 TSTUD= -G01CAF(PB,NDOF,IFAIL) IF(OMON.EQ.'Y') THEN IF(IFAIL.GT.0) THEN WRITE(NOUT,FMT=201) IFAIL PRINT 201,IFAIL ENDIF ENDIF IF(ITYPE.EQ.4) THEN ETQM= SQRT(VOR(2,2))*TSTUD EALS= SQRT(VOR(3,3))*TSTUD EAEM= SQRT(VOR(4,4))*TSTUD EBQM= SQRT(VOR(5,5))*TSTUD IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=180) ETQM,EALS,EAEM,EBQM PRINT 180,ETQM,EALS,EAEM,EBQM ENDIF ENDIF IF(ITYPE.EQ.5) THEN ETQM= SQRT(VOR(2,2))*TSTUD EHM= SQRT(VOR(3,3))*TSTUD EAEM= SQRT(VOR(4,4))*TSTUD EBQM= SQRT(VOR(5,5))*TSTUD IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=182) ETQM,EHM,EAEM,EBQM PRINT 182,ETQM,EHM,EAEM,EBQM ENDIF ENDIF IF(ITYPE.EQ.6) THEN ETQM= SQRT(VOR(2,2))*TSTUD EALS= SQRT(VOR(3,3))*TSTUD EAEM= SQRT(VOR(4,4))*TSTUD EBQM= SQRT(VOR(5,5))*TSTUD EHM= SQRT(VOR(6,6))*TSTUD IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=181) ETQM,EALS,EAEM,EBQM,EHM PRINT 181,ETQM,EALS,EAEM,EBQM,EHM ENDIF ENDIF * IF(OPE.EQ.'Y') THEN * IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=26) PRINT 26 ENDIF DO I=1,14 IGAMMA= I MSGLVL= 0 EPSRF= -1.D0 MODE= 0 DO J=1,N HFORW(J)= -1.D0 ENDDO IFAIL= 1 CALL E04XAF(MSGLVL,N,EPSRF,X,MODE,OBJFUN,NDEC,HFORW,OBJF, # OBJGRD,HCNTRL,HESIAN,IWARN,WORKP,IUSER,USER, # INFO,IFAIL) IF(IFAIL.NE.0.AND.IFAIL.NE.2) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=202) IFAIL PRINT 202,IFAIL ENDIF ENDIF VARG= 0.D0 DO L=1,N DO K=1,N VARG= VARG+OBJGRD(L)*OBJGRD(K)*V(L,K) ENDDO ENDDO * ERRG= SQRT(VARG)*TSTUD IF(I.EQ.1) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=40) 1.D3*ERRG PRINT 40,1.D3*ERRG ENDIF ELSE IF(I.EQ.2) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=41) ERRG PRINT 41,ERRG ENDIF ELSE IF(I.EQ.3) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=42) ERRG PRINT 42,ERRG ENDIF ELSE IF(I.EQ.4) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=43) ERRG PRINT 43,ERRG ENDIF ELSE IF(I.EQ.5) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=44) 1.D3*ERRG PRINT 44,1.D3*ERRG ENDIF ELSE IF(I.EQ.6) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=444) ERRG PRINT 444,ERRG ENDIF ELSE IF(I.EQ.7) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=445) ERRG PRINT 445,ERRG ENDIF ELSE IF(I.EQ.8) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=446) ERRG PRINT 446,ERRG ENDIF ELSE IF(I.EQ.9) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=447) ERRG PRINT 447,ERRG ENDIF ELSE IF(I.EQ.10) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=448) ERRG PRINT 448,ERRG ENDIF ELSE IF(I.EQ.11) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=449) ERRG PRINT 44,ERRG ENDIF ELSE IF(I.EQ.12) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=450) ERRG PRINT 450,ERRG ENDIF ELSE IF(I.EQ.13) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=451) ERRG PRINT 451,ERRG ENDIF ELSE IF(I.EQ.14) THEN IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=452) ERRG PRINT 452,ERRG ENDIF ENDIF ENDDO * ENDIF * 2 FORMAT(///' E04FCF, ERROR EXIT TYPE',I3) 3 FORMAT(///' ON EXIT THE SUM OF SQUARES IS',E20.5/ # ' WITH A TOTAL CHI^2 OF ',E20.7/) 4 FORMAT(' AT THE POINT',2E20.8) 5 FORMAT(/' CORRESPONDING TO '/) 66 FORMAT(/' TOP MASS (GEV) =',E20.7/ # ' ALPHA_S =',E20.7/ # ' Z MASS (GEV) =',E20.7) 67 FORMAT(/' TOP MASS (GEV) =',E20.7/ # ' Z MASS (GEV) =',E20.7/ # ' H MASS (GEV) =',E20.7) 68 FORMAT(/' TOP MASS (GEV) =',E20.7/ # ' ALPHA_S =',E20.7/ # ' Z MASS (GEV) =',E20.7/ # ' H MASS (GEV) =',E20.7) 60 FORMAT(/' 1/ALPHA_EM(L) =',E20.7/ # ' B MASS (GEV) =',E20.7) 70 FORMAT(//' RESIDUAL OF M_T PENALTY =',E20.5//) 23 FORMAT(//' THE RESIDUALS ARE: ') 8 FORMAT(/' M_Z ') 9 FORMAT(/' G_Z ') 10 FORMAT(/' SIGMA_0 ') 11 FORMAT(/' R_L ') 12 FORMAT(/' A_FB(L) ') 13 FORMAT(/' R_B ') 14 FORMAT(/' R_C ') 15 FORMAT(/' A_FB(B) ') 16 FORMAT(/' A_FB(C) ') 160 FORMAT(/' A_LR(B) ') 161 FORMAT(/' A_LR(C) ') 17 FORMAT(/' SIN^2(E) ') 18 FORMAT(/' A_TAU ') 19 FORMAT(/' M_W ') 20 FORMAT(/' SINW^2 ') 21 FORMAT(/' SIN^2(B) ') 22 FORMAT(' ',E20.7) 24 FORMAT(//' THE OBSERVABLES AT THE SOLUTION ARE: '/) 25 FORMAT(/' M_Z (GEV) =',E20.7,/ # ' G_Z (GEV) =',E20.7,/ # ' SIGMA_0 (NB) =',E20.7,/ # ' R_L =',E20.7,/ # ' A_FB(L) =',E20.7,/ # ' R_B =',E20.7,/ # ' R_C =',E20.7,/ # ' A_FB(B) =',E20.7,/ # ' A_FB(C) =',E20.7,/ # ' SIN^2(E) =',E20.7,/ # ' SIN^2(B) =',E20.7,/ # ' A_TAU =',E20.7,/ # ' M_W (GEV) =',E20.7,/ # ' SINW^2 =',E20.7,/ # ' A_LR(B) =',E20.7,/ # ' A_LR(C) =',E20.7,///) 200 FORMAT(///' E04YCF, ERROR EXIT TYPE',I3) 120 FORMAT(//' ESTIMATE OF THE VARIANCES OF REGRESSION COEFF. #(X-SPACE) # '/) 220 FORMAT(//' ESTIMATE OF THE VARIANCES OF REGRESSION COEFF. '/) 201 FORMAT(///' G01CAF, ERROR EXIT TYPE',I3) 26 FORMAT(//' 68% CL ERRORS '/) 202 FORMAT(///' E04XAF, ERROR EXIT TYPE',I3) 40 FORMAT(' GAMMA_L (MEV) =',E20.7) 41 FORMAT(' SIN^2(E) =',E20.7) 42 FORMAT(' SIN^2(B) =',E20.7) 43 FORMAT(' R_B =',E20.7) 44 FORMAT(' W MASS (MEV) =',E20.7) 444 FORMAT(' R_C =',E20.7) 445 FORMAT(' R_H =',E20.7) 446 FORMAT(' SIGMA_0 (NB) =',E20.7) 447 FORMAT(' A_FB(L) =',E20.7) 448 FORMAT(' A_FB(B) =',E20.7) 449 FORMAT(' A_FB(C) =',E20.7) 450 FORMAT(' A_LR(B) =',E20.7) 451 FORMAT(' A_LR(C) =',E20.7) 452 FORMAT(' A_LR(L) =',E20.7) 180 FORMAT(/' 68% CL INTERVALS ',// # ' TOP MASS (GEV) =',E20.7/ # ' ALPHA_S =',E20.7/ # ' 1/ALPHA_EM(L) =',E20.7/ # ' B MASS (GEV) =',E20.7) 182 FORMAT(/' 68% CL INTERVALS ',// # ' TOP MASS (GEV) =',E20.7/ # ' HIGGS MASS (GEV) =',E20.7/ # ' 1/ALPHA_EM(L) =',E20.7/ # ' B MASS (GEV) =',E20.7) 181 FORMAT(/' 68% CL INTERVALS ',// # ' TOP MASS (GEV) =',E20.7/ # ' ALPHA_S =',E20.7/ # ' 1/ALPHA_EM(L) =',E20.7/ # ' B MASS (GEV) =',E20.7/ # ' H MASS (GEV) =',E20.7) * RETURN END * *-----LSQFUN----------------------------------------------------------- * REQUIRED BY FITTER * SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVECC,IW,LIW,W,LW) IMPLICIT REAL*8(A-H,O-Z) * PARAMETER(NMAX1=5,IA1=NMAX1,IV1=NMAX1,NMAX2=6,IA2=NMAX2, # IV2=NMAX2,NIN=20,NOUT=21) * COMMON/ABM/ABQM COMMON/HIGGSM/HM COMMON/FIXED/AX(3) COMMON/AAEM/OIAEM COMMON/EXC/IUO,IUA(16) COMMON/RES1/VM1(IV1,NMAX1),R1(NMAX1) COMMON/RES2/VM2(IV2,NMAX2),R2(NMAX2) COMMON/TYPE/ITYPE,IPENT,IPENA,IPENAL COMMON/URES/FVECCL(NMAX1),FVECCQ(NMAX2) COMMON/SCALE/TQM0,DTQM,ALS0,DALS,ZM0,DZM,HM0,DHM COMMON/TH/ZMT,WTT,RLT,SHADT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ATAUT,WMT,SW2T,ALRBT,ALRCT,ST2BT COMMON/EXP/ZME,WTE,RLE,SHADE,ALE,GBGHE,GCGHE,ABE,ACE, # ST2EE,ATAUE,WME,SW2E,ALRBE,ALRCE,ST2BE COMMON/ERR/EEZM,EWT,ERL,ESHAD,EAL,EGBGH,EGCGH,EAB,EAC, # EST2E,EATAU,EWM,ESW2,EALRB,EALRC,EST2B * DIMENSION FVECC(M),W(LW),XC(N),IW(LIW) * CALL EWEAK(N,XC) * *-----THIS IS THE NORMALIZED CHI^2 (CHI^2/DOF) * MEFF= M-IUO CHIN= SQRT((MEFF-N)*1.D0) FVECCL(1)= (ZMT-ZME)/EEZM/CHIN FVECCL(2)= (WTT-WTE)/EWT/CHIN FVECCL(3)= (SHADT-SHADE)/ESHAD/CHIN FVECCL(4)= (RLT-RLE)/ERL/CHIN FVECCL(5)= (ALT-ALE)/EAL/CHIN FVECCQ(1)= (GBGHT-GBGHE)/EGBGH/CHIN FVECCQ(2)= (GCGHT-GCGHE)/EGCGH/CHIN FVECCQ(3)= (ABT-ABE)/EAB/CHIN FVECCQ(4)= (ACT-ACE)/EAC/CHIN FVECCQ(5)= (ALRBT-ALRBE)/EALRB/CHIN FVECCQ(6)= (ALRCT-ALRCE)/EALRC/CHIN * *-----THE CORRELATION MATRICES ARE TAKEN INTO ACCOUNT * ND1= NMAX1 DO I=1,ND1 FVECC(I)= 0.D0 DO J=1,ND1 FVECC(I)= FVECC(I)+VM1(J,I)*FVECCL(J) ENDDO FVECC(I)= FVECC(I)*SQRT(R1(I)) ENDDO ND2= NMAX1+NMAX2 DO I=ND1+1,ND2 I0= I-ND1 FVECC(I)= 0.D0 DO J=1,ND2 FVECC(I)= FVECC(I)+VM2(J,I0)*FVECCQ(J) ENDDO FVECC(I)= FVECC(I)*SQRT(R2(I0)) ENDDO * FVECC(12)= (ST2ET-ST2EE)/EST2E/CHIN FVECC(13)= (ATAUT-ATAUE)/EATAU/CHIN FVECC(14)= (WMT-WME)/EWM/CHIN FVECC(15)= (SW2T-SW2E)/ESW2/CHIN FVECC(16)= (ST2BT-ST2BE)/EST2B/CHIN FVECC(17)= (ABQM-4.7D0)/0.2D0/CHIN FVECC(18)= (OIAEM-128.89D0)/0.09D0/CHIN * DO I=1,16 FVECC(I)= IUA(I)*FVECC(I) ENDDO * *-----PENALTY FUNCTIONS * IF(IPENAL.EQ.0) RETURN * IF(IPENT.NE.0) THEN IF(ITYPE.EQ.2) THEN PTQM= AX(1) ELSE PTQM= TQM0+DTQM*XC(2) ENDIF ENDIF IF(IPENA.NE.0) THEN IF(ITYPE.EQ.2) THEN PALS= ALS0+DALS*XC(2) ELSE IF(ITYPE.EQ.6) THEN PALS= ALS0+DALS*XC(3) ELSE IF(ITYPE.EQ.3) THEN PALS= AX(2) ELSE IF(ITYPE.EQ.5) THEN PALS= AX(1) ENDIF ENDIF * IF(IPENT.NE.0.AND.IPENA.EQ.0) THEN FVECC(19)= (PTQM-175.D0)/6.D0/CHIN ENDIF IF(IPENT.EQ.0.AND.IPENA.NE.0) THEN FVECC(19)= (PALS-0.119D0)/0.006D0/CHIN ENDIF IF(IPENT.NE.0.AND.IPENA.NE.0) THEN FVECC(19)= (PTQM-175.D0)/6.D0/CHIN FVECC(20)= (PALS-0.119D0)/0.006D0/CHIN ENDIF * RETURN END * *-----LSQMON------------------------------------------------------------ * REQUIRED BY FITTER * SUBROUTINE LSQMON(M,N,XC,FVECC,FJACC,LJC,S,IGRADE,NITER,NF, # IW,LIW,W,LW) IMPLICIT REAL*8(A-H,O-Z) * PARAMETER(NDEC=6,NOUT=21) * DIMENSION FJACC(LJC,N),FVECC(M),S(N),W(LW),XC(N),IW(LIW), # G(NDEC) * EXTERNAL F06EAF,LSQGRD * FSUMQ= F06EAF(M,FVECC,1,FVECC,1) CALL LSQGRD(M,N,FVECC,FJACC,LJC,G) GTG= F06EAF(N,G,1,G,1) WRITE(NOUT,FMT=1) NITER,NF,FSUMQ,GTG,IGRADE PRINT 1,NITER,NF,FSUMQ,GTG,IGRADE WRITE(NOUT,FMT=2) PRINT 2 DO J=1,N WRITE(NOUT,FMT=3) XC(J),G(J),S(J) PRINT 3,XC(J),G(J),S(J) ENDDO * 1 FORMAT(///' ITNS F EVALS SUMSQ GTG ', * ' GRADE',/' ',I4,6X,I5,6X,1P,E13.5,6X,1P,E9.1,6X,I3) 2 FORMAT(/' X G SINGULAR VALU', * 'ES') 3 FORMAT(' ',1P,E13.5,10X,1P,E9.1,10X,1P,E9.1) * RETURN END * *-----LSQGRD-------------------------------------------------------------- * REQUIRED BY FITTER * SUBROUTINE LSQGRD(M,N,FVECC,FJACC,LJC,G) IMPLICIT REAL*8(A-H,O-Z) * DIMENSION FJACC(LJC,N),FVECC(M),G(N) * DO J=1,N SUM= 0.D0 DO I=1,M SUM= SUM+FJACC(I,J)*FVECC(I) ENDDO G(J)= SUM+SUM ENDDO RETURN END * *-----OBJFUN--------------------------------------------------------- * CALLS WIDTH IN ORDER TO COMPUTE THE GRADIENT VECTORS OF * SUBROUTINE OBJFUN(MODE,N,X,OBJF,OBJGRD,NSTATE,IUSER,USER) IMPLICIT REAL*8(A-H,O-Z) * COMMON/ABM/ABQM COMMON/HIGGSM/HM COMMON/AAEM/OIAEM COMMON/FIXED/AX(3) COMMON/CGAMMA/IGAMMA COMMON/PARAM/PI,PIS,DELTA COMMON/TYPE/ITYPE,IPENT,IPENA,IPENAL COMMON/SCALEP/AEMIL0,DAEMIL,BQM0,DBQM COMMON/SCALE/TQM0,DTQM,ALS0,DALS,ZM0,DZM,HM0,DHM * DIMENSION OBJGRD(N),USER(1),X(N),IUSER(1) * IF(ITYPE.EQ.2) THEN ZM= ZM0+DZM*X(1) TQM= AX(1) HM= AX(2) ALPHAS= ALS0+DALS*X(2) ABQM= BQM0+DBQM*X(4) OIAEM= AEMIL0+DAEMIL*X(3) ELSE IF(ITYPE.EQ.3) THEN ZM= ZM0+DZM*X(1) TQM= TQM0+DTQM*X(2) HM= AX(1) ALPHAS= AX(2) ABQM= BQM0+DBQM*X(4) OIAEM= AEMIL0+DAEMIL*X(3) ELSE IF(ITYPE.EQ.4) THEN ZM= ZM0+DZM*X(1) TQM= TQM0+DTQM*X(2) HM= AX(1) ALPHAS= ALS0+DALS*X(3) ABQM= BQM0+DBQM*X(5) OIAEM= AEMIL0+DAEMIL*X(4) ELSE IF(ITYPE.EQ.5) THEN ZM= ZM0+DZM*X(1) TQM= TQM0+DTQM*X(2) HM= HM0+DHM*X(3) ALPHAS= AX(1) ABQM= BQM0+DBQM*X(5) OIAEM= AEMIL0+DAEMIL*X(4) ELSE IF(ITYPE.EQ.6) THEN ZM= ZM0+DZM*X(1) TQM= TQM0+DTQM*X(2) HM= HM0+DHM*X(6) ALPHAS= ALS0+DALS*X(3) ABQM= BQM0+DBQM*X(5) OIAEM= AEMIL0+DAEMIL*X(4) ENDIF * CALL INITIAL * CALL WIDTH (ZM,TQM,HM,ALPHAS,WT,WH,WB,WC,RH,SHAD, # OWM,AFBEFF,AFBBEFF,ALREFF,AFBCEFF,ALRBEFF, # ALRCEFF,ST2E,ST2B) * IF(IGAMMA.EQ.1) THEN OBJF= WT ELSE IF(IGAMMA.EQ.2) THEN OBJF= ST2E ELSE IF(IGAMMA.EQ.3) THEN OBJF= ST2B ELSE IF(IGAMMA.EQ.4) THEN OBJF= WB/WH ELSE IF(IGAMMA.EQ.5) THEN OBJF= OWM ELSE IF(IGAMMA.EQ.6) THEN OBJF= WC/WH ELSE IF(IGAMMA.EQ.7) THEN OBJF= RH ELSE IF(IGAMMA.EQ.8) THEN OBJF= SHAD ELSE IF(IGAMMA.EQ.9) THEN OBJF= AFBEFF ELSE IF(IGAMMA.EQ.10) THEN OBJF= AFBBEFF ELSE IF(IGAMMA.EQ.11) THEN OBJF= AFBCEFF ELSE IF(IGAMMA.EQ.12) THEN OBJF= ALRBEFF ELSE IF(IGAMMA.EQ.13) THEN OBJF= ALRCEFF ELSE IF(IGAMMA.EQ.14) THEN OBJF= ALREFF ENDIF * RETURN END * *-----EWEAK------------------------------------------------------------- * WEAK CORRECTIONS ARE COMPUTED * SUBROUTINE EWEAK(N,XC) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*1 OMON,OPTM * PARAMETER(NOBS=10,NOUT=21) * COMMON/ECM/RS COMMON/ABM/ABQM COMMON/HIGGSM/HM COMMON/AAEM/OIAEM COMMON/FIXED/AX(3) COMMON/MONITOR/OMON COMMON/PARAM/PI,PIS,DELTA COMMON/TYPE/ITYPE,IPENT,IPENA,IPENAL COMMON/SCALEP/AEMIL0,DAEMIL,BQM0,DBQM COMMON/SCALE/TQM0,DTQM,ALS0,DALS,ZM0,DZM,HM0,DHM COMMON/TH/ZMT,WTT,RLT,SHADT,ALT,GBGHT,GCGHT,ABT,ACT, # ST2ET,ATAUT,WMT,SW2T,ALRBT,ALRCT,ST2BT * DIMENSION XC(N),SIGMA0(NOBS),C(20) * *-----INITIALIZES THE VARIABLES * IF(ITYPE.EQ.2) THEN ZM= ZM0+DZM*XC(1) TQM= AX(1) HM= AX(2) ALPHAS= ALS0+DALS*XC(2) ABQM= BQM0+DBQM*XC(4) OIAEM= AEMIL0+DAEMIL*XC(3) ELSE IF(ITYPE.EQ.3) THEN ZM= ZM0+DZM*XC(1) TQM= TQM0+DTQM*XC(2) HM= AX(1) ALPHAS= AX(2) ABQM= BQM0+DBQM*XC(4) OIAEM= AEMIL0+DAEMIL*XC(3) ELSE IF(ITYPE.EQ.4) THEN ZM= ZM0+DZM*XC(1) TQM= TQM0+DTQM*XC(2) HM= AX(1) ALPHAS= ALS0+DALS*XC(3) ABQM= BQM0+DBQM*XC(5) OIAEM= AEMIL0+DAEMIL*XC(4) ELSE IF(ITYPE.EQ.5) THEN ZM= ZM0+DZM*XC(1) TQM= TQM0+DTQM*XC(2) ALPHAS= AX(1) HM= HM0+DHM*XC(3) ABQM= BQM0+DBQM*XC(5) OIAEM= AEMIL0+DAEMIL*XC(4) ELSE IF(ITYPE.EQ.6) THEN ZM= ZM0+DZM*XC(1) TQM= TQM0+DTQM*XC(2) HM= HM0+DHM*XC(6) ALPHAS= ALS0+DALS*XC(3) ABQM= BQM0+DBQM*XC(5) OIAEM= AEMIL0+DAEMIL*XC(4) ENDIF ZMT= ZM * CALL INITIAL * IF(OMON.EQ.'Y') THEN WRITE(NOUT,FMT=6) ZM,ALPHAS,TQM,HM PRINT 6,ZM,ALPHAS,TQM,HM WRITE(NOUT,FMT=30) PRINT 30 WRITE(NOUT,FMT=9)(XC(I),I=1,N) PRINT 9, (XC(I),I=1,N) ENDIF * 6 FORMAT(/' CURRENT VALUES FOR THE PARAMETERS ARE: '/ # ' Z MASS (GEV) = ',E20.5,' ALPHA_S = ',E20.5/ # ' TOP MASS (GEV) = ',E20.5,' HIGGS MASS (GEV) = ',E20.5) 30 FORMAT(//' CURRENT FITTING (SCALED) PARAMETERS ARE: ',/) 9 FORMAT (E20.5) * *-----THE Z0 PARAMETERS ARE COMPUTED * CALL WIDTH (ZM,TQM,HM,ALPHAS,WT,WH,WB,WC,RH,SHAD, # OWM,AFBEFF,AFBBEFF,ALREFF,AFBCEFF,ALRBEFF, # ALRCEFF,ST2E,ST2B) * WMT= OWM SW2T= 1.D0-OWM*OWM/ZM/ZM GBGHT= WB/WH GCGHT= WC/WH ALT= AFBEFF ALRBT= ALRBEFF ALRCT= ALRCEFF ABT= AFBBEFF ACT= AFBCEFF RLT= RH SHADT= SHAD WTT= WT WHT= WH ATAUT= ALREFF ST2ET= ST2E ST2BT= ST2B OST2E= ST2E * RETURN END * *-----INITIAL------------------------------------------------------------- *-----MISCELLANEA * SUBROUTINE INITIAL IMPLICIT REAL*8 (A-H,O,P,R-Z) IMPLICIT REAL*16(Q) REAL*8 MM,NM,MM2,NM2 * COMMON/FN/NF COMMON/ABM/ABQM COMMON/PARAM/PI,PIS,DELTA COMMON/QVARIA/QALPHA,QSLLC COMMON/SHARE/AEXP,TAEXP,HOF,FPI COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/COUPLINGS/ALPHA,GF,ALS,RS0,CONV COMMON/FMASSES/EM,MM,TM,NM,UQM,DQM,CQM,SQM,BQM COMMON/QFMASSES/QEM,QMM,QTM,QNM,QUQM,QDQM,QCQM,QSQM,QBQM COMMON/FMASSES2/EM2,MM2,TM2,NM2,UQM2,DQM2,CQM2,SQM2,BQM2,TQM2 COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * *-----INITIALIZATION OF SHARABLE VARIABLES * QPIS= QPI*QPI PIS= PI*PI FPI= 4.D0*PI AEXP= ALPHA/FPI TAEXP= 2.D0*AEXP QSLLC= 1.Q0+3.Q0/4.Q0*QALPHA/QPI HOF= 19.D0-2.D0*PIS BQM= ABQM QBQM= BQM*1.D15*1.Q-15 * *-----FERMION MASSES * EM2= EM*EM MM2= MM*MM TM2= TM*TM NM2= NM*NM UQM2= UQM*UQM DQM2= DQM*DQM CQM2= CQM*CQM SQM2= SQM*SQM BQM2= BQM*BQM * QEM2= QEM*QEM QMM2= QMM*QMM QTM2= QTM*QTM QNM2= QNM*QNM QUQM2= QUQM*QUQM QDQM2= QDQM*QDQM QCQM2= QCQM*QCQM QSQM2= QSQM*QSQM QBQM2= QBQM*QBQM * RETURN END * SUBROUTINE WIDTH(ZM,TQM,HM,OALS,WT,WH,WB,WC,RH,SHAD, # OWM,AFBEFF,AFBBEFF,ALREFF,AFBCEFF,ALRBEFF, # ALRCEFF,ST2E,ST2B) IMPLICIT REAL*8 (A-H,O,P,R-Z) IMPLICIT REAL*16(Q) REAL*8 MM,NM,MM2,NM2 * CHARACTER*1 OAAS CHARACTER*1 OU1,OU2,OU3,OU4,OU5,OU6,OU7 * PARAMETER(NOUT=21,NPO=23) * COMMON/FN/NF COMMON/AAS/SC COMMON/AASO/OAAS COMMON/AFJTR/ALST COMMON/SSCAL/QSPSC COMMON/SUB/OSINT2,ODSWW COMMON/PARAM/PI,PIS,DELTA COMMON/MIX/QALST,QALSTZ,QALSTS COMMON/SHARE/AEXP,TAEXP,HOF,FPI COMMON/QNUM/BQL,BQN,BQUQ,BQDQ,ZIU,ZID COMMON/PAR/STH2,STH4,CTH2,ZM2,WM2,HM2 COMMON/MIXC/QV1,QA1,QF1,QV1P,QA1P,QV1I COMMON/THU/OU1,OU2,OU3,OU4,OU5,OU6,OU7 COMMON/COUPLINGS/ALPHA,GF,ALS,RS0,CONV COMMON/FMASSES/EM,MM,TLM,NM,UQM,DQM,CQM,SQM,BQM COMMON/QFMASSES/QEM,QMM,QTM,QNM,QUQM,QDQM,QCQM,QSQM,QBQM COMMON/FMASSES2/EM2,MM2,TLM2,NM2,UQM2,DQM2,CQM2,SQM2,BQM2,TQM2 COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 COMMON/QCDCORR/VCORQ,ACORU,ACORD,ACORB,RBM2,RCM2,VCMB,ACMB,VCMC, # ACMC,ACMT,ALSR,CAQCDB,CAQCDC,CAMB,CAMC,CAMT,ACMM, # QCDND * DIMENSION XA(17),YA(17),ZA(17) DIMENSION ZW(9),ZWQCDS(9) * DATA YA/0.D0,0.D0,-2.70D0,-4.13D0,-5.25D0,-6.18D0, # -6.95D0,-7.62D0,-8.19D0,-8.68D0,-9.11D0,-9.48D0, # -9.81D0,-10.1D0,-10.4D0,-10.6D0,-10.8D0/ DATA ZA/0.D0,0.D0,3.90D0,2.83D0,2.16D0,1.73D0, # 1.47D0,1.32D0,1.24D0,1.23D0,1.26D0,1.33D0, # 1.42D0,1.53D0,1.66D0,1.80D0,1.95D0/ DATA Z3/1.20205690315959428540D0/ DATA SCS2/0.260434137632162099D0/ DATA SCD3/-3.027009493987652020D0/ DATA SCB4/-1.762800087073770086D0/ DATA RZ5/1.03692775514336992633D0/ * *-----VARIA * Z2= PIS/6.D0 Z4= PIS*PIS/90.D0 DO I=1,9 ZWQCDS(I)= 0.D0 ENDDO * *-----MASSES * ZM2= ZM*ZM QZM= ZM*1.D15*1.Q-15 QZM2= QZM*QZM ZM3= ZM2*ZM TQM2= TQM*TQM QTQM= TQM*1.D15*1.Q-15 QTQM2= QTQM*QTQM HM2= HM*HM QHM= HM*1.D15*1.Q-15 QHM2= QHM*QHM XI= TQM2/ZM2 XI2= XI*XI ALS= OALS * *-----LO SIN^2 * ARG0= 1.D0-2.D0*PI*ALPHA/GF/ZM2 ST02= 0.5D0*(1.D0-SQRT(ARG0)) ST04= ST02*ST02 CT02= 1.D0-ST02 * GWEAK= GF*ZM2 GWEAK2= GWEAK*GWEAK GFC= GWEAK/2.D0/PIS * *-----QCD CORRECTION TO RHO * IF(TQM.GT.ZM) THEN RSQCD= TQM ELSE RSQCD= ZM ENDIF ALST= RALPHAS(RS0,RSQCD,ALS,NF) ALSTZ= RALPHAS(RS0,ZM,ALS,NF) SPSC= SC*TQM QSPSC= SPSC*1.D15*1.Q-15 ALSTS= RALPHAS(RS0,SPSC,ALS,NF) SCLU= 0.248D0*TQM ALSTU= RALPHAS(RS0,SCLU,ALS,NF) IF(OU7.EQ.'N') THEN RSCEFF= 0.D0 ELSE IF(OU7.EQ.'Y') THEN RSCEFF= 3.D0/4.D0*XI*2.86D0*(ALSTU-ALST)/PI ENDIF * QALST= ALST*1.D15*1.Q-15 QALSTZ= ALSTZ*1.D15*1.Q-15 QALSTS= ALSTS*1.D15*1.Q-15 * *-----IMPROVED SIN^2 FROM ALPHA,G_F AND M_Z * P2Z= -ZM2 QP2Z= -QZM2 WM2= ZM2*CT02 WM= SQRT(WM2) QWM= WM*1.D15*1.Q-15 QWM2= QWM*QWM * *-----THE ALPHA*ALPHA_S CORRECTIONS ARE INITIALIZED * CALL ALALS(QP2Z,QV1,QA1,QDF1,QV1P,QA1P,QV1I) * CALL PSELF(QP2Z,PGGFZ,PGGLQZ,PGGBZ,PGGNPZ,PPGGZ,PPGGNPZ, # PGGI,PGZI) OSINT2= ST02 CALL VBSELF0(PGGF0,PGGB0,SWW0F,SWW0B) PIFZ= PGGFZ-PGGF0+PGGNPZ AEXPH= AEXP/(1.D0-AEXP*PIFZ) AEXPH2= AEXPH*AEXPH ALPHAH= 4.D0*PI*AEXPH * ARG= 1.D0-2.D0*PI*ALPHA/GF/ZM2/(1.D0-AEXP*PIFZ) * STH2= 0.5D0*(1.D0-SQRT(ARG)) STH2M= STH2 CTH2= 1.D0-STH2 STH4= STH2*STH2 STH6= STH4*STH2 CTH4= CTH2*CTH2 CSTH2= CTH2*STH2 CMSTH2= CTH2-STH2 TSCTH2= 2.D0/CTH2 FSCTH2= 4.D0/CTH2 CTTH2= CTH2/STH2 FACT= 1.D0/CSTH2 * *-----GIVES THE W MASS TO APPEAR INSIDE LOOPS * WM2= ZM2*CTH2 WM= SQRT(WM2) QWM= WM*1.D15*1.Q-15 QWM2= QWM*QWM DG= 2.D0*(3.D0/STH2+(7.D0-4.D0*STH2)/4.D0/STH4*LOG(CTH2)) ODSWW= STH2*WM2*DG OSINT2= STH2 * DELW= DELTA-LOG(WM2) DELZ= DELTA-LOG(ZM2) * *-----COMPUTES THE BOSONIC PART OF PI WITH THE CORRECT W MASS * CALL PSELF(QP2Z,PGGFZ,PGGLQZ,PGGBZ,PGGNPZ,PPGGZ,PPGGNPZ, # PGGI,PGZI) CALL VBSELF0(PGGF0,PGGB0,SWW0F,SWW0B) PIBZ= PGGBZ-PGGB0 SWW0= SWW0F+SWW0B * *-----COMPUTES RHO * CALL VBSELF(QP2Z,S3GFZ,S33FZ,S3GBZ,S33BZ,SP3GZ,SP33Z) * RESS2= CMSTH2*(S3GFZ+S3GBZ+STH2*ZM2*(PGGFZ+PGGLQZ+PGGBZ)) DSIGF= (SWW0B+S3GBZ-S33BZ)/ZM2+CSTH2*PIBZ-RESS2/ZM2 SIGF= (SWW0F+S3GFZ-S33FZ)/ZM2 SIGFP= (SWW0F+S3GFZ-S33FZ)/ZM2+DSIGF * *-----O(M_T^4) CORRECTIONS TO RHO AND TAU, ANALYTICAL IF M_H > 3 M_T, * OTHERWISE INTERPOLATION * IF(HM.GT.3.D0*TQM) THEN RTH= TQM2/HM2 RTH2= RTH*RTH RTHL= LOG(RTH) RTHL2= RTHL*RTHL HOC= 49.D0/4.D0+PIS+27.D0/2.D0*RTHL+3.D0/2.D0*RTHL2+ # 1.D0/3.D0*RTH*(2.D0-12.D0*PIS+12.D0*RTHL-27.D0*RTHL2)+ # 1.D0/48.D0*RTH2*(1613.D0-240.D0*PIS-1500.D0*RTHL- # 720.D0*RTHL2) TCOR= 1.D0/144.D0*(311.D0+24.D0*PIS+282.D0*RTHL+90.D0*RTHL2- # 4.D0*RTH*(40.D0+6.D0*PIS+15.D0*RTHL+18.D0*RTHL2)+ # 3.D0/100.D0*RTH2*(24209.D0-6000.D0*PIS-45420.D0*RTHL- # 18000.D0*RTHL2)) ELSE * NINT= 17 XTH= ABS(HM)/TQM XA(1)= 0.D0 XA(2)= 0.02D0 YA(1)= -0.74D0 YA(2)= YA(1)-4.D0*PI*XA(2) DO I=3,17 XA(I)= 0.D0+0.2D0*(I-2) ENDDO CALL POLINT(XA,YA,NINT,XTH,HOC,DHOC) * NINT= 17 XTH= ABS(HM)/TQM XA(1)= 0.D0 XA(2)= 0.02D0 ZA(1)= 5.71D0 ZA(2)= ZA(1)+4.D0*PI*XA(2) DO I=3,17 XA(I)= 0.D0+0.2D0*(I-2) ENDDO CALL POLINT(XA,ZA,NINT,XTH,TCOR,DTCOR) ENDIF * IF(OU7.EQ.'N') THEN NFP1= NF+1 DQCD3= 157.D0/648.D0-3313.D0/162.D0*Z2-308.D0/27.D0*Z3+ # 1051.D0/90.D0*Z4-4.D0/3.D0*Z2*LOG(2.D0)+ # 441.D0/8.D0*SCS2-1.D0/9.D0*SCB4-1.D0/18.D0*SCD3- # (1.D0/18.D0-13.D0/9.D0*Z2+4.D0/9.D0*Z3)*NFP1- # (11.D0/6.D0-NF/9.D0)*(1.D0+2.D0*Z2)*2.D0*LOG(SC) DQCD3= DQCD3*ALSTS*ALSTS/PIS DRHO3= -3.D0/4.D0*XI*DQCD3 ELSE IF(OU7.EQ.'Y') THEN DRHO3= 0.D0 ENDIF SIGFC= SIGF+RSCEFF-3.D0/16.D0*GFC*HOC*XI2+DRHO3 SIGFPC= SIGFP+RSCEFF-3.D0/16.D0*GFC*HOC*XI2+DRHO3 RHO= 1.D0/(1.D0+GFC*SIGFC) RHOP= 1.D0/(1.D0+GFC*SIGFPC) * *----COMPUTES THE EFFECTIVE SIN^2 * ARGR= 1.D0-2.D0*PI*ALPHA/GF/ZM2/RHO* # 1.D0/(1.D0-AEXP*PIFZ) ARGRP= 1.D0-2.D0*PI*ALPHA/GF/ZM2/RHOP* # 1.D0/(1.D0-AEXP*PIFZ) * STR2= 0.5D0*(1.D0-SQRT(ARGRP)) CTR2= 1.D0-STR2 STR2M= 0.5D0*(1.D0-SQRT(ARGR)) CTR2M= 1.D0-STR2M * IF(OU5.EQ.'Y') THEN STH2= STR2 CTH2= CTR2 STH2M= STR2M ENDIF STH4= STH2*STH2 STH6= STH4*STH2 CTH4= CTH2*CTH2 CSTH2= CTH2*STH2 CMSTH2= CTH2-STH2 TSCTH2= 2.D0/CTH2 FSCTH2= 4.D0/CTH2 CTTH2= CTH2/STH2 FACT= 1.D0/CSTH2 SIGB= CMSTH2*(S3GFZ+S3GBZ+STH2*ZM2*(PGGFZ+PGGLQZ+ # PGGBZ))/ZM2 SIGBM= (SWW0B+S3GBZ-S33BZ)/ZM2+CSTH2*PIBZ AUX0= (SIGFP+RSCEFF+DRHO3)/CMSTH2 AUX1= (SIGFP+RSCEFF+DRHO3)/CMSTH2 AUX2= 3.D0/16.D0/CMSTH2*(16.D0/3.D0*AUX0*AUX0- # 1.D0/CSTH2*HOC*XI2) S33Z= S33FZ+S33BZ S3GZ= S3GFZ+S3GBZ PZ= PIFZ+PIBZ * *-----COMPUTES THE Z0 WAVE-FUNCTION FACTOR * DPIZ= (SWW0B-S33BZ-S3GFZ-ZM2*SP33Z+ # 2.D0*STH2*(S3GZ+ZM2*SP3GZ)+STH4*ZM2*ZM2* # (PPGGZ+PPGGNPZ))/ZM2 * WZZ= DPIZ/CSTH2 * *-----COMPUTES THE RESIDUAL CORRECTIONS TO SIN^2 * DST2= SIGB/CMSTH2 * *-----COMUPTES THE TOTAL Z0 WIDTH * *-----AUXILIARY FUNCTIONS FOR VERTICES * N1= 1 N2= 2 FF1= QRF(N1,QEM,QEM2,QZM2,QEM2) FF2= QRF(N1,QEM,QNM2,QWM2,QNM2) FF3= QRF(N2,QEM,QWM2,QNM2,QWM2) FF3P= 2.D0*FF3+DELW * *-----QCD & MASS CORRECTIONS * CALL CORRQCD(ZM) DO JC=1,9 * IF(OU1.EQ.'N') THEN AX= AEXP ELSE IF(OU1.EQ.'Y') THEN AX= AEXPH ENDIF * *-----LEPTONIC WIDTHS * IF(JC.LE.3) THEN IF(JC.EQ.1) THEN FM= EM ACM= 0.D0 ELSE IF(JC.EQ.2) THEN FM= MM ACM= ACMM ELSE IF(JC.EQ.3) THEN FM= TLM ACM= ACMT ENDIF * FM2= FM*FM VHL= ZID-2.D0*BQL*STH2 VHL2= VHL*VHL VRL= ZID-2.D0*BQL*STR2 VRL2= VRL*VRL CALL WFF(BQL,ZID,WVL,WAL) FRAD= 1.D0+ALPHAH/PI*(3.D0/4.D0) * *-----VERTEX CORRECTIONS * FVL= WVL*VHL+WAL*ZID-2.D0*VHL*(VHL2+3.D0/4.D0)* # FF1+4.D0*ZID*CTH2*FF2-CTH4*ZID*FF3P * FAL= WVL*ZID+WAL*VHL-2.D0*ZID*(3.D0*VHL2+0.25D0)* # FF1+4.D0*ZID*CTH2*FF2-CTH4*ZID*FF3P * *-----EFFECTIVE ASYMMETRIES * COMPUTED ARE ALSO ASYMMETRIES ACCORDING TO THE CONVENTIONAL * DEFINITION WHERE POWERS OF G_V AND G_A ARE NEVER EXPANDED * IF(JC.EQ.1) THEN DGV= 2.D0*FVL/CSTH2+2.D0*DST2-0.5D0*VHL*WZZ DGA= 2.D0*FAL/CSTH2+0.25D0*WZZ GV= VRL+AX*DGV GA= -0.5D0+AX*DGA GV2= VRL2+2.D0*AX*VHL*DGV GA2= 0.25D0-AX*DGA IF(OU2.EQ.'N') THEN AFBN= 0.25D0*VRL2+AX*(-VHL2*DGA+0.5D0*VHL*DGV) AFBD= (VRL2+0.25D0)*(VRL2+0.25D0)+2.D0*AX* # (VHL2+0.25D0)*(2.D0*VHL*DGV-DGA) AFBEFF= 3.D0*AFBN/AFBD ALREFFN= -0.5D0*VRL+AX*(-0.5D0*DGV+VHL*DGA) ALREFFD= VRL2+0.25D0+AX*(2.D0*VHL*DGV-DGA) ALREFF= 2.D0*ALREFFN/ALREFFD ELSE IF(OU2.EQ.'Y') THEN GVGA= GV/GA GVGA2= GVGA*GVGA AFBEFF= 3.D0*GVGA2/(1.D0+GVGA2)**2 ALREFF= 2.D0*GVGA/(1.D0+GVGA2) ENDIF IF(OU2.EQ.'N') THEN ST2E= STR2+0.5D0*AX*(DGV+(4.D0*STR2-1.D0)*DGA) ELSE IF(OU2.EQ.'Y') THEN ST2E= (STR2+0.5D0*AX*(DGV-DGA))/(1.D0-2.D0*AX*DGA) ENDIF ENDIF IF(OU2.EQ.'N') THEN ZW(JC)= GF*RHO*ZM3/6.D0/PI*((VRL2+0.25D0)*FRAD+ # GA2*ACM+AX*(2.D0*FACT*(2.D0*VHL*FVL-FAL)+ # 4.D0*VHL*DST2-(VHL2+0.25D0)*WZZ)) * ELSE IF(OU2.EQ.'Y') THEN ZW(JC)= GF*RHO*ZM3/6.D0/PI*(GV*GV+GA*GA*(1.D0+ACM))* # FRAD ENDIF * *-----HADRONIC WIDTHS * ELSE IF(JC.GT.3.AND.JC.LT.8) THEN SINGF= -0.5D0-2.D0/3.D0*STR2 VUP= 0.5D0-4.D0/3.D0*STR2 VDO= -0.5D0+2.D0/3.D0*STR2 ALS3= (ALSR/PI)**3 IF(JC.EQ.4) THEN FM= UQM FMD= DQM VCM= 0.D0 ACM= 0.D0 ACORQ= ACORU VCORQC= -0.41318D0*ALS3*SINGF/VUP ZI3= ZIU BQCH= BQUQ BQCHD= BQDQ ELSE IF(JC.EQ.5) THEN FM= DQM FMD= UQM VCM= 0.D0 ACM= 0.D0 ACORQ= ACORD VCORQC= -0.41318D0*ALS3*SINGF/VDO ZI3= ZID BQCH= BQDQ BQCHD= BQUQ ELSE IF(JC.EQ.6) THEN FM= SQM FMD= CQM VCM= 0.D0 ACM= 0.D0 ACORQ= ACORD VCORQC= -0.41318D0*ALS3*SINGF/VDO ZI3= ZID BQCH= BQDQ BQCHD= BQUQ ELSE IF(JC.EQ.7) THEN FM= CQM FMD= SQM VCM= VCMC ACM= ACMC ACORQ= ACORU VCORQC= -0.41318D0*ALS3*SINGF/VUP ZI3= ZIU BQCH= BQUQ BQCHD= BQDQ ENDIF FM2= FM*FM FMD2= FMD*FMD BQCH2= BQCH*BQCH FRADQ= 1.D0+ALPHAH/PI*BQCH2*(3.D0/4.D0-1.D0/4.D0*ALSTZ/PI) * *-----UP/DOWN QUARKS (NOT B) * VHQ= ZI3-2.D0*BQCH*STH2 VHQ2= VHQ*VHQ VQT= -ZI3-BQCHD*STH2 VRQ= ZI3-2.D0*BQCH*STR2 VRQ2= VRQ*VRQ VRQT= -ZI3-2.D0*BQCHD*STR2 CALL WFF(BQCH,ZI3,WVQ,WAQ) * *-----VERTEX CORRECTIONS * FVQ= WVQ*VHQ+WAQ*ZI3-2.D0*VHQ*(VHQ2+3.D0/4.D0)* # FF1-4.D0*VQT*CTH2*FF2-CTH4*ZI3*FF3P * FAQ= WVQ*ZI3+WAQ*VHQ-2.D0*ZI3*(3.D0*VHQ2+0.25D0)* # FF1-4.D0*VQT*CTH2*FF2-CTH4*ZI3*FF3P * ACMQ0= -6.D0*FM2/ZM2 DGVQ= 2.D0*FVQ/CSTH2-2.D0*BQCH*DST2-0.5D0*VHQ*WZZ DGAQ= 2.D0*FAQ/CSTH2-0.5*ZI3*WZZ GVQ0= VRQ GAQ0= ZI3 GVQC= AX*DGVQ GAQC= AX*DGAQ GVQ= GVQ0+GVQC GAQ= GAQ0+GAQC GVQS= GVQ*GVQ CMX1= -1.57D0/16.D0/STR2/CTR2 CMX2= -1.68D0/8.D0/STR2 CMX3= -0.93D0*3.D0/2.D0*CTR2/STR2 FMIX= ALPHAH*ALSTZ/PIS CMIX= (1.D0/8.D0+VRQ2*(3.D0+2.D0*VRQ2))*CMX1+ # (VRQ+ZI3)*(VRQT-ZI3)*CMX2+ZI3*(VRQ+ZI3)*CMX3 GQQFMIX= GF*RHO*ZM3/2.D0/PI*FMIX*CMIX * IF(OU6.EQ.'Y') THEN IF(OU2.EQ.'Y') THEN GQQF= GF*RHO*ZM3/2.D0/PI*(GVQ*GVQ*(FRADQ+VCORQ+QCDND+ # VCM)+GAQ*GAQ*(FRADQ+ACORQ+QCDND+ACM)) ELSE IF(OU2.EQ.'N') THEN GQQF= GF*RHO*ZM3/2.D0/PI*((VRQ2+2.D0*VHQ*GVQC)* # (FRADQ+VCORQ+QCDND+VCM)+GAQ0*(GAQ0+2.D0*GAQC)* # (FRADQ+ACORQ+QCDND+ACM)+GAQC*GAQC*ACMQ0) ENDIF GQQF= GQQF+GQQFMIX ZW(JC)= GQQF ELSE IF(OU6.EQ.'N') THEN IF(OU2.EQ.'N') THEN GQQQCD= GF*RHO*ZM3/2.D0/PI*(VRQ2*(VCORQ+QCDND+VCM)+0.25D0* # (ACORQ+QCDND+ACM)+2.D0*GAQ0*GAQC*ACMQ0) GQQEW= GF*RHO*ZM3/2.D0/PI*((VRQ2+0.25D0)*FRADQ+ # AX*(4.D0*FACT*(VHQ*FVQ+ZI3*FAQ)- # 4.D0*BQCH*VHQ*DST2-(VHQ2+0.25D0)*WZZ)) ELSE IF(OU2.EQ.'Y') THEN GQQQCD= GF*RHO*ZM3/2.D0/PI*(VRQ2*(VCORQ+QCDND+VCM)+ # 0.25D0*(ACORQ+QCDND+ACM)+(2.D0*GAQ0*GAQC+ # GAQC*GAQC)*ACMQ0) GQQEW= GF*RHO*ZM3/2.D0/PI*(GVQ*GVQ+GAQ*GAQ)*FRADQ ENDIF ZW(JC)= GQQEW+GQQQCD+GQQFMIX ENDIF ZWQCDS(JC)= GF*RHO*ZM3/2.D0/PI*VRQ2*VCORQC * *-----EFFECTIVE DECONVOLUTED C ASYMMETRY * IF(JC.EQ.7) THEN GVQ00= VHQ GVL0= VRL GVL00= VHL GVLC= AX*DGV GAL0= -0.5D0 GALC= AX*DGA BETC2= 1.D0-4.D0*CQM2/ZM2 BETC= SQRT(BETC2) BETCP= 0.5D0*(3.D0-BETC2) IF(OU2.EQ.'N') THEN AFBCN= GVL0*GAL0*GVQ0*GAQ0+GVLC*GAL0*GVQ00*GAQ0+ # GVL00*GALC*GVQ00*GAQ0+GVL00*GAL0*GVQC*GAQ0+ # GVL00*GAL0*GVQ00*GAQC AFBCN= 3.D0*BETC*AFBCN AFBCD= (GVL0*GVL0+GAL0*GAL0)*(BETCP*GVQ0*GVQ0+BETC2* # GAQ0*GAQ0)+2.D0*(GVL00*GVL00+GAL0*GAL0)*(BETCP* # GVQ00*GVQC+BETC2*GAQ0*GAQC)+(BETCP*GVQ00*GVQ00+ # BETC2*GAQ0*GAQ0)*2.D0*(GVL00*GVLC+GAL0*GALC) AFBCEFF= AFBCN/AFBCD ALRCN= GVQ0*GAQ0+GVQ00*GAQC+GAQ0*GVQC ALRCD= GVQ0*GVQ0+GAQ0*GAQ0+2.D0*(GVQ00*GVQC+GAQ0*GAQC) ALRCEFF= 2.D0*ALRCN/ALRCD ELSE IF(OU2.EQ.'Y') THEN XVL= GVL0+GVLC XAL= GAL0+GALC ETAL= XVL*XAL/(XVL*XVL+XAL*XAL) XVC= GVQ0+GVQC XAC= GAQ0+GAQC ETAC= XVC*XAC/(BETCP*XVC*XVC+BETC2*XAC*XAC) AFBCEFF= 3.D0*BETC*ETAL*ETAC ALRCEFF= 2.D0*XVC*XAC/(XVC*XVC+XAC*XAC) ENDIF ENDIF * *-----B QUARK * ELSE IF(JC.EQ.8) THEN * VHB= -0.5D0+2.D0/3.D0*STH2 VHB2= VHB*VHB VBT= -VHB VRB= -0.5D0+2.D0/3.D0*STR2 VRB2= VRB*VRB FRADB= 1.D0+ALPHAH/PI/9.D0*(3.D0/4.D0-1.D0/4.D0*ALSTZ/PI) * *-----VERTICES * DTW= TQM2-WM2 * B1= -0.5D0*(DELW+0.5D0)+0.5D0*TQM2/DTW*( # TQM2/DTW*LOG(TQM2/WM2)-1.D0) WCOM= 0.25D0*(1.D0+0.5D0*TQM2/WM2)*B1+1.D0/8.D0 WVB= -1.D0/32.D0*(4.D0*VHB2+1.D0)*(DELZ-0.5D0)+ # CTH2*WCOM * WAB= 1.D0/8.D0*VHB*(DELZ-0.5D0)+CTH2*WCOM * N3= 3 N4= 4 N5= 5 N6= 6 * FF4= QRF(N3,QBQM,QTQM2,QWM2,QTQM2) FF5= QRF(N4,QBQM,QTQM2,QWM2,QTQM2) FF6= QRF(N2,QBQM,QWM2,QTQM2,QWM2) FF7= QRF(N5,QBQM,QWM2,QTQM2,QWM2) FF8= QRF(N6,QBQM,QWM2,QTQM2,QWM2) FF9= QRF(N1,QBQM,QBQM2,QZM2,QBQM2) * TFAC= 4.D0*VHB*FF4-1.D0/3.D0*STH2*FF5+CTH2*FF6- # 0.25D0*CMSTH2*FF7+0.5D0*STH2*FF8+0.5D0*CTH2*DELW * FVB= VHB*WVB-0.5D0*WAB+(1.D0-2.D0*STH2+4.D0/3.D0*STH4- # 16.D0/27.D0*STH6)*FF9+CTH2*TFAC * FAB= -0.5D0*WVB+VHB*WAB+(1.D0-2.D0*STH2+4.D0/3.D0*STH4)* # FF9+CTH2*TFAC * IF(OU1.EQ.'N') THEN XIP= XI+(8.D0/3.D0*CTH2+1.D0/6.D0)*LOG(TQM2/WM2) ELSE IF(OU1.EQ.'Y') THEN XIP= 0.D0 ENDIF * * *-----CORRECTIONS O(G_F^2*M_T^4) * XVAR= GFC*XI/4.D0 XVAR2= XVAR*XVAR GBB4= GF*RHO*ZM3/2.D0/PI*XVAR2*((2.D0*VHB-1.D0)*TCOR+2.D0) SINGF= -0.5D0-2.D0/3.D0*STR2 VDO= -0.5D0+2.D0/3.D0*STR2 VCORQC= -0.41318D0*(ALSR/PI)**3*SINGF/VDO IF(OU4.EQ.'Y') THEN ACMB0= -6.D0*RBM2/ZM2 BET2= 1.D0-4.D0*RBM2/ZM2 ELSE IF(OU4.EQ.'N') THEN ACMB0= -6.D0*BQM2/ZM2 BET2= 1.D0-4.D0*BQM2/ZM2 ENDIF DQCDBV= 1.D0+VCORQ+QCDND+VCMB DQCDBA= 1.D0+ACORB+QCDND+ACMB DQCDBA0= 1.D0+ACMB0 * *-----CORRECTIONS O(ALS*G_F*M_T^2) * HVAR= GF/8.D0/PIS*TQM2*(ALSTZ/PI-PIS/3.D0*ALST/PI) HVARS= -GF/24.D0*TQM2*ALST/PI HVARSS= GF/8.D0/PIS*TQM2*ALSTZ/PI IF(OU2.EQ.'N') THEN GBBMIX= GF*RHO*ZM3/2.D0/PI*HVAR*(2.D0*VHB-1.D0) ELSE IF(OU2.EQ.'Y') THEN GBBMIX= GF*RHO*ZM3/2.D0/PI*HVARSS*(2.D0*VHB-1.D0) ENDIF * DGVB= 2.D0*FVB/CSTH2+2.D0/3.D0*DST2-0.5D0*VHB*WZZ DGAB= 2.D0*FAB/CSTH2+0.25D0*WZZ GVB0= VRB GAB0= -0.5D0 IF(OU1.EQ.'Y') THEN GVBC= GFC*(2.D0*FVB+CSTH2*(2.D0/3.D0*DST2-0.5D0*VHB*WZZ))+ # HVARS+XVAR2*TCOR GABC= GFC*(2.D0*FAB+0.25D0*CSTH2*WZZ)+HVARS+XVAR2*TCOR GVB= GVB0+GVBC GAB= GAB0+GABC ELSE IF(OU1.EQ.'N') THEN GVBC= 0.25D0*GFC*XIP+HVARS+XVAR2*TCOR+ # AEXP*(DGVB-0.25D0/CSTH2*XIP) GABC= 0.25D0*GFC*XIP+HVARS+XVAR2*TCOR+ # AEXP*(DGAB-0.25D0/CSTH2*XIP) GVB= GVB0+GVBC GAB= GAB0+GABC ENDIF * IF(OU6.EQ.'Y') THEN IF(OU2.EQ.'Y') THEN GBBF= GF*RHO*ZM3/2.D0/PI*(GVB*GVB*(FRADB+VCORQ+QCDND+VCMB)+ # GAB*GAB*(FRADB+ACORB+QCDND+ACMB)) ELSE IF(OU2.EQ.'N') THEN GBBF= GF*RHO*ZM3/2.D0/PI*((VRB2+2.D0*VHB*GVBC)* # (FRADB+VCORQ+QCDND+VCMB)+GAB0*(GAB0+2.D0*GABC)* # (FRADB+ACORB+QCDND+ACMB)+GABC*GABC*ACMB0) ENDIF ZW(JC)= GBBF ELSE IF(OU6.EQ.'N') THEN * IF(OU2.EQ.'Y') THEN GBBQCD= GF*RHO*ZM3/2.D0/PI*(VRB2*(VCORQ+QCDND+VCMB)+ # 0.25D0*(ACORB+QCDND+ACMB)+(2.D0*GAB0*GABC+ # GABC*GABC)*ACMB0) ELSE IF(OU2.EQ.'N') THEN GBBQCD= GF*RHO*ZM3/2.D0/PI*(VRB2*(VCORQ+QCDND+VCMB)+0.25D0* # (ACORB+QCDND+ACMB)+(2.D0*GAB0*GABC+XVAR2)*ACMB0) ENDIF * IF(OU1.EQ.'N'.AND.OU2.EQ.'N') THEN GBBEW= GF*RHO*ZM3/2.D0/PI*((VRB2+0.25D0)*FRADB+ # 0.25D0*GFC*XIP*(2.D0*VHB-1.D0)+AEXP*FACT* # (2.D0*VHB*(2.D0*FVB-0.25D0*XIP)-(2.D0*FAB- # 0.25D0*XIP))+AEXP*(4.D0/3.D0*VHB*DST2-(VHB2+ # 0.25D0)*WZZ)) ZW(JC)= GBBEW+GBBQCD+GBB4+GBBMIX ELSE IF(OU1.EQ.'Y'.AND.OU2.EQ.'N') THEN GBBEW= GF*RHO*ZM3/2.D0/PI*((VRB2+0.25D0)*FRADB+ # AEXPH*(2.D0*FACT*(2.D0*VHB*FVB-FAB)+4.D0/3.D0* # VHB*DST2-(VHB2+0.25D0)*WZZ)) ZW(JC)= GBBEW+GBBQCD+GBB4+GBBMIX ELSE IF(OU1.EQ.'Y'.AND.OU2.EQ.'Y') THEN ZW(JC)= GF*RHO*ZM3/2.D0/PI*FRADB*(GVB*GVB+GAB*GAB)+ # GBBQCD+GBBMIX ELSE IF(OU1.EQ.'N'.AND.OU2.EQ.'Y') THEN ZW(JC)= GF*RHO*ZM3/2.D0/PI*FRADB*(GVB*GVB+GAB*GAB)+ # GBBQCD+GBBMIX ENDIF ENDIF ZWQCDS(JC)= GF*RHO*ZM3/2.D0/PI*VRB2*VCORQC * *-----EFFECTIVE DECONVOLUTED B ASYMMETRY * GVB0= VRB GVB00= VHB GVL0= VRL GVL00= VHL GVLC= AX*DGV GAL0= -0.5D0 GALC= AX*DGA BET= SQRT(BET2) BETP= 0.5D0*(3.D0-BET2) IF(OU2.EQ.'N') THEN AFBBN= GVL0*GAL0*GVB0*GAB0+GVLC*GAL0*GVB00*GAB0+ # GVL00*GALC*GVB00*GAB0+GVL00*GAL0*GVBC*GAB0+ # GVL00*GAL0*GVB00*GABC AFBBN= 3.D0*BET*AFBBN AFBBD= (GVL0*GVL0+GAL0*GAL0)*(BETP*GVB0*GVB0+BET2*GAB0* # GAB0)+2.D0*(GVL00*GVL00+GAL0*GAL0)*(BETP*GVB00* # GVBC+BET2*GAB0*GABC)+(BETP*GVB00*GVB00+BET2* # GAB0*GAB0)*2.D0*(GVL00*GVLC+GAL0*GALC) AFBBEFF= AFBBN/AFBBD ABPOL= 2.D0*(-0.5D0*GVB0+GVB00*GABC-0.5D0*GVBC)/ # (GVB0*GVB0+0.25D0+2.D0*GVB00*GVBC-GABC) ALRBEFF= ABPOL ELSE IF(OU2.EQ.'Y') THEN XVL= GVL0+GVLC XAL= GAL0+GALC ETAL= XVL*XAL/(XVL*XVL+XAL*XAL) XVB= GVB0+GVBC XAB= GAB0+GABC ETAB= XVB*XAB/(BETP*XVB*XVB+BET2*XAB*XAB) AFBBEFF= 3.D0*BET*ETAL*ETAB ABPOL= 2.D0*GVB*GAB/(GVB*GVB+GAB*GAB) ALRBEFF= ABPOL ENDIF * IF(OU2.EQ.'N') THEN ST2B= STR2-1.5D0*((1.D0-4.D0/3.D0*STR2)*GABC-GVBC)+ # 4.D0*XVAR2*STR2 ELSE IF(OU2.EQ.'Y') THEN ST2B= (STR2+1.5D0*(GVBC-GABC))/(1.D0-2.D0*GABC) ENDIF * *-----INVISIBLE WIDTH * ELSE IF(JC.EQ.9) THEN * VNT= ZID-BQL*STH2 CALL WFF(BQN,ZIU,WVN,WAN) * *-----VERTEX CORRECTIONS * FF10= QRF(N1,QNM,QNM2,QZM2,QNM2) FF11= QRF(N1,QNM,QEM2,QWM2,QEM2) FF12= QRF(N2,QNM,QWM2,QEM2,QWM2) * FN= (WVN+WAN)*ZIU-2.D0*ZIU*FF10-4.D0*VNT*CTH2*FF11-2.D0*CTH4*ZIU* # FF12-CTH4*ZIU*DELW * ZW(JC)= GF*RHO*ZM3/3.D0/PI*(0.25D0+AX*(2.D0*FACT*FN-WZZ/4.D0)) * ENDIF ENDDO * WL= (ZW(1)+ZW(2)+ZW(3))/3.D0 WH= ZW(4)+ZW(7)+ZW(5)+ZW(6)+ZW(8)+ # ZWQCDS(4)+ZWQCDS(7)+ZWQCDS(5)+ZWQCDS(6)+ZWQCDS(8) WT= ZW(1)+ZW(2)+ZW(3)+3.D0*ZW(9)+WH WL= (ZW(1)+ZW(2)+ZW(3))/3.D0 WB= ZW(8) WC= ZW(7) PIC= PIS*PI SHAD= 12.D0*PI/ZM2*ZW(1)*WH/WT/WT*0.38937966D6 RH= WH/ZW(1) RD= ZW(5)/ZW(1) RB= ZW(8)/ZW(1) * *-----COMPUTES THE W MASS * QP2W= -QWM2 CALL ALALS(QP2W,QDV1,QDA1,QF1,QDV1P,QDA1P,QDV1I) CALL WSELF(QP2W,SWW,PWW) * DCZ= SIGBM/CMSTH2-S3GZ/ZM2-STH2M*PIBZ OWM2= (ZM2*RHO*CTR2M+AEXPH/STR2M*(SWW0-SWW))/ # (1.D0-AEXPH/STR2M*(PWW-DCZ)) OWM= SQRT(OWM2) SIN2W= 1.D0-OWM2/ZM2 * RETURN END * * *-----CORRQCD------------------------------------------------------------- *-----COMPUTES ALPHA_S(RUNNING) AND THE VARIOUS QCD CORRECTIONS * SUBROUTINE CORRQCD(SCAL) IMPLICIT REAL*8(A-H,O-Z) REAL*8 MM,NM,MM2,NM2 * COMMON/FN/NF COMMON/PARAM/PI,PIS,DELTA COMMON/COUPLINGS/ALPHA,GF,ALS,RS0,CONV COMMON/FMASSES/EM,MM,TLM,NM,UQM,DQM,CQM,SQM,BQM COMMON/FMASSES2/EM2,MM2,TLM2,NM2,UQM2,DQM2,CQM2,SQM2,BQM2,TQM2 COMMON/QCDCORR/VCORQ,ACORU,ACORD,ACORB,RBM2,RCM2,VCMB,ACMB,VCMC, # ACMC,ACMT,ALSR,CAQCDB,CAQCDC,CAMB,CAMC,CAMT,ACMM, # QCDND * DATA RZ3/1.20205690315959428540D0/,RA4/0.5174790617D0/ DATA RZ5/1.03692775514336992633D0/ * RZ4= PIS*PIS/90.D0 SCAL2= SCAL*SCAL TQM= SQRT(TQM2) ALSR= RALPHAS(RS0,SCAL,ALS,NF) AEXPS= ALSR/PI RVAR= 0.25D0*SCAL2/TQM2 PSC= RVAR*(-8.D0/135.D0*LOG(4.D0*RVAR)+176.D0/675.D0) * *-----POWER SUPPRESSED CORRECTIONS * A2= 1.40923D0+PSC VCORQ= AEXPS*(1.D0+AEXPS*(A2-12.76706D0*AEXPS)) CSI= SCAL/TQM/2.D0 CSI2= CSI*CSI CSI4= CSI2*CSI2 * AUXLN= LOG(SCAL/TQM) AUXLN2= LOG(TQM2/SCAL2) A3= -15.98773D0+AUXLN2*(-67.D0/18.D0+23.D0/12.D0*AUXLN2) ACORB= AEXPS*(1.D0+AEXPS*(A2-12.76706D0*AEXPS+2.D0*AUXLN- # 3.083D0+0.346D0*CSI2+0.211D0*CSI4+AEXPS*A3)) ACORD= ACORB ACORU= AEXPS*(1.D0+AEXPS*(A2-12.76706D0*AEXPS-2.D0*AUXLN+ # 3.083D0-0.346D0*CSI2-0.211D0*CSI4-AEXPS*A3)) * *-----COMPUTES THE RUNNING BOTTOM MASS * * *-----COMPUTES THE RUNNING CHARM MASS * IF(ALS.EQ.0.D0) THEN RBM= BQM RBM2= BQM2 RCM= CQM RCM2= CQM2 ELSE BLN= LOG(SCAL2/BQM2) CLN= LOG(SCAL2/CQM2) BCLN= LOG(BQM2/CQM2) CFB1= 4.D0/3.D0+BLN CFC1= 4.D0/3.D0+CLN+BCLN CFB2= 12.4D0-16.D0/9.D0+BLN*(341.D0/72.D0+ # 11.D0/24.D0*BLN) CFC2= 13.3D0-16.D0/9.D0+BCLN*(367.D0/72.D0+ # 13.D0/24.D0*BCLN)+CLN*(341.D0/72.D0+ # 11.D0/12.D0*BCLN+11.D0/24.D0*CLN) FN= 5.D0 B05= (11.D0-2.D0/3.D0*FN)/4.D0 B15= (102.D0-38.D0/3.D0*FN)/16.D0 B25= (2857.D0/2.D0-5033.D0/18.D0*FN+325.D0/54.D0* # FN*FN)/64.D0 G05= 1.D0 G15= (202.D0/3.D0-20.D0/9.D0*FN)/16.D0 G25= (1249.D0-(2216.D0/27.D0+160.D0/3.D0*RZ3)*FN- # 140.D0/81.D0*FN*FN)/64.D0 REX5= G05/B05 B05S= B05*B05 B05C= B05S*B05 B15S= B15*B15 * ALSZ= RALPHAS(RS0,SCAL,ALS,NF)/PI ALSB= RALPHAS(RS0,BQM,ALS,NF)/PI ALSZ2= ALSZ*ALSZ ALSB2= ALSB*ALSB BMM= BQM*(1.D0-4.D0/3.D0*ALSB+(16.D0/9.D0-12.4D0)*ALSB2) CFM15= G15/B05-B15*G05/B05S CFM25= G25/B05-B15*G15/B05S-B25*G05/B05S+B15S*G05/B05C ALSDF= ALSZ-ALSB RBM= BMM*(ALSZ/ALSB)**REX5*(1.D0+CFM15* # ALSDF+0.5D0*CFM15*CFM15*ALSDF*ALSDF+0.5D0*CFM25* # (ALSZ2-ALSB2)) RBM2= RBM*RBM * FN= 4.D0 B04= (11.D0-2.D0/3.D0*FN)/4.D0 B14= (102.D0-38.D0/3.D0*FN)/16.D0 B24= (2857.D0/2.D0-5033.D0/18.D0*FN+325.D0/54.D0* # FN*FN)/64.D0 G04= 1.D0 G14= (202.D0/3.D0-20.D0/9.D0*FN)/16.D0 G24= (1249.D0-(2216.D0/27.D0+160.D0/3.D0*RZ3)*FN- # 140.D0/81.D0*FN*FN)/64.D0 REX4= G04/B04 B04S= B04*B04 B04C= B04S*B04 B14S= B14*B14 * ALSC4= RALPHAS(RS0,CQM,ALS,NF)/PI ALSC42= ALSC4*ALSC4 XCM= 4.D0/3.D0*ALSC4+13.3D0*ALSC42 CMM4= CQM/(1.D0+XCM) CFM14= G14/B04-B14*G04/B04S CFM24= G24/B04-B14*G14/B04S-B24*G04/B04S+ # B14S*G04/B04C CMB= CMM4*(ALSB/ALSC4)**REX4*(1.D0+CFM14* # (ALSB-ALSC4)+0.5D0*CFM14*CFM14* # (ALSB-ALSC4)**2+0.5D0*CFM24*(ALSB2-ALSC42)) RCM= CMB*(ALSZ/ALSB)**REX5*(1.D0+CFM15* # (ALSZ-ALSB)+0.5D0*CFM15*CFM15* # (ALSZ-ALSB)**2+0.5D0*CFM25*(ALSZ2-ALSB2)) RCM2= RCM*RCM ENDIF * *-----MASS CORRECTIONS * * AV3= 89893.D0/54.D0-1645.D0/36.D0*PIS+820.D0/27.D0*RZ3- * # 36575.D0/54.D0*RZ5 * AA2= -(-2237.D0/8.D0+47.D0/6.D0*PIS+97.D0*RZ3)/6.D0 * AA2= AA2-3.D0 * AA3= -(-25024465.D0/7776.D0+15515.D0/108.D0*PIS+ * # 27545.D0/12.D0*RZ3+25.D0*RZ4-995.D0*RZ5)/6.D0 * FQCD1= 1.D0+AEXPS*(629.D0/6.D0+AV3*AEXPS)/12.D0 * FQCD2= 1.D0+AEXPS*(11.D0/3.D0+(AA2+LOG(4.D0*CSI2)+ * # AA3*AEXPS)*AEXPS) * FQCD2C= 1.D0+AEXPS*(11.D0/3.D0+(AA2+LOG(4.D0*CSI2)+ * # AA3*AEXPS)*AEXPS) * FQCD2T= AEXPS*AEXPS*(8.D0/81.D0-1.D0/27.D0*LOG(SCAL/TQM)) ** * VCMB= 12.D0*RBM2/SCAL2*AEXPS*FQCD1+RBM2*RBM2/SCAL2/SCAL2* * # (-6.D0+AEXPS*(-22.D0+AEXPS*(139.489D0-3.833D0* * # LOG(RBM2/SCAL2)))) * ACMB= -6.D0*RBM2/SCAL2*FQCD2-10.D0*RBM2/TQM2*FQCD2T+ * # RBM2*RBM2/SCAL2/SCAL2*(6.D0+AEXPS*(10.D0+AEXPS*( * # -217.728D0+26.833D0*LOG(RBM2/SCAL2)))) * VCMC= 12.D0*RCM2/SCAL2*AEXPS*FQCD1 * ACMC= -6.D0*RCM2/SCAL2*FQCD2C-10.D0*RCM2/TQM2*FQCD2T ** * CAMB= -4.D0*RBM2/SCAL2 * CAMC= -4.D0*RCM2/SCAL2 ** * CAQCDB= 4.D0/3.D0*ALSR*RBM/SCAL * CAQCDC= 4.D0/3.D0*ALSR*RCM/SCAL ** * ACMT= -6.D0*TLM2/SCAL2 * ACMM= -6.D0*MM2/SCAL2 * CAMT= -4.D0*TLM2/SCAL2 * * QCDND= (RBM2+RCM2)/SCAL2*AEXPS**3*(-560.D0/9.D0+ * # 140.D0/3.D0*RZ3) * FQCD1= 1.D0+AEXPS*(8.7D0+45.15D0*AEXPS) FQCD2= 1.D0+AEXPS*(11.D0/3.D0+(11.286D0+ # LOG(4.D0*CSI2))*AEXPS) FQCD2C= 1.D0+AEXPS*(11.D0/3.D0+(11.286D0+ # LOG(4.D0*CSI2))*AEXPS) FQCD2T= AEXPS*AEXPS*(8.D0/81.D0-1.D0/27.D0*LOG(SCAL/TQM)) * VCMB= 12.D0*RBM2/SCAL2*AEXPS*FQCD1 ACMB= -6.D0*RBM2/SCAL2*FQCD2-10.D0*RBM2/TQM2*FQCD2T VCMC= 12.D0*RCM2/SCAL2*AEXPS*FQCD1 ACMC= -6.D0*RCM2/SCAL2*FQCD2C-10.D0*RCM2/TQM2*FQCD2T * CAMB= -4.D0*RBM2/SCAL2 CAMC= -4.D0*RCM2/SCAL2 * CAQCDB= 4.D0/3.D0*ALSR*RBM/SCAL CAQCDC= 4.D0/3.D0*ALSR*RCM/SCAL * ACMT= -6.D0*TLM2/SCAL2 ACMM= -6.D0*MM2/SCAL2 CAMT= -4.D0*TLM2/SCAL2 * QCDND= 0.D0 RETURN END * *-----RALPHAS-------------------------------------------------------------- *-----COMPUTES ALPHA_S(RS) FOR NF FLAVORS GIVEN ALS= ALPHA_S(RS0) * REAL*8 FUNCTION RALPHAS(RS0,RS,ALS,NF) IMPLICIT REAL*8(A-H,O-Z) REAL*8 MM,NM * COMMON/PARAM/PI,PIS,DELTA COMMON/FMASSES/EM,MM,TLM,NM,UQM,DQM,CQM,SQM,BQM * *-----LIMITS FOR LAMBDA_5 ARE 1 MEV < LAMBDA_5 < 10 GEV * IF(ALS.EQ.0.D0) THEN RALPHAS= 0.D0 ELSE X1= 0.001D0 X2= 10.0D0 XACC= 1.D-12 QCDL= QCDLAM(NF,ALS,RS0,X1,X2,XACC) IF(RS.LT.BQM) THEN EX1= 2.D0/25.D0 EX2= 963.D0/14375.D0 RAT= BQM/QCDL QCDL= QCDL*RAT**EX1*(2.D0*LOG(RAT))**EX2 NFE= NF-1 ELSE NFE= NF ENDIF * QCDB0= 11.D0-2.D0/3.D0*NFE QCDB1= 102.D0-38.D0/3.D0*NFE QCDB2= 0.5D0*(2857.D0-5033.D0/9.D0*NFE+ # 325.D0/27.D0*NFE*NFE) QCDA= 2.D0*LOG(RS/QCDL) RALPHAS= 4.D0*PI/QCDB0/QCDA*(1.D0-QCDB1/QCDB0**2/QCDA* # LOG(QCDA)+(QCDB1/QCDB0**2/QCDA)**2*((LOG(QCDA)- # 0.5D0)**2+QCDB2*QCDB0/QCDB1**2-5.D0/4.D0)) * ENDIF RETURN END * *-----PQCDLAM----------------------------------------------------------------- *-----AUXILIARY ROUTINES TO PRINT THE FINAL LAMBA^NF_MSBAR * SUBROUTINE PQCDLAM(OALS,OQCDLAM) IMPLICIT REAL*8(A-H,O-Z) * COMMON/FN/NF COMMON/COUPLINGS/ALPHA,GF,ALS,RS0,CONV * X1= 0.01D0 X2= 2.0D0 XACC= 1.D-6 QALS= OALS OQCDLAM= QCDLAM(NF,QALS,RS0,X1,X2,XACC) * RETURN END * *-----QCDLAM----------------------------------------------------------------- *-----COMPUTES LAMBA^NF_MSBAR FROM ALPHA_S(RS0) * REAL*8 FUNCTION QCDLAM(NF,ALS,RS,X1,X2,XACC) IMPLICIT REAL*8(A-H,O-Z) * PARAMETER (JMAX=50,NOUT=21) * FMID= QCDSCALE(NF,ALS,RS,X2) F= QCDSCALE(NF,ALS,RS,X1) IF (F*FMID.GE.0.D0) THEN WRITE(NOUT,FMT=*) 'ROOT MUST BE BRACKETED FOR BISECTION' PRINT*,'ROOT MUST BE BRACKETED FOR BISECTION' WRITE(NOUT,FMT=1) ALS PRINT 1,ALS 1 FORMAT(/' ERROR DETECTED BY QCDLAM ',/ # ' CURRENT VALUE OF ALPHA_S = ',E20.5) STOP ENDIF IF (F.LT.0.D0) THEN QCDLAM= X1 DX= X2-X1 ELSE QCDLAM= X2 DX= X1-X2 ENDIF DO J=1,JMAX DX= DX*0.5D0 XMID= QCDLAM+DX FMID= QCDSCALE(NF,ALS,RS,XMID) IF (FMID.LE.0.D0) QCDLAM= XMID IF(ABS(DX).LT.XACC.OR.FMID.EQ.0.D0) RETURN ENDDO PAUSE 'TOO MANY BISECTIONS' END * *-----QCDSCALE------------------------------------------------------------- *-----COMPUTES LAMBA^NF_MSBAR FROM ALPHA_S(RS0) * REAL*8 FUNCTION QCDSCALE(NF,ALS,RS,X) IMPLICIT REAL*8(A-H,O-Z) * COMMON/PARAM/PI,PIS,DELTA * QCDB0= 11.D0-2.D0/3.D0*NF QCDB1= 102.D0-38.D0/3.D0*NF QCDB2= 0.5D0*(2857.D0-5033.D0/9.D0*NF+325.D0/27.D0*NF*NF) QCDA= 2.D0*LOG(RS/X) QCDSCALE= ALS-(4.D0*PI/QCDB0/QCDA*(1.D0-QCDB1/QCDB0**2/QCDA* # LOG(QCDA)+(QCDB1/QCDB0**2/QCDA)**2*((LOG(QCDA)- # 0.5D0)**2+QCDB2*QCDB0/QCDB1**2-5.D0/4.D0))) * RETURN END * * *-----WFF------------------------------------------------------------------ *-----FERMION WAVE-FUNCTION FACTORS (F NEQ. B, NON E.M.) * SUBROUTINE WFF(QF,ZIF,WV,WA) IMPLICIT REAL*8(A-H,O-Z) * COMMON/PARAM/PI,PIS,DELTA COMMON/PAR/STH2,STH4,CTH2,ZM2,WM2,HM2 * FZ= DELTA-LOG(ZM2)-0.5D0 FW= DELTA-LOG(WM2)-0.5D0 WV= -1.D0/32.D0*(4.D0*(ZIF-2.D0*QF*STH2)**2+1.D0)*FZ # -1.D0/8.D0*CTH2*FW WA= -1.D0/16.D0*(1.D0-8.D0*ZIF*QF*STH2)*FZ # -1.D0/8.D0*CTH2*FW * RETURN END * *-----QRF----------------------------------------------------------------- *-----VERTEX FORM FACTORS (REAL PART) * REAL*16 FUNCTION QRF(N,FM,RM12,RM22,RM32) IMPLICIT REAL*16(A-H,O-Z) * COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DIMENSION C0(2),C1(2,2),C2(2,4) * JFLAG= 1 P12= -FM*FM P22= -FM*FM S= -QZM2 * CALL CFF(JFLAG,P12,P22,S,RM12,RM22,RM32,C0,C1,C2) * IF (N.EQ.1) THEN QRF= -C2(1,4)/4.Q0+QZM2/8.Q0*(C1(1,1)+C2(1,3))+1.Q0/8.Q0 * ELSE IF(N.EQ.2) THEN QRF= -3.Q0*C2(1,4)+QZM2/2.Q0*(C0(1)+C1(1,1)+C2(1,3))+1.Q0/2.Q0 * ELSE IF(N.EQ.3) THEN QRF= -C2(1,4)/4.Q0+QZM2/8.Q0*(C1(1,1)+C2(1,3))+1.Q0/8.Q0- # RM12**2/RM22/16.Q0*C0(1) * ELSE IF(N.EQ.4) THEN QRF= RM12*(C0(1)-4.Q0/RM22*(-C2(1,4)/4.Q0+QZM2/8.Q0*(C1(1,2)+ # C2(1,3))+1.Q0/16.Q0)) * ELSE IF(N.EQ.5) THEN QRF= RM22/RM12*C2(1,4) * ELSE IF(N.EQ.6) THEN QRF= RM22*C0(1) * ENDIF * RETURN END * *-----QF------------------------------------------------------------------ *-----VERTEX FORM FACTORS (REAL AND IMAGINARY PARTS) * REAL*16 FUNCTION QF(N,I,P2,FM,RM12,RM22,RM32) IMPLICIT REAL*16(A-H,O-Z) * DIMENSION C0(2),C1(2,2),C2(2,4) * JFLAG= 1 P12= -FM*FM P22= -FM*FM * CALL CFF(JFLAG,P12,P22,P2,RM12,RM22,RM32,C0,C1,C2) * IF (N.EQ.1) THEN QF= -C2(I,4)/4.Q0-P2/8.Q0*(C1(I,1)+C2(I,3))+(2.Q0-I)/8.Q0 * ELSE IF(N.EQ.2) THEN QF= -3.Q0*C2(I,4)-P2/2.Q0*(C0(I)+C1(I,1)+C2(I,3))+ # (2.Q0-I)/2.Q0 * ELSE IF(N.EQ.3) THEN QF= -1.Q0/6.Q0*RM12/RM22*(2.Q0*C2(I,4)-(2.Q0-I)/2.Q0+ # P2*(C1(I,2)+C2(I,3))+RM12*C0(I)) * ELSE IF(N.EQ.4) THEN QF= RM22/2.Q0*(C0(I)+1.Q0/RM12*C2(I,4)) * ELSE IF (N.EQ.5) THEN QF= C0(I) * ELSE IF (N.EQ.6) THEN QF= -0.25Q0*C2(I,4)-P2/8.Q0*(C1(I,1)+C2(I,3)) # -1.Q0/16.Q0*RM12*RM12/RM22*C0(I)+(2.Q0-I)/8.Q0 * ELSE IF (N.EQ.7) THEN QF= RM12*C0(I)+RM12/RM22*(C2(I,4)+0.5Q0*P2*( # C1(I,2)+C2(I,3))-(2.Q0-I)/4.Q0) * ELSE IF (N.EQ.8) THEN QF= -3.Q0*C2(I,4)-P2/2.Q0*(C0(I)+C1(I,1)+C2(I,3)) # +(2.Q0-I)/2.Q0 * ELSE IF (N.EQ.9) THEN QF= RM22/RM12*C2(I,4) * ELSE IF (N.EQ.10) THEN QF= RM22*C0(I) * ENDIF * RETURN END * * *-----PSELF-------------------------------------------------------- * THE PHOTON TWO-POINT FUNCTION * SUBROUTINE PSELF(P2X,OPGGF,OPGGLQ,OPGGB,OPGGNP,OPPGG,OPPGGNP, # OIGG,OIGZ) IMPLICIT REAL*16(A-H,P-Z) IMPLICIT REAL*8(O) * COMMON/AAEM/OIAEM COMMON/SSCAL/QSPSC COMMON/QVARIA/QALPHA,QSLLC COMMON/MIX/QALST,QALSTZ,QALSTS COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/MIXC/QV1,QA1,QF1,QV1P,QA1P,QV1I COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DATA Z3/1.20205690315959428540Q0/ * AIAEM= OIAEM*1.D15*1.Q-15 CALL RBFF(P2X,QWM2,QWM2,B0WW,B1WW,B21WW) CALL RBFF(P2X,QEM2,QEM2,B0EE,B1EE,B21EE) CALL RBFF(P2X,QMM2,QMM2,B0MM,B1MM,B21MM) CALL RBFF(P2X,QTM2,QTM2,B0TAU,B1TAU,B21TAU) CALL RBFF(P2X,QTQM2,QTQM2,B0TT,B1TT,B21TT) * CALL RBFF(P2X,QUQM2,QUQM2,B0UU,B1UU,B21UU) CALL RBFF(P2X,QDQM2,QDQM2,B0DD,B1DD,B21DD) CALL RBFF(P2X,QCQM2,QCQM2,B0CC,B1CC,B21CC) CALL RBFF(P2X,QSQM2,QSQM2,B0SS,B1SS,B21SS) CALL RBFF(P2X,QBQM2,QBQM2,B0BB,B1BB,B21BB) * CALL RBPFF(P2X,QWM2,QWM2,E0WW,E1WW,E21WW) CALL RBPFF(P2X,QEM2,QEM2,E0EE,E1EE,E21EE) CALL RBPFF(P2X,QMM2,QMM2,E0MM,E1MM,E21MM) CALL RBPFF(P2X,QTM2,QTM2,E0TAU,E1TAU,E21TAU) CALL RBPFF(P2X,QTQM2,QTQM2,E0TT,E1TT,E21TT) * BFE= 2.Q0*B21EE-B0EE BFM= 2.Q0*B21MM-B0MM BFTAU= 2.Q0*B21TAU-B0TAU BFT= 2.Q0*B21TT-B0TT * BFU= 2.Q0*B21UU-B0UU BFD= 2.Q0*B21DD-B0DD BFC= 2.Q0*B21CC-B0CC BFS= 2.Q0*B21SS-B0SS BFB= 2.Q0*B21BB-B0BB BFUQ= BFU+BFC BFDQ= BFD+BFS+BFB * EFE= 2.Q0*E21EE-E0EE EFM= 2.Q0*E21MM-E0MM EFTAU= 2.Q0*E21TAU-E0TAU EFT= 2.Q0*E21TT-E0TT * QS2= QSPSC*QSPSC AEXPS= QALSTS/QPI AEXPSZ= QALSTZ/QPI RM= QTQM2/QS2 RL= LOG(RM) BX= -RL-4.Q0*Z3+55.Q0/12.Q0 RV= -0.25Q0*P2X/QTQM2 FLQI= 0.25Q0*QPI FB= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QS2) FLQ= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QZM2) * OPGGB= # -12.Q0*B21WW+7.Q0*B0WW+2.Q0/3.Q0 * OPGGF= # 4.Q0*(BFE+BFM+BFTAU)*QSLLC+16.Q0/3.Q0*BFT+ # 64.Q0/9.Q0*AEXPS*QTQM2/P2X*(RV*BX+QV1) * OPGGLQ= # 16.Q0/3.Q0*BFUQ+4.Q0/3.Q0*BFDQ- # 16.Q0/9.Q0*AEXPS*FB-160.Q0/9.Q0*AEXPSZ*FLQ * OIGG= 64.Q0/9.Q0*AEXPS*QTQM2/P2X*QV1I OIGZ= 8.Q0/3.Q0*AEXPS*QTQM2*QV1I-4.Q0/3.Q0*AEXPS* # P2X*FLQI-8.Q0*AEXPSZ*P2X*FLQI * OPPGG= # -(-12.Q0*E21WW+7.Q0*E0WW+4.Q0*(EFE+EFM+EFTAU)* # QSLLC+16.Q0/3.Q0*EFT)+4.Q0/9.Q0*AEXPS/QTQM2/ # RV/RV*(RV*QV1P-QV1) * P2Z= -QZM2 CALL RBFF(P2Z,QEM2,QEM2,B0EEZ,B1EEZ,B21EEZ) CALL RBFF(P2Z,QMM2,QMM2,B0MMZ,B1MMZ,B21MMZ) CALL RBFF(P2Z,QTM2,QTM2,B0TAUZ,B1TAUZ,B21TAUZ) CALL RBFF0(QEM2,QEM2,B0EE0,B1EE0,B21EE0) CALL RBFF0(QMM2,QMM2,B0MM0,B1MM0,B21MM0) CALL RBFF0(QTM2,QTM2,B0TAU0,B1TAU0,B21TAU0) BFEZ= 2.Q0*B21EEZ-B0EEZ BFMZ= 2.Q0*B21MMZ-B0MMZ BFTAUZ= 2.Q0*B21TAUZ-B0TAUZ BFE0= 2.Q0*B21EE0-B0EE0 BFM0= 2.Q0*B21MM0-B0MM0 BFTAU0= 2.Q0*B21TAU0-B0TAU0 PGGL= 4.Q0*(BFEZ+BFMZ+BFTAUZ)*QSLLC PGGL0= 4.Q0*(BFE0+BFM0+BFTAU0)*QSLLC B= 0.00299Q0 C= 1.Q0 AN= -B*LOG(1.Q0+C*QZM2)-QALPHA/4.Q0/QPI*(PGGL-PGGL0)+ # 1.Q0-AIAEM*QALPHA DA= AN-0.00165Q0 AP2X= ABS(P2X) PX= SQRT(AP2X) IF(PX.LT.0.3Q0) THEN A= 0.Q0 B= 0.00835Q0 C= 1.Q0 ELSE IF(PX.LT.3.Q0) THEN A= 0.Q0 B= 0.00238Q0 C= 3.927Q0 ELSE IF(PX.LT.100.Q0) THEN A= 0.00165Q0 B= 0.00299Q0 C= 1.Q0 ELSE IF(PX.GT.100.Q0) THEN A= 0.00221Q0 B= 0.00293Q0 C= 1.Q0 ENDIF * A= A+DA OPGGNP= 4.Q0*QPI/QALPHA*(A+B*LOG(1.Q0+C*AP2X)) OPPGGNP= -4.Q0*QPI/QALPHA*B/(1.Q0+QZM2) * RETURN END * *-----VBSELF0----------------------------------------------------------- * VECTOR-BOSONS TWO-POINT FUNCTIONS AT ZERO MOMENTUM * SUBROUTINE VBSELF0(OPGGF0,OPGGB0,OSWW0F,OSWW0B) IMPLICIT REAL*16(A-H,P-Z) IMPLICIT REAL*8(O) CHARACTER*1 OU1,OU2,OU3,OU4,OU5,OU6,OU7 * COMMON/SSCAL/QSPSC COMMON/SUB/OSINT2,ODSWW COMMON/QVARIA/QALPHA,QSLLC COMMON/THU/OU1,OU2,OU3,OU4,OU5,OU6,OU7 COMMON/MIX/QALST,QALSTZ,QALSTS COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DATA Z3/1.20205690315959428540Q0/ * PM= 1.Q-10 PM2= PM*PM * CALL RBFF0(QWM2,QWM2,B0WW0,B1WW0,B21WW0) CALL RBFF0(PM2,QWM2,B0PW0,B1PW0,B21PW0) CALL RBFF0(QEM2,QEM2,B0EE0,B1EE0,B21EE0) CALL RBFF0(QMM2,QMM2,B0MM0,B1MM0,B21MM0) CALL RBFF0(QTM2,QTM2,B0TAU0,B1TAU0,B21TAU0) CALL RBFF0(QTQM2,QTQM2,B0TT0,B1TT0,B21TT0) CALL RBFF0(QZM2,QWM2,B0ZW0,B1ZW0,B21ZW0) CALL RBFF0(QWM2,QHM2,B0WH0,B1WH0,B21WH0) * BFE= 2.Q0*B21EE0-B0EE0 BFM= 2.Q0*B21MM0-B0MM0 BFTAU= 2.Q0*B21TAU0-B0TAU0 BFT= 2.Q0*B21TT0-B0TT0 * QS2= QSPSC*QSPSC AEXPS= QALSTS/QPI RM= QTQM2/QS2 RL= LOG(RM) BX= -RL-4.Q0*Z3+55.Q0/12.Q0 BY= 3.Q0*RL*RL-11.Q0/2.Q0*RL+6.Q0*Z3+QPIS/2.Q0- # 11.Q0/8.Q0 F10= -3.Q0/2.Q0*Z3-QPIS/12.Q0+23.Q0/16.Q0 * OPGGB0= # -12.Q0*B21WW0+7.Q0*B0WW0+2.Q0/3.Q0 * OPGGF0= # 4.Q0*(BFE+BFM+BFTAU)*QSLLC+16.Q0/3.Q0*BFT- # 16.Q0/9.Q0*AEXPS*(BX+4.Q0*Z3-5.Q0/6.Q0) * CALL RBFF0(QNM2,QEM2,B0NE0,B1NE0,B21NE0) CALL RBFF0(QNM2,QMM2,B0NM0,B1NM0,B21NM0) CALL RBFF0(QNM2,QTM2,B0NT0,B1NT0,B21NT0) CALL RBFF0(QUQM2,QDQM2,B0UD0,B1UD0,B21UD0) CALL RBFF0(QCQM2,QSQM2,B0CS0,B1CS0,B21CS0) CALL RBFF0(QTQM2,QBQM2,B0TB0,B1TB0,B21TB0) * * OSWW0F= # QEM2*B1NE0+3.Q0*(QDQM2-QUQM2)*B1UD0-3.Q0*QUQM2* # B0UD0+QMM2*B1NM0+3.Q0*(QSQM2-QCQM2)*B1CS0- # 3.Q0*QCQM2*B0CS0+QTM2*B1NT0+3.Q0*QBQM2*B1TB0- # 3.Q0*QTQM2*(B0TB0+B1TB0)+4.Q0*AEXPS*QTQM2*( # 0.25D0*BY+F10) * OS0WW0B= # 9.Q0/2.Q0*(QZM2-QWM2)*B1ZW0+0.25Q0*(13.Q0*QZM2- # 21.Q0*QWM2)*B0ZW0+4.Q0*QWM2*B0WW0 OSH= 0.5Q0*(QWM2-QHM2)*B1WH0+0.25Q0*(5.Q0*QWM2-QHM2)* # B0WH0 IF(OU3.EQ.'N') THEN OS0WW0B= OS0WW0B+OSH ELSE IF(OU3.EQ.'Y') THEN OSWW0F= OSWW0F+OSH ENDIF * OS1WW0= # 2.Q0*(QWM2-QZM2)*(2.Q0*B1ZW0+B0ZW0)-2.Q0*QWM2* # (2.Q0*B1PW0+B0PW0) * OSWW0B= OS0WW0B+OSINT2*OS1WW0+ODSWW * RETURN END * *-----RBFF0----------------------------------------------------------------- * THE SCALAR FORM FACTORS B0,B1,B21 AT ZERO MOMENTUM * SUBROUTINE RBFF0(RM12,RM22,B0,B1,B21) IMPLICIT REAL*16(A-H,O-Z) * COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA * DIMENSION ARG(2),CLN(2),FR(3) * IF (RM12.EQ.RM22) THEN FACT= QDELTA-LOG(RM12) B0= FACT B1= -0.5Q0*FACT B21= FACT/3.Q0 RETURN ELSE N= 3 YR= RM12/(RM12-RM22) OMYR= -RM22/(RM12-RM22) YI= -QEPS/(RM12-RM22) CALL RCG(N,YR,YI,OMYR,FR) ARG(1)= OMYR ARG(2)= -YI CALL CQLNX(ARG,CLN) F1R= FR(1)-CLN(1) F2R= 2.Q0*FR(2)-CLN(1) F3R= 3.Q0*FR(3)-CLN(1) FACT= QDELTA-LOG(RM22) B0= FACT-F1R B1= 0.5Q0*(-FACT+F2R) B21= 1.Q0/3.Q0*(FACT-F3R) RETURN ENDIF END * *-----VBSELF------------------------------------------------------------ * Z^0 TWO-POINT FUNCTIONS AT ARBITRARY MOMENTUM * SUBROUTINE VBSELF(P2X,OS3GF,OS33F,OS3GB,OS33B,OSP3G,OSP33) IMPLICIT REAL*16(A-H,P-Z) IMPLICIT REAL*8(O) CHARACTER*1 OU1,OU2,OU3,OU4,OU5,OU6,OU7 * COMMON/SSCAL/QSPSC COMMON/QVARIA/QALPHA,QSLLC COMMON/MIX/QALST,QALSTZ,QALSTS COMMON/THU/OU1,OU2,OU3,OU4,OU5,OU6,OU7 COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/MIXC/QV1,QA1,QF1,QV1P,QA1P,QV1I COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DATA Z3/1.20205690315959428540Q0/ * CALL RBFF(P2X,QWM2,QWM2,B0WW,B1WW,B21WW) CALL RBFF(P2X,QEM2,QEM2,B0EE,B1EE,B21EE) CALL RBFF(P2X,QUQM2,QUQM2,B0UU,B1UU,B21UU) CALL RBFF(P2X,QDQM2,QDQM2,B0DD,B1DD,B21DD) CALL RBFF(P2X,QMM2,QMM2,B0MM,B1MM,B21MM) CALL RBFF(P2X,QCQM2,QCQM2,B0CC,B1CC,B21CC) CALL RBFF(P2X,QSQM2,QSQM2,B0SS,B1SS,B21SS) CALL RBFF(P2X,QTM2,QTM2,B0TAU,B1TAU,B21TAU) CALL RBFF(P2X,QTQM2,QTQM2,B0TT,B1TT,B21TT) CALL RBFF(P2X,QBQM2,QBQM2,B0BB,B1BB,B21BB) CALL RBFF(P2X,QNM2,QNM2,B0NN,B1NN,B21NN) CALL RBFF(P2X,QZM2,QHM2,B0ZH,B1ZH,B21ZH) * CALL RBFF0(QWM2,QWM2,B0WW0,B1WW0,B21WW0) * CALL RBPFF(P2X,QWM2,QWM2,E0WW,E1WW,E21WW) CALL RBPFF(P2X,QEM2,QEM2,E0EE,E1EE,E21EE) CALL RBPFF(P2X,QUQM2,QUQM2,E0UU,E1UU,E21UU) CALL RBPFF(P2X,QDQM2,QDQM2,E0DD,E1DD,E21DD) CALL RBPFF(P2X,QMM2,QMM2,E0MM,E1MM,E21MM) CALL RBPFF(P2X,QCQM2,QCQM2,E0CC,E1CC,E21CC) CALL RBPFF(P2X,QSQM2,QSQM2,E0SS,E1SS,E21SS) CALL RBPFF(P2X,QTM2,QTM2,E0TAU,E1TAU,E21TAU) CALL RBPFF(P2X,QTQM2,QTQM2,E0TT,E1TT,E21TT) CALL RBPFF(P2X,QBQM2,QBQM2,E0BB,E1BB,E21BB) CALL RBPFF(P2X,QNM2,QNM2,E0NN,E1NN,E21NN) CALL RBPFF(P2X,QZM2,QHM2,E0ZH,E1ZH,E21ZH) * BFE= 2.Q0*B21EE-B0EE BFM= 2.Q0*B21MM-B0MM BFTAU= 2.Q0*B21TAU-B0TAU BFU= 2.Q0*B21UU-B0UU BFD= 2.Q0*B21DD-B0DD BFC= 2.Q0*B21CC-B0CC BFS= 2.Q0*B21SS-B0SS BFT= 2.Q0*B21TT-B0TT BFB= 2.Q0*B21BB-B0BB BFN= 2.Q0*B21NN-B0NN BFL= BFE+BFM+BFTAU BFUQ= BFU+BFC+BFT BFDQ= BFD+BFS+BFB * EFE= 2.Q0*E21EE-E0EE EFM= 2.Q0*E21MM-E0MM EFTAU= 2.Q0*E21TAU-E0TAU EFU= 2.Q0*E21UU-E0UU EFD= 2.Q0*E21DD-E0DD EFC= 2.Q0*E21CC-E0CC EFS= 2.Q0*E21SS-E0SS EFT= 2.Q0*E21TT-E0TT EFB= 2.Q0*E21BB-E0BB EFN= 2.Q0*E21NN-E0NN EFL= EFE+EFM+EFTAU EFUQ= EFU+EFC+EFT EFDQ= EFD+EFS+EFB * QS2= QSPSC*QSPSC AEXPS= QALSTS/QPI AEXPSZ= QALSTZ/QPI RM= QTQM2/QS2 RL= LOG(RM) BX= -RL-4.Q0*Z3+55.Q0/12.Q0 BY= 3.Q0*RL*RL-11.Q0/2.Q0*RL+6.Q0*Z3+QPIS/2.Q0- # 11.Q0/8.Q0 RV= -0.25Q0*P2X/QTQM2 XV= -P2X/QTQM2 FB= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QS2) FBP= -Z3+43.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QS2) FLQ= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QZM2) FLQP= -Z3+43.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QZM2) * OS3GF= # P2X*(BFL+2.Q0*BFUQ+BFDQ)+8.Q0/3.Q0*AEXPS*QTQM2*( # RV*BX+QV1)-4.Q0/3.Q0*AEXPS*P2X*FB-8.Q0*AEXPSZ*P2X*FLQ * OS3GB= # P2X*(-10.Q0*B21WW+13.Q0/2.Q0*B0WW+2.Q0/3.Q0) # -2.Q0*QWM2*(B0WW-B0WW0) * PZZF= # +0.5Q0*(BFL+3.Q0*BFUQ+3.Q0*BFDQ+3.Q0*BFN) * PZZB= # -9.Q0*B21WW+25.Q0/4.Q0*B0WW+2.Q0/3.Q0 * PZZH= -B21ZH-B1ZH-0.25Q0*B0ZH IF(OU3.EQ.'N') THEN PZZB= PZZB+PZZH ELSE IF(OU3.EQ.'Y') THEN PZZF= PZZF+PZZH ENDIF * SIGZZF= # -0.5Q0*(QEM2*B0EE+3.Q0*QUQM2*B0UU+3.Q0*QDQM2* # B0DD+QMM2*B0MM+3.Q0*QCQM2*B0CC+3.Q0*QSQM2*B0SS+ # QTM2*B0TAU+3.Q0*QBQM2*B0BB)-3.Q0/2.Q0*QTQM2*B0TT * SIGZZB= # -2.Q0*QWM2*B0WW+4.Q0*QWM2*B0WW0 * SIGZZH= QZM2*(0.5Q0*B1ZH+5.Q0/4.Q0* # B0ZH)-0.5Q0*QHM2*(B1ZH+0.5Q0*B0ZH) IF(OU3.EQ.'N') THEN SIGZZB= SIGZZB+SIGZZH ELSE IF(OU3.EQ.'Y') THEN SIGZZF= SIGZZF+SIGZZH ENDIF * OS33F= P2X*PZZF+SIGZZF+AEXPS*QTQM2*(2.Q0*RV*BX+BY+QV1+QA1)- # 2.Q0*AEXPS*P2X*FB-8.Q0*AEXPSZ*P2X*FLQ * OS33B= P2X*PZZB+SIGZZB * OSP3G= # -10.Q0*B21WW+13.Q0/2.Q0*B0WW+2.Q0/3.Q0+ # BFL+2.Q0*BFUQ+BFDQ+P2X*(10.Q0*E21WW-13.Q0/2.Q0* # E0WW-EFL-2.Q0*EFUQ-EFDQ)+2.Q0*QWM2*E0WW- # AEXPS*(2.Q0/3.Q0*(BX+QV1P)+4.Q0/3.Q0*FBP)- # 8.Q0*AEXPSZ*FLQP * PZZ= PZZF+PZZB * PPZZ= # -9.Q0*E21WW+25.Q0/4.Q0*E0WW # +0.5Q0*(EFL+3.Q0*EFUQ+3.Q0*EFDQ+3.Q0*EFN) # -E21ZH-E1ZH-0.25Q0*E0ZH * SIGPZZ= # -2.Q0*QWM2*E0WW-0.5Q0*(QEM2*E0EE+3.Q0*QUQM2*E0UU # +3.Q0*QDQM2*E0DD+QMM2*E0MM+3.Q0*QCQM2*E0CC+3.Q0* # QSQM2*E0SS+QTM2*E0TAU+3.Q0*QBQM2*E0BB)-3.Q0/2.Q0* # QTQM2*E0TT+QZM2*(0.5Q0*E1ZH+5.Q0/4.Q0*E0ZH)-0.5Q0* # QHM2*(E1ZH+0.5Q0*E0ZH) * OSP33= -P2X*PPZZ+PZZ-SIGPZZ-AEXPS*(0.25Q0*(2.Q0*BX+QV1P+ # QA1P)+2.Q0*FBP)-8.Q0*AEXPSZ*FLQP * RETURN END * *-----WSELF------------------------------------------------------------ * W TWO-POINT FUNCTIONS AT ARBITRARY MOMENTUM * SUBROUTINE WSELF(P2X,OSWW,OPWW) IMPLICIT REAL*16(A-H,P-Z) IMPLICIT REAL*8(O) * COMMON/SSCAL/QSPSC COMMON/SUB/OSINT2,ODSWW COMMON/QVARIA/QALPHA,QSLLC COMMON/MIX/QALST,QALSTZ,QALSTS COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/MIXC/QV1,QA1,QF1,QV1P,QA1P,QV1I COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DATA Z3/1.20205690315959428540Q0/ * PM= 1.Q-10 PM2= PM*PM * CALL RBFF(P2X,QNM2,QEM2,B0NE,B1NE,B21NE) CALL RBFF(P2X,QNM2,QMM2,B0NM,B1NM,B21NM) CALL RBFF(P2X,QNM2,QTM2,B0NT,B1NT,B21NT) CALL RBFF(P2X,QUQM2,QDQM2,B0UD,B1UD,B21UD) CALL RBFF(P2X,QCQM2,QSQM2,B0CS,B1CS,B21CS) CALL RBFF(P2X,QTQM2,QBQM2,B0TB,B1TB,B21TB) CALL RBFF(P2X,QZM2,QWM2,B0ZW,B1ZW,B21ZW) CALL RBFF(P2X,QWM2,QHM2,B0WH,B1WH,B21WH) CALL RBFF(P2X,QPM2,QWM2,B0PW,B1PW,B21PW) * CALL RBFF0(QWM2,QWM2,B0WW0,B1WW0,B21WW0) * QS2= QSPSC*QSPSC AEXPS= QALSTS/QPI AEXPSZ= QALSTZ/QPI RM= QTQM2/QS2 RL= LOG(RM) BX= -RL-4.Q0*Z3+55.Q0/12.Q0 BY= 3.Q0*RL*RL-11.Q0/2.Q0*RL+6.Q0*Z3+QPIS/2.Q0- # 11.Q0/8.Q0 RV= -0.25Q0*P2X/QTQM2 XV= -P2X/QTQM2 FB= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QS2) FLQ= -Z3+55.Q0/48.Q0-0.25Q0*LOG(ABS(P2X)/QZM2) * SIG0WW= # 9.Q0/2.Q0*(QZM2-QWM2)*B1ZW+(13.Q0*QZM2-21.Q0*QWM2)/4.Q0* # B0ZW+0.5Q0*(QWM2-QHM2)*B1WH+0.25Q0*(5.Q0*QWM2-QHM2)*B0WH+ # QEM2*B1NE+QMM2*B1NM+QTM2*B1NT+3.Q0*((QDQM2-QUQM2)*B1UD- # QUQM2*B0UD+(QSQM2-QCQM2)*B1CS-QCQM2*B0CS+(QBQM2-QTQM2)* # B1TB-QTQM2*B0TB)+4.Q0*QWM2*B0WW0 * SIG1WW= # 2.Q0*(QWM2-QZM2)*(2.Q0*B1ZW+B0ZW)-2.Q0*QWM2*(2.Q0*B1PW+ # B0PW) * PI0WW= 2.Q0/3.Q0-9.Q0*B21ZW-9.Q0*B1ZW+7.Q0/4.Q0*B0ZW-B21WH-B1WH- # 0.25Q0*B0WH+2.Q0*(B21NE+B21NM+B21NT+B1NE+B1NM+B1NT+3.Q0*( # B21UD+B21CS+B21TB+B1UD+B1CS+B1TB)) * PI1WW= 8.Q0*B21ZW-2.Q0*B0ZW+8.Q0*B1ZW-8.Q0*B21PW-8.Q0*B1PW+ # 2.Q0*B0PW * OPWW= PI0WW+OSINT2*PI1WW-8.Q0*AEXPSZ*FLQ OSWW= SIG0WW+OSINT2*SIG1WW+4.Q0*AEXPS*QTQM2*(0.25Q0*(XV*BX+BY)+ # QF1) * RETURN END * *-----POLINT-------------------------------------------------------------------- * GIVES THE HIGHER ORDER CORRECTIONS TO RHO WHEN M_H < 2 M_T * (POLINOMIAL INTERPOLATION) * SUBROUTINE POLINT(XA,YA,N,X,Y,DY) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (NMAX=30) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) * NS= 1 DIF= ABS(X-XA(1)) DO 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) ENDDO Y= YA(NS) NS= NS-1 DO M=1,N-1 DO I=1,N-M H0= XA(I)-X HP= XA(I+M)-X W= C(I+1)-D(I) DEN= H0-HP IF(DEN.EQ.0.D0) PAUSE DEN= W/DEN D(I)= HP*DEN C(I)= H0*DEN ENDDO IF(2*NS.LT.N-M) THEN DY= C(NS+1) ELSE DY= D(NS) NS= NS-1 ENDIF Y= Y+DY ENDDO RETURN END * *-----RSPENCE----------------------------------------------------- * GIVES THE REAL PART OF LI_2(Z) IN DOUBLE PRECISION. ACCURACY * IS ABOUT 8 DIGITS * REAL*8 FUNCTION RSPENCE(XR,XI) IMPLICIT REAL*8(A,B,D-H,O-Z) IMPLICIT COMPLEX*16 (C) * COMMON/PARAM/PI,PIS,DELTA * DATA B0/1.D0/,B1/-0.25D0/,B2/0.277777777777778D-1/, # B4/-0.277777777777778D-3/,B6/0.472411186696901D-5/ * C0= CMPLX(0.D0) PIS6= PIS/6.D0 CX= CMPLX(XR,XI) COMX= 1.D0-CX IF (XR.LT.0.D0) THEN CY= 1.D0-CX SIGN1= -1.D0 CLNX= LOG(CX) CLNOMX= LOG(COMX) CADD1= PIS6-CLNX*CLNOMX ELSE CY= CX SIGN1= 1.D0 CADD1= C0 ENDIF COMY= 1.D0-CY YM2= CONJG(CY)*CY YM= SQRT(YM2) IF (YM.GT.1.D0) THEN CZ= CONJG(CY)/YM2 SIGN2= -1.D0 COY= -CY CLNOY= LOG(COY) CADD2= -PIS6-0.5D0*CLNOY*CLNOY ELSE CZ= CY SIGN2= 1.D0 CADD2= C0 ENDIF COMZ= 1.D0-CZ ZR= REAL(CZ) IF (ZR.GT.0.5D0) THEN CT= 1.D0-CZ COMT= 1.D0-CT SIGN3= -1.D0 CLNZ= LOG(CZ) CLNOMZ= LOG(COMZ) CADD3= PIS6-CLNZ*CLNOMZ ELSE CT= CZ COMT= 1.D0-CT SIGN3= 1.D0 CADD3= C0 ENDIF CPAR= LOG(COMT) CPAR2= CPAR*CPAR * CRES= -CPAR*(B0-CPAR*(B1-CPAR*(B2+CPAR2*(B4+B6*CPAR2)))) * CLI2= SIGN1*(SIGN2*(SIGN3*CRES+CADD3)+CADD2)+CADD1 * RSPENCE= REAL(CLI2) * RETURN END * * * *-----ELECTROWEAK LIBRARY QFORMFS * *------------------------------------------------------------------------ *-----QFORMFS-(WITHOUT 4-POINT FUNCTIONS)-------------------------------- *------------------------------------------------------------------------ * *-------------------------BFF-------------------------------------* * COMPUTES THE ONE-LOOP TWO-POINT FORM FACTORS /I*PI^2 * * B0,B1,B21 * * !!! B22 MUST BE COMPUTED SEPARATELY !!! * * INPUT PARAMETERS ARE P^2,M1^2,M2^2 * *-----------------------------------------------------------------* SUBROUTINE BFF(P2,RM12,RM22,B0,B1,B21) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA DIMENSION CMP2(2) DIMENSION CLNMP2(2) DIMENSION B0(2),B1(2),B21(2) DIMENSION GFPR(3),GFMR(3),GFPI(3),GFMI(3) CMP2(1)= -P2 CMP2(2)= -QEPS * M1^2= M2^2 ---> 0 CALL ROOTS(P2,RM12,RM22,RPR,RPI,RMR,RMI,OMRPR,OMRMR) * CALL CQLNX(CMP2,CLNMP2) N= 3 CALL CG(N,RPR,RPI,OMRPR,GFPR,GFPI) CALL CG(N,RMR,RMI,OMRMR,GFMR,GFMI) * AUXDEL= QDELTA-CLNMP2(1) B0(1)= AUXDEL-GFPR(1)-GFMR(1) B0(2)= -CLNMP2(2)-GFPI(1)-GFMI(1) B1(1)= -AUXDEL/2.Q0+GFPR(2)+GFMR(2) B1(2)= +CLNMP2(2)/2.Q0+GFPI(2)+GFMI(2) B21(1)= AUXDEL/3.Q0-GFPR(3)-GFMR(3) B21(2)= -CLNMP2(2)/3.Q0-GFPI(3)-GFMI(3) * RETURN END * *-------------------------RBFF------------------------------------* * COMPUTES THE REAL PART OF THE ONE-LOOP TWO-POINT * * FORM FACTORS /I*PI^2 * * B0,B1,B21 * * !!! B22 MUST BE COMPUTED SEPARATELY !!! * * INPUT PARAMETERS ARE P^2,M1^2,M2^2 * *-----------------------------------------------------------------* SUBROUTINE RBFF(P2,RM12,RM22,RB0,RB1,RB21) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA DIMENSION CMP2(2) DIMENSION CLNMP2(2) DIMENSION GFPR(3),GFMR(3) CMP2(1)= -P2 CMP2(2)= -QEPS * M1^2= M2^2 ---> 0 CALL ROOTS(P2,RM12,RM22,RPR,RPI,RMR,RMI,OMRPR,OMRMR) CALL CQLNX(CMP2,CLNMP2) N= 3 * CALL RCG(N,RPR,RPI,OMRPR,GFPR) CALL RCG(N,RMR,RMI,OMRMR,GFMR) * AUXDEL= QDELTA-CLNMP2(1) RB0= AUXDEL-GFPR(1)-GFMR(1) RB1= -AUXDEL/2.Q0+GFPR(2)+GFMR(2) RB21= AUXDEL/3.Q0-GFPR(3)-GFMR(3) * RETURN END * *-----------------BPFF--------------------------------------* * COMPUTES MINUS THE DERIVATIVE OF ONE-LOOP TWO-POINT * * FORM FACTORS /I*PI^2 * * B0P,B1P,B21P * * INPUT PARAMETERS ARE P^2,M1^2,M2^2 * *-----------------------------------------------------------* SUBROUTINE BPFF(P2,RM12,RM22,B0P,B1P,B21P) IMPLICIT REAL*16(A-H,O-Z) DIMENSION B0P(2),B1P(2),B21P(2) DIMENSION GFPR(4),GFPI(4),GFMR(4),GFMI(4) IF (RM12.EQ.RM22.AND.RM12.LT.1.Q-15) THEN B0P(1)= 1.Q0/P2 B0P(2)= 0.Q0 B1P(1)= -1.Q0/2.Q0/P2 B1P(2)= 0.Q0 B21P(1)= 1.Q0/3.Q0/P2 B21P(2)= 0.Q0 RETURN ELSE CALL ROOTS(P2,RM12,RM22,RPR,RPI,RMR,RMI,OMRPR,OMRMR) DRR= RPR-RMR DRI= RPI-RMI DRM2= DRR*DRR+DRI*DRI N= 4 CALL CG(N,RPR,RPI,OMRPR,GFPR,GFPI) CALL CG(N,RMR,RMI,OMRMR,GFMR,GFMI) EX0R= 2.Q0*(GFMR(2)-GFPR(2))-GFMR(1)+GFPR(1) EX0I= 2.Q0*(GFMI(2)-GFPI(2))-GFMI(1)+GFPI(1) EX1R= 2.Q0*(GFMR(2)-GFPR(2))-3.Q0*(GFMR(3)-GFPR(3)) EX1I= 2.Q0*(GFMI(2)-GFPI(2))-3.Q0*(GFMI(3)-GFPI(3)) EX2R= 4.Q0*(GFMR(4)-GFPR(4))-3.Q0*(GFMR(3)-GFPR(3)) EX2I= 4.Q0*(GFMI(4)-GFPI(4))-3.Q0*(GFMI(3)-GFPI(3)) B0P(1)= 1.Q0/DRM2/P2*(DRR*EX0R+DRI*EX0I) B0P(2)= 1.Q0/DRM2/P2*(DRR*EX0I-DRI*EX0R) B1P(1)= 1.Q0/DRM2/P2*(DRR*EX1R+DRI*EX1I) B1P(2)= 1.Q0/DRM2/P2*(DRR*EX1I-DRI*EX1R) B21P(1)= 1.Q0/DRM2/P2*(DRR*EX2R+DRI*EX2I) B21P(2)= 1.Q0/DRM2/P2*(DRR*EX2I-DRI*EX2R) RETURN ENDIF END *-----------------RBPFF-------------------------------------* * COMPUTES MINUS THE DERIVATIVE OF THE REAL PART OF * * ONE-LOOP TWO-POINT FORM FACTORS /I*PI^2 * * B0P,B1P,B21P * * INPUT PARAMETERS ARE P^2,M1^2,M2^2 * *-----------------------------------------------------------* SUBROUTINE RBPFF(P2,RM12,RM22,RB0P,RB1P,RB21P) IMPLICIT REAL*16(A-H,O-Z) DIMENSION GFPR(4),GFPI(4),GFMR(4),GFMI(4) IF (RM12.EQ.RM22.AND.RM12.LT.1.Q-15) THEN RB0P= 1.Q0/P2 RB1P= -1.Q0/2.Q0/P2 RB21P= 1.Q0/3.Q0/P2 RETURN ELSE CALL ROOTS(P2,RM12,RM22,RPR,RPI,RMR,RMI,OMRPR,OMRMR) DRR= RPR-RMR DRI= RPI-RMI DRM2= DRR*DRR+DRI*DRI N= 4 CALL CG(N,RPR,RPI,OMRPR,GFPR,GFPI) CALL CG(N,RMR,RMI,OMRMR,GFMR,GFMI) EX0R= 2.Q0*(GFMR(2)-GFPR(2))-GFMR(1)+GFPR(1) EX0I= 2.Q0*(GFMI(2)-GFPI(2))-GFMI(1)+GFPI(1) EX1R= 2.Q0*(GFMR(2)-GFPR(2))-3.Q0*(GFMR(3)-GFPR(3)) EX1I= 2.Q0*(GFMI(2)-GFPI(2))-3.Q0*(GFMI(3)-GFPI(3)) EX2R= 4.Q0*(GFMR(4)-GFPR(4))-3.Q0*(GFMR(3)-GFPR(3)) EX2I= 4.Q0*(GFMI(4)-GFPI(4))-3.Q0*(GFMI(3)-GFPI(3)) RB0P= 1.Q0/DRM2/P2*(DRR*EX0R+DRI*EX0I) RB1P= 1.Q0/DRM2/P2*(DRR*EX1R+DRI*EX1I) RB21P= 1.Q0/DRM2/P2*(DRR*EX2R+DRI*EX2I) RETURN ENDIF END *------------------------CFF-----------------------------------* * COMPUTES THE ONE-LOOP THREE-POINT FORM FACTORS/I*PI^2 * * C0,C11,C12,C21,C22,C23,C24 * * ACCORDING TO THE CONVENTION * * * * (Q^2+M1^2)((Q+P1)^2+M2^2)((Q+P1+P2)^2+M3^2) * * * * INPUT PARAMETERS ARE P1^2,P2^2, S=P5^2=(P1+P2)^2 * * M1^2,M2^2,M3^2 * * * * JFLAG = 1,2 SELECTS ONE OF THE ROOTS FOR ALPHA * * ! NOW ALPHA CAN BE COMPLEX * *--------------------------------------------------------------* SUBROUTINE CFF(JFLAG,P12,P22,S,RM12,RM22,RM32,C0,C1,C2) IMPLICIT REAL*16(A-H,O-Z) DIMENSION XMI(2,2),R1(2,2),TMP1(2,2),R2_13(2,2),TMP2_13(2,2) DIMENSION R2_32(2,2),TMP2_32(2,2),TMP0(2),TMP24(2) DIMENSION B0_12(2),B1_12(2),B21_12(2) DIMENSION B0_13(2),B1_13(2),B21_13(2) DIMENSION B0_23(2),B1_23(2),B21_23(2) DIMENSION C0(2),C1(2,2),C2(2,4) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA P1DP2= 0.5Q0*(S-P12-P22) DEN= -S*S/4.Q0*(1.Q0+((P12-P22)*(P12-P22)/S- # 2.Q0*(P12+P22))/S) DM12= RM12-RM22 DM13= RM12-RM32 DM23= RM22-RM32 RA= -P22 RB= -P12 RC= -2.Q0*P1DP2 RD= DM23+P22 RE= DM12+P12+2.Q0*P1DP2 RAA= -P22+DM23 RBB= P12+DM12 RCC= -2.Q0*P1DP2-P22+DM23 RDD= -P12+DM12 DISCP= +2.Q0*SQRT(ABS(DEN)) DISCM= -DISCP IF(DEN.LT.0.Q0) THEN CALL ALPHA(P12,P22,S,AP,AM,OMAP,OMAM) IF(JFLAG.EQ.1) THEN ALP= AP OMALP= OMAP DISC= DISCP ELSE IF(JFLAG.EQ.2) THEN ALP= AM OMALP= OMAM DISC= DISCM ELSE PRINT*,'JFLAG > 2' STOP ENDIF ALPI= 0.Q0 Y0R= -(RD+RE*ALP)/DISC Y0I= 0.Q0 Y1R= -(RAA+RBB*ALP)/DISC Y1I= Y0I OMY1R= (RCC+RDD*ALP)/DISC Y2R= Y0R/OMALP Y2I= Y0I/OMALP OMY2R= OMY1R/OMALP Y3R= -Y0R/ALP Y3I= -Y0I/ALP OMY3R= Y1R/ALP ELSE IF(DEN.GT.0) THEN ALP= (P12+P22-S)/2.Q0/P12 OMALP= (P12-P22+S)/2.Q0/P12 CALPI= SQRT(DEN)/P12 IF(JFLAG.EQ.1) THEN ALPI= -CALPI DISC= DISCP ELSE IF(JFLAG.EQ.2) THEN ALPI= CALPI DISC= DISCM ELSE PRINT*,'JFLAG > 2' STOP ENDIF Y0R= -RE*ALPI/DISC Y0I= (RD+RE*ALP)/DISC Y1R= -RBB*ALPI/DISC Y1I= (RAA+RBB*ALP)/DISC OMY1R= RDD*ALPI/DISC OMALPM2= OMALP*OMALP+ALPI*ALPI Y2R= (Y0R*OMALP-Y0I*ALPI)/OMALPM2 Y2I= (Y0R*ALPI+Y0I*OMALP)/OMALPM2 OMY2R= (OMY1R*OMALP+Y1I*ALPI)/OMALPM2 ALPM2= ALP*ALP+ALPI*ALPI Y3R= -(Y0R*ALP+Y0I*ALPI)/ALPM2 Y3I= -(-Y0R*ALPI+Y0I*ALP)/ALPM2 OMY3R= (Y1R*ALP+Y1I*ALPI)/ALPM2 ENDIF CALL ROOTS(P12,RM22,RM12,RP1R,RP1I,RM1R,RM1I,OMRP1R,OMRM1R) CALL ROOTS(S,RM32,RM12,RP2R,RP2I,RM2R,RM2I,OMRP2R,OMRM2R) CALL ROOTS(P22,RM32,RM22,RP3R,RP3I,RM3R,RM3I,OMRP3R,OMRM3R) A1= -P12 A2= -S A3= -P22 EB1R= -1.Q0-DM12/P12 EB1I= 0.Q0 EC1I= QEPS/P12 EB2R= -1.Q0-DM13/S EB2I= 0.Q0 EC2I= QEPS/S EB3R= -1.Q0-DM23/P22 EB3I= 0.Q0 EC3I= QEPS/P22 CALL S2(Y1R,Y1I,OMY1R,A1,EB1R,EB1I,EC1I,RP1R,RP1I,RM1R,RM1I, # OMRP1R,OMRM1R,S21R,S21I) CALL S2(Y2R,Y2I,OMY2R,A2,EB2R,EB2I,EC2I,RP2R,RP2I,RM2R,RM2I, # OMRP2R,OMRM2R,S22R,S22I) CALL S2(Y3R,Y3I,OMY3R,A3,EB3R,EB3I,EC3I,RP3R,RP3I,RM3R,RM3I, # OMRP3R,OMRM3R,S23R,S23I) * SCALAR 3-POINT FUNCTION C0 IF(DEN.LT.0.Q0) THEN TMP0(1)= (S21R-S22R+S23R)/DISC TMP0(2)= (S21I-S22I+S23I)/DISC ELSE IF(DEN.GT.0.Q0) THEN TMP0(1)= (S21I-S22I+S23I)/DISC TMP0(2)= -(S21R-S22R+S23R)/DISC ENDIF F1= DM12-P12 F2= DM23-S+P12 CALL BFF(P12,RM12,RM22,B0_12,B1_12,B21_12) CALL BFF(S,RM12,RM32,B0_13,B1_13,B21_13) CALL BFF(P22,RM22,RM32,B0_23,B1_23,B21_23) XMI(1,1)= P22/DEN XMI(1,2)= -P1DP2/DEN XMI(2,1)= XMI(1,2) XMI(2,2)= P12/DEN * C1 FORM FACTORS DO I= 1,2 R1(I,1)= 0.5Q0*(F1*TMP0(I)+B0_13(I)-B0_23(I)) R1(I,2)= 0.5Q0*(F2*TMP0(I)+B0_12(I)-B0_13(I)) ENDDO CALL MULTI2(TMP1,XMI,R1) * C2 FORM FACTORS TMP24(1)= 0.25Q0-0.5Q0*RM12*TMP0(1)+ # 0.25Q0*(B0_23(1)-F1*TMP1(1,1)-F2*TMP1(1,2)) TMP24(2)= -0.5Q0*RM12*TMP0(2)+ # 0.25Q0*(B0_23(2)-F1*TMP1(2,1)-F2*TMP1(2,2)) DO I= 1,2 R2_13(I,1)= 0.5Q0*(F1*TMP1(I,1)+B1_13(I)+B0_23(I))-TMP24(I) R2_13(I,2)= 0.5Q0*(F2*TMP1(I,1)+B1_12(I)-B1_13(I)) ENDDO CALL MULTI2(TMP2_13,XMI,R2_13) DO I= 1,2 R2_32(I,1)= 0.5Q0*(F1*TMP1(I,2)+B1_13(I)-B1_23(I)) R2_32(I,2)= 0.5Q0*(F2*TMP1(I,2)-B1_13(I))-TMP24(I) ENDDO CALL MULTI2(TMP2_32,XMI,R2_32) DO I= 1,2 C0(I)= TMP0(I) C1(I,1)= TMP1(I,1) C1(I,2)= TMP1(I,2) C2(I,1)= TMP2_13(I,1) C2(I,2)= TMP2_32(I,2) C2(I,3)= TMP2_13(I,2) C2(I,4)= TMP24(I) ENDDO RETURN END *------------------------------------------------------------* * COMPUTES THE FUNCTION S2 * *------------------------------------------------------------* SUBROUTINE S2(Y0R,Y0I,OMY0R,A,EBR,EBI,ECI,RPR,RPI,RMR,RMI, # OMRPR,OMRMR,S2R,S2I) IMPLICIT REAL*16(A-H,O-Z) DIMENSION CARG(2),COMARG(2),CLN(2) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA * RP AND RM ARE THE ROOTS OF A*X^2+B*X+C = 0 * EB= B/A AND EC= C/A ARE KEPT TO COMPUTE IM(RP*RM)= IM(EC) * AND IM{(Y0-RP)(Y0-RM)}= IM{Y0^2+EB*Y0+EC} A1= A*ECI A2= A*(2.Q0*Y0R*Y0I+EBR*Y0I+EBI*Y0R+ECI) IF (A1.LT.0.Q0) THEN DEL= 1.Q-90 ELSE DEL= -1.Q-90 ENDIF IF (A2.LT.0.Q0) THEN DELP= 1.Q-90 ELSE DELP= -1.Q-90 ENDIF * COMPUTES THE LOG WHICH OCCURS IN ASSOCIATION WITH ETA-FUNCTIONS Y0M2= Y0R*Y0R+Y0I*Y0I CARG(1)= Y0R/Y0M2 CARG(2)= -Y0I/Y0M2 COMARG(1)= -(OMY0R*Y0R-Y0I*Y0I)/Y0M2 COMARG(2)= Y0I/Y0M2 CALL CQLNOMX(CARG,COMARG,CLN) * IN CALLING ETA ONLY THE SIGN OF THE ARGUMENTS IS RELEVANT A11= -RPI A12= -RMI A1P= ECI A21= Y0I-RPI A22= Y0I-RMI A2P= 2.Q0*Y0R*Y0I+EBR*Y0I+EBI*Y0R+ECI A31= -DEL A32= DELP A3P= A*(DELP-DEL) CALL RLOG(Y0R,Y0I,OMY0R,RPR,RPI,OMRPR,RFPR,RFPI) CALL RLOG(Y0R,Y0I,OMY0R,RMR,RMI,OMRMR,RFMR,RFMI) ETA1= ETA(A11,A12,A1P) ETA2= ETA(A21,A22,A2P) ETA3= ETA(A31,A32,A3P) S2R= RFPR+RFMR-CLN(2)*(ETA1-ETA2-ETA3) S2I= RFPI+RFMI+CLN(1)*(ETA1-ETA2-ETA3) RETURN END *------------------------------------------------------------* * COMPUTES THE FUNCTION S1 * *------------------------------------------------------------* SUBROUTINE S1(Y0R,Y0I,OMY0R,B,RR,RI,S1R,S1I) IMPLICIT REAL*16(A-H,O-Z) DIMENSION CARG(2),COMARG(2),CLN(2) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA B1= -B*RI B2= B*(Y0I-RI) IF (B1.LT.0.Q0) THEN DEL= 1.Q-90 ELSE DEL= -1.Q-90 ENDIF IF (B2.LT.0.Q0) THEN DELP= 1.Q-90 ELSE DELP= -1.Q-90 ENDIF * COMPUTES THE LOG WHICH OCCURS IN ASSOCIATION WITH ETA-FUNCTIONS Y0I2= Y0I*Y0I Y0M2= Y0R*Y0R+Y0I2 SY0= Y0I/Y0M2 CARG(1)= Y0R/Y0M2 CARG(2)= -SY0 COMARG(1)= -(OMY0R*Y0R-Y0I2)/Y0M2 COMARG(2)= SY0 CALL CQLNOMX(CARG,COMARG,CLN) * IN CALLING ETA ONLY THE SIGN OF THE ARGUMENTS IS RELEVANT A1= -DEL A2= DELP AP= B*(DELP-DEL) OMRR= 1.Q0-RR CALL RLOG(Y0R,Y0I,OMY0R,RR,RI,OMRR,RFR,RFI) ETAS= ETA(A1,A2,AP) S1R= RFR+CLN(2)*ETAS S1I= RFI-CLN(1)*ETAS RETURN END *--------------------------------------------------* * COMPUTES THE FUNCTION R * *--------------------------------------------------* SUBROUTINE RLOG(Z0R,Z0I,OMZ0R,Z1R,Z1I,OMZ1R,RFUNR,RFUNI) IMPLICIT REAL*16(A-H,O-Z) DIMENSION CA1(2),CA2(2),CA3(2),ADD(2) DIMENSION CLN1(2),CLN2(2),CLN3(2),CT(10),SN(10) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA Z1M2= Z1R*Z1R+Z1I*Z1I Z1M= SQRT(Z1M2) Z0M2= Z0R*Z0R+Z0I*Z0I Z0M= SQRT(Z0M2) ZRT= Z1M/Z0M * |Z1| << |Z0| IF (ZRT.LT.1.Q-10) THEN CT(1)= (Z1R*Z0R+Z1I*Z0I)/Z1M/Z0M SN(1)= (Z0R*Z1I-Z1R*Z0I)/Z1M/Z0M DO K=2,10 CT(K)= CT(1)*CT(K-1)-SN(1)*SN(K-1) SN(K)= SN(1)*CT(K-1)+CT(1)*SN(K-1) ENDDO SUMR= CT(10) SUMI= SN(10) DO J=9,1,-1 SUMR= SUMR*ZRT+CT(J) SUMI= SUMI*ZRT+SN(J) ENDDO SUMR= SUMR*ZRT SUMI= SUMI*ZRT A1R= 1.Q0+SUMR A1I= SUMI A2R= (Z0I*Z0I-OMZ0R*Z0R)/Z0M2*A1R-Z0I/Z0M2*A1I A2I= (Z0I*Z0I-OMZ0R*Z0R)/Z0M2*A1I+Z0I/Z0M2*A1R OMA1R= -SUMR OMA2R= 1.Q0-A2R Z01R= Z0R-Z1R Z01I= Z0I-Z1I DEN= Z01R*Z01R+Z01I*Z01I ADD(1)= 0.Q0 ADD(2)= 0.Q0 SIGN= +1.Q0 ELSE * IF Z0R AND/OR Z1R ARE ROOTS VERY NEAR TO 1 * Z0R-Z1R = (1-Z1R)-(1-Z0R) AOMZ0R= ABS(OMZ0R) AOMZ1R= ABS(OMZ1R) AZ1I= ABS(Z1I) IF (AOMZ0R.LT.1.Q-10.OR.AOMZ1R.LT.1.Q-10) THEN Z01R= OMZ1R-OMZ0R ELSE Z01R= Z0R-Z1R ENDIF Z01I= Z0I-Z1I * IF Z0 AND Z1R ARE ROOTS VERY NEAR TO 1 * AND |Z1I| << 1 , Z0I = 0 IF (AOMZ0R.LT.1.Q-10.AND.AOMZ1R.LT.1.Q-10. # AND.AZ1I.LT.1.Q-10.AND.Z0I.EQ.0.Q0) THEN DEN= Z01R*Z01R+Z1I*Z1I A1R= Z01R/Z0R A1I= -Z1I/Z0R OMA1R= Z1R/Z0R CA3(1)= -A1R CA3(2)= -A1I CALL CQLNX(CA3,CLN3) ADD(1)= -QPIS/6.Q0-0.5Q0*(CLN3(1)*CLN3(1)-CLN3(2)*CLN3(2)) ADD(2)= -CLN3(1)*CLN3(2) SIGN= -1.Q0 ELSE DEN= Z01R*Z01R+Z01I*Z01I A1R= (Z0R*Z01R+Z0I*Z01I)/DEN A1I= (-Z0R*Z01I+Z0I*Z01R)/DEN OMA1R= -(Z01I*Z1I+Z1R*Z01R)/DEN ADD(1)= 0.Q0 ADD(2)= 0.Q0 SIGN= +1.Q0 ENDIF A2R= -(OMZ0R*Z01R-Z0I*Z01I)/DEN A2I= (OMZ0R*Z01I+Z0I*Z01R)/DEN OMA2R= (OMZ1R*Z01R-Z01I*Z1I)/DEN ENDIF * IN CALLING ETA ONLY THE SIGN OF THE ARGUMENTS IS RELEVANT A1= -Z1I A2= -Z01I API= Z1R*Z0I-Z0R*Z1I APII= -Z0I*OMZ1R+Z1I*OMZ0R CALL SPENCE(A1R,A1I,OMA1R,SP1R,SP1I) CALL SPENCE(A2R,A2I,OMA2R,SP2R,SP2I) CA1(1)= A1R CA1(2)= A1I CA2(1)= A2R CA2(2)= A2I CALL CQLNX(CA1,CLN1) CALL CQLNX(CA2,CLN2) ETAI= ETA(A1,A2,API) ETAII= ETA(A1,A2,APII) RFUNR= SIGN*SP1R-SP2R-ETAI*CLN1(2)+ETAII*CLN2(2)+ADD(1) RFUNI= SIGN*SP1I-SP2I+ETAI*CLN1(1)-ETAII*CLN2(1)+ADD(2) RETURN END *----------------------------------------------------------* * COMPUTES THE FUNTION ETA(A,B) * * LN(AB)= LN(A)+LN(B)+ETA(A,B) * *----------------------------------------------------------* REAL*16 FUNCTION ETA(Z1I,Z2I,ZPI) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA * ONLY THE SIGN OF THE ARGUMENTS IS RELEVANT IF (Z1I.LT.0.Q0.AND.Z2I.LT.0.Q0.AND.ZPI.GT.0.Q0) THEN ETA= 2.Q0*QPI RETURN ELSE IF (Z1I.GT.0.Q0.AND.Z2I.GT.0.Q0.AND.ZPI.LT.0.Q0) THEN ETA= -2.Q0*QPI RETURN ELSE ETA= 0.Q0 RETURN ENDIF END *----------------------------------------------------------* * COMPUTES THE FUNTION ETA(A,1/B) * * LN(A/B)= LN(A)-LN(B)+ETA(A,1/B) * *----------------------------------------------------------* REAL*16 FUNCTION ETAQ(Z1I,Z2I,ZQI) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA IF (Z1I.LT.0.Q0.AND.Z2I.GT.0.Q0.AND.ZQI.GT.0.Q0) THEN ETAQ= 2.Q0*QPI RETURN ELSE IF (Z1I.GT.0.Q0.AND.Z2I.LT.0.Q0.AND.ZQI.LT.0.Q0) THEN ETAQ= -2.Q0*QPI RETURN ELSE ETAQ= 0.Q0 RETURN ENDIF END * *-----SPENCE-------------------------------------------------------- * COMPUTES LI_2(X). ACCURACY IS ABOUT 16 DIGITS * SUBROUTINE SPENCE(XR,XI,OMXR,CLI2R,CLI2I) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA DIMENSION B(0:14),BF(0:14) DIMENSION X(2),OMX(2),Y(2),OY(2),OMY(2),Z(2),OMZ(2),T(2),OMT(2) DIMENSION CLNX(2),CLNOMX(2),CLNOY(2),CLNZ(2),CLNOMZ(2) DIMENSION ADD1(2),ADD2(2),ADD3(2),PAR(2),RES(2),CT(15),SN(15) X(1)= XR X(2)= XI OMX(1)= OMXR OMX(2)= -XI IF (XR.LT.0.Q0) THEN Y(1)= OMXR Y(2)= -XI SIGN1= -1.Q0 CALL CQLNX(X,CLNX) CALL CQLNOMX(X,OMX,CLNOMX) ADD1(1)= QPIS/6.Q0-CLNX(1)*CLNOMX(1)+CLNX(2)*CLNOMX(2) ADD1(2)= -CLNX(1)*CLNOMX(2)-CLNX(2)*CLNOMX(1) ELSE Y(1)= X(1) Y(2)= X(2) SIGN1= 1.Q0 ADD1(1)= 0.Q0 ADD1(2)= 0.Q0 ENDIF OMY(1)= 1.Q0-Y(1) OMY(2)= -Y(2) YM2= Y(1)*Y(1)+Y(2)*Y(2) YM= SQRT(YM2) IF (YM.GT.1.Q0) THEN Z(1)= Y(1)/YM2 Z(2)= -Y(2)/YM2 SIGN2= -1.Q0 OY(1)= -Y(1) OY(2)= -Y(2) CALL CQLNX(OY,CLNOY) ADD2(1)= -QPIS/6.Q0-0.5Q0*((CLNOY(1))**2-(CLNOY(2))**2) ADD2(2)= -CLNOY(1)*CLNOY(2) ELSE Z(1)= Y(1) Z(2)= Y(2) SIGN2= 1.Q0 ADD2(1)= 0.Q0 ADD2(2)= 0.Q0 ENDIF OMZ(1)= 1.Q0-Z(1) OMZ(2)= -Z(2) ZR= Z(1) IF (ZR.GT.0.5Q0) THEN T(1)= 1.Q0-Z(1) T(2)= -Z(2) OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) SIGN3= -1.Q0 CALL CQLNX(Z,CLNZ) CALL CQLNOMX(Z,OMZ,CLNOMZ) ADD3(1)= QPIS/6.Q0-CLNZ(1)*CLNOMZ(1)+CLNZ(2)*CLNOMZ(2) ADD3(2)= -CLNZ(1)*CLNOMZ(2)-CLNZ(2)*CLNOMZ(1) ELSE T(1)= Z(1) T(2)= Z(2) OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) SIGN3= 1.Q0 ADD3(1)= 0.Q0 ADD3(2)= 0.Q0 ENDIF CALL CQLNOMX(T,OMT,PAR) B(0)= 1.Q0 B(1)= -1.Q0/2.Q0 B(2)= 1.Q0/6.Q0 B(4)= -1.Q0/30.Q0 B(6)= 1.Q0/42.Q0 B(8)= -1.Q0/30.Q0 B(10)= 5.Q0/66.Q0 B(12)= -691.Q0/2730.Q0 B(14)= 7.Q0/6.Q0 FACT= 1.Q0 DO N=0,14 BF(N)= B(N)/FACT FACT= FACT*(N+2.Q0) ENDDO PARR= PAR(1) PARI= PAR(2) PARM2= PARR*PARR+PARI*PARI PARM= SQRT(PARM2) CT(1)= PARR/PARM SN(1)= PARI/PARM DO N=2,15 CT(N)= CT(1)*CT(N-1)-SN(1)*SN(N-1) SN(N)= SN(1)*CT(N-1)+CT(1)*SN(N-1) ENDDO RES(1)= -((((((((BF(14)*CT(15)*PARM2+BF(12)*CT(13))*PARM2+ # BF(10)*CT(11))*PARM2+BF(8)*CT(9))*PARM2+ # BF(6)*CT(7))*PARM2+BF(4)*CT(5))*PARM2+ # BF(2)*CT(3))*(-PARM)+BF(1)*CT(2))*(-PARM)+ # BF(0)*CT(1))*PARM RES(2)= -((((((((BF(14)*SN(15)*PARM2+BF(12)*SN(13))*PARM2+ # BF(10)*SN(11))*PARM2+BF(8)*SN(9))*PARM2+ # BF(6)*SN(7))*PARM2+BF(4)*SN(5))*PARM2+ # BF(2)*SN(3))*(-PARM)+BF(1)*SN(2))*(-PARM)+ # BF(0)*SN(1))*PARM CLI2R= SIGN1*(SIGN2*(SIGN3*RES(1)+ADD3(1))+ADD2(1))+ADD1(1) CLI2I= SIGN1*(SIGN2*(SIGN3*RES(2)+ADD3(2))+ADD2(2))+ADD1(2) RETURN END *--------------------------------------------------* * COMPUTES LN(Z) * *--------------------------------------------------* SUBROUTINE CQLNX(ARG,RES) IMPLICIT REAL*16(A-H,O-Z) DIMENSION ARG(2),AARG(2),RES(2) QPI= 3.141592653589793238462643Q0 DO I= 1,2 AARG(I)= ABS(ARG(I)) ENDDO ZM2= (ARG(1))**2+(ARG(2))**2 ZM= SQRT(ZM2) RES(1)= LOG(ZM) IF (ARG(1).EQ.0.Q0) THEN IF (ARG(2).GT.0.Q0) THEN TETA= QPI/2.Q0 ELSE TETA= -QPI/2.Q0 ENDIF RES(2)= TETA RETURN ELSE IF (ARG(2).EQ.0.Q0) THEN IF (ARG(1).GT.0.Q0) THEN TETA= 0.Q0 ELSE TETA= QPI ENDIF RES(2)= TETA RETURN ELSE TNTETA= AARG(2)/AARG(1) TETA= ATAN(TNTETA) SR= ARG(1)/AARG(1) SI= ARG(2)/AARG(2) IF (SR.GT.0.Q0) THEN RES(2)= SI*TETA ELSE RES(2)= SI*(QPI-TETA) ENDIF RETURN ENDIF END *--------------------------------------* * COMPUTES LN(1-X) * * USUALLY |X| << 1 * *--------------------------------------* SUBROUTINE CQLNOMX(ARG,OMARG,RES) IMPLICIT REAL*16(A-H,O-Z) DIMENSION ARG(2),OMARG(2),RES(2),ARES(2),CT(10),SN(10) ZR= ARG(1) ZI= ARG(2) ZM2= ZR*ZR+ZI*ZI ZM= SQRT(ZM2) IF (ZM.LT.1.Q-7) THEN CT(1)= ZR/ZM SN(1)= ZI/ZM DO N=2,10 CT(N)= CT(1)*CT(N-1)-SN(1)*SN(N-1) SN(N)= SN(1)*CT(N-1)+CT(1)*SN(N-1) ENDDO ARES(1)= CT(10)/10.Q0 ARES(2)= SN(10)/10.Q0 DO K=9,1,-1 ARES(1)= ARES(1)*ZM+CT(K)/K ARES(2)= ARES(2)*ZM+SN(K)/K ENDDO ARES(1)= -ARES(1)*ZM ARES(2)= -ARES(2)*ZM ELSE CALL CQLNX(OMARG,ARES) ENDIF DO I= 1,2 RES(I)= ARES(I) ENDDO RETURN END *-----------------------------------------------------* * COMPUTES THE FUNCTION G_N(X) * *-----------------------------------------------------* SUBROUTINE CG(N,ZR,ZI,OMZR,GFR,GFI) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA DIMENSION GFR(N),GFI(N) DIMENSION A(2,4),RA(16),Z(2),OMZ(2),OZ(2),CLNOMZ(2),CLNOZ(2) DIMENSION CA(2,10),AUX(2),ZP(2),CT(16),SN(16) DATA RA/-0.2Q+0,0.666666666666667Q-1, # -0.952380952380952Q-2,-0.396825396825397Q-3, # 0.317460317460317Q-3,-0.132275132275132Q-4, # -0.962000962000962Q-5,0.105218855218855Q-5, # 0.266488361726450Q-6,-0.488745528428000Q-7, # -0.675397500794000Q-8,0.190720263471000Q-8, # 0.153663007690000Q-9,-0.679697905790000Q-10, # -0.293683556000000Q-11,0.228836696000000Q-11/ Z(1)= ZR Z(2)= ZI OMZ(1)= OMZR OMZ(2)= -ZI IF (ZR.EQ.0.Q0.AND.ZI.EQ.0.Q0) THEN DO K=1,N GFR(K)= -1.Q0/K/K GFI(K)= 0.Q0 ENDDO RETURN ELSE IF (ZR.EQ.1.Q0.AND.ZI.EQ.0.Q0) THEN A(1,1)= -1.Q0 A(2,1)= QPI DO J=2,4 A(1,J)= ((J-1.Q0)*A(1,J-1)-1.Q0/J)/J A(2,J)= (J-1.Q0)/J*A(2,J-1) ENDDO DO K=1,N GFR(K)= A(1,K) GFI(K)= A(2,K) ENDDO RETURN ELSE ZMOD2= ZR*ZR+ZI*ZI ZMOD= SQRT(ZMOD2) CALL CQLNX(OMZ,CLNOMZ) * |Z| < 4 IF (ZMOD.LT.4.Q0) THEN OZ(1)= -Z(1) OZ(2)= -Z(2) CALL CQLNX(OZ,CLNOZ) CA(1,1)= OMZ(1)*CLNOMZ(1)-OMZ(2)*CLNOMZ(2)+Z(1)*CLNOZ(1)- # Z(2)*CLNOZ(2)-1.Q0 CA(2,1)= OMZ(1)*CLNOMZ(2)+OMZ(2)*CLNOMZ(1)+Z(1)*CLNOZ(2)+ # Z(2)*CLNOZ(1) IF (N.EQ.1) THEN GFR(1)= CA(1,1) GFI(1)= CA(2,1) RETURN ELSE DO J= 2,N JM1= J-1 CA(1,J)= ((J-1.Q0)*(Z(1)*CA(1,JM1)-Z(2)*CA(2,JM1))+ # OMZ(1)*CLNOMZ(1)-OMZ(2)*CLNOMZ(2)-1.Q0/J)/J CA(2,J)= ((J-1.Q0)*(Z(1)*CA(2,JM1)+Z(2)*CA(1,JM1))+ # OMZ(1)*CLNOMZ(2)+OMZ(2)*CLNOMZ(1))/J ENDDO DO K=1,N GFR(K)= CA(1,K) GFI(K)= CA(2,K) ENDDO RETURN ENDIF * |Z| > 4 ELSE AUX(1)= (-ZR*OMZR+ZI**2)/ZMOD2 AUX(2)= ZI/ZMOD2 CALL CQLNX(AUX,ZP) ZPM2= ZP(1)*ZP(1)+ZP(2)*ZP(2) ZPM= SQRT(ZPM2) CT(1)= ZP(1)/ZPM SN(1)= ZP(2)/ZPM DO K=2,16 CT(K)= CT(1)*CT(K-1)-SN(1)*SN(K-1) SN(K)= SN(1)*CT(K-1)+CT(1)*SN(K-1) ENDDO CA(1,4)= RA(16)*CT(16)*ZPM+RA(15)*CT(15) CA(2,4)= RA(16)*SN(16)*ZPM+RA(15)*SN(15) DO J=14,1,-1 CA(1,4)= CA(1,4)*ZPM+RA(J)*CT(J) CA(2,4)= CA(2,4)*ZPM+RA(J)*SN(J) ENDDO CA(1,4)= CA(1,4)*ZPM CA(2,4)= CA(2,4)*ZPM DO J= 3,1,-1 JP1= J+1 CA(1,J)= ((CA(1,JP1)+1.Q0/JP1)*Z(1)+CA(2,JP1)*Z(2))/ZMOD2 CA(2,J)= (CA(2,JP1)*Z(1)-(CA(1,JP1)+1.Q0/JP1)*Z(2))/ZMOD2 ENDDO DO K=1,N GFR(K)= (CA(1,K)+CLNOMZ(1))/K GFI(K)= (CA(2,K)+CLNOMZ(2))/K ENDDO RETURN ENDIF ENDIF END *-----------------------------------------------------* * COMPUTES THE FUNCTION RE G_N(X) * *-----------------------------------------------------* SUBROUTINE RCG(N,ZR,ZI,OMZR,GFR) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA DIMENSION GFR(N) DIMENSION A(4),RA(16),Z(2),OMZ(2),OZ(2),CLNOMZ(2),CLNOZ(2) DIMENSION CA(2,10),AUX(2),ZP(2),CT(16),SN(16) DATA RA/-0.2Q+0,0.666666666666667Q-1, # -0.952380952380952Q-2,-0.396825396825397Q-3, # 0.317460317460317Q-3,-0.132275132275132Q-4, # -0.962000962000962Q-5,0.105218855218855Q-5, # 0.266488361726450Q-6,-0.488745528428000Q-7, # -0.675397500794000Q-8,0.190720263471000Q-8, # 0.153663007690000Q-9,-0.679697905790000Q-10, # -0.293683556000000Q-11,0.228836696000000Q-11/ Z(1)= ZR Z(2)= ZI OMZ(1)= OMZR OMZ(2)= -ZI IF (ZR.EQ.0.Q0.AND.ZI.EQ.0.Q0) THEN DO K=1,N GFR(K)= -1.Q0/K/K ENDDO RETURN ELSE IF (ZR.EQ.1.Q0.AND.ZI.EQ.0.Q0) THEN A(1)= -1.Q0 DO J=2,4 A(J)= ((J-1.Q0)*A(J-1)-1.Q0/J)/J ENDDO DO K=1,N GFR(K)= A(K) ENDDO RETURN ELSE ZMOD2= ZR*ZR+ZI*ZI ZMOD= SQRT(ZMOD2) CALL CQLNX(OMZ,CLNOMZ) * |Z| < 4 IF (ZMOD.LT.4.Q0) THEN OZ(1)= -Z(1) OZ(2)= -Z(2) CALL CQLNX(OZ,CLNOZ) CA(1,1)= OMZ(1)*CLNOMZ(1)-OMZ(2)*CLNOMZ(2)+Z(1)*CLNOZ(1)- # Z(2)*CLNOZ(2)-1.Q0 CA(2,1)= OMZ(1)*CLNOMZ(2)+OMZ(2)*CLNOMZ(1)+Z(1)*CLNOZ(2)+ # Z(2)*CLNOZ(1) IF (N.EQ.1) THEN GFR(1)= CA(1,1) RETURN ELSE DO J= 2,N JM1= J-1 CA(1,J)= ((J-1.Q0)*(Z(1)*CA(1,JM1)-Z(2)*CA(2,JM1))+ # OMZ(1)*CLNOMZ(1)-OMZ(2)*CLNOMZ(2)-1.Q0/J)/J CA(2,J)= ((J-1.Q0)*(Z(1)*CA(2,JM1)+Z(2)*CA(1,JM1))+ # OMZ(1)*CLNOMZ(2)+OMZ(2)*CLNOMZ(1))/J ENDDO DO K=1,N GFR(K)= CA(1,K) ENDDO RETURN ENDIF * |Z| > 4 ELSE AUX(1)= (-ZR*OMZR+ZI**2)/ZMOD2 AUX(2)= ZI/ZMOD2 CALL CQLNX(AUX,ZP) ZPM2= ZP(1)*ZP(1)+ZP(2)*ZP(2) ZPM= SQRT(ZPM2) CT(1)= ZP(1)/ZPM SN(1)= ZP(2)/ZPM DO K=2,16 CT(K)= CT(1)*CT(K-1)-SN(1)*SN(K-1) SN(K)= SN(1)*CT(K-1)+CT(1)*SN(K-1) ENDDO CA(1,4)= RA(16)*CT(16)*ZPM+RA(15)*CT(15) CA(2,4)= RA(16)*SN(16)*ZPM+RA(15)*SN(15) DO J=14,1,-1 CA(1,4)= CA(1,4)*ZPM+RA(J)*CT(J) CA(2,4)= CA(2,4)*ZPM+RA(J)*SN(J) ENDDO CA(1,4)= CA(1,4)*ZPM CA(2,4)= CA(2,4)*ZPM DO J= 3,1,-1 JP1= J+1 CA(1,J)= ((CA(1,JP1)+1.Q0/JP1)*Z(1)+CA(2,JP1)*Z(2))/ZMOD2 CA(2,J)= (CA(2,JP1)*Z(1)-(CA(1,JP1)+1.Q0/JP1)*Z(2))/ZMOD2 ENDDO DO K=1,N GFR(K)= (CA(1,K)+CLNOMZ(1))/K ENDDO RETURN ENDIF ENDIF END *------------------------------------------------* * COMPUTES E(I)= XMI(I,J)*D(J) * * I,J=1,2 * *------------------------------------------------C SUBROUTINE MULTI2(EA,XMI,DA) IMPLICIT REAL*16(A-H,O-Z) DIMENSION EA(2,2),DA(2,2),XMI(2,2) DO J= 1,2 DO I= 1,2 EA(J,I)= XMI(I,1)*DA(J,1)+XMI(I,2)*DA(J,2) ENDDO ENDDO RETURN END *------------------------------------------------* * COMPUTES E(I)= XMI(I,J)*D(J) * * I,J=1,3 * *------------------------------------------------* SUBROUTINE MULTI3(EA,XMI,DA) IMPLICIT REAL*16(A-H,O-Z) DIMENSION EA(2,3),DA(2,3),XMI(3,3) DO J= 1,2 DO I= 1,3 EA(J,I)= XMI(I,1)*DA(J,1)+XMI(I,2)*DA(J,2)+XMI(I,3)*DA(J,3) ENDDO ENDDO RETURN END *-------------------------------------------------* * COMPUTES THE ROOTS OF THE QUADRATIC FORM * * -P2*X^2+(P2+M2^2-M1^2)*X+M1^2 = 0 * * WITH: * * REAL MASSES M^2=M^2-I*EPS * *-------------------------------------------------* SUBROUTINE ROOTS(P2,RM12,RM22,RPR,RPI,RMR,RMI,OMRPR,OMRMR) IMPLICIT REAL*16(A-H,O-Z) COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA IF (RM12.EQ.RM22) THEN DISC2= (P2+4.Q0*RM12)*P2 ELSE DISC2= (RM22+2.Q0*(P2-RM12))*RM22+(P2+RM12)*(P2+RM12) ENDIF IF (DISC2.GT.0.Q0) THEN DISC= SQRT(DISC2) RPI= +QEPS/DISC RMI= -QEPS/DISC R1= ABS(RM22/P2) R2= ABS(RM22/RM12) * ONE OF THE ROOTS VERY NEAR TO 1 IF (R1.LT.1.Q-10.OR.R2.LT.1.Q-10) THEN A= -P2 B= P2+RM12-RM22 C= RM22 IF (B.GT.0.Q0) THEN Q= -0.5Q0*(B+DISC) RMR= 1.Q0-C/Q RPR= 1.Q0-Q/A OMRPR= Q/A OMRMR= C/Q ELSE Q= -0.5Q0*(B-DISC) RMR= 1.Q0-Q/A RPR= 1.Q0-C/Q OMRPR= C/Q OMRMR= Q/A ENDIF ELSE A= -P2 B= P2+RM22-RM12 C= RM12 IF (B.GT.0.Q0) THEN Q= -0.5Q0*(B+DISC) RPR= C/Q RMR= Q/A OMRPR= 1.Q0-RPR OMRMR= 1.Q0-RMR ELSE Q= -0.5Q0*(B-DISC) RPR= Q/A RMR= C/Q OMRPR= 1.Q0-RPR OMRMR= 1.Q0-RMR ENDIF ENDIF ELSE DISC= SQRT(-DISC2) A= -P2 B= P2+RM22-RM12 RPR= -B/2.Q0/A RMR= -B/2.Q0/A OMRPR= 1.Q0-RPR OMRMR= 1.Q0-RMR RPI= +DISC/2.Q0/A RMI= -DISC/2.Q0/A ENDIF RETURN END *-------------------------------------------------* * COMPUTES THE ROOTS OF * * -P1^2*ALP^2+(P1^2+P2^2-S)*ALP-P2^2= 0 * *-------------------------------------------------* SUBROUTINE ALPHA(P12,P22,S,ALP,ALM,OMALP,OMALM) IMPLICIT REAL*16(A-H,O-Z) ARG= (P22-2.Q0*(P12+S))*P22+(S-P12)*(S-P12) RT= SQRT(ARG) S1= ABS(S/P12) S2= ABS(S/P22) * ONE OF THE ROOTS VERY NEAR TO 1 IF (S1.LT.1.Q-10.OR.S2.LT.1.Q-10) THEN B= P12-P22+S IF (B.GT.0.Q0) THEN Q= -0.5Q0*(B+RT) BEP= -S/Q BEM= -Q/P12 ELSE Q= -0.5Q0*(B-RT) BEP= -Q/P12 BEM= -S/Q ENDIF ALP= 1.Q0-BEM ALM= 1.Q0-BEP OMALP= BEM OMALM= BEP RETURN ELSE B= P12+P22-S IF (B.GT.0.Q0) THEN Q= -0.5Q0*(B+RT) ALP= -P22/Q ALM= -Q/P12 ELSE Q= -0.5Q0*(B-RT) ALP= -Q/P12 ALM= -P22/Q ENDIF OMALP= 1.Q0-ALP OMALM= 1.Q0-ALM RETURN ENDIF END * *-----ALALS * SUBROUTINE ALALS(P2X,V1,A1,F1,V1P,A1P,V1I) IMPLICIT REAL*16(A-H,O-Z) * COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA COMMON/QMASSES2/QEM2,QMM2,QTM2,QNM2,QUQM2,QDQM2,QCQM2,QSQM2, # QBQM2,QTQM2,QZM2,QWM2,QHM2 * DIMENSION R(2),X(2),XO(2),UMX(2),A(2),B(2),B2(2),AB(2), # RMR(2),RUMR(2),RP(2),RM(2),F(2),H(2),G(2),RM2(2), # RM4(2),RI(2),R2I(2),XI(2),X2I(2),RAR(2),RRAR(2), # XP(2),UMR(2),TRMR(2),TRUMR(2),RPP(2) * DATA Z3/1.20205690315959428540Q0/ * *-----THIS SUBROUTINE USES THE WRONG METRIC, THUS THE SIGN OF P^2 * MUST BE CHANGED * Z2= QPIS/6.Q0 R(1)= -0.25Q0*P2X/QTQM2 R(2)= 0.25Q0*QEPS/QTQM2 X(1)= -P2X/QTQM2 X(2)= QEPS/QTQM2 DO I=1,2 XO(I)= -X(I) ENDDO UMX(1)= 1.Q0-X(1) UMX(2)= -X(2) UMR(1)= 1.Q0-R(1) UMR(2)= -R(2) * CALL CQLNX(XO,A) CALL CQLNOMX(X,UMX,B) B2(1)= B(1)*B(1)-B(2)*B(2) B2(2)= 2.Q0*B(1)*B(2) AB(1)= A(1)*B(1)-A(2)*B(2) AB(2)= A(1)*B(2)+A(2)*B(1) * IF(R(1).GT.0.Q0) THEN RMR(1)= 0.5Q0*R(2)/SQRT(R(1)) RMR(2)= -SQRT(R(1)) ELSE IF(R(1).LT.0.Q0) THEN RMR(1)= SQRT(-R(1)) RMR(2)= -0.5Q0*R(2)/SQRT(-R(1)) ENDIF IF(UMR(1).GT.0.Q0) THEN RUMR(1)= SQRT(UMR(1)) RUMR(2)= -0.5Q0*R(2)/SQRT(UMR(1)) ELSE IF(UMR(1).LT.0.Q0) THEN RUMR(1)= 0.5Q0*R(2)/SQRT(-UMR(1)) RUMR(2)= -SQRT(-UMR(1)) ENDIF DO I=1,2 RP(I)= RUMR(I)+RMR(I) RM(I)= RUMR(I)-RMR(I) ENDDO DO I=1,2 TRMR(I)= RP(I)-RM(I) TRUMR(I)= RP(I)+RM(I) ENDDO CALL CQLNX(RP,F) CALL CQLNX(TRUMR,H) CALL CQLNX(TRMR,G) * RM2(1)= RM(1)*RM(1)-RM(2)*RM(2) RM2(2)= 2.Q0*RM(1)*RM(2) ARM= RM(1) ARM2= ARM*ARM ARM4= ARM2*ARM2 ZRM= RM(2)/RM(1) ZRM2= ZRM*ZRM ZRM4= ZRM2*ZRM2 BRM= RM(2) BRM2= BRM*BRM BRM4= BRM2*BRM2 URM= RM(1)/RM(2) URM2= URM*URM URM4= URM2*URM2 IF(ABS(ZRM).LT.1.Q0) THEN RM4(1)= ARM4*(1.Q0-6.Q0*ZRM2+ZRM4) ELSE IF(ABS(URM).LT.1.Q0) THEN RM4(1)= BRM4*(1.Q0-6.Q0*URM2+URM4) ENDIF RM4(2)= 4.Q0*RM(1)*RM(2)*RM2(1) R2= R(1)*R(1)+R(2)*R(2) RI(1)= R(1)/R2 RI(2)= -R(2)/R2 RR2= R(1)*R(1)-R(2)*R(2) R22= RR2*RR2+4.Q0*R(1)*R(1)*R(2)*R(2) R2I(1)= RR2/R22 R2I(2)= -2.Q0*R(1)*R(2)/R22 X2= X(1)*X(1)+X(2)*X(2) XI(1)= X(1)/X2 XI(2)= -X(2)/X2 XX2= X(1)*X(1)-X(2)*X(2) X22= XX2*XX2+4.Q0*X(1)*X(1)*X(2)*X(2) X2I(1)= XX2/X22 X2I(2)= -2.Q0*X(1)*X(2)/X22 * RAR(1)= 1.Q0-R(1)/R2 RAR(2)= R(2)/R2 IF(RAR(1).GT.0.Q0) THEN RRAR(1)= SQRT(RAR(1)) RRAR(2)= 0.5Q0*RAR(2)/SQRT(RAR(1)) ELSE IF(RAR(1).LT.0.Q0) THEN RRAR(1)= 0.5Q0*RAR(2)/SQRT(-RAR(1)) RRAR(2)= SQRT(-RAR(1)) ENDIF * XP2= (1.Q0-X(1))*(1.Q0-X(1))+X(2)*X(2) XP(1)= (1.Q0-X(1))/XP2 XP(2)= X(2)/XP2 RPP2= (1.Q0-R(1))*(1.Q0-R(1))+R(2)*R(2) RPP(1)= (1.Q0-R(1))/RPP2 RPP(2)= R(2)/RPP2 * AR1R= RM2(1) AR1I= RM2(2) UMAR1R= 1.Q0-RM2(1) AR2R= RM4(1) AR2I= RM4(2) UMAR2R= 1.Q0-RM4(1) AR3R= XP(1) AR3I= XP(2) UMAR3R= 1.Q0-XP(1) CALL SPENCE(AR1R,AR1I,UMAR1R,CLI2SR,CLI2SI) CALL SPENCE(AR2R,AR2I,UMAR2R,CLI2FR,CLI2FI) CALL SPENCE(AR3R,AR3I,UMAR3R,CLI2XR,CLI2XI) CALL CLI3(AR1R,AR1I,CLI3SR,CLI3SI) CALL CLI3(AR2R,AR2I,CLI3FR,CLI3FI) CALL CLI3(AR3R,AR3I,CLI3XR,CLI3XI) * COMB3R= 2.Q0*CLI3SR-CLI3FR COMB3I= 2.Q0*CLI3SI-CLI3FI COMB2R= CLI2SR-CLI2FR COMB2I= CLI2SI-CLI2FI F2R= F(1)*F(1)-F(2)*F(2) F2I= 2.Q0*F(1)*F(2) * R1P= R(1)+1.5Q0 AUX1R= -F(1)+G(1)/3.Q0+2.Q0/3.Q0*H(1) AUX1I= -F(2)+G(2)/3.Q0+2.Q0/3.Q0*H(2) AUX2R= -3.Q0*F(1)+2.Q0*G(1)+4.Q0*H(1) AUX2I= -3.Q0*F(2)+2.Q0*G(2)+4.Q0*H(2) AUX3R= COMB2R+F(1)*AUX2R-F(2)*AUX2I AUX3PR= -2.Q0*(R1P*F(1)-R(2)*F(2)) AUX3I= COMB2I+F(1)*AUX2I+F(2)*AUX2R AUX3PI= -2.Q0*(R1P*F(2)+R(2)*F(1)) AUX4R= COMB2R+F(1)*AUX2R-F(2)*AUX2I AUX4PR= -2.Q0*((R(1)-3.Q0+0.25Q0*RI(1))*F(1)- # (R(2)+0.25Q0*RI(2))*F(2)) AUX4I= COMB2I+F(1)*AUX2I+F(2)*AUX2R AUX4PI= -2.Q0*((R(1)-3.Q0+0.25Q0*RI(1))*F(2)+ # (R(2)+0.25Q0*RI(2))*F(1)) AUX5R= R(1)-1.Q0/6.Q0-7.Q0/48.Q0*RI(1) AUX5I= R(2)-7.Q0/48.Q0*RI(2) AUX6R= R(1)-11.Q0/12.Q0+5.Q0/48.Q0*RI(1)+1.Q0/32.Q0* # R2I(1) AUX6I= R(2)+5.Q0/48.Q0*RI(2)+1.Q0/32.Q0*R2I(2) * EX0VR= 4.Q0*(R(1)-0.25Q0*RI(1)) EX0VI= 4.Q0*(R(2)-0.25Q0*RI(2)) * EX0AR= 4.Q0*(R(1)-1.5Q0+0.5Q0*RI(1)) EX0AI= 4.Q0*(R(2)+0.5Q0*RI(2)) * EX0XR= X(1)-1.5Q0+0.5Q0*X2I(1) EX0XI= X(2)+0.5Q0*X2I(2) * EX1R= COMB3R+8.Q0/3.Q0*(F(1)*COMB2R-F(2)*COMB2I)+ # 4.Q0*(F2R*AUX1R-F2I*AUX1I) EX1I= COMB3I+8.Q0/3.Q0*(F(1)*COMB2I+F(2)*COMB2R)+ # 4.Q0*(F2R*AUX1I+F2I*AUX1R) * EX2R= 8.Q0/3.Q0*((R(1)+0.5Q0)*AUX3R-R(2)*AUX3I)+AUX3PR EX2I= 8.Q0/3.Q0*((R(1)+0.5Q0)*AUX3I+R(2)*AUX3R)+AUX3PI * EX3R= -8.Q0*(AUX5R*F2R-AUX5I*F2I)+13.Q0/6.Q0+Z3*RI(1) EX3I= -8.Q0*(AUX5R*F2I+AUX5I*F2R)+Z3*RI(2) * EX4R= 8.Q0/3.Q0*((R(1)-1.Q0)*AUX4R-R(2)*AUX4I)+AUX4PR EX4I= 8.Q0/3.Q0*((R(1)-1.Q0)*AUX4I+R(2)*AUX4R)+AUX4PI * EX5R= -8.Q0*(AUX6R*F2R-AUX6I*F2I)+ # 13.Q0/6.Q0-3.Q0*Z2+(0.25Q0-2.Q0*Z3)*RI(1) EX5I= -8.Q0*(AUX6R*F2I+AUX6I*F2R)+(0.25Q0-2.Q0*Z3)* # RI(2) * EX6R= CLI3XR+2.Q0/3.Q0*(B(1)*CLI2XR-B(2)*CLI2XI)-1.Q0/6.Q0* # (B2(1)*(A(1)-B(1))-B2(2)*(A(2)-B(2))) EX6I= CLI3XI+2.Q0/3.Q0*(B(1)*CLI2XI+B(2)*CLI2XR)-1.Q0/6.Q0* # (B2(1)*(A(2)-B(2))+B2(2)*(A(1)-B(1))) * EX7R= 1.Q0/3.Q0*((X(1)+0.5Q0-0.5Q0*XI(1))*(CLI2XR-AB(1))- # (X(2)-0.5Q0*XI(2))*(CLI2XI-AB(2)))+1.Q0/3.Q0*(B2(1)* # (X(1)-1.Q0/8.Q0-XI(1)+5.Q0/8.Q0*X2I(1))-B2(2)* # (X(2)-XI(2)+5.Q0/8.Q0*X2I(2)))-0.25Q0*(B(1)* # (X(1)-5.Q0/2.Q0+2.Q0/3.Q0*XI(1)+5.Q0/6.Q0*X2I(1))- # B(2)*(X(2)+2.Q0/3.Q0*XI(2)+5.Q0/6.Q0*X2I(2))) EX7I= 1.Q0/3.Q0*((X(1)+0.5Q0-0.5Q0*XI(1))*(CLI2XI-AB(2))+ # (X(2)-0.5Q0*XI(2))*(CLI2XR-AB(1)))+1.Q0/3.Q0*(B2(1)* # (X(2)-XI(2)+5.Q0/8.Q0*X2I(2))+B2(2)*(X(1)-1.Q0/8.Q0- # XI(1)+5.Q0/8.Q0*X2I(1)))-0.25Q0*(B(1)*(X(2)+2.Q0/3.Q0* # XI(2)+5.Q0/6.Q0*X2I(2))+B(2)*(X(1)-5.Q0/2.Q0+2.Q0/3.Q0* # XI(1)+5.Q0/6.Q0*X2I(1))) * EX8R= -3.Q0/4.Q0*Z2+13.Q0/12.Q0-5.Q0/24.Q0*XI(1)-0.5Q0*Z3*X2I(1) EX8I= -5.Q0/24.Q0*XI(2)-0.5Q0*Z3*X2I(2) * V1R= EX0VR*EX1R-EX0VI*EX1I+RRAR(1)*EX2R-RRAR(2)*EX2I+EX3R V1I= EX0VR*EX1I+EX0VI*EX1R+RRAR(1)*EX2I+RRAR(2)*EX2R+EX3I A1R= EX0AR*EX1R-EX0AI*EX1I+RRAR(1)*EX4R-RRAR(2)*EX4I+EX5R A1I= EX0AR*EX1I+EX0AI*EX1R+RRAR(1)*EX4I+RRAR(2)*EX4R+EX5I F1X= EX0XR*EX6R-EX0XI*EX6I+EX7R+EX8R F1XI= EX0XR*EX6I+EX0XI*EX6R+EX7I+EX8I * EX0VPR= 4.Q0*(1.Q0+0.25Q0*R2I(1)) EX0VPI= R2I(2) EX0APR= 4.Q0*(1.Q0-0.5Q0*R2I(1)) EX0API= -2.Q0*R2I(2) * BUX1R= -COMB2R-F(1)*AUX2R+F(2)*AUX2I BUX1I= -COMB2I-F(1)*AUX2I-F(2)*AUX2R BUX2R= 2.Q0/3.Q0*(RI(1)*BUX1R-RI(2)*BUX1I) BUX2I= 2.Q0/3.Q0*(RI(1)*BUX1I+RI(2)*BUX1R) BUX3R= BUX2R-(1.Q0-5.Q0/6.Q0*RI(1))*F(1)- # 5.Q0/6.Q0*RI(2)*F(2) BUX3I= BUX2I-(1.Q0-5.Q0/6.Q0*RI(1))*F(2)+ # 5.Q0/6.Q0*RI(2)*F(1) BUX4R= F2R*(4.Q0/3.Q0*RI(1)-11.Q0/6.Q0*R2I(1)- # 4.Q0*RPP(1))-F2I*(4.Q0/3.Q0*RI(2)-11.Q0/6.Q0* # R2I(2)-4.Q0*RPP(2))-1.Q0-1.5Q0*RI(1)-Z3*R2I(1) BUX4I= F2R*(4.Q0/3.Q0*RI(2)-11.Q0/6.Q0*R2I(2)- # 4.Q0*RPP(2))+F2I*(4.Q0/3.Q0*RI(1)-11.Q0/6.Q0* # R2I(1)-4.Q0*RPP(1))-1.5Q0*RI(2)-Z3*R2I(2) CUX2R= 4.Q0/3.Q0*(-RI(1)*BUX1R+RI(2)*BUX1I) CUX2I= 4.Q0/3.Q0*(-RI(1)*BUX1I-RI(2)*BUX1R) CUX3R= CUX2R-(1.Q0+13.Q0/6.Q0*RI(1)-0.5Q0*R2I(1))*F(1)+ # (13.Q0/6.Q0*RI(2)-0.5Q0*R2I(2))*F(2) CUX3I= CUX2I-(1.Q0+13.Q0/6.Q0*RI(1)-0.5Q0*R2I(1))*F(2)- # (13.Q0/6.Q0*RI(2)-0.5Q0*R2I(2))*F(1) CUX4R= F2R*RI(1)-F2I*RI(2) CUX4I= F2R*RI(2)+F2I*RI(1) CUX5R= CUX4R*(-20.Q0/3.Q0+13.Q0/6.Q0*RI(1)+0.5Q0*R2I(1))- # CUX4I*(13.Q0/6.Q0*RI(2)+0.5Q0*R2I(2))-1.Q0+3.Q0*RI(1)+ # (2.Q0*Z3-0.5Q0)*R2I(1) CUX5I= CUX4R*(13.Q0/6.Q0*RI(2)+0.5Q0*R2I(2))+ # CUX4I*(20.Q0/3.Q0+13.Q0/6.Q0*RI(1)+0.5Q0*R2I(1))+ # 3.Q0*RI(2)+(2.Q0*Z3-0.5Q0)*R2I(2) * V1PR= EX0VPR*EX1R-EX0VPI*EX1I+2.Q0*(RRAR(1)*BUX3R- # RRAR(2)*BUX3I)+BUX4R A1PR= EX0APR*EX1R-EX0API*EX1I+2.Q0*(RRAR(1)*CUX3R- # RRAR(2)*CUX3I)+CUX5R * V1= V1R A1= A1R F1= F1X V1P= V1PR A1P= A1PR * RETURN END * *-----CLI3 * COMPUTES LI_3(X) FOR COMPLEX X * SUBROUTINE CLI3(XR,XI,CLI3R,CLI3I) IMPLICIT REAL*16(A-H,O-Z) * COMMON/QPARAM/QPI,QPIS,QEPS,QDELTA * DIMENSION B(0:14),BF(0:14) DIMENSION X(2),Y(2),ADDX(2),OX(2),CLNX(2),PAR(2),CT(15),SN(15), # RES(2),U1(2),U2(2),CLNY(2),OMY(2),CLNOMY(2),ADDY(2), # PAR1(2),PAR2(2),CT1(15),SN1(15),CT2(15),SN2(15), # RES1(2),RES2(2),T(2),RESA(2),RESB(2),CLNT(2), # RES3(2),RES4(2),CLNOMT(2),ADDT(2),ADDT2(2),OMT(2), # OMU1(2),OMU2(2) * DATA B/1.Q0,-0.75Q0,0.236111111111111111111111111111111Q0, # -3.472222222222222222222222222222222Q-2, # 6.481481481481481481481481481481482Q-4, # 4.861111111111111111111111111111111Q-4, # -2.393550012597631645250692869740488Q-5, # -1.062925170068027210884353741496599Q-5, # 7.794784580498866213151927437641723Q-7, # 2.526087595532039976484420928865373Q-7, # -2.359163915200471237027273583310139Q-8, # -6.168132746415574698402981231264060Q-9, # 6.824456748981078267312315451125495Q-10, # 1.524285616929084572552216019859487Q-10, # -1.916909414174054295837274763110831Q-11/ DATA Z3/1.20205690315959428540Q0/ * Z2= QPIS/6.Q0 DO N=0,14 BF(N)= B(N)/(N+1.Q0) ENDDO * X(1)= XR X(2)= XI * XM2= X(1)*X(1)+X(2)*X(2) XM= SQRT(XM2) * *-----THE MODULUS OF X IS CHECKED * XTST= XM-1.Q0 IF(XTST.LE.1.Q-33) THEN Y(1)= X(1) Y(2)= X(2) ADDX(1)= 0.Q0 ADDX(2)= 0.Q0 ELSE IF(XM.GT.1.Q-33) THEN Y(1)= X(1)/XM2 Y(2)= -X(2)/XM2 OX(1)= -X(1) OX(2)= -X(2) CALL CQLNX(OX,CLNX) RLNX2= CLNX(1)*CLNX(1) AILNX2= CLNX(2)*CLNX(2) ADDX(1)= -CLNX(1)*(Z2+1.Q0/6.Q0*(RLNX2-3.Q0*AILNX2)) ADDX(2)= -CLNX(2)*(Z2+1.Q0/6.Q0*(3.Q0*RLNX2-AILNX2)) ENDIF * *-----ONCE X --> Y, |Y|<1 THE SIGN OF RE(Y) IS CHECKED * IF RE(Y)>0 A TRANSFORMATION IS REQUIRED FOR RE(Y)>1/2 * Y2R= Y(1)*Y(1)-Y(2)*Y(2) IF(Y(1).GE.0.Q0.OR.Y2R.LT.0.Q0) THEN YTST= Y(1)-0.5Q0 IF(YTST.LE.1.Q-33) THEN * *-----LI_3(Y) IS COMPUTED * OMY(1)= 1.Q0-Y(1) OMY(2)= -Y(2) CALL CQLNOMX(Y,OMY,PAR) PR= -PAR(1) PI= -PAR(2) P2= PR*PR+PI*PI PM= SQRT(P2) CT(1)= PR/PM SN(1)= PI/PM DO N=2,15 CT(N)= CT(1)*CT(N-1)-SN(1)*SN(N-1) SN(N)= SN(1)*CT(N-1)+CT(1)*SN(N-1) ENDDO RES(1)= PM*(BF(0)*CT(1)+PM*(BF(1)*CT(2)+PM* # (BF(2)*CT(3)+PM*(BF(3)*CT(4)+PM* # (BF(4)*CT(5)+PM*(BF(5)*CT(6)+PM* # (BF(6)*CT(7)+PM*(BF(7)*CT(8)+PM* # (BF(8)*CT(9)+PM*(BF(9)*CT(10)+PM* # (BF(10)*CT(11)+PM*(BF(11)*CT(12)+PM* # (BF(12)*CT(13)+PM*(BF(13)*CT(14)+PM* # (BF(14)*CT(15)))))))))))))))) RES(2)= PM*(BF(0)*SN(1)+PM*(BF(1)*SN(2)+PM* # (BF(2)*SN(3)+PM*(BF(3)*SN(4)+PM* # (BF(4)*SN(5)+PM*(BF(5)*SN(6)+PM* # (BF(6)*SN(7)+PM*(BF(7)*SN(8)+PM* # (BF(8)*SN(9)+PM*(BF(9)*SN(10)+PM* # (BF(10)*SN(11)+PM*(BF(11)*SN(12)+PM* # (BF(12)*SN(13)+PM*(BF(13)*SN(14)+PM* # (BF(14)*SN(15)))))))))))))))) CLI3R= RES(1)+ADDX(1) CLI3I= RES(2)+ADDX(2) RETURN ELSE IF(YTST.GT.1.Q-33) THEN YM2= Y(1)*Y(1)+Y(2)*Y(2) U1(1)= 1.Q0-Y(1)/YM2 U1(2)= Y(2)/YM2 U2(1)= 1.Q0-Y(1) U2(2)= -Y(2) CALL CQLNX(Y,CLNY) OMY(1)= 1.Q0-Y(1) OMY(2)= -Y(2) CALL CQLNOMX(Y,OMY,CLNOMY) ADDY(1)= Z3+Z2*CLNY(1)+1.Q0/6.D0*CLNY(1)* # (CLNY(1)*CLNY(1)-3.Q0*CLNY(2)*CLNY(2))- # 0.5Q0*CLNOMY(1)*(CLNY(1)*CLNY(1)-CLNY(2)* # CLNY(2))+CLNY(1)*CLNY(2)*CLNOMY(2) ADDY(2)= Z2*CLNY(2)+1.Q0/6.Q0*CLNY(2)*(3.Q0* # CLNY(1)*CLNY(1)-CLNY(2)*CLNY(2))-0.5Q0* # CLNOMY(2)*(CLNY(1)*CLNY(1)-CLNY(2)*CLNY(2))- # CLNY(1)*CLNOMY(1)*CLNY(2) * *-----LI_3(1-1/Y) IS COMPUTED * OMU1(1)= 1.Q0-U1(1) OMU1(2)= -U1(2) CALL CQLNOMX(U1,OMU1,PAR1) PR1= -PAR1(1) PI1= -PAR1(2) P12= PR1*PR1+PI1*PI1 PM1= SQRT(P12) CT1(1)= PR1/PM1 SN1(1)= PI1/PM1 DO N=2,15 CT1(N)= CT1(1)*CT1(N-1)-SN1(1)*SN1(N-1) SN1(N)= SN1(1)*CT1(N-1)+CT1(1)*SN1(N-1) ENDDO RES1(1)= PM1*(BF(0)*CT1(1)+PM1*(BF(1)*CT1(2)+PM1* # (BF(2)*CT1(3)+PM1*(BF(3)*CT1(4)+PM1* # (BF(4)*CT1(5)+PM1*(BF(5)*CT1(6)+PM1* # (BF(6)*CT1(7)+PM1*(BF(7)*CT1(8)+PM1* # (BF(8)*CT1(9)+PM1*(BF(9)*CT1(10)+PM1* # (BF(10)*CT1(11)+PM1*(BF(11)*CT1(12)+PM1* # (BF(12)*CT1(13)+PM1*(BF(13)*CT1(14)+PM1* # (BF(14)*CT1(15)))))))))))))))) RES1(2)= PM1*(BF(0)*SN1(1)+PM1*(BF(1)*SN1(2)+PM1* # (BF(2)*SN1(3)+PM1*(BF(3)*SN1(4)+PM1* # (BF(4)*SN1(5)+PM1*(BF(5)*SN1(6)+PM1* # (BF(6)*SN1(7)+PM1*(BF(7)*SN1(8)+PM1* # (BF(8)*SN1(9)+PM1*(BF(9)*SN1(10)+PM1* # (BF(10)*SN1(11)+PM1*(BF(11)*SN1(12)+PM1* # (BF(12)*SN1(13)+PM1*(BF(13)*SN1(14)+PM1* # (BF(14)*SN1(15)))))))))))))))) * *-----LI_3(1-Y) IS COMPUTED * OMU2(1)= 1.Q0-U2(1) OMU2(2)= -U2(2) CALL CQLNOMX(U2,OMU2,PAR2) PR2= -PAR2(1) PI2= -PAR2(2) P22= PR2*PR2+PI2*PI2 PM2= SQRT(P22) CT2(1)= PR2/PM2 SN2(1)= PI2/PM2 DO N=2,15 CT2(N)= CT2(1)*CT2(N-1)-SN2(1)*SN2(N-1) SN2(N)= SN2(1)*CT2(N-1)+CT2(1)*SN2(N-1) ENDDO RES2(1)= PM2*(BF(0)*CT2(1)+PM2*(BF(1)*CT2(2)+PM2* # (BF(2)*CT2(3)+PM2*(BF(3)*CT2(4)+PM2* # (BF(4)*CT2(5)+PM2*(BF(5)*CT2(6)+PM2* # (BF(6)*CT2(7)+PM2*(BF(7)*CT2(8)+PM2* # (BF(8)*CT2(9)+PM2*(BF(9)*CT2(10)+PM2* # (BF(10)*CT2(11)+PM2*(BF(11)*CT2(12)+PM2* # (BF(12)*CT2(13)+PM2*(BF(13)*CT2(14)+PM2* # (BF(14)*CT2(15)))))))))))))))) RES2(2)= PM2*(BF(0)*SN2(1)+PM2*(BF(1)*SN2(2)+PM2* # (BF(2)*SN2(3)+PM2*(BF(3)*SN2(4)+PM2* # (BF(4)*SN2(5)+PM2*(BF(5)*SN2(6)+PM2* # (BF(6)*SN2(7)+PM2*(BF(7)*SN2(8)+PM2* # (BF(8)*SN2(9)+PM2*(BF(9)*SN2(10)+PM2* # (BF(10)*SN2(11)+PM2*(BF(11)*SN2(12)+PM2* # (BF(12)*SN2(13)+PM2*(BF(13)*SN2(14)+PM2* # (BF(14)*SN2(15)))))))))))))))) CLI3R= -RES1(1)-RES2(1)+ADDX(1)+ADDY(1) CLI3I= -RES1(2)-RES2(2)+ADDX(2)+ADDY(2) RETURN ENDIF * *-----IF RE(Y)<0 A TRANSFORMATION IS REQUIRED IN TERMS OF T = -Y * AND OF T^2 * ELSE IF(Y(1).LT.0.Q0) THEN * *-----FIRST T * T(1)= -Y(1) T(2)= -Y(2) IF(T(1).LE.0.5Q0) THEN * *-----LI_3(T) IS COMPUTED * OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) CALL CQLNOMX(T,OMT,PAR) PR= -PAR(1) PI= -PAR(2) P2= PR*PR+PI*PI PM= SQRT(P2) CT(1)= PR/PM SN(1)= PI/PM DO N=2,15 CT(N)= CT(1)*CT(N-1)-SN(1)*SN(N-1) SN(N)= SN(1)*CT(N-1)+CT(1)*SN(N-1) ENDDO RESA(1)= PM*(BF(0)*CT(1)+PM*(BF(1)*CT(2)+PM* # (BF(2)*CT(3)+PM*(BF(3)*CT(4)+PM* # (BF(4)*CT(5)+PM*(BF(5)*CT(6)+PM* # (BF(6)*CT(7)+PM*(BF(7)*CT(8)+PM* # (BF(8)*CT(9)+PM*(BF(9)*CT(10)+PM* # (BF(10)*CT(11)+PM*(BF(11)*CT(12)+PM* # (BF(12)*CT(13)+PM*(BF(13)*CT(14)+PM* # (BF(14)*CT(15)))))))))))))))) RESA(2)= PM*(BF(0)*SN(1)+PM*(BF(1)*SN(2)+PM* # (BF(2)*SN(3)+PM*(BF(3)*SN(4)+PM* # (BF(4)*SN(5)+PM*(BF(5)*SN(6)+PM* # (BF(6)*SN(7)+PM*(BF(7)*SN(8)+PM* # (BF(8)*SN(9)+PM*(BF(9)*SN(10)+PM* # (BF(10)*SN(11)+PM*(BF(11)*SN(12)+PM* # (BF(12)*SN(13)+PM*(BF(13)*SN(14)+PM* # (BF(14)*SN(15)))))))))))))))) ELSE IF(T(1).GT.0.5Q0) THEN TM2= T(1)*T(1)+T(2)*T(2) U1(1)= 1.Q0-T(1)/TM2 U1(2)= T(2)/TM2 U2(1)= 1.Q0-T(1) U2(2)= -T(2) CALL CQLNX(T,CLNT) OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) CALL CQLNOMX(T,OMT,CLNOMT) ADDT(1)= Z3+Z2*CLNT(1)+1.Q0/6.D0*CLNT(1)* # (CLNT(1)*CLNT(1)-3.Q0*CLNT(2)*CLNT(2))- # 0.5Q0*CLNOMT(1)*(CLNT(1)*CLNT(1)-CLNT(2)* # CLNT(2))+CLNT(1)*CLNT(2)*CLNOMT(2) ADDT(2)= Z2*CLNT(2)+1.Q0/6.Q0*CLNT(2)*(3.Q0* # CLNT(1)*CLNT(1)-CLNT(2)*CLNT(2))-0.5Q0* # CLNOMT(2)*(CLNT(1)*CLNT(1)-CLNT(2)*CLNT(2))- # CLNT(1)*CLNOMT(1)*CLNT(2) * *-----LI3(1-1/T) IS COMPUTED * OMU1(1)= 1.Q0-U1(1) OMU1(2)= -U1(2) CALL CQLNOMX(U1,OMU1,PAR1) PR1= -PAR1(1) PI1= -PAR1(2) P12= PR1*PR1+PI1*PI1 PM1= SQRT(P12) CT1(1)= PR1/PM1 SN1(1)= PI1/PM1 DO N=2,15 CT1(N)= CT1(1)*CT1(N-1)-SN1(1)*SN1(N-1) SN1(N)= SN1(1)*CT1(N-1)+CT1(1)*SN1(N-1) ENDDO RES1(1)= PM1*(BF(0)*CT1(1)+PM1*(BF(1)*CT1(2)+PM1* # (BF(2)*CT1(3)+PM1*(BF(3)*CT1(4)+PM1* # (BF(4)*CT1(5)+PM1*(BF(5)*CT1(6)+PM1* # (BF(6)*CT1(7)+PM1*(BF(7)*CT1(8)+PM1* # (BF(8)*CT1(9)+PM1*(BF(9)*CT1(10)+PM1* # (BF(10)*CT1(11)+PM1*(BF(11)*CT1(12)+PM1* # (BF(12)*CT1(13)+PM1*(BF(13)*CT1(14)+PM1* # (BF(14)*CT1(15)))))))))))))))) RES1(2)= PM1*(BF(0)*SN1(1)+PM1*(BF(1)*SN1(2)+PM1* # (BF(2)*SN1(3)+PM1*(BF(3)*SN1(4)+PM1* # (BF(4)*SN1(5)+PM1*(BF(5)*SN1(6)+PM1* # (BF(6)*SN1(7)+PM1*(BF(7)*SN1(8)+PM1* # (BF(8)*SN1(9)+PM1*(BF(9)*SN1(10)+PM1* # (BF(10)*SN1(11)+PM1*(BF(11)*SN1(12)+PM1* # (BF(12)*SN1(13)+PM1*(BF(13)*SN1(14)+PM1* # (BF(14)*SN1(15)))))))))))))))) * *-----LI3(1-T) IS COMPUTED * OMU2(1)= 1.Q0-U2(1) OMU2(2)= -U2(2) CALL CQLNOMX(U2,OMU2,PAR2) PR2= -PAR2(1) PI2= -PAR2(2) P22= PR2*PR2+PI2*PI2 PM2= SQRT(P22) CT2(1)= PR2/PM2 SN2(1)= PI2/PM2 DO N=2,15 CT2(N)= CT2(1)*CT2(N-1)-SN2(1)*SN2(N-1) SN2(N)= SN2(1)*CT2(N-1)+CT2(1)*SN2(N-1) ENDDO RES2(1)= PM2*(BF(0)*CT2(1)+PM2*(BF(1)*CT2(2)+PM2* # (BF(2)*CT2(3)+PM2*(BF(3)*CT2(4)+PM2* # (BF(4)*CT2(5)+PM2*(BF(5)*CT2(6)+PM2* # (BF(6)*CT2(7)+PM2*(BF(7)*CT2(8)+PM2* # (BF(8)*CT2(9)+PM2*(BF(9)*CT2(10)+PM2* # (BF(10)*CT2(11)+PM2*(BF(11)*CT2(12)+PM2* # (BF(12)*CT2(13)+PM2*(BF(13)*CT2(14)+PM2* # (BF(14)*CT2(15)))))))))))))))) RES2(2)= PM2*(BF(0)*SN2(1)+PM2*(BF(1)*SN2(2)+PM2* # (BF(2)*SN2(3)+PM2*(BF(3)*SN2(4)+PM2* # (BF(4)*SN2(5)+PM2*(BF(5)*SN2(6)+PM2* # (BF(6)*SN2(7)+PM2*(BF(7)*SN2(8)+PM2* # (BF(8)*SN2(9)+PM2*(BF(9)*SN2(10)+PM2* # (BF(10)*SN2(11)+PM2*(BF(11)*SN2(12)+PM2* # (BF(12)*SN2(13)+PM2*(BF(13)*SN2(14)+PM2* # (BF(14)*SN2(15)))))))))))))))) RESA(1)= -RES1(1)-RES2(1)+ADDT(1) RESA(2)= -RES1(2)-RES2(2)+ADDT(2) ENDIF * *-----THEN T^2 * T(1)= Y(1)*Y(1)-Y(2)*Y(2) T(2)= 2.Q0*Y(1)*Y(2) IF(T(1).LE.0.5Q0) THEN * *-----LI_3(T^2) IS COMPUTED * OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) CALL CQLNOMX(T,OMT,PAR) PR= -PAR(1) PI= -PAR(2) P2= PR*PR+PI*PI PM= SQRT(P2) CT(1)= PR/PM SN(1)= PI/PM DO N=2,15 CT(N)= CT(1)*CT(N-1)-SN(1)*SN(N-1) SN(N)= SN(1)*CT(N-1)+CT(1)*SN(N-1) ENDDO RESB(1)= PM*(BF(0)*CT(1)+PM*(BF(1)*CT(2)+PM* # (BF(2)*CT(3)+PM*(BF(3)*CT(4)+PM* # (BF(4)*CT(5)+PM*(BF(5)*CT(6)+PM* # (BF(6)*CT(7)+PM*(BF(7)*CT(8)+PM* # (BF(8)*CT(9)+PM*(BF(9)*CT(10)+PM* # (BF(10)*CT(11)+PM*(BF(11)*CT(12)+PM* # (BF(12)*CT(13)+PM*(BF(13)*CT(14)+PM* # (BF(14)*CT(15)))))))))))))))) RESB(2)= PM*(BF(0)*SN(1)+PM*(BF(1)*SN(2)+PM* # (BF(2)*SN(3)+PM*(BF(3)*SN(4)+PM* # (BF(4)*SN(5)+PM*(BF(5)*SN(6)+PM* # (BF(6)*SN(7)+PM*(BF(7)*SN(8)+PM* # (BF(8)*SN(9)+PM*(BF(9)*SN(10)+PM* # (BF(10)*SN(11)+PM*(BF(11)*SN(12)+PM* # (BF(12)*SN(13)+PM*(BF(13)*SN(14)+PM* # (BF(14)*SN(15)))))))))))))))) ELSE IF(T(1).GT.0.5Q0) THEN TM2= T(1)*T(1)+T(2)*T(2) U1(1)= 1.Q0-T(1)/TM2 U1(2)= T(2)/TM2 U2(1)= 1.Q0-T(1) U2(2)= -T(2) CALL CQLNX(T,CLNT) OMT(1)= 1.Q0-T(1) OMT(2)= -T(2) CALL CQLNOMX(T,OMT,CLNOMT) ADDT2(1)= Z3+Z2*CLNT(1)+1.Q0/6.D0*CLNT(1)* # (CLNT(1)*CLNT(1)-3.Q0*CLNT(2)*CLNT(2))- # 0.5Q0*CLNOMT(1)*(CLNT(1)*CLNT(1)-CLNT(2)* # CLNT(2))+CLNT(1)*CLNT(2)*CLNOMT(2) ADDT2(2)= Z2*CLNT(2)+1.Q0/6.Q0*CLNT(2)*(3.Q0* # CLNT(1)*CLNT(1)-CLNT(2)*CLNT(2))-0.5Q0* # CLNOMT(2)*(CLNT(1)*CLNT(1)-CLNT(2)*CLNT(2))- # CLNT(1)*CLNOMT(1)*CLNT(2) * *-----LI_3(1-1/T^2) IS COMPUTED * OMU1(1)= 1.Q0-U1(1) OMU1(2)= -U1(2) CALL CQLNOMX(U1,OMU1,PAR1) PR1= -PAR1(1) PI1= -PAR1(2) P12= PR1*PR1+PI1*PI1 PM1= SQRT(P12) CT1(1)= PR1/PM1 SN1(1)= PI1/PM1 DO N=2,15 CT1(N)= CT1(1)*CT1(N-1)-SN1(1)*SN1(N-1) SN1(N)= SN1(1)*CT1(N-1)+CT1(1)*SN1(N-1) ENDDO RES3(1)= PM1*(BF(0)*CT1(1)+PM1*(BF(1)*CT1(2)+PM1* # (BF(2)*CT1(3)+PM1*(BF(3)*CT1(4)+PM1* # (BF(4)*CT1(5)+PM1*(BF(5)*CT1(6)+PM1* # (BF(6)*CT1(7)+PM1*(BF(7)*CT1(8)+PM1* # (BF(8)*CT1(9)+PM1*(BF(9)*CT1(10)+PM1* # (BF(10)*CT1(11)+PM1*(BF(11)*CT1(12)+PM1* # (BF(12)*CT1(13)+PM1*(BF(13)*CT1(14)+PM1* # (BF(14)*CT1(15)))))))))))))))) RES3(2)= PM1*(BF(0)*SN1(1)+PM1*(BF(1)*SN1(2)+PM1* # (BF(2)*SN1(3)+PM1*(BF(3)*SN1(4)+PM1* # (BF(4)*SN1(5)+PM1*(BF(5)*SN1(6)+PM1* # (BF(6)*SN1(7)+PM1*(BF(7)*SN1(8)+PM1* # (BF(8)*SN1(9)+PM1*(BF(9)*SN1(10)+PM1* # (BF(10)*SN1(11)+PM1*(BF(11)*SN1(12)+PM1* # (BF(12)*SN1(13)+PM1*(BF(13)*SN1(14)+PM1* # (BF(14)*SN1(15)))))))))))))))) * *-----LI_3(1-T^2) IS COMPUTED * OMU2(1)= 1.Q0-U2(1) OMU2(2)= -U2(2) CALL CQLNOMX(U2,OMU2,PAR2) PR2= -PAR2(1) PI2= -PAR2(2) P22= PR2*PR2+PI2*PI2 PM2= SQRT(P22) CT2(1)= PR2/PM2 SN2(1)= PI2/PM2 DO N=2,15 CT2(N)= CT2(1)*CT2(N-1)-SN2(1)*SN2(N-1) SN2(N)= SN2(1)*CT2(N-1)+CT2(1)*SN2(N-1) ENDDO RES4(1)= PM2*(BF(0)*CT2(1)+PM2*(BF(1)*CT2(2)+PM2* # (BF(2)*CT2(3)+PM2*(BF(3)*CT2(4)+PM2* # (BF(4)*CT2(5)+PM2*(BF(5)*CT2(6)+PM2* # (BF(6)*CT2(7)+PM2*(BF(7)*CT2(8)+PM2* # (BF(8)*CT2(9)+PM2*(BF(9)*CT2(10)+PM2* # (BF(10)*CT2(11)+PM2*(BF(11)*CT2(12)+PM2* # (BF(12)*CT2(13)+PM2*(BF(13)*CT2(14)+PM2* # (BF(14)*CT2(15)))))))))))))))) RES4(2)= PM2*(BF(0)*SN2(1)+PM2*(BF(1)*SN2(2)+PM2* # (BF(2)*SN2(3)+PM2*(BF(3)*SN2(4)+PM2* # (BF(4)*SN2(5)+PM2*(BF(5)*SN2(6)+PM2* # (BF(6)*SN2(7)+PM2*(BF(7)*SN2(8)+PM2* # (BF(8)*SN2(9)+PM2*(BF(9)*SN2(10)+PM2* # (BF(10)*SN2(11)+PM2*(BF(11)*SN2(12)+PM2* # (BF(12)*SN2(13)+PM2*(BF(13)*SN2(14)+PM2* # (BF(14)*SN2(15)))))))))))))))) RESB(1)= -RES3(1)-RES4(1)+ADDT2(1) RESB(2)= -RES3(2)-RES4(2)+ADDT2(2) ENDIF CLI3R= -RESA(1)+0.25Q0*RESB(1)+ADDX(1) CLI3I= -RESA(2)+0.25Q0*RESB(2)+ADDX(2) RETURN ENDIF END