c---------------------------------------------------------------- c simple_dmet.f 020707 HS c = cleaned version of ssbox_dmet_p9_new_ff_eh_2007p4.f c removed all viscosity etc. tabulations c but friction, adhesion retained c simple_dmet_rota.f 100707 HS c spin related lines uncommented c---------------------------------------------------------------- c Integrates particle orbits with RK4, including c central force, selfgravity and impact forces for overlapping pairs c To speed up calculations, selfgravity is not updated at every RK4 step c Instead, selfgravity acceleration A0 and its derivative dA0/dt are c returned by simple_gravi.f after each large timestep dt, while c the used RK4 step is dtstep=dt/alfstep (typically alfstep=40) c The list of potentially colliding pairs is also updated after each dt c In case of selfgravity this is done in simple_gravi.f, c in the nongravitating case in simple_nongravi.f c An impact can be in progress when a new list of potentially colliding c particles is calculated -> taken into account via a temporary list c of indices c-------------------------------------------------------------------------- c Equations of motion: c-------------------------------------------------------------------------- c 1) central force c ax=at(i,1)= w2*rpt(i,2)-factemp*rt(i,1) c ay=at(i,2)=-w2*rpt(i,1) c az=at(i,3)=-wz2*rt(i,3) c-------------------------------------------------------------------------- c 2) selfgravity c [ax,ay,az]= A= A0+ dA0/dt (t-t0) c where t-t0 is the time since the beginning of the large timestep dt c gmet=2 -> no self-gravity (A0=dA0/dt=0) c gmet=5 -> self-gravity constant in RK4 step (dA0/dt=0) c gmet=6 -> use dA0/dt=aagrav in RK4 step to improve gravity c-------------------------------------------------------------------------- c 3) impact forces c Normal component c (= in the direction joining the particle centers) c a_normal=-(alfome2*x + xp/alftau) c calculation of alfome2,alftau : c ometau=sqrt(pii**2/(2.*log(alfeps))**2+0.25) c alfome2=alfome**2 c alftau=ometau/alfome c Tangential component c (= in the direction of relative velocity to tangent plane of impact) c a_tangential=ft c lfric= coefficient of friction (negative) c fricmodel = model for friction (added 23.01.03) c fricmodel=0 -> as before c mu_fric=lfric, with beta limitation: c mu_lfric set to zero if gt/gt0 as above, but zigzag prevented by using ft*dt < gt c =2 -> mu_fric chosen to yield beta, ft=mu_fric*fn c =3 -> kmu_fric chosen to yield beta, ft=kmu_fric*gt*fn c Adhesion c F_adhesion= 3*pi*gamma_adh*r_eff c gamma_adh > 0 -> add force between all overlapping pairs c-------------------------------------------------------------------------- c Also possible to ignore certain type of impacts: c nob_imp = 1 -> ignore impacts over all boundaries c 2 -> ignore impacts over radial boundaries c 3 -> ignore impacts over tangential boundaries c nob_imp = 4 -> ignore all impacts c Normally nob_imp=0 --> calculate all impacts, taking correctly into c account the periodic bounday conditions c-------------------------------------------------------------------------- c total of 9*LKM integration variables c R,RP,ROTA c rotation related quantities commented out in the current version c-------------------------------------------------------------------------- c*************************************************************** SUBROUTINE DIMPACTS_RKHF4 C*************************************************************** INCLUDE 'ssbox.inc' c variables passed to RK4: dimension yy(nstar*9),dydx(nstar*9) c dimension yy(nstar*6),dydx(nstar*6) double precision lfriccol c store for collisions in progress common /rkhf/alfold(nix),taucol(nix), + epscol(nix),vcicol(nix),timcol(nix),vticol(nix), + lfriccol(nix) c divide each gravity step into ALFSTEP smaller integration steps c calculate impact force between overlapping particles call cpu(time1_dimp) cpuapu_ori=cpu_int+cpu_pota+cpu_sav c---------------------------------------------------------------- c FIRST ALFSTEP c---------------------------------------------------------------- c in the first step: c the current arrays holding the impacting particles need to be updated c since the INDIX-array (list of potential colliders = all nearby pairs) c has changed --> c 1) search II,JJ of present overlaps c 2) check the values IIOLD, JJOLD corresponding to II and JJ c from the table strored at prev. call to DIMPACTS c 3) transform the values: epssav(IIC) -> EPSCOL(ix(ii)) c taucol = stored alftau <-- tausav =spring parameter c alfold = prev alf <-- alfsav =overlap depth (pos for overlap) c iic <-- iisav =overlapping particle1 c iic <-- jjsav =overlapping particle2 c epscol = stored eps <-- epssav c vcicol = stored vcinit <-- vcisav v_nor initial normal velocity c vticol = stored vtinit <-- vtisav v_tan initial tangential velocity c timcol = stored tstart <-- timsav c NSAV = number of stored overlapping pairs c NSAV = 0 -> update unnecessary c except that alfodl(i) is set to a large neg. number call cpu(time1) do 5 i=1,ix alfold(i)=-1000. lfriccol(i)=lfric 5 continue if(NSAV.ne.0) then numimp=0 do 10 inew=1,ix if(eroix(inew).le.0) then ii=indix(inew,1) jj=indix(inew,2) ssum=s(ii)+s(jj) ero=eroix(inew) alf=ssum-ero c--------UPDATE do 20 ic=1,nsav iiold=iisav(ic) jjold=jjsav(ic) if(iiold.eq.ii.and.jjold.eq.jj) then numimp=numimp+1 alfold(inew)=alfsav(ic) taucol(inew)=tausav(ic) epscol(inew)=epssav(ic) vcicol(inew)=vcisav(ic) timcol(inew)=timsav(ic) vticol(inew)=vtisav(ic) lfriccol(inew)=lfricsav(ic) goto 21 endif 20 continue 21 continue endif 10 continue endif call cpu(time2) cpu_sav=cpu_sav+(time2-time1) c---------------------------------------------- c LOOP OVER ALFSTEP c---------------------------------------------- timesg=time dtstep=dt/alfstep alfmax=0. alfmean=0. novermean=0 lkm1=lkm lkm2=2*lkm1 lkm3=3*lkm1 lkm4=4*lkm1 lkm5=5*lkm1 lkm6=6*lkm1 lkm7=7*lkm1 lkm8=8*lkm1 c---------------------------------------------- c initialize r,rp,rota -> yy c---------------------------------------------- call cpu(time1) do 120 i=1,lkm1 yy(i)=r(i,1) yy(i+lkm1)=r(i,2) yy(i+lkm2)=r(i,3) yy(i+lkm3)=rp(i,1) yy(i+lkm4)=rp(i,2) yy(i+lkm5)=rp(i,3) yy(i+lkm6)=rota(i,1) yy(i+lkm7)=rota(i,2) yy(i+lkm8)=rota(i,3) 120 continue call cpu(time2) cpu_int=cpu_int+(time2-time1) c---------------------------------------------- c INTEGRATION OVER ALFSTEP small steps c---------------------------------------------- c idetect=1 -> detect starting/ending of impact c idetect=0 -> do not detect -"- c viscosity tabulation: weight=vflag c different in each RK-substep c 1./6.+2./6.+2./6.+1./6. do 130 istep=1,alfstep t=timesg+1.d0*(istep-1)*dtstep idetect=1 vflag=1./6. call derivs_hf4(t,yy,dydx,lkm1,idetect,vflag) call rk4_hf4(yy,dydx,lkm1,t,dtstep,yy) 130 continue c just to update forces in .oidl idetect=0 vflag=0. t=t+dtstep call derivs_hf4(t,yy,dydx,lkm1,idetect,vflag) c---------------------------------------------- c store back yy -> r,rp, rota c---------------------------------------------- call cpu(time1) do 140 i=1,lkm1 r(i,1)=yy(i) r(i,2)=yy(i+lkm1) r(i,3)=yy(i+lkm2) rp(i,1)=yy(i+lkm3) rp(i,2)=yy(i+lkm4) rp(i,3)=yy(i+lkm5) rota(i,1)=yy(i+lkm6) rota(i,2)=yy(i+lkm7) rota(i,3)=yy(i+lkm8) 140 continue call cpu(time2) cpu_int=cpu_int+(time2-time1) c---------------------------------------------------------------- c LAST ALFSTEP c---------------------------------------------------------------- c store arrays for currently overlapping particles c NSAV = number of overlaps c taucol = stored alftau --> tausav c alfold = prev alf --> alfsav c ii --> iisav c jj --> jjsav c epscol = eps(v) --> epssav c vci = initial vn c vti = initial vt c timsav = start of impact c lfricsav = call cpu(time1) nsav=0 do 200 i=1,ix if(alfold(i).gt.0) then nsav=nsav+1 iisav(nsav)=indix(i,1) jjsav(nsav)=indix(i,2) alfsav(nsav)=alfold(i) tausav(nsav)=taucol(i) epssav(nsav)=epscol(i) vcisav(nsav)=vcicol(i) vtisav(nsav)=vticol(i) timsav(nsav)=timcol(i) lfricsav(nsav)=lfriccol(i) endif 200 continue call cpu(time2) cpu_sav=cpu_sav+(time2-time1) c---------------------------------------------------------------- c current values during step: c alfmax = max overlap (relative to ssum) c alfmean = mean overlap (relative to ssum) c novermean = mean number of overlaps if(novermean.ne.0) alfmean=alfmean/novermean novermean=novermean/alfstep call cpu(time2_dimp) cpu_potb=cpu_potb+time2_dimp-time1_dimp- + (cpu_int+cpu_pota+cpu_sav-cpuapu_ori) c---------------------------------------------------------------- c update system time time=time+dt return end c****************************************************************** SUBROUTINE RK4_Hf4(Y,DYDX,N,X,H,YOUT) c****************************************************************** C N VARIABLES Y AND THEIR DERIVATES DYDX KNOWN AT X C TIMESTEP H, INCREMENTED VARIABLES YOUT c from numerical recipes INCLUDE 'ssbox.inc' PARAMETER (NMAX=9*nstar) c PARAMETER (NMAX=6*nstar) DIMENSION Y(*),DYDX(*),YOUT(*) DIMENSION YT(NMAX),DYT(NMAX),DYM(NMAX),DYS(NMAX) HH=.5d0*H H6=H/6.d0 XH=X+HH N6=9*N c N6=6*N c do not detect impacts here idetect=0 c viscosity tabulation: weight=vflag DO 11 I=1,N6 YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE vflag=2./6. CALL DERIVS_hf4(XH,YT,DYT,N,idetect,vflag) DO 12 I=1,N6 YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE vflag=2./6. CALL DERIVS_hf4(XH,YT,DYM,N,idetect,vflag) DO 13 I=1,N6 YT(I)=Y(I)+H*DYM(I) 13 CONTINUE vflag=1./6. CALL DERIVS_hf4(X+H,YT,DYS,N,idetect,vflag) DO 14 I=1,N6 YOUT(I)=Y(I)+H6*(DYDX(I)+2.d0*DYT(I)+2.d0*DYM(I)+DYS(I)) 14 CONTINUE RETURN END C*************************************************************** SUBROUTINE DERIVS_Hf4(T,YY,DYDX,N,idetect,vflag) C*************************************************************** c routine which constructs forces INCLUDE 'ssbox.inc' common /rkhf/alfold(nix),taucol(nix), + epscol(nix),vcicol(nix),timcol(nix),vticol(nix), + lfriccol(nix) dimension rt(nstar,3),rpt(nstar,3),at(nstar,3) dimension rotat(nstar,3) dimension yy(*),dydx(*) double precision cc(3),apuvec1(3),apuvec2(3) double precision vcoll(3),cvc(3),ccvc(3) double precision lfriccol double precision drota(nstar,3) c r,v,f for collision partners i,j double precision meff,fparti(3),fpartj(3) c idetect=1 -> detect starting/ending of impact c idetect=0 -> do not -"- alfome2=alfome**2 lkm1=n lkm2=2*lkm1 lkm3=3*lkm1 lkm4=4*lkm1 lkm5=5*lkm1 lkm6=6*lkm1 lkm7=7*lkm1 lkm8=8*lkm1 bvak=vvak*t bvak=-bvak bvak=bvak-int(bvak/2./yl)*2.*yl bvak=-bvak if(bvak.gt.yl) bvak=bvak-2.*yl if(bvak.lt.-yl) bvak=bvak+2.*yl dt_small=dt/alfstep c----------------------------------------------------------------- c insert yy to rt,rpt,rotat c----------------------------------------------------------------- call cpu(time1) do 99 i=1,lkm1 rt(i,1)=yy(i) rt(i,2)=yy(i+lkm1) rt(i,3)=yy(i+lkm2) rpt(i,1)=yy(i+lkm3) rpt(i,2)=yy(i+lkm4) rpt(i,3)=yy(i+lkm5) rotat(i,1)=yy(i+lkm6) rotat(i,2)=yy(i+lkm7) rotat(i,3)=yy(i+lkm8) 99 continue c----------------------------------------------------------------- c 0. c CONSTRUCT self-gravitational force c----------------------------------------------------------------- c gravi (gmet=5 or gmet=6) has returned: c AGRAV = gravitational force (far+near) c in the beginning of alfstep (t=timesg) c AAGRAV = derivative of nearby force c gmet=5 -> zero c gmet=6 -> contains all nearby pairs c here: c AGRAVP = gravity affecting each particle c gmet=5 -> set AGRAVP=AGRAV c gmet=6 -> set AGRAVP=AGRAV+AAGRAV*(t-timesg) c time since self-gravity was calculated dtgrav=t-timesg c for gmet=2 (nongravi) agrav and aagrav normally zero c Included only to allow for axisymmetric gravity c or some other later additions if(gmet.eq.2) then do 11 i=1,lkm1 agravp(i,1)=agrav(i,1)+aagrav(i,1)*dtgrav agravp(i,2)=agrav(i,2)+aagrav(i,2)*dtgrav agravp(i,3)=agrav(i,3)+aagrav(i,3)*dtgrav 11 continue endif if(gmet.eq.5) then do 1 i=1,lkm1 agravp(i,1)=agrav(i,1) agravp(i,2)=agrav(i,2) agravp(i,3)=agrav(i,3) 1 continue endif if(gmet.eq.6) then do 2 i=1,lkm1 agravp(i,1)=agrav(i,1)+aagrav(i,1)*dtgrav agravp(i,2)=agrav(i,2)+aagrav(i,2)*dtgrav agravp(i,3)=agrav(i,3)+aagrav(i,3)*dtgrav 2 continue endif c----------------------------------------------------------------- c 1. c CONSTRUCT DERIVATIVES DUE TO CENTRAL FORCE, MUTUAL GRAVITY c AND MOONLET c----------------------------------------------------------------- c----------------------------------- c 1.1 central force and gravity: c----------------------------------- do 100 i=1,lkm1 c uncomment the next lines if rotation is put back c (change is just due to rotation of the coordinate system) drota(i,1)=w0*rotat(i,2) drota(i,2)=-w0*rotat(i,1) drota(i,3)=0. c central force+gravity at(i,1)= w2*rpt(i,2)-factemp*rt(i,1)+agravp(i,1) at(i,2)=-w2*rpt(i,1)+agravp(i,2) at(i,3)=-wz2*rt(i,3)+agravp(i,3) 100 continue call cpu(time2) cpu_int=cpu_int+(time2-time1) c ignore impacts ? if(nob_imp.ge.4) goto 777 c----------------------------------------------------------------- c 2. DERIVATIVES DUE TO IMPACT FORCES c LOOP OVER LIST OF POTENTIALLY COLLIDING PAIRS c (these were constructed in simple_gravi.f or simple_nongravi.f) c----------------------------------------------------------------- call cpu(time1) nover=0 do 110 i=1,ix c---------------------------------------- c 2.0 choose correct images c---------------------------------------- ii=indix(i,1) jj=indix(i,2) mm=mmix(i) nn=nnix(i) ssum=s(jj)+s(ii) dx=rt(ii,1)-rt(jj,1)+nn*avak dy=rt(ii,2)-rt(jj,2)+nn*bvak+mm*cvak dz=rt(ii,3)-rt(jj,3) nncor=nn mmcor=mm if(dy.gt.yl) then dy=dy-2.*yl mmcor=mmcor-1 endif if(dy.lt.-yl) then dy=dy+2.*yl mmcor=mmcor+1 endif if(dx.gt.xl) then dx=dx-2.*xl nncor=nncor-1 endif if(dx.lt.-xl) then dx=dx+2.*xl nncor=nncor+1 endif ero2=dx*dx+dy*dy+dz*dz ero=sqrt(ero2) alf=ssum-ero c----------------------------------------------------------- c in case of no overlap skip impact calculation if(ero2.gt.ssum*ssum) goto 111 c----------------------------------------------------------------------- c 2.2 START TREATING OVERLAPPING PAIRS c----------------------------------------------------------------------- nover=nover+1 c 2.2.1 relative velocity dvx=rpt(ii,1)-rpt(jj,1) dvy=rpt(ii,2)-rpt(jj,2)+nncor*vvak dvz=rpt(ii,3)-rpt(jj,3) cv=dx*dvx+dy*dvy+dz*dvz alfp=-cv/ero apuvec1(1)=rotat(ii,1)+rotat(jj,1) apuvec1(2)=rotat(ii,2)+rotat(jj,2) apuvec1(3)=rotat(ii,3)+rotat(jj,3)-ero*w0 cc(1)=-dx/ero cc(2)=-dy/ero cc(3)=-dz/ero call dcross(apuvec1,cc,apuvec2) vcoll(1)=-dvx-apuvec2(1) vcoll(2)=-dvy-apuvec2(2) vcoll(3)=-dvz-apuvec2(3) call dcross(vcoll,cc,apuvec1) call dcross(cc,apuvec1,cvc) vtinit=sqrt(cvc(1)**2+cvc(2)**2+cvc(3)**2) if(vtinit.ne.0) then cvc(1)=cvc(1)/vtinit cvc(2)=cvc(2)/vtinit cvc(3)=cvc(3)/vtinit endif call dcross(cc,cvc,ccvc) c------------------------------------------------------- c 2.2.2 c NEW IMPACT STARTS ? c detected only if idetect=1 c------------------------------------------------------- c------------------------------------- c 2.2.2.1 c call DILCOL_RKHf4 if alfold < 0 c calculate eps c store values in the beginning of impact if(idetect.eq.1.and.alfold(i).lt.0.d0) then call DILCOL_RKHf4(I,dx,dy,dz,dvx,dvy,dvz, + alftau,eps,vcinit) c 180603 c ignore impacts over boundaries? if(nob_imp.ne.0) then if(nob_imp.eq.1) then if(nncor.ne.0) alftau=0. if(mmcor.ne.0) alftau=0. endif if(nob_imp.eq.2) then if(nncor.ne.0) alftau=0. endif if(nob_imp.eq.3) then if(mmcor.ne.0) alftau=0. endif endif c note: vcinit=abs(cv) taucol(i)=alftau epscol(i)=eps vcicol(i)=vcinit timcol(i)=t vticol(i)=vtinit lfriccol(i)=lfric c end: if(idetect.eq.1.and.alfold(i).lt.0.d0) endif c------------------------------------------------------- c 2.2.3 c CALCULATE FORCES BETWEEN IMPACTING PARTICLES c taucol(i)=0 -> skip this c------------------------------------------------------- alftau=taucol(i) if(alftau.eq.0) goto 111 c FIRST SOME STATISTICS: c alfmax = maximum overlap (relative to smaller particle) c ero1max= maximum pre-dimpact separation leading to impact smin=s(jj) if(s(ii).le.s(jj)) smin=s(ii) if(alf/smin.ge.alfmax) alfmax=alf/smin alfmean=alfmean+alf/smin if(eroix(i).ge.ero1max) ero1max=eroix(i) c------------------------------------------------------- c 2.2.3.1 c NORMAL IMPACT FORCE c------------------------------------------------------- c linear force model c remember: particle 1<->ii, 2<->jj c so: (dx,dy,dz)/ero= - C apu=alfome2*alf+alfp/alftau apu1=myy(jj)/(myy(ii)+myy(jj)) apu2=-myy(ii)/(myy(ii)+myy(jj)) meff=myy(ii)*myy(jj)/(myy(ii)+myy(jj)) do 800 k=1,3 fparti(k)=-apu*meff*cc(k) fpartj(k)= apu*meff*cc(k) 800 continue do 801 k=1,3 at(ii,k)=at(ii,k)+fparti(k)/myy(ii) at(jj,k)=at(jj,k)+fpartj(k)/myy(jj) 801 continue c------------------------------------------------ c 2.2.3.2 c ADHESIVE FORCE c------------------------------------------------ gamma_adh_use=0. if(gamma_adh.gt.0) gamma_adh_use=gamma_adh if(gamma_adh_use.gt.0) then r_eff=s(ii)*s(jj)/(s(ii)+s(jj)) apu_adh=-3.*pii*gamma_adh_use*r_eff/meff do 810 k=1,3 fparti(k)=-apu_adh*meff*cc(k) fpartj(k)= apu_adh*meff*cc(k) 810 continue do 811 k=1,3 at(ii,k)=at(ii,k)+fparti(k)/myy(ii) at(jj,k)=at(jj,k)+fpartj(k)/myy(jj) 811 continue endif c-------------------------------------------------------- c 2.2.3.3 c TANGENTIAL FORCE c------------------------------------------------------- c FRICTION: lfric corresponds to mu c limit friction so that eps_t > (1-beta) c vticol = tangential velocity in the beginning of impact c vtinit = current tangential velocity if(fricmodel.eq.0) then if(vtinit.lt.vticol(i)*(1.-beta)) lfriccol(i)=0. if(vticol(i).eq.0) lfriccol(i)=0. ft=0. if(lfriccol(i).ne.0.and.apu.gt.0) ft=-lfric*apu endif if(fricmodel.eq.1) then if(vtinit.lt.vticol(i)*(1.-beta)) lfriccol(i)=0. if(vticol(i).eq.0) lfriccol(i)=0. ft=0. if(lfriccol(i).ne.0.and.apu.gt.0) ft=-lfric*apu c prevent zigzag if(ft*dt_small.gt.vtinit) then ft=vtinit/dt_small*.5 endif endif if(fricmodel.eq.2) then ft=0. if(vcicol(i).ne.0) then lfriccol(i)=2.d0/7.d0*beta* + vticol(i)/vcicol(i)/(1.d0+epscol(i)) ft=lfriccol(i)*apu endif endif if(fricmodel.eq.3) then ft=0. if(vcicol(i).ne.0) then lfriccol(i)=-2./7.*dlog(1.d0-beta) + /vcicol(i)/(1.d0+epscol(i)) ft=lfriccol(i)*apu*vtinit endif endif c------------------------------------------------------ c INCLUDE TANGENTIAL FORCES c unless lfriccol = 0 c or apu<0 if(lfriccol(i).eq.0) goto 111 if(apu.le.0) goto 111 do 820 k=1,3 fparti(k)= meff*ft*cvc(k) fpartj(k)=-meff*ft*cvc(k) 820 continue do 821 k=1,3 at(ii,k)=at(ii,k)+fparti(k)/myy(ii) at(jj,k)=at(jj,k)+fpartj(k)/myy(jj) 821 continue c changes in spins, assuming alfa=2/5 c (J=alfa*m*R^2 moment of inertia) apu1=myy(jj)/(myy(jj)+myy(ii)) apu2=-myy(ii)/(myy(jj)+myy(ii)) c uncomment if needed drota(ii,1)=drota(ii,1)+2.5*apu1*ft*ccvc(1) drota(ii,2)=drota(ii,2)+2.5*apu1*ft*ccvc(2) drota(ii,3)=drota(ii,3)+2.5*apu1*ft*ccvc(3) drota(jj,1)=drota(jj,1)-2.5*apu2*ft*ccvc(1) drota(jj,2)=drota(jj,2)-2.5*apu2*ft*ccvc(2) drota(jj,3)=drota(jj,3)-2.5*apu2*ft*ccvc(3) 9797 continue c--------------------------------------------------------------------------- c END OVERLAP c OR NO OVERLAP c--------------------------------------------------------------------------- 111 continue c---------------------------------------------------------------------- c idetect=1 ? c set alfold(i) to current value if(idetect.eq.1) then alfold(i)=alf lfriccol(i)=lfric endif 110 continue call cpu(time2) cpu_pota=cpu_pota+(time2-time1) novermean=novermean+nover c jump here if no impacts 777 continue c-------------------------------------------------------------- C construct derivatives c-------------------------------------------------------------- call cpu(time1) do 120 i=1,lkm1 dydx(i)=rpt(i,1) dydx(i+lkm1)=rpt(i,2) dydx(i+lkm2)=rpt(i,3) dydx(i+lkm3)=at(i,1) dydx(i+lkm4)=at(i,2) dydx(i+lkm5)=at(i,3) dydx(i+lkm6)=drota(i,1) dydx(i+lkm7)=drota(i,2) dydx(i+lkm8)=drota(i,3) 120 continue call cpu(time2) cpu_int=cpu_int+(time2-time1) return end C******************************************************* SUBROUTINE DILCOL_RKHF4(icol,dx,dy,dz,dvx,dvy,dvz,alftau,eps, + vcinit) c******************************************************* C CALLED IN THE BEGINNING OF IMPACT C CALCULATE EPS(C.V) , TAU(ALFOME,EPS) C STORE IMPACT-STATISTICS INCLUDE 'ssbox.inc' double precision lamda_use ssmin=s(1) ssmax=s(lkm) c-------------------------------------------------- c RELATIVE POSITION AND VELOCITY c-------------------------------------------------- ii=indix(icol,1) jj=indix(icol,2) mm=mmix(icol) nn=nnix(icol) ssum=s(jj)+s(ii) ero2=dx*dx+dy*dy+dz*dz ero=sqrt(ero2) cv=(dx*dvx+dy*dvy+dz*dvz)/ero cv0=cv cv=abs(cv) size1=s(ii) size2=s(jj) ssum=size1+size2 is1=is(ii) is2=is(jj) c-------------------------------------------------- c IMPACT MODEL c-------------------------------------------------- C IMODEL=1 -> CONSTANT ALFA = LAMDA C IMODEL=2 -> ALFA(V)= (C.V/V0)**(-.234) BRIDGES C IMODEL=3 -> ALFA=LAMDA/(1+(C.V)/V0) C IMODEL=4 -> ALFA=LAMDA*exp(v0*v(cm/sec)) C IMODEL=5 -> ALFA=LAMDA*v**(v0*v(cm/sec)) c IMODEL=6 -> DILLEY I c IMODEL=7 -> DILLEY II c IMODEL=8 -> DILLEY III c IMODEL=9 -> DILLEY ksi0=lamda p=v0 c addition 17.3.99: c IMODEL=10 -> tabulated eps(v) imodel_use=imodel lamda_use=lamda v0_use=v0 if(imodel.eq.11) then lamda_use=alftab(is1,is2) imodel_use=1 endif if(imodel.eq.22) then v0_use=v0tab(is1,is2) imodel_use=2 endif call ssbox_imodel(imodel_use,lamda_use,v0_use, + cv,size1,size2,alfa, + ssmin,ssmax,is1,is2) if(coll.lt.10) then write(*,*) 'iscoll,is1,is2,alfa ',coll,is1,is2,alfa write(*,*) ' orb,ii, jj, cv ',orb,ii,jj,cv0 endif C ****** SLOW IMPACT -> ALFA=1 C STORE NUMBER OF SLOW IMPACTS cvcrit=vcrita*ssum*w0 IF(cv.LT.cvcrit) then alfa=1.d0 nhidas=nhidas+1 endif IF(cv.LT. 10.*cvcrit) NHIDAS10=NHIDAS10+1 IF(cv.LT. 5.0*cvcrit) NHIDAS5 =NHIDAS5+1 IF(cv.LT. 2.0*cvcrit) NHIDAS2 =NHIDAS2+1 IF(cv.LT. 0.5*cvcrit) NHIDAS05=NHIDAS05+1 c--------------------------------------------------- c EPS, ALFTAU, VCINIT c alfa=1 -> dissipation term set 1/alftau is zero eps=alfa vcinit=cv if(alfa.ge.1.) alfa=0.999d0 ometau=sqrt(pii**2/(2.*log(eps))**2+0.25) alftau=ometau/alfome IF(cv.gt.cvcrit) then COLL=COLL+1 COLLIT=COLLIT + 1 endif RETURN END C**************************************************** SUBROUTINE DCROSS(X,Y,Z) C**************************************************** c vector product Z = X * Y double precision X(3),Y(3),Z(3) Z(1)=X(2)*Y(3)-X(3)*Y(2) Z(2)=X(3)*Y(1)-X(1)*Y(3) Z(3)=X(1)*Y(2)-X(2)*Y(1) RETURN END C**************************************************** SUBROUTINE ssbox_imodel(imodel,lamda,v0,cv,size1,size2,alfa, + smin,smax,is1,is2) C**************************************************** c-------------------------------------------------------- c elasticity model specified by: imodel,lamda,v0 c calculates alfa=eps_n normal coefficient of restitution c as a function of cv= v \dot c c includes also size1 and size2 for models 6-9 c checked 26.03.02 HS c add back imodel=11 and imodel=22 c--------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) double precision lamda c-------------------------------------------------- c IMPACT MODEL c-------------------------------------------------- C IMODEL=1 -> CONSTANT ALFA = LAMDA C IMODEL=2 -> ALFA(V)= (C.V/V0)**(-.234) BRIDGES C IMODEL=3 -> ALFA=LAMDA/(1+(C.V)/V0) C IMODEL=4 -> ALFA=LAMDA*exp(v0*v(cm/sec)) C IMODEL=5 -> ALFA=LAMDA*v**(v0*v(cm/sec)) c IMODEL=6 -> DILLEY I c IMODEL=7 -> DILLEY II c IMODEL=8 -> DILLEY III c IMODEL=9 -> DILLEY ksi0=lamda p=v0 c new impact models c IMODEL=-11 BRIDGES eps1=0.32*v^(-0.24) FROSTED T=210 K r=2.5 c IMODEL=-12 HATZES eps2=0.48*v^(-0.20) FROSTED T=123 K r=2.5 c IMODEL=-13 HATZES compacted frost r=2.5 cm c eps3=0.90*exp(-0.395*v) +0.002*v^(-0.9) c IMODEL=-14 HATZES compacted frost r=20 cm c eps4=0.90*exp(-0.22*v) +0.01*v^(-0.6) c OBS in GG-book Model IV had -0.21 instead of -0.22 c OBS in GG-book Model IV had -0.385 instead of 0.395 C----------------------------------------------------- c 1 -> alfa=const=lamda IF(IMODEL.EQ.1) ALFA=LAMDA c 2 -> bridges84 IF(IMODEL.EQ.2) THEN VIMP=ABS(cv/V0) ALFA=1.D0 IF(VIMP.NE.0) ALFA=VIMP**(-.234D0) IF(ALFA.GT.1) ALFA=1.D0 IF(ALFA.LT..25) ALFA=.25D0 ENDIF c 3 -> SALO87 etc. IF(IMODEL.EQ.3) ALFA=LAMDA/(1.D0+cv/v0) c 4 -> hatzes88 eksponential lamda*exp(v0*v(cm/s)) IF(IMODEL.EQ.4) THEN VIMP=cv*100. ALFA=lamda*exp(v0*vimp) IF(ALFA.GT.1) ALFA=1.D0 IF(ALFA.LT..25) ALFA=.25D0 ENDIF c 5 -> hatzes88 power-law (frosty) lamda**(v0*v(cm/s)) IF(IMODEL.EQ.5) THEN VIMP=cv*100. ALFA=lamda*vimp**v0 IF(ALFA.GT.1) ALFA=1.D0 IF(ALFA.LT..25) ALFA=.25D0 ENDIF c 6,7,8,9 -> DILLEY92 IF(IMODEL.ge.6.and.imodel.le.9) THEN xsi0=lamda vexp=v0 if(imodel.eq.6) then xsi0=0.16 vexp=0.65 endif if(imodel.eq.7) then xsi0=0.25 vexp=0.43 endif if(imodel.eq.8) then xsi0=0.34 vexp=0.36 endif VIMP=cv*100 sizemin=size1 sizemax=size2 if(size1.gt.size2) then sizemin=size2 sizemax=size1 endif deltax=sizemin/sizemax xsi=xsi0*vimp**vexp*(1.+deltax)**0.2* + (sizemin/0.025)**(-0.2) alfa=exp(-3.1415*xsi/sqrt(1.-xsi*xsi)) if(alfa.ge.1) alfa=1. if(alfa.le.0) alfa=0.0001 endif c---------------------------------------------------------------------- c bridges84 power T=210 IF(IMODEL.EQ.-11) THEN VIMP=cv*100 ALFA=0.32*vimp**(-0.24) IF(ALFA.GT.1) ALFA=1.D0 END IF c hatzes88 power T=123 IF(IMODEL.EQ.-12) THEN VIMP=cv*100 ALFA=0.48*vimp**(-0.20) IF(ALFA.GT.1) ALFA=1.D0 END IF c hatzes88 for r=2.5cm IF(IMODEL.EQ.-13) THEN VIMP=cv*100 ALFA=0.9*exp(-0.395*vimp)+0.002*vimp**(-0.9) IF(ALFA.GT.1) ALFA=1.D0 END IF c hatzes88 for r=20cm IF(IMODEL.EQ.-14) THEN VIMP=cv*100 ALFA=0.9*exp(-0.22*vimp)+0.01*vimp**(-0.6) IF(ALFA.GT.1) ALFA=1.D0 END IF c---------------------------------------------------------------- c new models 160404 if(imodel.eq.20) then epsmin=0.1 epsmax=0.5 alfa1=epsmin+(size1-smin)/(smax-smin)*(epsmax-epsmin) alfa2=epsmin+(size2-smin)/(smax-smin)*(epsmax-epsmin) alfa=0.5*(alfa1+alfa2) endif if(imodel.eq.21) then epsmin=0.1 epsmax=0.6 alfa1=epsmin+(size1-smin)/(smax-smin)*(epsmax-epsmin) alfa2=epsmin+(size2-smin)/(smax-smin)*(epsmax-epsmin) alfa=0.5*(alfa1+alfa2) endif if(imodel.eq.30) then epsmin=0.5 epsmax=0.1 alfa1=epsmin+(size1-smin)/(smax-smin)*(epsmax-epsmin) alfa2=epsmin+(size2-smin)/(smax-smin)*(epsmax-epsmin) alfa=0.5*(alfa1+alfa2) endif if(imodel.eq.31) then epsmin=0.6 epsmax=0.1 alfa1=epsmin+(size1-smin)/(smax-smin)*(epsmax-epsmin) alfa2=epsmin+(size2-smin)/(smax-smin)*(epsmax-epsmin) alfa=0.5*(alfa1+alfa2) endif if(imodel.eq.40) then epsmin=lamda epsmax=v0 alfa1=epsmin+(size1-smin)/(smax-smin)*(epsmax-epsmin) alfa2=epsmin+(size2-smin)/(smax-smin)*(epsmax-epsmin) alfa=0.5*(alfa1+alfa2) endif c---------------------------------------------------------------------- RETURN END