src.acoustools.Utilities.Forward_models

  1import torch
  2from torch import Tensor
  3
  4from acoustools.Utilities.Boards import TRANSDUCERS
  5from acoustools.Utilities.Utilities import is_batched_points
  6from acoustools.Utilities.Setup import device, DTYPE
  7import acoustools.Constants as Constants
  8
  9
 10def forward_model(points:Tensor, transducers:Tensor|None = None, p_ref = Constants.P_ref, norms:Tensor=None, k=Constants.k, transducer_radius = Constants.radius) -> Tensor:
 11    '''
 12    Create the piston model forward propagation matrix for points and transducers\\
 13    :param points: Point position to compute propagation to \\
 14    :param transducers: The Transducer array, default two 16x16 arrays \\
 15    :param p_ref: The value to use for p_ref\\
 16    :param norms: Tensor of normals for transduers\\ 
 17    Returns forward propagation matrix \\
 18    '''
 19    if transducers is None:
 20        transducers = TRANSDUCERS
 21    
 22    if norms is None:
 23        norms = (torch.zeros_like(transducers) + torch.tensor([0,0,1], device=device)) * torch.sign(transducers[:,2].real).unsqueeze(1).to(DTYPE)
 24
 25    if is_batched_points(points):
 26        return forward_model_batched(points, transducers, p_ref=p_ref,norms=norms, k=k, transducer_radius=transducer_radius).to(DTYPE)
 27    else:
 28        return forward_model_unbatched(points, transducers, p_ref=p_ref,norms=norms, k=k, transducer_radius=transducer_radius).to(DTYPE)
 29
 30def forward_model_unbatched(points, transducers = TRANSDUCERS, p_ref = Constants.P_ref,norms:Tensor|None=None, k=Constants.k, transducer_radius = Constants.radius):
 31    '''
 32    @private
 33    Compute the piston model for acoustic wave propagation NOTE: Unbatched, use `forward_model_batched` for batched computation \\
 34    `points` Point position to compute propagation to \\
 35    `transducers` The Transducer array, default two 16x16 arrays \\
 36    :param p_ref: The value to use for p_ref\\
 37    :param norms: Tensor of normals for transduers\\ 
 38    Returns forward propagation matrix \\
 39    Written by Giorgos Christopoulos, 2022
 40    '''
 41
 42    if norms is None:
 43        norms = (torch.zeros_like(transducers) + torch.tensor([0,0,1], device=device)) * torch.sign(transducers[:,2].real).unsqueeze(1).to(DTYPE)
 44    
 45    m=points.size()[1]
 46    n=transducers.size()[0]
 47    
 48    transducers_x=torch.reshape(transducers[:,0],(n,1))
 49    transducers_y=torch.reshape(transducers[:,1],(n,1))
 50    transducers_z=torch.reshape(transducers[:,2],(n,1))
 51
 52
 53    points_x=torch.reshape(points[0,:],(m,1))
 54    points_y=torch.reshape(points[1,:],(m,1))
 55    points_z=torch.reshape(points[2,:],(m,1))
 56    
 57    dx = (transducers_x.T-points_x) **2
 58    dy = (transducers_y.T-points_y) **2
 59    dz = (transducers_z.T-points_z) **2
 60
 61    distance=torch.sqrt(dx+dy+dz)
 62
 63    distance_axis_sub = torch.sub(points,transducers).to(DTYPE)
 64    norms = norms.unsqueeze(3).expand( m, 3, n)
 65    sine_angle = torch.norm(torch.cross(distance_axis_sub, norms, dim=1),2, dim=1) / distance
 66
 67
 68    bessel_arg=k*transducer_radius*sine_angle #planar_dist / dist = sin(theta)
 69
 70    directivity=1/2-torch.pow(bessel_arg,2)/16+torch.pow(bessel_arg,4)/384
 71    
 72    p = 1j*k*distance
 73    phase = torch.exp(p)
 74    
 75    trans_matrix=2*p_ref*torch.multiply(torch.divide(phase,distance),directivity)
 76    return trans_matrix
 77
 78def forward_model_batched(points, transducers = TRANSDUCERS, p_ref = Constants.P_ref,norms:Tensor|None=None, k=Constants.k, transducer_radius=Constants.radius):
 79
 80    '''
 81    @private
 82    computed batched piston model for acoustic wave propagation
 83    `points` Point position to compute propagation to \\
 84    `transducers` The Transducer array, default two 16x16 arrays \\
 85    :param p_ref: The value to use for p_ref\\
 86    :param norms: Tensor of normals for transduers\\ 
 87    Returns forward propagation matrix \\
 88    '''
 89    B = points.shape[0]
 90    N = points.shape[2]
 91    M = transducers.shape[0]
 92
 93    if norms is None:
 94        norms = (torch.zeros_like(transducers) + torch.tensor([0,0,1], device=device)) * torch.sign(transducers[:,2].real).unsqueeze(1).to(DTYPE)
 95    
 96    # p = torch.permute(points,(0,2,1))
 97    transducers = torch.unsqueeze(transducers,2)
 98    transducers = transducers.expand((B,-1,-1,N))
 99    points = torch.unsqueeze(points,1)
100    points = points.expand((-1,M,-1,-1))
101
102    # distance_axis = (transducers - points) **2
103    # distance_axis_sub = transducers - points
104    distance_axis_sub = torch.sub(points,transducers).to(DTYPE)
105    distance_axis = distance_axis_sub * distance_axis_sub
106    distance = torch.sqrt(torch.sum(distance_axis,dim=2))
107    
108    # sine_angle = torch.divide(planar_distance,distance)
109    norms = norms.unsqueeze(0).unsqueeze(3).expand(B, M, 3, N)
110
111    sine_angle = torch.norm(torch.cross(distance_axis_sub, norms, dim=2),2, dim=2) / distance
112
113
114
115
116    bessel_arg=k*transducer_radius*sine_angle
117    directivity=1/2-torch.pow(bessel_arg,2)/16+torch.pow(bessel_arg,4)/384
118    
119    p = 1j*k*distance
120    phase = torch.exp(p)
121
122    trans_matrix=2*p_ref*torch.multiply(torch.divide(phase,distance),directivity)
123
124    return trans_matrix.permute((0,2,1)).to(DTYPE).to(device)
125
126def green_propagator(points:Tensor, board:Tensor, k:float=Constants.k) -> Tensor:
127    '''
128    Computes the Green's function propagation matrix from `board` to `points` \n
129    :param points: Points to use
130    :param board: transducers to use
131    :param k: Wavenumber of sound to use
132    :return: Green propagation matrix
133    '''
134
135    B = points.shape[0]
136    N = points.shape[2]
137    M = board.shape[0]
138    board = board.unsqueeze(0).unsqueeze_(3)
139    points = points.unsqueeze(1)
140    
141    # distances_axis = torch.abs(points-board)
142    distances_axis = (points-board)**2
143    distances = torch.sqrt(torch.sum(distances_axis, dim=2))
144
145
146    
147    # green = -1* (torch.exp(1j*k*distances)) / (4 * Constants.pi *distances)
148    green = -1* (torch.exp(1j*k*distances)) / (4 * Constants.pi *distances)
149
150
151    return green.mT