Extending diffraction functions to correctly calculate some rare scenarios

applications/SoundFieldSimulation/Diffraction/@itaInfiniteWedge/get_aperture_point_far_field.m

@@ -21,14 +21,31 @@ end
 % Auxilary plane spanned by source_receiver_dir and aux_plane_dir (closest
 % line between aperture vector and source-receiver-vector
 aux_plane_vec = cross( source_receiver_dir, aperture_dir );
-aux_plane_dir = aux_plane_vec / norm( aux_plane_vec );
-aux_plane_normal = cross( source_receiver_dir, aux_plane_dir );
-
-% Determine intersection of line (aperture) and auxiliary plane
-lambda_divisor = dot( aux_plane_normal, aperture_dir );
-assert( lambda_divisor ~= 0 )
-d = dot( aux_plane_normal, obj.location ); % Distance to origin
-lambda = ( d - dot( aux_plane_normal, obj.location ) ) / lambda_divisor;
-aperture_point = obj.location + lambda * aperture_dir;
+
+if norm( aux_plane_vec ) ~= 0
+    
+    % Directions are no parallel, a closest point must exist
+    aux_plane_dir = aux_plane_vec / norm( aux_plane_vec );
+    aux_plane_normal = cross( source_receiver_dir, aux_plane_dir );
+
+    % Determine intersection of line (aperture) and auxiliary plane
+    lambda_divisor = dot( aux_plane_normal, aperture_dir );
+    assert( lambda_divisor ~= 0 )
+    d = dot( aux_plane_normal, obj.location ); % Distance to origin
+    lambda = ( d - dot( aux_plane_normal, obj.location ) ) / lambda_divisor;
+    aperture_point = obj.location + lambda * aperture_dir;
+    
+else
+    
+    % Directions are parallel, project middle point between souce & target
+    % onto aperture
+    
+    % Project middle point onto aperture
+    lambda = dot( ( source_pos + source_receiver_vec ./ 2 )' - obj.location, aperture_dir );
+    aperture_point = obj.location + lambda * aperture_dir;
+        
+end
+
+assert( any( ~isnan( aperture_point ) ) )
 
 end
\ No newline at end of file

applications/SoundFieldSimulation/PropagationModels/@itaGeoPropagation/tf_diffraction.m

@@ -41,6 +41,17 @@ if size( effective_receiver_position, 1 ) == 4
     effective_receiver_position = effective_receiver_position( 1:3 )';
 end
 
+% Validate
+apex_point = w.get_aperture_point( effective_source_position( 1:3 )', effective_receiver_position( 1:3 )' );
+if all( apex_point == effective_source_position  ) || ...
+   all ( apex_point == effective_receiver_position )
+    warning( 'Skipping a path segment with a double anchor point' )
+    return
+elseif any( isnan( apex_point ) )
+    warning( 'Skipping a path segment with an invalid aperture point' )
+    return
+end
+
 % Return only the diffraction component and the subsequent
 % sphere-cylinder-shaped wave spreading loss factor to next effective
 % receiver (ignore effective source pos)
@@ -50,7 +61,6 @@ switch( diffraction_model )
         % Includes spreading loss after diffraction, but not phase shift
         linear_freq_data = obj.tf_diffraction_utd( w, effective_source_position, effective_receiver_position );
                 
-        apex_point = w.get_aperture_point( effective_source_position( 1:3 )', effective_receiver_position( 1:3 )' );
         eff_receiver_distance = norm( apex_point - effective_receiver_position );
         phase_delay_after_diffr = obj.phase_delay( eff_receiver_distance );
         
@@ -66,7 +76,6 @@ switch( diffraction_model )
         diffraction_dft = fft( btms_ir, obj.num_bins * 2 - 1 ); % odd DFT spectrum
         diffraction_hdft = diffraction_dft( 1:ceil( obj.num_bins ) );
         
-        apex_point = w.get_aperture_point( effective_source_position( 1:3 )', effective_receiver_position( 1:3 )' );
         eff_source_distance = norm( apex_point - effective_source_position );
         spreading_loss = ita_propagation_spreading_loss( eff_source_distance );
         phase_delay_after_diffr = obj.phase_delay( eff_source_distance );