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=5000)
c	parameter (nvar=16,ima=1)
c	parameter (nvar=80,ima=5)
c	parameter (nvar=480,ima=30)
	parameter (nvar=960,ima=60)
c	parameter (nvar=1920,ima=120)
	dimension pnu(ima)
	dimension fnue(ima),fnumu(ima),fnutau(ima)
	dimension fbarnue(ima),fbarnumu(ima),fbarnutau(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/momenta/pnu,ystep1
	common/density/densx,yex
	common/constants1/cdm21,cdm32,denscoef,neubcoef
	common/dmasses/dm21,dm32,rhonorm,barrhonorm
	common/angles1/sinalpha,cosalpha,singamma,cosgamma
	common/angles2/sinbeta,cosbeta
	common/ct/pi

c       both dm2 in units eV^2 and positive

c       atmospheric Delta m_32
        dm32=-3.1d-3

c	BEST FIT for ATM, 
	tanatm2=1.4d0

c       choose one of the following
c	BEST FIT for LMA, 
	dm21=4.5d-5
	tansolar2=0.41d0

c	BEST FIT for LOW, 
c	dm21=1.d-7
c	tansolar2=0.71d0

c	tansolar2=1.d-2

c	alpha=theta12, gamma=theta23

	sinalpha=dsqrt(tansolar2/(1.d0+tansolar2))
	cosalpha=dsqrt(1.d0/(1.d0+tansolar2))
	singamma=dsqrt(tanatm2/(1.d0+tanatm2))
	cosgamma=dsqrt(1.d0/(1.d0+tanatm2))

c       the third angle, beta=theta13

	tanthird2=0.d0
c       maximum according to hep-ph/0110307
c	tanthird2=0.065d0
	tanthird2=1.d-4
	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)
c	sinalpha=1.d0/dsqrt(2.d0)
c	cosalpha=1.d0/dsqrt(2.d0)
c	singamma=0.d0
c	cosgamma=1.d0
c	sinalpha=1.d-2
c	cosalpha=dsqrt(1.d0-sinalpha*sinalpha)

c       check with 2nu program here XXXXXXXXX !!!!
c	s2tsq=0.9722d0
c	singamma=0.6455d0
c	cosgamma=dsqrt(1.d0-singamma*singamma)
c	sinalpha=0.d0
c	cosalpha=1.d0
c	dm21=0.d0

 	open(unit=30,file='evol.dat',status='old')
c 	open(unit=40,file='fnu.dat',status='old')


c======= some parameters =========================

	pmed=3.15137d0
c	
c       neutrino sphere radius in cm
	Rnu=1.1d6


c       neutrino luminosity in each flavour in units 10^51 erg/sec
c       put zero if no neutrino background effect
c	Lneu=1.d4
	Lneu=1.d0
	off=0.d0
	
c       some constants (here Rnu in 10km!)
	cdm21=5.068d4*(Rnu/1.d6)*dm21/2.d0
	cdm32=5.068d4*(Rnu/1.d6)*dm32/2.d0
	denscoef=-2.548d4*(Rnu/1.d6)*dsqrt(2.d0)
c	denscoef=0.d0
	neubcoef=3.7d7*Lneu/(Rnu/1.d6)/(2.d0*pi*pi)

	rmin=1.d9/Rnu
	rmax=5.d10/Rnu

c       initial degeneracies (?)
	xie=0.d0
	xitau=0.d0
	ximu=-1.d-1

c=============================================================

c       fix average energies
	Enue=1.1d1
c	Ebarnue=1.6d1
	Ebarnue=2.5d1
	Enumu=2.5d1
	Ebarnumu=Enumu

c       evaluate FD neutrino distributions
	pmed=3.15137d0
	Tnue=Enue/pmed
	Tbarnue=Ebarnue/pmed
	Tnumu=Enumu/pmed
	Tbarnumu=Ebarnumu/pmed

c       momentum grid (linear)

c       set #2
	Emin=0.5d0
	Emax=1.d2
	ystep1=(Emax-Emin)/float(ima-1)
	
	pnu(1)=Emin
c	pnu(1)=1.5d1

	do i=2,ima
	   pnu(i)=pnu(i-1)+ystep1
c	   write(*,*)i,pnu(i)
	end do
c	stop

 13	 continue

