src.acoustools.BEM.Gradients.E_Gradient

  1import torch
  2from torch import Tensor
  3
  4from vedo import Mesh
  5
  6import hashlib, pickle
  7
  8from acoustools.Utilities import DTYPE, forward_model_grad, forward_model_second_derivative_unmixed, forward_model_second_derivative_mixed
  9from acoustools.BEM.Forward_models import get_cache_or_compute_H
 10from acoustools.Mesh import get_areas, get_centres_as_points, get_normals_as_points
 11import acoustools.Constants as Constants
 12
 13
 14 
 15# def get_G_partial(points:Tensor, scatterer:Mesh, board:None=None, return_components:bool=False) -> tuple[Tensor, Tensor, Tensor]:
 16#     '''
 17#     Computes gradient of the G matrix in BEM \n
 18#     :param points: Points to propagate to
 19#     :param scatterer: The mesh used (as a `vedo` `mesh` object)
 20#     :param board: Ignored
 21#     :param return_components: if true will return the subparts used to compute
 22#     :return: Gradient of the G matrix in BEM
 23#     '''
 24#     #Bk3. Pg. 26
 25
 26#     areas = get_areas(scatterer)
 27#     centres = get_centres_as_points(scatterer)
 28#     normals = get_normals_as_points(scatterer)
 29
 30#     N = points.shape[2]
 31#     M = centres.shape[2]
 32
 33#     points = points.unsqueeze(3).expand(-1,-1,-1,M)
 34#     centres = centres.unsqueeze(2).expand(-1,-1,N,-1)
 35
 36#     diff = centres - points    
 37#     diff_square = diff**2
 38#     diff_square_sum = torch.sum(diff_square, 1)
 39#     distances = torch.sqrt(torch.sum(diff_square, 1))
 40#     distances_expanded = distances.unsqueeze(1).expand((1,3,N,M))
 41#     distances_expanded_square = distances_expanded**2
 42#     distances_expanded_cube = distances_expanded**3
 43
 44#     # G  =  e^(ikd) / 4pi d
 45#     G = areas * torch.exp(1j * Constants.k * distances_expanded) / (4*3.1415*distances_expanded)
 46
 47#     #Ga =  [i*da * e^{ikd} * (kd+i) / 4pi d^2]
 48
 49#     #d = distance
 50#     #da = -(at - a)^2 / d
 51
 52#     da = -diff / distances_expanded
 53#     kd = Constants.k * distances_expanded
 54#     phase = torch.exp(1j*kd)
 55#     # Ga =  areas * ( (1j*diff*phase * (kd - 1j))/ (4*3.1415*distances_expanded_cube))
 56#     print(diff)
 57#     Ga = (areas * diff * phase * (1-1j*Constants.k*distances_expanded)) / (4*Constants.pi * distances_expanded_cube)
 58
 59#     #P = (ik - 1/d)
 60#     P = (1j*Constants.k - 1/distances_expanded)
 61#     #Pa = da / d^2
 62#     Pa = -da / distances_expanded_cube
 63
 64#     #C = (diff \cdot normals) / distances
 65
 66#     nx = normals[:,0]
 67#     ny = normals[:,1]
 68#     nz = normals[:,2]
 69
 70#     dx = diff[:,0,:]
 71#     dy = diff[:,1,:]
 72#     dz = diff[:,2,:]
 73
 74#     C = (nx*dx + ny*dy + nz*dz) / distances
 75
 76
 77#     distances_cubed = distances**3
 78#     distance_square = distances**2
 79
 80#     # Cx = (nx*(dy**2 + dz**2) - dx * (ny*dy + nz*dz)) / distances_cubed
 81#     # Cy = (ny*(dx**2 + dz**2) - dy * (nx*dx + nz*dz)) / distances_cubed
 82#     # Cz = (nz*(dx**2 + dy**2) - dz * (nx*dx + ny*dy)) / distances_cubed
 83
 84#     Cx = (dx * (nx*dx + ny*dy + nz*dz) - nx*(diff_square_sum) ) / distances_cubed
 85#     Cy = (dy * (nx*dx + ny*dy + nz*dz) - ny*(diff_square_sum) ) / distances_cubed
 86#     Cz = (dz * (nx*dx + ny*dy + nz*dz) - nz*(diff_square_sum) ) / distances_cubed
 87
 88#     Cx.unsqueeze_(1)
 89#     Cy.unsqueeze_(1)
 90#     Cz.unsqueeze_(1)
 91
 92#     Ca = torch.cat([Cx, Cy, Cz],axis=1)
 93
 94#     grad_G = Ga*P*C + G*P*Ca + G*Pa*C
 95
 96#     grad_G =  grad_G.to(DTYPE)
 97
 98#     if return_components:
 99#         return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:], G,P,C,Ga,Pa,Ca
