program antiwithneub
	implicit double precision (a-z)
	integer nstep,n1,n2,n3,i,j,jj,i0
	integer iask,ntheta,icount,ima,nask,kmed
c	parameter (nstep=50000)
c	parameter (nstep=10000)
*	parameter (nstep=30000)
c	parameter (nstep=5000)
	parameter (nstep=50000)
c	parameter (nvar=1044,ima=174)
	parameter (nvar=1800,ima=300)
*	parameter (nvar=600,ima=100)
c	parameter (nvar=456,ima=76)
c
	dimension pnu(ima),fnue(ima),fnumu(ima)
	dimension rdat(34),dens(34),ye(34),dummy(34)
	dimension fbarnue(ima),fbarnumu(ima)
	dimension fun(nvar),dfun(nvar)
	dimension P0(ima),Px(ima),Py(ima),Pz(ima)
	dimension Pbarx(ima),Pbary(ima),Pbarz(ima)

	common/constants1/st,ct,Rnu,neubcoef,ystep1
	common/constants2/osccoef,denscoef,off
	common/ener/Enue,Enumu,Tnue,Tnumu
	common/momenta/pnu,fnue,fnumu,fbarnue,fbarnumu
	common/density/densx,yex,yenow2
	common/neuback/intnux,intnuy,intnuz
	common/ct/pi

	double precision Lnue,Lnuebar,Lnumu

c       fix parameters here
cx	tan2t=1.d2
*	tan2t=1.0067d-2		! tan^2(theta) ???
	tan2t = 0.429		! solars
c       IH
*	tan2t=1.d0/tan2t

c       largest tan2 theta13
cc	tan2t=6.5d-2

c       dm2 in eV^2 (always positive)
c	dm2=1.d1
cx	dm2=1.d0
*	dm2=3.d-3
	dm2=8.d-5
cx	dm2=0.1701254d-1

c       atmospheric range
cx	dm2=0.8d-2
c	dm2=0.5d-2

c       calculate trigon. quantities
	st=2.d0*dsqrt(tan2t)/(1.d0+tan2t) ! sin(2*theta)? I think
	ct=(1.d0-tan2t)/(1.d0+tan2t) ! cos(2*theta)
*	write(*,*)tan2t,st,ct

*	stop

c       neutrino sphere radius in 10 km units
	Rnu=1.1d0
c       neutrino luminosity in each flavour in units 10^51 erg/sec
c       put zero if no neutrino background effect
c	Lneu=1.d4
c	Lneu=1.d2
c	Lneu=5.d2
	Lneu=1.d2		! 1.e51 erg/s
	
	Lnue = 1.d0		! Lnue/Lneu
	Lnuebar = 1.d0		! Lnuebar/Lneu
	Lnumu = 1.d0		! Lnumu/Lneu

	off=1.d0

ccc	write(*,*)'tan^2(theta)=',tan2t
ccc	write(*,*)'Dm^2=',dm2,' eV^2'
ccc	write(*,*)

c       some constants
	osccoef=5.068d4*Rnu*dm2/2.d0
	denscoef=2.548d4*Rnu*dsqrt(2.d0)
	neubcoef=3.7d7*Lneu/Rnu/(2.d0*pi*pi)

 	open(unit=40,file='snprofile_ricard.dat',status='old')
	do n2=1,34
c       radius in units Rnu, dens in 10^7 g/cm3
	   read(40,*)rdat(n2),dens(n2),ye(n2),dummy(n2)
	   rdat(n2)=rdat(n2)/1.d6/Rnu
	   dens(n2)=dens(n2)/1.d7 ! changed by a factor 1.d-20
c	   write(*,*)rdat(n2),dens(n2),ye(n2),dummy(n2)
	end do
	close(40)

c 	open(unit=30,file='evol.dat',status='old')
c 	open(unit=60,file='evolA.dat',status='old')
c 	open(unit=70,file='m1.dat',status='old')
c 	open(unit=80,file='m2.dat',status='old')
* 	open(unit=90,file='lnu.dat',status='old')
 	open(unit=90,file='lnu.dat')

c       new
	open(unit=11,file='P.dat',status='new')
	open(unit=12,file='Pbar.dat',status='new')

c       end new

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