c	neutrino density matrix

	 flux=0.5d0
	do j=1,ima

	   fnue(j)=pnu(j)**2/Tnue**4/(dexp(pnu(j)/Tnue)+1.d0)
	   fnumu(j)=flux*pnu(j)**2/Tnumu**4/(dexp(pnu(j)/Tnumu)+1.d0)
	   fnutau(j)=fnumu(j)
	   fbarnue(j)=pnu(j)**2/Tbarnue**4/(dexp(pnu(j)/Tbarnue)+1.d0)
	   fbarnumu(j)=flux*pnu(j)**2/Tbarnumu**4/(dexp(pnu(j)/Tbarnumu)+1.d0)
	   fbarnutau(j)=fbarnumu(j)

c       rhoee
	   rhoee(j)=fnue(j)
	   fun(j)=rhoee(j)

c       anti rhoee
	   barrhoee(j)=fbarnue(j)
	   fun(j+8*ima)=barrhoee(j)

c       rhomumu
	   rhomumu(j)=fnumu(j)
	   fun(j+ima)=rhomumu(j)

c       anti rhomumu
	   barrhomumu(j)=fbarnumu(j)
	   fun(j+9*ima)=barrhomumu(j)

c       rhotautau
	   rhotautau(j)=fnutau(j)

c       anti rhotautau
	   barrhotautau(j)=fbarnutau(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)


c	   write(50,'(5e15.7)')pnu(j),rhoee(j),rhomumu(j),rhotautau(j)
c	write(50,'(5e15.7)')pnu(j),barrhoee(j),barrhomumu(j),barrhotautau(j)
	end do

 111	continue

c       linear
	rstep=(rmax-rmin)/float(nstep-1)

c       logarithmic
	logstep=(dlog10(rmax)-dlog10(rmin))/float(nstep-1)

	r1=rmin

	icount=1

c       x-loop

	do i=1,nstep

c       linear
c	   r2=r1+rstep
c       log
	   r2=r1*1.d1**logstep

	   h0=rstep/1.d2

c       electron background
	   yex=0.5d0
	   densx=2.d-3*(1.d9/r2/Rnu)**3/yex

c	   write(50,*)r2*Rnu,yex*densx*1.d7

	   CALL EVOL1(r1,r2,fun,h0)


		 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

	      if(icount.eq.10)then

	      open(unit=20,file='fnu.dat',status='old')
	      open(unit=70,file='afnu.dat',status='old')
	      do n1=1,ima
		 fae=rhoee(n1)
		 famu=rhomumu(n1)
		 fatau=rhotautau(n1)
		 faantie=barrhoee(n1)
		 faantimu=barrhomumu(n1)
		 faantitau=barrhotautau(n1)
	    write(20,'(5e15.7)')pnu(n1),fae,famu,fatau
	    write(70,'(5e15.7)')pnu(n1),faantie,faantimu,faantitau
	      end do

	      close(20)
	      close(70)

	      icount=0
	   end if

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

	      write(*,'(5e15.7)')r2*Rnu,rhoee(1)/fnue(1)
c	      write(30,'(5e15.7)')r2,rhoee(2)/fnue(2),rhoee(3)/fnue(3)
	      write(30,'(5e15.7)')r2,rhoee(10)/fnue(10),rhoee(20)/fnue(20),rhoee(25)/fnue(25)
	      write(50,'(5e15.7)')r2,rhomumu(10)/fnumu(20),rhomumu(20)/fnumu(20),
     #          rhomumu(25)/fnumu(25)
	      write(60,'(5e15.7)')r2,rhotautau(10)/fnutau(20),rhotautau(20)/fnutau(20),
     #          rhotautau(25)/fnutau(25)
c	      write(30,'(5e15.7)')r2*Rnu,rhoee(1)/fnue(1),rhomumu(1)/fnumu(1),
c     #            rhotautau(1)/fnutau(1)
c     #            barrhoee(1)/fbarnue(1),barrhomumu(1)/fbarnumu(1)

c	      write(*,'(5e15.7)')r2,lnumu/lnu0,lnue/lnu0,lnutau/lnu0
c	      write(30,'(5e15.7)')r2,lnumu/lnu0,lnue/lnu0,lnutau/lnu0

c	      write(40,'(5e15.7)')r2,rhoe,rhomu

	   r1=r2

	   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=80,ima=5)
c	parameter (nvar=480,ima=30)
	parameter (nvar=960,ima=60)
c	parameter (nvar=1920,ima=120)

c       FOR CERN
c	dimension fun(nvar),w(10*nvar)

c       FOR NAG
	dimension fun(nvar),w(28+21*nvar)

	external d02cjf,d02cjx,d02cjw
	external dfundx

	ifail=0
	eps=1.d-6

