program threenu implicit double precision (a-z) integer nstep,n1,n2,n3,i,j,jj,off,nchange integer iask,ntheta,icount,ima,nask,kmed parameter (nstep=100000) c parameter (nvar=16,ima=1) c parameter (nvar=48,ima=3) c parameter (nvar=160,ima=10) c parameter (nvar=256,ima=16) cx parameter (nvar=480,ima=30) parameter (nvar=960,ima=60) c parameter (nvar=1920,ima=120) dimension yy(ima),feq(ima) dimension fxie(ima),fximu(ima),fxitau(ima) dimension faxie(ima),faximu(ima),faxitau(ima) dimension fun(nvar),dfun(nvar) dimension barrhoee(ima),barrhomumu(ima) dimension barrhotautau(ima),barrhonorm(ima) dimension barA(ima),barB(ima),barC(ima) dimension barD(ima),barG(ima),barH(ima) dimension rhoee(ima),rhomumu(ima),rhotautau(ima),rhonorm(ima) dimension A(ima),B(ima),C(ima) dimension D(ima),G(ima),H(ima) common/constants1/cfhx,cdm21,cdm32 common/constants2/cvehx,cvlhx,cenerhx,xw,off common/dmasses/dm21,dm32,rhonorm,barrhonorm common/angles1/sinalpha,cosalpha,singamma,cosgamma common/angles2/sinbeta,cosbeta common/ct/pi common/momenta/yy,feq,ystep1 common/meandmu/me,mmu,x1 common/rhoemuons/rhomu,rhoe external dadapt external frhoe external frhomu c self interactions off, put off=0? off=0 c change dump -> dump2 (for calculations with large Self) nchange=1 c both dm2 in units eV^2 and positive c atmospheric Delta m_32 dm32=3.1d-3 c choose one of the following c BEST FIT for LMA, dm21=4.5d-5 c BEST FIT for LOW, c dm21=1.d-7 c the third angle, beta=theta13 tanthird2=0.d0 c tanthird2=0.2d0 c tanthird2=0.02d0 c tanthird2=0.005d0 c maximum according to hep-ph/0110307 c tanthird2=0.065d0 sinbeta=dsqrt(tanthird2/(1.d0+tanthird2)) cosbeta=dsqrt(1.d0/(1.d0+tanthird2)) c put maximal mixing for atm and solar singamma=1.d0/dsqrt(2.d0) cosgamma=1.d0/dsqrt(2.d0) sinalpha=1.d0/dsqrt(2.d0) cosalpha=1.d0/dsqrt(2.d0) c sun mixing zero limit c dm21=0.d0 c sinalpha=0.d0 c cosalpha=1.d0 open(unit=30,file='asim.dat',status='old') c======= some parameters ========================= c start with dump, change to dump2 later dump=1.d4 dump2=1.d3 c muon mass in MeV mmu=1.05658d2 c electron mass in MeV me=0.511d0 pnumed=3.15137d0 mplanck=1.221d0 zeta3=1.20206d0 hub=dsqrt(8.d0*pi/3.d0) Gf=1.1664d0 xw=0.231d0 mw=80.42d0 c cvehx=8.d0*dsqrt(2.d0)*Gf*mplanck*1.d5*7.d0*pi*pi cvehx=cvehx/hub/1.8d2/mw/mw cvlhx=dsqrt(2.d0)*Gf*mplanck*1.d11 cvlhx=cvlhx/hub/dump if(off.eq.0)then cvlhx=0.d0 end if cenerhx=8.d0*dsqrt(2.d0)*Gf*mplanck*1.d5*(1.d0-xw) cenerhx=cenerhx/3.d0/hub/mw/mw cdm21=mplanck*dm21/hub/2.d0*1.d10 cdm32=mplanck*dm32/hub/2.d0*1.d10 cfhx=zeta3*mplanck*Gf*Gf*pnumed/hub/3.d0/pi/pi/pi c for pure vac c xini=0.005d0 c xmax=0.2d0 c x range c with potentials and collisions xini=0.05d0 c just solar c xini=0.12d0 c for LMA xmax=1.d0 c for LOW, if no self c xmax=4.d0 c xmax=9.d0 c LOW sun c xini=0.6d0 c xmax=5.d0 c initial degeneracies xie=0.d0 xitau=0.d0 ximu=3.d0 c ximu=0.d0 c xitau=1.d-1 c ximu=-5.d-2 c xitau=ximu c xie=-ximu c============================================================= c momentum grid (linear) cx ymax=1.2d1 ymax=2.d1 ymin=0.08d0 ystep1=(ymax-ymin)/float(ima) yy(1)=ymin do i=2,ima yy(i)=yy(i-1)+ystep1 c write(*,*)i,yy(i) end do c yy(1)=3.d0 c yy(2)=4.d0 c yy(3)=5.d0 c ystep1=1.d0 c neutrino density matrix do j=1,ima feq(j)=1.d0/(dexp(yy(j))+1.d0) fxie(j)=1.d0/(dexp(yy(j)-xie)+1.d0) fximu(j)=1.d0/(dexp(yy(j)-ximu)+1.d0) fxitau(j)=1.d0/(dexp(yy(j)-xitau)+1.d0) faxie(j)=1.d0/(dexp(yy(j)+xie)+1.d0) faximu(j)=1.d0/(dexp(yy(j)+ximu)+1.d0) faxitau(j)=1.d0/(dexp(yy(j)+xitau)+1.d0) c rhoee rhoee(j)=fxie(j)/feq(j) fun(j)=rhoee(j) c anti rhoee barrhoee(j)=faxie(j)/feq(j) fun(j+8*ima)=barrhoee(j) c rhomumu rhomumu(j)=fximu(j)/feq(j) fun(j+ima)=rhomumu(j) c anti rhomumu barrhomumu(j)=faximu(j)/feq(j) fun(j+9*ima)=barrhomumu(j) c rhotautau rhotautau(j)=fxitau(j)/feq(j) c anti rhotautau barrhotautau(j)=faxitau(j)/feq(j) rhonorm(j)=rhoee(j)+rhomumu(j)+rhotautau(j) barrhonorm(j)=barrhoee(j)+barrhomumu(j)+barrhotautau(j) c emu terms A(j)=0.d0 fun(j+2*ima)=A(j) B(j)=0.d0 fun(j+3*ima)=B(j) c anti emu terms barA(j)=0.d0 fun(j+10*ima)=barA(j) barB(j)=0.d0 fun(j+11*ima)=barB(j) c etau terms C(j)=0.d0 fun(j+4*ima)=C(j) D(j)=0.d0 fun(j+5*ima)=D(j) c anti etau terms barC(j)=0.d0 fun(j+12*ima)=barC(j) barD(j)=0.d0 fun(j+13*ima)=barD(j) c mutau terms G(j)=0.d0 fun(j+6*ima)=G(j) H(j)=0.d0 fun(j+7*ima)=H(j) c anti mutau terms barG(j)=0.d0 fun(j+14*ima)=barG(j) barH(j)=0.d0 fun(j+15*ima)=barH(j) end do c ASYMMETRIES lnu0=0.d0 do j=1,ima lnu0=lnu0+yy(j)*yy(j)*feq(j) end do lnu0=lnu0/(2.d0*pi*pi)*ystep1 c linear xstep=(xmax-xini)/float(nstep) c logarithmic logstep=(dlog10(xmax)-dlog10(xini))/float(nstep) x1=xini icount=1 c x-loop do i=1,nstep c linear x2=x1+xstep c log c x2=x1*1.d1**logstep h0=xstep/1.d2 if(nchange.eq.1)then if(x2.gt.2.5d-1)then nchange=0 cvlhx=cvlhx*dump/dump2 end if end if c calulate rho(mu+mubar) if(x1.gt.2.d-1)then rhomu=0.d0 else call dadapt(frhomu,0.d0,1.d2,1,0.d0,1.d-5,rhomu,err) rhomu=rhomu*2.d0/pi/pi end if c calulate rho(e+ebar) if(x1.gt.0.1d0)then call dadapt(frhoe,0.d0,1.d2,1,0.d0,1.d-5,rhoe,err) rhoe=rhoe*2.d0/pi/pi else rhoe=1.15145d0 end if CALL EVOL1(x1,x2,fun,h0) if(icount.eq.40)then do jj=1,ima rhoee(jj)=fun(jj) rhomumu(jj)=fun(jj+ima) rhotautau(jj)=rhonorm(jj)-rhoee(jj)-rhomumu(jj) barrhoee(jj)=fun(jj+8*ima) barrhomumu(jj)=fun(jj+9*ima) barrhotautau(jj)=barrhonorm(jj)-barrhoee(jj)-barrhomumu(jj) end do lnue=0.d0 lnumu=0.d0 lnutau=0.d0 do j=1,ima lnue=lnue+yy(j)*yy(j)*feq(j)*(rhoee(j)-barrhoee(j)) lnumu=lnumu+yy(j)*yy(j)*feq(j)*(rhomumu(j)-barrhomumu(j)) lnutau=lnutau+yy(j)*yy(j)*feq(j)*(rhotautau(j)-barrhotautau(j)) end do lnue=lnue/(2.d0*pi*pi)*ystep1 lnumu=lnumu/(2.d0*pi*pi)*ystep1 lnutau=lnutau/(2.d0*pi*pi)*ystep1 c calculate the value of degeneracies c better calculation Amu0=(-2.7d1*lnumu+ # dsqrt(3.d0*(2.43d2*lnumu*lnumu+pi*pi)))**(1.d0/3.d0) Amu1=(pi**4/3.d0)**(1.d0/3.d0) Amu2=(pi/3.d0)**(2.d0/3.d0) ximux=Amu1/Amu0-Amu2*Amu0 Ae0=(-2.7d1*lnue+ # dsqrt(3.d0*(2.43d2*lnue*lnue+pi*pi)))**(1.d0/3.d0) xiex=Amu1/Ae0-Amu2*Ae0 Atau0=(-2.7d1*lnutau+ # dsqrt(3.d0*(2.43d2*lnutau*lnutau+pi*pi)))**(1.d0/3.d0) xitaux=Amu1/Atau0-Amu2*Atau0 c approx for small ximu / xie /xitau c ximux=6.d0*lnumu c xiex=6.d0*lnue c xitaux=6.d0*lnutau write(*,'(5e15.7)')x2,ximux,xiex,xitaux c write(10,'(5e15.7)')x2,ximux,xiex,xitaux write(30,'(5e15.7)')x2,ximux,xiex,xitaux c write(*,'(5e15.7)')x2,lnumu/lnu0,lnue/lnu0,lnutau/lnu0 c write(30,'(5e15.7)')x2,lnumu/lnu0,lnue/lnu0,lnutau/lnu0 c write(40,'(5e15.7)')x2,rhoe,rhomu icount=0 end if x1=x2 icount=icount+1 end do close(30) c close(40) stop end c============================================================= subroutine evol1(x1,x2,fun,h0) implicit double precision (a-z) integer nstep,nvar,ifail,i,j,ima c parameter (nvar=16,ima=1) c parameter (nvar=48,ima=3) c parameter (nvar=160,ima=10) c parameter (nvar=256,ima=16) cx parameter (nvar=480,ima=30) parameter (nvar=960,ima=60) c parameter (nvar=1920,ima=120) c FOR CERN cx dimension fun(nvar),w(10*nvar) c FOR NAG dimension fun(nvar),w(28+21*nvar) external d02cjf,d02cjx,d02cjw external dfundx ifail=0 c eps=1.d-3 eps=1.d-6 cx CALL DDEQMR(nvar,x1,x2,fun,h0,eps,dfundx,w) call d02cjf(x1,x2,nvar,fun,dfundx,eps,'Default', # d02cjx,d02cjw,w,ifail) return end c============================================================= subroutine dfundx(x,fun,dfun) implicit double precision (a-z) integer nvar,i,j,jj,k,itest,ima,off c parameter (nvar=16,ima=1) c parameter (nvar=48,ima=3) c parameter (nvar=160,ima=10) c parameter (nvar=256,ima=16) cx parameter (nvar=480,ima=30) parameter (nvar=960,ima=60) c parameter (nvar=1920,ima=120) dimension fun(nvar),dfun(nvar) dimension yy(ima),feq(ima) dimension barrhoee(ima),barrhomumu(ima) dimension barrhotautau(ima),barrhonorm(ima) dimension barA(ima),barB(ima),barC(ima) dimension barD(ima),barG(ima),barH(ima) dimension rhoee(ima),rhomumu(ima),rhotautau(ima),rhonorm(ima) dimension A(ima),B(ima),C(ima) dimension D(ima),G(ima),H(ima) common/constants1/cfhx,cdm21,cdm32 common/constants2/cvehx,cvlhx,cenerhx,xw,off common/dmasses/dm21,dm32,rhonorm,barrhonorm common/angles1/sina,cosa,sing,cosg common/angles2/sinb,cosb common/ct/pi common/momenta/yy,feq,ystep1 common/rhoemuons/rhomu,rhoe xx=x*x x4=x*x*x*x drho=dsqrt(pi*pi/3.d1*7.25d0+rhoe+rhomu) c now there are 2 vacuum oscillation terms dm21phx=cdm21*xx/drho dm32phx=cdm32*xx/drho c rho(electron+positron) background vehx=cvehx/drho/x4*(rhoe/1.15145d0) c rho(muons+antimuons) background vmuhx=cvehx/drho/x4*(rhomu/1.15145d0) c (n_nu-n_antinu) background vlhx=cvlhx/drho/xx c rho(n_nu+n_antinu) background c vUhx=cenerhx/drho/x4 vUhx=0.d0 c collisions, first approximation cGamma=cfhx/drho/x4 c old values c dhx=12.43d0*cGamma c chx=1.21d0*cGamma c mckellar-thomson 1994 for numu-nutau c dmutauhx=(3.2d1*xw*xw-1.6d1*xw+1.6d1)*cGamma c cmutauhx=(1.6d1*xw*xw-8.d0*xw+1.4d1)*cGamma c our values for nue-nutau dmutauhx=(8.d0*xw*xw-4.d0*xw+4.d0)*cGamma c cmutauhx=-dmutauhx cmutauhx=0.d0 c new values emu or etau demuhx=12.21d0*cGamma c cemuhx=-2.42d0*cGamma cemuhx=0.d0 c write(*,*)'xw=',xw c stop c old values c demuhx=12.43d0*cGamma c cemuhx=1.21d0*cGamma c============================================================= Vee=0.d0 Vmm=0.d0 Vtt=0.d0 Remu=0.d0 Yemu=0.d0 Retau=0.d0 Yetau=0.d0 Rmutau=0.d0 Ymutau=0.d0 do jj=1,ima rhoee(jj)=fun(jj) rhomumu(jj)=fun(jj+ima) rhotautau(jj)=rhonorm(jj)-rhoee(jj)-rhomumu(jj) barrhoee(jj)=fun(jj+8*ima) barrhomumu(jj)=fun(jj+9*ima) barrhotautau(jj)=barrhonorm(jj)-barrhoee(jj) # -barrhomumu(jj) A(jj)=fun(jj+2*ima) B(jj)=fun(jj+3*ima) C(jj)=fun(jj+4*ima) D(jj)=fun(jj+5*ima) G(jj)=fun(jj+6*ima) H(jj)=fun(jj+7*ima) barA(jj)=fun(jj+10*ima) barB(jj)=fun(jj+11*ima) barC(jj)=fun(jj+12*ima) barD(jj)=fun(jj+13*ima) barG(jj)=fun(jj+14*ima) barH(jj)=fun(jj+15*ima) Vee=Vee+yy(jj)*yy(jj)*feq(jj)*(rhoee(jj)-barrhoee(jj)) Vmm=Vmm+yy(jj)*yy(jj)*feq(jj) # *(rhomumu(jj)-barrhomumu(jj)) Vtt=Vtt+yy(jj)*yy(jj)*feq(jj) # *(rhotautau(jj)-barrhotautau(jj)) Remu=Remu+yy(jj)*yy(jj)*feq(jj)*(A(jj)-barA(jj)) Yemu=Yemu+yy(jj)*yy(jj)*feq(jj)*(B(jj)-barB(jj)) Retau=Retau+yy(jj)*yy(jj)*feq(jj)*(C(jj)-barC(jj)) Yetau=Yetau+yy(jj)*yy(jj)*feq(jj)*(D(jj)-barD(jj)) Rmutau=Rmutau+yy(jj)*yy(jj)*feq(jj)*(G(jj)-barG(jj)) Ymutau=Ymutau+yy(jj)*yy(jj)*feq(jj)*(H(jj)-barH(jj)) c write(*,*)jj,Px(jj),Py(jj),Pz(jj) end do Vee=Vee/(2.d0*pi*pi)*ystep1 Vmm=Vmm/(2.d0*pi*pi)*ystep1 Vtt=Vtt/(2.d0*pi*pi)*ystep1 Remu=Remu/(2.d0*pi*pi)*ystep1 Yemu=Yemu/(2.d0*pi*pi)*ystep1 Retau=Retau/(2.d0*pi*pi)*ystep1 Yetau=Yetau/(2.d0*pi*pi)*ystep1 Rmutau=Rmutau/(2.d0*pi*pi)*ystep1 Ymutau=Ymutau/(2.d0*pi*pi)*ystep1 c write(*,*)Vee c write(*,*)Vmm c write(*,*)Vtt c write(*,*)Remu c write(*,*)Yemu c write(*,*)Retau c write(*,*)Yetau c write(*,*)Rmutau c write(*,*)Ymutau c write(*,*) c stop *************************************************************************** * evolution equations *************************************************************************** * dfun(i)= dfun/dx (i) *************************************************************************** * anyadir /yy(i) !!!! do i=1,ima c drhoee/dx keeB=-2.d0*cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i) 1 *sinb*sing + dm21phx/yy(i)*cosa**2*sinb*sing) keeD=-2.d0*cosb*(dm32phx/yy(i)*cosg*sinb + dm21phx/yy(i)*cosa**2 1 *cosg*sinb - dm21phx/yy(i)*cosa*sina*sing) dfun(i)=keeB*B(i)+keeD*D(i) 1 +2.d0*vlhx*(-Remu*B(i)-Retau*D(i)+Yemu*A(i)+Yetau*C(i)) c dbarrhoee/dx dfun(i+8*ima)=-keeB*barB(i)-keeD*barD(i) 1 +2.d0*vlhx*(-Remu*barB(i)-Retau*barD(i)+Yemu*barA(i)+Yetau*barC(i)) c drhomumu/dx kmumuB=2.d0*(dm21phx/yy(i)*cosa*cosb *cosg*sina + dm32phx/yy(i) 1 *cosb*sinb*sing + dm21phx/yy(i)*cosa**2*cosb *sinb*sing) kmumuH=2.d0*(dm21phx/yy(i)*cosa*(cosg**2-sing**2)*sina*sinb - 1 dm32phx/yy(i) *cosb**2*cosg*sing - dm21phx/yy(i)*cosg 1 *sina**2*sing + dm21phx/yy(i)*cosa**2 *cosg*sinb**2 2 *sing) dfun(i+ima)=kmumuB*B(i)+kmumuH*H(i) 1 +2.d0*vlhx*(Remu*B(i)-Rmutau*H(i)-Yemu*A(i)+Ymutau*G(i)) c dbarrhomumu/dx dfun(i+9*ima)=-kmumuB*barB(i)-kmumuH*barH(i) 1 +2.d0*vlhx*(Remu*barB(i)-Rmutau*barH(i)-Yemu*barA(i)+Ymutau*barG(i)) c dA/dx kaB=dm21phx/yy(i)*cosg**2*sina**2 + 2.d0*dm21phx/yy(i)*cosa*cosg 1 *sina*sinb*sing + dm32phx/yy(i)*(sinb**2 - cosb**2*sing**2) 1 + dm21phx/yy(i)*cosa**2*(-cosb**2 + sinb**2*sing**2) kaD=dm21phx/yy(i)*cosa*(cosg**2-sing**2)*sina*sinb - cosg 1 *(dm32phx/yy(i)*cosb**2 + dm21phx/yy(i)*sina**2)*sing 1 + dm21phx/yy(i)*cosa**2*cosg*sinb**2*sing kaH=-(cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg*sinb 1 - dm21phx/yy(i)*cosa*sina*sing)) dfun(i+2*ima)=kaB*B(i) + kaD*D(i) + kaH*H(i) 1 - (vehx-vmuhx)*yy(i)*B(i) 2 - demuhx*yy(i)*A(i) 3 - cemuhx*yy(i)*barA(i) 4 +vlhx*((Vee-Vmm)*B(i)-Retau*H(i)-Rmutau*D(i)+Yetau*G(i) 5 +Ymutau*C(i)-Yemu*(rhoee(i)-rhomumu(i))) c dbarA/dx dfun(i+10*ima)=-kaB*barB(i) - kaD*barD(i) - kaH*barH(i) 1 + (vehx-vmuhx)*yy(i)*barB(i) 2 - demuhx*yy(i)*barA(i) 3 - cemuhx*yy(i)*A(i) 4 +vlhx*((Vee-Vmm)*barB(i)-Retau*barH(i)-Rmutau*barD(i) 5 +Yetau*barG(i)+Ymutau*barC(i)-Yemu*(barrhoee(i)-barrhomumu(i))) c dB/dx kbA=-(dm21phx/yy(i)*cosg**2*sina**2) - 2.d0*dm21phx/yy(i)*cosa 1 *cosg*sina*sinb*sing - dm32phx/yy(i)*(sinb**2 - cosb**2 1 *sing**2) - dm21phx/yy(i)*cosa**2*(-cosb**2 + sinb**2 2 *sing**2) kbC=-(dm21phx/yy(i)*cosa*(cosg**2-sing**2)*sina*sinb) + cosg 1 *(dm32phx/yy(i)*cosb**2 + dm21phx/yy(i)*sina**2)*sing 1 - dm21phx/yy(i)*cosa**2*cosg*sinb**2*sing kbG=-(cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg*sinb 1 - dm21phx/yy(i)*cosa*sina*sing)) kbee=cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i)*sinb 1 *sing + dm21phx/yy(i)*cosa**2*sinb*sing) kbmumu=-(cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i) 1 *sinb*sing + dm21phx/yy(i)*cosa**2*sinb*sing)) dfun(i+3*ima)= kbA*A(i)+kbC*C(i)+kbG*G(i)+kbee*rhoee(i) 1 +kbmumu*rhomumu(i) + (vehx-vmuhx)*yy(i)*A(i) 2 - demuhx*yy(i)*B(i) 3 + cemuhx*yy(i)*barB(i) 4 +vlhx*(-(Vee-Vmm)*A(i)-Retau*G(i)+Rmutau*C(i)-Yetau*H(i) 5 +Ymutau*D(i)+Remu*(rhoee(i)-rhomumu(i))) c dbarB/dx dfun(i+11*ima)=-kbA*barA(i)-kbC*barC(i)-kbG*barG(i)-kbee*barrhoee(i) 1 -kbmumu*barrhomumu(i) - (vehx-vmuhx)*yy(i)*barA(i) 2 -demuhx*yy(i)*barB(i) 3 + cemuhx*yy(i)*B(i) 4 +vlhx*(-(Vee-Vmm)*barA(i)-Retau*barG(i)+Rmutau*barC(i) 5 -Yetau*barH(i)+Ymutau*barD(i)+Remu*(barrhoee(i)-barrhomumu(i))) c dC/dx kcB=dm21phx/yy(i)*cosa*(cosg**2-sing**2)*sina*sinb - cosg 1 *(dm32phx/yy(i)*cosb**2 + dm21phx/yy(i)*sina**2)*sing 1 + dm21phx/yy(i)*cosa**2*cosg*sinb**2*sing kcD=-(dm32phx/yy(i)*cosb**2*cosg**2) + dm32phx/yy(i)*sinb**2 1 + dm21phx/yy(i)*cosa**2*(-cosb**2 + cosg**2*sinb**2) 1 - 2.d0*dm21phx/yy(i)*cosa*cosg*sina*sinb*sing 2 + dm21phx/yy(i)*sina**2*sing**2 kcH=cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i)*sinb*sing 1 + dm21phx/yy(i)*cosa**2*sinb*sing) dfun(i+4*ima)=kcB*B(i)+kcD*D(i)+kcH*H(i) 1 - vehx*yy(i)*D(i) 2 - demuhx*yy(i)*C(i) 3 - cemuhx*yy(i)*barC(i) 4 +vlhx*((Vee-Vtt)*D(i)+Remu*H(i)-Rmutau*B(i)+Yemu*G(i) 5 -Ymutau*A(i)-Yetau*(rhoee(i)-rhotautau(i))) c dbarC/dx dfun(i+12*ima)=-kcB*barB(i)-kcD*barD(i)-kcH*barH(i) 1 + vehx*yy(i)*barD(i) 2 - demuhx*yy(i)*barC(i) 3 - cemuhx*yy(i)*C(i) 4 +vlhx*((Vee-Vtt)*barD(i)+Remu*barH(i)-Rmutau*barB(i) 5 +Yemu*barG(i)-Ymutau*barA(i)-Yetau*(barrhoee(i)-barrhotautau(i))) c dD/dx kdA=-(dm21phx/yy(i)*cosa*(cosg**2-sing**2)*sina*sinb) + cosg 1 *(dm32phx/yy(i)*cosb**2 + dm21phx/yy(i)*sina**2)*sing 1 - dm21phx/yy(i)*cosa**2*cosg*sinb**2*sing kdC=dm32phx/yy(i)*cosb**2*cosg**2 - dm32phx/yy(i)*sinb**2 1 + dm21phx/yy(i)*cosa**2*(cosb**2 - cosg**2*sinb**2) 1 + 2.d0*dm21phx/yy(i)*cosa*cosg*sina*sinb*sing 2 - dm21phx/yy(i)*sina**2*sing**2 kdG=-(cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i)*sinb 1 *sing + dm21phx/yy(i)*cosa**2*sinb*sing)) kdee=cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg*sinb 1 - dm21phx/yy(i)*cosa*sina*sing) kdtautau=-(cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg 1 *sinb - dm21phx/yy(i)*cosa*sina*sing)) dfun(i+5*ima)= kdA*A(i)+kdC*C(i)+kdG*G(i)+kdee*rhoee(i) 1 +kdtautau*rhotautau(i) + vehx*yy(i)*C(i) 2 - demuhx*yy(i)*D(i) 3 + cemuhx*yy(i)*barD(i) 4 +vlhx*(-(Vee-Vtt)*C(i)-Remu*G(i)+Rmutau*A(i)+Yemu*H(i) 5 -Ymutau*B(i)+Retau*(rhoee(i)-rhotautau(i))) c dbarD/dx dfun(i+13*ima)=-kdA*barA(i)-kdC*barC(i)-kdG*barG(i)-kdee*barrhoee(i) 1 -kdtautau*barrhotautau(i) - vehx*yy(i)*barC(i) 2 - demuhx*yy(i)*barD(i) 3 + cemuhx*yy(i)*D(i) 4 +vlhx*(-(Vee-Vtt)*barC(i)-Remu*barG(i)+Rmutau*barA(i) 5 +Yemu*barH(i)-Ymutau*barB(i)+Retau*(barrhoee(i)-barrhotautau(i))) c dG/dx kgB=cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg*sinb 1 - dm21phx/yy(i)*cosa*sina*sing) kgD=cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i)*sinb 1 *sing + dm21phx/yy(i)*cosa**2*sinb*sing) kgH=-(dm32phx/yy(i)*cosb**2*(cosg**2-sing**2)) + dm21phx/yy(i) 1 *(cosg**2*(-sina**2 + cosa**2*sinb**2) - 4.d0*cosa*cosg 1 *sina*sinb*sing + (sina**2 - cosa**2*sinb**2)*sing**2) dfun(i+6*ima)= kgB*B(i)+kgD*D(i)+kgH*H(i) 1 - vmuhx*yy(i)*H(i) 2 - dmutauhx*yy(i)*G(i) 3 - cmutauhx*yy(i)*barG(i) 4 +vlhx*((Vmm-Vtt)*H(i)+Remu*D(i)+Retau*B(i)-Yemu*C(i) 5 -Yetau*A(i)-Ymutau*(rhomumu(i)-rhotautau(i))) c dbarG/dx dfun(i+14*ima)=-kgB*barB(i)-kgD*barD(i)-kgH*barH(i) 1 + vmuhx*yy(i)*barH(i) 2 - dmutauhx*yy(i)*barG(i) 3 - cmutauhx*yy(i)*G(i) 4 +vlhx*((Vmm-Vtt)*barH(i)+Remu*barD(i)+Retau*barB(i) 5 -Yemu*barC(i)-Yetau*barA(i)-Ymutau*(barrhomumu(i)-barrhotautau(i))) c dH/dx khA=cosb*((dm32phx/yy(i) + dm21phx/yy(i)*cosa**2)*cosg*sinb 1 - dm21phx/yy(i)*cosa*sina*sing) khC=-(cosb*(dm21phx/yy(i)*cosa*cosg*sina + dm32phx/yy(i)*sinb 1 *sing + dm21phx/yy(i)*cosa**2*sinb*sing)) khG=dm32phx/yy(i)*cosb**2*(cosg**2-sing**2) + dm21phx/yy(i) 1 *(cosg**2*(sina**2 - cosa**2*sinb**2) - sina**2*sing**2 1 + cosa**2*sinb**2*sing**2 + 4.d0*sina*cosa*sinb 2 *sing*cosg) khmumu=dm32phx/yy(i)*cosb**2*cosg*sing - dm21phx/yy(i)*(cosa 1 *cosg*sinb - sina*sing)*(cosg*sina + cosa*sinb*sing) khtautau=-(dm32phx/yy(i)*cosb**2*cosg*sing) + dm21phx/yy(i) 1 *(cosa*cosg*sinb - sina*sing)*(cosg*sina + cosa*sinb 1 *sing) dfun(i+7*ima)= khA*A(i)+khC*C(i)+khG*G(i)+khmumu*rhomumu(i) 1 +khtautau*rhotautau(i) + vmuhx*yy(i)*G(i) 2 - dmutauhx*yy(i)*H(i) 3 + cmutauhx*yy(i)*barH(i) 4 +vlhx*(-(Vmm-Vtt)*G(i)-Remu*C(i)+Retau*A(i)-Yemu*D(i) 5 +Yetau*B(i)+Rmutau*(rhomumu(i)-rhotautau(i))) c dbarH/dx dfun(i+15*ima)=-khA*barA(i)-khC*barC(i)-khG*barG(i) 1 -khmumu*barrhomumu(i) 1 -khtautau*barrhotautau(i) - vmuhx*yy(i)*barG(i) 2 - dmutauhx*yy(i)*barH(i) 3 + cmutauhx*yy(i)*H(i) 4 +vlhx*(-(Vmm-Vtt)*barG(i)-Remu*barC(i)+Retau*barA(i) 5 -Yemu*barD(i)+Yetau*barB(i)+Rmutau*(barrhomumu(i)-barrhotautau(i))) c write(*,*)i c write(*,*)A(i) c write(*,*)B(i) c write(*,*)C(i) c write(*,*)D(i) c write(*,*)barA(i) c write(*,*)barB(i) c write(*,*)barC(i) c write(*,*)barD(i) c write(*,*)rhoee(i) c write(*,*)barrhoee(i) c write(*,*) c write(*,*)i c write(*,*)dfun(i+6*ima) c write(*,*)dfun(i+7*ima) c write(*,*)dfun(i+14*ima) c write(*,*)dfun(i+15*ima) c write(*,*)dfun(i) c write(*,*)dfun(i+ima) c write(*,*)dfun(i+8*ima) c write(*,*)dfun(i+9*ima) c write(*,*) end do c stop return end c============================================================= double precision function frhomu(y) implicit double precision (a-z) integer nstep,nvar,i,j,k common/meandmu/me,mmu,x emu=dsqrt(y*y+mmu*mmu*x*x) frhomu=y*y*emu/(dexp(emu)+1.d0) return end double precision function frhoe(y) implicit double precision (a-z) integer nstep,nvar,i,j,k common/meandmu/me,mmu,x ee=dsqrt(y*y+me*me*x*x) frhoe=y*y*ee/(dexp(ee)+1.d0) return end block data implicit double precision (a-z) common/ct/pi C Valore di pi greco. data pi/3.141592653589793d0/ end