function PDEPlotCut(action,flag)

persistent Npts;
persistent plotcut_fig plotcut_orient plotcut_draw;
persistent plotcut_xy plotcut_xy_label plotcut_yxmax plotcut_yxmax_label;
persistent plotcut_yxmin plotcut_yxmin_label;

if(nargin<1)
  action='initialize';
  Npts=129;
end

if(strcmp(lower(action),'initialize'))
  pde_ax=findobj(allchild(0),'Tag','PDEAxes');
  limits=get(pde_ax,'Xlim');


  plotcut_fig=findobj(allchild(0),'Tag','plotcut_fig');
  if(length(plotcut_fig))
    figure(plotcut_fig);
    position=get(plotcut_fig,'Position');
  else
    plotcut_fig=figure;
    position=get(plotcut_fig,'Position');
    position=[position(1) position(2)+(1-.5)*position(4) ... 
              .3*position(3) .5*position(4)];
  end

  set(plotcut_fig...
      ...%,'Visible','off'...
      ,'Tag','plotcut_fig'...
      ,'NumberTitle','off'...
      ,'Name','CrossSection'...
      ,'resize','on'...
      ,'MenuBar','none'...
      ,'Position',position...  
    );
 
  plotcut_orient=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_orient'...
	    ,'Style','popupmenu'...
	    ,'Units','normalized'...
	    ,'Position',[0 .8 1 .2]...
	    ,'String','Horizontal|Vertical'...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    ,'CallBack','PDEPlotCut(''plotcut_orient'');'...
	    );

  plotcut_xy=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_xy'...
	    ,'Style','edit'...
	    ,'Units','normalized'...
	    ,'Position',[.5 .6 .5 .2]...
	    ,'BackGroundColor',[1 1 1]...
	    ,'String',num2str(0) ...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    ,'CallBack','PDEPlotCut(''bound'',1)'...
	    );

  plotcut_xy_label=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_xy_label'...
	    ,'Style','text'...
	    ,'Units','normalized'...
	    ,'Position',[0 .6 .5 .2]...
	    ,'String','y ='...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    );

  plotcut_yxmax=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_yxmax'...
	    ,'Style','edit'...
	    ,'Units','normalized'...
	    ,'Position',[.5 .4 .5 .2]...
	    ,'BackGroundColor',[1 1 1]...
	    ,'String',num2str(limits(2)) ...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    ,'CallBack','PDEPlotCut(''bound'',2)'...
	    );

  plotcut_yxmax_label=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_yxmax_label'...
	    ,'Style','text'...
	    ,'Units','normalized'...
	    ,'Position',[0 .4 .5 .2]...
	    ,'String','xmax ='...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    );

  plotcut_yxmin=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_yxmin'...
	    ,'Style','edit'...
	    ,'Units','normalized'...
	    ,'Position',[.5 .2 .5 .2]...
	    ,'BackGroundColor',[1 1 1]...
	    ,'String',num2str(limits(1)) ...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    ,'CallBack','PDEPlotCut(''bound'',3)'...
	    );

  plotcut_yxmin_label=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_yxmin_label'...
	    ,'Style','text'...
	    ,'Units','normalized'...
	    ,'Position',[0 .2 .5 .2]...
	    ,'String','xmin ='...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    );

  plotcut_draw=...
  uicontrol(plotcut_fig...
	    ,'Tag','plotcut_draw'...
	    ,'Style','Push'...
	    ,'Units','normalized'...
	    ,'Position',[0 0 1 .2]...
	    ,'String','Draw'...
	    ,'FontUnits','normalized'...
	    ,'FontSize',.5 ...
	    ,'CallBack','PDEPlotCut(''plotcut_draw'');'...
	    );

elseif(strcmp(lower(action),'plotcut_orient'));
  pde_ax=findobj(allchild(0),'Tag','PDEAxes');
  HorizOrVert=get(plotcut_orient,'Value');
  switch(HorizOrVert)
    case 1 %horizontal cut coordinate is y
      limits=get(pde_ax,'Xlim');

      set(plotcut_xy_label,'String','y =');
      set(plotcut_xy,'String',num2str(0.5*sum(limits)));

      set(plotcut_yxmax_label,'String','xmax=');
      set(plotcut_yxmax,'String',num2str(limits(2)));

      set(plotcut_yxmin_label,'String','xmin=');
      set(plotcut_yxmin,'String',num2str(limits(1)));

    case 2 %vertical cut coordinate is x
      limits=get(pde_ax,'Ylim');

      set(plotcut_xy_label,'String','x =');
      set(plotcut_xy,'String',num2str(0.5*sum(limits)));

      set(plotcut_yxmax_label,'String','ymax=');
      set(plotcut_yxmax,'String',num2str(limits(2)));

      set(plotcut_yxmin_label,'String','ymin=');
      set(plotcut_yxmin,'String',num2str(limits(1)));
  end

