/[MITgcm]/MITgcm_contrib/high_res_cube/matlab-grid-generator/excess_of_quad.m
ViewVC logotype

Annotation of /MITgcm_contrib/high_res_cube/matlab-grid-generator/excess_of_quad.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.1.1.1 - (hide annotations) (download) (vendor branch)
Tue Nov 11 18:08:08 2003 UTC (21 years, 8 months ago) by cnh
Branch: MAIN, initial
CVS Tags: baseline, HEAD
Changes since 1.1: +0 -0 lines 
Checking in work done with Dimitri on high-resolution cube gridding and parallel 
communications. 
   o code is in a contrib experiment for now so we can continue collaborating
     on it. However most code is general and will be moved into main branch once 
     it is fully hardened.
   o There are README files in the contrib root and in the subdirectories that
     explain the contents
1 cnh 1.1 function [excess]=excess_of_quad(V1,V2,V3,V4)
2     % excess=excess_of_quad(V1,V2,V3,V4);
3     %
4     % where Vn are vectors
5    
6     plane1=plane_normal(V1,V2);
7     plane2=plane_normal(V2,V3);
8     plane3=plane_normal(V3,V4);
9     plane4=plane_normal(V4,V1);
10    
11     angle12=pi-angle_between_vectors(plane2,plane1);
12     angle23=pi-angle_between_vectors(plane3,plane2);
13     angle34=pi-angle_between_vectors(plane4,plane3);
14     angle41=pi-angle_between_vectors(plane1,plane4);
15    
16     excess=angle12+angle23+angle34+angle41-2*pi;
17    
18     %%%
19    
20     function [plane]=plane_normal(P1,P2);
21    
22     if ndims(P1)==4
23     plane(:,:,:,1)=P1(:,:,:,2).*P2(:,:,:,3)-P1(:,:,:,3).*P2(:,:,:,2);
24     plane(:,:,:,2)=P1(:,:,:,3).*P2(:,:,:,1)-P1(:,:,:,1).*P2(:,:,:,3);
25     plane(:,:,:,3)=P1(:,:,:,1).*P2(:,:,:,2)-P1(:,:,:,2).*P2(:,:,:,1);
26     mag=sqrt(plane(:,:,:,1).^2 + plane(:,:,:,2).^2 + plane(:,:,:,3).^2);
27     plane(:,:,:,1)=plane(:,:,:,1)./mag;
28     plane(:,:,:,2)=plane(:,:,:,2)./mag;
29     plane(:,:,:,3)=plane(:,:,:,3)./mag;
30     elseif ndims(P1)==3
31     plane(:,:,1)=P1(:,:,2).*P2(:,:,3)-P1(:,:,3).*P2(:,:,2);
32     plane(:,:,2)=P1(:,:,3).*P2(:,:,1)-P1(:,:,1).*P2(:,:,3);
33     plane(:,:,3)=P1(:,:,1).*P2(:,:,2)-P1(:,:,2).*P2(:,:,1);
34     mag=sqrt(plane(:,:,1).^2 + plane(:,:,2).^2 + plane(:,:,3).^2);
35     plane(:,:,1)=plane(:,:,1)./mag;
36     plane(:,:,2)=plane(:,:,2)./mag;
37     plane(:,:,3)=plane(:,:,3)./mag;
38     else
39     error('Wrong ndims')
40     end
41    
  ViewVC Help
Powered by ViewVC 1.1.22