c        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(rx,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=80,ima=5)
c	parameter (nvar=480,ima=30)
	parameter (nvar=960,ima=60)
c	parameter (nvar=1920,ima=120)

	dimension fun(nvar),dfun(nvar)
	dimension pnu(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/momenta/pnu,ystep1
	common/density/densx,yex
	common/constants1/cdm21,cdm32,denscoef,neubcoef
	common/dmasses/dm21,dm32,rhonorm,barrhonorm
	common/angles1/sina,cosa,sing,cosg
	common/angles2/sinb,cosb
	common/ct/pi

	yehere=yex

c       now there are 2 vacuum oscillation terms

	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)

c	   write(*,*)jj,Px(jj),Py(jj),Pz(jj)

	end do

cx	   write(*,*)cdm32,cdm21
cx	   write(*,*)cosa,cosb,cosg
cx	   write(*,*)sina,sinb,sing
cx	   write(*,*)denscoef*yehere*densx
cx	   stop

***************************************************************************
*                   evolution equations
***************************************************************************
*                  dfun(i)= dfun/dx (i)
***************************************************************************


	do i=1,ima

c       drhoee/dr

	   keeB=-2.d0*cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)
	1	*sinb*sing + cdm21/pnu(i)*cosa**2*sinb*sing)

	   keeD=-2.d0*cosb*(cdm32/pnu(i)*cosg*sinb + cdm21/pnu(i)*cosa**2
	1	*cosg*sinb - cdm21/pnu(i)*cosa*sina*sing)

	   dfun(i)=keeB*B(i)+keeD*D(i)
c	1	+2.d0*vlhx*(-Remu*B(i)-Retau*D(i)+Yemu*A(i)+Yetau*C(i))

c       dbarrhoee/dr

	   dfun(i+8*ima)=-keeB*barB(i)-keeD*barD(i)
c	1  +2.d0*vlhx*(-Remu*barB(i)-Retau*barD(i)+Yemu*barA(i)+Yetau*barC(i))

c       drhomumu/dr

	   kmumuB=2.d0*(cdm21/pnu(i)*cosa*cosb *cosg*sina + cdm32/pnu(i)
	1	*cosb*sinb*sing + cdm21/pnu(i)*cosa**2*cosb *sinb*sing)

	   kmumuH=2.d0*(cdm21/pnu(i)*cosa*(cosg**2-sing**2)*sina*sinb -
	1	cdm32/pnu(i) *cosb**2*cosg*sing - cdm21/pnu(i)*cosg
	1	*sina**2*sing + cdm21/pnu(i)*cosa**2 *cosg*sinb**2
	2	*sing)

	   dfun(i+ima)=kmumuB*B(i)+kmumuH*H(i)
c	1	+2.d0*vlhx*(Remu*B(i)-Rmutau*H(i)-Yemu*A(i)+Ymutau*G(i))

c       dbarrhomumu/dr

	   dfun(i+9*ima)=-kmumuB*barB(i)-kmumuH*barH(i)
c	1 +2.d0*vlhx*(Remu*barB(i)-Rmutau*barH(i)-Yemu*barA(i)+Ymutau*barG(i))

c       dA/dr

	   kaB=cdm21/pnu(i)*cosg**2*sina**2 + 2.d0*cdm21/pnu(i)*cosa*cosg
	1	*sina*sinb*sing + cdm32/pnu(i)*(sinb**2 - cosb**2*sing**2)
	1	+ cdm21/pnu(i)*cosa**2*(-cosb**2 + sinb**2*sing**2)

	   kaD=cdm21/pnu(i)*cosa*(cosg**2-sing**2)*sina*sinb - cosg
	1	*(cdm32/pnu(i)*cosb**2 + cdm21/pnu(i)*sina**2)*sing
	1	+ cdm21/pnu(i)*cosa**2*cosg*sinb**2*sing

	   kaH=-(cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg*sinb
	1	- cdm21/pnu(i)*cosa*sina*sing))

	   dfun(i+2*ima)=kaB*B(i) + kaD*D(i) + kaH*H(i)
	1	- denscoef*yehere*densx*B(i)
c	4       +vlhx*((Vee-Vmm)*B(i)-Retau*H(i)-Rmutau*D(i)+Yetau*G(i)
c	5       +Ymutau*C(i)-Yemu*(rhoee(i)-rhomumu(i)))

