function [IX,M] = make_geometry(ax,bx,xc,yc,R,GEOM) 
% function [IX,M] = make_geometry(ax,bx,xc,yc,R,GEOM)
% pre-processor for the interface structures
% INPUT: ax, bx: boundaries of the embedding box
%        xc, yc: origin of the inclusion
%        GEOM  : shape of the inclusion, 'circle' or 'flower'
% OUTPUT: interface structure IX with fields
%         coord: x and y coordinates of the intersections
%         type : 'x' or 'y'. For later purposes they have to be
%                 ordered, at first come all 'x'-type intersections
%                 and then all 'y'-type.
%         anch: indices of the achnor points   
%         n   : normal
%         t   : tangent
%               The pair (n,t) forms the positively oriented
%               coordinate system 
% based on the codes by Andreas Wiegmann
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global GRID
  
ay = ax; by = bx;

% function from ALW
[xp,yp,M,th,flow_k,flow_w,flow_c,flow_r0] = make_interface(GEOM,R,xc,yc,GRID);

% routines of ALW take only normalized geometries=>scale everything
l = bx - ax; % square already before assumed
hh = GRID.hx/l; % scaled mesh

pars = [hh,GRID.nx-1,GRID.ny-1];

% scale & shift also interface
xp = (xp-ax)/l; yp = (yp-ay)/l;

% interpolate also values of theta (parameter for circle & flower)
% inside/outside indicator (dummy1) has flipud problem, so 
% because of this do not use it. Use the analytically known
% M instead
[Inti,Intv,lo,ki] = geometry(pars,[xp yp],[th zeros(size(yp))],4); 
IX_theta = Intv(:,7);

% scale back the coordinates
ix = Intv(:,1)*l  + ax; iy = Intv(:,2)*l  + ay;
IX.coord = [ix iy];

% ALW provides lower left grid point as anchor
% transform to inner 
achx = Inti(:,2); 
achy = Inti(:,3);

IN = length(ix);

IX.type  = Inti(:,1);

for ctr = 1:IN  
  if(M(achx(ctr),achy(ctr))==0)
    if IX.type(ctr)==2
      % x-intersection
      achx(ctr) = achx(ctr)+1;
    else
      % y-intersection
      achy(ctr) = achy(ctr)+1;
    end
  end
end

IX.anch  = [achx achy];
IX.n     = [Intv(:,4) -Intv(:,3)];
IX.t     = [Intv(:,3) Intv(:,4)];


% finally, need to reorder everything
y_inters = find(IX.type==1);
x_inters = find(IX.type==2);
INx      = length(x_inters);


IX.coord = IX.coord([x_inters;y_inters],:);
IX.n     = IX.n([x_inters;y_inters],:);
IX.t     = IX.t([x_inters;y_inters],:);
IX.anch  = IX.anch([x_inters;y_inters],:); 
IX.type(1:INx) = 'x';
IX.type(INx+1:end) = 'y';
% put type to be a character!
IX.type = char(IX.type);


