src.acoustools.BEM.Gorkov

 1import torch
 2from torch import Tensor
 3
 4from vedo import Mesh
 5
 6from acoustools.Utilities import TRANSDUCERS, device
 7from acoustools.Mesh import load_scatterer
 8
 9from acoustools.BEM.Forward_models import compute_E
10from acoustools.BEM.Gradients import BEM_forward_model_grad
11
12import acoustools.Constants as Constants
13
14
15
16
17def BEM_gorkov_analytical(activations:Tensor,points:Tensor,scatterer:Mesh|None|str=None,
18                          board:Tensor|None=None,H:Tensor|None=None,E:Tensor|None=None, return_components:bool=False, dims='XYZ',
19                          V:float=Constants.V,
20                          **params) -> Tensor:
21    '''
22    Returns Gor'kov potential computed analytically from the BEM model\n
23    :param activations: Transducer hologram
24    :param points: Points to propagate to
25    :param scatterer: The mesh used (as a `vedo` `mesh` object) or string of path to mesh
26    :param board: Transducers to use 
27    :param H: Precomputed H - if None H will be computed
28    :param E: Precomputed E - if None E will be computed
29    :param return_components: 
30    :param dims: Dimensions to consider gradient in
31    :param V: Volume of particles
32    :return: Gor'kov potential at point U
33    '''
34    if board is None:
35        board = TRANSDUCERS
36    if type(scatterer) == str:
37            scatterer = load_scatterer(scatterer)
38    
39    path = params['path']
40    
41    if E is None:
42        E = compute_E(scatterer,points,board,H=H,path=path)
43
44    Ex, Ey, Ez = BEM_forward_model_grad(points,scatterer,board,H=H,path=path)
45
46    p = E@activations
47    
48    if 'X' in dims.upper():
49        px = Ex@activations
50    else:
51        px = torch.Tensor((0,)).to(device)
52    
53    if 'Y' in dims.upper():
54        py = Ey@activations
55    else:
56        py = torch.Tensor((0,)).to(device)
57    
58    if 'Z' in dims.upper():
59        pz = Ez@activations
60    else:
61        pz = torch.Tensor((0,)).to(device)
62    
63    K1 = V / (4*Constants.p_0*Constants.c_0**2)
64    K2 = 3*V / (4*(2*Constants.f**2 * Constants.p_0))
65
66    a = K1 * torch.abs(p)**2 
67    b = K2*(torch.abs(px)**2 + torch.abs(py)**2 + torch.abs(pz)**2)
68
69    U = a-b
70
71    if return_components:
72        return U, a ,b
73
74    return U