bin/BCVEGPY/p_lorentz.F

0001 c************************************************************

0002 c...some subroutines for p-wave calculation.

0003 c************************************************************
0004 c*** genpolar is used to generate the polarization vector ***
0005 c*** ep of the vector boson with momentum p and mass pm.  ***
0006 c************************************************************

0007       subroutine genpolar3(p,ep,pm)
0008       implicit double precision(a-h, o-z)
0009         implicit integer(i-n)
0010       dimension p(4),ep(4,2:4)
0011
0012       pxy=dsqrt(p(2)**2+p(3)**2)
0013       pxyz=dsqrt(p(2)**2+p(3)**2+p(4)**2)
0014
0015       if(pxy.lt.1.0d-16) then
0016          if(p(4).gt.0.0d0) then
0017            sign=1.0d0
0018          else
0019            sign=-1.0d0
0020          endif
0021 c******    ploarization vector 1 of momentum p *****
0022          ep(1,2)=0.0d0
0023          ep(2,2)=sign
0024          ep(3,2)=0.0d0
0025          ep(4,2)=0.0d0
0026 c******    ploarization vector 2 of momentum p *****
0027          ep(1,3)=0.0d0
0028          ep(2,3)=0.0d0
0029          ep(3,3)=sign
0030          ep(4,3)=0.0d0
0031 c******    ploarization vector 3 of momentum p *****
0032          ep(1,4)=p(4)/pm
0033          ep(2,4)=0.0d0
0034          ep(3,4)=0.0d0
0035          ep(4,4)=p(1)/pm
0036       else
0037 c******    ploarization vector 1 of momentum p *****
0038          ep(1,2)=0.0d0
0039          ep(2,2)=-p(3)/pxy
0040          ep(3,2)=p(2)/pxy
0041          ep(4,2)=0.0d0
0042 c******    ploarization vector 2 of momentum p *****
0043          ep(1,3)=0.0d0
0044          ep(2,3)=(p(2)*p(4))/(pxyz*pxy)
0045          ep(3,3)=(p(3)*p(4))/(pxyz*pxy)
0046          ep(4,3)=-pxy/pxyz
0047
0048 c******    ploarization vector 3 of momentum p *****
0049          ep(1,4)=pxyz/pm
0050          ep(2,4)=p(2)*p(1)/(pxyz*pm)
0051          ep(3,4)=p(3)*p(1)/(pxyz*pm)
0052          ep(4,4)=p(4)*p(1)/(pxyz*pm)
0053       endif
0054
0055       return
0056       end
0057 c*********************************************
0058       subroutine genpolar(p,ep)
0059       implicit double precision(a-h, o-z)
0060         implicit integer(i-n)
0061       dimension p(4),ep(4,2:3)
0062
0063       pxy=dsqrt(p(2)**2+p(3)**2)
0064       pxyz=dsqrt(p(2)**2+p(3)**2+p(4)**2)
0065
0066       if (pxy.lt.1.0d-16) then
0067          if(p(4).gt.0.0d0) then
0068            sign=1.0d0
0069          else
0070            sign=-1.0d0
0071          endif
0072 c******    ploarization vector 1 of momentum p *****
0073          ep(1,2)=0.0d0
0074          ep(2,2)=sign
0075          ep(3,2)=0.0d0
0076          ep(4,2)=0.0d0
0077 c******    ploarization vector 2 of momentum p *****
0078          ep(1,3)=0.0d0
0079          ep(2,3)=0.0d0
0080          ep(3,3)=sign
0081          ep(4,3)=0.0d0
0082       else
0083 c******    ploarization vector 1 of momentum p *****
0084          ep(1,2)=0.0d0
0085          ep(2,2)=-p(3)/pxy
0086          ep(3,2)=p(2)/pxy
0087          ep(4,2)=0.0d0
0088 c******    ploarization vector 2 of momentum p *****
0089          ep(1,3)=0.0d0
0090          ep(2,3)=-(p(2)*p(4))/(pxyz*pxy)
0091          ep(3,3)=-(p(3)*p(4))/(pxyz*pxy)
0092          ep(4,3)=pxy/pxyz
0093       endif
0094
0095       return
0096       end
0097 c*************************************************************

0098 c*** gentensor is used to generate the polarization tensor ***
0099 c*** ept with momentum p.                                  ***
0100 c...this is for the general results (m > 0).
0101       subroutine gentensor(p,ept,pm)
0102       implicit double precision(a-h, o-z)
0103         implicit integer(i-n)
0104       dimension p(4),ept(4,4,2:6),ep(4,2:4)
0105
0106       call genpolar3(p,ep,pm)
0107
0108       do i=1,4
0109        do j=1,4
0110         ept(i,j,2)=1.0d0/dsqrt(2.0d0)*(ep(i,2)*ep(j,3)+ep(i,3)*ep(j,2))
0111         ept(i,j,3)=1.0d0/dsqrt(2.0d0)*(ep(i,2)*ep(j,4)+ep(i,4)*ep(j,2))
0112         ept(i,j,4)=1.0d0/dsqrt(2.0d0)*(ep(i,4)*ep(j,3)+ep(i,3)*ep(j,4))
0113         ept(i,j,5)=1.0d0/dsqrt(2.0d0)*(ep(i,2)*ep(j,2)-ep(i,3)*ep(j,3))
0114         ept(i,j,6)=1.0d0/dsqrt(6.0d0)*(ep(i,2)*ep(j,2)+ep(i,3)*ep(j,3)-
0115      &            2.0d0*ep(i,4)*ep(j,4))
0116        enddo
0117       enddo
0118
0119       return
0120       end
0121 c********************************************************
0122       subroutine def_p_eq_p_dot_tp(n,k,e,j)
0123       implicit double precision(a-h, o-z)
0124         implicit integer(i-n)
0125       common/ppp/pp(4,40),guv(4)
0126       dimension e(4,4,2:6)
0127
0128       do i=1,4
0129         pp(i,n)=pp(1,k)*e(1,i,j)-pp(2,k)*e(2,i,j)-pp(3,k)*e(3,i,j)
0130      .         -pp(4,k)*e(4,i,j)
0131       enddo
0132
0133       return
0134       end
0135 c*************************************************

