c********************************************************************** c* This subroutine uses the 4th order Runge-Kutta method to * c* integrate the ODE. The equation of motion is defined in * c* Levine's paper except that we use (alpha*v - beta*v^3) instead * c* of Levine's (alpha*v_unit - beta*v). * c* The old configuration is input from "xold" and then the new * c* configuration after the integration will be output in "xnew". * c* (The content in "xold" will not be altered.) * c********************************************************************** subroutine method(t,dt,fm,lm,ifirst,ilast, & xold, xnew, & xm,vm,km) c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include "./linearad.i" c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer fm,lm,ifirst,ilast double precision t,dt double precision xold(fm:lm,1:4), xnew(fm:lm,1:4) double precision xm(fm:lm,1:2), vm(fm:lm,1:2), km(fm:lm,1:4) integer ie, je, irk double precision vlen, xlen, force(1:2) double precision cla,clr,erk c ****************************************************************** cla=ca/la clr=cr/lr c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Initiate the Runge-Kutta method c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ do ie=ifirst,ilast do je=1,4 xnew(ie,je)=xold(ie,je) km(ie,je)=0.d0 enddo enddo c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c Start the Runge-Kutta method (irk=1,2,3,4 - the 4 steps of RK4) c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ do irk=1,4 if (irk.eq.1) then erk=0.d0 else if(irk.eq.4) then erk=1.d0 else erk=0.5d0 endif do ie=ifirst, ilast xm(ie,1)=xold(ie,1)+erk*km(ie,1) xm(ie,2)=xold(ie,2)+erk*km(ie,2) vm(ie,1)=xold(ie,3)+erk*km(ie,3) vm(ie,2)=xold(ie,4)+erk*km(ie,4) enddo do ie=ifirst,ilast vlen=(vm(ie,1)**2+vm(ie,2)**2)**0.5d0 force(1)=0.d0 force(2)=0.d0 c Interaction forces do je=ifirst,ilast xlen=( (xm(ie,1)-xm(je,1))**2 & +(xm(ie,2)-xm(je,2))**2 )**0.5d0 if (xlen.ne.0.d0) then force(1) = force(1) & -cla*(xm(ie,1)-xm(je,1))/xlen*exp(-xlen/la) & +clr*(xm(ie,1)-xm(je,1))/xlen*exp(-xlen/lr) force(2) = force(2) & -cla*(xm(ie,2)-xm(je,2))/xlen*exp(-xlen/la) & +clr*(xm(ie,2)-xm(je,2))/xlen*exp(-xlen/lr) endif enddo c Calculate the increments km(ie,1)=dt*vm(ie,1) km(ie,2)=dt*vm(ie,2) km(ie,3)=dt*( alpha*vm(ie,1) & -beta*vm(ie,1)*vlen*vlen & +force(1) ) km(ie,4)=dt*( alpha*vm(ie,2) & -beta*vm(ie,2)*vlen*vlen & +force(2) ) enddo c Update the increments to "xnew" if ((irk.eq.1).or.(irk.eq.4)) then erk=1.d0/6.d0 else erk=1.d0/3.d0 endif do ie=ifirst,ilast do je=1,4 xnew(ie,je)=xnew(ie,je)+erk*km(ie,je) enddo enddo enddo return end c************************************************************************* c* Subroutine "method2" is an old subroutine that uses the averaged * c* self-propelling force defined in Levine's paper. * c************************************************************************* subroutine method2(t,dt,fm,lm,ifirst,ilast, & xold, xnew, & xm,vm,km) c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ include "./linearad.i" c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ integer fm,lm,ifirst,ilast double precision t,dt double precision xold(fm:lm,1:4), xnew(fm:lm,1:4) double precision xm(fm:lm,1:2), vm(fm:lm,1:2), km(fm:lm,1:4) integer ie, je, irk double precision vlen, xlen, force(1:2),univel(1:2),univlen double precision cla,clr,erk c ****************************************************************** c multiply fluxes by dt (ie, compute time integral over cell side) c ****************************************************************** cla=ca/la clr=cr/lr do ie=ifirst,ilast do je=1,4 xnew(ie,je)=xold(ie,je) km(ie,je)=0.d0 enddo enddo do irk=1,4 if (irk.eq.1) then erk=0.d0 else if(irk.eq.4) then erk=1.d0 else erk=0.5d0 endif do ie=ifirst, ilast c write(*,*) ie xm(ie,1)=xold(ie,1)+erk*km(ie,1) xm(ie,2)=xold(ie,2)+erk*km(ie,2) vm(ie,1)=xold(ie,3)+erk*km(ie,3) vm(ie,2)=xold(ie,4)+erk*km(ie,4) enddo do ie=ifirst,ilast vlen=(vm(ie,1)**2+vm(ie,2)**2)**0.5d0 force(1)=0.d0 force(2)=0.d0 univel(1)=0.d0 univel(2)=0.d0 c if (t.lt.0) then c if (ie.lt.ilast/2) then c do je=ifirst,ilast/2 c xlen=( (xm(ie,1)-xm(je,1))**2 c & +(xm(ie,2)-xm(je,2))**2 )**0.5d0 c force(1) = force(1) c & -cla*(xm(ie,1)-xm(je,1))*exp(-xlen/la) c & +clr*(xm(ie,1)-xm(je,1))*exp(-xlen/lr) c force(2) = force(2) c & -cla*(xm(ie,2)-xm(je,2))*exp(-xlen/la) c & +clr*(xm(ie,2)-xm(je,2))*exp(-xlen/lr) c enddo c else c do je=ilast/2,ilast c xlen=( (xm(ie,1)-xm(je,1))**2 c & +(xm(ie,2)-xm(je,2))**2 )**0.5d0 c force(1) = force(1) c & -cla*(xm(ie,1)-xm(je,1))*exp(-xlen/la) c & +clr*(xm(ie,1)-xm(je,1))*exp(-xlen/lr) c force(2) = force(2) c & -cla*(xm(ie,2)-xm(je,2))*exp(-xlen/la) c & +clr*(xm(ie,2)-xm(je,2))*exp(-xlen/lr) c enddo c endif c else do je=ifirst,ilast if(je.ne.ie) then xlen=( (xm(ie,1)-xm(je,1))**2 & +(xm(ie,2)-xm(je,2))**2 )**0.5d0 force(1) = force(1) & -cla*(xm(ie,1)-xm(je,1))/xlen*exp(-xlen/la) & +clr*(xm(ie,1)-xm(je,1))/xlen*exp(-xlen/lr) force(2) = force(2) & -cla*(xm(ie,2)-xm(je,2))/xlen*exp(-xlen/la) & +clr*(xm(ie,2)-xm(je,2))/xlen*exp(-xlen/lr) c endif c if(je.ne.ie) then univlen=(vm(je,1)**2+vm(je,2)**2)**0.5d0 if (univlen.ne.0.d0) then univel(1)= univel(1) & +vm(je,1)/univlen*exp(-xlen/lc) univel(2)= univel(2) & +vm(je,2)/univlen*exp(-xlen/lc) endif endif enddo c endif univlen=(univel(1)**2+univel(2)**2)**0.5d0 km(ie,1)=dt*vm(ie,1) km(ie,2)=dt*vm(ie,2) km(ie,3)=dt*( alpha*vm(ie,1) & -beta*vm(ie,1)*vm(ie,1)*vm(ie,1) & +force(1) ) km(ie,4)=dt*( alpha*vm(ie,2) & -beta*vm(ie,2)*vm(ie,2)*vm(ie,2) & +force(2) ) enddo if ((irk.eq.1).or.(irk.eq.4)) then erk=1.d0/6.d0 else erk=1.d0/3.d0 endif do ie=ifirst,ilast do je=1,4 xnew(ie,je)=xnew(ie,je)+erk*km(ie,je) enddo enddo enddo return end