c       evaluate FD neutrino distributions
	pmed=3.15137d0		! if eta=0 !!
	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

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

c	neutrino density matrix

	do j=1,ima

	   fnue(j)=pnu(j)**2/Tnue**4/(dexp(pnu(j)/Tnue)+1.d0)
	   fnue(j) = fnue(j)*Lnue !  ! 

	   fnumu(j)=pnu(j)**2/Tnumu**4/(dexp(pnu(j)/Tnumu)+1.d0)
	   fnumu(j) = fnumu(j)*Lnumu !"L-renormalization"

	   fbarnue(j)=pnu(j)**2/Tbarnue**4/(dexp(pnu(j)/Tbarnue)+1.d0)
	   fbarnue(j) = fbarnue(j)*Lnuebar ! "L-renormalization"

	   fbarnumu(j)=pnu(j)**2/Tbarnumu**4/(dexp(pnu(j)/Tbarnumu)+1.d0)
	   fbarnumu(j) = fbarnumu(j)*Lnumu ! "L-renormalization"

c	   write(*,*)pnu(j),fnue(j),fnumu(j)
c       P0
	   P0(j)=1.d0
c       Px
	   Px(j)=0.d0
	   fun(j)=Px(j)
c       Py
	   Py(j)=0.d0
	   fun(j+ima)=Px(j)
c       Pz
	   Pz(j)=(fnue(j)-fnumu(j))/(fnue(j)+fnumu(j))
	   fun(j+2*ima)=Pz(j)
c	   write(10,'(5e15.7)')pnu(j),fnue(j),fnumu(j),Pz(j)
c       Pxbar
	   Pbarx(j)=0.d0
	   fun(j+3*ima)=Pbarx(j)
c       Py
	   Pbary(j)=0.d0
	   fun(j+4*ima)=Pbary(j)
c       Pz
	   Pbarz(j)=fbarnue(j)-fbarnumu(j)
	   Pbarz(j)=Pbarz(j)/(fbarnue(j)+fbarnumu(j))
	   fun(j+5*ima)=Pbarz(j)
cxc	   write(50,*)pnu(j),dabs(Pz(j)),dabs(Pbarz(j))
cx	   write(*,*)pnu(j),Pz(j),Pbarz(j)
	end do

c	stop

	rnue=0.d0
	nnue=0.d0
	ranue=0.d0
	nanue=0.d0
	rnumu=0.d0
	nnumu=0.d0

	do j=1,ima
	   rnue=rnue+pnu(j)*0.5d0*(1.d0+Pz(j))*(fnue(j)+fnumu(j))
	   nnue=nnue+0.5d0*(1.d0+Pz(j))*(fnue(j)+fnumu(j))
	   ranue=ranue
     #    +pnu(j)*0.5d0*(1.d0+Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   nanue=nanue
     #    +0.5d0*(1.d0+Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   rnumu=rnumu+pnu(j)*0.5d0*(1.d0-Pz(j))*(fnue(j)+fnumu(j))
	   nnumu=nnumu+0.5d0*(1.d0-Pz(j))*(fnue(j)+fnumu(j))

	end do

	avnue=rnue/nnue
	avanue=ranue/nanue
	avnumu=rnumu/nnumu
	Lnuei=nnue*avnue*ystep1
	Lnuebari=nanue*avanue*ystep1
	Lnumui=nnumu*avnumu*ystep1

cx	write(*,*)avnue,avnumu,avanue
cx	write(*,*)Lnuei,Lnuebari,Lnumui
cx	write(*,*)avPx_nue,avPy_nue,avPz_nue
cx	write(*,*)avPx_barnue,avPy_barnue,avPz_barnue

c	stop
c       range in r (evolution)

c	rmin=1.05d0
	rmin=1.12d0! unitats en 1.e7 cm
c	rmin=1.08d0
cxx	rmin=2.d0

c	rmax=1.d2
c	rmax=1.d2/Rnu
c	rmax=1.d1/Rnu
c	rmax=3.18d0
c	rmin=1.5d0
c	rmax=3.d0/Rnu

c       for 1, 10
c	rmax=5.d0
c	rmax=1.d1

c       for 0.05, 0.005
c	rmax=55.d0

	rmax=10000.d0
*	rmax=500.d0


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