0136       double precision function sprpp(k1,k2)
0137       implicit double precision(a-h, o-z)
0138         implicit integer(i-n)
0139       common/ppp/pp(4,40),guv(4)
0140
0141       sprpp=pp(1,k1)*pp(1,k2)-pp(2,k1)*pp(2,k2)-pp(3,k1)*pp(3,k2)-
0142      .      pp(4,k1)*pp(4,k2)
0143
0144         return
0145         end
0146 c*************************************************

0147       double precision function epsp(n1,n2,n3,n4)
0148       implicit double precision(a-h, o-z)
0149         implicit integer(i-n)
0150       common/ppp/ pp(4,40),guv(4)
0151
0152       epsp=pp(1,n1)*(pp(2,n2)*(pp(3,n3)*pp(4,n4)-pp(4,n3)*pp(3,n4))
0153      .     -pp(3,n2)*(pp(2,n3)*pp(4,n4)-pp(4,n3)*pp(2,n4))
0154      .     +pp(4,n2)*(pp(2,n3)*pp(3,n4)-pp(3,n3)*pp(2,n4)))
0155      .     -pp(2,n1)*(pp(1,n2)*(pp(3,n3)*pp(4,n4)-pp(4,n3)*pp(3,n4))
0156      .     -pp(3,n2)*(pp(1,n3)*pp(4,n4)-pp(4,n3)*pp(1,n4))
0157      .     +pp(4,n2)*(pp(1,n3)*pp(3,n4)-pp(3,n3)*pp(1,n4)))
0158      .     -pp(3,n1)*(pp(2,n2)*(pp(1,n3)*pp(4,n4)-pp(4,n3)*pp(1,n4))
0159      .     -pp(1,n2)*(pp(2,n3)*pp(4,n4)-pp(4,n3)*pp(2,n4))
0160      .     +pp(4,n2)*(pp(2,n3)*pp(1,n4)-pp(1,n3)*pp(2,n4)))
0161      .     -pp(4,n1)*(pp(2,n2)*(pp(3,n3)*pp(1,n4)-pp(1,n3)*pp(3,n4))
0162      .     -pp(3,n2)*(pp(2,n3)*pp(1,n4)-pp(1,n3)*pp(2,n4))
0163      .     +pp(1,n2)*(pp(2,n3)*pp(3,n4)-pp(3,n3)*pp(2,n4)))
0164
0165       return
0166       end
0167 c********************************************************

