function R2adm(action)

global c g p e t w a1 a2   %c curve
                           %g holds the geometry specification, 
                           %p(1,:) and p(2,:) hold the x- and y-coordinates of                             the nodal points,
                           %e holds edge type specification
                           %t(1:3,i) holds node numbers of triangle i
                           %w current mesh function
                           %a1 and a2 handles of the two axes
                           %c

%  a1=findobj('Position',[0.1 0.53 0.55 0.3]);
%  a2=findobj('Position',[0.1 0.13 0.55 0.3])
  a1=findobj('Userdata','a1');
  a2=findobj('Userdata','a2');
%  set(gcbf,'WindowButtonDownFcn','')
  rotate3d off  


switch action


case 'curve'

%  delete(findobj(gcbf,'Tag','rotaObj'))
  set(findobj(gcbf,'Tag','x='),'Visible','on')
  set(findobj(gcbf,'Tag','x'),'Visible','on')
  set(findobj(gcbf,'Tag','y='),'Visible','on')
  set(findobj(gcbf,'Tag','y'),'Visible','on')
  set(findobj(gcbf,'Tag','ok'),'Visible','on')
%  g=[];
  axes(a2)
  cla
  set(a2,'XLim',[-.1 1.1],'YLim',[-.1 1.1])
  xlabel('x')
  ylabel('y')
  c=[];
  line([0 0],[0 0],'Tag','candidate','Color','r')
  line([0 0],[0 0],'Tag','curve','Color','g')
  text(0,0,'','Tag','coords')
  set(gcbf,'WindowButtonMotionFcn','R2adm(''WBMF'')')
  set(gcbf,'WindowButtonDownFcn','R2adm(''WBDF'')')


case 'ok'

  x=get(findobj(gcbf,'Tag','x'),'String');
  y=get(findobj(gcbf,'Tag','y'),'String');
  if length(x)>0 & length(y)>0 & isempty(c)
    for i=1:41
      s=(i-1)/40;
      c(i,1)=eval(x);
      c(i,2)=eval(y);
    end
    set(a2,'XLimMode','auto','YLimMode','auto')
    set(findobj(gcbf,'Tag','curve'),...
       'XData',c(:,1),'YData',c(:,2),'Userdata',c)
  end
  R2adm('closeup')


case 'domain'

  g=[];
  if isempty(c)
    set(findobj(gcbf,'Tag','curveint'),...
              'String','define boudary curve','Visible','on')
    pause(2)
    set(findobj(gcbf,'Tag','curveint'),'Visible','off')
  else
    for i=1:length(c(:,1))-1
      g(2:3,i)=c(i:i+1,1);
      g(4:5,i)=c(i:i+1,2);
    end
    g(6,:)=ones(size(g(2,:)));
    g(7,:)=0*g(6,:);
    g(1,:)=2*g(6,:);
 end
