c THE RATES COMPUTED BY THIS PROGRAM ALL NEED TO BE MULTIPLED BY 1/2. C THIS IS A MISTAKE THAT WAS CAUGHT LATE IN THE PROCESS SO I WON'T C CORRECT IT HERE c This program reads the results of modes_efield.f and computes a c rate. Implicit None Integer bj,i,j,ik,kn,nk,n,junk,nbands,ib,nok,nbad,ibik,it,nu Real*16 ki,delk,k,B2,pi,b,u1,m Real*16 sum,delt,rate,Ikz,ktest,p,rateex,ratedif Real*16 int,intrho,aout,aouttest,Itest,inttest,Itexact,diftest Real*16 Ar,Ai,Br,Bi,testw,B2test,difB2 Real*16 w0,u,uin,om0,Ip1,Ip2,Irho,Irate,Ip2tot,Ip1tot,Irhotot,Iratetot Real*16 intrate,itime,eps,h1,hmin,a,intp1,intp2 Real*16 sinOm,cosOm,rho,p1,p2,difaout,itimein Real*16 ditimedt,umatch,utime,delu,kf,uprint,test1,test2 Real*16 volscaled,Iratescaled Parameter (bj=10,kn=1000000) Dimension k(bj,kn),intrate(bj,kn),intrho(bj,kn),intp1(bj,kn),intp2(bj,kn) Dimension ki(bj),delk(bj),nk(bj),itime(bj),ditimedt(bj),Ar(bj,kn),Ai(bj,kn),Br(bj,kn),Bi(bj,kn) External Timeint,Bsstep Common/param/b,u1 OPEN(UNIT = 11, FILE = 'rate_ds.dat') OPEN(UNIT = 12, FILE = 'B2.1') OPEN(UNIT = 13, FILE = 'rate_ds.out') Open(unit = 14, file = 'rate_ds.tst') Open(unit = 15, file='pressure.integrand') Open(unit = 16, file='pressure.out') Open(unit = 17, file='pressure.integrand1') c This next file outputs the integral for the rate and the volume for c the rate and the rate and the energy density and the ave pressure Open(unit = 18, file = 'rate_ds.plt') Write(18,*)'0 b,u1,Iratescaled,volscaled,rate,rho,p1 ' READ(11,*) JUNK Read(11,*) nbands,m,u1,b,delu,nu,uprint Read(11,*) JUNK Read(11,*) eps,h1,hmin Read(11,*) JUNK aout = qexp(u1) c Let's test this and get the matching time from the file B2.1 Read(12,*) aouttest,umatch If(qabs((aout-aouttest)/aout).gt.1.q-6) then Print*,' Problem with aout. In terms of u1 we have aout = ',aout,' In terms of B2.1 we have aout = ',aouttest Stop Endif c I am computing the integral for the time part of the 4-volume by solving a c first order differential equation so I need the starting value at time c u = -u1 uin = - u1 c This is the scale factor at time u = -u1 a = qcosh(-2*b*u1)**(-1.q0/(2.q0*b)) c itime is the time part of the 4-volume that I want to divide by to c get the rate c This is wrong! c itime(1) = a**3 itime(1) = 0.q0 Call Odeint(itime,1,uin,u1,eps,h1,hmin,nok,nbad,Timeint,Bsstep) c This was the old version where the pressure was computed at only one c time c Now read in the time to do the computation of the pressure at. If it is 0 then we'll match at time t0 c Print*,'What time do you want to compute the pressure at? Saying 0 means the matching time.' c Read(*,*) utime c If(qabs(utime).lt.1.q-10) then c utime = umatch c Endif Do ib=1,nbands c Set up rate_ds.dat so that for the first band ki = 0. READ(11,*) ki(ib),delk(ib), nk(ib) kf = ki(ib)+delk(ib)*(nk(ib)-1) Print*,'ib = ',ib,'ki(ib) = ',ki(ib) Print*,' kf(ib) = ',kf Enddo pi = 2.q0*qasin(1.q0) Do ib=1,nbands ibik = 20+ib If(ib.eq.1) then Open(21,file='AB.1') Else if(ib.eq.2) then Open(22,file='AB.2') Else if(ib.eq.3) then Open(23,file='AB.3') Else if(ib.eq.4) then Open(24,file='AB.4') Else if(ib.eq.5) then Open(25,file='AB.5') Endif Do ik = 1,nk(ib) If(ik.eq.1) then k(ib,ik) = ki(ib) Else k(ib,ik) = k(ib,ik-1)+delk(ib) Endif c Print out the values at the boundaries If(ik.eq.1.or.ik.eq.nk(ib)) then Print*,'ib = ',ib,' ik = ',ik Print*,' k = ',k(ib,ik) Endif c We don't want to read anything for the very first value of k since it c is zero If(ib.eq.1.and.ik.eq.1) Then B2 = 0.q0 intrate(ib,ik) = 0.q0 intrho(ib,ik) = 0.q0 intp1(ib,ik) = 0.q0 intp2(ib,ik) = 0.q0 Go to 999 Endif 777 Read(ibik,*,end=998) ktest,Ar(ib,ik),Ai(ib,ik),Br(ib,ik),Bi(ib,ik) If(qabs((k(ib,ik)-ktest)/ktest).gt.1.q-10) then c Print*,' Wrong value of k from AB.ib, ib = ',ib,' ik = ',ik c Print*,' k = ',k(ib,ik),' ktest = ',ktest Go to 777 c Stop Endif B2 = Br(ib,ik)**2+Bi(ib,ik)**2 om0 = qsqrt(k(ib,ik)**2/aout**2 + m**2) c We set the phase to zero at the matching time umatch, so that at any time u so c we have cos(2*om0*(u-umatch)) and sin(2*om0*(u-umatch)) c In this version we'll compute the pressure at the matching time for c the first pass and then go to later times after that c cosOm = qcos(2.q0*om0*(utime-umatch)) c sinOm = qsin(2.q0*om0*(utime-umatch)) cosOm = 1.q0 sinOm = 0.q0 intrate(ib,ik) = k(ib,ik)**2*qlog(1.q0+B2) intrho(ib,ik) = k(ib,ik)**2*om0*B2 c Write(16,19) k(ib,ik),om0,B2,intrho(ib,ik) c The first one is the |B|^2 contribution to the pressure intp1(ib,ik) = 2.q0*k(ib,ik)**4*B2/(om0*aout**2) c The second is the oscillating part of the integrand for the pressure intp2(ib,ik) = -(m**2*k(ib,ik)**2/om0)*(2.q0*Ai(ib,ik)*Bi(ib,ik)*cosOm + - 2.q0*Ar(ib,ik)*Br(ib,ik)*cosOm - 2.q0*Ar(ib,ik)*Bi(ib,ik)*sinOm + - 2.q0*Ai(ib,ik)*Br(ib,ik)*sinOm) 999 Write(15,19) k(ib,ik),intrate(ib,ik),intrho(ib,ik),intp1(ib,ik),intp2(ib,ik) Enddo Enddo c Now compute the integrals. It is useful to have two sets of loops so c that in the first the integrand can be easily printed out and tests done c and in the second the integrals can be computed c c I'll compute the integrals for each band and also keep a running total Iratetot = 0.q0 Irhotot = 0.q0 Ip1tot = 0.q0 Ip2tot = 0.q0 Do ib=1,nbands Irate = 0.q0 Irho = 0.q0 Ip1 = 0.q0 Ip2 = 0.q0 Itest = 0.q0 Do ik = 1,nk(ib)-2,2 c This is to test the integration routine Itest = Itest+ (k(ib,ik)**2+4.q0*k(ib,ik+1)**2+k(ib,ik+2)**2)*delk(ib)/3.q0 Irate = Irate + (intrate(ib,ik)+4.q0*intrate(ib,ik+1)+intrate(ib,ik+2))*delk(ib)/3.q0 Irho = Irho + (intrho(ib,ik)+4.q0*intrho(ib,ik+1)+intrho(ib,ik+2))*delk(ib)/3.q0 Ip1 = Ip1 + (intp1(ib,ik)+4.q0*intp1(ib,ik+1)+intp1(ib,ik+2))*delk(ib)/3.q0 Ip2 = Ip2 + (intp2(ib,ik)+4.q0*intp2(ib,ik+1)+intp2(ib,ik+2))*delk(ib)/3.q0 Enddo Itexact = (1.q0/3.q0)*(k(ib,nk(ib))**3-k(ib,1)**3) diftest = (Itexact-Itest)/Itexact Write(14,*) ' Test integral for band number ',ib,' diftest = ',diftest c Print out the contribution to the rate for the band rate = Irate/(2.q0*pi**2*itime(1)) rho = Irho/(2.q0*pi**2*aout**3) Write(14,*) ' rate = ', rate, ' rho = ', rho Write(14,*) ' Ip1 = ', Ip1, ' Ip2 = ',Ip2 Iratetot = Iratetot+Irate Irhotot = Irhotot + Irho Ip1tot = Ip1tot + Ip1 Ip2tot = Ip2tot + Ip2 Enddo c I want a scaled volume and a scaled integral so I can get the u1 c dependence out volscaled = 2.q0*pi**2*itime(1)/aout**3 Iratescaled = Iratetot/aout**3 rate = Iratetot/(2.q0*pi**2*itime(1)) rho = Irhotot/(2.q0*pi**2*aout**3) p1 = Ip1tot/(12.q0*pi**2*aout**3) p2 = Ip2tot/(12.q0*pi**2*aout**3) p = p1+p2 Write(18,19) b,u1,Iratescaled,volscaled,rate,rho,p1 Write(13,*) 'Test integral results' Write(13,19) Itexact,Itest,diftest Write(13,*) 'Results for b = ',b,' and u1 = ',u1 Write(13,*) 'itime rate rho ptotal p from B2 p from oscillaions' Write(13,19)itime(1),rate,rho,p,p1,p2 c Now compute the pressure at the other times u = umatch-delu c u and t are the same essentially in the static region Do it = 1,nu+1 c I want to print out the actual time - not the time since umatch u = u+delu Do ib =1,nbands Do ik = 1,nk(ib) om0 = qsqrt(k(ib,ik)**2/aout**2 + m**2) cosOm = qcos(2.q0*om0*(u-umatch)) sinOm = qsin(2.q0*om0*(u-umatch)) intp2(ib,ik) = -(m**2*k(ib,ik)**2/om0)*(2.q0*Ai(ib,ik)*Bi(ib,ik)*cosOm + - 2.q0*Ar(ib,ik)*Br(ib,ik)*cosOm - 2.q0*Ar(ib,ik)*Bi(ib,ik)*sinOm + - 2.q0*Ai(ib,ik)*Br(ib,ik)*sinOm) If(qabs(u-uprint).lt.1.q-5) then test1 = Ai(ib,ik)*Bi(ib,ik)+Ar(ib,ik)*Br(ib,ik) test2 = -Ar(ib,ik)*Bi(ib,ik)+Ai(ib,ik)*Br(ib,ik) Write(17,19) k(ib,ik),intp2(ib,ik),om0,cosOm,sinOm,test1,test2 Endif Enddo Enddo Ip2tot = 0.q0 Do ib=1,nbands Ip2 = 0.q0 Do ik = 1,nk(ib)-2,2 Ip2 = Ip2 + (intp2(ib,ik)+4.q0*intp2(ib,ik+1)+intp2(ib,ik+2))*delk(ib)/3.q0 Enddo c Print out the contribution to the integral for the band Write(14,*) ' Ip2 = ',Ip2 Ip2tot = Ip2tot + Ip2 Enddo p2 = Ip2tot/(12.q0*pi**2*aout**3) p = p1+p2 Write(16,19) u,p1,p2,p Enddo 19 Format(1x,9(d14.7,1x)) Stop 998 Print*,'It got to the end of a file.' Print*,' Wrong value of k from AB.ib, ib = ',ib,' ik = ',ik Print*,' k = ',k(ib,ik),' ktest = ',ktest Stop End Subroutine Timeint(u,itime,ditimedt,nvar) Implicit None Integer nvar Real*16 u,u1,b,itime,ditimedt,a Dimension itime(nvar),ditimedt(nvar) Common/param/b,u1 c Print*,'Got to subroutine' a = qexp((-qlog(qcosh(b*(u-u1)))+qlog(qcosh(b*(u+u1))))/(2.q0*b)) ditimedt(1) = a**3 c Print*,'u = ',u,' u1 = ',u1,' b = ',b c Print*,'a = ',a,'ditimedt = ',ditimedt(1) Return End