Adding script to analyse angle errors evaluating Snell's law

matlab/propagation_path_eval_Snells_law_angle_error.m

+%% Load
+
+json_file_path = 'CombinedModelTest_Test.skp.json';
+txt = fileread( json_file_path );
+json_content = jsondecode( txt );
+
+N = numel( json_content.propagation_paths );
+fprintf( 'Found %i propagation paths\n', N )
+
+
+%% Eval
+
+error_angle_refl_rad = [];
+error_angle_diffr_rad = [];
+
+for n = 1:N
+    prop_path = json_content.propagation_paths( n );
+    assert( strcmpi( prop_path.class, 'propagation_path' ) )
+    
+    for a = 2:( numel( prop_path.propagation_anchors ) - 1 )
+        
+        if iscell( prop_path.propagation_anchors( a - 1 ) )
+            prop_anchor_prev = prop_path.propagation_anchors{ a - 1 };
+        else
+            prop_anchor_prev = prop_path.propagation_anchors( a - 1 );
+        end
+        
+        if iscell( prop_path.propagation_anchors( a ) )
+            prop_anchor_curr = prop_path.propagation_anchors{ a };
+        else
+            prop_anchor_curr = prop_path.propagation_anchors( a );
+        end
+        
+        if iscell( prop_path.propagation_anchors( a + 1 ) )
+            prop_anchor_next = prop_path.propagation_anchors{ a + 1 };
+        else
+            prop_anchor_next = prop_path.propagation_anchors( a + 1 );
+        end
+        
+        prop_segment_incoming_vec = prop_anchor_curr.interaction_point - prop_anchor_prev.interaction_point;
+        prop_segment_outgoing_vec = prop_anchor_next.interaction_point - prop_anchor_curr.interaction_point;
+        
+        if strcmpi( prop_anchor_curr.anchor_type, 'specular_reflection' ) && isfield( prop_anchor_curr, 'face_normal' )
+            a1 = dot( prop_segment_incoming_vec, prop_anchor_curr.face_normal ) / norm( prop_segment_incoming_vec ) / norm( prop_anchor_curr.face_normal );
+            a2 = ( -1 ) * dot( prop_segment_outgoing_vec, prop_anchor_curr.face_normal ) / norm( prop_segment_outgoing_vec ) / norm( prop_anchor_curr.face_normal );
+            error_angle_refl_rad( end + 1 ) = abs( acos( a2 ) - acos( a1 ) );
+        elseif strcmpi( prop_anchor_curr.anchor_type, 'outer_edge_diffraction' ) || strcmpi( prop_anchor_curr.anchor_type, 'inner_edge_diffraction' )
+            % @todo get edge normal, validate angles
+            vertex_dir = prop_anchor_curr.vertex_end -  prop_anchor_curr.vertex_start;
+            a1 = dot( prop_segment_incoming_vec, vertex_dir ) / norm( prop_segment_incoming_vec ) / norm( vertex_dir );
+            a2 = dot( prop_segment_outgoing_vec , vertex_dir ) / norm( prop_segment_outgoing_vec ) / norm( vertex_dir );
+            error_angle_diffr_rad( end + 1 ) = abs( acos( a2 ) - acos( a1 ) );
+        end
+    end
+end
+
+
+%% Statistics
+
+error_angle_rad = [ error_angle_refl_rad error_angle_diffr_rad ];
+fprintf( 'Propagation path angular error based on Snells law: %.1frad (%.1frad)\n', mean( error_angle_rad ), std2( error_angle_rad ) )
+fprintf( '\tspecular reflections only: %.3frad (%.3frad)\n', mean( error_angle_refl_rad ), std2( error_angle_refl_rad ) )
+fprintf( '\tspecular diffractions only: %.3frad (%.3frad)\n', mean( error_angle_diffr_rad ), std2( error_angle_diffr_rad ) )
+
+plot( error_angle_diffr_rad )
\ No newline at end of file