0168       subroutine genppp()
0169       implicit double precision(a-h, o-z)
0170         implicit integer(i-n)
0171 #include "inclppp.h"
0172 #include "inclcon.h"
0173       common/counter/ibcstate,nev
0174
0175         if(ibcstate.eq.3) then
0176         p1p2=sprpp(2,1)
0177         p1p3=sprpp(3,1)
0178         p1p4=sprpp(4,1)
0179         p1p5=sprpp(5,1)
0180         p1p6=sprpp(6,1)
0181         p1p7=sprpp(7,1)
0182         p1p8=sprpp(8,1)
0183         p2p3=sprpp(3,2)
0184         p2p4=sprpp(4,2)
0185         p2p5=sprpp(5,2)
0186         p2p6=sprpp(6,2)
0187         p2p7=sprpp(7,2)
0188         p2p8=sprpp(8,2)
0189         p3p4=sprpp(4,3)
0190         p3p5=sprpp(5,3)
0191         p3p6=sprpp(6,3)
0192         p3p7=sprpp(7,3)
0193         p3p8=sprpp(8,3)
0194         p4p5=sprpp(5,4)
0195         p4p6=sprpp(6,4)
0196         p4p7=sprpp(7,4)
0197         p4p8=sprpp(8,4)
0198         p5p6=sprpp(6,5)
0199         p5p7=sprpp(7,5)
0200         p5p8=sprpp(8,5)
0201         p6p6=sprpp(6,6)
0202         p6p7=sprpp(7,6)
0203         p6p8=sprpp(8,6)
0204         p7p7=sprpp(7,7)
0205         p7p8=sprpp(8,7)
0206         p8p8=sprpp(8,8)
0207         end if
0208
0209         if(ibcstate.eq.4) then
0210         p1p2=sprpp(2,1)
0211         p1p3=sprpp(3,1)
0212         p1p4=sprpp(4,1)
0213         p1p5=sprpp(5,1)
0214         p1p6=sprpp(6,1)
0215         p1p7=sprpp(7,1)
0216         p2p3=sprpp(3,2)
0217         p2p4=sprpp(4,2)
0218         p2p5=sprpp(5,2)
0219         p2p6=sprpp(6,2)
0220         p2p7=sprpp(7,2)
0221         p3p4=sprpp(4,3)
0222         p3p5=sprpp(5,3)
0223         p3p6=sprpp(6,3)
0224         p3p7=sprpp(7,3)
0225         p4p5=sprpp(5,4)
0226         p4p6=sprpp(6,4)
0227         p4p7=sprpp(7,4)
0228         p5p6=sprpp(6,5)
0229         p5p7=sprpp(7,5)
0230         p6p6=sprpp(6,6)
0231         p6p7=sprpp(7,6)
0232         p7p7=sprpp(7,7)
0233         p3p5s2=p3p5**2
0234         p2p3s2=p2p3**2
0235         p3p4s2=p3p4**2
0236         p1p3s2=p1p3**2
0237         p3p6s2=p3p6**2
0238         end if
0239
0240         if(ibcstate.eq.5) then
0241         p1p2=sprpp(2,1)
0242         p1p3=sprpp(3,1)
0243         p1p4=sprpp(4,1)
0244         p1p5=sprpp(5,1)
0245         p1p6=sprpp(6,1)
0246         p1p7=sprpp(7,1)
0247         p1p8=sprpp(8,1)
0248         p2p3=sprpp(3,2)
0249         p2p4=sprpp(4,2)
0250         p2p5=sprpp(5,2)
0251         p2p6=sprpp(6,2)
0252         p2p7=sprpp(7,2)
0253         p2p8=sprpp(8,2)
0254         p3p4=sprpp(4,3)
0255         p3p5=sprpp(5,3)
0256         p3p6=sprpp(6,3)
0257         p3p7=sprpp(7,3)
0258         p3p8=sprpp(8,3)
0259         p4p5=sprpp(5,4)
0260         p4p6=sprpp(6,4)
0261         p4p7=sprpp(7,4)
0262         p4p8=sprpp(8,4)
0263         p5p6=sprpp(6,5)
0264         p5p7=sprpp(7,5)
0265         p5p8=sprpp(8,5)
0266         p6p6=sprpp(6,6)
0267         p6p7=sprpp(7,6)
0268         p6p8=sprpp(8,6)
0269         p7p7=sprpp(7,7)
0270         p7p8=sprpp(8,7)
0271         p8p8=sprpp(8,8)
0272         p3p5s2=p3p5**2
0273         p2p3s2=p2p3**2
0274         p1p3s2=p1p3**2
0275         p3p4s2=p3p4**2
0276         end if
0277
0278       if(ibcstate.eq.6) then
0279         p1p2=sprpp(2,1)
0280         p1p3=sprpp(3,1)
0281         p1p4=sprpp(4,1)
0282         p1p5=sprpp(5,1)
0283         p1p6=sprpp(6,1)
0284         p1p7=sprpp(7,1)
0285         p1p8=sprpp(8,1)
0286         p1p9=sprpp(9,1)
0287         p1p10=sprpp(10,1)
0288         p1p11=sprpp(11,1)
0289         p1p12=sprpp(12,1)
0290         p1p13=sprpp(13,1)
0291         p2p3=sprpp(3,2)
0292         p2p4=sprpp(4,2)
0293         p2p5=sprpp(5,2)
0294         p2p6=sprpp(6,2)
0295         p2p7=sprpp(7,2)
0296         p2p8=sprpp(8,2)
0297         p2p9=sprpp(9,2)
0298         p2p10=sprpp(10,2)
0299         p2p11=sprpp(11,2)
0300         p2p12=sprpp(12,2)
0301         p2p13=sprpp(13,2)
0302         p3p4=sprpp(4,3)
0303         p3p5=sprpp(5,3)
0304         p3p6=sprpp(6,3)
0305         p3p7=sprpp(7,3)
0306         p3p8=sprpp(8,3)
0307         p3p9=sprpp(9,3)
0308         p3p10=sprpp(10,3)
0309         p3p11=sprpp(11,3)
0310         p3p12=sprpp(12,3)
0311         p3p13=sprpp(13,3)
0312         p4p5=sprpp(5,4)
0313         p4p6=sprpp(6,4)
0314         p4p7=sprpp(7,4)
0315         p4p8=sprpp(8,4)
0316         p4p9=sprpp(9,4)
0317         p4p10=sprpp(10,4)
0318         p4p11=sprpp(11,4)
0319         p4p12=sprpp(12,4)
0320         p4p13=sprpp(13,4)
0321         p5p6=sprpp(6,5)
0322         p5p7=sprpp(7,5)
0323         p5p8=sprpp(8,5)
0324         p5p9=sprpp(9,5)
0325         p5p10=sprpp(10,5)
0326         p5p11=sprpp(11,5)
0327         p5p12=sprpp(12,5)
0328         p5p13=sprpp(13,5)
0329         p6p6=sprpp(6,6)
0330         p6p7=sprpp(7,6)
0331         p6p8=sprpp(8,6)
0332         p6p9=sprpp(9,6)
0333         p6p10=sprpp(10,6)
0334         p6p11=sprpp(11,6)
0335         p6p12=sprpp(12,6)
0336         p6p13=sprpp(13,6)
0337         p7p7=sprpp(7,7)
0338         p7p8=sprpp(8,7)
0339         p7p9=sprpp(9,7)
0340         p7p10=sprpp(10,7)
0341         p7p11=sprpp(11,7)
0342         p7p12=sprpp(12,7)
0343         p7p13=sprpp(13,7)
0344         p8p8=sprpp(8,8)
0345         p8p9=sprpp(9,8)
0346         p8p10=sprpp(10,8)
0347         p8p11=sprpp(11,8)
0348         p8p12=sprpp(12,8)
0349         p8p13=sprpp(13,8)
0350         p9p9=sprpp(9,9)
0351         p9p10=sprpp(10,9)
0352         p9p11=sprpp(11,9)
0353         p9p12=sprpp(12,9)
0354         p9p13=sprpp(13,9)
0355         p10p10=sprpp(10,10)
0356         p10p11=sprpp(11,10)
0357         p10p12=sprpp(12,10)
0358         p10p13=sprpp(13,10)
0359         p11p11=sprpp(11,11)
0360         p11p12=sprpp(12,11)
0361         p11p13=sprpp(13,11)
0362         p12p12=sprpp(12,12)
0363         p12p13=sprpp(13,12)
0364         p13p13=sprpp(13,13)
0365         end if
0366
0367         w20=1.0d0/(2*p1p2+2*p1p3*ffmcfmb-2*p1p3+2*p2p3*ffmcfmb-2*p2p3+
0368      . ffmcfmb**2*hbcm2-2*ffmcfmb*hbcm2-fmb2+hbcm2)
0369       w19=1.0d0/(2*p1p2)
0370       w18=1.0d0/(2*p3p4*ffmcfmb+2*p3p5*ffmcfmb+2*p4p5+ffmcfmb**2*
0371      . hbcm2+fmc2)
0372       w17=1.0d0/(2*p1p2-2*p1p3*ffmcfmb-2*p2p3*ffmcfmb+ffmcfmb**2*
0373      . hbcm2-fmc2)
0374       w16=1.0d0/(2*p2p3*ffmcfmb-2*p2p3-2*p2p5-2*p3p5*ffmcfmb+2*p3p5+
0375      . ffmcfmb**2*hbcm2-2*ffmcfmb*hbcm2+fmb2+hbcm2)
0376       w15=1.0d0/(2*p1p3*ffmcfmb-2*p1p3-2*p1p5-2*p3p5*ffmcfmb+2*p3p5+
0377      . ffmcfmb**2*hbcm2-2*ffmcfmb*hbcm2+fmb2+hbcm2)
0378       w14=1.0d0/(2*p3p4-fmb2+fmc2+hbcm2)
0379       w13=1.0d0/(2*p1p3*ffmcfmb-2*p1p3+ffmcfmb**2*hbcm2-2*ffmcfmb*
0380      . hbcm2-fmb2+hbcm2)
0381       w12=(-1.0d0)/(2*p1p3*ffmcfmb-ffmcfmb**2*hbcm2+fmc2)
0382       w11=(-1.0d0)/(2*p2p3*ffmcfmb+2*p2p4-2*p3p4*ffmcfmb-ffmcfmb**2*
0383      . hbcm2-fmc2)
0384       w10=(-1.0d0)/(2*p2p4)
0385       w9=(-1.0d0)/(2*p1p5)
0386       w8=1.0d0/(2*p2p3*ffmcfmb-2*p2p3+ffmcfmb**2*hbcm2-2*ffmcfmb*
0387      . hbcm2-fmb2+hbcm2)
0388       w7=(-1.0d0)/(2*p2p3*ffmcfmb-ffmcfmb**2*hbcm2+fmc2)
0389       w6=1.0d0/(2*p3p5+fmb2-fmc2+hbcm2)
0390       w5=(-1.0d0)/(2*p1p3*ffmcfmb+2*p1p4-2*p3p4*ffmcfmb-ffmcfmb**2*
0391      . hbcm2-fmc2)
0392       w4=(-1.0d0)/(2*p1p4)
0393       w3=(-1.0d0)/(2*p2p5)
0394       w2=(-1.0d0)/(2*p3p5*ffmcfmb-2*p3p5-ffmcfmb**2*hbcm2+2*ffmcfmb*
0395      . hbcm2-fmb2-hbcm2)
0396       w1=1.0d0/(2*p3p4*ffmcfmb+ffmcfmb**2*hbcm2+fmc2)
0397
0398       return
0399       end
0400 c******************************************************