c       log
	rstep = (log10(rmax)-log10(rmin))/float(nstep-1) ! dlogr


	r1=rmin

	icount=1
	do i=1,nstep
c       linear
*	   r2=r1+rstep
*	   h0=rstep/1.d2

c       log
	   r2= log10(r1)+rstep
	   r2 = 10.d0**r2
	   h0=r2*rstep/1.d2



*** check
	   write(*,*)i,r1,r2,rstep,r1*Rnu,r2*Rnu
***

c       interpolate values of density and Ye for currrent r
	   do jj=1,34
	      if(rdat(jj).ge.r2)then
		 rratio=dlog(r2)-dlog(rdat(jj-1))
		 rratio=rratio/(dlog(rdat(jj))-dlog(rdat(jj-1)))

		 ldensx=dlog(dens(jj))-dlog(dens(jj-1))
		 ldensx=ldensx*rratio+dlog(dens(jj-1))
		 densx=dexp(ldensx)

		 lyex=dlog(ye(jj))-dlog(ye(jj-1))
		 lyex=lyex*rratio+dlog(ye(jj-1))
		 yex=dexp(lyex)

c		 densx=(dens(jj)-dens(jj-1))*(r2-rdat(jj-1))
c		 densx=densx/(rdat(jj)-rdat(jj-1))+dens(jj-1)
c		 yex=(ye(jj)-ye(jj-1))*(r2-rdat(jj-1))
c		 yex=yex/(rdat(jj)-rdat(jj-1))+ye(jj-1)
c		 write(10,'(5e15.7)')r2,densx,yex
		 goto 33
	      end if
	   end do
	   if(r2.gt.rdat(34))then
	      densx=dens(34)*(rdat(34)/r2)**3
	      yex=ye(34)
c	      write(*,*)'radius outside range!!'
c	      stop
	   end if

 33	   continue


*** check
*	   densx = 1.d1/1d7	!
*	   yex = 0.5d0
	   
	   write(91,*)r2,densx,yex
	   write(*,*)r2,densx,yex
***	   


c       calculate integrals over neutrino dist
	   intnuznue=0.d0
	   intnuzbarnue=0.d0
	   intnuznumu=0.d0

	   do jj=1,ima
c       just for check (113)
	      intnuznue=intnuznue+fnue(jj)
	      intnuznumu=intnuznumu+fnumu(jj)
	      intnuzbarnue=intnuzbarnue+fbarnue(jj)
	   end do


c       calculate geometrical factor
	   Frgeo=0.5d0*(1.d0-dsqrt(1.d0-1.d0/(r1*r1)))
	   Frgeo=Frgeo/r1/r1

c----------------------------------------------
c       Delta V diagonal check area

	   goto 113

	   intnuznue=intnuznue*Frgeo*neubcoef*ystep1
	   intnuzbarnue=intnuzbarnue*Frgeo*neubcoef*ystep1
	   intnuznumu=intnuznumu*Frgeo*neubcoef*ystep1

c       pmuestra shows example in MeV
	   pmuestra=1.d1
	   pmuestra=1.d0
	   dm2eff=denscoef*yex*densx*3.946d-5/Rnu*pmuestra
	   neueff=intnuz*3.946d-5/Rnu*pmuestra
	   neueff1=intnuznue*3.946d-5/Rnu*pmuestra
	   neueff2=intnuzbarnue*3.946d-5/Rnu*pmuestra
	   neueff3=intnuznumu*3.946d-5/Rnu*pmuestra
c	   write(50,*)r1*Rnu*1.d1,intnuz/r1/r1
c	   write(50,*)r1*Rnu*1.d1,dm2eff,neueff,
c     #          neueff1,neueff2,neueff3
	   write(50,*)r1*Rnu*1.d1,dm2eff,
     #          neueff1,neueff2,neueff3
c----------------------------------------------

 113	    continue
	    call evol1(r1,r2,fun,h0)

	   do j=1,ima

	      Px(j)=fun(j)
	      Py(j)=fun(j+ima)
	      Pz(j)=fun(j+2*ima)
	      Pbarx(j)=fun(j+3*ima)
	      Pbary(j)=fun(j+4*ima)
	      Pbarz(j)=fun(j+5*ima)

	   end do

	   if(icount.eq.3)then