100    
101#     return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:]
102
103
104def BEM_forward_model_grad(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, 
105                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
106                           path:str="Media") -> tuple[Tensor, Tensor, Tensor] | tuple[Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor, Tensor]:
107    '''
108    Computes the gradient of the forward propagation for BEM\n
109    :param scatterer: The mesh used (as a `vedo` `mesh` object)
110    :param transducers: Transducers to use, if `None` uses `acoustools.Utilities.TRANSDUCERS`
111    :param use_cache_H_grad: If true uses the cache system, otherwise computes `H` and does not save it
112    :param print_lines: if true prints messages detaling progress
113    :param H: Precomputed `H` - if `None` `H` will be computed
114    :param return_components: if true will return the subparts used to compute
115    :param path: path to folder containing `BEMCache/` 
116    :return: Ex, Ey, Ez
117    '''
118    if transducers is None:
119        transducers = TRANSDUCERS
120
121    B = points.shape[0]
122    if H is None:
123        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines)
124    
125    Fx, Fy, Fz  = forward_model_grad(points, transducers)
126    Gx, Gy, Gz = get_G_partial(points, scatterer, transducers)
127
128
129    Ex = Fx - Gx@H #Minus coming out of the PM grad -> Not entoerly sure where that comes from but seems to fix it?
130    Ey = Fy - Gy@H
131    Ez = Fz - Gz@H
132
133
134    if return_components:
135        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE), Fx, Fy, Fz, Gx, Gy, Gz, H
136    else:
137        return Ex.to(DTYPE), Ey.to(DTYPE), Ez.to(DTYPE)
138    
139
140def get_G_second_unmixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False) -> tuple[Tensor, Tensor, Tensor]:
141
142    #Bk.3 Pg.80
143    #P is much higher than the others?
144
145    areas = get_areas(scatterer)
146    centres = get_centres_as_points(scatterer)
147    normals = get_normals_as_points(scatterer).unsqueeze(2)
148
149    N = points.shape[2]
150    M = centres.shape[2]
151
152    points.requires_grad_()
153    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True)
154    # exit()
155
156    points  = points.unsqueeze(3)  # [B, 3, N, 1]
157    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
158    
159
160    diff = points - centres    
161    diff_square = diff**2
162    diff_square_sum = torch.sum(diff_square, 1)
163    distances = torch.sqrt(diff_square_sum)
164    distances_square = distances ** 2
165    distances_cube = distances ** 3
166    distances_five = distances ** 5
167
168    distances_expanded = distances.unsqueeze(1).expand((1,3,N,M))
169    distances_expanded_square = distances_expanded**2
170    distances_expanded_cube = distances_expanded**3
171    distances_expanded_four = distances_expanded**4
172    distances_expanded_five = distances_expanded**5
173
174    distance_xy = diff[:,0,:,:] **2 + diff[:,1,:,:] **2
175    distance_xz = diff[:,0,:,:] **2 + diff[:,2,:,:] **2
176    distance_yz = diff[:,1,:,:] **2 + diff[:,2,:,:] **2
177
178
179    da = diff / distances_expanded
180    distance_bs = torch.stack([distance_yz,distance_xz,distance_xy], dim =1)
181    daa = distance_bs / distances_expanded_cube
182
183    kd = Constants.k * distances_expanded
184
185    exp_ikd = torch.exp(1j*kd)
186    
187    # G  =  e^(ikd) / 4pi d
188    # Gaa = (areas/4*Constants.pi) * (-1/distances_expanded_cube * torch.exp(1j*kd) * (da**2 * (kd**2 + 2*1j*kd - 2) + distances_expanded *daa * (1-1j * kd)))
189    # Gaa = -areas / (4*Constants.pi * distances_expanded_cube) * exp_ikd * (da**2 * (kd**2 + 2j * kd - 2) + distances_expanded *daa * (1-1j * kd))
190  
191    # Gaa = areas * ((-Constants.k**2 * diff_square * exp_ikd) / diff_square_sum + (1j*Constants.k * exp_ikd) / distances - (1j * Constants.k * diff_square * exp_ikd)/ distances_expanded_cube) / (4 * Constants.pi * distances)
192    # Gaa += areas * (((3*diff_square)/distances_expanded_five - 1/(distances_expanded_cube)) * exp_ikd) / (4 * Constants.pi)
193    # Gaa -= (1j * Constants.k * areas * diff_square * exp_ikd) / (2 * Constants.pi * distances_expanded_four) #Same as Gaa above
194    Gaa = areas/(4*Constants.pi * 1j) *exp_ikd * ((1j*Constants.k**2 * diff_square)/ distances_expanded_cube + (1j*Constants.k)/distances_expanded_square - 1/distances_expanded_cube - (3*1j*Constants.k*diff_square)/distances_expanded_four + (3*diff_square)/distances_expanded_five)
195    Gxx = Gaa[:,0]
196    Gyy = Gaa[:,1]
197    Gzz = Gaa[:,2]
198    
199    #P = ik - 1/d
200    # Paa = (distances_expanded * daa - 2*da**2) / distances_expanded_cube
201    # Paa = (3 * diff_square - distances_expanded_square) / distances_expanded_five
202    Paa = 1/distances_expanded_cube - (3*diff_square) / distances_expanded_five  #Equal to lines above
203    Pxx = Paa[:,0]
204    Pyy = Paa[:,1]
205    Pzz = Paa[:,2]
206
207    #C = (diff dot c )/ d
208
209    dx = diff[:,0,:,:]
210    dy = diff[:,1,:,:]
211    dz = diff[:,2,:,:]
212
213
214    nx = normals[:,0,:]
215    ny = normals[:,1,:]
216    nz = normals[:,2,:]
217
218    nd = nx * dx + ny * dy + nz * dz
219
220    Cxx = (2 * nx * dx)/distances_cube - nd/distances_cube + (3*nd*dx**2)/distances_five
221    Cyy = (2 * ny * dy)/distances_cube - nd/distances_cube + (3*nd*dy**2)/distances_five
222    Czz = (2 * nz * dz)/distances_cube - nd/distances_cube + (3*nd*dz**2)/distances_five
223
224    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1)
225
226    # Cxx = (((3 * dax**2)/distances_five) - (1/distances_cube)) * nd - (2*nx*dx) / distances_cube
227    # Cyy = (((3 * day**2)/distances_five) - (1/distances_cube)) * nd - (2*ny*dy) / distances_cube 
228    # Czz = (((3 * daz**2)/distances_five) - (1/distances_cube)) * nd - (2*nz*dz) / distances_cube 
229    # Caa = torch.stack([Cxx, Cyy, Czz], dim=1) 
230
231    # grad_2_G_unmixed_old = 2 * Ca * (Ga * P + G * Pa) + C * (2 * Ga*Pa + Gaa * P + G * Paa) + G*P*Caa
232    # grad_2_G_unmixed = C * (2*Ga * Pa + Gaa * P + G*Paa) + P * (2*Ga * Ca + G*Caa) + 2*G * Pa * Ca
233
234    Cx = Ca[:,0]
235    Cy = Ca[:,1]
236    Cz = Ca[:,2]
237
238    Px = Pa[:,0]
239    Py = Pa[:,1]
240    Pz = Pa[:,2]
241
242    Gx = Ga[:,0]
243    Gy = Ga[:,1]
244    Gz = Ga[:,2]
245
246    G = G[:,0,:]
247    P = P[:,0,:] 
248    
249
250
251    grad_2_G_unmixed_xx = Gxx * P * C + G*Pxx*C + G*P*Cxx+  2*Gx*Px*C + 2*Gx*P*Cx + 2*G*Px*Cx
252    grad_2_G_unmixed_yy = Gyy * P * C + G*Pyy*C + G*P*Cyy+  2*Gy*Py*C + 2*Gy*P*Cy + 2*G*Py*Cy
253    grad_2_G_unmixed_zz = Gzz * P * C + G*Pzz*C + G*P*Czz+  2*Gz*Pz*C + 2*Gz*P*Cz + 2*G*Pz*Cz
254
255    if return_components:
256        return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz, G,P,C,Ga,Pa,Ca
257    return grad_2_G_unmixed_xx, grad_2_G_unmixed_yy, grad_2_G_unmixed_zz
258
259
260def get_G_second_mixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False) -> tuple[Tensor, Tensor, Tensor]:
261
262    #Bk.3 Pg.81
263    
264    areas = get_areas(scatterer).unsqueeze(1)
265    centres = get_centres_as_points(scatterer)
266    normals = get_normals_as_points(scatterer).unsqueeze(2)
267
268    N = points.shape[2]
269    M = centres.shape[2]
270
271    grad_Gx, grad_Gy, grad_Gz, G,P,C,Ga,Pa,Ca = get_G_partial(points, scatterer, board, True)
272
273    points  = points.unsqueeze(3)  # [B, 3, N, 1]
274    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
275
276    diff = points - centres    
277    diff_square = diff**2
278    distances = torch.sqrt(torch.sum(diff_square, 1))
279    distances_square = distances ** 2
280    distances_cube = distances ** 3
281    distances_four = distances ** 4
282    distances_five = distances ** 5
283    distances_six = distances ** 6
284    
285    kd = Constants.k * distances
286    ikd = 1j * kd
287
288    dx = diff[:,0,:,:]
289    dy = diff[:,1,:,:]
290    dz = diff[:,2,:,:]
291
292    # print(dx.shape, dy.shape, dz.shape, distances_cube.shape)
293    dxy = -dx*dy / distances_cube
294    dxz = -dx*dz / distances_cube
295    dyz = -dy*dz / distances_cube
296
297    # print(dx.shape, distances.shape)
298    dax = dx / distances
299    day = dy / distances
300    daz = dz / distances
301
302    exp_ikd = torch.exp(ikd)
303
304
305
306    # print(areas.shape, exp_ikd.shape, day.shape, dax.shape, distances.shape, distances_cube.shape)
307    Gxy = (-areas * exp_ikd * (day * dax * (kd**2 + 2*ikd - 2) + distances * dxy * (1 - ikd))) / (4*Constants.pi * distances_cube)
308    Gxz = (-areas * exp_ikd * (daz * dax * (kd**2 + 2*ikd - 2) + distances * dxz * (1 - ikd))) / (4*Constants.pi * distances_cube)
309    Gyz = (-areas * exp_ikd * (day * daz * (kd**2 + 2*ikd - 2) + distances * dyz * (1 - ikd))) / (4*Constants.pi * distances_cube)
310
311    # Gab = (areas/(Constants.pi*4j)) * exp_ikd * ((-1*Constants.k)/distances_cube - (3j * Constants.k)/distances_four + 3/distances_five )
312   
313    # print(exp_ikd.shape, distances_six.shape, distances.shape, dx.shape,dy.shape,dz.shape)
314    # Gab = (areas/(Constants.pi*4j)) * exp_ikd * (1/distances_six * (3-1j*Constants.k*distances))
315    # Gxy = (dx * dy) * Gab
316    # Gxz = (dx * dz) * Gab
317    # Gyz = (dz * dy) * Gab
318
319
320    nx = normals[:,0,:]
321    ny = normals[:,1,:]
322    nz = normals[:,2,:]
323
324    # print(nx.shape, dx.shape, ny.shape, dy.shape, nz.shape, dz.shape)
325    nd = nx * dx + ny * dy + nz * dz
326
327    # Cxy = (-ny*dx*distances_square - nx*dy*distances_square + 3 * dy*dx * nd) / distances_five
328    # Cxz = (-nz*dx*distances_square - nx*dz*distances_square + 3 * dz*dx * nd) / distances_five
329    # Cyz = (-ny*dz*distances_square - ny*dz*distances_square + 3 * dz*dy * nd) / distances_five
330
331    # print(nx.shape, dy.shape, ny.shape, dx.shape, distances_cube.shape,distances_square.shape )
332    Cxy = (-nx*dy - ny*dx + (3*nd*dx*dy/distances_cube)) / distances_square
333    Cxz = (-nx*dz - nz*dx + (3*nd*dx*dz/distances_cube)) / distances_square
334    Cyz = (-nz*dy - ny*dz + (3*nd*dz*dy/distances_cube)) / distances_square
335
336    # Cxy = (2*ny*dx - nx*dy) / distances_cube - (3*(ny*distances_square - nd*dy)* dx) / distances_five
337    # Cxz = (2*nz*dx - nx*dz) / distances_cube - (3*(nz*distances_square - nd*dz)* dx) / distances_five
338    # Cyz = (2*ny*dz - nz*dy) / distances_cube - (3*(ny*distances_square - nd*dy)* dz) / distances_five
339    
340   
341    # Pxy = -3 * (dx * dy) / distances_five
342    # Pxz = -3 * (dx * dz) / distances_five
343    # Pyz = -3 * (dz * dy) / distances_five
344
345    Pxy = (distances * dxy - 2 * day * dax) / distances_cube
346    Pxz = (distances * dxz - 2 * daz * dax) / distances_cube
347    Pyz = (distances * dyz - 2 * day * daz) / distances_cube
348
349    
350    G = G[:,0,:]
351    P = P[:,0,:] 
352    # C = C.unsqueeze(1).expand(-1,3,-1,-1)
353
354    Gx = Ga[:,0,:]
355    Gy = Ga[:,1,:]
356    Gz = Ga[:,2,:]
357
358    Px = Pa[:,0,:]
359    Py = Pa[:,1,:]
360    Pz = Pa[:,2,:]
361
362    Cx = Ca[:,0,:]
363    Cy = Ca[:,1,:]
364    Cz = Ca[:,2,:]
365
366
367    # print(Gx*Py*C, G*Py*Cx)
368    # print(Gxy*P*C , Gy*Px*C , Gy*P*Cx , Gx*Py*C , G*Pxy*C , G*Py*Cx , Gx*P*Cy , G*Px*Cy , G*P*Cxy, sep="\n\n")
369    # print(G.shape, Gx.shape, Gy.shape, Gz.shape,P.shape , Px.shape, Py.shape, Pz.shape, C.shape, Cx.shape, Cy.shape, Cz.shape )
370
371    # grad_2_G_mixed_xy = Gxy*P*C + Gy*Px*C + Gy*P*Cx + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*Px*Cy + G*P*Cxy
372    # grad_2_G_mixed_xz = Gxz*P*C + Gz*Px*C + Gz*P*Cx + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*Px*Cz + G*P*Cxz
373    # grad_2_G_mixed_yz = Gyz*P*C + Gy*Pz*C + Gy*P*Cz + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*Pz*Cy + G*P*Cyz
374
375    grad_2_G_mixed_xy = Gxy*P*C + Gx*Py*C + G*Pxy*C + G*Py*Cx + Gx*P*Cy + G*P*Cxy
376    grad_2_G_mixed_xz = Gxz*P*C + Gx*Pz*C + G*Pxz*C + G*Pz*Cx + Gx*P*Cz + G*P*Cxz
377    grad_2_G_mixed_yz = Gyz*P*C + Gz*Py*C + G*Pyz*C + G*Py*Cz + Gz*P*Cy + G*P*Cyz
378
379
380
381    return grad_2_G_mixed_xy, grad_2_G_mixed_xz, grad_2_G_mixed_yz
382
383
384def BEM_forward_model_second_derivative_unmixed(points:Tensor, scatterer:Mesh, transducers:Tensor=None, use_cache_H:bool=True, 
385                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
386                           path:str="Media"):
387    
388    
389    if transducers is None:
390        transducers = TRANSDUCERS
391
392    if H is None:
393        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines)
394
395    Fxx, Fyy, Fzz = forward_model_second_derivative_unmixed(points,transducers=transducers)
396    Gxx, Gyy, Gzz = get_G_second_unmixed(points, scatterer, transducers)
397
398    Exx = Fxx + Gxx@H
399    Eyy = Fyy + Gyy@H
400    Ezz = Fzz + Gzz@H
401
402    return Exx, Eyy, Ezz
403
404def BEM_forward_model_second_derivative_mixed(points:Tensor, scatterer:Mesh, transducers:Tensor|Mesh=None, use_cache_H:bool=True, 
405                           print_lines:bool=False, H:Tensor|None=None, return_components:bool=False,
406                           path:str="Media"):
407    
408       
409    if transducers is None:
410        transducers = TRANSDUCERS
411
412    if H is None:
413        H = get_cache_or_compute_H(scatterer,transducers,use_cache_H, path, print_lines)
414
415    Fxy, Fxz, Fyz = forward_model_second_derivative_mixed(points,transducers=transducers)
416    Gxy, Gxz, Gyz = get_G_second_mixed(points, scatterer, transducers)
417
418
419    Exy = Fxy + Gxy@H
420    Exz = Fxz + Gxz@H
421    Eyz = Fyz + Gyz@H
422
423    return Exy, Exz, Eyz
424
425
426
427
428
429# def get_G_second_mixed(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False) -> tuple[Tensor, Tensor, Tensor]:
430
431#     #Bk.3 Pg.81
432    
433#     areas = get_areas(scatterer)
434#     centres = get_centres_as_points(scatterer)
435#     normals = get_normals_as_points(scatterer).unsqueeze(2)
436
437#     N = points.shape[2]
438#     M = centres.shape[2]
439
440#     points = points.unsqueeze(3).expand(-1,-1,-1,M)
441#     centres = centres.unsqueeze(2).expand(-1,-1,N,-1)
442
443#     diff = points - centres    
444#     diff_square = diff**2
445#     distances = torch.sqrt(torch.sum(diff_square, 1))
446#     distances_square = distances ** 2
447#     distances_cube = distances ** 3
448#     distances_four = distances ** 4
449#     distances_five = distances ** 5
450
451    
452#     distances_expanded = distances.unsqueeze(1).expand((1,3,N,M)) 
453
454#     kd = Constants.k * distances
455#     exp_ikd = torch.exp(1j*kd)
456
457
458#     dx = diff[:,0,:,:]
459#     dy = diff[:,1,:,:]
460#     dz = diff[:,2,:,:]
461
462#     dxy = -dx*dy / distances_cube
463#     dxz = -dx*dz / distances_cube
464#     dyz = -dy*dz / distances_cube
465
466#     da = diff / distances_expanded
467#     dax = da[:,0,:,:]
468#     day = da[:,1,:,:]
469#     daz = da[:,2,:,:]
470 
471#     F = (areas * 1j * Constants.k * exp_ikd) / ( 4 * Constants.pi * distances)
472#     P = (-1 * areas * exp_ikd) / ( 4 * Constants.pi * distances_square)
473
474#     Fa = (-1 * Constants.k * areas * da * exp_ikd * (kd + 1j)) / (Constants.pi * 4 * distances_square)
475#     Pa = (areas * da * exp_ikd * (2 - 1j*kd)) / (Constants.pi * 4 * distances_cube)
476
477
478#     Fxy = (Constants.k * areas * exp_ikd)/(4 * Constants.pi * distances_cube) * (day * dax * (-1j * kd**2 + 2*kd + 2j) - distances*dxy * (kd + 1j))
479#     Fxz = (Constants.k * areas * exp_ikd)/(4 * Constants.pi * distances_cube) * (daz * dax * (-1j * kd**2 + 2*kd + 2j) - distances*dxz * (kd + 1j))
480#     Fyz = (Constants.k * areas * exp_ikd)/(4 * Constants.pi * distances_cube) * (day * daz * (-1j * kd**2 + 2*kd + 2j) - distances*dyz * (kd + 1j))
481
482#     Pxy = (areas * exp_ikd)/(4 * Constants.pi * distances_four) * (day * dax * (kd**2 + 4j*kd -6) -distances*dxy * (2-1j*kd))
483#     Pxz = (areas * exp_ikd)/(4 * Constants.pi * distances_four) * (daz * dax * (kd**2 + 4j*kd -6) -distances*dxz * (2-1j*kd))
484#     Pyz = (areas * exp_ikd)/(4 * Constants.pi * distances_four) * (day * daz * (kd**2 + 4j*kd -6) -distances*dyz * (2-1j*kd))
485
486
487#     nx = normals[:,0,:]
488#     ny = normals[:,1,:]
489#     nz = normals[:,2,:]
490
491#     nd = nx * dx + ny * dy + nz * dz
492
493#     C = (diff * normals).sum(dim=1) / distances
494
495#     Cxy = (-ny*dx*distances_square - nx*dy*distances_square - 3 * dy*dx * nd) / distances_five
496#     Cxz = (-nz*dx*distances_square - nx*dz*distances_square - 3 * dz*dx * nd) / distances_five
497#     Cyz = (-ny*dz*distances_square - ny*dz*distances_square - 3 * dz*dy * nd) / distances_five
498
499
500#     Cx = (nx*(dx**2 +dy**2 + dz**2) - dx * (ny*dy + nz*dz)) / distances_cube
501#     Cy = (ny*(dx**2 +dy**2 + dz**2) - dy * (nx*dx + nz*dz)) / distances_cube
502#     Cz = (nz*(dx**2 +dy**2 + dz**2) - dz * (nx*dx + ny*dy)) / distances_cube
503
504#     Fx = Fa[:,0,:]
505#     Fy = Fa[:,1,:]
506#     Fz = Fa[:,2,:]
507
508#     Px = Pa[:,0,:]
509#     Py = Pa[:,1,:]
510#     Pz = Pa[:,2,:]
511
512#     print(Fx.shape, Fxy.shape)
513#     print(Px.shape, Pxy.shape )
514#     print(Cx.shape, Cxy.shape, C.shape)
515
516
517#     grad_2_G_mixed_xy = Cy * (Fx + Px) + Cx * (Fy + Py) + C * (Fxy + Pxy) + Cxy * (F+P)
518#     grad_2_G_mixed_xz = Cz * (Fx + Px) + Cx * (Fz + Pz) + C * (Fxz + Pxz) + Cxz * (F+P)
519#     grad_2_G_mixed_yz = Cy * (Fz + Pz) + Cz * (Fy + Py) + C * (Fyz + Pyz) + Cyz * (F+P)
520
521
522
523#     return grad_2_G_mixed_xy, grad_2_G_mixed_xz, grad_2_G_mixed_yz
524
525def get_G_partial(points:Tensor, scatterer:Mesh, board:Tensor|None=None, return_components:bool=False) -> tuple[Tensor, Tensor, Tensor]:
526    '''
527    Computes gradient of the G matrix in BEM \n
528    :param points: Points to propagate to
529    :param scatterer: The mesh used (as a `vedo` `mesh` object)
530    :param board: Ignored
531    :param return_components: if true will return the subparts used to compute
532    :return: Gradient of the G matrix in BEM
533    '''
534    #Bk3. Pg. 26
535    # if board is None:
536    #     board = TRANSDUCERS
537
538    areas = get_areas(scatterer)
539    centres = get_centres_as_points(scatterer)
540    normals = get_normals_as_points(scatterer)
541
542
543    N = points.shape[2]
544    M = centres.shape[2]
545
546
547    # points = points.unsqueeze(3).expand(-1,-1,-1,M)
548    # centres = centres.unsqueeze(2).expand(-1,-1,N,-1)
549    points  = points.unsqueeze(3)  # [B, 3, N, 1]
550    centres = centres.unsqueeze(2)  # [B, 3, 1, M]
551
552    diff = points - centres    
553    diff_square = diff**2
554    distances = torch.sqrt(torch.sum(diff_square, 1))
555    distances_expanded = distances.unsqueeze(1)#.expand((1,3,N,M))
556    distances_expanded_square = distances_expanded**2
557    distances_expanded_cube = distances_expanded**3
558
559    # G  =  e^(ikd) / 4pi d
560    G = areas * torch.exp(1j * Constants.k * distances_expanded) / (4*3.1415*distances_expanded)
561
562    #Ga =  [i*da * e^{ikd} * (kd+i) / 4pi d^2]
563
564    #d = distance
565    #da = -(at - a)^2 / d
566    da = diff / distances_expanded
567    kd = Constants.k * distances_expanded
568    phase = torch.exp(1j*kd)
569    Ga =  areas * ( (1j*da*phase * (kd + 1j))/ (4*3.1415*distances_expanded_square))
570
571    #P = (ik - 1/d)
572    P = (1j*Constants.k - 1/distances_expanded)
573    #Pa = da / d^2
574    # Pa = da / distances_expanded_square
575    Pa = diff / distances_expanded_cube
576
577    #C = (diff \cdot normals) / distances
578
579    nx = normals[:,0]
580    ny = normals[:,1]
581    nz = normals[:,2]
582
583    dx = diff[:,0,:]
584    dy = diff[:,1,:]
585    dz = diff[:,2,:]
586
587    n_dot_d = nx*dx + ny*dy + nz*dz
588
589    dax = da[:,0,:]
590    day = da[:,1,:]
591    daz = da[:,2,:]
592
593    C = (nx*dx + ny*dy + nz*dz) / distances
594
595
596    distances_cubed = distances**3
597    distance_square = distances**2
598
599
600    Cx = 1/distance_square * (nx * distances - (n_dot_d * dx) / distances)
601    Cy = 1/distance_square * (ny * distances - (n_dot_d * dy) / distances)
602    Cz = 1/distance_square * (nz * distances - (n_dot_d * dz) / distances)
603
604    Cx.unsqueeze_(1)
605    Cy.unsqueeze_(1)
606    Cz.unsqueeze_(1)
607
608    Ca = torch.cat([Cx, Cy, Cz],axis=1)
609
610    grad_G = Ga*P*C + G*P*Ca + G*Pa*C
611
612    grad_G =  grad_G.to(DTYPE)
613
614    if return_components:
615        return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:], G,P,C,Ga,Pa, Ca
616    
617    return grad_G[:,0,:], grad_G[:,1,:], grad_G[:,2,:]