Print["\nFunctions dmp[] and prt0[]are equivalent Print[] statements"];

R1=0.795*10^(-3);  (* 1 is sphere *)
R2=Infinity;                 (* 2 is flat *)

R=radius=(1/R1+1/R2)^(-1);  (* Equivalent radius *)
dmp[{R1,R2,R}];

nu1=0.22; C1=1.295*Exp[0.736*nu1]; (* WC indentor - ball *)
nu2=0.30; C2=1.295*Exp[0.736*nu2];  (* Pure nickel - flat *)

E1=475*10^9;
E2=160*10^9;

EE=((1-nu1^2)/E1+(1-nu2^2)/E2)^(-1); (* equivalent E *)
dmp[EE];

Print["\nHardness is NEVER used in the JG model!\n"];
Sy1=1.*10^20; (* WC is a ceramic ; make Sy1 huge so it won't be in effect/used anyway *)
Sy2=159.*10^6; (* pure nickel in paper by Dr. *)

dmp[{nu1, nu2, C1, C2, Sy1, Sy2, C1*Sy1,C2*Sy2}];
CSy=Min[C1*Sy1, C2*Sy2]; dmp[CSy];
If[C1*Sy1<C2*Sy2, CC=C1;Sy=Sy1, CC=C2;Sy=Sy2];
dmp[{CC,Sy, CC*Sy}];

(* JG model *)
omegac=wc=(Pi*CSy/(2*EE))^2*R;
Pc=(4/3)*(R/EE)^2(Pi*CSy/2)^3;
dmp[wcjg=omegac];
dmp[Pcjg=Pc];

(*JG model - force vs. disp. *)
wstar=w/wcjg;
PF1jg=wstar^(3/2);dmp[PF1jg];  (* Hertzian solution - it stretches accurately to wstar=1.9 *)

ey=Sy/EE;
B=0.14*Exp[23*ey];dmp[B];
aoR=0.5*Pi*CC*ey*Sqrt[wstar*(wstar/1.9)^B]; dmp[aoR];

HGSyorg=2.84*(1-Exp[-0.82*aoR^(-0.7)]); dmp[HGSyorg];   (* original paper by JG *)
HGSy=2.84-0.92*(1-Cos[Pi*aoR]); dmp[HGSy];  (* newer expression, from COR paper to handle extreme deformation and palstcity *)

orgHGSy=0*1; (* use 0 to apply newer value; use 1 to apply the original orger expression; BOTH are very close anyways in up to moderate wstar*)
If[orgHGSy==1,HGSy=HGSyorg];

wmax=N[25*10^(-6)*1]; dmp[wmax]; (* maximum wstar *)

Print["\nHG/Sy"]
Plot[{HGSy,HGSyorg},{w,0,wmax}, GridLines->Automatic,  PlotLegends->{"HG/Sy newer","HG/Sy original "}]

PF2jg=Exp[-0.25*wstar^(5/12)]*wstar^(3/2)+(4*HGSy/CC)*(1-Exp[-wstar^(5/9)/25])*wstar; dmp[PF2jg];

pfjg=Pcjg*Which[w<1.9*wcjg,PF1jg,
                                  w≥1.9*wcjg,PF2jg];
dmp[pfjg];
dmp[Pcjg*PF1jg /. w-> wcjg]
dmp[Pcjg*PF2jg /. w->wmax]

ww=Table[w,       {w,0,wmax,wmax/100}]*1.*10^6;
ffjg=Table[pfjg,{w,0,wmax,wmax/100}];
dmp[{Min[ww],Max[ww],Min[ffjg],Max[ffjg]}];

dmp[wmax/wcjg]

Print["\na/R"]
Plot[aoR,{w,0,wmax}, GridLines->Automatic]

Print["\nForce vs. disp.: Linear-Linear"]
Plot[{pfjg},{w,0,wmax}, GridLines->Automatic, PlotLegends->{"Jackson-Green,JG"}]
Print["\nForce vs. disp.: Log-Log"]
LogLogPlot[{pfjg},{w,wcjg/10.0,wmax}, GridLines->Automatic,PlotLegends->{"Jackson-Green,JG"}]
Print["\nForce vs. disp.: Linear-Log"]
LogPlot[{pfjg},{w,0,wmax}, GridLines->Automatic,PlotLegends->{"Jackson-Green,JG"}]
Print["\nForce vs. disp.: Log-Linear"]
LogLinearPlot[{pfjg},{w,wcjg/10.0,wmax}, GridLines->Automatic, PlotLegends->{"Jackson-Green,JG"}]