c	   if(icount.eq.6)then
c	   if(icount.eq.15)then
c	   if(icount.eq.10)then
c	   if(icount.eq.20)then

c       find average energies
	      rnue=0.d0
	      nnue=0.d0
	      ranue=0.d0
	      nanue=0.d0
	      rnumu=0.d0
	      nnumu=0.d0
	      ranumu=0.d0
	      nanumu=0.d0
	      rnue2=0.d0
	      ranue2=0.d0
	      rnumu2=0.d0
	      
	      do j=1,ima
	   rnue2=rnue2+pnu(j)*pnu(j)*0.5d0*
     #          (1.d0+Pz(j))*(fnue(j)+fnumu(j))
	   rnue=rnue+pnu(j)*0.5d0*(1.d0+Pz(j))*(fnue(j)+fnumu(j))
	   nnue=nnue+0.5d0*(1.d0+Pz(j))*(fnue(j)+fnumu(j))
	   ranue2=ranue2+pnu(j)*pnu(j)*0.5d0*
     #          (1.d0+Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   ranue=ranue
     #    +pnu(j)*0.5d0*(1.d0+Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   nanue=nanue
     #    +0.5d0*(1.d0+Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   rnumu2=rnumu2+pnu(j)*pnu(j)*0.5d0*
     #          (1.d0-Pz(j))*(fnue(j)+fnumu(j))
	   rnumu=rnumu+pnu(j)*0.5d0*(1.d0-Pz(j))*(fnue(j)+fnumu(j))
	   nnumu=nnumu+0.5d0*(1.d0-Pz(j))*(fnue(j)+fnumu(j))

	   ranumu=ranumu
     #    +pnu(j)*0.5d0*(1.d0-Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	   nanumu=nanumu
     #    +0.5d0*(1.d0-Pbarz(j))*(fbarnue(j)+fbarnumu(j))
	      end do

	      avnue2=rnue2/nnue
	      avanue2=ranue2/nanue
	      avnumu2=rnumu2/nnumu

	      avnue=rnue/nnue
	      avanue=ranue/nanue
	      avnumu=rnumu/nnumu
	      avanumu=ranumu/nanumu

	      lnue=nnue*avnue*ystep1/Lnuei
	      lnuebar=nanue*avanue*ystep1/Lnuebari
	      lnumu=nnumu*avnumu*ystep1/Lnumui
	      lnumubar=nanumu*avanumu*ystep1/Lnumui

c	      yenow=1.d0/(1.d0+Ebarnue/avnue)
	      yenow=1.d0/(1.d0+lnuebar*avanue/lnue/avnue)
	      yenow2=1.d0/(1.d0+nanue*avanue2/nnue/avnue2)

ccc	      write(*,'(5e15.7)')r2,yenow,yenow2
c	      write(50,'(5e15.7)')r2,avnue,avnumu,yenow,avanue,avanumu
	      write(50,'(5e15.7)')r2,avnue,avnumu,avanue,avanumu
c	      write(40,'(5e15.7)')r2,avnue2,avnumu2,yenow2,avanue2
	      write(40,'(5e15.7)')r2,yenow,yenow2
c	      write(*,'(5e15.7)')r2,lnue,lnuebar,lnumu,lnumubar
	      write(90,'(5e15.7)')r2,lnue,lnuebar,lnumu,lnumubar

c       for 76 modes (13,30,39 MeV) set #2
c	      write(*,'(5e15.7)')r2,Pz(10),Pz(23),Pz(30)
cx	      write(30,'(5e15.7)')r2,Pz(10),Pz(23),Pz(30)
cx	      write(60,'(5e15.7)')r2,Pbarz(10),Pbarz(23),Pbarz(30)
cx	      write(70,'(5e15.7)')r2,Px(10),Py(10),Pz(10)
cx	      write(80,'(5e15.7)')r2,Px(30),Py(30),Pz(30)
c	      write(50,'(5e15.7)')r2,intnux/Frgeo/neubcoef,
c     #        intnuy/Frgeo/neubcoef,intnuz/Frgeo/neubcoef

