Solve the Briliantov, Poschel, et al visco-elastic impact model in dimensionless form. 
MYL 11/5/06

References:
C.Thornton& Z.Ning:"a theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres," Powder Tech.99,154-162 (1998).
TN_theory_dim-.nb
Ram\'{\i}rez, et al [R. Ram\'{\i}rez, T. P\"oschel, N. V. Brilliantov and T Schwager: 
 ``Coefficient of restitution of colliding viscoelastic spheres,'' PRE {\bf 60}, 4465-4472 (1999)],
 N.V. Brilliantov, F. Spahn, J.-M. Hertzsch and T. P\"oschel:``Model for collisions in granular gases,'' PRE 53, 5382-5392 (1996).

Instructions:
	implement statements in red
	adjust performance using statements in blue

Variables:
h	dim(-) interpenetration relative to maximum Hertz value hmax
Dvdw	dimensionless parameter 
	= Sqrt[8π/3] dd^(1/2)*gam^(1/2)*rr^(1/4)/hmax^(3/4)
dd	= (3/2) (1-nu^2)/E with nu = Poisson ratio; E = Young's modulus (Pa)
gam	surface adhesion energy (J/m^2)
rr	sphere radius (m)
a	dim(-) contact patch radius relative to Sqrt[rr*hmax]

