matlab_class/gcmfaces_misc/gcmfaces_interp_coeffs.m

function [prof_interp,tile_corners]=gcmfaces_interp_coeffs(prof_lon,prof_lat,varargin);
%[prof_interp]=gcmfaces_interp_coeffs(prof_lon,prof_lat);
%object:    compute bilinear interpolation weights for prof_lon, prof_lat
%inputs:    prof_lon, prof_lat are column vectors
%(optional) ni,nj is the MITgcm tile size (2 numbers total)
%outputs:   prof_interp contains face and tile numbers,
%           indices within tiles (within 0:ni+1,0:nj+1)
%           and interpolation weights (between 0 and 1)
%           of the four neighboring grid points.
%(optional) tile_corners contains XC11,XCNINJ,YC11,YCNINJ (for MITprof)
%
%note: pathological cases (e.g. at edges) remain to be treated.
%example:
%prof=MITprof_load('ctd_feb2013.nc');
%[prof_interp]=gcmfaces_interp_coeffs(prof.prof_lon,prof.prof_lat);

gcmfaces_global;

doDisplay=0;
doVerbose=0; if myenv.verbose>=1; doVerbose=myenv.verbose; end;
%set doVerbose to display points that could not be triangulated (if any)

%set-up tile information (domain decomposition to ni,nj blocs)
if nargin<=2;
    ni=30; nj=30;
else;
    ni=varargin{1};
    nj=varargin{2};
end;

map_tile=gcmfaces_loc_tile(ni,nj);
loc_tile=gcmfaces_loc_tile(ni,nj,prof_lon,prof_lat);
prof_tile=loc_tile.tileNo;
list_tile=unique(prof_tile);

%initialize output:
prof_interp.point=loc_tile.point;  
prof_interp.face=NaN*prof_lon;
prof_interp.tile=NaN*prof_lon;
prof_interp.i=NaN*repmat(prof_lon,[1 4]);
prof_interp.j=NaN*repmat(prof_lon,[1 4]);
prof_interp.w=NaN*repmat(prof_lon,[1 4]);
prof_interp.XC=NaN*repmat(prof_lon,[1 4]);
prof_interp.YC=NaN*repmat(prof_lon,[1 4]);
%
tile_corners.XC11=NaN*prof_lon;
tile_corners.YC11=NaN*prof_lon;
tile_corners.XCNINJ=NaN*prof_lon;
tile_corners.YCNINJ=NaN*prof_lon;