c       for 100 modes (6.5,9.5,11.5,15,5 MeV)
c	      write(*,'(5e15.7)')r2,Pz(7),Pz(10),Pz(16)
	      write(30,'(5e15.7)')r2,
     #        Pz(7)*(fnue(7)+fnumu(7))/(fnue(7)-fnumu(7)),
     #        Pz(10)*(fnue(10)+fnumu(10))/(fnue(10)-fnumu(10)),
     #        Pz(12)*(fnue(12)+fnumu(12))/(fnue(12)-fnumu(12)),
     #        Pz(16)*(fnue(16)+fnumu(16))/(fnue(16)-fnumu(16))
	      write(60,'(5e15.7)')r2,
     #  Pbarz(7)*(fbarnue(7)+fbarnumu(7))/(fbarnue(7)-fbarnumu(7)),
     #  Pbarz(10)*(fbarnue(10)+fbarnumu(10))/(fbarnue(10)-fbarnumu(10)),
     #  Pbarz(12)*(fbarnue(12)+fbarnumu(12))/(fbarnue(12)-fbarnumu(12)),
     #  Pbarz(16)*(fbarnue(16)+fbarnumu(16))/(fbarnue(16)-fbarnumu(16))

*	      open(unit=20,file='fnu.dat',status='old')
	      open(unit=20,file='fnu.dat')
	      do n1=1,ima
		 fae=0.5d0*(1.d0+Pz(n1))
		 fae=fae*(fnue(n1)+fnumu(n1))
		 famu=0.5d0*(1.d0-Pz(n1))
		 famu=famu*(fnue(n1)+fnumu(n1))
		 faantie=0.5d0*(1.d0+Pbarz(n1))
		 faantie=faantie*(fbarnue(n1)+fbarnumu(n1))
		 faantimu=0.5d0*(1.d0-Pbarz(n1))
		 faantimu=faantimu*(fbarnue(n1)+fbarnumu(n1))
	    write(20,'(5e15.7)')pnu(n1),fae,famu,faantie,faantimu
	 end do
	 close(20)


c new average of Pi
	   avPx_nue=0.d0
	   avPy_nue=0.d0
	   avPz_nue=0.d0
	   avPx_barnue=0.d0
	   avPy_barnue=0.d0
	   avPz_barnue=0.d0


	   do j=1,ima
	   check=(fnue(j)+fnumu(j))/(fnue(j)-fnumu(j))/float(ima)
	   avPx_nue=avPx_nue+Px(j)*check
	   avPy_nue=avPy_nue+Py(j)*check
	   avPz_nue=avPz_nue+Pz(j)*check

	   check2=(fbarnue(j)+fbarnumu(j))/(fbarnue(j)-fbarnumu(j))
	   check2=check2/float(ima)
	   avPx_barnue=avPx_barnue+Pbarx(j)*check2
	   avPy_barnue=avPy_barnue+Pbary(j)*check2
	   avPz_barnue=avPz_barnue+Pbarz(j)*check2
	   
	   end do
	   write(11,'(5e15.7)')r2,avPx_nue,avPy_nue,avPz_nue
	   write(12,'(5e15.7)')r2,avPx_barnue,avPy_barnue,avPz_barnue
c	   write(*,*)r2,avPx_nue,avPy_nue,avPz_nue
	    
c end new

	 goto 115

	   open(unit=10,file='ratios.dat',status='old')
	      do n1=1,ima
		 c1=denscoef*yex*densx !*Py(n1)
		 c2=osccoef/pnu(n1)*ct !*Py(n1)
cx		 c3=intnuy*Pz(n1)-intnuz*Py(n1)
		 c3=intnuz !*Py(n1)
	    write(10,'(5e15.7)')pnu(n1),c1,c2,c3
	 end do
	 close(10)

 115	 continue

	      icount=0
	   end if

 111	   r1=r2	      

	   icount=icount+1

	end do
cx	close(30)
cx	close(60)
cx	close(70)
cx	close(80)
	close(90)

c        new
	close(11)
	close(12)
c       end new
	end

c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 	subroutine evol1(x1,x2,fun,h0)

	implicit double precision (a-z)
	integer nstep,nvar,ifail,i,j,ima
c	parameter (nvar=1044,ima=174)
	parameter (nvar=1800,ima=300)
*	parameter (nvar=600,ima=100)
c	parameter (nvar=456,ima=76)

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	eps=1.d-8
c	eps=1.d-2