% set(a2,'Userdata',g)
  [p,e,t]=initmesh(g);
  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  axes(a2)
  pdesurf(p,t,zeros(size(p(1,:)))')
  set(gca,'Userdata','a2')
  axis('equal')
  line(c(:,1),c(:,2),'Color','g','Tag','curve','Userdata',c)
  set(a2,'XLim',XLim,'YLim',YLim)
  view(2)
  h=findobj(gcbf,'Tag','rotaObj');
  delete(h)
XLim=get(a2,'XLim');
YLim=get(a2,'YLim');
if min(p(1,:))<XLim(1)|max(p(1,:))>XLim(2)|min(p(2,:))<YLim(1)|max(p(2,:))>YLim(2)
XLim=[min(p(1,:)) max(p(1,:))];
YLim=[min(p(2,:)) max(p(2,:))];
set(a2,'XLim',XLim,'YLim',YLim)
end

case 'WBMF'

  axes(a2)
  cp=get(a2,'CurrentPoint');
  x2=cp(1,1:2);
  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  dX=XLim(2)-XLim(1);
  dY=YLim(2)-YLim(1);
  if XLim(1)<x2(1) & x2(1)<XLim(2) & YLim(1)<x2(2) & x2(2)<YLim(2)
    if x2(1)-XLim(1)<0.02*dX
      XLim(1)=XLim(1)-.1*dX;
    end
    if XLim(2)-x2(1)<0.02*dX
      XLim(2)=XLim(2)+.1*dX; 
    end
    if x2(2)-YLim(1)<0.02*dY
      YLim(1)=YLim(1)-.1*dY;
    end
    if YLim(2)-x2(2)<0.02*dY
      YLim(2)=YLim(2)+.1*dY;
    end
    set(a2,'XLim',XLim,'YLim',YLim)
    xy=['(' num2str(x2(1)) ',' num2str(x2(2)) ')'];
    set(findobj(gcbf,'Tag','coords'),'String',xy,'Position',x2)
  end
  if isempty(c)
  else
    l=length(c(:,1));
    x1=c(l,:);
    set(findobj(gcbf,'Tag','candidate'),'XData',[x1(1) x2(1)],'YData',[x1(2) x2(2)])
  end


case 'WBDF'

  cp=get(a2,'CurrentPoint');
  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  if XLim(1) < cp(1,1) & cp(1,1) < XLim(2) & ...
     YLim(1) < cp(1,2) & cp(1,2) < YLim(2)
     if isempty(c)
       c=cp(1,1:2);
     else
       l=length(c(:,1));
       x0=c(1,:);
       x1=c(l,:);
       x2=cp(1,1:2);
       if norm(x2-x0)<.1
         x2=x0;
         R2adm('closeup')
       end
       c(l+1,:)=x2;
     end
     set(findobj(gcbf,'Tag','curve'),'XData',c(:,1),'YData',c(:,2))
  else
    disp('closing curve input')
    R2adm('closeup')
  end
    

case 'closeup'

  set(findobj(gcbf,'Tag','candidate'),'Visible','off')
  set(gcbf,'WindowButtonMotionFcn','')
  set(gcbf,'WindowButtonDownFcn','')
  set(findobj(gcbf,'Tag','coords'),'String','')
  set(findobj(gcbf,'Tag','x='),'Visible','off')
  set(findobj(gcbf,'Tag','x'),'Visible','off')
  set(findobj(gcbf,'Tag','y='),'Visible','off')
  set(findobj(gcbf,'Tag','y'),'Visible','off')
  set(findobj(gcbf,'Tag','ok'),'Visible','off')


case 'initiate'

  data=get(gcf,'Userdata');
  c=data{1};
  g=data{2};
  p=data{3};
  e=data{4};
  t=data{5};
  [x,y]=getxy(p);
  wxy='sin(12*x.*y)';
  w=eval(wxy);
%  R2adm('uv')
  

case 'close'

  data{1}=c;
  data{2}=g;
  data{3}=p;
  data{4}=e;
  data{5}=t;
%  set(gcf,'Userdata',data)


case 'modify'

  sp=get(a2,'CurrentPoint');
  index=0;
  dist=100;
  for j=1:length(c(:,1))
    newdist=norm(c(j,:)-sp(1,1:2));
    if newdist<dist
      index=j;
      dist=newdist;
    end
  end
  set(gcf,'Userdata',index)
  set(gcbf,'WindowButtonMotionFcn','R2adm(''move'')')
  set(gcbf,'WindowButtonUpFcn','R2adm(''stop'')')
  delete(findobj(gcbf,'Tag','curve'))
  line(c(:,1),c(:,2),'Color','g','Tag','curve','Userdata',c)


case 'move'

  index=get(gcf,'Userdata');
  cp=get(a2,'CurrentPoint');    
  x2=cp(1,1:2);
  xy=['(' num2str(x2(1)) ',' num2str(x2(2)) ')'];
  set(findobj(gcbf,'Tag','coords'),...
             'Position',x2,'String',xy,'Visible','on')
  l=length(c(:,1));
  if c(1,:)==c(l,:) & (index==1 | index==l)
    c(1,:)=cp(1,1:2);
    c(l,:)=c(1,:);
  else
    c(index,:)=cp(1,1:2);
  end
  set(findobj(gcbf,'Tag','curve'),'XData',c(:,1),'YData',c(:,2))


case 'stop'

  set(gcbf,'WindowButtonMotionFcn','')
  set(gcbf,'WindowButtonUpFcn','')
  set(findobj(gcbf,'Tag','coords'),'Visible','off')
  set(findobj(gcbf,'Tag','coords'),'String','')


case 'refine'

  axes(a2)
  [p,e,t]=refinemesh(g,p,e,t);
  pdemesh(p,e,t)
  set(gca,'Userdata','a2')
  axis('equal')
  line(c(:,1),c(:,2),'Color','g','Tag','curve')


case 'gallery'

  string=get(gcbo,'String');
  value=get(gcbo,'Value');
  w=get(gcbo,'Userdata');
  if value>1
    set(findobj(gcbf,'Tag',w),'String',string{value})
  end
  set(gcbo,'Value',1)


case 'add to gallery'

  wxy=get(gcbo,'String');
  str=get(findobj(gcbf,'Tag','gallery'),'String');
  l=length(str);  
  found='no ';
  for i=1:l
    stri=str{i};
    l1=length(wxy);
    l2=length(stri);
    if wxy(1:min(l1,l2))==stri(1:min(l1,l2))
      found='yes';
    end
  end
  if length(wxy)>6 & found=='no '
    str{l+1}=wxy;
    set(findobj(gcbf,'Tag','gallery'),'String',str)
    disp(['function ' wxy 'added to gallery'])
  end


case 'uv'

%if cbo==findobj(gcbf,'Tag','u')|cbo==findobj(gcbf,'Tag','v')
  [x,y]=getxy(p);
  wxy=get(gcbo,'String');
  w=eval(wxy);
  axes(a1)
  pdesurf(p,t,w')
  set(gca,'Userdata','a1')
 % axis('equal')
  R2adm('add to gallery')
  R2adm('add plotvalue')
  w_tag=get(gcbo,'Tag');
  title([w_tag '(x,y)'])
%end


case 'div'

  h=0.000001;
  [ux,uy,vx,vy]=getduv(p,h);
  div=ux+vy;
  axes(a1)
%  CP=get(a1,'CameraPosition');
%  CT=get(a1,'CameraTarget');
%  CVA=get(a1,'CameraViewAngle');
  pdesurf(p,t,div')
  set(gca,'Userdata','a1')
%  axis('equal')
%  set(a1,'CameraPosition',CP)
%  set(a1,'CameraTarget',CT)
%  set(a1,'CameraViewAngle',CVA)
  w=div;
  R2adm('add plotvalue')
  title('div(u,v)')


case 'rot'

  h=0.000001;
  [ux,uy,vx,vy]=getduv(p,h);
  rot=vx-uy;
  axes(a1)
  pdesurf(p,t,rot')
  set(gca,'Userdata','a1')
%  axis('equal')
  w=rot;
  R2adm('add plotvalue')
  title('rot(u,v)')


case 'div(grad(u))'

  uxy=get(findobj(gcbf,'Tag','u'),'String');
  [x,y]=getxy(p);
  h=0.0001;
  ux=diffq(x,y,uxy,h,0);
  uxx=(ux+diffq(x,y,uxy,-h,0))/h;
  uy=diffq(x,y,uxy,0,h);
  uyy=(uy+diffq(x,y,uxy,0,-h))/h;
  Lapl=uxx+uyy;
  axes(a1)
  pdesurf(p,t,Lapl')
  set(gca,'Userdata','a1')
%  axis('equal')
  w=Lapl;
  R2adm('add plotvalue')
  title('div(grad(u))')


case 'add plotvalue'

  axes(a1)
  ZLim=get(a1,'ZLim');
  hold on
  pdesurf(p,t,min(0,ZLim(1))*ones(size(w))')
  set(gca,'Userdata','a1')
%  axis('equal')
  hold off
  grid on


case 'integrate'

  axes(a2)
  hold on
  pdemesh(p,e,t)
  set(gca,'Userdata','a2')
  axis('equal')
  line(c(:,1),c(:,2),'Color','g','Tag','curve')
  hold off
  axes(a1)
  hold on
  pdemesh(p,e,t)
  set(gca,'Userdata','a1')
%  axis('equal')
  line(c(:,1),c(:,2),'Color','g','Tag','curve')
  hold off
  x=p(1,:);
  y=p(2,:);
  s1x=p(1,t(1,:))-p(1,t(3,:));
  s2x=p(1,t(2,:))-p(1,t(3,:));
  s1y=p(2,t(1,:))-p(2,t(3,:));
  s2y=p(2,t(2,:))-p(2,t(3,:));
  areas=abs(s1x.*s2y-s2x.*s1y)/2;
  w_means=(w(t(1,:))+w(t(2,:))+w(t(3,:)))/3;
  int=sum(w_means.*areas);
  set(findobj(gcbf,'Tag','integral'),...
       'String',['integral = ' num2str(int)],'Visible','on')
  pause(2)
  set(findobj(gcbf,'Tag','integral'),'Visible','off')


case 'tangent plane'

  [x,y]=getxy(p);
  uxy=get(findobj(gcbf,'Tag','u'),'String');
  u=eval(uxy);
  axes(a1)
  pdesurf(p,t,u')
  set(gca,'Userdata','a1')
%  axis('equal')
  R2adm('add plotvalue')
  title('u(x,y)')
  cp=get(a2,'CurrentPoint');
  x=cp(1,1);
  y=cp(1,2);
  u0=eval(uxy);
  h=0.000001;
  ux=diffq(x,y,uxy,h,0);
  uy=diffq(x,y,uxy,0,h);
  x0=x;
  y0=y;
  [x,y]=getxy(p);
  u=u0+ux*(x-x0)+uy*(y-y0);
  axes(a1)
  ZLim=get(a1,'ZLim');
  oldtps=findobj('Type','patch','Userdata','tp');
  delete(oldtps)
  hold on
  h=pdesurf(p,t,u');
  set(gca,'Userdata','a1')
%  axis('equal')
  set(h,'FaceColor','y','Userdata','tp')
  hold off
  set(a1,'ZLim',ZLim)
  set(gcf,'WindowButtonDownFcn','')
  set(gcf,'WindowButtonUpFcn','')

  
case 'mapped domain'

  if ~isempty(g)
    axes(a2)
    R2adm('domain')
    R2adm('mesh')
    title('xy domain')
    [x,y]=getxy(p);
    [uxy,vxy]=getuv(0);
    u=eval(uxy);
    v=eval(vxy);
    R2adm('Jac')
    axes(a1)
    pdesurf([u;v],t,(sign(w)-1)')
    set(gca,'Userdata','a1')
    axis('equal')
    view(2)    
    hold on
    pdemesh([u;v],e,t)
    set(gca,'Userdata','a1')
    axis('equal')
    hold off
    title('uv domain')
    set(findobj(gcbf,'Tag','contours'),'Value',0)
    set(findobj(gcbf,'Tag','gradients'),'Value',0)
 end
 set(findobj(gcbf,'Tag','jacobian'),...
    'String','jacobian','Callback','R2adm(''Jac'')')


case 'uv field'

  if isempty(g)
    XLim=get(a2,'XLim');
    dX=diff(XLim);
    YLim=get(a2,'YLim');
    dY=diff(YLim);
    x=[];
    y=[];
    z=(1:20)/21;
    for i=1:20
      x=[x XLim(1)+z*dX];
      y=[y YLim(1)+i*ones(size(z))/21*dY];
    end
  else
    x=(p(1,t(1,:))+p(1,t(2,:))+p(1,t(3,:)))/3;
    y=(p(2,t(1,:))+p(2,t(2,:))+p(2,t(3,:)))/3;
  end
  uxy=get(findobj(gcbf,'Tag','u'),'String');
  vxy=get(findobj(gcbf,'Tag','v'),'String');
  u=eval(uxy);
  v=eval(vxy);
  axes(a2)
  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  quiver(x,y,u,v)
  set(gca,'Userdata','a2')
  line(c(:,1),c(:,2),'Color','g','Tag','curve')
  set(a2,'XLim',XLim,'YLim',YLim)
  title('uv field')
  %line(c(:,1),c(:,2),'Color','g');
%  if gcbo==findobj(gcbf,'Tag','uv field')
%    pause(2)
%    R2adm('domain')
%    R2adm('mesh')
%    title('xy domain')
%    R2adm('mapped domain')
%  end


case 'Jac'
   
  h=0.000001;
  [ux,uy,vx,vy]=getduv(p,h);
  Jac=ux.*vy-uy.*vx;
  axes(a1)
  pdesurf(p,t,Jac')
  set(gca,'Userdata','a1')
%  axis('equal')
  w=Jac;
  R2adm('add plotvalue')
  title('d(u,v)/d(x,y)')
  set(findobj(gcbf,'Tag','jacobian'),...
     'String','uv map','Callback','R2adm(''mapped domain'')')


case 'contours'

  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  axes(a2)
  cla
  tit='';
  value=get(findobj(gcbf,'Tag','contours'),'Value');
  if value==1
     pdecont(p,t,w')
     set(gca,'Userdata','a2')
     axis('equal')
     line(c(:,1),c(:,2),'Color','g','Tag','curve')
    tit='levels';
  end
  value=get(findobj(gcbf,'Tag','gradients'),'Value');
  if value==1
    [wx,wy]=pdegrad(p,t,w');
    line(c(:,1),c(:,2),'Color','g','Tag','curve')
    x=(p(1,t(1,:))+p(1,t(2,:))+p(1,t(3,:)))/3;
    y=(p(2,t(1,:))+p(2,t(2,:))+p(2,t(3,:)))/3;
    hold on
    quiver(x,y,wx,wy)
    set(gca,'Userdata','a2')
%line(c(:,1),c(:,2),'Color','g','Tag','curve')
    hold off
    tit=[tit 'gradients'];
  end
  title(tit)
  set(a2,'XLim',XLim,'YLim',YLim)
  line(c(:,1),c(:,2),'Color','g','Tag','curve')


case 'gradients'

  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  axes(a2)
  cla
  tit='';
  value=get(findobj(gcbf,'Tag','contours'),'Value');
  if value==1
    pdecont(p,t,w')
    set(gca,'Userdata','a2')
    axis('equal')
    line(c(:,1),c(:,2),'Color','g','Tag','curve')
    tit='levels';
  end
  value=get(findobj(gcbf,'Tag','gradients'),'Value');
  if value==1
    [wx,wy]=pdegrad(p,t,w');
    x=(p(1,t(1,:))+p(1,t(2,:))+p(1,t(3,:)))/3;
    y=(p(2,t(1,:))+p(2,t(2,:))+p(2,t(3,:)))/3;
    hold on
    quiver(x,y,wx,wy)
    set(gca,'Userdata','a2')
    hold off
    tit=[tit 'gradients'];
  end
  title(tit)
  set(a2,'XLim',XLim,'YLim',YLim)
  line(c(:,1),c(:,2),'Color','g','Tag','curve')


case 'mesh'

  XLim=get(a2,'XLim');
  YLim=get(a2,'YLim');
  axes(a2)
  hold on
  pdemesh(p,e,t)
  set(gca,'Userdata','a2')
  axis('equal')
  line(c(:,1),c(:,2),'Color','g','Tag','curve')
  hold off
  set(a2,'XLim',XLim,'YLim',YLim)


case 'u ds'

size(c)
  dx=diff(c(:,1));
  dy=diff(c(:,2));
  ds=sqrt(dx.^2+dy.^2);
  a=ones(size(ds));
  b=zeros(size(ds));
  int=sumup(c(:,1),c(:,2),a,b,ds);

%  R2adm('vis bdry int')

  set(findobj(gcbf,'Tag','curveint'),...
        'String',['sum(uds) = ' num2str(int)],'Visible','on')
  pause(2)
  set(findobj(gcbf,'Tag','curveint'),'Visible','off')


case 'vis bdry int'

  axes(a1)
  uxy=get(findobj(gcbf,'Tag','u'),'String');
  X=c(:,1);
  Y=c(:,2);
  dx=diff(X);
  dy=diff(Y);
  ds=sqrt(dx.^2+dy.^2);
  ph=patch([0 0],[0 0],[0 0],'y');
  for i=1:length(X)-1
    for j=.2:.2:1
    x=X(i)+(j-.2)*dx(i);
    y=Y(i)+(j-.2)*dy(i);
    u1=eval(uxy);
    x1=x;
    y1=y;
    x=X(i)+j*dx(i);
    y=Y(i)+j*dy(i);
    u2=eval(uxy);
    x2=x;
    y2=y;
    set(ph,'XData',[x1 x2 x2 x1],...
           'YData',[y1 y2 y2 y1],...
           'ZData',[0 0 u2 u1])
    drawnow
    end
  end
  delete(ph)


case 'work'

  dx=diff(c(:,1));
  dy=diff(c(:,2));
  ds=sqrt(dx.^2+dy.^2);
  a=dx./ds;
  b=dy./ds;
  int=sumup(c(:,1),c(:,2),a,b,ds);
  R2adm('uv field')
  set(findobj(gcbf,'Tag','curveint'),...
        'String',['sum(udx+vdy) = ' num2str(int)],'Visible','on')
  pause(4)
  set(findobj(gcbf,'Tag','curveint'),'Visible','off')



case 'crossflow'

  dx=diff(c(:,1));
  dy=diff(c(:,2));
  ds=sqrt(dx.^2+dy.^2);
  a=dy./ds;
  b=-dx./ds;
  int=sumup(c(:,1),c(:,2),a,b,ds);
  R2adm('uv field')
  set(findobj(gcbf,'Tag','curveint'),...
        'String',['sum(udy-vdx) = ' num2str(int)],'Visible','on')
  pause(4)
  set(findobj(gcbf,'Tag','curveint'),'Visible','off')



end




function dw=diffq(x,y,wxy,h,k)

  w=eval(wxy);
  x=x+h;
  y=y+k;
  dw=(eval(wxy)-w)/sqrt(h^2+k^2);


function int=sumup(X,Y,a,b,ds)

  uxy=get(findobj(gcbf,'Tag','u'),'String');
  vxy=get(findobj(gcbf,'Tag','v'),'String');
  dx=diff(X);
  dy=diff(Y);
  int=0;
  z=(0:10)/10;
  for i=1:length(ds)
    x=X(i)+dx(i)*z;
    y=Y(i)+dy(i)*z;
    u=eval(uxy);
    v=eval(vxy);
    u=(sum(u(1:10))+sum(u(2:11)))/20;
    v=(sum(v(1:10))+sum(v(2:11)))/20;
    int=int+(u*a(i)+v*b(i))*ds(i);
%    patch([X(i) X(i)+dx(i) X(i)+dx(i) X(i)],[Y(i) Y(i)+dy(i) Y(i)+dy(i) Y(i)],[0 0 1 1])
  end


function [ux,uy,vx,vy]=getduv(p,h)

  [uxy,vxy]=getuv(0);
  [x,y]=getxy(p);
  ux=diffq(x,y,uxy,h,0);
  vx=diffq(x,y,vxy,h,0);
  uy=diffq(x,y,uxy,0,h);
  vy=diffq(x,y,vxy,0,h);


function [uxy,vxy]=getuv(dummy)

  uxy=get(findobj(gcbf,'Tag','u'),'String');
  vxy=get(findobj(gcbf,'Tag','v'),'String');


function [x,y]=getxy(p)

  x=p(1,:);
  y=p(2,:);
