| 1215 | # kk: list of curvature fields |
| 1216 | # ll: list of labels |
| 1217 | def FigAng2(vf, kk, ll, po): |
| 1218 | dim = GetDim(vf.shape) |
| 1219 | x1, y1, z1, hx, hy, hz = GetGeom(vf.shape) |
| 1220 | x, y, z = GetMesh(x1, y1, z1) |
| 1221 | |
| 1222 | # interface cells |
| 1223 | th = 0 |
| 1224 | ii = np.where((vf > th) & (vf < 1. - th)) |
| 1225 | x = x[ii] |
| 1226 | y = y[ii] |
| 1227 | z = z[ii] |
| 1228 | |
| 1229 | # hires angle |
| 1230 | anh = np.linspace(-np.pi, np.pi, 200) |
| 1231 | degh = np.degrees(anh) |
| 1232 | # hires points |
| 1233 | cx, cy, cz, rx, ry, rz = LoadBub() |
| 1234 | r = (rx * ry) / ((ry * np.cos(anh))**2 + (rx * np.sin(anh))**2)**0.5 |
| 1235 | xh = cx + r * np.cos(anh) |
| 1236 | yh = cy + r * np.sin(anh) |
| 1237 | # hires curvature |
| 1238 | keh = GetExactK(dim, xh, yh, cz) |
| 1239 | # average curvature |
| 1240 | kea = keh.mean() |
| 1241 | |
| 1242 | dx = x - cx |
| 1243 | dy = y - cy |
| 1244 | an = np.arctan2(dy, dx) |
| 1245 | deg = np.degrees(an) |
| 1246 | s = np.argsort(an) |
| 1247 | |
| 1248 | fig, ax = PlotInit() |
| 1249 | |
| 1250 | # plot estimates |
| 1251 | for k, l in zip(kk, ll): |
| 1252 | if k is not None: |
| 1253 | ax.plot(deg[s], k[ii][s] / kea, label=l) |
| 1254 | |
| 1255 | # plot lowres exact curvature |
| 1256 | ax.plot(deg[s], GetExactK(dim, x, y, z)[s] / kea, label="exact") |
| 1257 | |
| 1258 | # plot exact curvature |
| 1259 | ax.plot(degh, keh / kea, label="exact", c="0.5", ls='--') |
| 1260 | |
| 1261 | ax.set_xlabel(r"angle [deg]") |
| 1262 | ax.set_ylabel(r"normalized curvature") |
| 1263 | q = 180. |
| 1264 | ax.set_xticks(np.arange(-q, q * 1.1, 90.)) |
| 1265 | ax.set_xlim(-q, q) |
| 1266 | ax.legend() |
| 1267 | ax.set_ylim(0., 2.) |
| 1268 | ax.grid() |
| 1269 | PlotSave(fig, ax, po) |
| 1270 | |
| 1271 | |
| 1272 | def Univel(): |