aa4[h_,d_]:=(1/12) ((Sqrt[3]) √(8 h+((32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+6 √(((4 h)/3)-((8 (2^(2/3)) (h^2))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))-((((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))/(6 (2^(2/3))))+((Sqrt[3]) (d^2))/(√(8 h+((32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))));force[a_,d_]:=2^(3/2)*3*Pi/5*(Gamma[2/5]/Gamma[9/10])^2*(4*a^3/3 - d*a^(3/2));forceh2[h_,d_]:=force[aa4[h,d],d];hhmin[d_]:=-3/8*d^(4/3);fvisc[h_,d_]:=(1/(5 ((Gamma[(9/10)])^2)))(6 (Sqrt[2]) (-(1/(16 (Sqrt[3])))(d (((Sqrt[3]) (8-(32 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3)))+(64 (2^(2/3)) h)/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(1536 (h^2)+((13824 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(3 (((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3)))))/(2 √(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))+(3 ((4/3)+(8 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(9 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3)))-(768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(18 (2^(2/3)) (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3)))-(16 (2^(2/3)) h)/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))-((Sqrt[3]) (d^2) (8-(32 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3)))+(64 (2^(2/3)) h)/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(1536 (h^2)+((13824 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(3 (((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3)))))/(2 (((8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))^(3/2)))))/(√(((4 h)/3)-(8 (2^(2/3)) (h^2))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))-(1/(6 (2^(2/3))))((((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+((Sqrt[3]) (d^2))/(√(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))))) √((Sqrt[3]) √(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+6 √(((4 h)/3)-(8 (2^(2/3)) (h^2))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))-(1/(6 (2^(2/3))))((((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+((Sqrt[3]) (d^2))/(√(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))))))+(1/432) (((Sqrt[3]) (8-((32 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3))))+((64 (2^(2/3)) h)/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+((1536 (h^2)+((13824 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(3 (((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3))))))/(2 √(8 h+((32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))+(3 ((4/3)+((8 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(9 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3))))-(768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(18 (2^(2/3)) (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3)))-((16 (2^(2/3)) h)/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))-((Sqrt[3]) (d^2) (8-(32 (2^(2/3)) (h^2) (768 (h^2)+((6912 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)]))))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(4/3)))+(64 (2^(2/3)) h)/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(1536 (h^2)+((13824 (d^4) (h^2))/(Sqrt[81 (d^8)+1536 (d^4) (h^3)])))/(3 (((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(2/3)))))/(2 (((8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))^(3/2)))))/(√(((4 h)/3)-((8 (2^(2/3)) (h^2))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))-((((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))/(6 (2^(2/3))))+((Sqrt[3]) (d^2))/(√(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))))))) ((((Sqrt[3]) √(8 h+((32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))+6 √(((4 h)/3)-(8 (2^(2/3)) (h^2))/(3 (((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))-((((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))/(6 (2^(2/3))))+((Sqrt[3]) (d^2))/(√(8 h+(32 (2^(2/3)) (h^2))/(((27 (d^4)+256 (h^3)+3 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3))+(((54 (d^4)+512 (h^3)+6 (Sqrt[81 (d^8)+1536 (d^4) (h^3)])))^(1/3)))))))^2)) π ((Gamma[(2/5)])^2));

Prepare a large table of data for the optimization fit:

if monitoring is important, implement
Print[{"v0/vc = ",vvc," e= ",e," tf1 = ",tf1}];

Use with fit_poschel.m

[rr0,vvc]=meshgrid(3:(3/39):6,-6:.05:0);

clear all;[tau,gam]=fit_poschel('data','restit_JKR_visco',7.27579E-12,0.0015,7920,1.4E-7,92);

p1=0.0013018987235`;p2=1.16627716109425`;q1=-3.46624928174111`;vvcm[r_]:=((p1 r+p2)/(r+q1));
Hz=(4/5) π (((((Gamma[(2/5)])/(Gamma[(9/10)]))))^2);
v0=((4 (Sqrt[π]) Gamma[(2/5)])/(5 Gamma[(9/10)]));
lnrr0max=6;lnrr0min=3;numk=40;
For[k=1,k≤numk,k=k+1,RR0=(E^(lnrr0min+(((lnrr0max-lnrr0min) (k-1))/(numk-1))));Print[{"k = ",k," RR0 = ",RR0}];tf1=1.6`;minlnv=-6;maxlnv=0;incr=0.05`;numpts=1+((maxlnv-minlnv)/incr);mini=IntegerPart[1+((Log[vvcm[RR0] 1.08`]-minlnv)/incr)];For[i=mini,i≤numpts,i=i+1,vvc=(E^(minlnv+(((i-1) (maxlnv-minlnv))/(numpts-1))));Dvdw=(1/((vvc^(3/5)) (Sqrt[RR0])));Tb=(1/((RR0^(7/6)) (Dvdw^(1/3))));tf=tf1 1.`;solution=NDSolve[{v'[t]==-Re[forceh2[h[t],Dvdw]+Tb v[t] fvisc[h[t],Dvdw]],h'[t]==v[t],v[0]==v0,h[0]==0},{v,h},{t,0,tf}];tf1=Last[Extract[FindRoot[Extract[Re[h[t]]/.solution,1]==hhmin[Dvdw],{t,tf}],1]];e=Extract[Re[-((v[tf1])/(v[0]))]/.solution,1];varr[i]=Log[vvc];earr[i]=e];varray=Table[varr[j],{j,mini,i-1}];earray=Table[earr[j],{j,mini,i-1}];vpad=Table[Log[(E^(minlnv+(((j-1) (maxlnv-minlnv))/(numpts-1))))],{j,1,mini-1}];epad=Table[0,{j,1,mini-1}];eres[k]=Join[epad,earray];result=Table[{Extract[Join[vpad,varray],j],Extract[Join[epad,earray],j]},{j,1,Extract[Dimensions[Join[vpad,varray]],1]}];ListPlot[result,PlotLabel->{"R/R0 = ",RR0},AxesLabel->{"v0/vc","e"},Joined->True,PlotRange->All]];

twhole=Table[{Extract[eres[1],j],Extract[eres[2],j],Extract[eres[3],j],Extract[eres[4],j],Extract[eres[5],j],Extract[eres[6],j],Extract[eres[7],j],Extract[eres[8],j],Extract[eres[9],j],Extract[eres[10],j],Extract[eres[11],j],Extract[eres[12],j],Extract[eres[13],j],Extract[eres[14],j],Extract[eres[15],j],Extract[eres[16],j],Extract[eres[17],j],Extract[eres[18],j],Extract[eres[19],j],Extract[eres[20],j],Extract[eres[21],j],Extract[eres[22],j],Extract[eres[23],j],Extract[eres[24],j],Extract[eres[25],j],Extract[eres[26],j],Extract[eres[27],j],Extract[eres[28],j],Extract[eres[29],j],Extract[eres[30],j],Extract[eres[31],j],Extract[eres[32],j],Extract[eres[33],j],Extract[eres[34],j],Extract[eres[35],j],Extract[eres[36],j],Extract[eres[37],j],Extract[eres[38],j],Extract[eres[39],j],Extract[eres[40],j]},{j,1,Extract[Dimensions[eres[1]],1]}];Export["restit_viscoelastic_JKR.xls",twhole]