IX_theta = IX_theta([x_inters;y_inters]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DOUBLE INTERSECTION CHECK - those must be skipped
to_skip = zeros(IN,1);
for p = 1:INx  
  hminus = IX.coord(p,1) - GRID.x(IX.anch(p,1));
  
  if ((hminus>1e-10)&(hminus<GRID.hx-1e-10))&(M(IX.anch(p,1)+1,IX.anch(p,2))==1)
    % interface is to the right, but the next
    % neighbour to the right is still in the domain
    % => double intersection
    disp('skip the double x-intersection!')
    to_skip(p) = 1;
  elseif ((hminus<-1e-10)&(hminus>-GRID.hx+1e-10))&(M(IX.anch(p,1)-1,IX.anch(p,2))==1)
    disp('skip the double x-intersection!')
    to_skip(p) = 1;    
  end
  
  if abs(hminus)<1e-10
    % set intersection point to be exactly at the grid point
    IX.coord(p,1) = GRID.x(IX.anch(p,1));
    % set the grid point to be inside
    M(IX.anch(p,1),IX.anch(p,2))=1;
  end
  
  if hminus>GRID.hx-(1e-10)
    % set intersection point to be exactly at the NEXT grid point
    IX.coord(p,1) = GRID.x(IX.anch(p,1)+1);
    % set the grid point to be inside
    M(IX.anch(p,1)+1,IX.anch(p,2))=1;
    % reassign the anchor:
    IX.anch(p,1) = IX.anch(p,1)+1;
  elseif hminus<-GRID.hx+(1e-10)
    % set intersection point to be exactly at the PREVIOUS grid point
    IX.coord(p,1) = GRID.x(IX.anch(p,1)-1);
    % set the grid point to be inside
    M(IX.anch(p,1)-1,IX.anch(p,2))=1;
    % reassign the anchor:
    IX.anch(p,1) = IX.anch(p,1)-1;   
  end
    
  % x-intersections
  if M(IX.anch(p,1),IX.anch(p,2))==0
    % anchor is outside!
    to_skip(p) = 1;
  end  
  
end

for p = INx+1:IN
  hminus = IX.coord(p,2) - GRID.y(IX.anch(p,2));
  
  if ((hminus>1e-10)&(hminus<GRID.hy-1e-10))&(M(IX.anch(p,1),IX.anch(p,2)+1)==1)
    % interface is to the right, but the next
    % neighbour to the right is still in the domain
    % => double intersection
    disp('skip the double y-intersection!')
    to_skip(p) = 1;    
  elseif ((hminus<-1e-10)&(hminus>-GRID.hy+1e-10))&(M(IX.anch(p,1),IX.anch(p,2)-1)==1)
    disp('skip the double y-intersection!')
    to_skip(p) = 1;    
  end  

  if abs(hminus)<1e-10
    % set intersection point to be exactly at the grid point
    IX.coord(p,2) = GRID.y(IX.anch(p,2));
    % set the grid point to be inside
    M(IX.anch(p,1),IX.anch(p,2))=1;
  end  
  
  if hminus>GRID.hy-1e-10
    % set the intersection to be at the NEXT grid point
    IX.coord(p,2) = GRID.y(IX.anch(p,2)+1);
    % set the grid point to be inside
    M(IX.anch(p,1),IX.anch(p,2)+1)=1;
    % reassign the anchor value
    IX.anch(p,2) = IX.anch(p,2)+1;
  elseif hminus<-GRID.hy+1e-10
    % set the intersection to be at the PREVIOUS grid point
    IX.coord(p,2) = GRID.y(IX.anch(p,2)-1);
    % set the grid point to be inside
    M(IX.anch(p,1),IX.anch(p,2)-1)=1;
    % reassign the anchor value
    IX.anch(p,2) = IX.anch(p,2)-1;
  end
    
  if M(IX.anch(p,1),IX.anch(p,2))==0
    % anchor is outside!
    to_skip(p) = 1;
  end  
end

to_skip = find(to_skip==1);

IX.coord(to_skip,:) = [];
IX.anch(to_skip,:)  = [];
IX.n(to_skip,:)     = [];
IX.t(to_skip,:)     = [];
IX.type(to_skip,:)  = [];

IX_theta(to_skip) = [];

ix(to_skip) = [];
iy(to_skip) = [];


IN = length(IX.type);

% rewrite values of IX.n and IX.t by analytic values
% possible for parametrically given interfaces
if isequal(GEOM,'circle')
  n1 = (IX.coord(:,1)-xc)./R;
  n2 = (IX.coord(:,2)-yc)./R;    
  
  IX.n = [n1 n2];
  IX.t = [-n2 n1];
  
elseif isequal(GEOM,'flower')
  dx = -flow_k*(flow_r0+flow_c*sin(flow_w*IX_theta)).*sin(IX_theta)+...
       flow_k*flow_c*flow_w*cos(flow_w*IX_theta).*cos(IX_theta);
  dy = flow_k*(flow_r0+flow_c*sin(flow_w*IX_theta)).*cos(IX_theta)+...
       flow_k*flow_c*flow_w*cos(flow_w*IX_theta).*sin(IX_theta);
  
  the_length = sqrt(dx.^2+dy.^2);
  
  IX.t(:,1) = dx./the_length;
  IX.t(:,2) = dy./the_length;
  
  IX.n(:,1) = IX.t(:,2);
  IX.n(:,2) = -IX.t(:,1);
else
  warning('No analytic geometry calculations!')
end

%draw_interface
% part of test_arc
% plot the interface

figure('name',['Interface, n=',num2str(GRID.n)]);
subplot(1,2,1)
hold on
plot(IX.coord(:,1),IX.coord(:,2),'r.');
plot(GRID.x(IX.anch(:,1)),GRID.y(IX.anch(:,2)),'bo');
quiver(IX.coord(:,1),IX.coord(:,2),IX.n(:,1),IX.n(:,2),'g');
quiver(IX.coord(:,1),IX.coord(:,2),IX.t(:,1),IX.t(:,2),'m');
plot([ax bx bx ax ax],[ay ay by by ay],'k');
axis equal
grid on
title('interface')

subplot(1,2,2)
PP = mesh(GRID.X,GRID.Y,double(M));
set(PP,'linestyle','none')
set(PP,'marker','*')
set(PP,'facecolor','none')
title('indicator')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xp,yp,M,th,flow_k,flow_w,flow_c,flow_r0] = ...
      make_interface(GEOM,R,xc,yc,GRID)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isequal(GEOM,'circle')
  % make the data points, based on which splines will be later fit
  n1  = 5201; % number of points along the boundary    
  
  % circle  
  t0 =  1e-8; % arg shift to avoid intersection trouble    
  th  = [0:2*pi/(n1-1):2*pi]';    
  xp  = xc + R*cos(th +t0);    
  yp  = yc + R*sin(th +t0);      
  
  % indicator for '+' and '-' sides
  % set M to be 1 in '-' side (i.e., to the left)
  M = (((GRID.X-xc).^2+(GRID.Y-yc).^2)<=R^2);  
  
  % the flower variables are set to dummies to 
  % run this toy code also for Matlab 7
  flow_k = 0;
  flow_w = 0;
  flow_c = 0;
  flow_r0 = 0;
  
