src.acoustools.Gorkov 

  1import torch
  2from acoustools.Utilities import device, propagate, forward_model_batched, forward_model_grad, TRANSDUCERS, DTYPE
  3import acoustools.Constants as c
  4
  5from torch import Tensor
  6from types import FunctionType
  7
  8
  9def gorkov(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
 10                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, **params) -> Tensor:
 11    '''
 12    Use to compute the gorkov potential at a point. Alias for `src.acoustools.Gorkov.gorkov_analytical`
 13    '''
 14    return gorkov_analytical(activations, points, board, axis, V=V, p_ref=p_ref, k=k, transducer_radius=transducer_radius, 
 15                             medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density,particle_speed=particle_speed , **params)
 16
 17def gorkov_analytical(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", board_norms = None,
 18                        V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, 
 19                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p , **params) -> Tensor:
 20    '''
 21    Computes the Gorkov potential using analytical derivative of the piston model \n
 22    :param activation: The transducer activations to use 
 23    :param points: The points to compute the potential at
 24    :param board: The transducer boards to use
 25    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
 26    :param V: particle Volume
 27    :return: gorkov potential at each point
 28    ```Python
 29    from acoustools.Utilities import create_points, add_lev_sig
 30    from acoustools.Solvers import wgs
 31    from acoustools.Gorkov import gorkov_analytical
 32
 33    N=1
 34    B=1
 35    points = create_points(N,B)
 36    x = wgs(points)
 37    x = add_lev_sig(x)
 38    
 39    U_a = gorkov_analytical(x,points)
 40
 41    print("Analytical",U_a.data.squeeze())
 42    ```
 43    '''
 44
 45    if board is None:
 46        board = TRANSDUCERS
 47
 48    Fx, Fy, Fz = forward_model_grad(points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius)
 49    F = forward_model_batched(points,board, p_ref=p_ref, k=k, transducer_radius=transducer_radius)
 50    
 51    p = torch.abs(F@activations)**2
 52
 53    px = (Fx@activations) if 'X' in axis else 0
 54    py = (Fy@activations) if 'Y' in axis else 0
 55    pz = (Fz@activations) if 'Z' in axis else 0
 56
 57    grad  = torch.cat((px,py,pz),dim=2)
 58
 59    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
 60    g = (torch.sum(torch.abs(grad)**2, dim=2, keepdim=True))
 61    
 62    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
 63    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
 64    
 65    U = K1*p - K2*g
 66    return U
 67
 68
 69def get_gorkov_constants(V=c.V, p_0 = c.p_0, p_p=c.p_p, c_0=c.c_0, c_p=c.c_p, angular_frequency=c.angular_frequency ):
 70    '''
 71    Returns K1 and K2 for use in Gorkov computations, Uses the form shown in `Holographic acoustic elements for manipulation of levitated objects` \n
 72    :param: V: Particle Volume
 73    :param p_0: Density of medium
 74    :param p_p: Density of particle
 75    :param c_0: Speed of sound in medium
 76    :param c_p: speed of sound in particle
 77    :param angular_frequency: The angular frequency
 78    :returns K1, K2:
 79
 80    '''
 81    #Derived Bk.3 Pg.91
 82    K1 = 1/4*V*(1/(c_0**2*p_0) - 1/(c_p**2*p_p)) 
 83    K2 = 3/4 * V * ((p_p - p_0) / (angular_frequency**2 * p_0 * (p_0 + 2*p_p))) 
 84
 85    # K1 = V / (4*p_0*c_0**2)
 86    # K2 = 3*V / (4*(2*f**2 * p_0))
 87
 88    return K1, K2
 89
 90def gorkov_autograd(activations:Tensor, points:Tensor, K1:float|None=None, K2:float|None=None, V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
 91                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,
 92                     retain_graph:bool=False,board:Tensor|None=None, board_norms = None,**params) -> Tensor:
 93    '''
 94    Computes the Gorkov potential using pytorch's autograd system\n
 95    :param activation: The transducer activations to use 
 96    :param points: The points to compute the potential at. if `None` will use `acoustools.Utilities.TRANSDUCERS`
 97    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
 98    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
 99    :param board: The transducer boards to use
100    :param retain_graph: Value will be passed to autograd
101    :return: gorkov potential at each point
102
103    ```Python
104    from acoustools.Utilities import create_points, add_lev_sig
105    from acoustools.Solvers import wgs
106    from acoustools.Gorkov import gorkov_autograd
107
108    N=1
109    B=1
110    points = create_points(N,B)
111    x = wgs(points)
112    x = add_lev_sig(x)
113    
114    U_ag = gorkov_autograd(x,points)
115
116    print("Autograd", U_ag.data.squeeze())
117    ```
118    '''
119
120    if board is None:
121        board = TRANSDUCERS
122
123    var_points = torch.autograd.Variable(points.data, requires_grad=True).to(device).to(DTYPE)
124
125    B = points.shape[0]
126    N = points.shape[2]
127    
128    if len(activations.shape) < 3:
129        activations.unsqueeze_(0)    
130    
131    pressure = propagate(activations.to(DTYPE),var_points,board=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=board_norms)
132    pressure.backward(torch.ones((B,N))+0j, inputs=var_points, retain_graph=retain_graph)
133    grad_pos = var_points.grad
134
135    if K1 is None or K2 is None:
136        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
137        if K1 is None:
138            K1 = K1_
139        if K2 is None:
140            K2 = K2_
141
142    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
143    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
144
145
146    gorkov = K1 * torch.abs(pressure) **2 - K2 * torch.sum((torch.abs(grad_pos)**2),1)
147    return gorkov
148
149
150def gorkov_fin_diff(activations: Tensor, points:Tensor, axis:str="XYZ", stepsize:float = 0.000135156253,K1:float|None=None, K2:float|None=None,
151                    prop_function:FunctionType=propagate,prop_fun_args:dict={}, board:Tensor|None=None, board_norms=None,
152                    V=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
153                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,) -> Tensor:
154    '''
155    Computes the Gorkov potential using finite differences to compute derivatives \n
156    :param activation: The transducer activations to use 
157    :param points: The points to compute the potential at
158    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
159    :param stepsize: The distance aroud points to add, default 0.000135156253
160    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
161    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
162    :param prop_function: Function to use to compute pressure
163    :param prop_fun_args: Arguments to pass to `prop_function`
164    :param board: The transducer boards to use if `None` use `acoustools.Utilities.TRANSDUCERS`
165    :param V: Particle Volume
166    :return: gorkov potential at each point
167
168    ```Python
169    from acoustools.Utilities import create_points, add_lev_sig
170    from acoustools.Solvers import wgs
171    from acoustools.Gorkov import gorkov_fin_diff
172
173    N=1
174    B=1
175    points = create_points(N,B)
176    x = wgs(points)
177    x = add_lev_sig(x)
178    
179    U_fd = gorkov_fin_diff(x,points)
180
181    print("Finite Differences",U_fd.data.squeeze())
182    ```
183    '''
184    # torch.autograd.set_detect_anomaly(True)
185    if board is None:
186        board = TRANSDUCERS
187    B = points.shape[0]
188    D = len(axis)
189    N = points.shape[2]
190
191    
192    if len(activations.shape) < 3:
193        activations = torch.unsqueeze(activations,0).clone().to(device)
194
195    fin_diff_points = get_finite_diff_points_all_axis(points, axis, stepsize)
196
197    pressure_points = prop_function(activations, fin_diff_points,board=board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms = board_norms,**prop_fun_args)
198    # if len(pressure_points.shape)>1:
199    # pressure_points = torch.squeeze(pressure_points,2)
200
201    pressure = pressure_points[:,:N]
202    pressure_fin_diff = pressure_points[:,N:]
203
204    split = torch.reshape(pressure_fin_diff,(B,2, ((2*D))*N // 2))
205    
206    grad = (split[:,0,:] - split[:,1,:]) / (2*stepsize)
207    
208    grad = torch.reshape(grad,(B,D,N))
209    grad_abs_square = torch.pow(torch.abs(grad),2)
210    grad_term = torch.sum(grad_abs_square,dim=1)
211
212    if K1 is None or K2 is None:
213        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
214        if K1 is None:
215            K1 = K1_
216        if K2 is None:
217            K2 = K2_
218    
219    # p_in =  torch.abs(pressure)
220    p_in = torch.sqrt(torch.real(pressure) **2 + torch.imag(pressure)**2)
221    if len(p_in.shape) > 2:
222        p_in.squeeze_(2)
223    # p_in = torch.squeeze(p_in,2)
224
225    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
226    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
227
228    U = K1 * p_in**2 - K2 *grad_term
229    
230    return U
231
232
233def get_finite_diff_points(points:Tensor , axis:Tensor, stepsize:float = 0.000135156253) -> Tensor:
234    '''
235    Gets points for finite difference calculations in one axis\n
236    :param points: Points around which to find surrounding points
237    :param axis: The axis to add points in
238    :param stepsize: The distance aroud points to add, default 0.000135156253
239    :return: points 
240    '''
241    #points = Bx3x4
242    points_h = points.clone()
243    points_neg_h = points.clone()
244    points_h[:,axis,:] = points_h[:,axis,:] + stepsize
245    points_neg_h[:,axis,:] = points_neg_h[:,axis,:] - stepsize
246
247    return points_h, points_neg_h
248
249def get_finite_diff_points_all_axis(points: Tensor,axis: str="XYZ", stepsize:float = 0.000135156253) -> Tensor:
250    '''
251    Gets points for finite difference calculations\\
252    :param points: Points around which to find surrounding points\\
253    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\\
254    :param stepsize: The distance aroud points to add, default 0.000135156253\\
255    :return: Points
256    '''
257    B = points.shape[0]
258    D = len(axis)
259    N = points.shape[2]
260    fin_diff_points=  torch.zeros((B,3,((2*D)+1)*N)).to(device).to(DTYPE)
261    fin_diff_points[:,:,:N] = points.clone()
262
263    i = 2
264    if "X" in axis:
265        points_h, points_neg_h = get_finite_diff_points(points, 0, stepsize)
266        fin_diff_points[:,:,N:i*N] = points_h
267        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
268
269        i += 1
270
271    
272    if "Y" in axis:
273        points_h, points_neg_h = get_finite_diff_points(points, 1, stepsize)
274        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
275        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
276        i += 1
277    
278    if "Z" in axis:
279        points_h, points_neg_h = get_finite_diff_points(points, 2, stepsize)
280        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
281        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
282        i += 1
283    
284    return fin_diff_points
def gorkov( activations: torch.Tensor, points: torch.Tensor, board: torch.Tensor | None = None, axis: str = 'XYZ', V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, **params) -> torch.Tensor: 
10def gorkov(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
11                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p, **params) -> Tensor:
12    '''
13    Use to compute the gorkov potential at a point. Alias for `src.acoustools.Gorkov.gorkov_analytical`
14    '''
15    return gorkov_analytical(activations, points, board, axis, V=V, p_ref=p_ref, k=k, transducer_radius=transducer_radius, 
16                             medium_density=medium_density, medium_speed=medium_speed, particle_density=particle_density,particle_speed=particle_speed , **params)

Use to compute the gorkov potential at a point. Alias for src.acoustools.Gorkov.gorkov_analytical

def gorkov_analytical( activations: torch.Tensor, points: torch.Tensor, board: torch.Tensor | None = None, axis: str = 'XYZ', board_norms=None, V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, **params) -> torch.Tensor: 
18def gorkov_analytical(activations: Tensor, points: Tensor,board:Tensor|None=None, axis:str="XYZ", board_norms = None,
19                        V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius, 
20                        medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p , **params) -> Tensor:
21    '''
22    Computes the Gorkov potential using analytical derivative of the piston model \n
23    :param activation: The transducer activations to use 
24    :param points: The points to compute the potential at
25    :param board: The transducer boards to use
26    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
27    :param V: particle Volume
28    :return: gorkov potential at each point
29    ```Python
30    from acoustools.Utilities import create_points, add_lev_sig
31    from acoustools.Solvers import wgs
32    from acoustools.Gorkov import gorkov_analytical
33
34    N=1
35    B=1
36    points = create_points(N,B)
37    x = wgs(points)
38    x = add_lev_sig(x)
39    
40    U_a = gorkov_analytical(x,points)
41
42    print("Analytical",U_a.data.squeeze())
43    ```
44    '''
45
46    if board is None:
47        board = TRANSDUCERS
48
49    Fx, Fy, Fz = forward_model_grad(points, transducers=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius)
50    F = forward_model_batched(points,board, p_ref=p_ref, k=k, transducer_radius=transducer_radius)
51    
52    p = torch.abs(F@activations)**2
53
54    px = (Fx@activations) if 'X' in axis else 0
55    py = (Fy@activations) if 'Y' in axis else 0
56    pz = (Fz@activations) if 'Z' in axis else 0
57
58    grad  = torch.cat((px,py,pz),dim=2)
59
60    K1, K2 = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
61    g = (torch.sum(torch.abs(grad)**2, dim=2, keepdim=True))
62    
63    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
64    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
65    
66    U = K1*p - K2*g
67    return U

Computes the Gorkov potential using analytical derivative of the piston model

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at
  • board: The transducer boards to use
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
  • V: particle Volume
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_analytical

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_a = gorkov_analytical(x,points)

print("Analytical",U_a.data.squeeze())
def get_gorkov_constants( V=4.188790204666667e-09, p_0=1.2, p_p=29.36, c_0=343, c_p=1052, angular_frequency=251327.41228000002): 
70def get_gorkov_constants(V=c.V, p_0 = c.p_0, p_p=c.p_p, c_0=c.c_0, c_p=c.c_p, angular_frequency=c.angular_frequency ):
71    '''
72    Returns K1 and K2 for use in Gorkov computations, Uses the form shown in `Holographic acoustic elements for manipulation of levitated objects` \n
73    :param: V: Particle Volume
74    :param p_0: Density of medium
75    :param p_p: Density of particle
76    :param c_0: Speed of sound in medium
77    :param c_p: speed of sound in particle
78    :param angular_frequency: The angular frequency
79    :returns K1, K2:
80
81    '''
82    #Derived Bk.3 Pg.91
83    K1 = 1/4*V*(1/(c_0**2*p_0) - 1/(c_p**2*p_p)) 
84    K2 = 3/4 * V * ((p_p - p_0) / (angular_frequency**2 * p_0 * (p_0 + 2*p_p))) 
85
86    # K1 = V / (4*p_0*c_0**2)
87    # K2 = 3*V / (4*(2*f**2 * p_0))
88
89    return K1, K2

Returns K1 and K2 for use in Gorkov computations, Uses the form shown in Holographic acoustic elements for manipulation of levitated objects

Parameters
  • V: Particle Volume
  • p_0: Density of medium
  • p_p: Density of particle
  • c_0: Speed of sound in medium
  • c_p: speed of sound in particle
  • angular_frequency: The angular frequency :returns K1, K2:
def gorkov_autograd( activations: torch.Tensor, points: torch.Tensor, K1: float | None = None, K2: float | None = None, V: float = 4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052, retain_graph: bool = False, board: torch.Tensor | None = None, board_norms=None, **params) -> torch.Tensor: 
 91def gorkov_autograd(activations:Tensor, points:Tensor, K1:float|None=None, K2:float|None=None, V:float=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
 92                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,
 93                     retain_graph:bool=False,board:Tensor|None=None, board_norms = None,**params) -> Tensor:
 94    '''
 95    Computes the Gorkov potential using pytorch's autograd system\n
 96    :param activation: The transducer activations to use 
 97    :param points: The points to compute the potential at. if `None` will use `acoustools.Utilities.TRANSDUCERS`
 98    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
 99    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
100    :param board: The transducer boards to use
101    :param retain_graph: Value will be passed to autograd
102    :return: gorkov potential at each point
103
104    ```Python
105    from acoustools.Utilities import create_points, add_lev_sig
106    from acoustools.Solvers import wgs
107    from acoustools.Gorkov import gorkov_autograd
108
109    N=1
110    B=1
111    points = create_points(N,B)
112    x = wgs(points)
113    x = add_lev_sig(x)
114    
115    U_ag = gorkov_autograd(x,points)
116
117    print("Autograd", U_ag.data.squeeze())
118    ```
119    '''
120
121    if board is None:
122        board = TRANSDUCERS
123
124    var_points = torch.autograd.Variable(points.data, requires_grad=True).to(device).to(DTYPE)
125
126    B = points.shape[0]
127    N = points.shape[2]
128    
129    if len(activations.shape) < 3:
130        activations.unsqueeze_(0)    
131    
132    pressure = propagate(activations.to(DTYPE),var_points,board=board, p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms=board_norms)
133    pressure.backward(torch.ones((B,N))+0j, inputs=var_points, retain_graph=retain_graph)
134    grad_pos = var_points.grad
135
136    if K1 is None or K2 is None:
137        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
138        if K1 is None:
139            K1 = K1_
140        if K2 is None:
141            K2 = K2_
142
143    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
144    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
145
146
147    gorkov = K1 * torch.abs(pressure) **2 - K2 * torch.sum((torch.abs(grad_pos)**2),1)
148    return gorkov

Computes the Gorkov potential using pytorch's autograd system

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at. if None will use acoustools.Utilities.TRANSDUCERS
  • K1: The value for K1 in the Gorkov equation, if None will use c.V / (4*c.p_0*c.c_0**2)
  • K2: The value for K2 in the Gorkov equation, if None will use 3*c.V / (4*(2*c.f**2 * c.p_0))
  • board: The transducer boards to use
  • retain_graph: Value will be passed to autograd
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_autograd

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_ag = gorkov_autograd(x,points)

print("Autograd", U_ag.data.squeeze())
def gorkov_fin_diff( activations: torch.Tensor, points: torch.Tensor, axis: str = 'XYZ', stepsize: float = 0.000135156253, K1: float | None = None, K2: float | None = None, prop_function: function = <function propagate>, prop_fun_args: dict = {}, board: torch.Tensor | None = None, board_norms=None, V=4.188790204666667e-09, p_ref=3.4000000000000004, k=732.7329804081634, transducer_radius=0.0045, medium_density=1.2, medium_speed=343, particle_density=29.36, particle_speed=1052) -> torch.Tensor: 
151def gorkov_fin_diff(activations: Tensor, points:Tensor, axis:str="XYZ", stepsize:float = 0.000135156253,K1:float|None=None, K2:float|None=None,
152                    prop_function:FunctionType=propagate,prop_fun_args:dict={}, board:Tensor|None=None, board_norms=None,
153                    V=c.V, p_ref=c.P_ref, k=c.k, transducer_radius = c.radius,
154                     medium_density=c.p_0, medium_speed = c.c_0, particle_density = c.p_p, particle_speed = c.c_p,) -> Tensor:
155    '''
156    Computes the Gorkov potential using finite differences to compute derivatives \n
157    :param activation: The transducer activations to use 
158    :param points: The points to compute the potential at
159    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
160    :param stepsize: The distance aroud points to add, default 0.000135156253
161    :param K1: The value for K1 in the Gorkov equation, if `None` will use `c.V / (4*c.p_0*c.c_0**2)`
162    :param K2: The value for K2 in the Gorkov equation, if `None` will use `3*c.V / (4*(2*c.f**2 * c.p_0))`
163    :param prop_function: Function to use to compute pressure
164    :param prop_fun_args: Arguments to pass to `prop_function`
165    :param board: The transducer boards to use if `None` use `acoustools.Utilities.TRANSDUCERS`
166    :param V: Particle Volume
167    :return: gorkov potential at each point
168
169    ```Python
170    from acoustools.Utilities import create_points, add_lev_sig
171    from acoustools.Solvers import wgs
172    from acoustools.Gorkov import gorkov_fin_diff
173
174    N=1
175    B=1
176    points = create_points(N,B)
177    x = wgs(points)
178    x = add_lev_sig(x)
179    
180    U_fd = gorkov_fin_diff(x,points)
181
182    print("Finite Differences",U_fd.data.squeeze())
183    ```
184    '''
185    # torch.autograd.set_detect_anomaly(True)
186    if board is None:
187        board = TRANSDUCERS
188    B = points.shape[0]
189    D = len(axis)
190    N = points.shape[2]
191
192    
193    if len(activations.shape) < 3:
194        activations = torch.unsqueeze(activations,0).clone().to(device)
195
196    fin_diff_points = get_finite_diff_points_all_axis(points, axis, stepsize)
197
198    pressure_points = prop_function(activations, fin_diff_points,board=board,p_ref=p_ref, k=k, transducer_radius=transducer_radius, norms = board_norms,**prop_fun_args)
199    # if len(pressure_points.shape)>1:
200    # pressure_points = torch.squeeze(pressure_points,2)
201
202    pressure = pressure_points[:,:N]
203    pressure_fin_diff = pressure_points[:,N:]
204
205    split = torch.reshape(pressure_fin_diff,(B,2, ((2*D))*N // 2))
206    
207    grad = (split[:,0,:] - split[:,1,:]) / (2*stepsize)
208    
209    grad = torch.reshape(grad,(B,D,N))
210    grad_abs_square = torch.pow(torch.abs(grad),2)
211    grad_term = torch.sum(grad_abs_square,dim=1)
212
213    if K1 is None or K2 is None:
214        K1_, K2_ = get_gorkov_constants(V=V, c_0=medium_speed, c_p=particle_speed, p_0=medium_density, p_p=particle_density)
215        if K1 is None:
216            K1 = K1_
217        if K2 is None:
218            K2 = K2_
219    
220    # p_in =  torch.abs(pressure)
221    p_in = torch.sqrt(torch.real(pressure) **2 + torch.imag(pressure)**2)
222    if len(p_in.shape) > 2:
223        p_in.squeeze_(2)
224    # p_in = torch.squeeze(p_in,2)
225
226    # K1 = 1/4*V*(1/(c.c_0**2*c.p_0) - 1/(c.c_p**2*c.p_p))
227    # K2 = 3/4 * V * ((c.p_0 - c.p_p) / (c.angular_frequency**2 * c.p_0 * (c.p_0 * 2*c.p_p)))
228
229    U = K1 * p_in**2 - K2 *grad_term
230    
231    return U

Computes the Gorkov potential using finite differences to compute derivatives

Parameters
  • activation: The transducer activations to use
  • points: The points to compute the potential at
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis
  • stepsize: The distance aroud points to add, default 0.000135156253
  • K1: The value for K1 in the Gorkov equation, if None will use c.V / (4*c.p_0*c.c_0**2)
  • K2: The value for K2 in the Gorkov equation, if None will use 3*c.V / (4*(2*c.f**2 * c.p_0))
  • prop_function: Function to use to compute pressure
  • prop_fun_args: Arguments to pass to prop_function
  • board: The transducer boards to use if None use acoustools.Utilities.TRANSDUCERS
  • V: Particle Volume
Returns

gorkov potential at each point

from acoustools.Utilities import create_points, add_lev_sig
from acoustools.Solvers import wgs
from acoustools.Gorkov import gorkov_fin_diff

N=1
B=1
points = create_points(N,B)
x = wgs(points)
x = add_lev_sig(x)

U_fd = gorkov_fin_diff(x,points)

print("Finite Differences",U_fd.data.squeeze())
def get_finite_diff_points( points: torch.Tensor, axis: torch.Tensor, stepsize: float = 0.000135156253) -> torch.Tensor: 
234def get_finite_diff_points(points:Tensor , axis:Tensor, stepsize:float = 0.000135156253) -> Tensor:
235    '''
236    Gets points for finite difference calculations in one axis\n
237    :param points: Points around which to find surrounding points
238    :param axis: The axis to add points in
239    :param stepsize: The distance aroud points to add, default 0.000135156253
240    :return: points 
241    '''
242    #points = Bx3x4
243    points_h = points.clone()
244    points_neg_h = points.clone()
245    points_h[:,axis,:] = points_h[:,axis,:] + stepsize
246    points_neg_h[:,axis,:] = points_neg_h[:,axis,:] - stepsize
247
248    return points_h, points_neg_h

Gets points for finite difference calculations in one axis

Parameters
  • points: Points around which to find surrounding points
  • axis: The axis to add points in
  • stepsize: The distance aroud points to add, default 0.000135156253
Returns

points

def get_finite_diff_points_all_axis( points: torch.Tensor, axis: str = 'XYZ', stepsize: float = 0.000135156253) -> torch.Tensor: 
250def get_finite_diff_points_all_axis(points: Tensor,axis: str="XYZ", stepsize:float = 0.000135156253) -> Tensor:
251    '''
252    Gets points for finite difference calculations\\
253    :param points: Points around which to find surrounding points\\
254    :param axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\\
255    :param stepsize: The distance aroud points to add, default 0.000135156253\\
256    :return: Points
257    '''
258    B = points.shape[0]
259    D = len(axis)
260    N = points.shape[2]
261    fin_diff_points=  torch.zeros((B,3,((2*D)+1)*N)).to(device).to(DTYPE)
262    fin_diff_points[:,:,:N] = points.clone()
263
264    i = 2
265    if "X" in axis:
266        points_h, points_neg_h = get_finite_diff_points(points, 0, stepsize)
267        fin_diff_points[:,:,N:i*N] = points_h
268        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
269
270        i += 1
271
272    
273    if "Y" in axis:
274        points_h, points_neg_h = get_finite_diff_points(points, 1, stepsize)
275        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
276        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
277        i += 1
278    
279    if "Z" in axis:
280        points_h, points_neg_h = get_finite_diff_points(points, 2, stepsize)
281        fin_diff_points[:,:,(i-1)*N:i*N] = points_h
282        fin_diff_points[:,:,D*N+(i-1)*N:D*N+i*N] = points_neg_h
283        i += 1
284    
285    return fin_diff_points

Gets points for finite difference calculations\

Parameters
  • points: Points around which to find surrounding points\
  • axis: The axes to add points in as a string containing 'X', 'Y' and/or 'Z' eg 'XYZ' will use all three axis but 'YZ' will only add points in the YZ axis\
  • stepsize: The distance aroud points to add, default 0.000135156253\
Returns

Points