function csg=make_cs_segments_and_patches(csg)
% 
% Add line segments and polygon patches to CS grid object.
%

clf
csg
xg=csg.gridarr(:,:,csg.xgpos);
yg=csg.gridarr(:,:,csg.ygpos);
iind=csg.index_i_center;
jind=csg.index_j_center;
clear xc yc;
nc=0;
np=0;
ng=csg.ng;
xcoord=cast(ones(2,ng*ng*6*4),class(xg));
ycoord=cast(ones(2,ng*ng*6*4),class(xg));

fprintf('Beginning segment calculation\n');

for t=1:numel(csg.fx)
 ilo=csg.ilog(t);
 jlo=csg.jlog(t);
 ihi=csg.ilog(t)+csg.fx(t);
 jhi=csg.jlog(t)+csg.fy(t);
 nlo=nc+1;
 nhi=nc+(ihi-ilo)*(jhi-jlo)*4;
 
 % Set i  ,j   -> i+1,j   segment
 r1=[ilo  :ihi-1];r2=[jlo:jhi-1];
 r3=[ilo+1:ihi  ];r4=[jlo:jhi-1];
 a=[xg(r1,r2)]; a=a(:);
 b=[xg(r3,r4)]; b=b(:);
 i1=find(a== 180 & b < 0);
 i2=find(a==-180 & b > 0);
 a(i1)= -180;
 a(i2)=  180;
 i3=find(b== 180 & a < 0);
 i4=find(b==-180 & a > 0);
 b(i3)= -180;
 b(i4)=  180;
 c=[yg(r1,r2)]; c=c(:);
 d=[yg(r3,r4)]; d=d(:);
 i1=find(c== 90 );
 i2=find(d== 90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 i1=find(c== -90 );
 i2=find(d== -90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 
 xcoord(1,nlo:4:nhi)=a;
 xcoord(2,nlo:4:nhi)=b;
 ycoord(1,nlo:4:nhi)=c;
 ycoord(2,nlo:4:nhi)=d;
 
 % Store associated index
 ai=[iind(r1,r4)]; ai=ai(:);
 bi=[jind(r1,r4)]; bi=bi(:);
 xcoordi((nlo-1)/4+1:nhi/4)=ai;
 ycoordi((nlo-1)/4+1:nhi/4)=bi;
 
 % Set i+1,j   -> i+1,j+1 segment
 r1=[ilo+1:ihi  ];r2=[jlo  :jhi-1];
 r3=[ilo+1:ihi  ];r4=[jlo+1:jhi  ];
 a=[xg(r1,r2)]; a=a(:);
 b=[xg(r3,r4)]; b=b(:);
 i1=find(a== 180 & b < 0);
 i2=find(a==-180 & b > 0);
 a(i1)= -180;
 a(i2)=  180;
 i3=find(b== 180 & a < 0);
 i4=find(b==-180 & a > 0);
 b(i3)= -180;
 b(i4)=  180;
 c=[yg(r1,r2)]; c=c(:);
 d=[yg(r3,r4)]; d=d(:);
 i1=find(c== 90 );
 i2=find(d== 90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 i1=find(c== -90 );
 i2=find(d== -90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 
 xcoord(1,nlo+1:4:nhi+1)=a;
 xcoord(2,nlo+1:4:nhi+1)=b;
 ycoord(1,nlo+1:4:nhi+1)=c;
 ycoord(2,nlo+1:4:nhi+1)=d;
  
 % Set i+1,j+1 -> i  ,j+1 segment
 r1=[ilo+1:ihi  ];r2=[jlo+1:jhi  ];
 r3=[ilo  :ihi-1];r4=[jlo+1:jhi  ];
 a=[xg(r1,r2)]; a=a(:);
 b=[xg(r3,r4)]; b=b(:);
 i1=find(a== 180 & b < 0);
 i2=find(a==-180 & b > 0);
 a(i1)= -180;
 a(i2)=  180;
 i3=find(b== 180 & a < 0);
 i4=find(b==-180 & a > 0);
 b(i3)= -180;
 b(i4)=  180;
 c=[yg(r1,r2)]; c=c(:);
 d=[yg(r3,r4)]; d=d(:);
 i1=find(c== 90 );
 i2=find(d== 90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 i1=find(c== -90 );
 i2=find(d== -90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 
 xcoord(1,nlo+2:4:nhi+2)=a;
 xcoord(2,nlo+2:4:nhi+2)=b;
 ycoord(1,nlo+2:4:nhi+2)=c;
 ycoord(2,nlo+2:4:nhi+2)=d;
  
 
 % Set i  ,j+1 -> i  ,j   segment
 r1=[ilo  :ihi-1];r2=[jlo+1:jhi  ];
 r3=[ilo  :ihi-1];r4=[jlo  :jhi-1];
 a=[xg(r1,r2)]; a=a(:);
 b=[xg(r3,r4)]; b=b(:);
 i1=find(a== 180 & b < 0);
 i2=find(a==-180 & b > 0);
 a(i1)= -180;
 a(i2)=  180;
 i3=find(b== 180 & a < 0);
 i4=find(b==-180 & a > 0);
 b(i3)= -180;
 b(i4)=  180;
 c=[yg(r1,r2)]; c=c(:);
 d=[yg(r3,r4)]; d=d(:);
 i1=find(c== 90 );
 i2=find(d== 90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 i1=find(c== -90 );
 i2=find(d== -90 );
 a(i1)=b(i1);
 b(i2)=a(i2);
 
 xcoord(1,nlo+3:4:nhi+3)=a;
 xcoord(2,nlo+3:4:nhi+3)=b;
 ycoord(1,nlo+3:4:nhi+3)=c;
 ycoord(2,nlo+3:4:nhi+3)=d;
  
 nc=nc+(ihi-ilo)*(jhi-jlo)*4;
 
 fprintf('Segment calculation for face %d\n',t);
 
end

csg.xcoord=xcoord;
csg.ycoord=ycoord;


px=cast(ones(5,nc/4),class(xcoord));
py=cast(ones(5,nc/4),class(ycoord));
r=1:4:nc;
px(1,:)=xcoord(1,r  );
px(2,:)=xcoord(2,r  );
px(3,:)=xcoord(2,r+1);
px(4,:)=xcoord(2,r+2);
px(5,:)=xcoord(2,r+3);
py(1,:)=ycoord(1,r  );
py(2,:)=ycoord(2,r  );
py(3,:)=ycoord(2,r+1);
py(4,:)=ycoord(2,r+2);
py(5,:)=ycoord(2,r+3);

% For any polygon with a negative longitude set longitude 180
% to -180 where the polygon lonigtude is exactly 180.
[  r0,  c0]=find(px< 0  );
[r180,c180]=find(px==180);
[cfix,ia,ib]=intersect(c180,c0);
rfix=r180(ia);
foo=px(:,cfix);
foo(foo==180)=-180;
px(:,cfix)=foo;

% Set longitude to 0 for all points that are at the pole.
px(py== 90)=0;
px(py==-90)=0;

csg.px=px;
csg.py=py;
csg.xcoordi=xcoordi;
csg.ycoordi=ycoordi;

return