*** check
*        write(*,*)'entrem en d02cj ...'
***

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


	subroutine dfundx(rx,fun,dfun)
	implicit double precision (a-z)
	integer nvar,i,j,jj,k,itest,ima
c	parameter (nvar=1044,ima=174)
	parameter (nvar=1800,ima=300)
*	parameter (nvar=600,ima=100)
c	parameter (nvar=456,ima=76)
	dimension pnu(ima),fnue(ima),fnumu(ima)
	dimension fbarnue(ima),fbarnumu(ima)
	dimension fun(nvar),dfun(nvar)
	dimension P0(ima),Px(ima),Py(ima),Pz(ima)
	dimension Pbarx(ima),Pbary(ima),Pbarz(ima)

	common/constants1/st,ct,Rnu,neubcoef,ystep1
	common/constants2/osccoef,denscoef,off
	common/ct/pi
	common/momenta/pnu,fnue,fnumu,fbarnue,fbarnumu
	common/density/densx,yex,yenow2
	common/neuback/intnux,intnuy,intnuz

	do jj=1,ima
	   Px(jj)=fun(jj)
	   Py(jj)=fun(jj+ima)
	   Pz(jj)=fun(jj+2*ima)
	   Pbarx(jj)=fun(jj+3*ima)
	   Pbary(jj)=fun(jj+4*ima)
	   Pbarz(jj)=fun(jj+5*ima)
	end do

 100	continue

c       calculate integrals over neutrino dist
	   intnux=0.d0
	   intnuy=0.d0
	   intnuz=0.d0

	   do jj=1,ima
	      intnux=intnux+(fnue(jj)+fnumu(jj))*Px(jj)
     #            -(fbarnue(jj)+fbarnumu(jj))*Pbarx(jj)
	      intnuy=intnuy+(fnue(jj)+fnumu(jj))*Py(jj)
     #            -(fbarnue(jj)+fbarnumu(jj))*Pbary(jj)
	      intnuz=intnuz+(fnue(jj)+fnumu(jj))*Pz(jj)
     #            -(fbarnue(jj)+fbarnumu(jj))*Pbarz(jj)

	   end do


c       calculate geometrical factor
	   Frgeo=0.5d0*(1.d0-dsqrt(1.d0-1.d0/(rx*rx)))
	   Frgeo=Frgeo/rx/rx

	   intnux=intnux*Frgeo*neubcoef*ystep1
	   intnuy=intnuy*Frgeo*neubcoef*ystep1
	   intnuz=intnuz*Frgeo*neubcoef*ystep1

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

c	rkm=rx*Rnu*1.d1
c	if(rkm.le.1.2485d1)then
	   yehere=yex
c	else
c	   yehere=yenow2
c	end if

	do k=1,ima

c       dPx/dt
	   dfun(k)= -(denscoef*yehere*densx-osccoef/pnu(k)*ct)*Py(k)
     #          +off*(intnuy*Pz(k)-intnuz*Py(k))

c       dPy/dt
	   dfun(k+ima)= (denscoef*yehere*densx-osccoef/pnu(k)*ct)*Px(k)
     #          -osccoef/pnu(k)*st*Pz(k)
     #          +off*(intnuz*Px(k)-intnux*Pz(k))

c       dPz/dt
	   dfun(k+2*ima)= osccoef/pnu(k)*st*Py(k)
     #          +off*(intnux*Py(k)-intnuy*Px(k))

c       dPbarx/dt
	dfun(k+3*ima)= (-denscoef*yehere*densx-osccoef/pnu(k)*ct)*Pbary(k)
     #          +off*(intnuy*Pbarz(k)-intnuz*Pbary(k))

c       dPbary/dt
	dfun(k+4*ima)= (denscoef*yehere*densx+osccoef/pnu(k)*ct)*Pbarx(k)
     #          +osccoef/pnu(k)*st*Pbarz(k)
     #          +off*(intnuz*Pbarx(k)-intnux*Pbarz(k))

c       dPbarz/dt
	   dfun(k+5*ima)= -osccoef/pnu(k)*st*Pbary(k)
     #          +off*(intnux*Pbary(k)-intnuy*Pbarx(k))

	end do

	return 
        end

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

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

	end