c       dbarA/dr

	   dfun(i+10*ima)=-kaB*barB(i) - kaD*barD(i) - kaH*barH(i)
	1	- denscoef*yehere*densx*barB(i)
c	4       +vlhx*((Vee-Vmm)*barB(i)-Retau*barH(i)-Rmutau*barD(i)
c	5       +Yetau*barG(i)+Ymutau*barC(i)-Yemu*(barrhoee(i)-barrhomumu(i)))

c       dB/dr
	   
	   kbA=-(cdm21/pnu(i)*cosg**2*sina**2) - 2.d0*cdm21/pnu(i)*cosa
	1	*cosg*sina*sinb*sing - cdm32/pnu(i)*(sinb**2 - cosb**2
	1	*sing**2) - cdm21/pnu(i)*cosa**2*(-cosb**2 + sinb**2
	2	*sing**2)

	   kbC=-(cdm21/pnu(i)*cosa*(cosg**2-sing**2)*sina*sinb) + cosg
	1	*(cdm32/pnu(i)*cosb**2 + cdm21/pnu(i)*sina**2)*sing
	1	- cdm21/pnu(i)*cosa**2*cosg*sinb**2*sing

	   kbG=-(cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg*sinb
	1	- cdm21/pnu(i)*cosa*sina*sing))

	   kbee=cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)*sinb
	1	*sing + cdm21/pnu(i)*cosa**2*sinb*sing)

	   kbmumu=-(cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)
	1	*sinb*sing + cdm21/pnu(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) + denscoef*yehere*densx*A(i)
c	4       +vlhx*(-(Vee-Vmm)*A(i)-Retau*G(i)+Rmutau*C(i)-Yetau*H(i)
c	5       +Ymutau*D(i)+Remu*(rhoee(i)-rhomumu(i)))


c       dbarB/dr

	   dfun(i+11*ima)=-kbA*barA(i)-kbC*barC(i)-kbG*barG(i)-kbee*barrhoee(i)
	1	-kbmumu*barrhomumu(i) + denscoef*yehere*densx*barA(i)
c	4       +vlhx*(-(Vee-Vmm)*barA(i)-Retau*barG(i)+Rmutau*barC(i)
c	5       -Yetau*barH(i)+Ymutau*barD(i)+Remu*(barrhoee(i)-barrhomumu(i)))


c       dC/dr

	   kcB=cdm21/pnu(i)*cosa*(cosg**2-sing**2)*sina*sinb - cosg
	1	*(cdm32/pnu(i)*cosb**2 + cdm21/pnu(i)*sina**2)*sing
	1	+ cdm21/pnu(i)*cosa**2*cosg*sinb**2*sing

	   kcD=-(cdm32/pnu(i)*cosb**2*cosg**2) + cdm32/pnu(i)*sinb**2
	1	+ cdm21/pnu(i)*cosa**2*(-cosb**2 + cosg**2*sinb**2)
	1	- 2.d0*cdm21/pnu(i)*cosa*cosg*sina*sinb*sing
	2	+ cdm21/pnu(i)*sina**2*sing**2

	   kcH=cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)*sinb*sing
	1	+ cdm21/pnu(i)*cosa**2*sinb*sing)

	   dfun(i+4*ima)=kcB*B(i)+kcD*D(i)+kcH*H(i)
	1	- denscoef*yehere*densx*D(i)
c	4       +vlhx*((Vee-Vtt)*D(i)+Remu*H(i)-Rmutau*B(i)+Yemu*G(i)
c	5       -Ymutau*A(i)-Yetau*(rhoee(i)-rhotautau(i)))

c       dbarC/dr

	   dfun(i+12*ima)=-kcB*barB(i)-kcD*barD(i)-kcH*barH(i)
	1	- denscoef*yehere*densx*barD(i)
c	4    +vlhx*((Vee-Vtt)*barD(i)+Remu*barH(i)-Rmutau*barB(i)
c	5    +Yemu*barG(i)-Ymutau*barA(i)-Yetau*(barrhoee(i)-barrhotautau(i)))