elseif isequal(GEOM,'flower')
  % flower-shape inclusion from Li98
  n1  = 3000; % number of points along the boundary    
  
  % circle  
  t0 =  1e-8; % arg shift to avoid intersection trouble    
  th  = [0:2*pi/(n1-1):2*pi]';    

  th = th+t0;
  
  flow_r0 = R;
  
  % perturbation parameters for the flower
  flow_k = 0.55;
  flow_w = 5;
  flow_c = 0.15;
  
  %r0 = 0.55;
  %r0 = 0.35+1e-5;    
  ww = flow_w;
  r = flow_r0 + flow_c*sin(ww*th);
  xp = flow_k*r.*cos(th) + xc;
  yp = flow_k*r.*sin(th) + yc;
 
  % polar coordinates
  RR = sqrt((GRID.Y-yc).^2 + (GRID.X-xc).^2);

  % switch off the warnings, the RR can be zero if the 
  % origin belongs to the grid
  warning off
  TH = sign(GRID.Y-yc).*acos((GRID.X-xc)./RR) + 2*pi*((GRID.Y-yc)<0);
  warning on 
  
  M = (RR.^2<=(0.55^2)*(flow_r0 + 0.15*sin(ww*TH)).^2);  
  
  % set the origin (RR=0) also to belong to M
  ORIG    = find(RR==0);
  M(ORIG) = 1;
  
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Inti,Intv,lo,ki] = geometry(pars,datapoints,datavalues,point); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1996-1999.
  
XS                    = datapoints(:,1); 
YS                    = datapoints(:,2); 
 