elseif(strcmp(lower(action),'bound'));
  pde_ax=findobj(allchild(0),'Tag','PDEAxes');

  HorizOrVert=get(plotcut_orient,'Value');
  switch(HorizOrVert)
  case 1 %horizontal cut
    limits=get(pde_ax,'Xlim');
  case 2 %vertical cut
    limits=get(pde_ax,'YLim');
  end

  switch(flag) %which edit field was just used

  case 1 %plotcut_xy
    xy=str2num(get(plotcut_xy,'String'));
    if(length(xy)==0); xy=0; end
    switch(HorizOrVert)
    case 1 %horizontal cut
      xy=min(max(xy,limits(1)),limits(2));
    case 2 %vertical cut
      xy=min(max(xy,limits(1)),limits(2));
    end          
    set(plotcut_xy,'String',num2str(xy));
         
  case 2 %plotcut_yxmax
    yxmax=str2num(get(plotcut_yxmax,'String'));
    if(length(yxmax)==0); yxmax=0; end
    switch(HorizOrVert)
    case 1 %horizontal cut
      yxmax=min(max(yxmax,limits(1)),limits(2));
    case 2 %vertical cut
      yxmax=min(max(yxmax,limits(1)),limits(2));
    end          
    set(plotcut_yxmax,'String',num2str(yxmax));

  case 3 %plotcut_yxmin
    yxmin=str2num(get(plotcut_yxmin,'String'));
    if(length(yxmin)==0); yxmin=0; end
    switch(HorizOrVert)
    case 1 %horizontal cut
      yxmin=min(max(yxmin,limits(1)),limits(2));
    case 2 %vertical cut
      yxmin=min(max(yxmin,limits(1)),limits(2));
    end
    set(plotcut_yxmin,'String',num2str(yxmin));
  end 