c       dD/dr

	   kdA=-(cdm21/pnu(i)*cosa*(cosg**2-sing**2)*sina*sinb) + cosg
	1	*(cdm32/pnu(i)*cosb**2 + cdm21/pnu(i)*sina**2)*sing
	1	- cdm21/pnu(i)*cosa**2*cosg*sinb**2*sing

	   kdC=cdm32/pnu(i)*cosb**2*cosg**2 - cdm32/pnu(i)*sinb**2
	1	+ cdm21/pnu(i)*cosa**2*(cosb**2 - cosg**2*sinb**2)
	1	+ 2.d0*cdm21/pnu(i)*cosa*cosg*sina*sinb*sing
	2	- cdm21/pnu(i)*sina**2*sing**2

	   kdG=-(cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)*sinb
	1	*sing + cdm21/pnu(i)*cosa**2*sinb*sing))

	   kdee=cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg*sinb
	1	- cdm21/pnu(i)*cosa*sina*sing)

	   kdtautau=-(cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg
	1	*sinb - cdm21/pnu(i)*cosa*sina*sing))

	   dfun(i+5*ima)= kdA*A(i)+kdC*C(i)+kdG*G(i)+kdee*rhoee(i)
	1	+kdtautau*rhotautau(i) + denscoef*yehere*densx*C(i)
c	4       +vlhx*(-(Vee-Vtt)*C(i)-Remu*G(i)+Rmutau*A(i)+Yemu*H(i)
c	5       -Ymutau*B(i)+Retau*(rhoee(i)-rhotautau(i)))

c       dbarD/dr

	   dfun(i+13*ima)=-kdA*barA(i)-kdC*barC(i)-kdG*barG(i)-kdee*barrhoee(i)
	1	-kdtautau*barrhotautau(i) + denscoef*yehere*densx*barC(i)
c	4    +vlhx*(-(Vee-Vtt)*barC(i)-Remu*barG(i)+Rmutau*barA(i)
c	5    +Yemu*barH(i)-Ymutau*barB(i)+Retau*(barrhoee(i)-barrhotautau(i)))

c       dG/dr

	   kgB=cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg*sinb
	1	- cdm21/pnu(i)*cosa*sina*sing)

	   kgD=cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)*sinb
	1	*sing + cdm21/pnu(i)*cosa**2*sinb*sing)

	   kgH=-(cdm32/pnu(i)*cosb**2*(cosg**2-sing**2)) + cdm21/pnu(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)
c	4       +vlhx*((Vmm-Vtt)*H(i)+Remu*D(i)+Retau*B(i)-Yemu*C(i)
c	5       -Yetau*A(i)-Ymutau*(rhomumu(i)-rhotautau(i)))

c       dbarG/dr

	   dfun(i+14*ima)=-kgB*barB(i)-kgD*barD(i)-kgH*barH(i)
c	4   +vlhx*((Vmm-Vtt)*barH(i)+Remu*barD(i)+Retau*barB(i)
c	5   -Yemu*barC(i)-Yetau*barA(i)-Ymutau*(barrhomumu(i)-barrhotautau(i)))

c       dH/dr

	   khA=cosb*((cdm32/pnu(i) + cdm21/pnu(i)*cosa**2)*cosg*sinb
	1	- cdm21/pnu(i)*cosa*sina*sing)

	   khC=-(cosb*(cdm21/pnu(i)*cosa*cosg*sina + cdm32/pnu(i)*sinb
	1	*sing + cdm21/pnu(i)*cosa**2*sinb*sing))

	   khG=cdm32/pnu(i)*cosb**2*(cosg**2-sing**2) + cdm21/pnu(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=cdm32/pnu(i)*cosb**2*cosg*sing - cdm21/pnu(i)*(cosa
	1	*cosg*sinb - sina*sing)*(cosg*sina + cosa*sinb*sing)

	   khtautau=-(cdm32/pnu(i)*cosb**2*cosg*sing) + cdm21/pnu(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)
c	4       +vlhx*(-(Vmm-Vtt)*G(i)-Remu*C(i)+Retau*A(i)-Yemu*D(i)
c	5       +Yetau*B(i)+Rmutau*(rhomumu(i)-rhotautau(i)))

c       dbarH/dr

	   dfun(i+15*ima)=-khA*barA(i)-khC*barC(i)-khG*barG(i)
	1	-khmumu*barrhomumu(i)-khtautau*barrhotautau(i)
c	4   +vlhx*(-(Vmm-Vtt)*barG(i)-Remu*barC(i)+Retau*barA(i)
c	5   -Yemu*barD(i)+Yetau*barB(i)+Rmutau*(barrhomumu(i)-barrhotautau(i)))

	end do

c	write(*,*)dfun
c	stop

	return 
	end

c=============================================================

	block data
	implicit double precision (a-z)
	common/ct/pi

C	Valore di pi greco.
	data pi/3.141592653589793d0/

	end