[Intv,Inti,lo,ki]     = intersections(datapoints,datavalues,pars);
XP                    = Intv(:,1); 
YP                    = Intv(:,2); 
[Intx,Inty,IntD,Inti] = interval(Intv(:,1:2),Inti,pars,point);  
h                     = pars(1); 
nx                    = pars(2); 
ny                    = pars(3); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [intv,inti,lo,ki] = intersections(xy_in,p_in,pars); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1995-1999.

global myeps 
myeps   = 1.0e-12; % Used for fudging in xloop, yloop, also  
                   % in Cardano's formula. 
   
h       = pars(1); 
nx      = pars(2); 
ny      = pars(3); 
xo      = xy_in(:,1)'; 
yo      = xy_in(:,2)'; 
n1      = size(xo,2); 
do(1,:) = diff(xo); 
do(2,:) = diff(yo); 
lo      = [0,cumsum(sqrt(do(1,:).^2+do(2,:).^2))]; 

[brx,coefx,lx,kx] = unmkpp( spline(lo,xo)); 
[bry,coefy,ly,ky] = unmkpp( spline(lo,yo)); 

Sol = zeros(4*(nx+ny),9); % x,y,x,y,k 
                          % t: parameter value, auxiliary 
                          % i: index lower breakpoint, auxiliary 
[minx,maxx,derivx] = xextrema(lx,coefx,coefy,brx,bry); 
[miny,maxy,derivy] = yextrema(lx,coefx,coefy,brx,bry); 
 
mxsol = 1; 

%xloop
for i = 1:lx 
  % fudging to avoid double counting data points on the mesh 
  % On all intervals, always include starting point, not endpoint. 
  % On last interval, include endpoint if not periodic. This should be  
  % fixed, maybe using the information obtained from interval to get  
  % rid of double intersections. Possibly both, if "spike" in graph.  
  fdg1 =   - myeps; 
  if (i == lx)
    fdg2 =   myeps; 
  else 
    fdg2 = - myeps; 
  end 
 
  nmin = round(minx(i)/h);  
  nmin = nmin + (nmin*h<minx(i)); 
  nmin = max(nmin,0); 
  nmax = round(maxx(i)/h); 
  nmax = nmax - (nmax*h>maxx(i)); 
  nmax = min(nmax,nx); 
   
  if xo(i) <  xo(i+1) 
    jvals = nmin: 1:nmax; 
  else 
    jvals = nmax:-1:nmin; 
  end 
  for j = jvals 
    xv = j*h; 
    xs = fcardan(coefx(i,:)-[0,0,0,xv]); 
    ns = 0; 
    for jj = 1:xs(5) 
      ts = xs(jj); 
      if ts > -myeps & ts < brx(i+1)-brx(i)+myeps 
        yv    = polyval(coefy(i,:),ts); 
        if 0 <= yv & yv <= h*ny, 
          xt    = polyval(derivx(i,:),ts); 
          yt    = polyval(derivy(i,:),ts); 
          nrmt  = sqrt(xt^2+yt^2); 
          xtt   = polyval([6*coefx(i,1),2*coefx(i,2)],ts); 
          ytt   = polyval([6*coefy(i,1),2*coefy(i,2)],ts); 
          kappa = (xt*ytt-yt*xtt)/nrmt^3; 
          xt    = xt/nrmt; 
          yt    = yt/nrmt; 
           
          if j > 0 & j < nx , 
            % interior intersection 
            Sol(mxsol,1:9) = [xv,yv,xt,yt,kappa,ts+brx(i),i,1,j+1];  
          else 
            % intersection on boundary 
            Sol(mxsol,1:9) = [xv,yv,xt,yt,kappa,ts+brx(i),i,3,j+1];  
          end 
         
          ns    = ns+1; 
          if mxsol >1 
            if norm(Sol(mxsol,1:2)'-Sol(mxsol-1,1:2)',2) > 100*myeps 
              mxsol = mxsol+1; 
               
            end 
          else 
            mxsol=mxsol+1; 
          end 
        end 
      else 
        % ts skipped 
      end 
    end 
    if ns == 3,  
      error('mesh cannot resolve interface: 3 intersections')  
    end;     
  end 
