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

Contents 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 - (show annotations) (download)
Tue Mar 15 18:05:55 2005 UTC (19 years ago) by cnh
Branch: MAIN
CVS Tags: HEAD
Adding cube face plotting utility + t32 example that it uses.

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