1 |
2 (* Clean register. *) |
3 Clear[“Global’*”] |
4 |
5 (* Define first surface function *) |
6 za[ra_] := Module[{h}, h = Cos[ra 0.5]]; |
7 |
8 (* Define variables of Eq. (5) *) |
9 Phi[ra_] := Module[{h}, |
10 h = Sqrt[1 − (ra + (−ta + za[ra]) Derivative[1][za][ra])^2/ |
11 (n^2 (ra^2 + (ta − za[ra])^2) (1 + (Derivative[1][za][ra])^2))]/ |
12 Sqrt[1 + (Derivative[1][za][ra])^2]] |
13 |
14 ri[ra_] := Module[{h}, |
15 h = (ra + (−ta + za[ra]) Derivative[1][za][ra])/( |
16 n Sqrt[ra^2 + (ta − za[ra])^2] (1 + (Derivative[1][za][ra])^2)) − |
17 Derivative[1][za][ra] Phi[ra]] |
18 |
19 zi[ra_] := Module[{h}, |
20 h = (Derivative[1][za][ |
21 ra] (ra + (−ta + za[ra]) Derivative[1][za][ra]))/( |
22 n Sqrt[ra^2 + (ta − |
23 za[ra])^2] (1 + (Derivative[1][za][ra])^2)) +Phi[ |
24 ra]]; |
25 |
26 fi[ra_] := Module[{h}, h = ta − tb − Sign[ta] Sqrt[ra^2 + (ta − za[ra])^2]]; |
27 |
28 h0[ra_] := Module[{h}, |
29 h = ri[ra]^2 za[ra] + fi[ra] n zi[ra] − |
30 ra ri[ra] zi[ra] + (t + tb) zi[ra]^2 − n^2 (za[ra] + t zi[ra])]; |
31 |
32 h1[ra_] := Module[{h}, |
33 h = ra^2 + 2 ra ri[ra] t + (tb − za[ra])^2 + |
34 t^2 (ri[ra]^2 + (−1 + zi[ra])^2) − |
35 2 t (tb − za[ra]) (−1 + zi[ra])]; |
36 |
37 zb[ra_] := Module[{h}, |
38 h = (h0[ra] + |
39 s1 Sqrt[zi[ra]^2 (fi[ra]^2 − |
40 2 fi[ra] n (ra ri[ra] + ri[ra]^2 t + |
41 zi[ra] (t (zi[ra] − 1) − tb + za[ra])) + |
42 h1[ra] n^2 − (ra zi[ra] + |
43 ri[ra] (t + tb − za[ra]))^2)])/(−n^2 + 1)]; |
44 |
45 rb[ra_] := Module[{h}, h = ra + (ri[ra] (−za[ra] + zb[ra]))/zi[ra]]; |
46 |
47 (* Define system configuration’s *) |
48 s1 = 1; |
49 n = 1.5; |
50 t = 8; |
51 ta = −60; |
52 tb = 70; |
53 basez = ta; |
54 topez = tb + t; |
55 rmax = 16.8; (* Iterative value *) |
56 |
57 (* System plot *) |
58 g1 = Plot[{za[ra], topez, basez}, {ra, −rmax, rmax}, |
59 PlotStyle −> {Black, Transparent, Transparent}, |
60 AspectRatio −> Automatic, AxesLabel −> {r, z}]; |
61 g2 = ParametricPlot[{rb[ra], zb[ra]}, {ra, −rmax, rmax}, |
62 AspectRatio −> 1/1, AxesOrigin −> {0, 0}, PlotStyle −> {Black}, |
63 Mesh −> None]; |
64 Show[g1, g2] |