0401       subroutine coupling()
0402       implicit double precision(a-h, o-z)
0403         implicit integer(i-n)
0404 #include "inclcon.h"
0405       common/cccc/c36,c35,c34,c33,c32,c31,c30,c29
0406      . ,c28,c27,c26,c25,c24,c23,c22,c21,c20,c19,c18,c17
0407      . ,c16,c15,c14,c13,c12,c11,c10,c9,c8,c7,c6,c5,c4,c3
0408      . ,c2,c1s2,c1s1
0409       common/counter/ibcstate,nev
0410
0411         c0 = 0
0412         if(ibcstate.eq.3) c0=cbc1p1
0413       if(ibcstate.eq.4) c0=cbcp0
0414       if(ibcstate.eq.5) c0=cbcp1
0415         if(ibcstate.eq.6) c0=cbcp2
0416
0417         c36=-c0
0418       c35=-c0
0419       c34=-c0
0420       c33=c0
0421       c32=c0
0422       c31=c0
0423       c30=c0
0424       c29=c0
0425       c28=c0
0426       c27=c0
0427       c26=c0
0428       c25=c0
0429       c24=c0
0430       c23=c0
0431       c22=c0
0432       c21=c0
0433       c20=c0
0434       c19=c0
0435       c18=c0
0436       c17=c0
0437       c16=c0
0438       c15=c0
0439       c14=c0
0440       c13=c0
0441       c12=c0
0442       c11=c0
0443       c10=c0
0444       c9=c0
0445       c8=c0
0446       c7=c0
0447       c6=c0
0448       c5=c0
0449       c4=c0
0450       c3=c0
0451       c2=c0
0452       c1s2=-c0
0453       c1s1=-c0
0454
0455       return
0456       end
0457 c************************