%loop over tiles
for ii=1:length(list_tile);
    %1) determine face of current tile ...
    tmp1=1*(map_tile==list_tile(ii));
    tmp11=sum(sum(tmp1,1),2); tmp12=[];
    for ff=1:tmp11.nFaces; tmp12=[tmp12,tmp11{ff}]; end;
    iiFace=find(tmp12);
    %... and its index range within face ...
    tmp1=tmp1{iiFace};
    tmp11=sum(tmp1,2);
    iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
    tmp11=sum(tmp1,1);
    jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
    %... as well as the list of profiles in tile
    ii_prof=find(prof_tile==list_tile(ii));
    %tile corners
    XC11=mygrid.XC{iiFace}(iiMin,jjMin);
    YC11=mygrid.YC{iiFace}(iiMin,jjMin);
    XCNINJ=mygrid.XC{iiFace}(iiMax,jjMax);
    YCNINJ=mygrid.YC{iiFace}(iiMax,jjMax);

    clear tmp*;
    
    %2) stereographic projection to current tile center:
    ii0=floor((iiMin+iiMax)/2);
    jj0=floor((jjMin+jjMax)/2);
    XC0=mygrid.XC{iiFace}(ii0,jj0);
    YC0=mygrid.YC{iiFace}(ii0,jj0);
    %for grid locations:
    [xx,yy]=gcmfaces_stereoproj(XC0,YC0);
    %for profile locations
    [prof_x,prof_y]=gcmfaces_stereoproj(XC0,YC0,prof_lon,prof_lat);
    
    % nrm=sqrt(prof_x.^2+prof_y.^2);
    %ii_prof=find(nrm<tan(pi/4/2));%points inside of pi/4 cone
    
    %3) form array of grid cell quadrilaterals
    xxx=exch_T_N(xx); yyy=exch_T_N(yy);
    xc=exch_T_N(mygrid.XC); yc=exch_T_N(mygrid.YC);
    x_quad=[]; y_quad=[]; xc_quad=[]; yc_quad=[]; i_quad=[]; j_quad=[];
    for pp=1:4;
        switch pp;
            case 1; di=0; dj=0;
            case 2; di=1; dj=0;
            case 3; di=1; dj=1;
            case 4; di=0; dj=1;
        end;
        %note the shift in indices due to exchange above
        tmpx=xxx{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpx=tmpx(:); x_quad=[x_quad tmpx];
        tmpy=yyy{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpy=tmpy(:); y_quad=[y_quad tmpy];
        %
        tmpx=xc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpx=tmpx(:); xc_quad=[xc_quad tmpx];
        tmpy=yc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
        tmpy=tmpy(:); yc_quad=[yc_quad tmpy];
        %
        tmpi=[0+di:iiMax-iiMin+1+di]'*ones(1,jjMax-jjMin+2);
        tmpi=tmpi(:); i_quad=[i_quad tmpi];
        tmpj=ones(jjMax-jjMin+2,1)*[0+dj:jjMax-jjMin+1+dj];
        tmpj=tmpj(:); j_quad=[j_quad tmpj];
    end;
    
    %4) associate profile locations with quadrilaterals
    [angsum]=gcmfaces_polygonangle(x_quad,y_quad,prof_x(ii_prof),prof_y(ii_prof));
    [II,JJ]=find(abs(angsum)>179);%+-360 for an interior point (+-180 for an edge point)
    if length(unique(JJ))~=length(JJ)&doVerbose;
        n0=num2str(length(JJ)-length(unique(JJ)));
        warning(['multiple polygons (' n0 ')']);
        %store indices of such instances for display
        [a,b]=hist(JJ,unique(JJ));
        KK=find(a>1);
        kk_prof=ii_prof(KK);
        kk_quad={};
        for kk=1:length(KK);
          kk_quad{kk}=II(find(JJ==KK(kk)));
        end
    else;
        kk_prof=[];
    end;
    if length(unique(JJ))<length(ii_prof)&doVerbose;
        n0=num2str(length(ii_prof)-length(unique(JJ)));
        n1=num2str(length(ii_prof));
        warning(['no polygon for ' n0 ' / ' n1]);
        %the following will then remove the corresponding profiles form ii_prof
    end;
    [C,IA,IC] = unique(JJ);
    %
    ii_prof0=ii_prof;
    ii_prof=ii_prof(C);%treated profiles
    jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles
    %
    ii_quad=II(IA);

    if length(kk_prof)>0&doVerbose>=3;
      for kk=1:length(kk_prof);
        figureL;
        tmpx=x_quad(kk_quad{kk},[1:4 1])';
        tmpy=y_quad(kk_quad{kk},[1:4 1])';
        plot(tmpx,tmpy,'k.-','MarkerSize',36); hold on;
        plot(prof_x(kk_prof(kk)),prof_y(kk_prof(kk)),'r.','MarkerSize',36)
        aa=axis;
        aa(1:2)=aa(1:2)+abs(diff(aa(1:2)))*[-0.1 0.1];
        aa(3:4)=aa(3:4)+abs(diff(aa(3:4)))*[-0.1 0.1];
        axis(aa);
        keyboard;
      end;
    end;

    if length(jj_prof)>0&doVerbose>=2;
        figureL;
        %
        subplot(2,1,1);
        plot(prof_x(ii_prof),prof_y(ii_prof),'.');
        hold on; plot(x_quad(:),y_quad(:),'r.');
        plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',36);
        %
        subplot(2,1,2);
        tmpx=convert2vector(mygrid.XC);
        tmpy=convert2vector(mygrid.YC);
        tmpi=convert2vector(map_tile);
        tmpi=find(tmpi==ii);
        plot(tmpx(:),tmpy(:),'r.'); hold on; 
        plot(tmpx(tmpi),tmpy(tmpi),'c.');
        plot(prof_lon(ii_prof),prof_lat(ii_prof),'.');
        plot(prof_lon(jj_prof),prof_lat(jj_prof),'k.','MarkerSize',36);
        %
        tmp1=prof_lon([ii_prof;jj_prof]);
        tmp2=prof_lat([ii_prof;jj_prof]);
        axis([min(tmp1) max(tmp1) min(tmp2) max(tmp2)]);
        %
        keyboard;
    end;
    
    if doDisplay;
        figureL;
        xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
        yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
        pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on;
        cc=caxis; cc(2)=cc(2)*2; caxis(cc);
        plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20);
        plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60);
        axis([-0.6 0.6 -0.6 0.6]/6);
    end;
    
    if ~isempty(ii_prof);
        %5) determine bi-linear interpolation weights:
        px=x_quad(ii_quad,:);
        py=y_quad(ii_quad,:);
        ox=prof_x(ii_prof);
        oy=prof_y(ii_prof);
        [ow]=gcmfaces_quadmap(px,py,ox,oy);
        
        %to double check interpolation
        % pw=squeeze(ow);
        % oxInterp=sum(pw.*px,2);
        % oyInterp=sum(pw.*py,2);

        %round up coefficient to 4th digit (also to avoid slight negatives)
        test1=~isempty(find(ow(:)<-1e-5));
        if test1; error('interp weight < 0 -- something went wrong'); end;
        test1=~isempty(find(ow(:)>1+1e-5));
        if test1; error('interp weight >1 -- something went wrong'); end;
        %
        ow=1e-4*round(ow*1e4);
        sumw=repmat(sum(ow,3),[1 1 4]);
        ow=ow./sumw;
        
        %6) output interpolation specs:
        prof_interp.face(ii_prof,1)=iiFace*(1+0*ii_quad);
        prof_interp.tile(ii_prof,1)=list_tile(ii)*(1+0*ii_quad);
        prof_interp.i(ii_prof,:)=i_quad(ii_quad,:);
        prof_interp.j(ii_prof,:)=j_quad(ii_quad,:);
        prof_interp.w(ii_prof,:)=squeeze(ow);
        prof_interp.XC(ii_prof,:)=xc_quad(ii_quad,:);
        prof_interp.YC(ii_prof,:)=yc_quad(ii_quad,:);
        %
        tile_corners.XC11(ii_prof)=XC11;
        tile_corners.YC11(ii_prof)=YC11;
        tile_corners.XCNINJ(ii_prof)=XCNINJ;
        tile_corners.YCNINJ(ii_prof)=YCNINJ;
    end;
    