end 
 
mysol = mxsol; 
for i = 1:ly
  % fudging to avoid double counting data points on the mesh
  % On all intervals, always include starting point, not endpoint.
  % On last interval, include endpoint if not periodic. This should be 
  % fixed, maybe using the information obtained from interval to get 
  % rid of double intersections. Possibly both, if "spike" in graph.  
  fdg1 =   - myeps;
  if (i == ly) 
    fdg2 =   myeps;
  else
    fdg2 = - myeps;
  end

  nmin = round(miny(i)/h); 
  nmin = nmin + (nmin*h<miny(i));
  nmin = max(nmin,0);
  nmax = round(maxy(i)/h);
  nmax = nmax - (nmax*h>maxy(i));
  nmax = min(nmax,ny);
  
  if yo(i) <  yo(i+1)
    jvals = nmin: 1:nmax;
  else
    jvals = nmax:-1:nmin;
  end
  for j = jvals
    yv = j*h;
    ys = fcardan([coefy(i,1:3),coefy(i,4)-yv]);
    ns = 0;
    for jj = 1:ys(5)
      ts = ys(jj);
      if ts  > -myeps & ts < bry(i+1)-bry(i)+myeps
        xv    = polyval(coefx(i,:),ts);
        if 0 <= xv & xv <= h*nx,
          xt    = polyval(derivx(i,:),ts);
          yt    = polyval(derivy(i,:),ts);
          nrmt  = sqrt(xt^2+yt^2);
          xtt   = polyval([6*coefx(i,1),2*coefx(i,2)],ts);
          ytt   = polyval([6*coefy(i,1),2*coefy(i,2)],ts);
          kappa = (xt*ytt-yt*xtt)/nrmt^3;
          xt    = xt/nrmt;
          yt    = yt/nrmt;
          
          if j > 0 & j < ny,
            % interior intersection
            Sol(mysol,1:9) = [xv,yv,xt,yt,kappa,ts+bry(i),i,2,j+1]; 
          else
            % intersection on boundary
            Sol(mysol,1:9) = [xv,yv,xt,yt,kappa,ts+bry(i),i,4,j+1]; 
          end
        
          ns    = ns+1;
          if mysol >1
            if norm(Sol(mysol,1:2)'-Sol(mysol-1,1:2)',2) > 100*myeps
              mysol = mysol+1;
            end
          else
            mysol=mysol+1;
          end
        end 
      else
        % ts skipped
      end
    end
    if ns == 3, 
      error('mesh cannot resolve interface: 3 intersections') 
    end;
  end
end

Sol    = Sol(1:mysol-1,:); 
[Ss,I] = sort(Sol(1:mysol-1,6)); 
Sol    = Sol(I,:); 

%int_findp
np = size(p_in);
if sum(np == [0,0]) == 2,
  pout = [Sol(1:mysol-1,6)];
else
  if np(1) ~= n1,
    error('p input columns need to be same length as xy columns')
  end
  pout = zeros(mysol-1,1+3*np(2));
  pout(:,1) = [Sol(1:mysol-1,6)];
  for j = 1:np(2)
    po = p_in(:,j);
    [brp,coefp,lp,kp] = unmkpp( spline(lo,po));
    jj = (j-1)*4;
    for i = 1:mysol -1
      k = Sol(i,7);
      t = Sol(i,6)-brp(k);
      pout(i,2+jj) = polyval(coefp(k,:),t);
      pout(i,3+jj) = polyval([3*coefp(k,1),2*coefp(k,2),coefp(k,3)],t);
      pout(i,4+jj) = polyval([6*coefp(k,1),2*coefp(k,2)],t);
      pout(i,5+jj) = polyval([6*coefp(k,1)],t);
    end
  end   
end