0458       subroutine genpppn(n)
0459       implicit double precision(a-h, o-z)
0460         implicit integer(i-n)
0461 #include "inclppp.h"
0462 #include "inclamp.h"
0463       common/ppp/pp(4,40),guv(4)
0464       common/counter/ibcstate,nev
0465
0466         if(ibcstate.eq.3) then
0467       do i=1,4
0468          pp(i,12)=0.d0
0469          pp(i,11)=0.d0
0470          pp(i,10)=0.d0
0471          pp(i,9)=0.d0
0472       enddo
0473         call wa0(3,9,c(2,n))
0474         call wa0(6,9,c(5,n))
0475         call wa0(7,9,c(8,n))
0476         call wa0(8,9,c(9,n))
0477         call wa0(2,9,c(31,n))
0478         call wa1(3,8,11,12,c(1,n))
0479         call wa1(3,2,11,12,c(28,n))
0480         call wa1(6,3,11,12,c(3,n))
0481         call wa1(6,8,11,12,c(4,n))
0482         call wa1(6,7,11,12,c(13,n))
0483         call wa1(6,2,11,12,c(18,n))
0484         call wa1(7,3,11,12,c(6,n))
0485         call wa1(7,8,11,12,c(7,n))
0486         call wa1(2,7,11,12,c(15,n))
0487         call wa1(2,8,11,12,c(27,n))
0488         call wa2(3,8,7,9,10,c(26,n))
0489         call wa2(3,2,7,9,10,c(16,n))
0490         call wa2(6,3,7,9,10,c(11,n))
0491         call wa2(6,3,8,9,10,c(25,n))
0492         call wa2(6,3,2,9,10,c(20,n))
0493         call wa2(6,8,7,9,10,c(14,n))
0494         call wa2(6,8,2,9,10,c(23,n))
0495         call wa2(6,2,7,9,10,c(17,n))
0496         call wa2(8,2,7,9,10,c(24,n))
0497         call wa2(2,3,8,9,10,c(32,n))
0498         call wa3(3,8,2,7,11,12,c(29,n))
0499         call wa3(6,3,8,7,11,12,c(12,n))
0500         call wa3(6,3,2,7,11,12,c(19,n))
0501         call wa3(6,3,2,8,11,12,c(30,n))
0502         call wa3(6,8,2,7,11,12,c(22,n))
0503         call wa4(6,3,8,2,7,9,10,c(21,n))
0504         epspn0=epsp(10,4,5,9)
0505         epspn1=epsp3(10,11,4)
0506         epspn2=epsp3(10,11,5)
0507         epspn3=epsp3(12,4,9)
0508         epspn4=epsp3(12,5,9)
0509         epspn5=epsp(11,12,4,5)
0510         epspn6=epsp3(11,12,4)
0511         epspn7=epsp3(12,4,5)
0512         p12p12=sprpp(12,12)
0513         p11p12=sprpp(11,12)
0514         p11p11=sprpp(11,11)
0515         p10p12=sprpp(10,12)
0516         p10p11=sprpp(10,11)
0517         p10p10=sprpp(10,10)
0518         p9p12=sprpp(9,12)
0519         p9p11=sprpp(9,11)
0520         p9p10=sprpp(9,10)
0521         p9p9=sprpp(9,9)
0522         p5p12=sprpp(5,12)
0523         p5p11=sprpp(5,11)
0524         p5p10=sprpp(5,10)
0525         p5p9=sprpp(5,9)
0526         p4p12=sprpp(4,12)
0527         p4p11=sprpp(4,11)
0528         p4p10=sprpp(4,10)
0529         p4p9=sprpp(4,9)
0530         p0p12=pp(1,12)
0531         p0p11=pp(1,11)
0532         p0p10=pp(1,10)
0533         p0p9=pp(1,9)
0534         p0p5=pp(1,5)
0535         p0p4=pp(1,4)
0536         end if
0537
0538         if(ibcstate.eq.5) then
0539         do i=1,4
0540          pp(i,12)=0.d0
0541          pp(i,11)=0.d0
0542          pp(i,10)=0.d0
0543          pp(i,9)=0.d0
0544         enddo
0545         call wa0(3,11,c(3,n))
0546         call wa0(6,11,c(7,n))
0547         call wa0(7,11,c(10,n))
0548         call wa0(8,11,c(11,n))
0549         call wa0(2,11,c(30,n))
0550         call wa1(3,8,9,10,c(2,n))
0551         call wa1(3,2,9,10,c(31,n))
0552         call wa1(6,3,9,10,c(5,n))
0553         call wa1(6,8,9,10,c(6,n))
0554         call wa1(6,7,9,10,c(15,n))
0555         call wa1(6,2,9,10,c(27,n))
0556         call wa1(7,3,9,10,c(8,n))
0557         call wa1(7,8,9,10,c(9,n))
0558         call wa1(2,7,9,10,c(28,n))
0559         call wa1(2,8,9,10,c(29,n))
0560         call wa2(3,8,7,11,12,c(1,n))
0561         call wa2(3,8,2,11,12,c(32,n))
0562         call wa2(3,2,7,11,12,c(17,n))
0563         call wa2(6,3,8,11,12,c(4,n))
0564         call wa2(6,3,7,11,12,c(13,n))
0565         call wa2(6,3,2,11,12,c(21,n))
0566         call wa2(6,8,7,11,12,c(16,n))
0567         call wa2(6,8,2,11,12,c(25,n))
0568         call wa2(6,2,7,11,12,c(19,n))
0569         call wa2(8,2,7,11,12,c(26,n))
0570         call wa3(3,8,2,7,9,10,c(18,n))
0571         call wa3(6,3,8,7,9,10,c(14,n))
0572         call wa3(6,3,8,2,9,10,c(23,n))
0573         call wa3(6,3,2,7,9,10,c(20,n))
0574         call wa3(6,8,2,7,9,10,c(24,n))
0575         call wa4(6,3,8,2,7,11,12,c(22,n))
0576         epspn0=epsp(10,4,5,9)
0577         epspn1=epsp3(10,4,9)
0578         epspn2=epsp3(10,11,4)
0579         epspn3=epsp3(10,11,5)
0580         epspn4=epsp3(10,4,5)
0581         epspn5=epsp3(12,4,9)
0582         epspn6=epsp3(12,5,9)
0583         epspn7=epsp(11,12,4,5)
0584         p12p12=sprpp(12,12)
0585         p11p12=sprpp(11,12)
0586         p11p11=sprpp(11,11)
0587         p10p12=sprpp(10,12)
0588         p10p11=sprpp(10,11)
0589         p10p10=sprpp(10,10)
0590         p9p12=sprpp(9,12)
0591         p9p11=sprpp(9,11)
0592         p9p10=sprpp(9,10)
0593         p9p9=sprpp(9,9)
0594         p5p12=sprpp(5,12)
0595         p5p11=sprpp(5,11)
0596         p5p10=sprpp(5,10)
0597         p5p9=sprpp(5,9)
0598         p4p12=sprpp(4,12)
0599         p4p11=sprpp(4,11)
0600         p4p10=sprpp(4,10)
0601         p4p9=sprpp(4,9)
0602         p0p12=pp(1,12)
0603         p0p11=pp(1,11)
0604         p0p10=pp(1,10)
0605         p0p9=pp(1,9)
0606         p0p5=pp(1,5)
0607         p0p4=pp(1,4)
0608         end if
0609
0610         if(ibcstate.eq.6) then
0611         do i=1,4
0612          pp(i,17)=0.d0
0613          pp(i,16)=0.d0
0614          pp(i,15)=0.d0
0615          pp(i,14)=0.d0
0616         enddo
0617         call wa0(3,14,c(33,n))
0618         call wa0(6,14,c(6,n))
0619         call wa0(7,14,c(10,n))
0620         call wa0(8,14,c(11,n))
0621         call wa0(9,14,c(12,n))
0622         call wa0(1,14,c(66,n))
0623         call wa0(10,14,c(29,n))
0624         call wa0(11,14,c(40,n))
0625         call wa0(12,14,c(62,n))
0626         call wa0(13,14,c(83,n))
0627         call wa1(3,8,16,17,c(1,n))
0628         call wa1(3,9,16,17,c(2,n))
0629         call wa1(3,11,16,17,c(41,n))
0630         call wa1(3,1,16,17,c(43,n))
0631         call wa1(3,12,16,17,c(48,n))
0632         call wa1(3,10,16,17,c(64,n))
0633         call wa1(3,13,16,17,c(84,n))
0634         call wa1(6,3,16,17,c(3,n))
0635         call wa1(6,8,16,17,c(4,n))
0636         call wa1(6,9,16,17,c(5,n))
0637         call wa1(6,1,16,17,c(38,n))
0638         call wa1(6,10,16,17,c(70,n))
0639         call wa1(6,11,16,17,c(51,n))
0640         call wa1(6,7,16,17,c(76,n))
0641         call wa1(6,12,16,17,c(92,n))
0642         call wa1(7,3,16,17,c(7,n))
0643         call wa1(7,8,16,17,c(8,n))
0644         call wa1(7,9,16,17,c(9,n))
0645         call wa1(7,1,16,17,c(77,n))
0646         call wa1(1,10,16,17,c(63,n))
0647         call wa1(1,8,16,17,c(65,n))
0648         call wa1(1,11,16,17,c(39,n))
0649         call wa1(1,9,16,17,c(42,n))
0650         call wa1(1,13,16,17,c(82,n))
0651         call wa1(10,7,16,17,c(28,n))
0652         call wa1(11,7,16,17,c(44,n))
0653         call wa1(12,1,16,17,c(93,n))
0654         call wa1(13,7,16,17,c(85,n))
0655         call wa2(3,8,1,14,15,c(74,n))
0656         call wa2(3,9,7,14,15,c(32,n))
0657         call wa2(3,11,1,14,15,c(73,n))
0658         call wa2(3,11,7,14,15,c(46,n))
0659         call wa2(3,10,1,14,15,c(81,n))
0660         call wa2(3,13,7,14,15,c(87,n))
0661         call wa2(3,6,12,14,15,c(94,n))
0662         call wa2(6,3,8,14,15,c(91,n))
0663         call wa2(6,3,9,14,15,c(80,n))
0664         call wa2(6,3,11,14,15,c(90,n))
0665         call wa2(6,3,10,14,15,c(88,n))
0666         call wa2(6,3,7,14,15,c(75,n))
0667         call wa2(6,8,7,14,15,c(15,n))
0668         call wa2(6,9,7,14,15,c(16,n))
0669         call wa2(6,1,10,14,15,c(20,n))
0670         call wa2(6,1,7,14,15,c(21,n))
0671         call wa2(6,1,8,14,15,c(26,n))
0672         call wa2(6,1,3,14,15,c(35,n))
0673         call wa2(6,1,9,14,15,c(37,n))
0674         call wa2(6,1,11,14,15,c(49,n))
0675         call wa2(6,1,12,14,15,c(61,n))
0676         call wa2(6,10,7,14,15,c(22,n))
0677         call wa2(6,11,7,14,15,c(59,n))
0678         call wa2(7,3,1,14,15,c(78,n))
0679         call wa2(8,3,7,14,15,c(31,n))
0680         call wa2(1,10,7,14,15,c(17,n))
0681         call wa2(1,8,7,14,15,c(23,n))
0682         call wa2(1,11,7,14,15,c(58,n))
0683         call wa2(1,9,7,14,15,c(60,n))
0684         call wa2(1,3,9,14,15,c(79,n))
0685         call wa2(10,3,7,14,15,c(27,n))
0686         call wa2(13,1,7,14,15,c(89,n))
0687         call wa3(3,8,1,7,16,17,c(68,n))
0688         call wa3(3,9,1,7,16,17,c(47,n))
0689         call wa3(3,11,1,7,16,17,c(45,n))
0690         call wa3(3,10,1,7,16,17,c(67,n))
0691         call wa3(3,13,1,7,16,17,c(86,n))
0692         call wa3(6,3,8,7,16,17,c(13,n))
0693         call wa3(6,3,9,7,16,17,c(14,n))
0694         call wa3(6,3,11,7,16,17,c(56,n))
0695         call wa3(6,3,10,7,16,17,c(72,n))
0696         call wa3(6,1,10,7,16,17,c(19,n))
0697         call wa3(6,1,8,7,16,17,c(25,n))
0698         call wa3(6,1,3,9,16,17,c(52,n))
0699         call wa3(6,1,3,11,16,17,c(50,n))
0700         call wa3(6,1,3,7,16,17,c(53,n))
0701         call wa3(6,1,3,12,16,17,c(57,n))
0702         call wa3(6,1,3,10,16,17,c(69,n))
0703         call wa3(6,1,3,8,16,17,c(71,n))
0704         call wa3(6,1,9,7,16,17,c(36,n))
0705         call wa3(6,1,11,7,16,17,c(54,n))
0706         call wa4(6,1,10,3,7,14,15,c(18,n))
0707         call wa4(6,1,8,3,7,14,15,c(24,n))
0708         call wa4(6,1,3,9,7,14,15,c(34,n))
0709         call wa4(6,1,3,11,7,14,15,c(55,n))
0710         epspn0=epsp(14,15,4,5)
0711         epspn1=epsp3(15,16,4)
0712         epspn2=epsp3(15,16,5)
0713         epspn3=epsp3(14,17,4)
0714         epspn4=epsp3(14,17,5)
0715         epspn5=epsp(16,17,4,5)
0716         epspn6=epsp3(16,17,4)
0717         epspn7=epsp3(17,4,5)
0718         p17p17=sprpp(17,17)
0719         p16p17=sprpp(16,17)
0720         p16p16=sprpp(16,16)
0721         p15p17=sprpp(15,17)
0722         p15p16=sprpp(15,16)
0723         p15p15=sprpp(15,15)
0724         p14p17=sprpp(14,17)
0725         p14p16=sprpp(14,16)
0726         p14p15=sprpp(14,15)
0727         p14p14=sprpp(14,14)
0728         p5p17=sprpp(5,17)
0729         p5p16=sprpp(5,16)
0730         p5p15=sprpp(5,15)
0731         p5p14=sprpp(5,14)
0732         p4p17=sprpp(4,17)
0733         p4p16=sprpp(4,16)
0734         p4p15=sprpp(4,15)
0735         p4p14=sprpp(4,14)
0736         p0p17=pp(1,17)
0737         p0p16=pp(1,16)
0738         p0p15=pp(1,15)
0739         p0p14=pp(1,14)
0740         p0p5=pp(1,5)
0741         p0p4=pp(1,4)
0742       end if
0743
0744         return
0745       end
0746 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0747 c%    the subroutine is used to do calcualtion for and plus the results     %