elseif(strcmp(lower(action),'plotcut_draw'))
  %just in case user changed bounds since cut selected
  PDEPlotCut('bound',1);
  PDEPlotCut('bound',2);
  PDEPlotCut('bound',3);

  %find the current axes... old one might have died
  pde_ax=findobj(allchild(0),'Tag','PDEAxes');
  h=findobj(get(pde_ax,'Children'),'Type','patch');
  for i=1:length(h)
    if(length(get(h(i),'ZData')))
      pde_ax_points=h(i);
      break
    end
  end

  HorizOrVert=get(plotcut_orient,'Value');
  %get interpolation points from cut line
  switch(HorizOrVert)
  case 1 %horizontal cut
    yi=str2num(get(plotcut_xy,'String'));
    xmin=str2num(get(plotcut_yxmin,'String'));
    xmax=str2num(get(plotcut_yxmax,'String'));
    if(xmax<xmin)
      temp=xmax;
      xmax=xmin;
      xmin=temp;
      set(plotcut_yxmin,'String',num2str(xmin));
      set(plotcut_yxmax,'String',num2str(xmax));
    end
  case 2 %vertical cut
    xi=str2num(get(plotcut_xy,'String'));
    ymin=str2num(get(plotcut_yxmin,'String'));
    ymax=str2num(get(plotcut_yxmax,'String'));
    if(ymax<ymin)
      temp=ymax;
      ymax=ymin;
      ymin=temp;
      set(plotcut_yxmin,'String',num2str(ymin));
      set(plotcut_yxmax,'String',num2str(ymax));
    end
  end


  %get the x, y, z, and triangulation data from the picture
  x=get(pde_ax_points,'XData');  x=x(:);
  y=get(pde_ax_points,'YData');  y=y(:);
  z=get(pde_ax_points,'CData');  z=z(:);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %reduce number of triangles for tsearch to search to  %
  %only those that the cross section line passes through%
  %also eliminate edge patches                          %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  tri0=get(pde_ax_points,'faces');

  switch(HorizOrVert)
  case 1 %horizontal cut
    yminmax=minmax(y(tri0));
    itri=find((yminmax(:,1)<=yi)&(yi<=yminmax(:,2))...
	      &(yminmax(:,1)~=yminmax(:,2)));
  case 2 %vertical cut
    xminmax=minmax(x(tri0));
    itri=find((xminmax(:,1)<=xi)&(xi<=xminmax(:,2))...
	      &(xminmax(:,1)~=xminmax(:,2)));
  end    

  tri1=tri0(itri,:);
  tri2=zeros(size(tri1));

  ioldpointskept=unique(tri1);

  x=x(ioldpointskept);
  y=y(ioldpointskept);
  z=z(ioldpointskept);

  for i=1:length(ioldpointskept)
    tri2(find(tri1==ioldpointskept(i)))=i;
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %DANGER picture contains duplicate %
  %points and triangles.  This screws% 
  %up tsearch so have to fix it...   %
  %         FIX STARTS HERE          %
  %even though tsearch is no longer  %
  %used I've decided to keep this    %
  %just in case                      %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  xyz=[x y z];
  [XYZ,iunique,junique]=unique(xyz(:,1:2),'rows');
  XYZ=xyz(iunique,:);

  tri3=tri2;    

  for i=1:length(iunique)
    j=find((x==XYZ(i,1))&(y==XYZ(i,2)));
    tri3(find(ismember(tri2,j)))=i;
  end
    
  tri4=unique(sort(tri3')','rows');
  x=XYZ(:,1); y=XYZ(:,2); z=XYZ(:,3);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %          FIX ENDS HERE           %
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  pdetitlestring=get(get(pde_ax,'Title'),'String');
  icolor=findstr(lower(pdetitlestring),'color: ');
  if(length(icolor))
    pdetitlestring=pdetitlestring((icolor+7):length(pdetitlestring));
  end

  figure;
  switch(HorizOrVert)
  case 1 %Horizontal
    [xi,zi]=myinterp(x,y,z,tri4,yi);
    ax=[max(xi(1),xmin) min(xi(length(xi)),xmax) min(zi) max(zi)];
    axis(ax); hold on;
    hax=plot(xi,zi,'-');
    title([pdetitlestring sprintf(': cut at y=%g',yi)]);
    xlabel('x');
  case 2 %vertical
    [yi,zi]=myinterp(y,x,z,tri4,xi);
    ax=[max(yi(1),ymin) min(yi(length(yi)),ymax) min(zi) max(zi)];
    axis(ax); hold on;
    hax=plot(yi,zi,'-');
    title([pdetitlestring sprintf(': cut at x=%g',xi)]);
    xlabel('y');
  end
  hold off;
  grid on;
end

return

function b=minmax(a)
b=[min(a')' max(a')'];

return

function [xi,zi]=myinterp(x,y,z,tri,yi)
  %assumes a horizontal cut if veritcal switch x and y
  N=length(tri);
  [ysort,I]=sort(y(tri)');
  xtemp=x(tri)';
  ztemp=z(tri)';
  xsort=zeros([3 N]);
  zsort=xsort;
  for i=1:N
    xsort(:,i)=xtemp(I(:,i),i);
     zsort(:,i)=ztemp(I(:,i),i);
  end
  xsort=xsort';
  ysort=ysort';
  zsort=zsort';

  %the line allways passes through the triangle edge connecting
  %the points with the triangle's maximum and minimum y 
  yada=(yi-ysort(:,1))./(ysort(:,3)-ysort(:,1));
  x1=yada.*(xsort(:,3)-xsort(:,1))+xsort(:,1);
  z1=yada.*(zsort(:,3)-zsort(:,1))+zsort(:,1);

  %part of time, the line passes through the triangle edge
  %connecting the points with the triangle's maximum and middle
  %value
  i=find(yi>=ysort(:,2));
  x2=zeros([N 1]);
  z2=x2;
  yada=(yi-ysort(i,2))./(ysort(i,3)-ysort(i,2));
  x2(i)=yada.*(xsort(i,3)-xsort(i,2))+xsort(i,2);
  z2(i)=yada.*(zsort(i,3)-zsort(i,2))+zsort(i,2);

  %part of time, the line passes through the triangle edge
  %connecting the points with the triangle's middle and minimum
  %value
  i=find(yi<ysort(:,2));
  yada=(yi-ysort(i,1))./(ysort(i,2)-ysort(i,1));
  x2(i)=yada.*(xsort(i,2)-xsort(i,1))+xsort(i,1);
  z2(i)=yada.*(zsort(i,2)-zsort(i,1))+zsort(i,1);

  %reorder so has history of line intersections with triangle
  %edges
  [xtemp,i]=sort([x1 x2],2);
  ztemp=[z1 z2];
  for j=1:length(z1)
    z12(j,:)=ztemp(j,i(j,:));
  end

  [xz,i,j]=unique([xtemp z12],'rows');
  x12=xz(:,1:2);
  z12=xz(:,3:4);


  NN=size(x12,1);  
  xi=zeros([NN+1 1]);
  zi=xi;
  xi(1)=x12(1,1);
  zi(1)=z12(1,1);
  j=1;
  for i=1:NN
    %remove duplicate points
    if(x12(i,1)~=x12(i,2))
      j=j+1;
      xi(j)=x12(i,2);
      zi(j)=z12(i,2);
    end

    %insert NaN's in holes between triangles
    if(i<NN)
      if(x12(i,2)~=x12(i+1,1))
        j=j+1;
        xi(j)=0.5*(x12(i,2)+x12(i+1,1));
        zi(j)=NaN;
        j=j+1;
        xi(j)=x12(i+1,1);
        zi(j)=z12(i+1,1);
      end
    end
  end

  xi=xi(1:j);
  zi=zi(1:j);

return