intv = Sol(1:mysol-1,1:5); % strip zeros and auxiliary variables 
inti = Sol(1:mysol-1,8:9); 
ki   = Sol(1:mysol-1,7); 
intv = [intv,pout]; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Intx,Inty,IntD,Inti] = interval(xy_in,Inti,pars,point);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1995-1999
h    = pars(1);
nx   = pars(2); 
ny   = pars(3);
np   = size(xy_in,1);
x    = xy_in(:,1)';
y    = xy_in(:,2)';
Inti = Inti(:,1:2);

Intx = zeros(ny,nx+1);
Inty = zeros(ny+1,nx);
IntD = zeros(np,4);
IntP =1;

% in case an intersection is also a lattice point, 
% where the curve touches the mesh we are in trouble!

for i = 1:np
  if any([1,3] == Inti(i,1)),
    xi(i) = Inti(i,2);
    yi(i) = round(y(i)/h);
    yi(i) = yi(i) - (yi(i)*h > y(i)) +1;
    Inti(i,3) = yi(i);
  elseif any([2,4] == Inti(i,1)),
    yi(i) = Inti(i,2);
    xi(i) = round(x(i)/h);
    xi(i) = xi(i) - (xi(i)*h > x(i)) +1;
    Inti(i,2) = xi(i);
    Inti(i,3) = yi(i);
  else
    Inti(i,1)
    error('illegal input in interval')
  end
end