end;%for ii=1:length(list_tile);

if doVerbose;
  n1=sum(~isnan(prof_interp.face));
  n2=sum(isnan(prof_interp.face));
  fprintf(['interpolated points: ' num2str(n1) '\n']);
  fprintf(['un-treated points:   ' num2str(n2) '\n']);
end;

1	gforget	1.2	function [prof_interp,tile_corners]=gcmfaces_interp_coeffs(prof_lon,prof_lat,varargin);
2	gforget	1.1	%[prof_interp]=gcmfaces_interp_coeffs(prof_lon,prof_lat);
3			%object: compute bilinear interpolation weights for prof_lon, prof_lat
4			%inputs: prof_lon, prof_lat are column vectors
5			%(optional) ni,nj is the MITgcm tile size (2 numbers total)
6			%outputs: prof_interp contains face and tile numbers,
7			% indices within tiles (within 0:ni+1,0:nj+1)
8			% and interpolation weights (between 0 and 1)
9			% of the four neighboring grid points.
10	gforget	1.2	%(optional) tile_corners contains XC11,XCNINJ,YC11,YCNINJ (for MITprof)
11	gforget	1.1	%
12			%note: pathological cases (e.g. at edges) remain to be treated.
13			%example:
14			%prof=MITprof_load('ctd_feb2013.nc');
15			%[prof_interp]=gcmfaces_interp_coeffs(prof.prof_lon,prof.prof_lat);
16
17			gcmfaces_global;
18
19			doDisplay=0;
20	gforget	1.4	doVerbose=0; if myenv.verbose>=1; doVerbose=myenv.verbose; end;
21	gforget	1.5	%set doVerbose to display points that could not be triangulated (if any)
22	gforget	1.1
23			%set-up tile information (domain decomposition to ni,nj blocs)
24			if nargin<=2;
25			ni=30; nj=30;
26			else;
27			ni=varargin{1};
28			nj=varargin{2};
29			end;
30
31			map_tile=gcmfaces_loc_tile(ni,nj);
32			loc_tile=gcmfaces_loc_tile(ni,nj,prof_lon,prof_lat);
33			prof_tile=loc_tile.tileNo;
34			list_tile=unique(prof_tile);
35
36			%initialize output:
37	gforget	1.6	prof_interp.point=loc_tile.point;
38	gforget	1.1	prof_interp.face=NaN*prof_lon;
39			prof_interp.tile=NaN*prof_lon;
40			prof_interp.i=NaN*repmat(prof_lon,[1 4]);
41			prof_interp.j=NaN*repmat(prof_lon,[1 4]);
42			prof_interp.w=NaN*repmat(prof_lon,[1 4]);
43	gforget	1.6	prof_interp.XC=NaN*repmat(prof_lon,[1 4]);
44			prof_interp.YC=NaN*repmat(prof_lon,[1 4]);
45	gforget	1.2	%
46			tile_corners.XC11=NaN*prof_lon;
47			tile_corners.YC11=NaN*prof_lon;
48			tile_corners.XCNINJ=NaN*prof_lon;
49			tile_corners.YCNINJ=NaN*prof_lon;
50	gforget	1.1
51			%loop over tiles
52			for ii=1:length(list_tile);
53			%1) determine face of current tile ...
54			tmp1=1*(map_tile==list_tile(ii));
55			tmp11=sum(sum(tmp1,1),2); tmp12=[];
56			for ff=1:tmp11.nFaces; tmp12=[tmp12,tmp11{ff}]; end;
57			iiFace=find(tmp12);
58			%... and its index range within face ...
59			tmp1=tmp1{iiFace};
60			tmp11=sum(tmp1,2);
61			iiMin=min(find(tmp11)); iiMax=max(find(tmp11));
62			tmp11=sum(tmp1,1);
63			jjMin=min(find(tmp11)); jjMax=max(find(tmp11));
64			%... as well as the list of profiles in tile
65			ii_prof=find(prof_tile==list_tile(ii));
66	gforget	1.2	%tile corners
67			XC11=mygrid.XC{iiFace}(iiMin,jjMin);
68			YC11=mygrid.YC{iiFace}(iiMin,jjMin);
69			XCNINJ=mygrid.XC{iiFace}(iiMax,jjMax);
70			YCNINJ=mygrid.YC{iiFace}(iiMax,jjMax);
71
72	gforget	1.1	clear tmp*;
73
74			%2) stereographic projection to current tile center:
75			ii0=floor((iiMin+iiMax)/2);
76			jj0=floor((jjMin+jjMax)/2);
77			XC0=mygrid.XC{iiFace}(ii0,jj0);
78			YC0=mygrid.YC{iiFace}(ii0,jj0);
79			%for grid locations:
80			[xx,yy]=gcmfaces_stereoproj(XC0,YC0);
81			%for profile locations
82			[prof_x,prof_y]=gcmfaces_stereoproj(XC0,YC0,prof_lon,prof_lat);
83
84			% nrm=sqrt(prof_x.^2+prof_y.^2);
85			%ii_prof=find(nrm<tan(pi/4/2));%points inside of pi/4 cone
86
87			%3) form array of grid cell quadrilaterals
88			xxx=exch_T_N(xx); yyy=exch_T_N(yy);
89	gforget	1.6	xc=exch_T_N(mygrid.XC); yc=exch_T_N(mygrid.YC);
90			x_quad=[]; y_quad=[]; xc_quad=[]; yc_quad=[]; i_quad=[]; j_quad=[];
91	gforget	1.1	for pp=1:4;
92			switch pp;
93			case 1; di=0; dj=0;
94			case 2; di=1; dj=0;
95			case 3; di=1; dj=1;
96			case 4; di=0; dj=1;
97			end;
98			%note the shift in indices due to exchange above
99			tmpx=xxx{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
100			tmpx=tmpx(:); x_quad=[x_quad tmpx];
101			tmpy=yyy{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
102			tmpy=tmpy(:); y_quad=[y_quad tmpy];
103			%
104	gforget	1.6	tmpx=xc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
105			tmpx=tmpx(:); xc_quad=[xc_quad tmpx];
106			tmpy=yc{iiFace}(iiMin+di:iiMax+1+di,jjMin+dj:jjMax+1+dj);
107			tmpy=tmpy(:); yc_quad=[yc_quad tmpy];
108			%
109	gforget	1.1	tmpi=[0+di:iiMax-iiMin+1+di]'*ones(1,jjMax-jjMin+2);
110			tmpi=tmpi(:); i_quad=[i_quad tmpi];
111			tmpj=ones(jjMax-jjMin+2,1)*[0+dj:jjMax-jjMin+1+dj];
112			tmpj=tmpj(:); j_quad=[j_quad tmpj];
113			end;
114
115			%4) associate profile locations with quadrilaterals
116			[angsum]=gcmfaces_polygonangle(x_quad,y_quad,prof_x(ii_prof),prof_y(ii_prof));
117			[II,JJ]=find(abs(angsum)>179);%+-360 for an interior point (+-180 for an edge point)
118			if length(unique(JJ))~=length(JJ)&doVerbose;
119			n0=num2str(length(JJ)-length(unique(JJ)));
120			warning(['multiple polygons (' n0 ')']);
121	gforget	1.4	%store indices of such instances for display
122			[a,b]=hist(JJ,unique(JJ));
123			KK=find(a>1);
124			kk_prof=ii_prof(KK);
125			kk_quad={};
126			for kk=1:length(KK);
127			kk_quad{kk}=II(find(JJ==KK(kk)));
128			end
129			else;
130			kk_prof=[];
131	gforget	1.1	end;
132			if length(unique(JJ))<length(ii_prof)&doVerbose;
133			n0=num2str(length(ii_prof)-length(unique(JJ)));
134			n1=num2str(length(ii_prof));
135			warning(['no polygon for ' n0 ' / ' n1]);
136			%the following will then remove the corresponding profiles form ii_prof
137			end;
138			[C,IA,IC] = unique(JJ);
139			%
140			ii_prof0=ii_prof;
141			ii_prof=ii_prof(C);%treated profiles
142			jj_prof=setdiff(ii_prof0,ii_prof);%un-treated profiles
143			%
144			ii_quad=II(IA);
145	gforget	1.4
146			if length(kk_prof)>0&doVerbose>=3;
147			for kk=1:length(kk_prof);
148			figureL;
149			tmpx=x_quad(kk_quad{kk},[1:4 1])';
150			tmpy=y_quad(kk_quad{kk},[1:4 1])';
151			plot(tmpx,tmpy,'k.-','MarkerSize',36); hold on;
152			plot(prof_x(kk_prof(kk)),prof_y(kk_prof(kk)),'r.','MarkerSize',36)
153			aa=axis;
154			aa(1:2)=aa(1:2)+abs(diff(aa(1:2)))*[-0.1 0.1];
155			aa(3:4)=aa(3:4)+abs(diff(aa(3:4)))*[-0.1 0.1];
156			axis(aa);
157			keyboard;
158			end;
159			end;
160
161			if length(jj_prof)>0&doVerbose>=2;
162			figureL;
163			%
164			subplot(2,1,1);
165			plot(prof_x(ii_prof),prof_y(ii_prof),'.');
166			hold on; plot(x_quad(:),y_quad(:),'r.');
167			plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',36);
168			%
169			subplot(2,1,2);
170			tmpx=convert2vector(mygrid.XC);
171			tmpy=convert2vector(mygrid.YC);
172			tmpi=convert2vector(map_tile);
173			tmpi=find(tmpi==ii);
174			plot(tmpx(:),tmpy(:),'r.'); hold on;
175			plot(tmpx(tmpi),tmpy(tmpi),'c.');
176			plot(prof_lon(ii_prof),prof_lat(ii_prof),'.');
177			plot(prof_lon(jj_prof),prof_lat(jj_prof),'k.','MarkerSize',36);
178			%
179			tmp1=prof_lon([ii_prof;jj_prof]);
180			tmp2=prof_lat([ii_prof;jj_prof]);
181			axis([min(tmp1) max(tmp1) min(tmp2) max(tmp2)]);
182			%
183			keyboard;
184			end;
185	gforget	1.1
186			if doDisplay;
187			figureL;
188			xx_tile=xxx{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
189			yy_tile=yyy{iiFace}(iiMin:iiMax+2,jjMin:jjMax+2);
190			pcolor(xx_tile,yy_tile,sqrt(xx_tile.^2+yy_tile.^2)); hold on;
191			cc=caxis; cc(2)=cc(2)*2; caxis(cc);
192			plot(prof_x(ii_prof),prof_y(ii_prof),'r.','MarkerSize',20);
193			plot(prof_x(jj_prof),prof_y(jj_prof),'k.','MarkerSize',60);
194			axis([-0.6 0.6 -0.6 0.6]/6);
195			end;
196
197			if ~isempty(ii_prof);
198			%5) determine bi-linear interpolation weights:
199			px=x_quad(ii_quad,:);
200			py=y_quad(ii_quad,:);
201			ox=prof_x(ii_prof);
202			oy=prof_y(ii_prof);
203			[ow]=gcmfaces_quadmap(px,py,ox,oy);
204
205			%to double check interpolation
206			% pw=squeeze(ow);
207			% oxInterp=sum(pw.*px,2);
208			% oyInterp=sum(pw.*py,2);
209	gforget	1.3
210			%round up coefficient to 4th digit (also to avoid slight negatives)
211			test1=~isempty(find(ow(:)<-1e-5));
212			if test1; error('interp weight < 0 -- something went wrong'); end;
213			test1=~isempty(find(ow(:)>1+1e-5));
214			if test1; error('interp weight >1 -- something went wrong'); end;
215			%
216			ow=1e-4round(ow1e4);
217			sumw=repmat(sum(ow,3),[1 1 4]);
218			ow=ow./sumw;
219	gforget	1.1
220			%6) output interpolation specs:
221			prof_interp.face(ii_prof,1)=iiFace(1+0ii_quad);
222			prof_interp.tile(ii_prof,1)=list_tile(ii)(1+0ii_quad);
223			prof_interp.i(ii_prof,:)=i_quad(ii_quad,:);
224			prof_interp.j(ii_prof,:)=j_quad(ii_quad,:);
225			prof_interp.w(ii_prof,:)=squeeze(ow);
226	gforget	1.6	prof_interp.XC(ii_prof,:)=xc_quad(ii_quad,:);
227			prof_interp.YC(ii_prof,:)=yc_quad(ii_quad,:);
228	gforget	1.2	%
229			tile_corners.XC11(ii_prof)=XC11;
230			tile_corners.YC11(ii_prof)=YC11;
231			tile_corners.XCNINJ(ii_prof)=XCNINJ;
232			tile_corners.YCNINJ(ii_prof)=YCNINJ;
233	gforget	1.1	end;
234
235			end;%for ii=1:length(list_tile);
236
237	gforget	1.5	if doVerbose;
238			n1=sum(~isnan(prof_interp.face));
239			n2=sum(isnan(prof_interp.face));
240			fprintf(['interpolated points: ' num2str(n1) '\n']);
241			fprintf(['un-treated points: ' num2str(n2) '\n']);
242			end;
243	gforget	1.1