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]
Improve this question
$\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, 2022 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$
    – Gopal Verma
    Mar 1, 2022 at 14:30

1 Answer 1

Reset to default
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

Improve this answer
$\endgroup$
12
  • $\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$
    – Gopal Verma
    Mar 3, 2022 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$
    – Alex Trounev
    Mar 3, 2022 at 5:26
  • $\begingroup$ Thanks@Alex, It is working. $\endgroup$
    – Gopal Verma
    Mar 3, 2022 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$
    – Gopal Verma
    Mar 4, 2022 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$
    – Alex Trounev
    Mar 4, 2022 at 5:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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