0748 c%  wa1 -- c* p1_hat * p0_hat                                               %

0749 c%  wa2 -- c* p2_hat * p1_hat * p0_hat                                      %

0750 c%  wa3 -- c* p3_hat * p2_hat * p1_hat * p0_hat                             %

0751 c%  wa4 -- c* p4_hat * p3_hat * p2_hat * p1_hat * p0_hat                    %

0752 c%  wa5 -- c* p5_hat * p4_hat * p3_hat * p2_hat * p1_hat * p0_hat           %

0753 c%  wa6 -- c* p6_hat * p5_hat * p4_hat * p3_hat * p2_hat * p1_hat * p0_hat  %

0754 c%         ........                                                         %

0755 c%  wan -- c* pn_hat * .......................................... * p0_hat  %

0756 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0757       subroutine wa0(n0,na,c)
0758       implicit double precision(a-h, o-z)
0759         implicit integer(i-n)
0760       common/ppp/pp(4,40),guv(4)
0761
0762       pp(1,na)=pp(1,na)+c*pp(1,n0)
0763       pp(2,na)=pp(2,na)+c*pp(2,n0)
0764       pp(3,na)=pp(3,na)+c*pp(3,n0)
0765       pp(4,na)=pp(4,na)+c*pp(4,n0)
0766
0767       return
0768       end
0769 c*********************

