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