* *-----WTOXSSW------------------------------------------------------------ * real*8 function wtoxssw(ndim,x) implicit real*8(a-h,o-z) * common/wtkount/ik common/wtistrf/isf common/wtipt/ifz(51) common/wtfsw/tfactsw common/wtnpr/ipr,ipr0 common/wtcutsw/tcs,s0cut common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,g8,tfact,pih,alw, # eta,feta,beta,g2,tfacth common/wtapar/ars,s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,sth2,cth2, # hsth2,tsth2,scth2,asth2,tth2,rzm,rzm2,rzg,rzmg,szg, # szgs,opszgs,sth4,cth4,ve,vf,vfp,rbqm2,rszw,rszw2, # s0w,s0z * dimension x(ndim),ww(5) * if(ipr.eq.4) then ipro= 1 else if(ipr.eq.5) then ipro= 2 endif do ii=1,5 ww(ii)= 0.d0 enddo em2= em*em em4= em2*em2 em6= em4*em2 s2= s*s s3= s2*s rem= em/sqrt(s) rwms= rwm2 rwmq= rwms*rwms rwgs= rwg*rwg rwgm= rwg*rwm rwgms= rwgs*rwms rwmpgs= rwms+rwgs s0= s0cut * if(ndim.eq.4) then ymx= x(1) ypx= x(2) zx= x(3) xx= x(4) else if(ndim.eq.6) then uvx= x(1) vvx= x(2) ymx= x(3) ypx= x(4) zx= x(5) xx= x(6) endif * iz= 1 ik= ik+1 * if(ndim.eq.4) then ueps= 0.d0 uv= 1.d0 uvs= uv*uv ujc= 1.d0 veps= 0.d0 vv= 1.d0 vjc= 1.d0 else if(uvx.gt.1.d0) then iz= 0 ifz(1)= ifz(1)+1 go to 2 endif vvl= 0.25d0*(rem+sqrt(rem*rem+4.d0*s0))* # (rem+sqrt(rem*rem+4.d0*s0)) omul= 1.d0-vvl ueps= omul*(1.d0-uvx)**hbeti uv= 1.d0-ueps ujc= omul**hbet uvs= uv*uv * uvl= uv-vvl veps= uvl*(1.d0-vvx)**hbeti vv= uv-veps vjc0= 1.d0-vvl/uv if(vjc0.le.0.d0) then iz= 0 ifz(2)= ifz(2)+1 go to 2 else vjc= vjc0**hbet endif endif * if(vv.lt.0.d0) then iz= 0 ifz(2)= ifz(2)+1 go to 2 endif svv= sqrt(vv) xm= uv xp= vv/uv xmop= xm/xp sh= vv*s shs= sh*sh ssh= sqrt(sh) vvs= vv*vv vvi= 1.d0/vv vvis= vvi*vvi vvic= vvis*vvi vviq= vvic*vvi * opxp= 1.d0+xp opxm= 1.d0+xm omxp= veps/uv omxm= ueps if(isf.eq.0) then stfp= 1.d0 stfm= 1.d0 else if(isf.gt.0) then if(omxp.eq.0) then stfp= d0gl else rcpx= 0.25d0*opxp*opxp rcpy= xp iflp= 1 rclp= s21baf(rcpx,rcpy,iflp) stfp= d0gl+eob*omxp**omhb*(-0.5d0*opxp+ # feta*(-4.d0*opxp*log(omxp)+ # 3.d0*opxp*log(xp)+4.d0*rclp-5.d0-xp)) endif if(omxm.eq.0) then stfm= d0gl else rcmx= 0.25d0*opxm*opxm rcmy= xm iflm= 1 rclm= s21baf(rcmx,rcmy,iflm) stfm= d0gl+eob*omxm**omhb*(-0.5d0*opxm+ # feta*(-4.d0*opxm*log(omxm)+ # 3.d0*opxm*log(xm)+4.d0*rclm-5.d0-xm)) endif endif * stf= stfp*stfm * *-----ym is initialized * yml= (4.d0*s0+tcs)/(4.d0+tcs)/vv ymu= 1.d0-em/ssh if(yml.gt.ymu) then iz= 0 ifz(3)= ifz(3)+1 go to 2 endif * ymjc= ymu-yml ym= ymjc*ymx+yml yms= ym*ym ymvv= ym*vv ymvvs= ymvv*ymvv ymvvc= ymvvs*ymvv ymvvq= ymvvc*ymvv * omym= (1.d0-s0/vv)*(1.d0-ymx)+em/ssh*ymx omym= 1.d0-ym omyms= omym*omym omymc= omyms*omym * *-----yp is initialized * ypl= s0/ymvv+zv ypu= 1.d0 if(ypl.gt.ypu) then iz= 0 ifz(4)= ifz(4)+1 go to 2 endif ypjc= ypu-ypl yp= ypjc*ypx+ypl ypi= 1.d0/yp ypis= ypi*ypi ypic= ypis*ypi ypiq= ypic*ypi ypiv= ypiq*ypi ypivi= ypiv*ypi omyp= 1.d0-yp omyps= omyp*omyp omypc= omyps*omyp * qh0s= em2*(xm*ym+omxm)*ym/omym qhcs= qh0s+0.25d0*xm*xm*omym*tcs*s * z0= qh0s/ym/sh zc0= qhcs/ym/sh zc1= 1.d0-s0/ymvv zc2= yp-s0/ymvv zc= dmin1(zc0,zc1,zc2) if(zc.eq.zc0) then dzc= 0.25d0*xm*xm*omym*tcs*s/ym/sh else dzc= zc-z0 endif * if(z0.gt.zc) then iz= 0 ifz(5)= ifz(5)+1 go to 2 endif * *-----Loop for integration starts here * z is initialized * do ii=1,3 * if(ii.eq.1) then zjac= dzc/zc/z0/ym/sh zv= z0*zc/(dzc*zx+z0) else if(ii.eq.2) then zjac= log(zc/z0) zv= exp(zjac*zx+log(z0)) else if(ii.eq.3) then zjac= dzc zv= zjac*zx+z0 else if(ii.eq.4) then zjac= dzc/zc/z0/(ym*sh)**2 zv= z0*zc/(dzc*zx+z0) else if(ii.eq.5) then zjac= log(zc/z0)/ym/sh zv= exp(zjac*zx+log(z0)) endif * *-----x is initialized * if(ipro.eq.2) then ums= uqm*uqm/s dms= dqm*dqm/s rum= sqrt(ums) rdm= sqrt(dms) tha= (rum-rdm)*(rum-rdm) thn= (rum+rdm)*(rum+rdm) else if(ipro.eq.1) then dms= rmm*rmm/s endif * xl1= s0/ymvv+zv xl2= (1.d0-0.5d0/zv)*yp xl3= thn/ymvv+zv xl= dmax1(xl1,xl2,xl3) xu= yp * xu1= yp * xu2= 2.5d3/s/ymvv+zv * xu= dmin1(xu1,xu2) if(xl.gt.xu) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif * atma= ((xl-zv)*ymvv-rwms)/rwgm atmb= ((xu-zv)*ymvv-rwms)/rwgm * if(atma.gt.1.d0.and.atmb.gt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (-atmb+atma)/rwgm else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (-pi+atmb+atma)/rwgm else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (-pih+atmb+atma)/rwgm else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (pi-atmb-atma)/rwgm else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (atmb-atma)/rwgm else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm xjac0= (pih-atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm xjac0= (-pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= atan(atmb) zmbt= atmb/rwgm xjac0= (atmb-atma)/rwgm endif * xjac= xjac0/ymvv argx= rwgm*(xjac0*xx+zmat) iftn= 1 s07= s07aaf(argx,iftn) if(iftn.ne.0) then print*,' call to s07aaf, IFAIL = ',iftn endif xv= rwms/ymvv+zv+rwgm/ymvv*s07 xvs= xv*xv xvc= xvs*xv xvq= xvc*xv xvv= xvq*xv xvvi= xvv*xv * pwtm2= ((yp-zv)*ymvv-rwgm*s07)**2+rwms*rwgs * qms= zv*ymvv qmss= qms*qms qps= (yp-xv)*ymvv qpm= -0.5d0*yp*ymvv qpms= qpm*qpm qra= qps/qpm qras= qra*qra ophqra= 1.d0+0.5d0*qra opqra= 1.d0+qra tpqra= 2.d0+qra ophqras= ophqra*ophqra bqs0= qps+2.d0*qpm bqs= bqs0*(1.d0+qms/bqs0) qex= qms/qpm qexs= qex*qex qexx= qms*qex qt= qps*qms/qpms grm2= qpms*(1.d0-qt) grm= sqrt(grm2) opgrm= 1.d0+grm/qpm omgrm= 1.d0-grm/qpm ophqrae= ophqra+0.25d0*qras * if(ipro.eq.2) then umds= ums-dms dmus= dms-ums upds= ums+dms * if(ums.ne.0.d0.and.dms.ne.0.d0) then if(qt.lt.1.d-4) then denu= vvs*(bqs0+qms-umds)/(4.d0*ums*qpms- # bqs0*qps*qms-qps*qmss-2.d0*qps*qms* # upds+2.d0*qpm*qms*umds) dend= vvs*(bqs0+qms+umds)/(4.d0*dms*qpms- # bqs0*qps*qms-qps*qmss-2.d0*qps*qms* # upds-2.d0*qpm*qms*umds) denu= (1.d0+2.d0*bqs*ums/(bqs+umds)**2)*denu dend= (1.d0+2.d0*bqs*dms/(bqs-umds)**2)*dend * unum= qpm*tpqra+2.d0*dms-ums+ # qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra* # qexx*ophqrae-0.5d0*ums* # qex+0.5d0*dms*(1.d0-qra)*qex uden= -ums+0.5d0*qms*qra*ophqra+0.25d0*qra* # qexx*ophqrae-0.5d0*ums* # qex+0.5d0*dms*opqra*qex rlncu= log(unum/uden) dnum= qpm*tpqra+2.d0*ums-dms+ # qms*(1.d0-0.5d0*qra*ophqra)-0.25d0*qra* # qexx*ophqrae-0.5d0*dms* # qex+0.5d0*ums*(1.d0-qra)*qex dden= -dms+0.5d0*qms*qra*ophqra+0.25d0*qra* # qexx*ophqrae-0.5d0*dms* # qex+0.5d0*ums*opqra*qex rlncd= log(dnum/dden) else denu= vvs*(bqs-umds)/(4.d0*ums*qpms-qms*qpm*( # qra*tpqra*qpm+2.d0*opqra*dmus+qra*qms)) dend= vvs*(bqs+umds)/(4.d0*dms*qpms-qms*qpm*( # qra*tpqra*qpm+2.d0*opqra*umds+qra*qms)) unum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+ # dmus/qpm)+dms*opgrm-ums uden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+ # dmus/qpm)+dms*omgrm-ums dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+ # umds/qpm)+ums*opgrm-dms dden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+ # umds/qpm)+ums*omgrm-dms rlncu= -log(unum/uden) rlncd= -log(dnum/dden) endif else if(qt.lt.1.d-4) then denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) unum= qpm*tpqra+qms*(1.d0-0.5d0*qra*ophqra)- # 0.25d0*qexx*qra*ophqras uden= 0.5d0*qms*qra*(ophqra+0.5d0*qex*ophqras+ # qra/8.d0*qexs*(opqra*ophqra)) if((unum/uden).le.0.d0) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif rlncu= log(unum/uden) rlncd= rlncd else denu= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) dend= vvs*bqs/(-qms*qpm*(qra*tpqra*qpm+qra*qms)) unum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm uden= ophqra*(qpm-grm)+0.5d0*qms*omgrm dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*opgrm dden= ophqra*(qpm-grm)+0.5d0*qms*omgrm if((unum/uden).le.0.d0) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif rlncu= -log(unum/uden) rlncd= -log(dnum/dden) endif endif * else if(ipro.eq.1) then umds= -dms dmus= dms upds= dms * if(qt.lt.1.d-4) then dend= vvs*(bqs0+qms+umds)/(4.d0*dms*qpms- # bqs0*qps*qms-qps*qmss-2.d0*qps*qms* # upds-2.d0*qpm*qms*umds) dend= (1.d0+2.d0*bqs*dms/(bqs-umds)**2)*dend * dnum= qpm*tpqra-dms+qms*(1.d0-0.5d0*qra*ophqra)- # 0.25d0*qra*qexx*ophqrae-0.5d0*dms* # qex+0.5d0*ums*(1.d0-qra)*qex dden= -dms+0.5d0*qms*qra*ophqra+0.25d0*qra* # qexx*ophqrae-0.5d0*dms*qex+0.5d0*ums*opqra*qex rlncd= log(dnum/dden) else dend= vvs*(bqs+umds)/(4.d0*dms*qpms-qms*qpm*( # qra*tpqra*qpm+2.d0*opqra*umds+qra*qms)) dnum= ophqra*qpm+grm*ophqra+0.5d0*qms*(opgrm+ # umds/qpm)-dms dden= ophqra*(qpm-grm)+0.5d0*qms*(omgrm+ # umds/qpm)-dms rlncd= -log(dnum/dden) endif endif * bd= -xp*xm+(xp+xm)*em2/s qyjc= 1.d0-em4/s2*(1.d0+(xp+xm*omxm)/bd)+em6/s3/bd*(xp+xm) * if(ii.eq.1) then psf= em2*(2.d0*ymvvs+omxm*(2.d0*ym*vvs-ymvvs)) else if(ii.eq.2) then psf= -ymvvs+2.d0*vvs*(ym-1.d0) else if(ii.ge.3) then psf= ymvvs endif * tjac0= ujc*vjc*ymjc*ypjc*zjac*xjac/qyjc tjac= tjac0*stf*psf/sh rmu2= -bqs rlk= (rmu2-tha)*(rmu2-thn) tjac= tjac*sqrt(rlk)/rmu2 * if(ipro.eq.2.and.ii.le.2) then * wcomu1= ypic*rwms*(-8./9.*xvq) wcomu1= wcomu1+ypis*rwms*(16./9.*xvc) wcomu1= wcomu1+ypi*rwms*(-4./3.*xvs) wcomu1= wcomu1+rwms*(4./9.*xv) wcomu1= wcomu1+ymvv*ypic*(4./9.*xvv) wcomu1= wcomu1+ymvv*ypis*(-8./9.*xvq) wcomu1= wcomu1+ymvv*ypi*(2./3.*xvc) wcomu1= wcomu1+ymvv*(-2./9.*xvs) wcomu1= ypi*ymvv*wcomu1 * wcomu2= ypic*(4./9.*xvc) wcomu2= wcomu2+ypis*(-8./9.*xvs) wcomu2= wcomu2+ypi*(2./3.*xv) wcomu2= wcomu2-2./9. wcomu2= ypi*rwms*rwmpgs*wcomu2 * wcomu= wcomu1+wcomu2 * wcomd1= ypiq*rwms*(-2./9.*xvq) wcomd1= wcomd1+ypic*rwms*xvc*(4./9.+4./9.*xv) wcomd1= wcomd1+ypis*rwms*xvs*(-1./3.-8./9.*xv-2./9.*xvs) wcomd1= wcomd1+ypi*rwms*xv*(1./9.+2./3.*xv+4./9.*xvs) wcomd1= wcomd1+rwms*omyp*(-1./9.*xv) wcomd1= wcomd1+rwms*xv*(-1./9.-1./3.*xv) wcomd1= wcomd1+ymvv*ypiq*(1./9.*xvv) wcomd1= wcomd1+ymvv*ypic*xvq*(-2./9.-2./9.*xv) wcomd1= wcomd1+ymvv*ypis*xvc*(1./6.+4./9.*xv+1./9.*xvs) wcomd1= wcomd1+ymvv*ypi*xvs*(-1./18.-1./3.*xv-2./9.*xvs) wcomd1= wcomd1+ymvv*omyp*(1./18.*xvs) wcomd1= wcomd1+ymvv*xvs*(1./18.+1./6.*xv) wcomd1= ymvv*wcomd1 * wcomd2= ypiq*(1./9.*xvc) wcomd2= wcomd2+ypic*xvs*(-2./9.-2./9.*xv) wcomd2= wcomd2+ypis*xv*(1./6.+4./9.*xv+1./9.*xvs) wcomd2= wcomd2+ypi*(-1./18.-1./3.*xv-2./9.*xvs) wcomd2= wcomd2+omyp*(1./18.) wcomd2= wcomd2+(1./18.+1./6.*xv) wcomd2= rwms*rwmpgs*wcomd2 * wcomd= wcomd1+wcomd2 * wcom01= ypiq*rwms*(14./3.*xvq) wcom01= wcom01+ypic*rwms*(-28./3.*xvc-10./3.*xvq) wcom01= wcom01+ypis*rwms*xvs*(62./9.+20./3.*xv+11./9.*xvs) wcom01= wcom01+ypi*rwms*xv*(-20./9.-50./9.* # xv-22./9.*xvs-2./3.*xvc) wcom01= wcom01+rwms*omyp*xv*(11./18.+xv) wcom01= wcom01+rwms*omyps*(1./3.*xv) wcom01= wcom01+rwms*xv*(23./18.+3./2.*xv+4./3.*xvs) wcom01= wcom01+ymvv*ypiq*(-7./3.*xvv) wcom01= wcom01+ymvv*ypic*xvq*(5.+5./3.*xv) wcom01= wcom01+ymvv*ypis*xvc*(-71./18.-11./3.*xv-11./18.*xvs) wcom01= wcom01+ymvv*ypi*xvs*(29./18.+59./18.*xv # +14./9.*xvs+1./3.*xvc) wcom01= wcom01+ymvv*omyp*xv*(-1./6.-5./36.*xv-7./6.*xvs) wcom01= wcom01+ymvv*omyps*(-2./3.*xvs) wcom01= wcom01+ymvv*omypc*(-1./6.*xv) wcom01= wcom01+ymvv*xvs*(-41./36.-3./4.*xv-xvs) wcom01= ymvv*wcom01 * wcom02= ypiq*(-7./3.*xvc) wcom02= wcom02+ypic*xvs*(13./3.+5./3.*xv) wcom02= wcom02+ypis*xv*(-53./18.-3*xv-11./18.*xvs) wcom02= wcom02+ypi*(11./18.+41./18.*xv+8./9.*xvs+1./3.*xvc) wcom02= wcom02+omyp*(-5./36.-1./6.*xv) wcom02= wcom02+(-5./36.-3./4.*xv-1./3.*xvs) wcom02= wcom02*rwms*rwmpgs * wcom0= wcom01+wcom02 * denum= denu*ums/vv dendm= dend*dms/vv * wcomum= ypis/vv*rwms*rwmpgs*(4./9.*xvs) wcomum= wcomum+ypi/vv*rwms*rwmpgs*(-8./9.*xv) wcomum= wcomum+1./vv*rwms*rwmpgs*(4./9.) wcomum= wcomum+ymvv*ypis/vv*rwms*(-8./9.*xvc) wcomum= wcomum+ymvv*ypi/vv*rwms*(16./9.*xvs) wcomum= wcomum+ymvv/vv*rwms*(-8./9.*xv) wcomum= wcomum+ymvvs*ypis/vv*(4./9.*xvq) wcomum= wcomum+ymvvs*ypi/vv*(-8./9.*xvc) wcomum= wcomum+ymvvs/vv*(4./9.*xvs) * wcomdm= ypis/vv*rwms*rwmpgs*(1./9.*xvs) wcomdm= wcomdm+ypi/vv*rwms*rwmpgs*xv*(-2./9.-2./9.*xv) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyp*(2./9.*xv) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyps*(1./9.) wcomdm= wcomdm+1./vv*rwms*rwmpgs*xv*(2./9.+1./9.*xv) wcomdm= wcomdm+ymvv*ypis/vv*rwms*(-2./9.*xvc) wcomdm= wcomdm+ymvv*ypi/vv*rwms*xvs*(4./9.+4./9.*xv) wcomdm= wcomdm+ymvv/vv*rwms*omyp*(-4./9.*xvs) wcomdm= wcomdm+ymvv/vv*rwms*omyps*(-2./9.*xv) wcomdm= wcomdm+ymvv/vv*rwms*xvs*(-4./9.-2./9.*xv) wcomdm= wcomdm+ymvvs*ypis/vv*(1./9.*xvq) wcomdm= wcomdm+ymvvs*ypi/vv*xvc*(-2./9.-2./9.*xv) wcomdm= wcomdm+ymvvs/vv*omyp*(2./9.*xvc) wcomdm= wcomdm+ymvvs/vv*omyps*(1./9.*xvs) wcomdm= wcomdm+ymvvs/vv*xvc*(2./9.+1./9.*xv) * wcom= wcomu*rlncu+wcomd*rlncd+wcom0+ # ymvv*(wcomum*denum+wcomdm*dendm) * else if(ipro.eq.2.and.ii.ge.3) then * wcomu1= ypic*rwms*(-8./9.*xvq) wcomu1= wcomu1+ypis*rwms*(16./9.*xvc) wcomu1= wcomu1+ypi*rwms*(-4./3.*xvs) wcomu1= wcomu1+rwms*(4./9.*xv) wcomu1= wcomu1+ymvv*ypic*(4./9.*xvv) wcomu1= wcomu1+ymvv*ypis*(-8./9.*xvq) wcomu1= wcomu1+ymvv*ypi*(2./3.*xvc) wcomu1= wcomu1+ymvv*(-2./9.*xvs) wcomu1= ypi*ymvv*wcomu1 * wcomu2= ypic*(4./9.*xvc) wcomu2= wcomu2+ypis*(-8./9.*xvs) wcomu2= wcomu2+ypi*(2./3.*xv) wcomu2= wcomu2-2./9. wcomu2= ypi*rwms*rwmpgs*wcomu2 * wcomu= wcomu1+wcomu2 * wcomd1= ypiq*rwms*(-2./9.*xvq) wcomd1= wcomd1+ypic*rwms*xvc*(4./9.+4./9.*xv) wcomd1= wcomd1+ypis*rwms*xvs*(-1./3.-8./9.*xv-2./9.*xvs) wcomd1= wcomd1+ypi*rwms*xv*(1./9.+2./3.*xv+4./9.*xvs) wcomd1= wcomd1+rwms*omyp*(-1./9.*xv) wcomd1= wcomd1+rwms*xv*(-1./9.-1./3.*xv) wcomd1= wcomd1+ymvv*ypiq*(1./9.*xvv) wcomd1= wcomd1+ymvv*ypic*xvq*(-2./9.-2./9.*xv) wcomd1= wcomd1+ymvv*ypis*xvc*(1./6.+4./9.*xv+1./9.*xvs) wcomd1= wcomd1+ymvv*ypi*xvs*(-1./18.-1./3.*xv-2./9.*xvs) wcomd1= wcomd1+ymvv*omyp*(1./18.*xvs) wcomd1= wcomd1+ymvv*xvs*(1./18.+1./6.*xv) wcomd1= ymvv*wcomd1 * wcomd2= ypiq*(1./9.*xvc) wcomd2= wcomd2+ypic*xvs*(-2./9.-2./9.*xv) wcomd2= wcomd2+ypis*xv*(1./6.+4./9.*xv+1./9.*xvs) wcomd2= wcomd2+ypi*(-1./18.-1./3.*xv-2./9.*xvs) wcomd2= wcomd2+omyp*(1./18.) wcomd2= wcomd2+(1./18.+1./6.*xv) wcomd2= rwms*rwmpgs*wcomd2 * wcomd= wcomd1+wcomd2 * wcom01= ypiq*rwms*(14./3.*xvq) wcom01= wcom01+ypic*rwms*(-28./3.*xvc-10./3.*xvq) wcom01= wcom01+ypis*rwms*xvs*(62./9.+20./3.*xv+11./9.*xvs) wcom01= wcom01+ypi*rwms*xv*(-20./9.-50./9.* # xv-22./9.*xvs-2./3.*xvc) wcom01= wcom01+rwms*omyp*xv*(11./18.+xv) wcom01= wcom01+rwms*omyps*(1./3.*xv) wcom01= wcom01+rwms*xv*(23./18.+3./2.*xv+4./3.*xvs) wcom01= wcom01+ymvv*ypiq*(-7./3.*xvv) wcom01= wcom01+ymvv*ypic*xvq*(5.+5./3.*xv) wcom01= wcom01+ymvv*ypis*xvc*(-71./18.-11./3.*xv-11./18.*xvs) wcom01= wcom01+ymvv*ypi*xvs*(29./18.+59./18.*xv # +14./9.*xvs+1./3.*xvc) wcom01= wcom01+ymvv*omyp*xv*(-1./6.-5./36.*xv-7./6.*xvs) wcom01= wcom01+ymvv*omyps*(-2./3.*xvs) wcom01= wcom01+ymvv*omypc*(-1./6.*xv) wcom01= wcom01+ymvv*xvs*(-41./36.-3./4.*xv-xvs) wcom01= ymvv*wcom01 * wcom02= ypiq*(-7./3.*xvc) wcom02= wcom02+ypic*xvs*(13./3.+5./3.*xv) wcom02= wcom02+ypis*xv*(-53./18.-3*xv-11./18.*xvs) wcom02= wcom02+ypi*(11./18.+41./18.*xv+8./9.*xvs+1./3.*xvc) wcom02= wcom02+omyp*(-5./36.-1./6.*xv) wcom02= wcom02+(-5./36.-3./4.*xv-1./3.*xvs) wcom02= wcom02*rwms*rwmpgs * wcom0= wcom01+wcom02 * denum= denu*ums/vv dendm= dend*dms/vv * wcomum= ypis/vv*rwms*rwmpgs*(4./9.*xvs) wcomum= wcomum+ypi/vv*rwms*rwmpgs*(-8./9.*xv) wcomum= wcomum+1./vv*rwms*rwmpgs*(4./9.) wcomum= wcomum+ymvv*ypis/vv*rwms*(-8./9.*xvc) wcomum= wcomum+ymvv*ypi/vv*rwms*(16./9.*xvs) wcomum= wcomum+ymvv/vv*rwms*(-8./9.*xv) wcomum= wcomum+ymvvs*ypis/vv*(4./9.*xvq) wcomum= wcomum+ymvvs*ypi/vv*(-8./9.*xvc) wcomum= wcomum+ymvvs/vv*(4./9.*xvs) * wcomdm= ypis/vv*rwms*rwmpgs*(1./9.*xvs) wcomdm= wcomdm+ypi/vv*rwms*rwmpgs*xv*(-2./9.-2./9.*xv) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyp*(2./9.*xv) wcomdm= wcomdm+1./vv*rwms*rwmpgs*omyps*(1./9.) wcomdm= wcomdm+1./vv*rwms*rwmpgs*xv*(2./9.+1./9.*xv) wcomdm= wcomdm+ymvv*ypis/vv*rwms*(-2./9.*xvc) wcomdm= wcomdm+ymvv*ypi/vv*rwms*xvs*(4./9.+4./9.*xv) wcomdm= wcomdm+ymvv/vv*rwms*omyp*(-4./9.*xvs) wcomdm= wcomdm+ymvv/vv*rwms*omyps*(-2./9.*xv) wcomdm= wcomdm+ymvv/vv*rwms*xvs*(-4./9.-2./9.*xv) wcomdm= wcomdm+ymvvs*ypis/vv*(1./9.*xvq) wcomdm= wcomdm+ymvvs*ypi/vv*xvc*(-2./9.-2./9.*xv) wcomdm= wcomdm+ymvvs/vv*omyp*(2./9.*xvc) wcomdm= wcomdm+ymvvs/vv*omyps*(1./9.*xvs) wcomdm= wcomdm+ymvvs/vv*xvc*(2./9.+1./9.*xv) * wcom= wcomu*rlncu+wcomd*rlncd+wcom0+ # ymvv*(wcomum*denum+wcomdm*dendm) * dw1= ymvv*ypivi*rlncu*vvis*rwms*(56./3.*xvv) dw1= dw1+ymvv*ypivi*rlncu*rwms*(64./3.*xvv) dw1= dw1+ymvv*ypivi*rlncd*vvis*rwms*(16./3.*xvv) dw1= dw1+ymvv*ypivi*rlncd*rwms*(8*xvv) dw1= dw1+ymvv*ypivi*vvis*rwms*(-500./9.*xvv) dw1= dw1+ymvv*ypivi*rwms*(-232./3.*xvv) dw1= dw1+ymvv*ypiv*rlncu*vvis*rwms*(-48*xvq-136./9.*xvv) dw1= dw1+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv) dw1= dw1+ymvv*ypiv*rlncd*vvis*rwms*(-40./3.* # xvq-104./9.*xvv) dw1= dw1+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv) dw1= dw1+ymvv*ypiv*vvis*rwms*(400./3.*xvq+680./9.*xvv) dw1= dw1+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv) dw1= dw1+ymvv*ypiq*rlncu*vvis*rwms*(362./9.* # xvc+364./9.*xvq) dw1= dw1+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq) dw1= dw1+ymvv*ypiq*rlncd*vvis*rwms*(95./9.* # xvc+278./9.*xvq+8*xvv) dw1= dw1+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./ # 3.*xvq+28./3.*xvv) dw1= dw1+ymvv*ypiq*vvis*rwms*(-101*xvc- # 1754./9.*xvq-28*xvv) dw1= dw1+ymvv*ypiq*rwms*(-2411./9.*xvc-758./ # 3.*xvq-254./9.*xvv) dw1= dw1+ymvv*ypic*rlncu*vvis*rwms*(-100./9. # *xvs-328./9.*xvc) dw1= dw1+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc) dw1= dw1+ymvv*ypic*rlncd*vvis*rwms*(-22./9.* # xvs-250./9.*xvc-212./9.*xvq-16./9.*xvv) dw1= dw1+ymvv*ypic*rlncd*rwms*(-142./9.*xvs- # 496./9.*xvc-80./3.*xvq-4./3.*xvv) dw1= dw1+ymvv*ypic*vvis*rwms*(199./9.*xvs+ # 505./3.*xvc+242./3.*xvq+40./9.*xvv) dw1= dw1+ymvv*ypic*rwms*(1285./9.*xvs+2519./9. # *xvc+236./3.*xvq+10./3.*xvv) dw1= dw1+ymvv*ypis*rlncu*vvis*rwms*(2./9.*xv+12*xvs) dw1= dw1+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs) dw1= dw1+ymvv*ypis*rlncd*vvis*rwms*(-1./9.* # xv+82./9.*xvs+223./9.*xvc+6*xvq) dw1= dw1+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9. # *xvs+265./9.*xvc+10./3.*xvq) dw1= dw1+ymvv*ypis*vvis*rwms*(29./18.*xv-941. # /18.*xvs-497./6.*xvc-15*xvq-2./3.*xvv) dw1= dw1+ymvv*ypis*rwms*(-593./18.*xv-1270./9. # *xvs-169./2.*xvc-25./3.*xvq) dw1= dw1+ymvv*ypi*rlncu*vvis*rwms*(-8./9.*xv) dw1= dw1+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv) dw1= dw1+ymvv*ypi*rlncd*vvis*rwms*(-2./3.* # xv-98./9.*xvs-68./9.*xvc) dw1= dw1+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.* # xv-46./3.*xvs-8./3.*xvc) dw1= dw1+ymvv*ypi*vvis*rwms*(-1./3.+25./9. # *xv+319./9.*xvs+173./9.*xvc+8./3.*xvq) dw1= dw1+ymvv*ypi*rwms*(23./9.+277./9.*xv+383. # /9.*xvs+29./3.*xvc-2./3.*xvq) dw1= dw1+ymvv*rlncd*vvis*rwms*omyp*(8./9.*xv) dw1= dw1+ymvv*rlncd*vvis*rwms*(7./9.*xv+38./9.*xvs) dw1= dw1+ymvv*rlncd*rwms*omyp*(1./9.) dw1= dw1+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs) dw1= dw1+ymvv*vvis*rwms*omyp*(-11./9.*xv-8./3.*xvs) dw1= dw1+ymvv*vvis*rwms*omyps*(-2./3.*xv) dw1= dw1+ymvv*vvis*rwms*(1./3.-61./18.*xv-77./ # 9.*xvs-4*xvc) dw1= dw1+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs) dw1= dw1+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw1= dw1+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc) dw1= dw1+ymvvs*ypivi*rlncu*vvis*(-28./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncu*(-32./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncd*vvis*(-8./3.*xvvi) dw1= dw1+ymvvs*ypivi*rlncd*(-4*xvvi) dw1= dw1+ymvvs*ypivi*vvis*(250./9.*xvvi) dw1= dw1+ymvvs*ypivi*(116./3.*xvvi) dw1= dw1+ymvvs*ypiv*rlncu*vvis*(24.*xvv+68./9.*xvvi) dw1= dw1+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*vvis*(20./3.*xvv+52./9.*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi) dw1= dw1+ymvvs*ypiv*vvis*(-200./3.*xvv-340./9.*xvvi) dw1= dw1+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi) dw1= dw1+ymvvs*ypiq*rlncu*vvis*(-178./9.*xvq-188./9.*xvv) dw1= dw1+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv) dw1= dw1+ymvvs*ypiq*rlncd*vvis*(-46./9.* # xvq-142./9.*xvv-4*xvvi) dw1= dw1+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./ # 3.*xvv-14./3.*xvvi) dw1= dw1+ymvvs*ypiq*vvis*(146./3.*xvq+910./9.*xvv+14*xvvi) dw1= dw1+ymvvs*ypiq*(1495./9.*xvq+428./3.* # xvv+127./9.*xvvi) dw1= dw1+ymvvs*ypic*rlncu*vvis*(44./9.*xvc+176./9.*xvq) dw1= dw1+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq) dw1= dw1+ymvvs*ypic*rlncd*vvis*(8./9.*xvc+ # 131./9.*xvq+112./9.*xvv+8./9.*xvvi) dw1= dw1+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.* # xvq+46./3.*xvv+2./3.*xvvi) dw1= dw1+ymvvs*ypic*vvis*(-68./9.*xvc-91 # *xvq-130./3.*xvv-20./9.*xvvi) dw1= dw1+ymvvs*ypic*(-908./9.*xvc-1612./9.* # xvq-137./3.*xvv-5./3.*xvvi) dw1= dw1+ymvvs*ypis*rlncu*vvis*(2./9.*xvs-20./3.*xvc) dw1= dw1+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc) dw1= dw1+ymvvs*ypis*rlncd*vvis*(2./9.*xvs- # 44./9.*xvc-128./9.*xvq-10./3.*xvv) dw1= dw1+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3. # *xvc-58./3.*xvq-2*xvv) dw1= dw1+ymvvs*ypis*vvis*(-23./9.*xvs+ # 523./18.*xvc+49*xvq+25./3.*xvv+1./3.*xvvi) dw1= dw1+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv) dw1= dw1+ymvvs*ypi*rlncu*vvis*(4./9.*xvs) dw1= dw1+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs) dw1= dw1+ymvvs*ypi*rlncd*vvis*(1./3.*xvs+ # 64./9.*xvc+43./9.*xvq) dw1= dw1+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs # +101./9.*xvc+2*xvq) dw1= dw1+ymvvs*ypi*vvis*(1./3.*xv-19./18.* # xvs-215./9.*xvc-221./18.*xvq-5./3.*xvv) dw1= dw1+ymvvs*ypi*(-32./9.*xv-163./6.*xvs # -203./6.*xvc-20./3.*xvq+1./3.*xvv) dw1= dw1+ymvvs*rlncd*vvis*omyp*(-7./9.*xvs) dw1= dw1+ymvvs*rlncd*vvis*(-5./9.*xvs-28./9.*xvc) dw1= dw1+ymvvs*rlncd*omyp*(-1./9.*xv) dw1= dw1+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc) dw1= dw1+ymvvs*vvis*omyp*(-1./24.+5./12.*xv- # 7./18.*xvs+10./3.*xvc) dw1= dw1+ymvvs*vvis*omyps*(-2./3.*xv+5./3.*xvs) dw1= dw1+ymvvs*vvis*omypc*(1./3.*xv) dw1= dw1+ymvvs*vvis*(-5./12.*xv+109./36.*xvs # +49./9.*xvc+10./3.*xvq) dw1= dw1+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc) dw1= dw1+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw1= dw1+ymvvs*omypc*(1./6.-1./6.*xv) dw1= dw1+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq) dw1= dw1+ypivi*rlncu*vvis*rwms*rwmpgs*(-28./3.*xvq) dw1= dw1+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq) dw1= dw1+ypivi*rlncd*vvis*rwms*rwmpgs*(-8./3.*xvq) dw1= dw1+ypivi*rlncd*rwms*rwmpgs*(-4*xvq) dw1= dw1+ypivi*vvis*rwms*rwmpgs*(250./9.*xvq) dw1= dw1+ypivi*rwms*rwmpgs*(116./3.*xvq) dw1= dw1+ypiv*rlncu*vvis*rwms*rwmpgs*(24.*xvc+68./9.*xvq) dw1= dw1+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq) dw1= dw1+ypiv*rlncd*vvis*rwms*rwmpgs*(20./3.*xvc+52./9.*xvq) dw1= dw1+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq) dw1= dw1+ypiv*vvis*rwms*rwmpgs*(-200./3.* # xvc-340./9.*xvq) dw1= dw1+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq) dw1= dw1+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./ # 9.*xvs-176./9.*xvc) dw1= dw1+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc) dw1= dw1+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./9. # *xvs-136./9.*xvc-4*xvq) dw1= dw1+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs # -62./3.*xvc-14./3.*xvq) dw1= dw1+ypiq*vvis*rwms*rwmpgs*(157./3.*xvs # +844./9.*xvc+14*xvq) dw1= dw1+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq) dw1= dw1+ypic*rlncu*vvis*rwms*rwmpgs*(56./9.* # xv+152./9.*xvs) dw1= dw1+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs) dw1= dw1+ypic*rlncd*vvis*rwms*rwmpgs*(14./9.* # xv+119./9.*xvs+100./9.*xvc+8./9.*xvq) dw1= dw1+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.* # xvs+34./3.*xvc+2./3.*xvq) dw1= dw1+ypic*vvis*rwms*rwmpgs*(-131./9.*xv # -232./3.*xvs-112./3.*xvc-20./9.*xvq) dw1= dw1+ypic*rwms*rwmpgs*(-383./9.*xv-907./9. # *xvs-33*xvc-5./3.*xvq) dw1= dw1+ypis*rlncu*vvis*rwms*rwmpgs*(-4./9.-16./3.*xv) dw1= dw1+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv) dw1= dw1+ypis*rlncd*vvis*rwms*rwmpgs*(-1./9. # -38./9.*xv-95./9.*xvs-8./3.*xvc) dw1= dw1+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9. # *xv-91./9.*xvs-4./3.*xvc) dw1= dw1+ypis*vvis*rwms*rwmpgs*(17./18.+209./ # 9.*xv+203./6.*xvs+20./3.*xvc+1./3.*xvq) dw1= dw1+ypis*rwms*rwmpgs*(107./18.+341./9.*xv # +499./18.*xvs+10./3.*xvc) dw1= dw1+ypi*rlncu*vvis*rwms*rwmpgs*(4./9.) dw1= dw1+ypi*rlncu*rwms*rwmpgs*(4./9.) dw1= dw1+ypi*rlncd*vvis*rwms*rwmpgs*(1./3.+ # 34./9.*xv+25./9.*xvs) dw1= dw1+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs) dw1= dw1+ypi*vvis*rwms*rwmpgs*(-31./18.- # 104./9.*xv-125./18.*xvs-xvc) dw1= dw1+ypi*rwms*rwmpgs*(-83./18.-169./18.* # xv-3*xvs+1./3.*xvc) dw1= dw1+rlncd*vvis*rwms*rwmpgs*omyp*(-1./9.) dw1= dw1+rlncd*vvis*rwms*rwmpgs*(-2./9.-10./9.*xv) dw1= dw1+rlncd*rwms*rwmpgs*(-2./3.) dw1= dw1+vvis*rwms*rwmpgs*omyp*(5./18.+1./3.*xv) dw1= dw1+vvis*rwms*rwmpgs*(7./9.+22./9.*xv+xvs) dw1= dw1+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw1= dw1+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs) dw1= dw1+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./9.*xvq) dw1= dw1+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4./3.*xvc) dw1= dw1+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4./3.*xvs) dw1= dw1+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs) dw1= dw1+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./9.*xv) dw1= dw1+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv) dw1= dw1+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.) dw1= dw1+denu*ymvvc*ypiq*vviq*rwms*(8./9.*xvv) dw1= dw1+denu*ymvvc*ypic*vviq*rwms*(-8./3.*xvq) dw1= dw1+denu*ymvvc*ypis*vviq*rwms*(8./3.*xvc) dw1= dw1+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc) dw1= dw1+denu*ymvvc*ypi*vviq*rwms*(-8./9.*xvs) dw1= dw1+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs) dw1= dw1+denu*ymvvc*vvis*rwms*(4./9.*xv) dw1= dw1+denu*ymvvq*ypiq*vviq*(-4./9.*xvvi) dw1= dw1+denu*ymvvq*ypic*vviq*(4./3.*xvv) dw1= dw1+denu*ymvvq*ypis*vviq*(-4./3.*xvq) dw1= dw1+denu*ymvvq*ypis*vvis*(-2./9.*xvq) dw1= dw1+denu*ymvvq*ypi*vviq*(4./9.*xvc) dw1= dw1+denu*ymvvq*ypi*vvis*(4./9.*xvc) dw1= dw1+denu*ymvvq*vvis*(-2./9.*xvs) dw1= dw1+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./9.*xvq) dw1= dw1+dend*ymvvs*ypic*vviq*rwms*rwmpgs # *(1./3.*xvc+2./9.*xvq) dw1= dw1+dend*ymvvs*ypis*vviq*rwms*rwmpgs # *(-1./3.*xvs-2./3.*xvc-1./9.*xvq) dw1= dw1+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs) dw1= dw1+dend*ymvvs*ypi*vviq*rwms*rwmpgs # *(1./9.*xv+2./3.*xvs+1./3.*xvc) dw1= dw1+dend*ymvvs*ypi*vvis*rwms*rwmpgs # *(1./9.*xv+1./9.*xvs) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./9.*xv) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*(-1./9. # *xv-1./3.*xvs) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs) dw1= dw1+dend*ymvvc*ypiq*vviq*rwms*(2./9.*xvv) dw1= dw1+dend*ymvvc*ypic*vviq*rwms*(-2./3.*xvq-4./9.*xvv) dw1= dw1+dend*ymvvc*ypis*vviq*rwms*(2./3.* # xvc+4./3.*xvq+2./9.*xvv) dw1= dw1+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc) dw1= dw1+dend*ymvvc*ypi*vviq*rwms*(-2./ # 9.*xvs-4./3.*xvc-2./3.*xvq) dw1= dw1+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc) dw1= dw1+dend*ymvvc*vviq*rwms*omyp*(2./9.*xvs) dw1= dw1+dend*ymvvc*vviq*rwms*(2./9.*xvs+2./3.*xvc) dw1= dw1+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs) dw1= dw1+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv) dw1= dw1+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc) dw1= dw1+dend*ymvvq*ypiq*vviq*(-1./9.*xvvi) dw1= dw1+dend*ymvvq*ypic*vviq*(1./3.*xvv+2./9.*xvvi) dw1= dw1+dend*ymvvq*ypis*vviq*(-1./3.* # xvq-2./3.*xvv-1./9.*xvvi) dw1= dw1+dend*ymvvq*ypis*vvis*(-1./18.*xvq) dw1= dw1+dend*ymvvq*ypi*vviq*(1./9.*xvc # +2./3.*xvq+1./3.*xvv) dw1= dw1+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq) dw1= dw1+dend*ymvvq*vviq*omyp*(-1./9.*xvc) dw1= dw1+dend*ymvvq*vviq*(-1./9.*xvc-1./3.*xvq) dw1= dw1+dend*ymvvq*vvis*omyp*(-1./9.*xvc) dw1= dw1+dend*ymvvq*vvis*omyps*(-1./18.*xvs) dw1= dw1+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq) * dw2= ymvv*ypivi*rlncu*vvis*rwms*(56*xvv) dw2= dw2+ymvv*ypivi*rlncu*rwms*(64./3.*xvv) dw2= dw2+ymvv*ypivi*rlncd*vvis*rwms*(16.*xvv) dw2= dw2+ymvv*ypivi*rlncd*rwms*(8*xvv) dw2= dw2+ymvv*ypivi*vvis*rwms*(-500./3.*xvv) dw2= dw2+ymvv*ypivi*rwms*(-232./3.*xvv) dw2= dw2+ymvv*ypiv*rlncu*vvis*rwms*(-144*xvq-136./3.*xvv) dw2= dw2+ymvv*ypiv*rlncu*rwms*(-64*xvq-40./3.*xvv) dw2= dw2+ymvv*ypiv*rlncd*vvis*rwms*(-40*xvq-104./3.*xvv) dw2= dw2+ymvv*ypiv*rlncd*rwms*(-24.*xvq-16.*xvv) dw2= dw2+ymvv*ypiv*vvis*rwms*(400.*xvq+680./3.*xvv) dw2= dw2+ymvv*ypiv*rwms*(232.*xvq+260./3.*xvv) dw2= dw2+ymvv*ypiq*rlncu*vvis*rwms*(362./3.* # xvc+364./3.*xvq) dw2= dw2+ymvv*ypiq*rlncu*rwms*(226./3.*xvc+116./3.*xvq) dw2= dw2+ymvv*ypiq*rlncd*vvis*rwms*(95./3.* # xvc+278./3.*xvq+24.*xvv) dw2= dw2+ymvv*ypiq*rlncd*rwms*(85./3.*xvc+142./ # 3.*xvq+28./3.*xvv) dw2= dw2+ymvv*ypiq*vvis*rwms*(-303*xvc- # 1754./3.*xvq-84*xvv) dw2= dw2+ymvv*ypiq*rwms*(-2411./9.*xvc-758./ # 3.*xvq-254./9.*xvv) dw2= dw2+ymvv*ypic*rlncu*vvis*rwms*(-100./3. # *xvs-328./3.*xvc) dw2= dw2+ymvv*ypic*rlncu*rwms*(-376./9.*xvs-380./9.*xvc) dw2= dw2+ymvv*ypic*rlncd*vvis*rwms*(-22./3.* # xvs-250./3.*xvc-212./3.*xvq-16./3.*xvv) dw2= dw2+ymvv*ypic*rlncd*rwms*(-142./9.*xvs- # 496./9.*xvc-80./3.*xvq-4./3.*xvv) dw2= dw2+ymvv*ypic*vvis*rwms*(199./3.*xvs+ # 505*xvc+242.*xvq+40./3.*xvv) dw2= dw2+ymvv*ypic*rwms*(1285./9.*xvs+2519./9. # *xvc+236./3.*xvq+10./3.*xvv) dw2= dw2+ymvv*ypis*rlncu*vvis*rwms*(2./3.*xv+36*xvs) dw2= dw2+ymvv*ypis*rlncu*rwms*(86./9.*xv+184./9.*xvs) dw2= dw2+ymvv*ypis*rlncd*vvis*rwms*(-1./3.* # xv+82./3.*xvs+223./3.*xvc+18*xvq) dw2= dw2+ymvv*ypis*rlncd*rwms*(32./9.*xv+274./9. # *xvs+265./9.*xvc+10./3.*xvq) dw2= dw2+ymvv*ypis*vvis*rwms*(29./6.*xv-941./ # 6.*xvs-497./2.*xvc-45*xvq-2*xvv) dw2= dw2+ymvv*ypis*rwms*(-593./18.*xv-1270./9. # *xvs-169./2.*xvc-25./3.*xvq) dw2= dw2+ymvv*ypi*rlncu*vvis*rwms*(-8./3.*xv) dw2= dw2+ymvv*ypi*rlncu*rwms*(-4./9.-32./9.*xv) dw2= dw2+ymvv*ypi*rlncd*vvis*rwms*(-2*xv- # 98./3.*xvs-68./3.*xvc) dw2= dw2+ymvv*ypi*rlncd*rwms*(-1./9.-62./9.* # xv-46./3.*xvs-8./3.*xvc) dw2= dw2+ymvv*ypi*vvis*rwms*(-1+25./3.*xv # +319./3.*xvs+173./3.*xvc+8*xvq) dw2= dw2+ymvv*ypi*rwms*(23./9.+277./9.*xv+383. # /9.*xvs+29./3.*xvc-2./3.*xvq) dw2= dw2+ymvv*rlncd*vvis*rwms*omyp*(8./3.*xv) dw2= dw2+ymvv*rlncd*vvis*rwms*(7./3.*xv+38./3.*xvs) dw2= dw2+ymvv*rlncd*rwms*omyp*(1./9.) dw2= dw2+ymvv*rlncd*rwms*(1./9.+10./3.*xv+2./3.*xvs) dw2= dw2+ymvv*vvis*rwms*omyp*(-11./3.*xv-8*xvs) dw2= dw2+ymvv*vvis*rwms*omyps*(-2*xv) dw2= dw2+ymvv*vvis*rwms*(1.-61./6.*xv-77./3.*xvs-12*xvc) dw2= dw2+ymvv*rwms*omyp*(-11./18.-7./3.*xv+2./3.*xvs) dw2= dw2+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw2= dw2+ymvv*rwms*(-13./9.-73./9.*xv-20./3.*xvs+xvc) dw2= dw2+ymvvs*ypivi*rlncu*vvis*(-28*xvvi) dw2= dw2+ymvvs*ypivi*rlncu*(-32./3.*xvvi) dw2= dw2+ymvvs*ypivi*rlncd*vvis*(-8*xvvi) dw2= dw2+ymvvs*ypivi*rlncd*(-4*xvvi) dw2= dw2+ymvvs*ypivi*vvis*(250./3.*xvvi) dw2= dw2+ymvvs*ypivi*(116./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncu*vvis*(72*xvv+68./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncu*(104./3.*xvv+20./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*vvis*(20*xvv+52./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*(40./3.*xvv+8*xvvi) dw2= dw2+ymvvs*ypiv*vvis*(-200*xvv-340./3.*xvvi) dw2= dw2+ymvvs*ypiv*(-386./3.*xvv-130./3.*xvvi) dw2= dw2+ymvvs*ypiq*rlncu*vvis*(-178./3.*xvq-188./3.*xvv) dw2= dw2+ymvvs*ypiq*rlncu*(-400./9.*xvq-64./3.*xvv) dw2= dw2+ymvvs*ypiq*rlncd*vvis*(-46./3.* # xvq-142./3.*xvv-12*xvvi) dw2= dw2+ymvvs*ypiq*rlncd*(-157./9.*xvq-80./ # 3.*xvv-14./3.*xvvi) dw2= dw2+ymvvs*ypiq*vvis*(146*xvq+910./3.*xvv+42.*xvvi) dw2= dw2+ymvvs*ypiq*(1495./9.*xvq+428./3.*xvv+127./9.*xvvi) dw2= dw2+ymvvs*ypic*rlncu*vvis*(44./3.*xvc+176./3.*xvq) dw2= dw2+ymvvs*ypic*rlncu*(244./9.*xvc+232./9.*xvq) dw2= dw2+ymvvs*ypic*rlncd*vvis*(8./3.*xvc+ # 131./3.*xvq+112./3.*xvv+8./3.*xvvi) dw2= dw2+ymvvs*ypic*rlncd*(97./9.*xvc+313./9.* # xvq+46./3.*xvv+2./3.*xvvi) dw2= dw2+ymvvs*ypic*vvis*(-68./3.*xvc- # 273*xvq-130*xvv-20./3.*xvvi) dw2= dw2+ymvvs*ypic*(-908./9.*xvc-1612./9.* # xvq-137./3.*xvv-5./3.*xvvi) dw2= dw2+ymvvs*ypis*rlncu*vvis*(2./3.*xvs-20*xvc) dw2= dw2+ymvvs*ypis*rlncu*(-64./9.*xvs-128./9.*xvc) dw2= dw2+ymvvs*ypis*rlncd*vvis*(2./3.*xvs- # 44./3.*xvc-128./3.*xvq-10*xvv) dw2= dw2+ymvvs*ypis*rlncd*(-25./9.*xvs-64./3. # *xvc-58./3.*xvq-2*xvv) dw2= dw2+ymvvs*ypis*vvis*(-23./3.*xvs+ # 523./6.*xvc+147*xvq+25*xvv+xvvi) dw2= dw2+ymvvs*ypis*(28*xvs+935./9.*xvc+511./9.*xvq+5*xvv) dw2= dw2+ymvvs*ypi*rlncu*vvis*(4./3.*xvs) dw2= dw2+ymvvs*ypi*rlncu*(4./9.*xv+28./9.*xvs) dw2= dw2+ymvvs*ypi*rlncd*vvis*(xvs+64./3.*xvc+43./3.*xvq) dw2= dw2+ymvvs*ypi*rlncd*(1./9.*xv+49./9.*xvs # +101./9.*xvc+2*xvq) dw2= dw2+ymvvs*ypi*vvis*(xv-19./6.*xvs- # 215./3.*xvc-221./6.*xvq-5*xvv) dw2= dw2+ymvvs*ypi*(-32./9.*xv-163./6.*xvs # -203./6.*xvc-20./3.*xvq+1./3.*xvv) dw2= dw2+ymvvs*rlncd*vvis*omyp*(-7./3.*xvs) dw2= dw2+ymvvs*rlncd*vvis*(-5./3.*xvs-28./3.*xvc) dw2= dw2+ymvvs*rlncd*omyp*(-1./9.*xv) dw2= dw2+ymvvs*rlncd*(-1./9.*xv-8./3.*xvs-2./3.*xvc) dw2= dw2+ymvvs*vvis*omyp*(-1./8.+5./4.*xv-7./6.*xvs+10*xvc) dw2= dw2+ymvvs*vvis*omyps*(-2*xv+5*xvs) dw2= dw2+ymvvs*vvis*omypc*(xv) dw2= dw2+ymvvs*vvis*(-5./4.*xv+109./12.*xvs # +49./3.*xvc+10*xvq) dw2= dw2+ymvvs*omyp*(5./24.+1./36.*xv+25./6.*xvs-5./6.*xvc) dw2= dw2+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw2= dw2+ymvvs*omypc*(1./6.-1./6.*xv) dw2= dw2+ymvvs*(85./36.*xv+251./36.*xvs+13./2.*xvc-5./6.*xvq) dw2= dw2+ypivi*rlncu*vvis*rwms*rwmpgs*(-28*xvq) dw2= dw2+ypivi*rlncu*rwms*rwmpgs*(-32./3.*xvq) dw2= dw2+ypivi*rlncd*vvis*rwms*rwmpgs*(-8*xvq) dw2= dw2+ypivi*rlncd*rwms*rwmpgs*(-4.*xvq) dw2= dw2+ypivi*vvis*rwms*rwmpgs*(250./3.*xvq) dw2= dw2+ypivi*rwms*rwmpgs*(116./3.*xvq) dw2= dw2+ypiv*rlncu*vvis*rwms*rwmpgs*(72*xvc+68./3.*xvq) dw2= dw2+ypiv*rlncu*rwms*rwmpgs*(88./3.*xvc+20./3.*xvq) dw2= dw2+ypiv*rlncd*vvis*rwms*rwmpgs*(20*xvc+52./3.*xvq) dw2= dw2+ypiv*rlncd*rwms*rwmpgs*(32./3.*xvc+8*xvq) dw2= dw2+ypiv*vvis*rwms*rwmpgs*(-200*xvc-340./3.*xvq) dw2= dw2+ypiv*rwms*rwmpgs*(-310./3.*xvc-130./3.*xvq) dw2= dw2+ypiq*rlncu*vvis*rwms*rwmpgs*(-184./ # 3.*xvs-176./3.*xvc) dw2= dw2+ypiq*rlncu*rwms*rwmpgs*(-278./9.*xvs-52./3.*xvc) dw2= dw2+ypiq*rlncd*vvis*rwms*rwmpgs*(-49./3. # *xvs-136./3.*xvc-12*xvq) dw2= dw2+ypiq*rlncd*rwms*rwmpgs*(-98./9.*xvs # -62./3.*xvc-14./3.*xvq) dw2= dw2+ypiq*vvis*rwms*rwmpgs*(157*xvs+844./3.*xvc+42.*xvq) dw2= dw2+ypiq*rwms*rwmpgs*(916./9.*xvs+110*xvc+127./9.*xvq) dw2= dw2+ypic*rlncu*vvis*rwms*rwmpgs*(56./3.*xv+152./3.*xvs) dw2= dw2+ypic*rlncu*rwms*rwmpgs*(44./3.*xv+148./9.*xvs) dw2= dw2+ypic*rlncd*vvis*rwms*rwmpgs*(14./3.* # xv+119./3.*xvs+100./3.*xvc+8./3.*xvq) dw2= dw2+ypic*rlncd*rwms*rwmpgs*(5*xv+61./3.* # xvs+34./3.*xvc+2./3.*xvq) dw2= dw2+ypic*vvis*rwms*rwmpgs*(-131./3.*xv # -232.*xvs-112*xvc-20./3.*xvq) dw2= dw2+ypic*rwms*rwmpgs*(-383./9.*xv-907./9. # *xvs-33*xvc-5./3.*xvq) dw2= dw2+ypis*rlncu*vvis*rwms*rwmpgs*(-4./3.-16.*xv) dw2= dw2+ypis*rlncu*rwms*rwmpgs*(-22./9.-56./9.*xv) dw2= dw2+ypis*rlncd*vvis*rwms*rwmpgs*(-1./3. # -38./3.*xv-95./3.*xvs-8*xvc) dw2= dw2+ypis*rlncd*rwms*rwmpgs*(-7./9.-82./9. # *xv-91./9.*xvs-4./3.*xvc) dw2= dw2+ypis*vvis*rwms*rwmpgs*(17./6.+209./ # 3.*xv+203./2.*xvs+20*xvc+xvq) dw2= dw2+ypis*rwms*rwmpgs*(107./18.+341./9.*xv # +499./18.*xvs+10./3.*xvc) dw2= dw2+ypi*rlncu*vvis*rwms*rwmpgs*(4./3.) dw2= dw2+ypi*rlncu*rwms*rwmpgs*(4./9.) dw2= dw2+ypi*rlncd*vvis*rwms*rwmpgs*(1+34./3.*xv+25./3.*xvs) dw2= dw2+ypi*rlncd*rwms*rwmpgs*(13./9.+37./9.*xv+2./3.*xvs) dw2= dw2+ypi*vvis*rwms*rwmpgs*(-31./6.- # 104./3.*xv-125./6.*xvs-3*xvc) dw2= dw2+ypi*rwms*rwmpgs*(-83./18.-169./18.*xv-3*xvs+1./3.*xvc) dw2= dw2+rlncd*vvis*rwms*rwmpgs*omyp*(-1./3.) dw2= dw2+rlncd*vvis*rwms*rwmpgs*(-2./3.-10./3.*xv) dw2= dw2+rlncd*rwms*rwmpgs*(-2./3.) dw2= dw2+vvis*rwms*rwmpgs*omyp*(5./6.+xv) dw2= dw2+vvis*rwms*rwmpgs*(7./3.+22./3.*xv+3*xvs) dw2= dw2+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw2= dw2+rwms*rwmpgs*(11./9.+5./6.*xv-1./6.*xvs) dw2= dw2+denu*ymvvs*ypiq*vviq*rwms*rwmpgs*(-4./3.*xvq) dw2= dw2+denu*ymvvs*ypic*vviq*rwms*rwmpgs*(4.*xvc) dw2= dw2+denu*ymvvs*ypis*vviq*rwms*rwmpgs*(-4.*xvs) dw2= dw2+denu*ymvvs*ypis*vvis*rwms*rwmpgs*(-2./9.*xvs) dw2= dw2+denu*ymvvs*ypi*vviq*rwms*rwmpgs*(4./3.*xv) dw2= dw2+denu*ymvvs*ypi*vvis*rwms*rwmpgs*(4./9.*xv) dw2= dw2+denu*ymvvs*vvis*rwms*rwmpgs*(-2./9.) dw2= dw2+denu*ymvvc*ypiq*vviq*rwms*(8./3.*xvv) dw2= dw2+denu*ymvvc*ypic*vviq*rwms*(-8*xvq) dw2= dw2+denu*ymvvc*ypis*vviq*rwms*(8*xvc) dw2= dw2+denu*ymvvc*ypis*vvis*rwms*(4./9.*xvc) dw2= dw2+denu*ymvvc*ypi*vviq*rwms*(-8./3.*xvs) dw2= dw2+denu*ymvvc*ypi*vvis*rwms*(-8./9.*xvs) dw2= dw2+denu*ymvvc*vvis*rwms*(4./9.*xv) dw2= dw2+denu*ymvvq*ypiq*vviq*(-4./3.*xvvi) dw2= dw2+denu*ymvvq*ypic*vviq*(4.*xvv) dw2= dw2+denu*ymvvq*ypis*vviq*(-4.*xvq) dw2= dw2+denu*ymvvq*ypis*vvis*(-2./9.*xvq) dw2= dw2+denu*ymvvq*ypi*vviq*(4./3.*xvc) dw2= dw2+denu*ymvvq*ypi*vvis*(4./9.*xvc) dw2= dw2+denu*ymvvq*vvis*(-2./9.*xvs) dw2= dw2+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-1./3.*xvq) dw2= dw2+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(xvc+2./3.*xvq) dw2= dw2+dend*ymvvs*ypis*vviq*rwms*rwmpgs*(-xvs-2*xvc-1./3.*xvq) dw2= dw2+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./18.*xvs) dw2= dw2+dend*ymvvs*ypi*vviq*rwms*rwmpgs*(1./3.*xv+2*xvs+xvc) dw2= dw2+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(1./9.*xv+1./9.*xvs) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-1./3.*xv) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*(-1./3.*xv-xvs) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-1./9.*xv) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./18.) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*(-1./9.*xv-1./18.*xvs) dw2= dw2+dend*ymvvc*ypiq*vviq*rwms*(2./3.*xvv) dw2= dw2+dend*ymvvc*ypic*vviq*rwms*(-2*xvq-4./3.*xvv) dw2= dw2+dend*ymvvc*ypis*vviq*rwms*(2*xvc+4.*xvq+2./3.*xvv) dw2= dw2+dend*ymvvc*ypis*vvis*rwms*(1./9.*xvc) dw2= dw2+dend*ymvvc*ypi*vviq*rwms*(-2./3.*xvs-4.*xvc-2*xvq) dw2= dw2+dend*ymvvc*ypi*vvis*rwms*(-2./9.*xvs-2./9.*xvc) dw2= dw2+dend*ymvvc*vviq*rwms*omyp*(2./3.*xvs) dw2= dw2+dend*ymvvc*vviq*rwms*(2./3.*xvs+2*xvc) dw2= dw2+dend*ymvvc*vvis*rwms*omyp*(2./9.*xvs) dw2= dw2+dend*ymvvc*vvis*rwms*omyps*(1./9.*xv) dw2= dw2+dend*ymvvc*vvis*rwms*(2./9.*xvs+1./9.*xvc) dw2= dw2+dend*ymvvq*ypiq*vviq*(-1./3.*xvvi) dw2= dw2+dend*ymvvq*ypic*vviq*(xvv+2./3.*xvvi) dw2= dw2+dend*ymvvq*ypis*vviq*(-xvq-2*xvv-1./3.*xvvi) dw2= dw2+dend*ymvvq*ypis*vvis*(-1./18.*xvq) dw2= dw2+dend*ymvvq*ypi*vviq*(1./3.*xvc+2*xvq+xvv) dw2= dw2+dend*ymvvq*ypi*vvis*(1./9.*xvc+1./9.*xvq) dw2= dw2+dend*ymvvq*vviq*omyp*(-1./3.*xvc) dw2= dw2+dend*ymvvq*vviq*(-1./3.*xvc-xvq) dw2= dw2+dend*ymvvq*vvis*omyp*(-1./9.*xvc) dw2= dw2+dend*ymvvq*vvis*omyps*(-1./18.*xvs) dw2= dw2+dend*ymvvq*vvis*(-1./9.*xvc-1./18.*xvq) * if(ii.eq.3) then wcom= wcom*(2.d0-ym)**2/yms-dw1+2.d0/ym*(1.d0-1.d0/ym)*dw2 else if(ii.eq.4) then wcom= 0.5d0*em4*omxm*omxm*(dw1-dw2) else if(ii.eq.5) then wcom= em2*((2.d0+0.5d0*omxm)*dw1+omxm*(2.d0/ym- # 3.d0/2.d0)*dw2) endif * else if(ipro.eq.1.and.ii.le.2) then * wcomd1= ypiq*rwms*(-2*xvq) wcomd1= wcomd1+ypic*rwms*(4*xvc+4*xvq) wcomd1= wcomd1+ypis*rwms*(-3*xvs-8*xvc-2*xvq) wcomd1= wcomd1+ypi*rwms*(xv+6*xvs+4*xvc) wcomd1= wcomd1+rwms*omyp*(-xv) wcomd1= wcomd1+rwms*(-xv-3*xvs) wcomd1= wcomd1+ymvv*ypiq*(xvv) wcomd1= wcomd1+ymvv*ypic*(-2*xvq-2*xvv) wcomd1= wcomd1+ymvv*ypis*(3./2.*xvc+4*xvq+xvv) wcomd1= wcomd1+ymvv*ypi*(-1./2.*xvs-3*xvc-2*xvq) wcomd1= wcomd1+ymvv*omyp*(1./2.*xvs) wcomd1= wcomd1+ymvv*(1./2.*xvs+3./2.*xvc) wcomd1= ymvv*wcomd1 * wcomd2= ypiq*(xvc) wcomd2= wcomd2+ypic*(-2*xvs-2*xvc) wcomd2= wcomd2+ypis*(3./2.*xv+4*xvs+xvc) wcomd2= wcomd2+ypi*(-1./2.-3*xv-2*xvs) wcomd2= wcomd2+omyp*(1./2.) wcomd2= wcomd2+(1./2.+3./2.*xv) wcomd2= rwms*rwmpgs*wcomd2 * wcom01= ypiq*rwms*(22./3.*xvq) wcom01= wcom01+ypic*rwms*(-44./3.*xvc-14*xvq) wcom01= wcom01+ypis*rwms*(10*xvs+28*xvc+7*xvq) wcom01= wcom01+ypi*rwms*(-8./3.*xv-20*xvs-14*xvc-2./3.*xvq) wcom01= wcom01+rwms*omyp*(17./6.*xv+xvs) wcom01= wcom01+rwms*omyps*(1./3.*xv) wcom01= wcom01+rwms*(17./6.*xv+19./2.*xvs+4./3.*xvc) wcom01= wcom01+ymvv*ypiq*(-11./3.*xvv) wcom01= wcom01+ymvv*ypic*(23./3.*xvq+7*xvv) wcom01= wcom01+ymvv*ypis*(-11./2.*xvc-15*xvq-7./2.*xvv) wcom01= wcom01+ymvv*ypi*(11./6.*xvs+23./2.*xvc+8*xvq+1./3.*xvv) wcom01= wcom01+ymvv*omyp*(-1./6.*xv-19./12.*xvs-7./6.*xvc) wcom01= wcom01+ymvv*omyps*(-2./3.*xvs) wcom01= wcom01+ymvv*omypc*(-1./6.*xv) wcom01= wcom01+ymvv*(-23./12.*xvs-23./4.*xvc-xvq) wcom01= ymvv*wcom01 * wcom02= ypiq*(-11./3.*xvc) wcom02= wcom02+ypic*(7*xvs+7*xvc) wcom02= wcom02+ypis*(-9./2.*xv-13*xvs-7./2.*xvc) wcom02= wcom02+ypi*(5./6.+17./2.*xv+6*xvs+1./3.*xvc) wcom02= wcom02+omyp*(-11./12.-1./6.*xv) wcom02= wcom02+(-11./12.-15./4.*xv-1./3.*xvs) wcom02= rwms*rwmpgs*wcom02 * wcom0= wcom01+wcom02 * wcomdm= ypis*vvi*rwms*rwmpgs*xvs wcomdm= wcomdm+ypi*vvi*rwms*rwmpgs*(-2*xv-2*xvs) wcomdm= wcomdm+vvi*rwms*rwmpgs*omyp*(2*xv) wcomdm= wcomdm+vvi*rwms*rwmpgs*omyps wcomdm= wcomdm+vvi*rwms*rwmpgs*(2*xv+xvs) wcomdm= wcomdm+ymvv*ypis*vvi*rwms*(-2*xvc) wcomdm= wcomdm+ymvv*ypi*vvi*rwms*(4*xvs+4*xvc) wcomdm= wcomdm+ymvv*vvi*rwms*omyp*(-4*xvs) wcomdm= wcomdm+ymvv*vvi*rwms*omyps*(-2*xv) wcomdm= wcomdm+ymvv*vvi*rwms*(-4*xvs-2*xvc) wcomdm= wcomdm+ymvvs*ypis*vvi*xvq wcomdm= wcomdm+ymvvs*ypi*vvi*(-2*xvc-2*xvq) wcomdm= wcomdm+ymvvs*vvi*omyp*(2*xvc) wcomdm= wcomdm+ymvvs*vvi*omyps*xvs wcomdm= wcomdm+ymvvs*vvi*(2*xvc+xvq) dendm= dend*dms/vv wcomdm= ymvv*dendm*wcomdm * wcom= wcomd*rlncd+wcom0+wcomdm * else if(ipro.eq.1.and.ii.ge.3) then * wcomd1= ypiq*rwms*(-2*xvq) wcomd1= wcomd1+ypic*rwms*(4*xvc+4*xvq) wcomd1= wcomd1+ypis*rwms*(-3*xvs-8*xvc-2*xvq) wcomd1= wcomd1+ypi*rwms*(xv+6*xvs+4*xvc) wcomd1= wcomd1+rwms*omyp*(-xv) wcomd1= wcomd1+rwms*(-xv-3*xvs) wcomd1= wcomd1+ymvv*ypiq*(xvv) wcomd1= wcomd1+ymvv*ypic*(-2*xvq-2*xvv) wcomd1= wcomd1+ymvv*ypis*(3./2.*xvc+4*xvq+xvv) wcomd1= wcomd1+ymvv*ypi*(-1./2.*xvs-3*xvc-2*xvq) wcomd1= wcomd1+ymvv*omyp*(1./2.*xvs) wcomd1= wcomd1+ymvv*(1./2.*xvs+3./2.*xvc) wcomd1= ymvv*wcomd1 * wcomd2= ypiq*(xvc) wcomd2= wcomd2+ypic*(-2*xvs-2*xvc) wcomd2= wcomd2+ypis*(3./2.*xv+4*xvs+xvc) wcomd2= wcomd2+ypi*(-1./2.-3*xv-2*xvs) wcomd2= wcomd2+omyp*(1./2.) wcomd2= wcomd2+(1./2.+3./2.*xv) wcomd2= rwms*rwmpgs*wcomd2 * wcom01= ypiq*rwms*(22./3.*xvq) wcom01= wcom01+ypic*rwms*(-44./3.*xvc-14*xvq) wcom01= wcom01+ypis*rwms*(10*xvs+28*xvc+7*xvq) wcom01= wcom01+ypi*rwms*(-8./3.*xv-20*xvs-14*xvc-2./3.*xvq) wcom01= wcom01+rwms*omyp*(17./6.*xv+xvs) wcom01= wcom01+rwms*omyps*(1./3.*xv) wcom01= wcom01+rwms*(17./6.*xv+19./2.*xvs+4./3.*xvc) wcom01= wcom01+ymvv*ypiq*(-11./3.*xvv) wcom01= wcom01+ymvv*ypic*(23./3.*xvq+7*xvv) wcom01= wcom01+ymvv*ypis*(-11./2.*xvc-15*xvq-7./2.*xvv) wcom01= wcom01+ymvv*ypi*(11./6.*xvs+23./2.*xvc+8*xvq+1./3.*xvv) wcom01= wcom01+ymvv*omyp*(-1./6.*xv-19./12.*xvs-7./6.*xvc) wcom01= wcom01+ymvv*omyps*(-2./3.*xvs) wcom01= wcom01+ymvv*omypc*(-1./6.*xv) wcom01= wcom01+ymvv*(-23./12.*xvs-23./4.*xvc-xvq) wcom01= ymvv*wcom01 * wcom02= ypiq*(-11./3.*xvc) wcom02= wcom02+ypic*(7*xvs+7*xvc) wcom02= wcom02+ypis*(-9./2.*xv-13*xvs-7./2.*xvc) wcom02= wcom02+ypi*(5./6.+17./2.*xv+6*xvs+1./3.*xvc) wcom02= wcom02+omyp*(-11./12.-1./6.*xv) wcom02= wcom02+(-11./12.-15./4.*xv-1./3.*xvs) wcom02= rwms*rwmpgs*wcom02 * wcom0= wcom01+wcom02 * wcomdm= ypis*vvi*rwms*rwmpgs*xvs wcomdm= wcomdm+ypi*vvi*rwms*rwmpgs*(-2*xv-2*xvs) wcomdm= wcomdm+vvi*rwms*rwmpgs*omyp*(2*xv) wcomdm= wcomdm+vvi*rwms*rwmpgs*omyps wcomdm= wcomdm+vvi*rwms*rwmpgs*(2*xv+xvs) wcomdm= wcomdm+ymvv*ypis*vvi*rwms*(-2*xvc) wcomdm= wcomdm+ymvv*ypi*vvi*rwms*(4*xvs+4*xvc) wcomdm= wcomdm+ymvv*vvi*rwms*omyp*(-4*xvs) wcomdm= wcomdm+ymvv*vvi*rwms*omyps*(-2*xv) wcomdm= wcomdm+ymvv*vvi*rwms*(-4*xvs-2*xvc) wcomdm= wcomdm+ymvvs*ypis*vvi*xvq wcomdm= wcomdm+ymvvs*ypi*vvi*(-2*xvc-2*xvq) wcomdm= wcomdm+ymvvs*vvi*omyp*(2*xvc) wcomdm= wcomdm+ymvvs*vvi*omyps*xvs wcomdm= wcomdm+ymvvs*vvi*(2*xvc+xvq) dendm= dend*dms/vv wcomdm= ymvv*dendm*wcomdm * wcom= wcomd*rlncd+wcom0+wcomdm * dw1= ymvv*ypivi*rlncd*vvis*rwms*(40*xvv) dw1= dw1+ymvv*ypivi*rlncd*rwms*(40*xvv) dw1= dw1+ymvv*ypivi*vvis*rwms*(-260./3.*xvv) dw1= dw1+ymvv*ypivi*rwms*(-296./3.*xvv) dw1= dw1+ymvv*ypiv*rlncd*vvis*rwms*(-104*xvq-80*xvv) dw1= dw1+ymvv*ypiv*rlncd*rwms*(-120*xvq-72*xvv) dw1= dw1+ymvv*ypiv*vvis*rwms*(640./3.*xvq+520./3.*xvv) dw1= dw1+ymvv*ypiv*rwms*(296*xvq+172*xvv) dw1= dw1+ymvv*ypiq*rlncd*vvis*rwms*(89*xvc # +218*xvq+48*xvv) dw1= dw1+ymvv*ypiq*rlncd*rwms*(141*xvc+214*xvq+36*xvv) dw1= dw1+ymvv*ypiq*vvis*rwms*(-511./3.*xvc # -1366./3.*xvq-108*xvv) dw1= dw1+ymvv*ypiq*rwms*(-339*xvc-1526./3.*xvq-86*xvv) dw1= dw1+ymvv*ypic*rlncd*vvis*rwms*(-26* # xvs-202*xvc-140*xvq-8*xvv) dw1= dw1+ymvv*ypic*rlncd*rwms*(-78*xvs-252* # xvc-104*xvq-4*xvv) dw1= dw1+ymvv*ypic*vvis*rwms*(133./3.*xvs+ # 1217./3.*xvc+310*xvq+64./3.*xvv) dw1= dw1+ymvv*ypic*rwms*(535./3.*xvs+1753./3.* # xvc+244*xvq+34./3.*xvv) dw1= dw1+ymvv*ypis*rlncd*vvis*rwms*(xv+70* # xvs+145*xvc+26*xvq) dw1= dw1+ymvv*ypis*rlncd*rwms*(18*xv+142*xvs # +119*xvc+10*xvq) dw1= dw1+ymvv*ypis*vvis*rwms*(-1./6.*xv- # 805./6.*xvs-635./2.*xvc-209./3.*xvq-2./3.*xvv) dw1= dw1+ymvv*ypis*rwms*(-81./2.*xv-316*xvs # -545./2.*xvc-27*xvq) dw1= dw1+ymvv*ypi*rlncd*vvis*rwms*(-6*xv-62*xvs-32*xvc) dw1= dw1+ymvv*ypi*rlncd*rwms*(-1-34*xv-66*xvs-8*xvc) dw1= dw1+ymvv*ypi*vvis*rwms*(-1./3.+29./3. # *xv+135*xvs+257./3.*xvc+8./3.*xvq) dw1= dw1+ymvv*ypi*rwms*(3+229./3.*xv+147* # xvs+25*xvc-2./3.*xvq) dw1= dw1+ymvv*rlncd*vvis*rwms*omyp*(4*xv) dw1= dw1+ymvv*rlncd*vvis*rwms*(5*xv+18*xvs) dw1= dw1+ymvv*rlncd*rwms*omyp*(1) dw1= dw1+ymvv*rlncd*rwms*(1+16*xv+2*xvs) dw1= dw1+ymvv*vvis*rwms*omyp*(1./3.-9*xv-8./3.*xvs) dw1= dw1+ymvv*vvis*rwms*omyps*(-2./3.*xv) dw1= dw1+ymvv*vvis*rwms*(1./3.-19./2.*xv-45*xvs-4*xvc) dw1= dw1+ymvv*rwms*omyp*(-19./6.-3*xv+2./3.*xvs) dw1= dw1+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw1= dw1+ymvv*rwms*(-10./3.-34*xv-12*xvs+xvc) dw1= dw1+ymvvs*ypivi*rlncd*vvis*(-20*xvvi) dw1= dw1+ymvvs*ypivi*rlncd*(-20*xvvi) dw1= dw1+ymvvs*ypivi*vvis*(130./3.*xvvi) dw1= dw1+ymvvs*ypivi*(148./3.*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*vvis*(52*xvv+40*xvvi) dw1= dw1+ymvvs*ypiv*rlncd*(64*xvv+36*xvvi) dw1= dw1+ymvvs*ypiv*vvis*(-320./3.*xvv-260./3.*xvvi) dw1= dw1+ymvvs*ypiv*(-482./3.*xvv-86*xvvi) dw1= dw1+ymvvs*ypiq*rlncd*vvis*(-44*xvq-110*xvv-24*xvvi) dw1= dw1+ymvvs*ypiq*rlncd*(-81*xvq-116*xvv-18*xvvi) dw1= dw1+ymvvs*ypiq*vvis*(250./3.*xvq+694./ # 3.*xvv+54*xvvi) dw1= dw1+ymvvs*ypiq*(203*xvq+844./3.*xvv+43*xvvi) dw1= dw1+ymvvs*ypic*rlncd*vvis*(12*xvc+103 # *xvq+72*xvv+4*xvvi) dw1= dw1+ymvvs*ypic*rlncd*(49*xvc+149*xvq+58*xvv+2*xvvi) dw1= dw1+ymvvs*ypic*vvis*(-56./3.*xvc- # 629./3.*xvq-162*xvv-32./3.*xvvi) dw1= dw1+ymvvs*ypic*(-364./3.*xvc-1084./3.* # xvq-139*xvv-17./3.*xvvi) dw1= dw1+ymvvs*ypis*rlncd*vvis*(-36*xvc-78*xvq-14*xvv) dw1= dw1+ymvvs*ypis*rlncd*(-13*xvs-92*xvc-74*xvq-6*xvv) dw1= dw1+ymvvs*ypis*vvis*(-5./3.*xvs+419. # /6.*xvc+177*xvq+115./3.*xvv+1./3.*xvvi) dw1= dw1+ymvvs*ypis*(100./3.*xvs+223*xvc+ # 177*xvq+49./3.*xvv) dw1= dw1+ymvvs*ypi*rlncd*vvis*(3*xvs+36*xvc+19*xvq) dw1= dw1+ymvvs*ypi*rlncd*(xv+25*xvs+45*xvc+6*xvq) dw1= dw1+ymvvs*ypi*vvis*(1./3.*xv-25./6.* # xvs-83*xvc-319./6.*xvq-5./3.*xvv) dw1= dw1+ymvvs*ypi*(-4*xv-391./6.*xvs-221. # /2.*xvc-18*xvq+1./3.*xvv) dw1= dw1+ymvvs*rlncd*vvis*omyp*(-3*xvs) dw1= dw1+ymvvs*rlncd*vvis*(-3*xvs-12*xvc) dw1= dw1+ymvvs*rlncd*omyp*(-xv) dw1= dw1+ymvvs*rlncd*(-xv-12*xvs-2*xvc) dw1= dw1+ymvvs*vvis*omyp*(-1./24.+1./12.*xv+ # 35./6.*xvs+10./3.*xvc) dw1= dw1+ymvvs*vvis*omyps*(-2./3.*xv+5./3.*xvs) dw1= dw1+ymvvs*vvis*omypc*(1./3.*xv) dw1= dw1+ymvvs*vvis*(-5./12.*xv+25./4.*xvs # +31*xvc+10./3.*xvq) dw1= dw1+ymvvs*omyp*(5./24.+13./4.*xv+29./6.*xvs # -5./6.*xvc) dw1= dw1+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw1= dw1+ymvvs*omypc*(1./6.-1./6.*xv) dw1= dw1+ymvvs*(17./4.*xv+365./12.*xvs+67./6.* # xvc-5./6.*xvq) dw1= dw1+ypivi*rlncd*vvis*rwms*rwmpgs*(-20*xvq) dw1= dw1+ypivi*rlncd*rwms*rwmpgs*(-20*xvq) dw1= dw1+ypivi*vvis*rwms*rwmpgs*(130./3.*xvq) dw1= dw1+ypivi*rwms*rwmpgs*(148./3.*xvq) dw1= dw1+ypiv*rlncd*vvis*rwms*rwmpgs*(52*xvc+40*xvq) dw1= dw1+ypiv*rlncd*rwms*rwmpgs*(56*xvc+36*xvq) dw1= dw1+ypiv*vvis*rwms*rwmpgs*(-320./3.* # xvc-260./3.*xvq) dw1= dw1+ypiv*rwms*rwmpgs*(-406./3.*xvc-86*xvq) dw1= dw1+ypiq*rlncd*vvis*rwms*rwmpgs*(-45* # xvs-108*xvc-24*xvq) dw1= dw1+ypiq*rlncd*rwms*rwmpgs*(-60*xvs-98* # xvc-18*xvq) dw1= dw1+ypiq*vvis*rwms*rwmpgs*(87*xvs+224*xvc+54*xvq) dw1= dw1+ypiq*rwms*rwmpgs*(136*xvs+682./3.*xvc+43*xvq) dw1= dw1+ypic*rlncd*vvis*rwms*rwmpgs*(14*xv+ # 99*xvs+68*xvc+4*xvq) dw1= dw1+ypic*rlncd*rwms*rwmpgs*(29*xv+103*xvs # +46*xvc+2*xvq) dw1= dw1+ypic*vvis*rwms*rwmpgs*(-77./3.*xv # -196*xvs-148*xvc-32./3.*xvq) dw1= dw1+ypic*rwms*rwmpgs*(-173./3.*xv-223* # xvs-105*xvc-17./3.*xvq) dw1= dw1+ypis*rlncd*vvis*rwms*rwmpgs*(-1- # 34*xv-67*xvs-12*xvc) dw1= dw1+ypis*rlncd*rwms*rwmpgs*(-5-50*xv-45*xvs-4*xvc) dw1= dw1+ypis*vvis*rwms*rwmpgs*(11./6.+193./ # 3.*xv+281./2.*xvs+94./3.*xvc+1./3.*xvq) dw1= dw1+ypis*rwms*rwmpgs*(49./6.+95*xv+191./2. # *xvs+32./3.*xvc) dw1= dw1+ypi*rlncd*vvis*rwms*rwmpgs*(3+26*xv+13*xvs) dw1= dw1+ypi*rlncd*rwms*rwmpgs*(9+21*xv+2*xvs) dw1= dw1+ypi*vvis*rwms*rwmpgs*(-11./2.-52 # *xv-65./2.*xvs-xvc) dw1= dw1+ypi*rwms*rwmpgs*(-85./6.-77./2.*xv # -7*xvs+1./3.*xvc) dw1= dw1+rlncd*vvis*rwms*rwmpgs*omyp*(-1) dw1= dw1+rlncd*vvis*rwms*rwmpgs*(-2-6*xv) dw1= dw1+rlncd*rwms*rwmpgs*(-4) dw1= dw1+vvis*rwms*rwmpgs*omyp*(11./6.+1./3.*xv) dw1= dw1+vvis*rwms*rwmpgs*(11./3.+40./3.*xv+xvs) dw1= dw1+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw1= dw1+rwms*rwmpgs*(17./3.+3./2.*xv-1./6.*xvs) dw1= dw1+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-xvq) dw1= dw1+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(3*xvc+2*xvq) dw1= dw1+dend*ymvvs*ypis*vviq*rwms*rwmpgs* # (-3*xvs-6*xvc-xvq) dw1= dw1+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./2.*xvs) dw1= dw1+dend*ymvvs*ypi*vviq*rwms*rwmpgs* # (xv+6*xvs+3*xvc) dw1= dw1+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(xv+xvs) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-xv) dw1= dw1+dend*ymvvs*vviq*rwms*rwmpgs*(-xv-3*xvs) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-xv) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./2.) dw1= dw1+dend*ymvvs*vvis*rwms*rwmpgs*(-xv-1./2.*xvs) dw1= dw1+dend*ymvvc*ypiq*vviq*rwms*(2*xvv) dw1= dw1+dend*ymvvc*ypic*vviq*rwms*(-6*xvq-4*xvv) dw1= dw1+dend*ymvvc*ypis*vviq*rwms*(6*xvc+12*xvq+2*xvv) dw1= dw1+dend*ymvvc*ypis*vvis*rwms*xvc dw1= dw1+dend*ymvvc*ypi*vviq*rwms*(-2*xvs-12*xvc-6*xvq) dw1= dw1+dend*ymvvc*ypi*vvis*rwms*(-2*xvs-2*xvc) dw1= dw1+dend*ymvvc*vviq*rwms*omyp*(2*xvs) dw1= dw1+dend*ymvvc*vviq*rwms*(2*xvs+6*xvc) dw1= dw1+dend*ymvvc*vvis*rwms*omyp*(2*xvs) dw1= dw1+dend*ymvvc*vvis*rwms*omyps*(xv) dw1= dw1+dend*ymvvc*vvis*rwms*(2*xvs+xvc) dw1= dw1+dend*ymvvq*ypiq*vviq*(-xvvi) dw1= dw1+dend*ymvvq*ypic*vviq*(3*xvv+2*xvvi) dw1= dw1+dend*ymvvq*ypis*vviq*(-3*xvq-6*xvv-xvvi) dw1= dw1+dend*ymvvq*ypis*vvis*(-1./2.*xvq) dw1= dw1+dend*ymvvq*ypi*vviq*(xvc+6*xvq+3*xvv) dw1= dw1+dend*ymvvq*ypi*vvis*(xvc+xvq) dw1= dw1+dend*ymvvq*vviq*omyp*(-xvc) dw1= dw1+dend*ymvvq*vviq*(-xvc-3*xvq) dw1= dw1+dend*ymvvq*vvis*omyp*(-xvc) dw1= dw1+dend*ymvvq*vvis*omyps*(-1./2.*xvs) dw1= dw1+dend*ymvvq*vvis*(-xvc-1./2.*xvq) dw2= ymvv*ypivi*rlncd*vvis*rwms*(120*xvv) dw2= dw2+ymvv*ypivi*rlncd*rwms*(40*xvv) dw2= dw2+ymvv*ypivi*vvis*rwms*(-260*xvv) dw2= dw2+ymvv*ypivi*rwms*(-296./3.*xvv) dw2= dw2+ymvv*ypiv*rlncd*vvis*rwms*(-312*xvq-240*xvv) dw2= dw2+ymvv*ypiv*rlncd*rwms*(-120*xvq-72*xvv) dw2= dw2+ymvv*ypiv*vvis*rwms*(640*xvq+520*xvv) dw2= dw2+ymvv*ypiv*rwms*(296*xvq+172*xvv) dw2= dw2+ymvv*ypiq*rlncd*vvis*rwms*(267*xvc # +654*xvq+144*xvv) dw2= dw2+ymvv*ypiq*rlncd*rwms*(141*xvc+214*xvq+36*xvv) dw2= dw2+ymvv*ypiq*vvis*rwms*(-511*xvc- # 1366*xvq-324*xvv) dw2= dw2+ymvv*ypiq*rwms*(-339*xvc-1526./3.*xvq-86*xvv) dw2= dw2+ymvv*ypic*rlncd*vvis*rwms*(-78* # xvs-606*xvc-420*xvq-24*xvv) dw2= dw2+ymvv*ypic*rlncd*rwms*(-78*xvs-252* # xvc-104*xvq-4*xvv) dw2= dw2+ymvv*ypic*vvis*rwms*(133*xvs+1217 # *xvc+930*xvq+64*xvv) dw2= dw2+ymvv*ypic*rwms*(535./3.*xvs+1753./3.* # xvc+244*xvq+34./3.*xvv) dw2= dw2+ymvv*ypis*rlncd*vvis*rwms*(3*xv+210 # *xvs+435*xvc+78*xvq) dw2= dw2+ymvv*ypis*rlncd*rwms*(18*xv+142*xvs # +119*xvc+10*xvq) dw2= dw2+ymvv*ypis*vvis*rwms*(-1./2.*xv- # 805./2.*xvs-1905./2.*xvc-209*xvq-2*xvv) dw2= dw2+ymvv*ypis*rwms*(-81./2.*xv-316*xvs # -545./2.*xvc-27*xvq) dw2= dw2+ymvv*ypi*rlncd*vvis*rwms*(-18*xv-186*xvs-96*xvc) dw2= dw2+ymvv*ypi*rlncd*rwms*(-1-34*xv-66*xvs-8*xvc) dw2= dw2+ymvv*ypi*vvis*rwms*(-1+29*xv+ # 405*xvs+257*xvc+8*xvq) dw2= dw2+ymvv*ypi*rwms*(3+229./3.*xv+147* # xvs+25*xvc-2./3.*xvq) dw2= dw2+ymvv*rlncd*vvis*rwms*omyp*(12*xv) dw2= dw2+ymvv*rlncd*vvis*rwms*(15*xv+54*xvs) dw2= dw2+ymvv*rlncd*rwms*omyp*(1) dw2= dw2+ymvv*rlncd*rwms*(1+16*xv+2*xvs) dw2= dw2+ymvv*vvis*rwms*omyp*(1-27*xv-8*xvs) dw2= dw2+ymvv*vvis*rwms*omyps*(-2*xv) dw2= dw2+ymvv*vvis*rwms*(1-57./2.*xv-135*xvs-12*xvc) dw2= dw2+ymvv*rwms*omyp*(-19./6.-3*xv+2./3.*xvs) dw2= dw2+ymvv*rwms*omyps*(-1./3.+1./3.*xv) dw2= dw2+ymvv*rwms*(-10./3.-34*xv-12*xvs+xvc) dw2= dw2+ymvvs*ypivi*rlncd*vvis*(-60*xvvi) dw2= dw2+ymvvs*ypivi*rlncd*(-20*xvvi) dw2= dw2+ymvvs*ypivi*vvis*(130*xvvi) dw2= dw2+ymvvs*ypivi*(148./3.*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*vvis*(156*xvv+120*xvvi) dw2= dw2+ymvvs*ypiv*rlncd*(64*xvv+36*xvvi) dw2= dw2+ymvvs*ypiv*vvis*(-320*xvv-260*xvvi) dw2= dw2+ymvvs*ypiv*(-482./3.*xvv-86*xvvi) dw2= dw2+ymvvs*ypiq*rlncd*vvis*(-132*xvq # -330*xvv-72*xvvi) dw2= dw2+ymvvs*ypiq*rlncd*(-81*xvq-116*xvv-18*xvvi) dw2= dw2+ymvvs*ypiq*vvis*(250*xvq+694* # xvv+162*xvvi) dw2= dw2+ymvvs*ypiq*(203*xvq+844./3.*xvv+43*xvvi) dw2= dw2+ymvvs*ypic*rlncd*vvis*(36*xvc+309 # *xvq+216*xvv+12*xvvi) dw2= dw2+ymvvs*ypic*rlncd*(49*xvc+149*xvq+58*xvv+2*xvvi) dw2= dw2+ymvvs*ypic*vvis*(-56*xvc-629* # xvq-486*xvv-32*xvvi) dw2= dw2+ymvvs*ypic*(-364./3.*xvc-1084./3.* # xvq-139*xvv-17./3.*xvvi) dw2= dw2+ymvvs*ypis*rlncd*vvis*(-108*xvc # -234*xvq-42*xvv) dw2= dw2+ymvvs*ypis*rlncd*(-13*xvs-92*xvc-74*xvq-6*xvv) dw2= dw2+ymvvs*ypis*vvis*(-5*xvs+419./2. # *xvc+531*xvq+115*xvv+xvvi) dw2= dw2+ymvvs*ypis*(100./3.*xvs+223*xvc+ # 177*xvq+49./3.*xvv) dw2= dw2+ymvvs*ypi*rlncd*vvis*(9*xvs+108*xvc+57*xvq) dw2= dw2+ymvvs*ypi*rlncd*(xv+25*xvs+45*xvc+6*xvq) dw2= dw2+ymvvs*ypi*vvis*(xv-25./2.*xvs- # 249*xvc-319./2.*xvq-5*xvv) dw2= dw2+ymvvs*ypi*(-4*xv-391./6.*xvs-221. # /2.*xvc-18*xvq+1./3.*xvv) dw2= dw2+ymvvs*rlncd*vvis*omyp*(-9*xvs) dw2= dw2+ymvvs*rlncd*vvis*(-9*xvs-36*xvc) dw2= dw2+ymvvs*rlncd*omyp*(-xv) dw2= dw2+ymvvs*rlncd*(-xv-12*xvs-2*xvc) dw2= dw2+ymvvs*vvis*omyp*(-1./8.+1./4.*xv+35. # /2.*xvs+10*xvc) dw2= dw2+ymvvs*vvis*omyps*(-2*xv+5*xvs) dw2= dw2+ymvvs*vvis*omypc*(xv) dw2= dw2+ymvvs*vvis*(-5./4.*xv+75./4.*xvs+ # 93*xvc+10*xvq) dw2= dw2+ymvvs*omyp*(5./24.+13./4.*xv+29./6.*xvs-5./6.*xvc) dw2= dw2+ymvvs*omyps*(3./2.*xv-1./2.*xvs) dw2= dw2+ymvvs*omypc*(1./6.-1./6.*xv) dw2= dw2+ymvvs*(17./4.*xv+365./12.*xvs+67./6.* # xvc-5./6.*xvq) dw2= dw2+ypivi*rlncd*vvis*rwms*rwmpgs*(-60*xvq) dw2= dw2+ypivi*rlncd*rwms*rwmpgs*(-20*xvq) dw2= dw2+ypivi*vvis*rwms*rwmpgs*(130*xvq) dw2= dw2+ypivi*rwms*rwmpgs*(148./3.*xvq) dw2= dw2+ypiv*rlncd*vvis*rwms*rwmpgs*(156*xvc+120*xvq) dw2= dw2+ypiv*rlncd*rwms*rwmpgs*(56*xvc+36*xvq) dw2= dw2+ypiv*vvis*rwms*rwmpgs*(-320*xvc-260*xvq) dw2= dw2+ypiv*rwms*rwmpgs*(-406./3.*xvc-86*xvq) dw2= dw2+ypiq*rlncd*vvis*rwms*rwmpgs*(-135*xvs-324*xvc-72*xvq) dw2= dw2+ypiq*rlncd*rwms*rwmpgs*(-60*xvs-98*xvc-18*xvq) dw2= dw2+ypiq*vvis*rwms*rwmpgs*(261*xvs+ # 672*xvc+162*xvq) dw2= dw2+ypiq*rwms*rwmpgs*(136*xvs+682./3.*xvc+43*xvq) dw2= dw2+ypic*rlncd*vvis*rwms*rwmpgs*(42*xv+ # 297*xvs+204*xvc+12*xvq) dw2= dw2+ypic*rlncd*rwms*rwmpgs*(29*xv+103*xvs+46*xvc+2*xvq) dw2= dw2+ypic*vvis*rwms*rwmpgs*(-77*xv-588 # *xvs-444*xvc-32*xvq) dw2= dw2+ypic*rwms*rwmpgs*(-173./3.*xv-223* # xvs-105*xvc-17./3.*xvq) dw2= dw2+ypis*rlncd*vvis*rwms*rwmpgs*(-3- # 102*xv-201*xvs-36*xvc) dw2= dw2+ypis*rlncd*rwms*rwmpgs*(-5-50*xv-45*xvs-4*xvc) dw2= dw2+ypis*vvis*rwms*rwmpgs*(11./2.+193* # xv+843./2.*xvs+94*xvc+xvq) dw2= dw2+ypis*rwms*rwmpgs*(49./6.+95*xv+191./2. # *xvs+32./3.*xvc) dw2= dw2+ypi*rlncd*vvis*rwms*rwmpgs*(9+78*xv+39*xvs) dw2= dw2+ypi*rlncd*rwms*rwmpgs*(9+21*xv+2*xvs) dw2= dw2+ypi*vvis*rwms*rwmpgs*(-33./2.- # 156*xv-195./2.*xvs-3*xvc) dw2= dw2+ypi*rwms*rwmpgs*(-85./6.-77./2.*xv # -7*xvs+1./3.*xvc) dw2= dw2+rlncd*vvis*rwms*rwmpgs*omyp*(-3) dw2= dw2+rlncd*vvis*rwms*rwmpgs*(-6-18*xv) dw2= dw2+rlncd*rwms*rwmpgs*(-4) dw2= dw2+vvis*rwms*rwmpgs*omyp*(11./2.+xv) dw2= dw2+vvis*rwms*rwmpgs*(11+40*xv+3*xvs) dw2= dw2+rwms*rwmpgs*omyp*(1./6.-1./6.*xv) dw2= dw2+rwms*rwmpgs*(17./3.+3./2.*xv-1./6.*xvs) dw2= dw2+dend*ymvvs*ypiq*vviq*rwms*rwmpgs*(-3*xvq) dw2= dw2+dend*ymvvs*ypic*vviq*rwms*rwmpgs*(9*xvc+6*xvq) dw2= dw2+dend*ymvvs*ypis*vviq*rwms*rwmpgs*(-9*xvs-18*xvc-3*xvq) dw2= dw2+dend*ymvvs*ypis*vvis*rwms*rwmpgs*(-1./2.*xvs) dw2= dw2+dend*ymvvs*ypi*vviq*rwms*rwmpgs*(3*xv+18*xvs+9*xvc) dw2= dw2+dend*ymvvs*ypi*vvis*rwms*rwmpgs*(xv+xvs) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*omyp*(-3*xv) dw2= dw2+dend*ymvvs*vviq*rwms*rwmpgs*(-3*xv-9*xvs) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyp*(-xv) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*omyps*(-1./2.) dw2= dw2+dend*ymvvs*vvis*rwms*rwmpgs*(-xv-1./2.*xvs) dw2= dw2+dend*ymvvc*ypiq*vviq*rwms*(6*xvv) dw2= dw2+dend*ymvvc*ypic*vviq*rwms*(-18*xvq-12*xvv) dw2= dw2+dend*ymvvc*ypis*vviq*rwms*(18*xvc+36*xvq+6*xvv) dw2= dw2+dend*ymvvc*ypis*vvis*rwms*(xvc) dw2= dw2+dend*ymvvc*ypi*vviq*rwms*(-6*xvs-36*xvc-18*xvq) dw2= dw2+dend*ymvvc*ypi*vvis*rwms*(-2*xvs-2*xvc) dw2= dw2+dend*ymvvc*vviq*rwms*omyp*(6*xvs) dw2= dw2+dend*ymvvc*vviq*rwms*(6*xvs+18*xvc) dw2= dw2+dend*ymvvc*vvis*rwms*omyp*(2*xvs) dw2= dw2+dend*ymvvc*vvis*rwms*omyps*(xv) dw2= dw2+dend*ymvvc*vvis*rwms*(2*xvs+xvc) dw2= dw2+dend*ymvvq*ypiq*vviq*(-3*xvvi) dw2= dw2+dend*ymvvq*ypic*vviq*(9*xvv+6*xvvi) dw2= dw2+dend*ymvvq*ypis*vviq*(-9*xvq-18*xvv-3*xvvi) dw2= dw2+dend*ymvvq*ypis*vvis*(-1./2.*xvq) dw2= dw2+dend*ymvvq*ypi*vviq*(3*xvc+18*xvq+9*xvv) dw2= dw2+dend*ymvvq*ypi*vvis*(xvc+xvq) dw2= dw2+dend*ymvvq*vviq*omyp*(-3*xvc) dw2= dw2+dend*ymvvq*vviq*(-3*xvc-9*xvq) dw2= dw2+dend*ymvvq*vvis*omyp*(-xvc) dw2= dw2+dend*ymvvq*vvis*omyps*(-1./2.*xvs) dw2= dw2+dend*ymvvq*vvis*(-xvc-1./2.*xvq) * if(ii.eq.3) then wcom= wcom*(2.d0-ym)**2/yms-dw1+2.d0/ym*(1.d0-1.d0/ym)*dw2 else if(ii.eq.4) then wcom= 0.5d0*em4*omxm*omxm*(dw1-dw2) else if(ii.eq.5) then wcom= em2*((2.d0+0.5d0*omxm)*dw1+omxm*(2.d0/ym- # 3.d0/2.d0)*dw2) endif * endif * www= wcom/pwtm2 * 1 if(iz.eq.0) then ww(ii)= 0.d0 iz= 1 else ww(ii)= www*tjac*tfactsw endif * enddo * 2 if(iz.eq.0) then wwt= 0.d0 else wwt= ww(1)+ww(2)+ww(3)+ww(4)+ww(5) endif * if(wwt.lt.0.d0) then iz= 0 ifz(7)= ifz(7)+1 wtoxssw= 0.d0 else wtoxssw= wwt endif * return end * *------------------------------------------------------------------------- * real*8 function wtoxsswc(ndim,x) implicit real*8(a-h,o-z) character*1 opeak * common/wtps/opeak common/wtkount/ik common/wtorti/indz common/wtistrf/isf common/wtipt/ifz(51) common/wtfsw/tfactsw common/wtnpr/ipr,ipr0 common/wticuts/iac(4) common/wtcutsw/tcs,s0cut common/wtfmass/em,rmm,tm,rnm,uqm,dqm,cqm,sqm,bqm,tqm,dmy common/wtort/zca,zcb,zcc,z0,beta1s,ztmp,zl,zu,zln,zun,sqz common/wtbpar/wm,zm,zg,gf,pi,pis,cfct,fcnt,ge,alphai,alwi common/wtcpar/alpha,hbet,hbeti,omhb,eob,d0gl,g8,tfact,pih,alw, # eta,feta,beta,g2,tfacth common/wtapar/ars,s,rwm,rwm2,rwg,rwmg,swg,swgs,opswgs,sth2,cth2, # hsth2,tsth2,scth2,asth2,tth2,rzm,rzm2,rzg,rzmg,szg, # szgs,opszgs,sth4,cth4,ve,vf,vfp,rbqm2,rszw,rszw2, # s0w,s0z common/wtcuts/aim(6),bim(6),ae(4),asa(4),bsa(4),afsa(6),bfsa(6), # ombsa(4),opbsa(4),teq,rae(4),omasa(4),opasa(4), # sg12,cg12,sg13,cg13,sg14,cg14,sg23,cg23,sg24, # cg24,sg34,cg34,sct120,sct130,sct140,sct230, # sct240,sct340,sgam(4),cgam(4) * dimension vu(4),vut(4) dimension vd(4),vdt(4) dimension vn(4),vnt(4) dimension x(ndim),ww(3,5),wwt(3),tjac(3,5),tjac0(3,5) * external c05adf,s09aaf * *-----order of integration is: uv,vv; y; x2,x1,tau/ups,z,phi_u/phi_d * ik= ik+1 * if(ipr.eq.4) then ipro= 1 else if(ipr.eq.5) then ipro= 2 endif do ii=1,3 wwt(ii)= 0.d0 do jj=1,5 ww(ii,jj)= 0.d0 tjac(ii,jj)= 0.d0 tjac0(ii,jj)= 0.d0 enddo enddo * em2= em*em em4= em2*em2 em6= em4*em2 s2= s*s s3= s2*s rem= em/sqrt(s) rwms= rwm2 rwmq= rwms*rwms rwgs= rwg*rwg rwgm= rwg*rwm rwgms= rwgs*rwms rwmgs= rwmg*rwmg rwmpgs= rwms+rwgs if(ipro.eq.2) then ums= uqm*uqm/s dms= dqm*dqm/s else if(ipro.eq.1) then ums= 0.d0 dms= rmm*rmm/s endif um= sqrt(ums) dm= sqrt(dms) s0= s0cut * if(ndim.eq.7) then yx= x(1) x2x= x(2) x1x= x(3) xx= x(4) taux= x(5) zx= x(6) phx= x(7) else if(ndim.eq.9) then uvx= x(1) vvx= x(2) yx= x(3) x2x= x(4) x1x= x(5) xx= x(6) taux= x(7) zx= x(8) phx= x(9) endif * iz= 1 * if(ndim.eq.7) then ueps= 0.d0 uv= 1.d0 uvs= uv*uv ujc= 1.d0 veps= 0.d0 vv= 1.d0 vjc= 1.d0 else if(uvx.gt.1.d0) then iz= 0 ifz(1)= ifz(1)+1 go to 1 endif vcut= (4.d0*s0+tcs)/(4.d0+tcs) vvl= 0.25d0*(rem+sqrt(rem*rem+4.d0*vcut))**2 omul= 1.d0-vvl ueps= omul*(1.d0-uvx)**hbeti uv= 1.d0-ueps ujc= omul**hbet uvs= uv*uv * uvl= uv-vvl veps= uvl*(1.d0-vvx)**hbeti vv= uv-veps vjc0= 1.d0-vvl/uv if(vjc0.le.0.d0) then iz= 0 ifz(2)= ifz(2)+1 go to 1 else vjc= vjc0**hbet endif endif * if(vv.lt.0.d0) then iz= 0 ifz(3)= ifz(3)+1 go to 1 endif svv= sqrt(vv) vvs= vv*vv vvc= vvs*vv vvq= vvc*vv vvi= 1.d0/vv rdms= dms*vv rums= ums*vv xm= uv xms= xm*xm xp= vv/uv xps= xp*xp xmop= xm/xp sh= vv*s ssh= sqrt(sh) * opxp= 1.d0+xp opxm= 1.d0+xm omxp= veps/uv omxm= ueps if(isf.eq.0) then stfp= 1.d0 stfm= 1.d0 else if(isf.gt.0) then if(omxp.eq.0) then stfp= d0gl else rcpx= 0.25d0*opxp*opxp rcpy= xp iflp= 1 rclp= s21baf(rcpx,rcpy,iflp) stfp= d0gl+eob*omxp**omhb*(-0.5d0*opxp+ # feta*(-4.d0*opxp*log(omxp)+ # 3.d0*opxp*log(xp)+4.d0*rclp-5.d0-xp)) endif if(omxm.eq.0) then stfm= d0gl else rcmx= 0.25d0*opxm*opxm rcmy= xm iflm= 1 rclm= s21baf(rcmx,rcmy,iflm) stfm= d0gl+eob*omxm**omhb*(-0.5d0*opxm+ # feta*(-4.d0*opxm*log(omxm)+ # 3.d0*opxm*log(xm)+4.d0*rclm-5.d0-xm)) endif endif * stf= stfp*stfm * *-----y is initialized * yl= (4.d0*s0+tcs)/(4.d0+tcs)/vv yu= 1.d0-em/ssh if(yl.gt.yu) then iz= 0 ifz(4)= ifz(4)+1 go to 1 endif yjc= yu-yl y= yjc*yx+yl ys= y*y yc= ys*y yq= yc*y yvv= y*vv yvvs= yvv*yvv omy= (1.d0-s0/vv)*(1.d0-yx)+em/ssh*yx omys= omy*omy omyc= omys*omy dmsh= dms/yvv umsh= ums/yvv dmh= sqrt(dmsh) umh= sqrt(umsh) umshs= umsh*umsh dmshs= dmsh*dmsh * *-----x2 is initialized * x2l1= (um+dm)*(um+dm)/yvv x2l2= s0/yvv x2l= dmax1(x2l1,x2l2) x2u= 1.d0 if(x2l.gt.x2u) then iz= 0 ifz(5)= ifz(5)+1 go to 1 endif vwmg= rwmg*yvv zma= yvv*x2l zmb= yvv*x2u zmas= zma-rwm2 zmbs= zmb-rwm2 atma= zmas/rwmg atmb= zmbs/rwmg * if(opeak.eq.'y') then if(atma.gt.1.d0.and.atmb.gt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm x2jc0= (-atmb+atma)/rwgm else if(atma.gt.1.d0.and.atmb.lt.-1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm x2jc0= (-pi+atmb+atma)/rwgm else if(atma.gt.1.d0.and.abs(atmb).lt.1.d0) then atma= 1.d0/atma atma= atan(atma) zmat= (pih-atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm x2jc0= (-pih+atmb+atma)/rwgm else if(atma.lt.-1.d0.and.atmb.gt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm x2jc0= (pi-atmb-atma)/rwgm else if(atma.lt.-1.d0.and.atmb.lt.-1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm x2jc0= (atmb-atma)/rwgm else if(atma.lt.-1.d0.and.abs(atmb).lt.1.d0) then atma= -1.d0/atma atma= atan(atma) zmat= (-pih+atma)/rwgm atmb= atan(atmb) zmbt= atmb/rwgm x2jc0= (pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.gt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= 1.d0/atmb atmb= atan(atmb) zmbt= (pih-atmb)/rwgm x2jc0= (pih-atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.atmb.lt.-1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= -1.d0/atmb atmb= atan(atmb) zmbt= (-pih+atmb)/rwgm x2jc0= (-pih+atmb-atma)/rwgm else if(abs(atma).lt.1.d0.and.abs(atmb).lt.1.d0) then atma= atan(atma) zmat= atma/rwgm atmb= atan(atmb) zmbt= atmb/rwgm x2jc0= (atmb-atma)/rwgm endif * x2jc= x2jc0/yvv argx2= rwgm*(x2jc0*x2x+zmat) iftn= 1 s07= s07aaf(argx2,iftn) x2= (rwms+rwgm*s07)/yvv pmjc= 1.d0 * else if(opeak.eq.'n') then smjc0= zmb-zma x2= (smjc0*x2x+zma)/yvv pmjc= 1.d0/((yvv*x2-rwm2)**2+rwmgs) x2jc= smjc0 endif omx2= 1.d0-x2 omx2s= omx2*omx2 x2s= x2*x2 x2c= x2s*x2 x2q= x2c*x2 pwtm2= ((1.d0+xv-x2+z)*yvv+rwms)**2+rwms*rwgs * *-----x1 is initialized * xlk= x2s-2.d0*(dmsh+umsh)*x2+(dmsh-umsh)*(dmsh-umsh) xlk= sqrt(xlk) x1l1= dmsh+2.d0*umh x1l2= x2 x1l3= dmsh+(umsh-dmsh+(1.d0+umsh-dmsh)*x2+x2s-omx2*xlk)/x2 x1l= dmax1(x1l1,x1l2,x1l3) x1u1= 1.d0 x1u2= (1.d0-dmh)**2+x2 x1u3= 1.d0-dmsh*omx2/x2 x1u4= dmsh+(umsh-dmsh+(1.d0+umsh-dmsh)*x2+x2s+omx2*xlk)/x2 x1u= dmin1(x1u1,x1u2,x1u3,x1u4) if(x1l.gt.x1u) then iz= 0 ifz(6)= ifz(6)+1 go to 1 endif x1jc= x1u-x1l x1= x1jc*x1x+x1l omx1= 1.d0-x1 omx1s= omx1*omx1 x1mo= x1-1.d0 x1mos= x1mo*x1mo x1mx2= x1-x2 x2mx1= x2-x1 x2mx1s= x2mx1*x2mx1 x1s= x1*x1 x1c= x1s*x1 x1q= x1c*x1 x1f= x1q*x1 x12= x1*x2 x12s= x1*x2s x12c= x1*x2c x12q= x1*x2q x1s2s= x1s*x2s edm= 1.d0+x2-x1 edms= edm*edm edmc= edms*edm edmq= edmc*edm edmf= edmq*edm edmsx= edmf*edm reu= x1-dmsh red= 1.d0+x2-x1+dmsh betus= 1.d0-4.d0*umsh/reu betds= 1.d0-4.d0*dmsh/red betu= sqrt(betus) betd= sqrt(betds) * *-----X is initialized * qh0s= em2*(xm*y+omxm)*y/omy qhcs= qh0s+0.25d0*xm*xm*omy*tcs*s * x01= qh0s/y/sh x02= -1.d0+dms/x1/yvv x0= dmax1(x01,x02) xc0= qhcs/y/sh xc1= 1.d0-s0/yvv xc= dmin1(xc0,xc1) if(xc.eq.xc0) then dxc= 0.25d0*xm*xm*omy*tcs/yvv else dxc= xc-x0 endif * if(x0.gt.xc) then iz= 0 ifz(7)= ifz(7)+1 go to 1 endif * *-----first loop starts here * do ii=1,3 * if(ii.eq.1) then xjac= dxc/xc/x0/y/sh xv= x0*xc/(dxc*xx+x0) else if(ii.eq.2) then xjac= log(xc/x0) xv= exp(xjac*xx+log(x0)) else if(ii.eq.3) then xjac= dxc xv= xjac*xx+x0 endif omxv= 1.d0-xv opxv= 1.d0+xv opxvs= opxv*opxv xvs= xv*xv * if(ii.eq.1) then psf= em2*(2.d0+omxm*(2.d0/y-1.d0)) else if(ii.eq.2) then psf= -1.d0+2.d0/y*(1.d0-1.d0/y) else if(ii.eq.3) then psf= 1.d0 endif * *-----second loop starts here * do jj=1,5 * if(jj.le.3) then psjc= 1.d0/reu/betu/opxv else psjc= 1.d0/red/betd/opxv endif * *-----tau/upsilon is initialized * *-----Here is the tau branch * if(jj.le.3) then if((umsh/reu).lt.1.d-5) then taul= -x1-xv+umsh+dmsh+umsh*opxv/reu tauu= -xv*(1.d0+x1+dmsh)+2.d0*dmsh+umsh*(1.d0-opxv/reu) else taul= -x1-xv+umsh+dmsh+0.5d0*(1.d0-betu)*opxv*reu tauu= -x1-xv+umsh+dmsh+0.5d0*(1.d0+betu)*opxv*reu endif phta= umsh*dmsh*(-1+x2) phta= phta+umsh*(1-x1*x2-x1+x2s) phta= phta+umshs*(1) phta= phta+dmsh*(-x1*x2+x1-x2+x2s) phta= phta-x2*x1mo*x2mx1 * phtb= xv*umsh*(2-2*x1*x2-x1*x2s-3*x1+x1s*x2+x1s+2*x2s) phtb= phtb+xv*dmsh*(-2*x1*x2-2*x1*x2s+2*x1+2*x1s*x2- # x1s-2*x2+3*x2s) phtb= phtb+xv*(x1*x2*x1mo*x2mx1-2*x2*x1mo*x2mx1) phtb= phtb+umsh*(2*x1*x2+3*x1*x2s+x1-3*x1s*x2-x1s-2*x2s) phtb= phtb+dmsh*(2*x1*x2s-2*x1s*x2+x1s-x2s) phtb= phtb-x1*x2*x1mo*x2mx1 * phtc= xv*umsh*(2*x1*x2s+x1-2*x1s*x2-2*x1s*x2s- # 2*x1s+2*x1c*x2+x1c) phtc= phtc+xv*dmsh*(4*x1*x2s-3*x1s*x2-3*x1s*x2s+ # x1s+3*x1c*x2-x1c-x2s) phtc= phtc+xv*(-x1*x2*x1mo*x2mx1+x1s*x2*x1mo*x2mx1) phtc= phtc+xvs*umsh*(1-2*x1*x2-2*x1*x2s-2*x1+2*x1s* # x2+x1s+2*x2s) phtc= phtc+xvs*dmsh*(-x1*x2-2*x1*x2s+x1+2*x1s*x2- # x1s-x2+2*x2s) phtc= phtc+xvs*(x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1) phtc= phtc+umsh*(-x1*x2-2*x1*x2s+2*x1s*x2+x1s* # x2s-x1c*x2+x2s) * phtdisc= umsh*opxvs*(8*x1*x2c+8*x1*x2q-2*x1s*x2s- # 16*x1s*x2c-2*x1s*x2q-2*x1c*x2+4*x1c* # x2s+4*x1c*x2c-2*x1c*x2q+4*x1q*x2+4*x1q* # x2c-2*x1f*x2-2*x1f*x2s-4*x2q) phtdisc= phtdisc+dmsh*opxvs*(-2*x1*x2q+2*x1s*x2c+6*x1s* # x2q+2*x1c*x2s-10*x1c*x2c-4*x1c*x2q-2* # x1q*x2+2*x1q*x2s+8*x1q*x2c+2*x1f*x2-4*x1f*x2s) phtdisc= phtdisc+opxvs*(x2s*x1mos*x2mx1s*x1s) if(phtb.gt.0.d0) then phtq= -0.5d0*(phtb+phtdisc) else if(phtb.lt.0.d0) then phtq= -0.5d0*(phtb-phtdisc) endif ataul= phtq/phta atauu= phtc/phtq taupl= dmin1(ataul,atauu) taupu= dmax1(ataul,atauu) taul= dmax1(taul,taupl) tauu= dmin1(tauu,taupu) if(taul.gt.tauu) then iz= 0 ifz(8)= ifz(8)+1 go to 3 endif arglnl= umsh-taul arglnu= umsh-tauu * *-----Here is the ups branch * else if(jj.ge.4) then upsl= -1.d0-x2+x1-xv+0.5d0*(1.d0-betd)*opxv*red upsu= -1.d0-x2+x1-xv+0.5d0*(1.d0+betd)*opxv*red phta= umsh*dmsh*(-1+x2) phta= phta+umsh*(1-x1*x2-x1+x2s) phta= phta+umshs*(1) phta= phta+dmsh*(-x1*x2+x1-x2+x2s) phta= phta-x2*x1mo*x2mx1 * phtb= xv*umsh*(1+2*x1*x2s-x1s*x2-x1s-x2+x2s-x2c) phtb= phtb+xv*dmsh*(-2*x1*x2+3*x1*x2s+x1-2*x1s*x2+ # x1s-x2+x2s-x2c) phtb= phtb+xv*(-x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1+ # x2s*x1mo*x2mx1) phtb= phtb+umsh*(1-2*x1*x2-2*x1*x2s-2*x1+x1s*x2+ # x1s+x2+x2s+x2c) phtb= phtb+dmsh*(2*x1*x2-x1*x2s+x1-x1s-x2-x2s+x2c) phtb= phtb+x1*x2*x1mo*x2mx1-x2*x1mo*x2mx1-x2s*x1mo*x2mx1 * phtc= xv*umsh*(3*x1*x2+3*x1*x2s+3*x1*x2c+x1-3* # x1s*x2-3*x1s*x2s-2*x1s+x1c*x2+x1c-x2- # x2s-x2c-x2q) phtc= phtc+xv*dmsh*(-2*x1*x2+x1*x2s+4*x1*x2c+x1s* # x2-5*x1s*x2s+x1s+2*x1c*x2-x1c+x2s-x2c-x2q) phtc= phtc+xv*(-x1*x2*x1mo*x2mx1-2*x1*x2s*x1mo* # x2mx1+x1s*x2*x1mo*x2mx1+x2s*x1mo*x2mx1+ # x2c*x1mo*x2mx1) phtc= phtc+xvs*umsh*(x1*x2+2*x1*x2s+x1-x1s*x2- # x1s-x2-x2c) phtc= phtc+xvs*dmsh*(-2*x1*x2+2*x1*x2s-x1s*x2+ # x1s+x2s-x2c) phtc= phtc+xvs*(-x1*x2*x1mo*x2mx1+x2s*x1mo*x2mx1) phtc= phtc+dmsh*(2*x1*x2s+x1*x2c-x1s*x2-2*x1s* # x2s+x1c*x2-x2c) * phtdisc= umsh*(-24*x1*xv*edms+60*x1*xv*edmc- # 56*x1*xv*edmq+24*x1*xv*edmf-4*x1*xv* # edmsx-12*x1*xvs*edms+30*x1*xvs*edmc # -28*x1*xvs*edmq+12*x1*xvs*edmf-2*x1* # xvs*edmsx-12*x1*edms+30*x1*edmc-28* # x1*edmq+12*x1*edmf-2*x1*edmsx+24*x1s # *xv*edms-48*x1s*xv*edmc+32*x1s*xv* # edmq-8*x1s*xv*edmf+12*x1s*xvs* # edms-24*x1s*xvs*edmc+16*x1s*xvs* # edmq-4*x1s*xvs*edmf+12*x1s*edms # -24*x1s*edmc+16*x1s*edmq-4*x1s* # edmf-8*x1c*xv*edms+12*x1c*xv*edmc # -4*x1c*xv*edmq-4*x1c*xvs*edms+6*x1c # *xvs*edmc-2*x1c*xvs*edmq-4*x1c* # edms+6*x1c*edmc-2*x1c*edmq+8*xv*edms) phtdisc= phtdisc+umsh*(-24*xv*edmc+28*xv*edmq-16*xv* # edmf+4*xv*edmsx+4*xvs*edms-12*xvs* # edmc+14*xvs*edmq-8*xvs*edmf+2* # xvs*edmsx+4*edms-12*edmc+14* # edmq-8*edmf+2*edmsx) phtdisc= phtdisc+dmsh*(-4-128*x1*xv*edm+180*x1*xv*edms # -92*x1*xv*edmc-8*x1*xv*edmq+20*x1*xv* # edmf-4*x1*xv*edmsx+32*x1*xv-64*x1*xvs* # edm+90*x1*xvs*edms-46*x1*xvs*edmc-4 # *x1*xvs*edmq+10*x1*xvs*edmf-2*x1*xvs* # edmsx+16*x1*xvs-64*x1*edm+90*x1*edms # -46*x1*edmc-4*x1*edmq+10*x1*edmf-2* # x1*edmsx+16*x1+168*x1s*xv*edm-208*x1s*xv* # edms+100*x1s*xv*edmc-8*x1s*xv*edmq # -4*x1s*xv*edmf-48*x1s*xv+84*x1s*xvs* # edm-104*x1s*xvs*edms+50*x1s*xvs* # edmc-4*x1s*xvs*edmq-2*x1s*xvs* # edmf-24*x1s*xvs+84*x1s*edm-104*x1s* # edms+50*x1s*edmc-4*x1s*edmq-2* # x1s*edmf-24*x1s-96*x1c*xv*edm+100*x1c # *xv*edms) phtdisc= phtdisc+dmsh*(-40*x1c*xv*edmc+4*x1c*xv* # edmq+32*x1c*xv-48*x1c*xvs*edm+50*x1c* # xvs*edms-20*x1c*xvs*edmc+2*x1c*xvs* # edmq+16*x1c*xvs-48*x1c*edm+50*x1c* # edms-20*x1c*edmc+2*x1c*edmq+16* # x1c+20*x1q*xv*edm-16*x1q*xv*edms+4* # x1q*xv*edmc-8*x1q*xv+10*x1q*xvs*edm-8 # *x1q*xvs*edms+2*x1q*xvs*edmc-4*x1q* # xvs+10*x1q*edm-8*x1q*edms+2*x1q* # edmc-4*x1q+36*xv*edm-56*xv*edms+28* # xv*edmc+12*xv*edmq-16*xv*edmf+4*xv* # edmsx-8*xv+18*xvs*edm-28*xvs*edms+ # 14*xvs*edmc+6*xvs*edmq-8*xvs*edmf # +2*xvs*edmsx-4*xvs+18*edm-28*edms # +14*edmc+6*edmq-8*edmf+2*edmsx) phtdisc= phtdisc+2*xv*edms*x1mos*x2s-4*xv*edmc* # x1mos*x2s+2*xv*edmq*x1mos*x2s+xvs* # edms*x1mos*x2s-2*xvs*edmc*x1mos* # x2s+xvs*edmq*x1mos*x2s+edms* # x1mos*x2s-2*edmc*x1mos*x2s+edmq*x1mos*x2s * if(phtb.gt.0.d0) then phtq= -0.5d0*(phtb+phtdisc) else if(phtb.lt.0.d0) then phtq= -0.5d0*(phtb-phtdisc) endif aupsl= phtq/phta aupsu= phtc/phtq upspl= dmin1(aupsl,aupsu) upspu= dmax1(aupsl,aupsu) upsl= dmax1(upsl,upspl) upsu= dmin1(upsu,upspu) if(upsl.gt.upsu) then iz= 0 ifz(8)= ifz(8)+1 go to 3 endif arglnl= dmsh-upsl arglnu= dmsh-upsu endif * *-----Quark propagators are mapped away (for no angular cuts) * if(asa(3).eq.-1.d0.and.asa(4).eq.-1.d0.and. # bsa(3).eq.1.d0.and.bsa(4).eq.1.d0) then if(jj.eq.1) then taujc= tauu-taul tau= taujc*taux+taul taus= tau*tau optau= 1.d0+tau else if(jj.eq.2) then taujc= log(arglnu/arglnl) tau= umsh-exp((taujc*taux+log(arglnl))) taus= tau*tau optau= 1.d0+tau else if(jj.eq.3) then taujc= 1.d0/arglnu-1.d0/arglnl tau= umsh-1.d0/(taujc*taux+1.d0/arglnl) taus= tau*tau optau= 1.d0+tau else if(jj.eq.4) then upsjc= log(arglnu/arglnl) ups= dmsh-exp((upsjc*taux+log(arglnl))) upss= ups*ups opups= 1.d0+ups else if(jj.eq.5) then upsjc= 1.d0/arglnu-1.d0/arglnl ups= dmsh-1.d0/(upsjc*taux+1.d0/arglnl) upss= ups*ups opups= 1.d0+ups endif else if(jj.le.3) then taujc= tauu-taul tau= taujc*taux+taul taus= tau*tau optau= 1.d0+tau else if(jj.ge.4) then upsjc= upsu-upsl ups= upsjc*taux+upsl upss= ups*ups opups= 1.d0+ups endif if(jj.eq.2) then taujc= taujc/(tau-umsh) else if(jj.eq.3) then taujc= taujc/(tau-umsh)/(tau-umsh) else if(jj.eq.4) then upsjc= upsjc/(ups-dmsh) else if(jj.eq.5) then upsjc= upsjc/(ups-dmsh)/(ups-dmsh) endif endif * *-----z is initialized * if(omx2.eq.0.d0.or.x1.eq.0.d0) then iz= 0 ifz(9)= ifz(9)+1 go to 3 else if(jj.le.3) then ct= (-omxv-2.d0*(tau+xv-umsh)/reu)/opxv/betu cp= (1.d0-2.d0*(x1-x2-umsh)/omx2/reu)/betu else if(jj.ge.4) then ct= (-omxv-2.d0*(ups+xv-dmsh)/red)/opxv/betd cp= (1.d0-2.d0*(omx1+umsh)/omx2/red)/betd endif if(ct.lt.-1.d0.or.ct.gt.1.d0) then iz= 0 ifz(10)= ifz(10)+1 go to 3 else st= sqrt(1.d0-ct*ct) endif if(cp.lt.-1.d0.or.cp.gt.1.d0) then iz= 0 ifz(11)= ifz(11)+1 go to 3 else sp= sqrt(1.d0-cp*cp) endif endif * zca= -1.d0 zcb= -xv-0.5d0*omx2*(omxv+opxv*ct*cp) zcc= -xvs+xv*omx2s-xv*omxv*omx2-0.25d0*omx2s*opxvs* # (ct*ct+cp*cp)-xv*opxv*omx2*ct*cp-0.5d0*omxv*opxv* # omx2s*ct*cp zdisc= 0.5d0*opxv*omx2*st*sp if(zcb.gt.0.d0) then zcq= -(zcb+zdisc) else if(zcb.lt.0.d0) then zcq= -(zcb-zdisc) endif azl= -zcq azu= zcc/zcq zl= dmin1(azl,azu) zu= dmax1(azl,azu) if(jj.le.3) then zl1= -tau-x1+x2-xv+umsh+dmsh-0.5d0*(1.d0+betd)*opxv*red zu1= -tau-x1+x2-xv+umsh+dmsh-0.5d0*(1.d0-betd)*opxv*red zl2= -1.d0+x2-xv zu2= -xv*x2 * phca= umsh*(16) phca= phca+dmsh*(8*x1) phca= phca+dmshs*(-4) phca= phca-4*x1s * phcb= xv*umsh*(32-8*x1+16*x2) phcb= phcb+xv*dmsh*(-8+8*x1+16*x2) phcb= phcb+xv*(-16*x1*x2-8*x1+16*x2) phcb= phcb+umsh*(16+8*x1*x2+16*x1-32*x2+16*tau) phcb= phcb+dmsh*(8*x1+8*x2*tau-8*x2-8*tau) phcb= phcb-8*x1*x2*tau+8*x1*x2-8*x1*tau-8*x1s+16*x2*tau * phcc= xv*umsh*(16+32*x1*x2+8*x1+8*x2*tau-32*x2-16*x2s+8*tau) phcc= phcc+xv*dmsh*(16*x1*x2-8*x1+8*x2*tau+8*x2-16*x2s-8*tau) phcc= phcc+xv*(-16*x1*x2*tau+16*x1*x2+16*x1*x2s-8*x1- # 16*x1s*x2+24*x2*tau+8*x2-16*x2s-8*tau) phcc= phcc+xvs*umsh*(8+16*x2) phcc= phcc+xvs*dmsh*(-8+16*x2) phcc= phcc+xvs*(-4-16*x1*x2+16*x2) phcc= phcc+umsh*(-8*x1*x2+16*x1-24*x2*tau-16*x2+8* # x2s*tau+8*x2s+16*tau) phcc= phcc+8*x1*x2*tau+8*x1*x2-8*x1*tau-4*x1s+8*x2*tau+ # 8*x2*taus-8*x2s*tau-4*x2s*taus-4*x2s-4*taus * phct= umsh*dmsh*(-1+x2) phct= phct+umsh*(1-x1*x2-x1+x2s) phct= phct+umshs*(1) phct= phct+dmsh*(-x1*x2+x1-x2+x2s) phct= phct-x2*x1mo*x2mx1 phdisc= 256.d0*phct*(tau-taupl)*(tau-taupu) if(phcb.gt.0.d0) then phcq= -0.5d0*(phcb+sqrt(phdisc)) else if(phcb.lt.0.d0) then phcq= -0.5d0*(phcb-sqrt(phdisc)) endif azl3= phcq/phca azu3= phcc/phcq zl3= dmin1(azl3,azu3) zu3= dmax1(azl3,azu3) zln= dmax1(zl,zl1,zl2,zl3) zun= dmin1(zu,zu1,zu2,zu3) else if(jj.ge.4) then if((umsh/reu).lt.1.d-5) then zl1= -1.d0-ups+opxv*dmsh-xv*x1+umsh*opxv/reu zu1= -1.d0+x1-ups-xv-umsh*opxv/reu else zl1= -1.d0+x1-ups-xv-0.5d0*(1.d0+betu)*opxv*reu zu1= -1.d0+x1-ups-xv-0.5d0*(1.d0-betu)*opxv*reu endif zl2= -1.d0+x2-xv zu2= -xv*x2 phca= dmsh*(8+8*x1-8*x2) phca= phca+dmshs*(-4) phca= phca-4+8*x1*x2+8*x1-4*x1s-8*x2-4*x2s * phcb= xv*umsh*(-8-8*x1+8*x2) phcb= phcb+xv*dmsh*(16+8*x1-8*x2) phcb= phcb+xv*(-8+16*x1*x2+8*x1-8*x2-16*x2s) phcb= phcb+umsh*(-8+8*x1-8*x2-16*ups) phcb= phcb+dmsh*(16-8*x1*x2-8*x2*ups-16*x2+8*x2s+8*ups) phcb= phcb-8+8*x1*x2*ups+8*x1*x2+8*x1*ups+16*x1-8*x1s- # 8*x2-8*x2s*ups-8*ups * phcc= xv*umsh*(-16+16*x1*x2+8*x1-8*x2*ups-16*x2s-8*ups) phcc= phcc+xv*dmsh*(16-8*x1-8*x2*ups-8*x2+8*ups) phcc= phcc+xv*(-8+16*x1*x2*ups+16*x1*x2+16*x1*x2s+8* # x1-16*x1s*x2+8*x2*ups-16*x2s*ups-16*x2s-8*ups) phcc= phcc+xvs*umsh*(-8) phcc= phcc+xvs*dmsh*(8) phcc= phcc+xvs*(-4+16*x1*x2-16*x2s) phcc= phcc+umsh*(-8+8*x1+8*x2*ups-8*ups) phcc= phcc+dmsh*(8+8*x1*x2-8*x1-16*x2*ups-8*x2+8*x2s* # ups+8*ups) phcc= phcc-4-8*x1*x2*ups+8*x1*ups+8*x1-4*x1s+8*x2*ups+ # 8*x2*upss-4*x2s*upss-8*ups-4*upss * phdisc= 256.d0*phta*(ups-upspl)*(ups-upspu) if(phcb.gt.0.d0) then phcq= -0.5d0*(phcb+sqrt(phdisc)) else if(phcb.lt.0.d0) then phcq= -0.5d0*(phcb-sqrt(phdisc)) endif azl3= phcq/phca azu3= phcc/phcq zl3= dmin1(azl3,azu3) zu3= dmax1(azl3,azu3) zln= dmax1(zl,zl1,zl2,zl3) zun= dmin1(zu,zu1,zu2,zu3) endif * if(zln.gt.zun) then iz= 0 ifz(12)= ifz(12)+1 go to 3 endif * *-----the overall sqrt is mapped away * sqz= 0.5d0*opxv*omx2*st*sp if(zl.eq.zln) then tl= -pi/2.d0 else ifas= 1 aa= (zcb-zln)/sqz tl= -s09aaf(aa,ifas) if(ifas.ne.0) then iz= 0 ifz(13)= ifz(13)+1 go to 3 endif endif if(zu.eq.zun) then tu= pi/2.d0 else ifas= 1 bb= (zcb-zun)/sqz tu= -s09aaf(bb,ifas) if(ifas.ne.0) then iz= 0 ifz(14)= ifz(14)+1 go to 3 endif endif * arg= (tu-tl)*zx+tl z= zcb+sqz*sin(arg) zjc= tu-tl zs= z*z zc= zs*z * if(st.ne.0.d0.and.sp.ne.0.d0) then acphi= z+xv+0.5d0*omx2*(opxv*ct*cp+omxv) bcphi= 0.5d0*omx2*opxv*st*sp if((acphi+bcphi).lt.0.d0.or.(acphi-bcphi).lt.0.d0) then iz= 0 ifz(15)= ifz(15)+1 go to 3 endif endif * *-----phi is initialized * if(jj.le.3) then cpu= cos(2.d0*pi*phx) spu= sin(2.d0*pi*phx) else if(jj.ge.4) then cpd= cos(2.d0*pi*phx) spd= sin(2.d0*pi*phx) endif * *-----here one starts with u,d variables and then n variables * if(jj.le.3) then eu= x1-dmsh ctu= (1.d0-2.d0*(tau+x1+xv-umsh-dmsh)/opxv/reu)/betu if(abs(ctu).gt.1.d0) then iz= 0 ifz(16)= ifz(16)+1 go to 3 endif epsu= 1.d0-abs(ctu) if(epsu.lt.1.d-5) then stu= sqrt(2.d0*epsu)*(1.d0-epsu/4.d0-epsu*epsu/32.d0) else stu= sqrt(1.d0-ctu*ctu) endif if(stu.eq.0.d0) then iz= 0 ifz(17)= ifz(17)+1 go to 3 endif * ed= 1.d0+x2-x1+dmsh ctd= (1.d0+2.d0*(tau+x1-x2+z+xv-umsh-dmsh)/opxv/ # red)/betd if(abs(ctd).gt.1.d0) then iz= 0 ifz(18)= ifz(18)+1 go to 3 endif epsd= 1.d0-abs(ctd) if(epsd.lt.1.d-5) then std= sqrt(2.d0*epsd)*(1.d0-epsd/4.d0-epsd*epsd/32.d0) else std= sqrt(1.d0-ctd*ctd) endif aca= betu*betd*stu*std*cpu acb= betu*betd*stu*std*spu acca= -1.d0+2.d0*(x2-umsh-dmsh)/reu/red acc= betu*betd*ctu*ctd+acca discaz= sqrt(phca)/opxv/reu/red*sqrt((z-zl3)*(z-zu3)) byp= (-acb+discaz)/(acc-aca) bym= (-acb-discaz)/(acc-aca) cpdp= (1.d0-byp*byp)/(1.d0+byp*byp) cpdm= (1.d0-bym*bym)/(1.d0+bym*bym) spdp= 2.d0*byp/(1.d0+byp*byp) spdm= 2.d0*bym/(1.d0+bym*bym) * en= omx2 ctn= 1.d0-2.d0*(1.d0-x2+z+xv)/opxv/omx2 if(abs(ctn).gt.1.d0) then iz= 0 ifz(19)= ifz(19)+1 go to 3 endif epsn= 1.d0-abs(ctn) if(epsn.lt.1.d-5) then stn= sqrt(2.d0*epsn)*(1.d0-epsn/4.d0-epsn*epsn/32.d0) else stn= sqrt(1.d0-ctn*ctn) endif cpnp= -(eu*stu*cpu+ed*std*cpdp)/en/stn cpnm= -(eu*stu*cpu+ed*std*cpdm)/en/stn spnp= -(eu*stu*spu+ed*std*spdp)/en/stn spnm= -(eu*stu*spu+ed*std*spdm)/en/stn tstnp= 1.d0-cpnp*cpnp-spnp*spnp tstnm= 1.d0-cpnm*cpnm-spnm*spnm * *-----here one starts with d,u variables and then n variables * else if(jj.ge.4) then ed= 1.d0+x2-x1+dmsh ctd= (1.d0-2.d0*(1.d0+x2-x1+ups+xv)/opxv/ # red)/betd if(abs(ctd).gt.1.d0) then iz= 0 ifz(16)= ifz(16)+1 go to 3 endif epsd= 1.d0-abs(ctd) if(epsd.lt.1.d-5) then std= sqrt(2.d0*epsd)*(1.d0-epsd/4.d0-epsd*epsd/32.d0) else std= sqrt(1.d0-ctd*ctd) endif if(std.eq.0.d0) then iz= 0 ifz(17)= ifz(17)+1 go to 3 endif * eu= x1-dmsh ctu= (1.d0+2.d0*(omx1+ups+z+xv)/opxv/reu)/betu if(abs(ctu).gt.1.d0) then iz= 0 ifz(18)= ifz(18)+1 go to 3 endif epsu= 1.d0-abs(ctu) if(epsu.lt.1.d-5) then stu= sqrt(2.d0*epsu)*(1.d0-epsu/4.d0-epsu*epsu/32.d0) else stu= sqrt(1.d0-ctu*ctu) endif aca= betu*betd*stu*std*cpd acb= betu*betd*stu*std*spd acca= -1.d0+2.d0*(x2-umsh-dmsh)/reu/red acc= betu*betd*ctu*ctd+acca discaz= sqrt(phca)/opxv/reu/red*sqrt((z-zl3)*(z-zu3)) byp= (-acb+sqrt(discaz))/(acc-aca) bym= (-acb-sqrt(discaz))/(acc-aca) cpup= (1.d0-byp*byp)/(1.d0+byp*byp) cpum= (1.d0-bym*bym)/(1.d0+bym*bym) spup= 2.d0*byp/(1.d0+byp*byp) spum= 2.d0*bym/(1.d0+bym*bym) * en= omx2 ctn= 1.d0-2.d0*(1.d0-x2+z+xv)/opxv/omx2 if(abs(ctn).gt.1.d0) then iz= 0 ifz(19)= ifz(19)+1 go to 3 endif epsn= 1.d0-abs(ctn) if(epsn.lt.1.d-5) then stn= sqrt(2.d0*epsn)*(1.d0-epsn/4.d0-epsn*epsn/32.d0) else stn= sqrt(1.d0-ctn*ctn) endif cpnp= -(eu*stu*cpup+ed*std*cpd)/en/stn cpnm= -(eu*stu*cpum+ed*std*cpd)/en/stn spnp= -(eu*stu*spup+ed*std*spd)/en/stn spnm= -(eu*stu*spum+ed*std*spd)/en/stn tstnp= 1.d0-cpnp*cpnp-spnp*spnp tstnm= 1.d0-cpnm*cpnm-spnm*spnm endif * if(jj.le.3) then if(abs(cpdp).le.1.d0.and.abs(cpnp).le.1.d0.and. # abs(spnp).le.1.d0.and.abs(tstnp).le.1.d-10) then cpd= cpdp spd= spdp cpn= cpnp spn= spnp else if(abs(cpdm).le.1.d0.and.abs(cpnm).le.1.d0.and. # abs(spnm).le.1.d0.and.abs(tstnm).le.1.d-10) then cpd= cpdm spd= spdm cpn= cpnm spn= spnm else iz= 0 ifz(20)= ifz(20)+1 go to 3 endif else if(jj.ge.4) then if(abs(cpup).le.1.d0.and.abs(cpnp).le.1.d0.and. # abs(spnp).le.1.d0.and.abs(tstnp).le.1.d-10) then cpu= cpup spu= spup cpn= cpnp spn= spnp else if(abs(cpum).le.1.d0.and.abs(cpnm).le.1.d0.and. # abs(spnm).le.1.d0.and.abs(tstnm).le.1.d-10) then cpu= cpum spu= spum cpn= cpnm spn= spnm else iz= 0 ifz(20)= ifz(20)+1 go to 3 endif endif * vfact= xp+xm*y vfacts= vfact*vfact velp= 1.d0-2.d0*xm*y/vfact+xv*xm*y/vfact*(-6.d0+2.d0*y+ # 6.d0*xm*y/vfact-2.d0*xm*ys/vfact)+xvs*xm*y/vfact*( # -8.d0+10.d0*y-2.d0*ys+26.d0*xm*y/vfact-22.d0*xm*ys/vfact+ # 4.d0*xm*yc/vfact-18.d0*xms*ys/vfacts+12.d0*xms*yc/vfacts- # 2.d0*xms*yq/vfacts) if(abs(velp).gt.1.d0) then iz= 0 ifz(21)= ifz(21)+1 go to 3 endif betpl= 2.d0*sqrt(xp*xm*y)/vfact*(1.d0+ # xv*(3.d0/2.d0-3.d0*xm*y/vfact+xm*ys/vfact-0.5d0*y)+ # xvs*(7.d0/8.d0-17.d0/2.d0*xm*y/vfact+8.d0*xm*ys/vfact- # 3.d0/2.d0*xm*yc/vfact+9.d0*xms*ys/vfacts-6.d0*xms*yc/ # vfacts+xms*yq/vfacts-7.d0/4.d0*y+3.d0/8.d0*ys)) * if(abs(betpl).gt.1.d0) then iz= 0 ifz(22)= ifz(22)+1 go to 3 endif betl= 1.d0-0.5d0*xv*omy+xvs/8.d0*(-1.d0+6.d0*y-ys) if(abs(betl).gt.1.d0) then iz= 0 ifz(23)= ifz(23)+1 go to 3 endif vels= xv*(omy-xv*y) if(vels.lt.0.d0) then vel= sqrt(xv*omy)*(1.d0-0.5d0*xv*y/omy) else vel= sqrt(vels) endif if(abs(vel).gt.1.d0) then iz= 0 ifz(24)= ifz(24)+1 go to 3 endif * betli= 1.d0/betl betpli= 1.d0/betpl * ysn= sqrt(y*sh) * vu(1)= 0.5d0*eu*stu*cpu*ysn vu(2)= 0.5d0*eu*stu*spu*ysn vu(3)= 0.5d0*eu*ctu*ysn vu(4)= 0.5d0*eu*ysn * vd(1)= 0.5d0*ed*std*cpd*ysn vd(2)= 0.5d0*ed*std*spd*ysn vd(3)= 0.5d0*ed*ctd*ysn vd(4)= 0.5d0*ed*ysn * vn(1)= 0.5d0*en*stn*cpn*ysn vn(2)= 0.5d0*en*stn*spn*ysn vn(3)= 0.5d0*en*ctn*ysn vn(4)= 0.5d0*en*ysn * vut(1)= vel*(vu(4)-vu(3))-vu(1) vut(2)= vu(2) vut(3)= (vel*(1.d0+velp)*vu(1)+(vel*vel-1.d0)*vu(3)- # (vel*vel+velp)*vu(4))/betl/betpl vut(4)= (-vel*(1.d0+velp)*vu(1)+velp*(1.d0-vel*vel)*vu(3)+ # (1.d0+vel*vel*velp)*vu(4))/betl/betpl * vdt(1)= vel*(vd(4)-vd(3))-vd(1) vdt(2)= vd(2) vdt(3)= (vel*(1.d0+velp)*vd(1)+(vel*vel-1.d0)*vd(3)- # (vel*vel+velp)*vd(4))/betl/betpl vdt(4)= (-vel*(1.d0+velp)*vd(1)+velp*(1.d0-vel*vel)*vd(3)+ # (1.d0+vel*vel*velp)*vd(4))/betl/betpl * vnt(1)= vel*(vn(4)-vn(3))-vn(1) vnt(2)= vn(2) vnt(3)= (vel*(1.d0+velp)*vn(1)+(vel*vel-1.d0)*vn(3)- # (vel*vel+velp)*vn(4))/betl/betpl vnt(4)= (-vel*(1.d0+velp)*vn(1)+velp*(1.d0-vel*vel)*vn(3)+ # (1.d0+vel*vel*velp)*vn(4))/betl/betpl * eut= vut(4) ctut= vut(3)/eut edt= vdt(4) ctdt= vdt(3)/edt ent= vnt(4) ctnt= vnt(3)/ent * * if(eut.lt.ae(3).and.edt.lt.ae(4)) then * iz= 0 * ifz(25)= ifz(25)+1 * go to 3 * endif * if(eut.lt.ae(3)) then * iz= 0 * ifz(26)= ifz(26)+1 * go to 3 * endif * if(edt.lt.ae(4)) then * iz= 0 * ifz(27)= ifz(27)+1 * go to 3 * endif * if(ctut.lt.asa(3).or.ctut.gt.bsa(3)) then * iz= 0 * ifz(28)= ifz(28)+1 * go to 3 * endif * if(ctdt.lt.asa(4).or.ctdt.gt.bsa(4)) then * iz= 0 * ifz(29)= ifz(29)+1 * go to 3 * endif * if(jj.le.3) then tjac0(ii,jj)= ujc*vjc*yjc*x2jc*x1jc*xjac*taujc*zjc* # psjc*pmjc else if(jj.ge.4) then tjac0(ii,jj)= ujc*vjc*yjc*x2jc*x1jc*xjac*upsjc*zjc* # psjc*pmjc endif tjac(ii,jj)= tjac0(ii,jj)*stf*psf/sh * if(ipro.eq.2.and.ii.le.2) then if(jj.eq.1) then * wcom= y*rwms*rwmpgs*(-vvs*tau*x1+1./3.*vvs* # tau*x2+2./3.*vvs*tau+2./3.*vvs*x1*x2+1./3.*vvs*x1 # *z+5./3.*vvs*x1-2*vvs*x1s-1./3.*vvs*x2*z-1./9. # *vvs*x2-2./9.*vvs*z-2./9.*vvs) wcom= wcom+ys*rwms*(4*vvc*tau*x1*x2- # vvc*tau*x1*z-3*vvc*tau*x1+1./3.*vvc*tau*x2*z-2./3. # *vvc*tau*x2-2./3.*vvc*tau*x2s+1./3.*vvc*tau*z+1. # /3.*vvc*tau+vvc*taus*x2-vvc*taus-4*vvc*x1* # x2-4./3.*vvc*x1*x2s+4./3.*vvc*x1*z+1./3.*vvc*x1 # *zs+3*vvc*x1+6*vvc*x1s*x2-2*vvc*x1s*z-4* # vvc*x1s+1./9.*vvc*x2*z-1./3.*vvc*x2*zs+13./9. # *vvc*x2+2./3.*vvc*x2s*z+2./9.*vvc*x2s-2./3.* # vvc*z-1./3.*vvc*zs-4./3.*vvc) wcom= wcom+yc*(3*vvq*tau*x1*x2*z+3*vvq # *tau*x1*x2-3*vvq*tau*x1*x2s-1./3.*vvq*tau*x2*z-1./ # 3.*vvq*tau*x2-1./3.*vvq*tau*x2s*z+1./3.*vvq*tau* # x2c+vvq*taus*x2*z+vvq*taus*x2-vvq*taus* # x2s-10./3.*vvq*x1*x2*z-1./3.*vvq*x1*x2*zs-3* # vvq*x1*x2-1./3.*vvq*x1*x2s*z+7./3.*vvq*x1*x2s # +2./3.*vvq*x1*x2c+4*vvq*x1s*x2*z+4*vvq*x1s* # x2-4*vvq*x1s*x2s+5./3.*vvq*x2*z+1./3.*vvq*x2* # zs+4./3.*vvq*x2+1./9.*vvq*x2s*z+1./3.*vvq* # x2s*zs-11./9.*vvq*x2s-1./3.*vvq*x2c*z-1./9. # *vvq*x2c) * else if(jj.eq.2) then * wcom= y*rwms*rwmpgs*(-4./9.*vvs*x1*x2-14./9.* # vvs*x1+4./9.*vvs*x1s*x2+20./9.*vvs*x1s-4./3.* # vvs*x1c+2./9.*vvs*x2+4./9.*vvs) wcom= wcom+ys*rwms*(16./9.*vvc*x1*x2+8. # /9.*vvc*x1*x2s-2./3.*vvc*x1*z-2./3.*vvc*x1-28./ # 9.*vvc*x1s*x2-8./9.*vvc*x1s*x2s+4./3.*vvc* # x1s*z+4./3.*vvc*x1s+8./3.*vvc*x1c*x2-4./3.* # vvc*x1c*z-4./3.*vvc*x1c-2./9.*vvc*x2-4./9.* # vvc*x2s) wcom= wcom+yc*(2./3.*vvq*x1*x2*z+2./3.* # vvq*x1*x2-2./9.*vvq*x1*x2s-4./9.*vvq*x1*x2c-4. # /3.*vvq*x1s*x2*z-4./3.*vvq*x1s*x2+8./9.*vvq* # x1s*x2s+4./9.*vvq*x1s*x2c+4./3.*vvq*x1c*x2* # z+4./3.*vvq*x1c*x2-4./3.*vvq*x1c*x2s-2./9.* # vvq*x2s+2./9.*vvq*x2c) * else if(jj.eq.3) then * wcom= y*rums*rwms*(-32./9.*vvs*x1*x2+16./9.* # vvs*x1s*x2+16./9.*vvs*x2) wcom= wcom+ys*rums*(16./9.*vvc*x1*x2s # -8./9.*vvc*x1s*x2s-8./9.*vvc*x2s) wcom= wcom+rums*rwms*rwmpgs*(16./9.*vv*x1 # -8./9.*vv*x1s-8./9.*vv) * else if(jj.eq.4) then * wcom= y*rwms*rwmpgs*(2./3.*vvs*x1*x2*z+4./9.* # vvs*x1*x2+10./9.*vvs*x1*z+1./3.*vvs*x1*zs+14./9. # *vvs*x1-2./9.*vvs*x1s*x2-2./3.*vvs*x1s*z-16./ # 9.*vvs*x1s+2./3.*vvs*x1c-4./9.*vvs*x2*z-1./3. # *vvs*x2*zs-2./9.*vvs*x2-4./9.*vvs*z-2./9.*vvs # *zs-4./9.*vvs) wcom= wcom+ys*rwms*(-8./9.*vvc*x1*x2*z # -28./9.*vvc*x1*x2-4./3.*vvc*x1*x2s*z-8./9.*vvc* # x1*x2s+10./3.*vvc*x1*z+5./3.*vvc*x1*zs+1./3.* # vvc*x1*zc+2*vvc*x1+2./3.*vvc*x1s*x2*z+32./9.* # vvc*x1s*x2+4./9.*vvc*x1s*x2s-8./3.*vvc*x1s* # z-2./3.*vvc*x1s*zs-2*vvc*x1s-4./3.*vvc* # x1c*x2+2./3.*vvc*x1c*z+2./3.*vvc*x1c+2./9.* # vvc*x2*z-2./9.*vvc*x2*zs-1./3.*vvc*x2*zc+8./9. # *vvc*x2+8./9.*vvc*x2s*z+2./3.*vvc*x2s*zs+4./ # 9.*vvc*x2s-4./3.*vvc*z-vvc*zs-1./3.*vvc* # zc-2./3.*vvc) wcom= wcom+yc*(-10./3.*vvq*x1*x2*z-5./ # 3.*vvq*x1*x2*zs-1./3.*vvq*x1*x2*zc-2*vvq*x1*x2 # -2./9.*vvq*x1*x2s*z-1./3.*vvq*x1*x2s*zs+14./9. # *vvq*x1*x2s+2./3.*vvq*x1*x2c*z+4./9.*vvq*x1* # x2c+8./3.*vvq*x1s*x2*z+2./3.*vvq*x1s*x2*zs+ # 2*vvq*x1s*x2-16./9.*vvq*x1s*x2s-2./9.*vvq* # x1s*x2c-2./3.*vvq*x1c*x2*z-2./3.*vvq*x1c*x2 # +2./3.*vvq*x1c*x2s+4./3.*vvq*x2*z+vvq*x2*zs # +1./3.*vvq*x2*zc+2./3.*vvq*x2+2./9.*vvq*x2s*z # +4./9.*vvq*x2s*zs+1./3.*vvq*x2s*zc-4./9.* # vvq*x2s-4./9.*vvq*x2c*z-1./3.*vvq*x2c*zs- # 2./9.*vvq*x2c) * else if(jj.eq.5) then * wcom= y*rdms*rwms*(-8./9.*vvs*x1*x2*z-8./9.* # vvs*x1*x2+4./9.*vvs*x1s*x2+8./9.*vvs*x2*z+4./9. # *vvs*x2*zs+4./9.*vvs*x2) wcom= wcom+ys*rdms*(4./9.*vvc*x1*x2s*z # +4./9.*vvc*x1*x2s-2./9.*vvc*x1s*x2s-4./9.* # vvc*x2s*z-2./9.*vvc*x2s*zs-2./9.*vvc*x2s) wcom= wcom+rdms*rwms*rwmpgs*(4./9.*vv*x1*z # +4./9.*vv*x1-2./9.*vv*x1s-4./9.*vv*z-2./9.*vv*zs # -2./9.*vv) * endif else if(ipro.eq.2.and.ii.eq.3) then if(jj.eq.1) then * wcom= y*rwms*rwmpgs*(-vvs*tau*x1+1./3.*vvs* # tau*x2+2./3.*vvs*tau+2./3.*vvs*x1*x2+1./3.*vvs*x1 # *z+5./3.*vvs*x1-2*vvs*x1s-1./3.*vvs*x2*z-1./9. # *vvs*x2-2./9.*vvs*z-2./9.*vvs) wcom= wcom+ys*rwms*(4*vvc*tau*x1*x2- # vvc*tau*x1*z-3*vvc*tau*x1+1./3.*vvc*tau*x2*z-2./3. # *vvc*tau*x2-2./3.*vvc*tau*x2s+1./3.*vvc*tau*z+1. # /3.*vvc*tau+vvc*taus*x2-vvc*taus-4*vvc*x1* # x2-4./3.*vvc*x1*x2s+4./3.*vvc*x1*z+1./3.*vvc*x1 # *zs+3*vvc*x1+6*vvc*x1s*x2-2*vvc*x1s*z-4* # vvc*x1s+1./9.*vvc*x2*z-1./3.*vvc*x2*zs+13./9. # *vvc*x2+2./3.*vvc*x2s*z+2./9.*vvc*x2s-2./3.* # vvc*z-1./3.*vvc*zs-4./3.*vvc) wcom= wcom+yc*(3*vvq*tau*x1*x2*z+3*vvq # *tau*x1*x2-3*vvq*tau*x1*x2s-1./3.*vvq*tau*x2*z-1./ # 3.*vvq*tau*x2-1./3.*vvq*tau*x2s*z+1./3.*vvq*tau* # x2c+vvq*taus*x2*z+vvq*taus*x2-vvq*taus* # x2s-10./3.*vvq*x1*x2*z-1./3.*vvq*x1*x2*zs-3* # vvq*x1*x2-1./3.*vvq*x1*x2s*z+7./3.*vvq*x1*x2s # +2./3.*vvq*x1*x2c+4*vvq*x1s*x2*z+4*vvq*x1s* # x2-4*vvq*x1s*x2s+5./3.*vvq*x2*z+1./3.*vvq*x2* # zs+4./3.*vvq*x2+1./9.*vvq*x2s*z+1./3.*vvq* # x2s*zs-11./9.*vvq*x2s-1./3.*vvq*x2c*z-1./9. # *vvq*x2c) * dw1= y*rwms*rwmpgs*(8./3.-4./3.*vvs*x1+ # vvs*x1s-2./3.*vvs*x2+7./9.*vvs+2*tau*x1-2./3. # *tau*x2-4./3.*tau-8./3.*x1*x2-2./3.*x1*z-8*x1+6* # x1s+2./3.*x2*z+16./9.*x2+2./9.*x2s+4./9.*z) dw1= dw1+ys*rwms*(-12*vv*tau*x1*x2+2 # *vv*tau*x1*z+8*vv*tau*x1-2./3.*vv*tau*x2*z+5*vv*tau*x2 # +8./3.*vv*tau*x2s-2./3.*vv*tau*z-11./3.*vv*tau-2*vv # *taus*x2+2*vv*taus+4./3.*vv*x1*x2*z+23*vv*x1*x2+8 # *vv*x1*x2s-13./3.*vv*x1*z-2./3.*vv*x1*zs-10*vv*x1 # -22*vv*x1s*x2+4*vv*x1s*z+10*vv*x1s+1./9.*vv*x2* # z+2./3.*vv*x2*zs-6*vv*x2-8./3.*vv*x2s*z-44./9.*vv # *x2s-4./9.*vv*x2c+5./3.*vv*z+2./3.*vv*zs+7./3.* # vv+1./3.*vvc*tau*x2-2./3.*vvc*tau+4./3.*vvc*x1*x2 # -vvc*x1*z-1./3.*vvc*x1-2*vvc*x1s*x2+2*vvc* # x1s*z-vvc*x2*z+1./9.*vvc*x2+4./3.*vvc*x2s-2. # /3.*vvc*z-vvc) dw1= dw1+yc*(-6*vvs*tau*x1*x2*z-15* # vvs*tau*x1*x2+14*vvs*tau*x1*x2s+vvs*tau*x1*z+5./ # 2.*vvs*tau*x1+2./3.*vvs*tau*x2*z+14./3.*vvs*tau*x2 # +2./3.*vvs*tau*x2s*z-11./3.*vvs*tau*x2s-2*vvs # *tau*x2c-1./2.*vvs*tau-2*vvs*taus*x2*z-5*vvs* # taus*x2+4*vvs*taus*x2s+vvs*taus+13./3.* # vvs*x1*x2*z+2./3.*vvs*x1*x2*zs+11*vvs*x1*x2-2./ # 3.*vvs*x1*x2s*z-15*vvs*x1*x2s-16./3.*vvs*x1* # x2c-1./2.*vvs*x1-6*vvs*x1s*x2*z-14*vvs*x1s* # x2+18*vvs*x1s*x2s+vvs*x1s*z+3./2.*vvs*x1s # -5./3.*vvs*x2*z-2./3.*vvs*x2*zs-7./3.*vvs*x2- # 5./9.*vvs*x2s*z-2./3.*vvs*x2s*zs+10./3.*vvs* # x2s+2*vvs*x2c*z+28./9.*vvs*x2c+2./9.*vvs* # x2q+2*vvq*tau*x1*x2+vvq*tau*x1*z+1./2.*vvq*tau* # x1+2./3.*vvq*tau*x2-1./3.*vvq*tau*x2s-1./2.*vvq # *tau) dw1= dw1+yc*(vvq*taus*x2-vvq*x1*x2 # *z-5./3.*vvq*x1*x2-vvq*x1*z-1./2.*vvq*x1-2* # vvq*x1s*x2*z+2*vvq*x1s*x2+vvq*x1s*x2s+2* # vvq*x1s*z+vvq*x1s*zs+1./2.*vvq*x1s+8./3.* # vvq*x2*z+2*vvq*x2+vvq*x2s*z-8./9.*vvq*x2s # -2./3.*vvq*x2c) dw2= y*rwms*rwmpgs*(8-4./3.*vvs*x1+vvs* # x1s-2./3.*vvs*x2+7./9.*vvs+6*tau*x1-2*tau*x2- # 4*tau-8*x1*x2-2*x1*z-24*x1+18*x1s+2*x2*z+16./3. # *x2+2./3.*x2s+4./3.*z) dw2= dw2+ys*rwms*(-36*vv*tau*x1*x2+6 # *vv*tau*x1*z+24*vv*tau*x1-2*vv*tau*x2*z+15*vv*tau*x2+ # 8*vv*tau*x2s-2*vv*tau*z-11*vv*tau-6*vv*taus*x2+6* # vv*taus+4*vv*x1*x2*z+69*vv*x1*x2+24*vv*x1*x2s-13* # vv*x1*z-2*vv*x1*zs-30*vv*x1-66*vv*x1s*x2+12*vv* # x1s*z+30*vv*x1s+1./3.*vv*x2*z+2*vv*x2*zs-18*vv* # x2-8*vv*x2s*z-44./3.*vv*x2s-4./3.*vv*x2c+5*vv*z # +2*vv*zs+7*vv+1./3.*vvc*tau*x2-2./3.*vvc*tau+ # 4./3.*vvc*x1*x2-vvc*x1*z-1./3.*vvc*x1-2*vvc* # x1s*x2+2*vvc*x1s*z-vvc*x2*z+1./9.*vvc*x2+4./ # 3.*vvc*x2s-2./3.*vvc*z-vvc) dw2= dw2+yc*(-18*vvs*tau*x1*x2*z-45* # vvs*tau*x1*x2+42*vvs*tau*x1*x2s+3*vvs*tau*x1*z+ # 15./2.*vvs*tau*x1+2*vvs*tau*x2*z+14*vvs*tau*x2+2* # vvs*tau*x2s*z-11*vvs*tau*x2s-6*vvs*tau*x2c- # 3./2.*vvs*tau-6*vvs*taus*x2*z-15*vvs*taus*x2+ # 12*vvs*taus*x2s+3*vvs*taus+13*vvs*x1*x2*z+2 # *vvs*x1*x2*zs+33*vvs*x1*x2-2*vvs*x1*x2s*z-45* # vvs*x1*x2s-16*vvs*x1*x2c-3./2.*vvs*x1-18* # vvs*x1s*x2*z-42*vvs*x1s*x2+54*vvs*x1s*x2s # +3*vvs*x1s*z+9./2.*vvs*x1s-5*vvs*x2*z-2* # vvs*x2*zs-7*vvs*x2-5./3.*vvs*x2s*z-2*vvs* # x2s*zs+10*vvs*x2s+6*vvs*x2c*z+28./3.*vvs* # x2c+2./3.*vvs*x2q+2*vvq*tau*x1*x2+vvq*tau*x1* # z+1./2.*vvq*tau*x1+2./3.*vvq*tau*x2-1./3.*vvq*tau # *x2s-1./2.*vvq*tau+vvq*taus*x2-vvq*x1*x2*z- # 5./3.*vvq*x1*x2) dw2= dw2+yc*(-vvq*x1*z-1./2.*vvq* # x1-2*vvq*x1s*x2*z+2*vvq*x1s*x2+vvq*x1s* # x2s+2*vvq*x1s*z+vvq*x1s*zs+1./2.*vvq* # x1s+8./3.*vvq*x2*z+2*vvq*x2+vvq*x2s*z-8./9. # *vvq*x2s-2./3.*vvq*x2c) * else if(jj.eq.2) then * wcom= y*rwms*rwmpgs*(-4./9.*vvs*x1*x2-14./9.* # vvs*x1+4./9.*vvs*x1s*x2+20./9.*vvs*x1s-4./3.* # vvs*x1c+2./9.*vvs*x2+4./9.*vvs) wcom= wcom+ys*rwms*(16./9.*vvc*x1*x2+8. # /9.*vvc*x1*x2s-2./3.*vvc*x1*z-2./3.*vvc*x1-28./ # 9.*vvc*x1s*x2-8./9.*vvc*x1s*x2s+4./3.*vvc* # x1s*z+4./3.*vvc*x1s+8./3.*vvc*x1c*x2-4./3.* # vvc*x1c*z-4./3.*vvc*x1c-2./9.*vvc*x2-4./9.* # vvc*x2s) wcom= wcom+yc*(2./3.*vvq*x1*x2*z+2./3.* # vvq*x1*x2-2./9.*vvq*x1*x2s-4./9.*vvq*x1*x2c-4. # /3.*vvq*x1s*x2*z-4./3.*vvq*x1s*x2+8./9.*vvq* # x1s*x2s+4./9.*vvq*x1s*x2c+4./3.*vvq*x1c*x2* # z+4./3.*vvq*x1c*x2-4./3.*vvq*x1c*x2s-2./9.* # vvq*x2s+2./9.*vvq*x2c) * dw1= y*rwms*rwmpgs*(-8./9.+4./9.*vvs*x1*x2 # -8./9.*vvs*x1-8./9.*vvs*x1s*x2+8./9.*vvs*x1s # +2./9.*vvs-4./9.*x1*x2+40./9.*x1+4./3.*x1s*x2- # 52./9.*x1s-8./9.*x1c*x2+4./3.*x1c+8./9.*x1q) dw1= dw1+ys*rwms*(-92./9.*vv*x1*x2+8. # /9.*vv*x1*x2s+4./3.*vv*x1+140./9.*vv*x1s*x2-8./3.* # vv*x1s*x2s-4./3.*vv*x1s*z-4*vv*x1s-16./3.*vv* # x1c*x2+16./9.*vv*x1c*x2s+4./3.*vv*x1c*z+8./3.* # vv*x1c-16./9.*vv*x1q*x2+16./9.*vv*x2+4./9.*vvc*x1 # *x2-8./9.*vvc*x1*x2s-16./9.*vvc*x1s*x2+16./9.* # vvc*x1s*x2s+4./3.*vvc*x1s-4./3.*vvc*x1c+8. # /9.*vvc*x2-2./3.*vvc) dw1= dw1+yc*(-4./3.*vvs*x1*x2+52./9. # *vvs*x1*x2s-4./9.*vvs*x1*x2c+4./3.*vvs*x1s*x2 # *z+4*vvs*x1s*x2-88./9.*vvs*x1s*x2s+4./3.* # vvs*x1s*x2c-4./3.*vvs*x1c*x2*z-8./3.*vvs* # x1c*x2+4*vvs*x1c*x2s-8./9.*vvs*x1c*x2c+8./ # 9.*vvs*x1q*x2s-8./9.*vvs*x2s+4./9.*vvq*x1* # x2s+4./9.*vvq*x1*x2c-4./3.*vvq*x1s*x2+8./9.* # vvq*x1s*x2s-8./9.*vvq*x1s*x2c+4./3.*vvq* # x1c*x2+2./3.*vvq*x2-10./9.*vvq*x2s) dw2= y*rwms*rwmpgs*(-8./3.+4./9.*vvs*x1*x2 # -8./9.*vvs*x1-8./9.*vvs*x1s*x2+8./9.*vvs*x1s # +2./9.*vvs-4./3.*x1*x2+40./3.*x1+4*x1s*x2-52./3. # *x1s-8./3.*x1c*x2+4*x1c+8./3.*x1q) dw2= dw2+ys*rwms*(-92./3.*vv*x1*x2+8. # /3.*vv*x1*x2s+4*vv*x1+140./3.*vv*x1s*x2-8*vv*x1s* # x2s-4*vv*x1s*z-12*vv*x1s-16*vv*x1c*x2+16./3.* # vv*x1c*x2s+4*vv*x1c*z+8*vv*x1c-16./3.*vv*x1q* # x2+16./3.*vv*x2+4./9.*vvc*x1*x2-8./9.*vvc*x1*x2s # -16./9.*vvc*x1s*x2+16./9.*vvc*x1s*x2s+4./3.* # vvc*x1s-4./3.*vvc*x1c+8./9.*vvc*x2-2./3.* # vvc) dw2= dw2+yc*(-4*vvs*x1*x2+52./3.* # vvs*x1*x2s-4./3.*vvs*x1*x2c+4*vvs*x1s*x2*z+ # 12*vvs*x1s*x2-88./3.*vvs*x1s*x2s+4*vvs*x1s* # x2c-4*vvs*x1c*x2*z-8*vvs*x1c*x2+12*vvs* # x1c*x2s-8./3.*vvs*x1c*x2c+8./3.*vvs*x1q* # x2s-8./3.*vvs*x2s+4./9.*vvq*x1*x2s+4./9.* # vvq*x1*x2c-4./3.*vvq*x1s*x2+8./9.*vvq*x1s* # x2s-8./9.*vvq*x1s*x2c+4./3.*vvq*x1c*x2+2./3. # *vvq*x2-10./9.*vvq*x2s) * else if(jj.eq.3) then * wcom= y*rums*rwms*(-32./9.*vvs*x1*x2+16./9.* # vvs*x1s*x2+16./9.*vvs*x2) wcom= wcom+ys*rums*(16./9.*vvc*x1*x2s # -8./9.*vvc*x1s*x2s-8./9.*vvc*x2s) wcom= wcom+rums*rwms*rwmpgs*(16./9.*vv*x1 # -8./9.*vv*x1s-8./9.*vv) * dw1= y*rums*rwms*(8./9.*vvs*x1*x2-8./9.* # vvs*x2+16./9.*x1*x2-32./9.*x1s*x2+16./9.*x1c*x2) dw1= dw1+y*rwms*rwmpgs*(-8./9.*vvs* # x1+4./9.*vvs*x1s+4./9.*vvs-8./9.*x1+8./3.*x1s # -8./3.*x1c+8./9.*x1q) dw1= dw1+ys*rums*(-8./9.*vv*x1*x2s # +16./9.*vv*x1s*x2s-8./9.*vv*x1c*x2s-4./9.*vvc # *x1*x2s+4./9.*vvc*x2s) dw1= dw1+ys*rwms*(16./9.*vv*x1*x2-16. # /3.*vv*x1s*x2+16./3.*vv*x1c*x2-16./9.*vv*x1q*x2+ # 16./9.*vvc*x1*x2-8./9.*vvc*x1s*x2-8./9.*vvc*x2) dw1= dw1+yc*(-8./9.*vvs*x1*x2s+8. # /3.*vvs*x1s*x2s-8./3.*vvs*x1c*x2s+8./9.*vvs # *x1q*x2s-8./9.*vvq*x1*x2s+4./9.*vvq*x1s*x2s # +4./9.*vvq*x2s) dw1= dw1+rums*rwms*rwmpgs*(-8./9.* # vvi*x1+16./9.*vvi*x1s-8./9.*vvi*x1c- # 4./9.*vv*x1+4./9.*vv) dw2= y*rums*rwms*(8./9.*vvs*x1*x2-8./9.* # vvs*x2+16./3.*x1*x2-32./3.*x1s*x2+16./3.*x1c*x2) dw2= dw2+y*rwms*rwmpgs*(-8./9.*vvs* # x1+4./9.*vvs*x1s+4./9.*vvs-8./3.*x1+8*x1s-8 # *x1c+8./3.*x1q) dw2= dw2+ys*rums*(-8./3.*vv*x1*x2s # +16./3.*vv*x1s*x2s-8./3.*vv*x1c*x2s-4./9.*vvc # *x1*x2s+4./9.*vvc*x2s) dw2= dw2+ys*rwms*(16./3.*vv*x1*x2-16 # *vv*x1s*x2+16*vv*x1c*x2-16./3.*vv*x1q*x2+16./9.* # vvc*x1*x2-8./9.*vvc*x1s*x2-8./9.*vvc*x2) dw2= dw2+yc*(-8./3.*vvs*x1*x2s+8 # *vvs*x1s*x2s-8*vvs*x1c*x2s+8./3.*vvs*x1q* # x2s-8./9.*vvq*x1*x2s+4./9.*vvq*x1s*x2s+4./9. # *vvq*x2s) dw2= dw2+rums*rwms*rwmpgs*(-8./3.* # vvi*x1+16./3.*vvi*x1s-8./3.*vvi*x1c- # 4./9.*vv*x1+4./9.*vv) * else if(jj.eq.4) then * wcom= y*rwms*rwmpgs*(2./3.*vvs*x1*x2*z+4./9.* # vvs*x1*x2+10./9.*vvs*x1*z+1./3.*vvs*x1*zs+14./9. # *vvs*x1-2./9.*vvs*x1s*x2-2./3.*vvs*x1s*z-16./ # 9.*vvs*x1s+2./3.*vvs*x1c-4./9.*vvs*x2*z-1./3. # *vvs*x2*zs-2./9.*vvs*x2-4./9.*vvs*z-2./9.*vvs # *zs-4./9.*vvs) wcom= wcom+ys*rwms*(-8./9.*vvc*x1*x2*z # -28./9.*vvc*x1*x2-4./3.*vvc*x1*x2s*z-8./9.*vvc* # x1*x2s+10./3.*vvc*x1*z+5./3.*vvc*x1*zs+1./3.* # vvc*x1*zc+2*vvc*x1+2./3.*vvc*x1s*x2*z+32./9.* # vvc*x1s*x2+4./9.*vvc*x1s*x2s-8./3.*vvc*x1s* # z-2./3.*vvc*x1s*zs-2*vvc*x1s-4./3.*vvc* # x1c*x2+2./3.*vvc*x1c*z+2./3.*vvc*x1c+2./9.* # vvc*x2*z-2./9.*vvc*x2*zs-1./3.*vvc*x2*zc+8./9. # *vvc*x2+8./9.*vvc*x2s*z+2./3.*vvc*x2s*zs+4./ # 9.*vvc*x2s-4./3.*vvc*z-vvc*zs-1./3.*vvc* # zc-2./3.*vvc) wcom= wcom+yc*(-10./3.*vvq*x1*x2*z-5./ # 3.*vvq*x1*x2*zs-1./3.*vvq*x1*x2*zc-2*vvq*x1*x2 # -2./9.*vvq*x1*x2s*z-1./3.*vvq*x1*x2s*zs+14./9. # *vvq*x1*x2s+2./3.*vvq*x1*x2c*z+4./9.*vvq*x1* # x2c+8./3.*vvq*x1s*x2*z+2./3.*vvq*x1s*x2*zs+ # 2*vvq*x1s*x2-16./9.*vvq*x1s*x2s-2./9.*vvq* # x1s*x2c-2./3.*vvq*x1c*x2*z-2./3.*vvq*x1c*x2 # +2./3.*vvq*x1c*x2s+4./3.*vvq*x2*z+vvq*x2*zs # +1./3.*vvq*x2*zc+2./3.*vvq*x2+2./9.*vvq*x2s*z # +4./9.*vvq*x2s*zs+1./3.*vvq*x2s*zc-4./9.* # vvq*x2s-4./9.*vvq*x2c*z-1./3.*vvq*x2c*zs- # 2./9.*vvq*x2c) * dw1= y*rwms*rwmpgs*(8./9.+22./9.*vvs*x1*x2 # +2./3.*vvs*x1*z+2./9.*vvs*x1-8./9.*vvs*x1s*x2 # -4./9.*vvs*x1s-4./3.*vvs*x2*z-14./9.*vvs*x2+2. # /9.*vvs-28./9.*x1*x2*z-56./9.*x1*x2-4./9.*x1*x2s- # 4*x1*z-2./3.*x1*zs-14./3.*x1+14./3.*x1s*x2+8./3.* # x1s*z+68./9.*x1s-8./9.*x1c*x2-14./3.*x1c+8./9. # *x1q+8./3.*x2*z+2./3.*x2*zs+22./9.*x2+4./9.*x2s # *z+4./9.*x2s+4./3.*z+4./9.*zs) dw1= dw1+ys*rwms*(25./3.*vv*x1*x2*z+4./ # 3.*vv*x1*x2*zs+28./3.*vv*x1*x2+80./9.*vv*x1*x2s*z+ # 136./9.*vv*x1*x2s+8./9.*vv*x1*x2c-19./3.*vv*x1*z-11. # /3.*vv*x1*zs-2./3.*vv*x1*zc-10./3.*vv*x1-22./3.*vv* # x1s*x2*z-154./9.*vv*x1s*x2-32./3.*vv*x1s*x2s+13. # /3.*vv*x1s*z+4./3.*vv*x1s*zs+11./3.*vv*x1s+32./ # 3.*vv*x1c*x2+16./9.*vv*x1c*x2s-2./3.*vv*x1c*z-4. # /3.*vv*x1c-16./9.*vv*x1q*x2-vv*x2*z+7./9.*vv*x2* # zs+2./3.*vv*x2*zc-10./9.*vv*x2-8*vv*x2s*z-8./3. # *vv*x2s*zs-56./9.*vv*x2s-8./9.*vv*x2c*z-8./9.* # vv*x2c+8./3.*vv*z+7./3.*vv*zs+2./3.*vv*zc+vv+ # 4./3.*vvc*x1*x2*z+8./9.*vvc*x1*x2-44./9.*vvc*x1* # x2s+5*vvc*x1*z+vvc*x1*zs+14./3.*vvc*x1+2./9. # *vvc*x1s*x2+16./9.*vvc*x1s*x2s-2*vvc*x1s*z # -10./3.*vvc*x1s+2./3.*vvc*x1c-8./3.*vvc*x2*z # -5./3.*vvc*x2*zs) dw1= dw1+ys*rwms*(-10./9.*vvc*x2+8. # /3.*vvc*x2s*z+28./9.*vvc*x2s-3*vvc*z-4./3.* # vvc*zs-2*vvc) dw1= dw1+yc*(19./3.*vvs*x1*x2*z+11./3. # *vvs*x1*x2*zs+2./3.*vvs*x1*x2*zc+10./3.*vvs*x1* # x2-13./3.*vvs*x1*x2s*z-2./3.*vvs*x1*x2s*zs-14. # /3.*vvs*x1*x2s-52./9.*vvs*x1*x2c*z-80./9.*vvs* # x1*x2c-4./9.*vvs*x1*x2q-13./3.*vvs*x1s*x2*z-4. # /3.*vvs*x1s*x2*zs-11./3.*vvs*x1s*x2+14./3.* # vvs*x1s*x2s*z+86./9.*vvs*x1s*x2s+6*vvs* # x1s*x2c+2./3.*vvs*x1c*x2*z+4./3.*vvs*x1c*x2 # -6*vvs*x1c*x2s-8./9.*vvs*x1c*x2c+8./9.* # vvs*x1q*x2s-8./3.*vvs*x2*z-7./3.*vvs*x2*zs- # 2./3.*vvs*x2*zc-vvs*x2-1./3.*vvs*x2s*z-11./9. # *vvs*x2s*zs-2./3.*vvs*x2s*zc+2./9.*vvs* # x2s+16./3.*vvs*x2c*z+2*vvs*x2c*zs+34./9.* # vvs*x2c+4./9.*vvs*x2q*z+4./9.*vvs*x2q-5* # vvq*x1*x2*z-vvq*x1*x2*zs-14./3.*vvq*x1*x2-2* # vvq*x1*x2s*z) dw1= dw1+yc*(-10./9.*vvq*x1*x2s+22. # /9.*vvq*x1*x2c+2*vvq*x1s*x2*z+10./3.*vvq*x1s* # x2+2./9.*vvq*x1s*x2s-8./9.*vvq*x1s*x2c-2./3. # *vvq*x1c*x2+3*vvq*x2*z+4./3.*vvq*x2*zs+2* # vvq*x2+8./3.*vvq*x2s*z+5./3.*vvq*x2s*zs+8./ # 9.*vvq*x2s-4./3.*vvq*x2c*z-14./9.*vvq*x2c) dw2= y*rwms*rwmpgs*(8./3.+22./9.*vvs*x1*x2 # +2./3.*vvs*x1*z+2./9.*vvs*x1-8./9.*vvs*x1s*x2 # -4./9.*vvs*x1s-4./3.*vvs*x2*z-14./9.*vvs*x2+2. # /9.*vvs-28./3.*x1*x2*z-56./3.*x1*x2-4./3.*x1*x2s- # 12*x1*z-2*x1*zs-14*x1+14*x1s*x2+8*x1s*z+68./3. # *x1s-8./3.*x1c*x2-14*x1c+8./3.*x1q+8*x2*z+2 # *x2*zs+22./3.*x2+4./3.*x2s*z+4./3.*x2s+4*z+4./ # 3.*zs) dw2= dw2+ys*rwms*(25*vv*x1*x2*z+4*vv* # x1*x2*zs+28*vv*x1*x2+80./3.*vv*x1*x2s*z+136./3.*vv* # x1*x2s+8./3.*vv*x1*x2c-19*vv*x1*z-11*vv*x1*zs-2 # *vv*x1*zc-10*vv*x1-22*vv*x1s*x2*z-154./3.*vv*x1s* # x2-32*vv*x1s*x2s+13*vv*x1s*z+4*vv*x1s*zs+11 # *vv*x1s+32*vv*x1c*x2+16./3.*vv*x1c*x2s-2*vv* # x1c*z-4*vv*x1c-16./3.*vv*x1q*x2-3*vv*x2*z+7./3. # *vv*x2*zs+2*vv*x2*zc-10./3.*vv*x2-24*vv*x2s*z-8 # *vv*x2s*zs-56./3.*vv*x2s-8./3.*vv*x2c*z-8./3.* # vv*x2c+8*vv*z+7*vv*zs+2*vv*zc+3*vv+4./3.* # vvc*x1*x2*z+8./9.*vvc*x1*x2-44./9.*vvc*x1*x2s+5 # *vvc*x1*z+vvc*x1*zs+14./3.*vvc*x1+2./9.*vvc* # x1s*x2+16./9.*vvc*x1s*x2s-2*vvc*x1s*z-10./3. # *vvc*x1s+2./3.*vvc*x1c-8./3.*vvc*x2*z-5./3.* # vvc*x2*zs-10./9.*vvc*x2+8./3.*vvc*x2s*z+28./9. # *vvc*x2s) dw2= dw2+ys*rwms*(-3*vvc*z-4./3.* # vvc*zs-2*vvc) dw2= dw2+yc*(19*vvs*x1*x2*z+11*vvs* # x1*x2*zs+2*vvs*x1*x2*zc+10*vvs*x1*x2-13*vvs* # x1*x2s*z-2*vvs*x1*x2s*zs-14*vvs*x1*x2s-52./ # 3.*vvs*x1*x2c*z-80./3.*vvs*x1*x2c-4./3.*vvs*x1* # x2q-13*vvs*x1s*x2*z-4*vvs*x1s*x2*zs-11* # vvs*x1s*x2+14*vvs*x1s*x2s*z+86./3.*vvs*x1s* # x2s+18*vvs*x1s*x2c+2*vvs*x1c*x2*z+4*vvs* # x1c*x2-18*vvs*x1c*x2s-8./3.*vvs*x1c*x2c+8. # /3.*vvs*x1q*x2s-8*vvs*x2*z-7*vvs*x2*zs-2* # vvs*x2*zc-3*vvs*x2-vvs*x2s*z-11./3.*vvs* # x2s*zs-2*vvs*x2s*zc+2./3.*vvs*x2s+16* # vvs*x2c*z+6*vvs*x2c*zs+34./3.*vvs*x2c+4./ # 3.*vvs*x2q*z+4./3.*vvs*x2q-5*vvq*x1*x2*z- # vvq*x1*x2*zs-14./3.*vvq*x1*x2-2*vvq*x1*x2s*z- # 10./9.*vvq*x1*x2s+22./9.*vvq*x1*x2c+2*vvq*x1s # *x2*z) dw2= dw2+yc*(10./3.*vvq*x1s*x2+2./9. # *vvq*x1s*x2s-8./9.*vvq*x1s*x2c-2./3.*vvq* # x1c*x2+3*vvq*x2*z+4./3.*vvq*x2*zs+2*vvq*x2+ # 8./3.*vvq*x2s*z+5./3.*vvq*x2s*zs+8./9.*vvq* # x2s-4./3.*vvq*x2c*z-14./9.*vvq*x2c) * else if(jj.eq.5) then * wcom= y*rdms*rwms*(-8./9.*vvs*x1*x2*z-8./9.* # vvs*x1*x2+4./9.*vvs*x1s*x2+8./9.*vvs*x2*z+4./9. # *vvs*x2*zs+4./9.*vvs*x2) wcom= wcom+ys*rdms*(4./9.*vvc*x1*x2s*z # +4./9.*vvc*x1*x2s-2./9.*vvc*x1s*x2s-4./9.* # vvc*x2s*z-2./9.*vvc*x2s*zs-2./9.*vvc*x2s) wcom= wcom+rdms*rwms*rwmpgs*(4./9.*vv*x1*z # +4./9.*vv*x1-2./9.*vv*x1s-4./9.*vv*z-2./9.*vv*zs # -2./9.*vv) * dw1= y*rdms*rwms*(-2./3.*vvs*x1*x2+2./3.* # vvs*x2*z+2./3.*vvs*x2+28./9.*x1*x2*z+8./9.*x1*x2* # zs+20./9.*x1*x2+8./3.*x1*x2s*z+28./9.*x1*x2s+4./ # 9.*x1*x2c-20./9.*x1s*x2*z-28./9.*x1s*x2-16./9.* # x1s*x2s+4./3.*x1c*x2-8./9.*x2*z-4./9.*x2*zs-4. # /9.*x2-20./9.*x2s*z-8./9.*x2s*zs-4./3.*x2s-4./ # 9.*x2c*z-4./9.*x2c) dw1= dw1+y*rwms*rwmpgs*(-2./9.*vvs* # x1*z-2./9.*vvs*x1+1./9.*vvs*x1s+2./9.*vvs*z+1. # /9.*vvs*zs+1./9.*vvs-4./3.*x1*x2*z-4./9.*x1*x2* # zs-8./9.*x1*x2-4./9.*x1*x2s*z-4./9.*x1*x2s-4./9. # *x1*z-2./9.*x1*zs-2./9.*x1+8./9.*x1s*x2*z+10./9.* # x1s*x2+2./9.*x1s*x2s+8./9.*x1s*z+2./9.*x1s* # zs+2./3.*x1s-4./9.*x1c*x2-4./9.*x1c*z-2./3.* # x1c+2./9.*x1q+4./9.*x2*z+2./9.*x2*zs+2./9.*x2 # +4./9.*x2s*z+2./9.*x2s*zs+2./9.*x2s) dw1= dw1+ys*rdms*(-14./9.*vv*x1*x2s* # z-4./9.*vv*x1*x2s*zs-10./9.*vv*x1*x2s-4./3.*vv*x1 # *x2c*z-14./9.*vv*x1*x2c-2./9.*vv*x1*x2q+10./9.*vv # *x1s*x2s*z+14./9.*vv*x1s*x2s+8./9.*vv*x1s*x2c # -2./3.*vv*x1c*x2s+4./9.*vv*x2s*z+2./9.*vv*x2s* # zs+2./9.*vv*x2s+10./9.*vv*x2c*z+4./9.*vv*x2c* # zs+2./3.*vv*x2c+2./9.*vv*x2q*z+2./9.*vv*x2q+1. # /3.*vvc*x1*x2s-1./3.*vvc*x2s*z-1./3.*vvc*x2s # ) dw1= dw1+ys*rwms*(8./9.*vv*x1*x2*z+4. # /9.*vv*x1*x2*zs+4./9.*vv*x1*x2+8./3.*vv*x1*x2s*z+8./ # 9.*vv*x1*x2s*zs+16./9.*vv*x1*x2s+8./9.*vv*x1*x2c* # z+8./9.*vv*x1*x2c-16./9.*vv*x1s*x2*z-4./9.*vv*x1s # *x2*zs-4./3.*vv*x1s*x2-16./9.*vv*x1s*x2s*z-20./ # 9.*vv*x1s*x2s-4./9.*vv*x1s*x2c+8./9.*vv*x1c*x2* # z+4./3.*vv*x1c*x2+8./9.*vv*x1c*x2s-4./9.*vv*x1q # *x2-8./9.*vv*x2s*z-4./9.*vv*x2s*zs-4./9.*vv*x2s # -8./9.*vv*x2c*z-4./9.*vv*x2c*zs-4./9.*vv*x2c+ # 4./9.*vvc*x1*x2*z+4./9.*vvc*x1*x2-2./9.*vvc*x1s* # x2-4./9.*vvc*x2*z-2./9.*vvc*x2*zs-2./9.*vvc*x2 # ) dw1= dw1+yc*(-4./9.*vvs*x1*x2s*z # -2./9.*vvs*x1*x2s*zs-2./9.*vvs*x1*x2s-4./3.* # vvs*x1*x2c*z-4./9.*vvs*x1*x2c*zs-8./9.*vvs*x1 # *x2c-4./9.*vvs*x1*x2q*z-4./9.*vvs*x1*x2q+8./9. # *vvs*x1s*x2s*z+2./9.*vvs*x1s*x2s*zs+2./3.* # vvs*x1s*x2s+8./9.*vvs*x1s*x2c*z+10./9.*vvs* # x1s*x2c+2./9.*vvs*x1s*x2q-4./9.*vvs*x1c* # x2s*z-2./3.*vvs*x1c*x2s-4./9.*vvs*x1c*x2c # +2./9.*vvs*x1q*x2s+4./9.*vvs*x2c*z+2./9.* # vvs*x2c*zs+2./9.*vvs*x2c+4./9.*vvs*x2q*z+ # 2./9.*vvs*x2q*zs+2./9.*vvs*x2q-2./9.*vvq*x1* # x2s*z-2./9.*vvq*x1*x2s+1./9.*vvq*x1s*x2s+2./ # 9.*vvq*x2s*z+1./9.*vvq*x2s*zs+1./9.*vvq*x2s # ) dw1= dw1+rdms*rwms*rwmpgs*(-4./3.* # vvi*x1*x2*z-14./9.*vvi*x1*x2-2./9.*vvi*x1* # x2s-14./9.*vvi*x1*z-4./9.*vvi*x1*zs-10./9. # *vvi*x1+8./9.*vvi*x1s*x2+10./9.*vvi* # x1s*z+14./9.*vvi*x1s-2./3.*vvi*x1c+10./ # 9.*vvi*x2*z+4./9.*vvi*x2*zs+2./3.*vvi*x2 # +2./9.*vvi*x2s*z+2./9.*vvi*x2s+4./9.* # vvi*z+2./9.*vvi*zs+2./9.*vvi+1./3.*vv* # x1-1./3.*vv*z-1./3.*vv) dw2= y*rdms*rwms*(-2./3.*vvs*x1*x2+2./3.* # vvs*x2*z+2./3.*vvs*x2+28./3.*x1*x2*z+8./3.*x1*x2* # zs+20./3.*x1*x2+8*x1*x2s*z+28./3.*x1*x2s+4./3.* # x1*x2c-20./3.*x1s*x2*z-28./3.*x1s*x2-16./3.*x1s # *x2s+4*x1c*x2-8./3.*x2*z-4./3.*x2*zs-4./3.*x2 # -20./3.*x2s*z-8./3.*x2s*zs-4*x2s-4./3.*x2c* # z-4./3.*x2c) dw2= dw2+y*rwms*rwmpgs*(-2./9.*vvs* # x1*z-2./9.*vvs*x1+1./9.*vvs*x1s+2./9.*vvs*z+1. # /9.*vvs*zs+1./9.*vvs-4*x1*x2*z-4./3.*x1*x2*zs # -8./3.*x1*x2-4./3.*x1*x2s*z-4./3.*x1*x2s-4./3.*x1 # *z-2./3.*x1*zs-2./3.*x1+8./3.*x1s*x2*z+10./3.* # x1s*x2+2./3.*x1s*x2s+8./3.*x1s*z+2./3.*x1s* # zs+2*x1s-4./3.*x1c*x2-4./3.*x1c*z-2*x1c+2. # /3.*x1q+4./3.*x2*z+2./3.*x2*zs+2./3.*x2+4./3.* # x2s*z+2./3.*x2s*zs+2./3.*x2s) dw2= dw2+ys*rdms*(-14./3.*vv*x1*x2s* # z-4./3.*vv*x1*x2s*zs-10./3.*vv*x1*x2s-4*vv*x1* # x2c*z-14./3.*vv*x1*x2c-2./3.*vv*x1*x2q+10./3.*vv* # x1s*x2s*z+14./3.*vv*x1s*x2s+8./3.*vv*x1s*x2c # -2*vv*x1c*x2s+4./3.*vv*x2s*z+2./3.*vv*x2s*zs # +2./3.*vv*x2s+10./3.*vv*x2c*z+4./3.*vv*x2c*zs # +2*vv*x2c+2./3.*vv*x2q*z+2./3.*vv*x2q+1./3.* # vvc*x1*x2s-1./3.*vvc*x2s*z-1./3.*vvc*x2s) dw2= dw2+ys*rwms*(8./3.*vv*x1*x2*z+4. # /3.*vv*x1*x2*zs+4./3.*vv*x1*x2+8*vv*x1*x2s*z+8./3.* # vv*x1*x2s*zs+16./3.*vv*x1*x2s+8./3.*vv*x1*x2c*z # +8./3.*vv*x1*x2c-16./3.*vv*x1s*x2*z-4./3.*vv*x1s* # x2*zs-4*vv*x1s*x2-16./3.*vv*x1s*x2s*z-20./3.*vv # *x1s*x2s-4./3.*vv*x1s*x2c+8./3.*vv*x1c*x2*z+4 # *vv*x1c*x2+8./3.*vv*x1c*x2s-4./3.*vv*x1q*x2-8./ # 3.*vv*x2s*z-4./3.*vv*x2s*zs-4./3.*vv*x2s-8./3.* # vv*x2c*z-4./3.*vv*x2c*zs-4./3.*vv*x2c+4./9.* # vvc*x1*x2*z+4./9.*vvc*x1*x2-2./9.*vvc*x1s*x2-4./ # 9.*vvc*x2*z-2./9.*vvc*x2*zs-2./9.*vvc*x2) dw2= dw2+yc*(-4./3.*vvs*x1*x2s*z # -2./3.*vvs*x1*x2s*zs-2./3.*vvs*x1*x2s-4*vvs # *x1*x2c*z-4./3.*vvs*x1*x2c*zs-8./3.*vvs*x1* # x2c-4./3.*vvs*x1*x2q*z-4./3.*vvs*x1*x2q+8./3. # *vvs*x1s*x2s*z+2./3.*vvs*x1s*x2s*zs+2*vvs # *x1s*x2s+8./3.*vvs*x1s*x2c*z+10./3.*vvs*x1s # *x2c+2./3.*vvs*x1s*x2q-4./3.*vvs*x1c*x2s*z # -2*vvs*x1c*x2s-4./3.*vvs*x1c*x2c+2./3.* # vvs*x1q*x2s+4./3.*vvs*x2c*z+2./3.*vvs*x2c* # zs+2./3.*vvs*x2c+4./3.*vvs*x2q*z+2./3.*vvs* # x2q*zs+2./3.*vvs*x2q-2./9.*vvq*x1*x2s*z-2./ # 9.*vvq*x1*x2s+1./9.*vvq*x1s*x2s+2./9.*vvq* # x2s*z+1./9.*vvq*x2s*zs+1./9.*vvq*x2s) dw2= dw2+rdms*rwms*rwmpgs*(-4*vvi # *x1*x2*z-14./3.*vvi*x1*x2-2./3.*vvi*x1*x2s- # 14./3.*vvi*x1*z-4./3.*vvi*x1*zs-10./3.* # vvi*x1+8./3.*vvi*x1s*x2+10./3.*vvi*x1s # *z+14./3.*vvi*x1s-2*vvi*x1c+10./3.* # vvi*x2*z+4./3.*vvi*x2*zs+2*vvi*x2+2./3. # *vvi*x2s*z+2./3.*vvi*x2s+4./3.*vvi*z # +2./3.*vvi*zs+2./3.*vvi+1./3.*vv*x1-1./3. # *vv*z-1./3.*vv) * endif * wcom= wcom*(2.d0-y)**2/ys-dw1+2.d0/y*(1.d0-1.d0/y)*dw2 * else if(ipro.eq.1.and.ii.le.2) then if(jj.eq.1) then * wcom= y*rwms*rwmpgs*(-vvs*tau*x1+vvs*tau*x2 # +2*vvs*x1*x2+vvs*x1*z+vvs*x1-2*vvs*x1s- # vvs*x2*z-vvs*x2) wcom= wcom+ys*rwms*(4*vvc*tau*x1*x2- # vvc*tau*x1*z-3*vvc*tau*x1+vvc*tau*x2*z-2*vvc* # tau*x2s+vvc*tau*z+vvc*tau+vvc*taus*x2-vvc # *taus-4*vvc*x1*x2-4*vvc*x1*x2s+4*vvc*x1*z+ # vvc*x1*zs+5*vvc*x1+6*vvc*x1s*x2-2*vvc*x1s # *z-4*vvc*x1s-vvc*x2*z-vvc*x2*zs+vvc*x2+ # 2*vvc*x2s*z+2*vvc*x2s-2*vvc*z-vvc*zs-2* # vvc) wcom= wcom+yc*(3*vvq*tau*x1*x2*z+3*vvq # *tau*x1*x2-3*vvq*tau*x1*x2s-vvq*tau*x2*z-vvq* # tau*x2-vvq*tau*x2s*z+vvq*tau*x2c+vvq*taus* # x2*z+vvq*taus*x2-vvq*taus*x2s-6*vvq*x1*x2*z # -vvq*x1*x2*zs-5*vvq*x1*x2-vvq*x1*x2s*z+3* # vvq*x1*x2s+2*vvq*x1*x2c+4*vvq*x1s*x2*z+4* # vvq*x1s*x2-4*vvq*x1s*x2s+3*vvq*x2*z+vvq* # x2*zs+2*vvq*x2+vvq*x2s*z+vvq*x2s*zs- # vvq*x2s-vvq*x2c*z-vvq*x2c) * else if(jj.eq.2) then wcom= 0.d0 * else if(jj.eq.3) then wcom= 0.d0 * else if(jj.eq.4) then * wcom= y*rwms*rwmpgs*(2*vvs*x1*x2*z+4*vvs*x1* # x2+2*vvs*x1*z+vvs*x1*zs+2*vvs*x1-2*vvs* # x1s*x2-2*vvs*x1s*z-4*vvs*x1s+2*vvs*x1c- # 2*vvs*x2*z-vvs*x2*zs-2*vvs*x2) wcom= wcom+ys*rwms*(-4*vvc*x1*x2-4* # vvc*x1*x2s*z-8*vvc*x1*x2s+10*vvc*x1*z+5*vvc # *x1*zs+vvc*x1*zc+6*vvc*x1+2*vvc*x1s*x2*z+ # 8*vvc*x1s*x2+4*vvc*x1s*x2s-8*vvc*x1s*z-2* # vvc*x1s*zs-6*vvc*x1s-4*vvc*x1c*x2+2*vvc # *x1c*z+2*vvc*x1c-2*vvc*x2*z-2*vvc*x2*zs- # vvc*x2*zc+4*vvc*x2s*z+2*vvc*x2s*zs+4* # vvc*x2s-4*vvc*z-3*vvc*zs-vvc*zc-2*vvc # ) wcom= wcom+yc*(-10*vvq*x1*x2*z-5*vvq # *x1*x2*zs-vvq*x1*x2*zc-6*vvq*x1*x2-2*vvq*x1* # x2s*z-vvq*x1*x2s*zs+2*vvq*x1*x2s+2*vvq*x1 # *x2c*z+4*vvq*x1*x2c+8*vvq*x1s*x2*z+2*vvq* # x1s*x2*zs+6*vvq*x1s*x2-4*vvq*x1s*x2s-2* # vvq*x1s*x2c-2*vvq*x1c*x2*z-2*vvq*x1c*x2+2 # *vvq*x1c*x2s+4*vvq*x2*z+3*vvq*x2*zs+vvq* # x2*zc+2*vvq*x2+2*vvq*x2s*z+2*vvq*x2s*zs # +vvq*x2s*zc-2*vvq*x2c*z-vvq*x2c*zs-2* # vvq*x2c) * else if(jj.eq.5) then * wcom= y*rdms*rwms*(-8*vvs*x1*x2*z-8*vvs*x1* # x2+4*vvs*x1s*x2+8*vvs*x2*z+4*vvs*x2*zs+4* # vvs*x2) wcom= wcom+ys*rdms*(4*vvc*x1*x2s*z+4 # *vvc*x1*x2s-2*vvc*x1s*x2s-4*vvc*x2s*z-2* # vvc*x2s*zs-2*vvc*x2s) wcom= wcom+rdms*rwms*rwmpgs*(4*vv*x1*z+4 # *vv*x1-2*vv*x1s-4*vv*z-2*vv*zs-2*vv) * endif else if(ipro.eq.1.and.ii.eq.3) then if(jj.eq.1) then * wcom= y*rwms*rwmpgs*(-vvs*tau*x1+vvs*tau*x2 # +2*vvs*x1*x2+vvs*x1*z+vvs*x1-2*vvs*x1s- # vvs*x2*z-vvs*x2) wcom= wcom+ys*rwms*(4*vvc*tau*x1*x2- # vvc*tau*x1*z-3*vvc*tau*x1+vvc*tau*x2*z-2*vvc* # tau*x2s+vvc*tau*z+vvc*tau+vvc*taus*x2-vvc # *taus-4*vvc*x1*x2-4*vvc*x1*x2s+4*vvc*x1*z+ # vvc*x1*zs+5*vvc*x1+6*vvc*x1s*x2-2*vvc*x1s # *z-4*vvc*x1s-vvc*x2*z-vvc*x2*zs+vvc*x2+ # 2*vvc*x2s*z+2*vvc*x2s-2*vvc*z-vvc*zs-2* # vvc) wcom= wcom+yc*(3*vvq*tau*x1*x2*z+3*vvq # *tau*x1*x2-3*vvq*tau*x1*x2s-vvq*tau*x2*z-vvq* # tau*x2-vvq*tau*x2s*z+vvq*tau*x2c+vvq*taus* # x2*z+vvq*taus*x2-vvq*taus*x2s-6*vvq*x1*x2*z # -vvq*x1*x2*zs-5*vvq*x1*x2-vvq*x1*x2s*z+3* # vvq*x1*x2s+2*vvq*x1*x2c+4*vvq*x1s*x2*z+4* # vvq*x1s*x2-4*vvq*x1s*x2s+3*vvq*x2*z+vvq* # x2*zs+2*vvq*x2+vvq*x2s*z+vvq*x2s*zs- # vvq*x2s-vvq*x2c*z-vvq*x2c) * dw1= y*rwms*rwmpgs*(vvs*x1s-2*vvs*x2+ # vvs+2*tau*x1-2*tau*x2-8*x1*x2-2*x1*z-4*x1+6* # x1s+2*x2*z+4*x2+2*x2s) dw1= dw1+ys*rwms*(-12*vv*tau*x1*x2+2 # *vv*tau*x1*z+8*vv*tau*x1-2*vv*tau*x2*z-vv*tau*x2+8*vv # *tau*x2s-2*vv*tau*z-3*vv*tau-2*vv*taus*x2+2*vv* # taus+4*vv*x1*x2*z+9*vv*x1*x2+24*vv*x1*x2s-9*vv*x1 # *z-2*vv*x1*zs-10*vv*x1-22*vv*x1s*x2+4*vv*x1s*z # +10*vv*x1s+3*vv*x2*z+2*vv*x2*zs+2*vv*x2-8*vv* # x2s*z-12*vv*x2s-4*vv*x2c+5*vv*z+2*vv*zs+3* # vv+vvc*tau*x2+vvc*x1*z+3*vvc*x1-2*vvc*x1s* # x2+2*vvc*x1s*z-3*vvc*x2*z-vvc*x2+4*vvc* # x2s-2*vvc*z-3*vvc) dw1= dw1+yc*(-6*vvs*tau*x1*x2*z-15* # vvs*tau*x1*x2+14*vvs*tau*x1*x2s+vvs*tau*x1*z+5./ # 2.*vvs*tau*x1+2*vvs*tau*x2*z+4*vvs*tau*x2+2*vvs # *tau*x2s*z+vvs*tau*x2s-6*vvs*tau*x2c-1./2.* # vvs*tau-2*vvs*taus*x2*z-5*vvs*taus*x2+4*vvs # *taus*x2s+vvs*taus+9*vvs*x1*x2*z+2*vvs*x1* # x2*zs+11*vvs*x1*x2-2*vvs*x1*x2s*z-5*vvs*x1* # x2s-16*vvs*x1*x2c-1./2.*vvs*x1-6*vvs*x1s*x2 # *z-14*vvs*x1s*x2+18*vvs*x1s*x2s+vvs*x1s*z # +3./2.*vvs*x1s-5*vvs*x2*z-2*vvs*x2*zs-3* # vvs*x2-3*vvs*x2s*z-2*vvs*x2s*zs-2*vvs* # x2s+6*vvs*x2c*z+8*vvs*x2c+2*vvs*x2q+2* # vvq*tau*x1*x2+vvq*tau*x1*z+1./2.*vvq*tau*x1-vvq # *tau*x2s-1./2.*vvq*tau+vvq*taus*x2-3*vvq*x1* # x2*z-5*vvq*x1*x2-vvq*x1*z-1./2.*vvq*x1-2*vvq* # x1s*x2*z) dw1= dw1+yc*(2*vvq*x1s*x2+vvq* # x1s*x2s+2*vvq*x1s*z+vvq*x1s*zs+1./2.* # vvq*x1s+4*vvq*x2*z+4*vvq*x2+3*vvq*x2s*z-2 # *vvq*x2c) dw2= y*rwms*rwmpgs*(vvs*x1s-2*vvs*x2+ # vvs+6*tau*x1-6*tau*x2-24*x1*x2-6*x1*z-12*x1+18* # x1s+6*x2*z+12*x2+6*x2s) dw2= dw2+ys*rwms*(-36*vv*tau*x1*x2+6 # *vv*tau*x1*z+24*vv*tau*x1-6*vv*tau*x2*z-3*vv*tau*x2+ # 24*vv*tau*x2s-6*vv*tau*z-9*vv*tau-6*vv*taus*x2+6* # vv*taus+12*vv*x1*x2*z+27*vv*x1*x2+72*vv*x1*x2s-27 # *vv*x1*z-6*vv*x1*zs-30*vv*x1-66*vv*x1s*x2+12*vv* # x1s*z+30*vv*x1s+9*vv*x2*z+6*vv*x2*zs+6*vv*x2- # 24*vv*x2s*z-36*vv*x2s-12*vv*x2c+15*vv*z+6*vv* # zs+9*vv+vvc*tau*x2+vvc*x1*z+3*vvc*x1-2* # vvc*x1s*x2+2*vvc*x1s*z-3*vvc*x2*z-vvc*x2+ # 4*vvc*x2s-2*vvc*z-3*vvc) dw2= dw2+yc*(-18*vvs*tau*x1*x2*z-45* # vvs*tau*x1*x2+42*vvs*tau*x1*x2s+3*vvs*tau*x1*z+ # 15./2.*vvs*tau*x1+6*vvs*tau*x2*z+12*vvs*tau*x2+6* # vvs*tau*x2s*z+3*vvs*tau*x2s-18*vvs*tau*x2c- # 3./2.*vvs*tau-6*vvs*taus*x2*z-15*vvs*taus*x2+ # 12*vvs*taus*x2s+3*vvs*taus+27*vvs*x1*x2*z+6 # *vvs*x1*x2*zs+33*vvs*x1*x2-6*vvs*x1*x2s*z-15* # vvs*x1*x2s-48*vvs*x1*x2c-3./2.*vvs*x1-18* # vvs*x1s*x2*z-42*vvs*x1s*x2+54*vvs*x1s*x2s # +3*vvs*x1s*z+9./2.*vvs*x1s-15*vvs*x2*z-6* # vvs*x2*zs-9*vvs*x2-9*vvs*x2s*z-6*vvs*x2s* # zs-6*vvs*x2s+18*vvs*x2c*z+24*vvs*x2c+6* # vvs*x2q+2*vvq*tau*x1*x2+vvq*tau*x1*z+1./2.* # vvq*tau*x1-vvq*tau*x2s-1./2.*vvq*tau+vvq* # taus*x2-3*vvq*x1*x2*z-5*vvq*x1*x2-vvq*x1*z-1./ # 2.*vvq*x1) dw2= dw2+yc*(-2*vvq*x1s*x2*z+2* # vvq*x1s*x2+vvq*x1s*x2s+2*vvq*x1s*z+vvq* # x1s*zs+1./2.*vvq*x1s+4*vvq*x2*z+4*vvq*x2+ # 3*vvq*x2s*z-2*vvq*x2c) * else if(jj.eq.2) then wcom= 0.d0 dw1= 0.d0 dw2= 0.d0 * else if(jj.eq.3) then wcom= 0.d0 dw1= 0.d0 dw2= 0.d0 * else if(jj.eq.4) then * wcom= y*rwms*rwmpgs*(2*vvs*x1*x2*z+4*vvs*x1* # x2+2*vvs*x1*z+vvs*x1*zs+2*vvs*x1-2*vvs* # x1s*x2-2*vvs*x1s*z-4*vvs*x1s+2*vvs*x1c- # 2*vvs*x2*z-vvs*x2*zs-2*vvs*x2) wcom= wcom+ys*rwms*(-4*vvc*x1*x2-4* # vvc*x1*x2s*z-8*vvc*x1*x2s+10*vvc*x1*z+5*vvc # *x1*zs+vvc*x1*zc+6*vvc*x1+2*vvc*x1s*x2*z+ # 8*vvc*x1s*x2+4*vvc*x1s*x2s-8*vvc*x1s*z-2* # vvc*x1s*zs-6*vvc*x1s-4*vvc*x1c*x2+2*vvc # *x1c*z+2*vvc*x1c-2*vvc*x2*z-2*vvc*x2*zs- # vvc*x2*zc+4*vvc*x2s*z+2*vvc*x2s*zs+4* # vvc*x2s-4*vvc*z-3*vvc*zs-vvc*zc-2*vvc # ) wcom= wcom+yc*(-10*vvq*x1*x2*z-5*vvq # *x1*x2*zs-vvq*x1*x2*zc-6*vvq*x1*x2-2*vvq*x1* # x2s*z-vvq*x1*x2s*zs+2*vvq*x1*x2s+2*vvq*x1 # *x2c*z+4*vvq*x1*x2c+8*vvq*x1s*x2*z+2*vvq* # x1s*x2*zs+6*vvq*x1s*x2-4*vvq*x1s*x2s-2* # vvq*x1s*x2c-2*vvq*x1c*x2*z-2*vvq*x1c*x2+2 # *vvq*x1c*x2s+4*vvq*x2*z+3*vvq*x2*zs+vvq* # x2*zc+2*vvq*x2+2*vvq*x2s*z+2*vvq*x2s*zs # +vvq*x2s*zc-2*vvq*x2c*z-vvq*x2c*zs-2* # vvq*x2c) * dw1= y*rwms*rwmpgs*(6*vvs*x1*x2+2*vvs*x1*z # +2*vvs*x1-4*vvs*x1s-4*vvs*x2*z-6*vvs*x2+2 # *vvs*z+2*vvs-12*x1*x2*z-16*x1*x2-4*x1*x2s-8* # x1*z-2*x1*zs-6*x1+10*x1s*x2+8*x1s*z+12*x1s # -6*x1c+8*x2*z+2*x2*zs+6*x2+4*x2s*z+4*x2s # ) dw1= dw1+ys*rwms*(17*vv*x1*x2*z+4*vv* # x1*x2*zs+12*vv*x1*x2+32*vv*x1*x2s*z+40*vv*x1*x2s # +8*vv*x1*x2c-19*vv*x1*z-11*vv*x1*zs-2*vv*x1*zc # -10*vv*x1-22*vv*x1s*x2*z-30*vv*x1s*x2-24*vv*x1s # *x2s+13*vv*x1s*z+4*vv*x1s*zs+11*vv*x1s+16* # vv*x1c*x2-2*vv*x1c*z-4*vv*x1c+5*vv*x2*z+5*vv*x2 # *zs+2*vv*x2*zc+2*vv*x2-24*vv*x2s*z-8*vv*x2s* # zs-16*vv*x2s-8*vv*x2c*z-8*vv*x2c+8*vv*z+7* # vv*zs+2*vv*zc+3*vv+4*vvc*x1*x2*z-12*vvc*x1* # x2s+15*vvc*x1*z+3*vvc*x1*zs+14*vvc*x1+6* # vvc*x1s*x2-6*vvc*x1s*z-10*vvc*x1s+2*vvc* # x1c-12*vvc*x2*z-5*vvc*x2*zs-6*vvc*x2+8* # vvc*x2s*z+12*vvc*x2s-9*vvc*z-4*vvc*zs-6 # *vvc) dw1= dw1+yc*(19*vvs*x1*x2*z+11*vvs* # x1*x2*zs+2*vvs*x1*x2*zc+10*vvs*x1*x2-9*vvs*x1 # *x2s*z-2*vvs*x1*x2s*zs-6*vvs*x1*x2s-20* # vvs*x1*x2c*z-24*vvs*x1*x2c-4*vvs*x1*x2q-13* # vvs*x1s*x2*z-4*vvs*x1s*x2*zs-11*vvs*x1s*x2 # +14*vvs*x1s*x2s*z+18*vvs*x1s*x2s+14*vvs* # x1s*x2c+2*vvs*x1c*x2*z+4*vvs*x1c*x2-10* # vvs*x1c*x2s-8*vvs*x2*z-7*vvs*x2*zs-2*vvs* # x2*zc-3*vvs*x2-5*vvs*x2s*z-5*vvs*x2s*zs # -2*vvs*x2s*zc-2*vvs*x2s+16*vvs*x2c*z+6* # vvs*x2c*zs+10*vvs*x2c+4*vvs*x2q*z+4*vvs # *x2q-15*vvq*x1*x2*z-3*vvq*x1*x2*zs-14*vvq*x1* # x2-6*vvq*x1*x2s*z-2*vvq*x1*x2s+6*vvq*x1*x2c # +6*vvq*x1s*x2*z+10*vvq*x1s*x2-2*vvq*x1s* # x2s-2*vvq*x1c*x2+9*vvq*x2*z+4*vvq*x2*zs+6 # *vvq*x2) dw1= dw1+yc*(10*vvq*x2s*z+5*vvq* # x2s*zs+4*vvq*x2s-4*vvq*x2c*z-6*vvq*x2c) dw2= y*rwms*rwmpgs*(6*vvs*x1*x2+2*vvs*x1*z # +2*vvs*x1-4*vvs*x1s-4*vvs*x2*z-6*vvs*x2+2 # *vvs*z+2*vvs-36*x1*x2*z-48*x1*x2-12*x1*x2s-24 # *x1*z-6*x1*zs-18*x1+30*x1s*x2+24*x1s*z+36* # x1s-18*x1c+24*x2*z+6*x2*zs+18*x2+12*x2s*z # +12*x2s) dw2= dw2+ys*rwms*(51*vv*x1*x2*z+12*vv* # x1*x2*zs+36*vv*x1*x2+96*vv*x1*x2s*z+120*vv*x1*x2s # +24*vv*x1*x2c-57*vv*x1*z-33*vv*x1*zs-6*vv*x1*zc # -30*vv*x1-66*vv*x1s*x2*z-90*vv*x1s*x2-72*vv*x1s # *x2s+39*vv*x1s*z+12*vv*x1s*zs+33*vv*x1s+48* # vv*x1c*x2-6*vv*x1c*z-12*vv*x1c+15*vv*x2*z+15*vv # *x2*zs+6*vv*x2*zc+6*vv*x2-72*vv*x2s*z-24*vv* # x2s*zs-48*vv*x2s-24*vv*x2c*z-24*vv*x2c+24* # vv*z+21*vv*zs+6*vv*zc+9*vv+4*vvc*x1*x2*z-12* # vvc*x1*x2s+15*vvc*x1*z+3*vvc*x1*zs+14*vvc* # x1+6*vvc*x1s*x2-6*vvc*x1s*z-10*vvc*x1s+2* # vvc*x1c-12*vvc*x2*z-5*vvc*x2*zs-6*vvc*x2+ # 8*vvc*x2s*z+12*vvc*x2s-9*vvc*z-4*vvc*zs # -6*vvc) dw2= dw2+yc*(57*vvs*x1*x2*z+33*vvs* # x1*x2*zs+6*vvs*x1*x2*zc+30*vvs*x1*x2-27*vvs* # x1*x2s*z-6*vvs*x1*x2s*zs-18*vvs*x1*x2s-60* # vvs*x1*x2c*z-72*vvs*x1*x2c-12*vvs*x1*x2q-39 # *vvs*x1s*x2*z-12*vvs*x1s*x2*zs-33*vvs*x1s* # x2+42*vvs*x1s*x2s*z+54*vvs*x1s*x2s+42*vvs # *x1s*x2c+6*vvs*x1c*x2*z+12*vvs*x1c*x2-30* # vvs*x1c*x2s-24*vvs*x2*z-21*vvs*x2*zs-6* # vvs*x2*zc-9*vvs*x2-15*vvs*x2s*z-15*vvs* # x2s*zs-6*vvs*x2s*zc-6*vvs*x2s+48*vvs* # x2c*z+18*vvs*x2c*zs+30*vvs*x2c+12*vvs* # x2q*z+12*vvs*x2q-15*vvq*x1*x2*z-3*vvq*x1*x2* # zs-14*vvq*x1*x2-6*vvq*x1*x2s*z-2*vvq*x1*x2s # +6*vvq*x1*x2c+6*vvq*x1s*x2*z+10*vvq*x1s*x2 # -2*vvq*x1s*x2s-2*vvq*x1c*x2+9*vvq*x2*z+4* # vvq*x2*zs) dw2= dw2+yc*(6*vvq*x2+10*vvq*x2s*z # +5*vvq*x2s*zs+4*vvq*x2s-4*vvq*x2c*z-6* # vvq*x2c) * else if(jj.eq.5) then * wcom= y*rdms*rwms*(-8*vvs*x1*x2*z-8*vvs*x1* # x2+4*vvs*x1s*x2+8*vvs*x2*z+4*vvs*x2*zs+4* # vvs*x2) wcom= wcom+ys*rdms*(4*vvc*x1*x2s*z+4 # *vvc*x1*x2s-2*vvc*x1s*x2s-4*vvc*x2s*z-2* # vvc*x2s*zs-2*vvc*x2s) wcom= wcom+rdms*rwms*rwmpgs*(4*vv*x1*z+4 # *vv*x1-2*vv*x1s-4*vv*z-2*vv*zs-2*vv) * dw1= y*rdms*rwms*(-6*vvs*x1*x2+6*vvs*x2* # z+6*vvs*x2+28*x1*x2*z+8*x1*x2*zs+20*x1*x2+24*x1 # *x2s*z+28*x1*x2s+4*x1*x2c-20*x1s*x2*z-28* # x1s*x2-16*x1s*x2s+12*x1c*x2-8*x2*z-4*x2*zs # -4*x2-20*x2s*z-8*x2s*zs-12*x2s-4*x2c*z- # 4*x2c) dw1= dw1+y*rwms*rwmpgs*(-2*vvs*x1*z # -2*vvs*x1+vvs*x1s+2*vvs*z+vvs*zs+vvs # -12*x1*x2*z-4*x1*x2*zs-8*x1*x2-4*x1*x2s*z-4*x1* # x2s-4*x1*z-2*x1*zs-2*x1+8*x1s*x2*z+10*x1s* # x2+2*x1s*x2s+8*x1s*z+2*x1s*zs+6*x1s-4* # x1c*x2-4*x1c*z-6*x1c+2*x1q+4*x2*z+2*x2*zs # +2*x2+4*x2s*z+2*x2s*zs+2*x2s) dw1= dw1+ys*rdms*(-14*vv*x1*x2s*z- # 4*vv*x1*x2s*zs-10*vv*x1*x2s-12*vv*x1*x2c*z-14* # vv*x1*x2c-2*vv*x1*x2q+10*vv*x1s*x2s*z+14*vv* # x1s*x2s+8*vv*x1s*x2c-6*vv*x1c*x2s+4*vv* # x2s*z+2*vv*x2s*zs+2*vv*x2s+10*vv*x2c*z+4*vv # *x2c*zs+6*vv*x2c+2*vv*x2q*z+2*vv*x2q+3* # vvc*x1*x2s-3*vvc*x2s*z-3*vvc*x2s) dw1= dw1+ys*rwms*(8*vv*x1*x2*z+4*vv* # x1*x2*zs+4*vv*x1*x2+24*vv*x1*x2s*z+8*vv*x1*x2s* # zs+16*vv*x1*x2s+8*vv*x1*x2c*z+8*vv*x1*x2c-16* # vv*x1s*x2*z-4*vv*x1s*x2*zs-12*vv*x1s*x2-16*vv* # x1s*x2s*z-20*vv*x1s*x2s-4*vv*x1s*x2c+8*vv* # x1c*x2*z+12*vv*x1c*x2+8*vv*x1c*x2s-4*vv*x1q* # x2-8*vv*x2s*z-4*vv*x2s*zs-4*vv*x2s-8*vv*x2c # *z-4*vv*x2c*zs-4*vv*x2c+4*vvc*x1*x2*z+4*vvc # *x1*x2-2*vvc*x1s*x2-4*vvc*x2*z-2*vvc*x2*zs- # 2*vvc*x2) dw1= dw1+yc*(-4*vvs*x1*x2s*z-2* # vvs*x1*x2s*zs-2*vvs*x1*x2s-12*vvs*x1*x2c*z # -4*vvs*x1*x2c*zs-8*vvs*x1*x2c-4*vvs*x1* # x2q*z-4*vvs*x1*x2q+8*vvs*x1s*x2s*z+2*vvs* # x1s*x2s*zs+6*vvs*x1s*x2s+8*vvs*x1s*x2c* # z+10*vvs*x1s*x2c+2*vvs*x1s*x2q-4*vvs* # x1c*x2s*z-6*vvs*x1c*x2s-4*vvs*x1c*x2c+2 # *vvs*x1q*x2s+4*vvs*x2c*z+2*vvs*x2c*zs+2 # *vvs*x2c+4*vvs*x2q*z+2*vvs*x2q*zs+2*vvs # *x2q-2*vvq*x1*x2s*z-2*vvq*x1*x2s+vvq*x1s* # x2s+2*vvq*x2s*z+vvq*x2s*zs+vvq*x2s) dw1= dw1+rdms*rwms*rwmpgs*(-12* # vvi*x1*x2*z-14*vvi*x1*x2-2*vvi*x1*x2s- # 14*vvi*x1*z-4*vvi*x1*zs-10*vvi*x1+8* # vvi*x1s*x2+10*vvi*x1s*z+14*vvi*x1s # -6*vvi*x1c+10*vvi*x2*z+4*vvi*x2*zs # +6*vvi*x2+2*vvi*x2s*z+2*vvi*x2s+4* # vvi*z+2*vvi*zs+2*vvi+3*vv*x1-3*vv*z # -3*vv) dw2= y*rdms*rwms*(-6*vvs*x1*x2+6*vvs*x2* # z+6*vvs*x2+84*x1*x2*z+24*x1*x2*zs+60*x1*x2+72* # x1*x2s*z+84*x1*x2s+12*x1*x2c-60*x1s*x2*z-84* # x1s*x2-48*x1s*x2s+36*x1c*x2-24*x2*z-12*x2* # zs-12*x2-60*x2s*z-24*x2s*zs-36*x2s-12* # x2c*z-12*x2c) dw2= dw2+y*rwms*rwmpgs*(-2*vvs*x1*z # -2*vvs*x1+vvs*x1s+2*vvs*z+vvs*zs+vvs # -36*x1*x2*z-12*x1*x2*zs-24*x1*x2-12*x1*x2s*z-12 # *x1*x2s-12*x1*z-6*x1*zs-6*x1+24*x1s*x2*z+30* # x1s*x2+6*x1s*x2s+24*x1s*z+6*x1s*zs+18* # x1s-12*x1c*x2-12*x1c*z-18*x1c+6*x1q+12*x2 # *z+6*x2*zs+6*x2+12*x2s*z+6*x2s*zs+6*x2s) dw2= dw2+ys*rdms*(-42*vv*x1*x2s*z- # 12*vv*x1*x2s*zs-30*vv*x1*x2s-36*vv*x1*x2c*z-42* # vv*x1*x2c-6*vv*x1*x2q+30*vv*x1s*x2s*z+42*vv* # x1s*x2s+24*vv*x1s*x2c-18*vv*x1c*x2s+12*vv* # x2s*z+6*vv*x2s*zs+6*vv*x2s+30*vv*x2c*z+12* # vv*x2c*zs+18*vv*x2c+6*vv*x2q*z+6*vv*x2q+3* # vvc*x1*x2s-3*vvc*x2s*z-3*vvc*x2s) dw2= dw2+ys*rwms*(24*vv*x1*x2*z+12* # vv*x1*x2*zs+12*vv*x1*x2+72*vv*x1*x2s*z+24*vv*x1* # x2s*zs+48*vv*x1*x2s+24*vv*x1*x2c*z+24*vv*x1* # x2c-48*vv*x1s*x2*z-12*vv*x1s*x2*zs-36*vv*x1s* # x2-48*vv*x1s*x2s*z-60*vv*x1s*x2s-12*vv*x1s* # x2c+24*vv*x1c*x2*z+36*vv*x1c*x2+24*vv*x1c*x2s # -12*vv*x1q*x2-24*vv*x2s*z-12*vv*x2s*zs-12*vv* # x2s-24*vv*x2c*z-12*vv*x2c*zs-12*vv*x2c+4* # vvc*x1*x2*z+4*vvc*x1*x2-2*vvc*x1s*x2-4*vvc*x2 # *z-2*vvc*x2*zs-2*vvc*x2) dw2= dw2+yc*(-12*vvs*x1*x2s*z-6* # vvs*x1*x2s*zs-6*vvs*x1*x2s-36*vvs*x1*x2c*z # -12*vvs*x1*x2c*zs-24*vvs*x1*x2c-12*vvs*x1* # x2q*z-12*vvs*x1*x2q+24*vvs*x1s*x2s*z+6* # vvs*x1s*x2s*zs+18*vvs*x1s*x2s+24*vvs* # x1s*x2c*z+30*vvs*x1s*x2c+6*vvs*x1s*x2q- # 12*vvs*x1c*x2s*z-18*vvs*x1c*x2s-12*vvs* # x1c*x2c+6*vvs*x1q*x2s+12*vvs*x2c*z+6* # vvs*x2c*zs+6*vvs*x2c+12*vvs*x2q*z+6*vvs # *x2q*zs+6*vvs*x2q-2*vvq*x1*x2s*z-2*vvq*x1 # *x2s+vvq*x1s*x2s+2*vvq*x2s*z+vvq*x2s* # zs+vvq*x2s) dw2= dw2+rdms*rwms*rwmpgs*(-36* # vvi*x1*x2*z-42*vvi*x1*x2-6*vvi*x1*x2s- # 42*vvi*x1*z-12*vvi*x1*zs-30*vvi*x1+24* # vvi*x1s*x2+30*vvi*x1s*z+42*vvi*x1s # -18*vvi*x1c+30*vvi*x2*z+12*vvi*x2*zs # +18*vvi*x2+6*vvi*x2s*z+6*vvi*x2s+ # 12*vvi*z+6*vvi*zs+6*vvi+3*vv*x1-3*vv # *z-3*vv) * endif wcom= wcom*(2.d0-y)**2/ys-dw1+2.d0/y*(1.d0-1.d0/y)*dw2 * endif www= y*wcom/pwtm2 * 3 if(iz.eq.0) then ww(ii,jj)= 0.d0 iz= 1 else ww(ii,jj)= www*tjac(ii,jj)*tfactsw endif * enddo * 2 if(iz.eq.0) then wwt(ii)= 0.d0 iz= 1 else wwt(ii)= ww(ii,1)+ww(ii,2)+ww(ii,3)+ww(ii,4)+ww(ii,5) endif * enddo * 1 if(iz.eq.0) then wwtt= 0.d0 else wwtt= wwt(1)+wwt(2)+wwt(3) endif * if(wwtt.lt.0.d0) then iz= 0 ifz(30)= ifz(30)+1 wtoxsswc= 0.d0 else wtoxsswc= wwtt endif * return end