0770       subroutine wa1(n1,n0,na,nb,c)
0771       implicit double precision(a-h, o-z)
0772         implicit integer(i-n)
0773       common/ppp/pp(4,40),guv(4)
0774
0775       call pabp(pp(1,n1),pp(1,n0),pp(1,na),pp(1,nb),c)
0776
0777       return
0778       end
0779 c********************

0780       subroutine wa2(n2,n1,n0,na,nb,c)
0781       implicit double precision(a-h, o-z)
0782         implicit integer(i-n)
0783       common/ppp/pp(4,40),guv(4)
0784       dimension pa(4),pb(4)
0785
0786       call    pab(pp(1,n1),pp(1,n0),pa,pb)
0787       call paabbp(pp(1,n2),pa,pb,pp(1,na),pp(1,nb),c)
0788
0789       return
0790       end
0791 c************************

0792       subroutine wa3(n3,n2,n1,n0,na,nb,c)
0793       implicit double precision(a-h, o-z)
0794         implicit integer(i-n)
0795       common/ppp/pp(4,40),guv(4)
0796       dimension pa(4),pb(4),qa(4),qb(4)
0797
0798       call    pab(pp(1,n1),pp(1,n0),pa,pb)
0799       call  paabb(pp(1,n2),pa,pb,qa,qb)
0800       call paabbp(pp(1,n3),qa,qb,pp(1,na),pp(1,nb),c)
0801
0802       return
0803       end
0804 c*********************

0805       subroutine wa4(n4,n3,n2,n1,n0,na,nb,c)
0806       implicit double precision(a-h, o-z)
0807         implicit integer(i-n)
0808       common/ppp/pp(4,40),guv(4)
0809       dimension pa(4),pb(4),qa(4),qb(4)
0810
0811       call    pab(pp(1,n1),pp(1,n0),pa,pb)
0812       call  paabb(pp(1,n2),pa,pb,qa,qb)
0813       call  paabb(pp(1,n3),qa,qb,pa,pb)
0814       call paabbp(pp(1,n4),pa,pb,pp(1,na),pp(1,nb),c)
0815
0816       return
0817       end
0818 c********************

0819       subroutine wa5(n5,n4,n3,n2,n1,n0,na,nb,c)
0820       implicit double precision(a-h, o-z)
0821         implicit integer(i-n)
0822       common/ppp/pp(4,40),guv(4)
0823       dimension pa(4),pb(4),qa(4),qb(4)
0824
0825       call    pab(pp(1,n1),pp(1,n0),pa,pb)
0826       call  paabb(pp(1,n2),pa,pb,qa,qb)
0827       call  paabb(pp(1,n3),qa,qb,pa,pb)
0828       call  paabb(pp(1,n4),pa,pb,qa,qb)
0829       call paabbp(pp(1,n5),qa,qb,pp(1,na),pp(1,nb),c)
0830
0831       return
0832       end
0833 c***************

0834       subroutine wa6(n6,n5,n4,n3,n2,n1,n0,na,nb,c)
0835       implicit double precision(a-h, o-z)
0836         implicit integer(i-n)
0837       common/ppp/pp(4,40),guv(4)
0838       dimension pa(4),pb(4),qa(4),qb(4)
0839
0840       call    pab(pp(1,n1),pp(1,n0),pa,pb)
0841       call  paabb(pp(1,n2),pa,pb,qa,qb)
0842       call  paabb(pp(1,n3),qa,qb,pa,pb)
0843       call  paabb(pp(1,n4),pa,pb,qa,qb)
0844       call  paabb(pp(1,n5),qa,qb,pa,pb)
0845       call paabbp(pp(1,n6),pa,pb,pp(1,na),pp(1,nb),c)
0846
0847       return
0848       end
0849 c*****************

0850       double precision function epsp3(n1,n2,n3)
0851       implicit double precision(a-h, o-z)
0852         implicit integer(i-n)
0853       common/ppp/pp(4,40),guv(4)
0854
0855       epsp3=pp(2,n1)*(pp(3,n2)*pp(4,n3)-pp(4,n2)*pp(3,n3))
0856      .     +pp(3,n1)*(pp(4,n2)*pp(2,n3)-pp(2,n2)*pp(4,n3))
0857      .     +pp(4,n1)*(pp(2,n2)*pp(3,n3)-pp(3,n2)*pp(2,n3))
0858
0859       return
0860       end
0861 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0862 c%      this is written by j. x. wnag; june 20, 2004           %

0863 c%      there are   16 times,  9 plus   for subroutinue pab    %

0864 c%      there are   32 times, 24 plus   for subroutinue paabb  %

0865 c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0866       subroutine pab(p1,q,p1a,p1b)
0867       implicit double precision(a-h, o-z)
0868         implicit integer(i-n)
0869       dimension p1(4),q(4),p1a(4),p1b(4)
0870 c*******************************************************************

0871 c***     p1a(1)= p1.q

0872 c***     p1a_v=q(1)*p1_v - p1(1)*q_v

0873 c*******************************************************************