for i = 1:np
  if any([1,3] == Inti(i,1)),
    P = Intx(ny+1-yi(i),xi(i));
    if P == 0,
      Intx(ny+1-yi(i),xi(i)) = IntP;
      IntD(IntP,1) = 1;
      IntD(IntP,2) = i;
      IntP = IntP +1;
    elseif IntD(P,1) > 1, 
			[[1:np]',IntD],P,IntP
      error('Too many intersections in interval 1');
    else
      IntD(P,1) = 2;
      IntD(P,3) = i;
    end
  elseif any([2,4] == Inti(i,1)),
    P = Inty(ny+2-yi(i),xi(i));
    if P == 0,
      Inty(ny+2-yi(i),xi(i)) = IntP;
      IntD(IntP,1) = 1;
      IntD(IntP,2) = i;
      IntP = IntP +1;
    elseif IntD(P,1) > 1, 
			[[1:np]',IntD],P,IntP
      error('Too many intersections in interval 2');
    else
      IntD(P,1) = 2;
      IntD(P,3) = i;
    end
  else
    error('  illegal input in interval')
  end
end

% For each interval, choose the interior(4) or 
% exterior(2) irregular point to develop around it.

IntD = IntD(1:IntP-1,:);

for i = 1: IntP-1
  if IntD(i,1) == 1, 
    IntD(i,4) = point;
  else 
    IntD(i,4) = 0;
  end
  for j = 1:IntD(i,1)
    Inti(IntD(i,j+1),4) = i;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [miny,maxy,derivy] = yextrema(ly,coefx,coefy,brx,bry)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1995-1999
  
derivy = zeros(ly,4);
miny   = zeros(1,ly);
maxy   = zeros(1,ly);


coefx = [coefx; 0,0,0,polyval(coefx(ly,:),brx(ly+1)-brx(ly))];
coefy = [coefy; 0,0,0,polyval(coefy(ly,:),brx(ly+1)-brx(ly))];
%yloop
for i = 1:ly
  evs = coefy(i:i+1,4);
  exs = coefx(i:i+1,4);  
  if evs(1) < evs(2),
    mny = evs(1);
    mxy = evs(2);
    mnx = exs(1);
    mxx = exs(2);
  else
    mxy = evs(1);
    mny = evs(2);
    mxx = exs(1);
    mnx = exs(2);
  end    

  derivy(i,:) = [0,3*coefy(i,1),2*coefy(i,2),coefy(i,3)];
  extr        = fcardan(derivy(i,:));
  
  if extr(4) > 0 & extr(4) < Inf,
    for j = 1:extr(4)
      if extr(j) >= 0 & extr(j) <= bry(i+1)-bry(i)
        v = polyval(coefy(i,:),extr(j));
        x = polyval(coefx(i,:),extr(j));
        if v < mny, mny = v; mnx = x; end
        if v > mxy, mxy = v; mxx = x; end
      end
    end
  end
  miny(i) = mny;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [minx,maxx,derivx] = xextrema(lx,coefx,coefy,brx,bry)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1995-1999

derivx = zeros(lx,4);
minx   = zeros(1,lx);
maxx   = zeros(1,lx);

coefx = [coefx; 0,0,0,polyval(coefx(lx,:),brx(lx+1)-brx(lx))];
coefy = [coefy; 0,0,0,polyval(coefy(lx,:),brx(lx+1)-brx(lx))];
for i = 1:lx % = ly
  evs = coefx(i:i+1,4);
  eys = coefy(i:i+1,4);  
  if evs(1) < evs(2),
    mnx = evs(1);
    mxx = evs(2);
    mny = eys(1);
    mxy = eys(2);
  else
    mxx = evs(1);
    mnx = evs(2);
    mxy = eys(1);
    mny = eys(2);
  end    

  derivx(i,:) = [0,3*coefx(i,1),2*coefx(i,2),coefx(i,3)];
  extr        = fcardan(derivx(i,:));
  
  if extr(4) > 0 & extr(4) < Inf,
    for j = 1:extr(4)
      if extr(j) >= 0 & extr(j) <= brx(i+1)-brx(i)
        v = polyval(coefx(i,:),extr(j));
        y = polyval(coefy(i,:),extr(j));
        if v < mnx, mnx = v; mny = y; end
        if v > mxx, mxx = v; mxy = y; end
      end
    end
  end

  minx(i) = mnx;
  maxx(i) = mxx;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = fcardan(coef)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ALW, 1995-1999

global myeps 

if abs(coef(1)) / max(abs(coef(2:4))) > myeps
  coef = round(coef/myeps)*myeps;
  r = coef(2)/coef(1);
  s = coef(3)/coef(1);
  t = coef(4)/coef(1);
  
  p = s - r^2/3;
  q = 2*r^3/27-r*s/3+t;
  
  D = (p/3)^3 + (q/2)^2;
  u = (-q/2+sqrt(D))^(1/3);
  if (p==0) & (q==0),
    v = 0;
  else
    v = -p/(3*u);
  end

  e1= -0.5+ 0.86602540378444i;
  e2= -0.5- 0.86602540378444i;
  
  x1 =    u +    v - r/3; if abs(imag(x1)) < myeps, x1 = real(x1); end
  x2 = e1*u + e2*v - r/3; if abs(imag(x2)) < myeps, x2 = real(x2); end
  x3 = e2*u + e1*v - r/3; if abs(imag(x3)) < myeps, x3 = real(x3); end
  if sum(abs(imag([x1,x2,x3])))==0,
    out = [sort([x1,x2,x3]),3,3];
  elseif imag(x1) == 0,
    out = [x1,x2,x3,3,1];
  elseif imag(x2) == 0;
    out = [x2,x1,x3,3,1];
  else
    out = [x3,x1,x2,3,1];
  end
else 
  if abs(coef(2)) < myeps
    if abs(coef(3)) < myeps
      if abs(coef(4)) < myeps
        out = [NaN,NaN,NaN,Inf,Inf];
      else
        out = [NaN,NaN,NaN,0,0];
      end	
    else 
      out = [-coef(4)/coef(3),0,0,1,1];
    end;
  else
    D = sqrt(coef(3)^2-4*coef(2)*coef(4));
    x1 = (-coef(3)-D)/(2*coef(2)); if abs(imag(x1)) < myeps, x1 = real(x1); end
    x2 = (-coef(3)+D)/(2*coef(2)); if abs(imag(x2)) < myeps, x2 = real(x2); end
    if imag(D) == 0,
      out = [min(x1,x2),max(x1,x2),0,2,2];
    else
      out = [x1,x2,0,2,0];
    end
  end;
end;