/[MITgcm]/MITgcm_contrib/cnh_cs_plot/make_cs_segments_and_patches.m
ViewVC logotype

Annotation of /MITgcm_contrib/cnh_cs_plot/make_cs_segments_and_patches.m

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


Revision 1.1 - (hide annotations) (download)
Tue Mar 15 18:05:55 2005 UTC (19 years, 2 months ago) by cnh
Branch: MAIN
CVS Tags: HEAD
Adding cube face plotting utility + t32 example that it uses.

1 cnh 1.1 function csg=make_cs_segments_and_patches(csg)
2     %
3     % Add line segments and polygon patches to CS grid object.
4     %
5    
6     clf
7     csg
8     xg=csg.gridarr(:,:,csg.xgpos);
9     yg=csg.gridarr(:,:,csg.ygpos);
10     iind=csg.index_i_center;
11     jind=csg.index_j_center;
12     clear xc yc;
13     nc=0;
14     np=0;
15     ng=csg.ng;
16     xcoord=cast(ones(2,ng*ng*6*4),class(xg));
17     ycoord=cast(ones(2,ng*ng*6*4),class(xg));
18    
19     fprintf('Beginning segment calculation\n');
20    
21     for t=1:numel(csg.fx)
22     ilo=csg.ilog(t);
23     jlo=csg.jlog(t);
24     ihi=csg.ilog(t)+csg.fx(t);
25     jhi=csg.jlog(t)+csg.fy(t);
26     nlo=nc+1;
27     nhi=nc+(ihi-ilo)*(jhi-jlo)*4;
28    
29     % Set i ,j -> i+1,j segment
30     r1=[ilo :ihi-1];r2=[jlo:jhi-1];
31     r3=[ilo+1:ihi ];r4=[jlo:jhi-1];
32     a=[xg(r1,r2)]; a=a(:);
33     b=[xg(r3,r4)]; b=b(:);
34     i1=find(a== 180 & b < 0);
35     i2=find(a==-180 & b > 0);
36     a(i1)= -180;
37     a(i2)= 180;
38     i3=find(b== 180 & a < 0);
39     i4=find(b==-180 & a > 0);
40     b(i3)= -180;
41     b(i4)= 180;
42     c=[yg(r1,r2)]; c=c(:);
43     d=[yg(r3,r4)]; d=d(:);
44     i1=find(c== 90 );
45     i2=find(d== 90 );
46     a(i1)=b(i1);
47     b(i2)=a(i2);
48     i1=find(c== -90 );
49     i2=find(d== -90 );
50     a(i1)=b(i1);
51     b(i2)=a(i2);
52    
53     xcoord(1,nlo:4:nhi)=a;
54     xcoord(2,nlo:4:nhi)=b;
55     ycoord(1,nlo:4:nhi)=c;
56     ycoord(2,nlo:4:nhi)=d;
57    
58     % Store associated index
59     ai=[iind(r1,r4)]; ai=ai(:);
60     bi=[jind(r1,r4)]; bi=bi(:);
61     xcoordi((nlo-1)/4+1:nhi/4)=ai;
62     ycoordi((nlo-1)/4+1:nhi/4)=bi;
63    
64     % Set i+1,j -> i+1,j+1 segment
65     r1=[ilo+1:ihi ];r2=[jlo :jhi-1];
66     r3=[ilo+1:ihi ];r4=[jlo+1:jhi ];
67     a=[xg(r1,r2)]; a=a(:);
68     b=[xg(r3,r4)]; b=b(:);
69     i1=find(a== 180 & b < 0);
70     i2=find(a==-180 & b > 0);
71     a(i1)= -180;
72     a(i2)= 180;
73     i3=find(b== 180 & a < 0);
74     i4=find(b==-180 & a > 0);
75     b(i3)= -180;
76     b(i4)= 180;
77     c=[yg(r1,r2)]; c=c(:);
78     d=[yg(r3,r4)]; d=d(:);
79     i1=find(c== 90 );
80     i2=find(d== 90 );
81     a(i1)=b(i1);
82     b(i2)=a(i2);
83     i1=find(c== -90 );
84     i2=find(d== -90 );
85     a(i1)=b(i1);
86     b(i2)=a(i2);
87    
88     xcoord(1,nlo+1:4:nhi+1)=a;
89     xcoord(2,nlo+1:4:nhi+1)=b;
90     ycoord(1,nlo+1:4:nhi+1)=c;
91     ycoord(2,nlo+1:4:nhi+1)=d;
92    
93     % Set i+1,j+1 -> i ,j+1 segment
94     r1=[ilo+1:ihi ];r2=[jlo+1:jhi ];
95     r3=[ilo :ihi-1];r4=[jlo+1:jhi ];
96     a=[xg(r1,r2)]; a=a(:);
97     b=[xg(r3,r4)]; b=b(:);
98     i1=find(a== 180 & b < 0);
99     i2=find(a==-180 & b > 0);
100     a(i1)= -180;
101     a(i2)= 180;
102     i3=find(b== 180 & a < 0);
103     i4=find(b==-180 & a > 0);
104     b(i3)= -180;
105     b(i4)= 180;
106     c=[yg(r1,r2)]; c=c(:);
107     d=[yg(r3,r4)]; d=d(:);
108     i1=find(c== 90 );
109     i2=find(d== 90 );
110     a(i1)=b(i1);
111     b(i2)=a(i2);
112     i1=find(c== -90 );
113     i2=find(d== -90 );
114     a(i1)=b(i1);
115     b(i2)=a(i2);
116    
117     xcoord(1,nlo+2:4:nhi+2)=a;
118     xcoord(2,nlo+2:4:nhi+2)=b;
119     ycoord(1,nlo+2:4:nhi+2)=c;
120     ycoord(2,nlo+2:4:nhi+2)=d;
121    
122    
123     % Set i ,j+1 -> i ,j segment
124     r1=[ilo :ihi-1];r2=[jlo+1:jhi ];
125     r3=[ilo :ihi-1];r4=[jlo :jhi-1];
126     a=[xg(r1,r2)]; a=a(:);
127     b=[xg(r3,r4)]; b=b(:);
128     i1=find(a== 180 & b < 0);
129     i2=find(a==-180 & b > 0);
130     a(i1)= -180;
131     a(i2)= 180;
132     i3=find(b== 180 & a < 0);
133     i4=find(b==-180 & a > 0);
134     b(i3)= -180;
135     b(i4)= 180;
136     c=[yg(r1,r2)]; c=c(:);
137     d=[yg(r3,r4)]; d=d(:);
138     i1=find(c== 90 );
139     i2=find(d== 90 );
140     a(i1)=b(i1);
141     b(i2)=a(i2);
142     i1=find(c== -90 );
143     i2=find(d== -90 );
144     a(i1)=b(i1);
145     b(i2)=a(i2);
146    
147     xcoord(1,nlo+3:4:nhi+3)=a;
148     xcoord(2,nlo+3:4:nhi+3)=b;
149     ycoord(1,nlo+3:4:nhi+3)=c;
150     ycoord(2,nlo+3:4:nhi+3)=d;
151    
152     nc=nc+(ihi-ilo)*(jhi-jlo)*4;
153    
154     fprintf('Segment calculation for face %d\n',t);
155    
156     end
157    
158     csg.xcoord=xcoord;
159     csg.ycoord=ycoord;
160    
161    
162     px=cast(ones(5,nc/4),class(xcoord));
163     py=cast(ones(5,nc/4),class(ycoord));
164     r=1:4:nc;
165     px(1,:)=xcoord(1,r );
166     px(2,:)=xcoord(2,r );
167     px(3,:)=xcoord(2,r+1);
168     px(4,:)=xcoord(2,r+2);
169     px(5,:)=xcoord(2,r+3);
170     py(1,:)=ycoord(1,r );
171     py(2,:)=ycoord(2,r );
172     py(3,:)=ycoord(2,r+1);
173     py(4,:)=ycoord(2,r+2);
174     py(5,:)=ycoord(2,r+3);
175    
176     % For any polygon with a negative longitude set longitude 180
177     % to -180 where the polygon lonigtude is exactly 180.
178     [ r0, c0]=find(px< 0 );
179     [r180,c180]=find(px==180);
180     [cfix,ia,ib]=intersect(c180,c0);
181     rfix=r180(ia);
182     foo=px(:,cfix);
183     foo(foo==180)=-180;
184     px(:,cfix)=foo;
185    
186     % Set longitude to 0 for all points that are at the pole.
187     px(py== 90)=0;
188     px(py==-90)=0;
189    
190     csg.px=px;
191     csg.py=py;
192     csg.xcoordi=xcoordi;
193     csg.ycoordi=ycoordi;
194    
195     return

  ViewVC Help
Powered by ViewVC 1.1.22