0874       p1a(1)=p1(1)*q(1)-p1(2)*q(2)-p1(3)*q(3)-p1(4)*q(4)
0875       p1a(2)=p1(2)*q(1)-p1(1)*q(2)
0876       p1a(3)=p1(3)*q(1)-p1(1)*q(3)
0877       p1a(4)=p1(4)*q(1)-p1(1)*q(4)
0878 c*******************************************************************

0879 c***     p1b(1)=0

0880 c***     p1b_v= p1_v x q_v

0881 c*******************************************************************

0882       p1b(1)=0.0d0
0883       p1b(2)=p1(3)*q(4)-p1(4)*q(3)
0884       p1b(3)=p1(4)*q(2)-p1(2)*q(4)
0885       p1b(4)=p1(2)*q(3)-p1(3)*q(2)
0886
0887       return
0888       end
0889 c******************************

0890       subroutine paabb(p2,pa,pb,p2a,p2b)
0891       implicit double precision(a-h, o-z)
0892         implicit integer(i-n)
0893       dimension p2(4),pa(4),pb(4),p2a(4),p2b(4)
0894 c*******************************************************************

0895 c***     p2a(1)= p2.pa

0896 c***     p2a_v=pa(1)*p2_v - p2(1)*pa_v + p2_v x pb_v

0897 c*******************************************************************

0898       p2a(1)=p2(1)*pa(1)-p2(2)*pa(2)-p2(3)*pa(3)-p2(4)*pa(4)
0899       p2a(2)=p2(2)*pa(1)-p2(1)*pa(2)+p2(3)*pb(4)-p2(4)*pb(3)
0900       p2a(3)=p2(3)*pa(1)-p2(1)*pa(3)+p2(4)*pb(2)-p2(2)*pb(4)
0901       p2a(4)=p2(4)*pa(1)-p2(1)*pa(4)+p2(2)*pb(3)-p2(3)*pb(2)
0902 c*******************************************************************

0903 c***     p2b(1)=-p2.pb

0904 c***     p2b_v= p2(1)*pb_v -pb(1)*p2_v + p2_v x pa_v

0905 c*******************************************************************

0906       p2b(1)=p2(2)*pb(2)+p2(3)*pb(3)+p2(4)*pb(4)-p2(1)*pb(1)
0907       p2b(2)=p2(1)*pb(2)-p2(2)*pb(1)+p2(3)*pa(4)-p2(4)*pa(3)
0908       p2b(3)=p2(1)*pb(3)-p2(3)*pb(1)+p2(4)*pa(2)-p2(2)*pa(4)
0909       p2b(4)=p2(1)*pb(4)-p2(4)*pb(1)+p2(2)*pa(3)-p2(3)*pa(2)
0910
0911       return
0912       end
0913 c******************************

0914       subroutine pabp(p1,q,p1a,p1b,c)
0915       implicit double precision(a-h, o-z)
0916         implicit integer(i-n)
0917       dimension p1(4),q(4),p1a(4),p1b(4)
0918 c*******************************************************************

0919 c***     p1a(1)= p1.q

0920 c***     p1a_v=q(1)*p1_v - p1(1)*q_v

0921 c*******************************************************************

0922       p1a(1)=p1a(1)+c*(p1(1)*q(1)-p1(2)*q(2)-p1(3)*q(3)-p1(4)*q(4))
0923       p1a(2)=p1a(2)+c*(p1(2)*q(1)-p1(1)*q(2))
0924       p1a(3)=p1a(3)+c*(p1(3)*q(1)-p1(1)*q(3))
0925       p1a(4)=p1a(4)+c*(p1(4)*q(1)-p1(1)*q(4))
0926 c*******************************************************************

0927 c***     p1b(1)=0

0928 c***     p1b_v= p1_v x q_v

0929 c*******************************************************************

0930 c      p1b(1)=0.0d0+pb1(1)

0931       p1b(2)=p1b(2)+c*(p1(3)*q(4)-p1(4)*q(3))
0932       p1b(3)=p1b(3)+c*(p1(4)*q(2)-p1(2)*q(4))
0933       p1b(4)=p1b(4)+c*(p1(2)*q(3)-p1(3)*q(2))
0934
0935       return
0936       end
0937 c**********************

0938       subroutine paabbp(p2,pa,pb,p2a,p2b,c)
0939       implicit double precision(a-h, o-z)
0940         implicit integer(i-n)
0941       dimension p2(4),pa(4),pb(4),p2a(4),p2b(4)
0942 c*******************************************************************

0943 c***     p2a(1)=p2a(1)+c*(p2.pa)

0944 c***     p2a_v=p2a_v+c*(pa(1)*p2_v - p2(1)*pa_v + p2_v x pb_v)

0945 c*******************************************************************

0946       p2a(1)=p2a(1)+c*(p2(1)*pa(1)-p2(2)*pa(2)-p2(3)*pa(3)-p2(4)*pa(4))
0947       p2a(2)=p2a(2)+c*(p2(2)*pa(1)-p2(1)*pa(2)+p2(3)*pb(4)-p2(4)*pb(3))
0948       p2a(3)=p2a(3)+c*(p2(3)*pa(1)-p2(1)*pa(3)+p2(4)*pb(2)-p2(2)*pb(4))
0949       p2a(4)=p2a(4)+c*(p2(4)*pa(1)-p2(1)*pa(4)+p2(2)*pb(3)-p2(3)*pb(2))
0950 c*******************************************************************

0951 c***     p2b(1)=p2b(1)-c*p2.pb

0952 c***     p2b_v= p2b_v+c*(p2(1)*pb_v -pb(1)*p2_v + p2_v x pa_v)

0953 c*******************************************************************

0954       p2b(1)=p2b(1)+c*(p2(2)*pb(2)+p2(3)*pb(3)+p2(4)*pb(4)-p2(1)*pb(1))
0955       p2b(2)=p2b(2)+c*(p2(1)*pb(2)-p2(2)*pb(1)+p2(3)*pa(4)-p2(4)*pa(3))
0956       p2b(3)=p2b(3)+c*(p2(1)*pb(3)-p2(3)*pb(1)+p2(4)*pa(2)-p2(2)*pa(4))
0957       p2b(4)=p2b(4)+c*(p2(1)*pb(4)-p2(4)*pb(1)+p2(2)*pa(3)-p2(3)*pa(2))
0958
0959       return
0960       end
0961