1
$\begingroup$

I tried to calculate mean curvature of laser induced axisymmetric deformed surface, which radial dependence of height is given by:

γw = 70.0*10^-3; (*SR in N/m*)
ρ = 1000;
c = 3*10^8;
g = 9.8;
size = 600;
u1 = γw*k^2 + ρ*g;
we1 = 7*10^-6; (*beam waist*)
n = 1.33;
P0 = 4.0;
rng = 120*10^-6; P1 = (P0/(c*Pi))*((n - 1)/(n + 1));

f70[r_] := P1*Exp[-(we1^2*k^2)/8]*k*BesselJ[0, r*k]/u1
h70[r_?NumericQ] := 
 Quiet[NIntegrate[f70[r], {k, 0, ∞}, 
   Method -> {"GlobalAdaptive"}, MinRecursion -> 4, 
   MaxRecursion -> 100]]

hp70 = Plot[10^9*(h70[r]), {r, -rng, rng}, Axes -> False, 
  PlotStyle -> {Blue, Thickness[0.006]}, PlotRange -> All, 
  ImageSize -> size, LabelStyle -> Directive[Blue, 18], Frame -> True]
$\endgroup$
2
  • $\begingroup$ Can you be more explicit in your question? What does the plot represent? Do you have a definition of the quantity you are trying to calculate? What is the relationship between the variables you introduced and that quantity? $\endgroup$
    – MarcoB
    Mar 1 at 13:23
  • $\begingroup$ Thanks @Marcob, plot represent the height as function of radial coordinate. I want to calculate principal curvatures so that I can calculate the mean curvature. $\endgroup$ Mar 1 at 14:30

1 Answer 1

1
$\begingroup$

The mean curvature in a case of axial symmetric surface h[y] can be defined in a form

H[y_] := (h''[y]/(1 + h'[y]^2)^(3/2) + h'[y]/Sqrt[1 + h'[y]^2]/y)/2;

Using numeric integral we can compute interpolation function described surface h as follows

\[Gamma]w = 70.0*10^-3;(*SR in N/m*)\[Rho] = 1000;
c = 3*10^8;
g = 9.8;
size = 600;
u1 = \[Gamma]w*k^2 + \[Rho]*g;
we1 = 7*10^-6;(*beam waist*)n = 1.33;
P0 = 4.0;
rng = 120*10^-6; P1 = (P0/(c*Pi))*((n - 1)/(n + 1));

f70[r_] := P1*Exp[-(we1^2*k^2)/8]*k*BesselJ[0, r*k]/u1
h70[r_?NumericQ] := 
 Quiet[NIntegrate[f70[r], {k, 0, \[Infinity]}, 
   Method -> {"GlobalAdaptive"}, MinRecursion -> 4, 
   MaxRecursion -> 100]]
lst = Table[{r, h70[r]}, {r, 10^-6, rng, 10^-6}];

h = Interpolation[lst, InterpolationOrder -> 4]; 

At r->0 function H[r] has singularity, therefore we need to regularize h by using series

f7 = Series[BesselJ[0, r*k], {r, 0, 7}] // Normal

(*Out[]= 1 - (k^2 r^2)/4 + (k^4 r^4)/64 - (k^6 r^6)/2304*)

 h7 = 
 Integrate[P1*Exp[-(we1^2*k^2)/8]*k*f7/u1, {k, 0, Infinity}, 
  Assumptions -> r > 0]

(*Out[]= 5.74997*10^-8 - 175.246 r^2 + 1.78824*10^12 r^4 - 
 1.62199*10^22 r^6 *)

Finally we define H7 and point r=rmin to join H7 and H

H7 = (D[h7, {r, 2}]/(1 + D[h7, r]^2)^(3/2) + 
    D[h7, r]/Sqrt[1 + D[h7, r]]/r)/2

 rmin = 
 r /. FindMinimum[Abs[H7 - H[r]], {r, 2 10^-6, 10^-6, 10^-5}][[2]] // 
  Quiet

(*Out[]= 1.79474*10^-6 *)

Visualization

{Plot[If[r < rmin , H7, H[r]], {r, 0, 10^-5}, PlotRange -> All, 
  PlotPoints -> 200], 
 Plot[If[r < rmin, H7, H[r]], {r, 0, rng}, PlotRange -> All, 
  PlotPoints -> 200]}

Figure 1

$\endgroup$
5
  • $\begingroup$ Thanks@Alex, It is working. Please also suggest that Can I straight forward calculate the mean curvature derivative? What about the small slope limit approximation? $\endgroup$ Mar 3 at 4:06
  • $\begingroup$ @GopalVerma It could be better to use piecewise function If[r < rmin , H7, H[r]] for approximation mean curvature only. For derivative we need some smooth differentiable function based on it. For instance, we can prepare new function as follows Hnew = Interpolation[ Table[{r, If[r < rmin, H7, H[r]]}, {r, 0, rng, rng/1000}], InterpolationOrder -> 3] $\endgroup$ Mar 3 at 5:26
  • $\begingroup$ Thanks@Alex, It is working. $\endgroup$ Mar 3 at 14:33
  • $\begingroup$ This code Hnew = Interpolation[ Table[{r, If[r < rmin, H7, H[r]]}, {r, 0, rng, rng/1000}], InterpolationOrder -> 3] showing 1/0 expression, but if I replaced 0 with small number e.g. 10^-7,then we can calculate derivative and plot it, but plotted graph is not smooth and does not have negative values (it should have in some region). Please suggest full code regarding derivative. $\endgroup$ Mar 4 at 3:07
  • $\begingroup$ @GopalVerma To avoid singularity 1/0 we can use H7 = (D[h7, {r, 2}]/(1 + D[h7, r]^2)^(3/2) + D[h7, r]/Sqrt[1 + D[h7, r]]/r)/2 // FullSimplify $\endgroup$ Mar 4 at 5:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.