#Knihovna mvcalp kveten 2001
#autor Roman Plch
#upraveno z verze pro Maple V R3 na verzi pro Maple 7
#Pokud je text napovedy v anglictine, je procedura prejata z knihovny
#mvcal (http://fmwww.bc.edu/MT/gopher/maple/mvcal/dir.html),
#pokud je text napovedy v cestine, je procedura nove naprogramovana.
#Veskery text napovedy z puvodni verze byl prepsan
#do Maple Worksheetu verze Maple V R5.
#Nove napsana byla napoveda pro proceduru mvcalp[takereal].


mvcalp:=module()

  export difer, difern, drawpoint, drawplane, dzd1, dzd2, dzdd1, dzdd2, dzdd12,
         GraphTan, ParamTan, ImplicitTan, graf_t,
         implicitdiff, implicitdiffb, kmen, lagrange, levelcurve, levelsurface,
         mvextrem, sing, takereal;

  option package;

takereal:= proc(list)
  local i,j,result,real;
    result := {};
      for i to nops(list) do
        real := true;
        for j to nops(op(i,list)) do
          if not type(op(2,op(j,op(i,list))),realcons) then real := false fi
        od;
        if real = true then result := result union {op(i,list)} fi
      od;
  RETURN(result)
end:

sing:= proc(f)
  local cp;
  cp:={solve({diff(f,x)=0, diff(f,y)=0}, {x,y})};
  RETURN(cp)
end:

mvextrem:= proc(f)
  local zxx, zyy, zxy, D, i, p2, pom;
  zxx:=diff(f,x,x);
  zyy:=diff(f,y,y);
  zxy:=diff(f,x,y);
  pom:=map(allvalues, sing(f));
  pom:=takereal(pom);
  for i from 1 to nops(pom) do
    p2:=op(i,pom);
    D:=evalf(subs(p2,zxx)*subs(p2,zyy)-subs(p2,zxy)^2);
    if D=0 then print( p2, ` nelze rozhodnout`);
    elif D<0 then print( p2, ` extrem nenastava`);
      else
        if evalf(subs(p2,zxx)) > 0 then print( p2, ` lokalni minimum`);
        else print( p2, ` lokalni maximun`);
        fi;
    fi;
  od;
end:

GraphTan:= proc(f,xrange,yrange,pt)

  local xmin,xmax,ymin,ymax,x0,y0,z0,dx,dy,xsour,ysour,
        tanfunc,gpha,gphb,tanpt,optio,rovnice;


  x0 := op(2,pt)[1]:
  y0 := op(2,pt)[2]:
  xmin := op(1,op(2,xrange)):
  xmax := op(2,op(2,xrange)):
  ymin := op(1,op(2,yrange)):
  ymax := op(2,op(2,yrange)):
  optio:=args[5..nargs];

  dx := subs(x=x0,y=y0,diff(f,x)):
  dy := subs(x=x0,y=y0,diff(f,y)):

  z0 := subs(x=x0,y=y0,f):

  tanfunc:=z0+dx*(x-x0)+dy*(y-y0):

  tanpt:=plots[pointplot3d]({[x0,y0,z0]},color=red):
  gpha := plot3d(f,xrange,yrange,optio):
  xsour:=abs(xmax-xmin)/4;
  ysour:=abs(ymax-ymin)/4;
  gphb:=plot3d(tanfunc,x=x0-xsour..x0+xsour,y=y0-ysour..y0+ysour,optio):

  rovnice:=z=tanfunc:
  print(`Tecna rovina ma rovnici `);
  print(rovnice);

  plots[display]([gpha, gphb, tanpt]);
end:

difer:=proc()

  local derx,dery,dif;

  if nargs=1 then
    #print('Diferencial v obecnem bode ');
    print(diff(args[1],x)*h+diff(args[1],y)*k);
  fi;
  if nargs=5 then
    #print('Diferencial funkce v danem bode a s danymi diferencemi je');
    derx:=subs(x=args[2],y=args[3],diff(args[1],x));
    dery:=subs(x=args[2],y=args[3],diff(args[1],y));
    dif:=derx*args[4]+dery*args[5];
    RETURN(dif);
  fi;
  if nargs<>1 and nargs <>5 then print ('Spatne_zadano') fi;
end:

difern:=proc(funkce,m)
  local j;
  RETURN(sum(binomial(m,j)*diff(funkce, x$j,y$(m-j))*h^j*k^(m-j), j=0..m));
end:

ParamTan:= proc (f,srange,trange,pt)
  #Define local variables
  local
  planefunc,s0,t0,smin,smax,tmin,tmax,pama,dfs,dft,z0,pamb,pamc,tanpt,optio;


  #Extract specific variables we need from the larger variables define as input from user
  s0:= op (2,pt)[1]:
  t0:= op (2,pt)[2]:
  smin := op(1,op(2,srange)):
  smax := op(2,op(2,srange)):
  tmin := op(1,op(2,trange)):
  tmax := op(2,op(2,trange)):
  optio:= args[5..nargs];

  #Plot the parametric surface.  Note: f is a 3-vector as defined by user
  pama:= plot3d(f,srange,trange,style=HIDDEN,axes=BOXED):

  #Find the tanget vectors by taking partials with respect to s and then to t
  dfs := subs(s=s0,t=t0,diff(f,s)):
  dft := subs(s=s0,t=t0,diff(f,t)):

  #Find the point of tangency
  z0  := subs (s=s0,t=t0,f);

  #Plot point of tangency under variable name tanpt
  tanpt := plots[pointplot3d] ({z0},color=black):

  #Parameterize plane using tangent vectors and k,l as parameters, add z0 to give proper translation
  planefunc := evalm (k*dfs+l*dft+z0):

  #Plot the tangent plane under variable name pamb, limit parameters so as not to obstruct viewing
  pamb:=plot3d(planefunc,k=-(smax-smin)/8..(smax-smin)/8,l=-(tmax-tmin)/8..(tmax-tmin)/8,style=PATCHNOGRID):

  #Display the plots
  plots[display] ([pama,pamb,tanpt],optio);
end:

ImplicitTan:= proc (F,xrange,yrange,zrange,pt)
  #Define local variable names
  local f,K,gradat,tanpt,x0,y0,z0,impa,impb,dx,dy,dz,vec0,planefunc,optio;



  #Extract different variables we need from larger variable given by user input
  x0 := op (2,pt)[1]:
  y0 := op (2,pt)[2]:
  z0 := op (2,pt)[3]:
  optio:=args[6..nargs];

  #Plot implicit surface under variable name impa
  impa := plots[implicitplot3d](F,xrange,yrange,zrange,style=HIDDEN,axes=BOXED):

  #If F is just an expression rename it f
  f := F:

  #If F is an equation, extract the expression as f and the constant as K
  if type(F,`=`) = 'true' then
    f:= lhs (F):
    K := rhs (F):
  fi;

  #Evaluate partials of f with repect to x, y, and z.
  dx := subs(x=x0,y=y0,z=z0,diff(f,x)):
  dy := subs(x=x0,y=y0,z=z0,diff(f,y)):
  dz := subs(x=x0,y=y0,z=z0,diff(f,z)):

  #Create the gradient vector at [x0,y0,z0] Which we know is normal to suface
  gradat := [dx,dy,dz]:

  #Find the fuction for tangent plane (implicitly defined)
  #Using fact that plane is given by--  grad(x0,y0,z0) dot ([x-x0,y-y0,z-z0]) = 0
  vec0:=  (evalm([x,y,z]-[x0,y0,z0])):
  planefunc:= innerprod(gradat,vec0):

  #Plot tangent plane under variable name impb
  impb := implicitplot3d(planefunc,xrange,yrange,zrange):

  #Plot point of tangency
  tanpt := pointplot3d({[x0,y0,z0]},color=red):

  #Display surface, tangent plane, tangent point
  display ([impa,impb,tanpt],optio);
end:

kmen:= proc(m,n)
  if diff(m,y)=diff(n,x) then
    simplify(integrate(m,x)+integrate(n-diff(integrate(m,x),y),y));
  else print(`Zadany vyraz neni diferencialem zadne funkce`);
  fi
end:

dzd1:= proc(z,uu,vv,u,v,x,y)
  diff(z,u)*diff(uu,x)+diff(z,v)*diff(vv,x);
end:

dzd2:= proc(z,uu,vv,u,v,x,y)
  diff(z,u)*diff(uu,y)+diff(z,v)*diff(vv,y);
end:

dzdd1:= proc(z,uu,vv,u,v,x,y)
  diff(z,u,u)*diff(uu,x)^2+2*diff(z,u,v)*diff(vv,x)*diff(uu,x)+
  diff(z,v,v)*diff(vv,x)^2+diff(z,u)*diff(uu,x,x)+
  diff(z,v)*diff(vv,x,x);
end:

dzdd2:= proc(z,uu,vv,u,v,x,y)
  diff(z,u,u)*diff(uu,y)^2+2*diff(z,u,v)*diff(vv,y)*diff(uu,y)+
  diff(z,v,v)*diff(vv,y)^2+diff(z,u)*diff(uu,y,y)+
  diff(z,v)*diff(vv,y,y);
end:

dzdd12:=proc(z,uu,vv,u,v,x,y)
  diff(z,u,u)*diff(uu,x)*diff(uu,y)+diff(z,u,v)*diff(vv,y)*diff(uu,x)+
  diff(z,u,v)*diff(vv,x)*diff(uu,y)+diff(z,v,v)*diff(vv,x)*diff(vv,y)+
  diff(z,u)*diff(uu,x,y)+diff(z,v)*diff(vv,x,y);
end:

implicitdiff:= proc(g)
  local tmp,DIFFg,DIFFy,DIFFy0,p1:
  DIFFg:= diff(g,x):
  DIFFy:=simplify(solve(subs(diff(y(x),x)=p1,DIFFg)=0,p1));
end:

implicitdiffb:= proc(x0,y0,g)
  local tmp,DIFFg,DIFFy,DIFFy0,p1:
  tmp:=subs(y(x)=y0,g):
  if (simplify(subs(x=x0,tmp)) <> 0) then
    ERROR(` x0,y0 appear not to be on the curve`):
  fi:
  DIFFg:= diff(g,x):
  DIFFy:=simplify(solve(subs(diff(y(x),x)=p1,DIFFg)=0,p1)):
  DIFFy0:= simplify(subs(x=x0,y(x0)=y0,DIFFy)):
  DIFFy0
end:

graf_t:= proc() local a,b,c,u,v,k;
  a:=diff(args[1],x);
  b:=diff(args[1],y);
  u:=op(1,args[2]);
  v:=op(2,args[2]);
  c:=eval(subs({x=u,y=v},args[1]));
  if c=0 then
    k:=(subs({x=u,y=v},a)*(x-u)+subs({x=u,y=v},b)*(y-v));
    print(`Rovnice tecny v bode`,args[2],`je `,k=0);
    if nargs(graf_t)=6 then
      RETURN (plots[implicitplot]({args[1],k},x=args[3]..args[4],y=args[5]..args[6]));
    fi;
  fi;
  if c<>0 then
    print(`Bod`,args[2],`nelezi na krivce `,args[1]=0);
  fi;
end:

drawpoint:= proc( point)
  local u1 , u2, u3, u4;
  u1 := point[1];
  u2 := point[2];
  u3 := point[3];
  u4 := max( abs(u1),abs(u2),abs(u3));
  plot3d( {[u1+t, u2+t, u3], [u1, u2+t, u3+t],[u1+t, u2, u3+t], [u1, u2,s*u3],
  [s*u1, u2, 0], [u1, s*u2,0]}, t=-0.1*u4-0.01..0.1*u4+0.01, s=0..1,
  grid=[2,2], view=[-u4-1..u4+1,
  -u4-1..u4+1, -u4-1..u4+1], axes= NORMAL);
end:

lagrange:= proc(f,c,g,h,v)
  local curve;
  curve := solve(f=c, y);
  plot({curve, solve(g,y)},h,v);
end:

drawplane:= proc()
  local a ,b , c, d, temp ;
  if (nops(args[1]) <> 2) then
    ERROR(`input equation is not a plane` )
  fi;
  temp := op(1, args[1])- op(2, args[1]):
  a:= diff(temp,x) : b:=diff(temp,y):  c:= diff(temp, z): d:=subs(x=0,y=0,z=0, temp):
  if not type (evalf(a), numeric) or not type (evalf(b), numeric) or not type
              (evalf(c), numeric) or not type (evalf(d), numeric)  then
    ERROR(`input equation is not a plane` )
  fi;
  if (op(1, args[2]) = 'x') or (op(1, args[3]) = 'x') then
    if (op(1, args[2]) = 'y') or (op(1, args[3]) = 'y')  then
      if (evalf(c) =0) then
        ERROR(`the given expression cannot draw a plane`)
      else
        plot3d( (-d-a*x-b*y)/c, args[2..nargs], grid=[5,5], axes=framed);
      fi;
    else
      if (evalf(b) =0) then
        ERROR(`the given expression cannot draw a plane`)
      else
       plot3d( [x, (-d-a*x-c*z)/b, z], args[2..nargs], grid=[5,5], axes=framed);
      fi;
    fi;
  else
    if (evalf(a) =0) then
      ERROR(`the given expression cannot draw a plane`)
    else
      plot3d([(-d-b*y-c*z)/a, y, z], args[2..nargs], grid=[5,5], axes=framed);
    fi;
  fi;
end:

levelcurve:= proc()
  local a1, a2, a3, a4, b, b1, b2, b3, b4, c1, e, f, F, c,
  opts, x, y, xrange,
  yrange, curve, r, i, ans, j;
  F:=args[1]:
  c:=args[2]:
  x:=lhs(args[3]):
  y:=lhs(args[4]):
  xrange:=args[3]:
  yrange:=args[4]:
  opts:=args[5..nargs]:
  if type(c, constant) then
    c1 := [c]
  else c1 := c
  fi;
  if not type(c1, list(constant)) then
    ERROR(`illegal constant`)
  fi;
  curve := {};
  b3 := diff(F,y$3);
  a3 := diff(F,x$3);
  if b3 =0 then
    if a3 =0 then
      for i from 1 to nops(c1) do
        ans := {solve(F=c1[i],x)};
        for j from 1 to nops(ans) do
          curve := curve union { [subs(y=t,ans[j]), t, t =op(2,yrange)]};
        od;
      od;
      for i from 1 to nops(c1) do
        ans := {solve(F=c1[i],y)};
        for j from 1 to nops(ans) do
          curve := curve union { [t, subs(x=t, ans[j]), t=op(2,xrange)]};
        od;
      od;
      if nargs<5 then plot(curve, xrange, yrange);
      else plot(curve, xrange, yrange, opts);
      fi;
    else
      for i from 1 to nops(c1) do
        curve := curve union {solve(F=c1[i],y)};
      od;
      if nargs<5 then plot(curve, xrange, yrange);
      else plot(curve, xrange, yrange, opts);
      fi;
    fi;
  else
    if a3 =0 then
      for i from 1 to nops(c1) do
        ans := {solve(F=c1[i],x)};
        for j from 1 to nops(ans) do
          curve := curve union { [ans[j], y, yrange]};
        od;
      od;
    if nargs < 5 then plot(curve, xrange, yrange);
    else plot(curve, xrange, yrange, opts);
    fi;
    else
      if diff( subs(x^4=r, subs(x^2=r, f)),x) = 0 then
        for i from 1 to nops(c1) do
          ans := {solve(F=c1[i],x)};
          for j from 1 to nops(ans) do
            curve := curve union { [ans[j], y,yrange]};
          od;
        od;
        if nargs<5 then plot(curve, xrange, yrange);
        else plot(curve, xrange, yrange, opts);
        fi;
        else
          for i from 1 to nops(c1) do
            curve := curve union {solve(F=c1[i],y)};
          od;
        if nargs<5 then plot(curve, xrange, yrange);
        else plot(curve, xrange, yrange,opts);
        fi;
      fi;
    fi;
  fi;
end:

levelsurface:= proc(f, c, h, v)
  local surface, i, j, ans;
  if nops(c) = 1 then
    if not (diff(f, z) = 0) then
      surface := solve(f = c, z);
      plot3d({surface}, h, v, grid=[12,12]);
    else
      surface :={ } ;
      ans := {solve(f=c, y) };
      for j from 1 to nops(ans) do
        surface := surface union{[x, ans[j], z]};
      od;
      plot3d(surface, h, z=-2..2, grid=[12,12]);
    fi;
  else
    surface := {};
    if not (diff(f, z) = 0) then
      for i from 1 to nops(c) do
        surface := surface union {solve(f = c[i], z)};
      od;
      plot3d(surface, h, v, grid=[12,12]);
    else
      for i from 1 to nops(c) do
        ans := {solve(f=c[i], y) };
        for j from 1 to nops(ans) do
          surface := surface union{[x, ans[j], z]};
        od;
      od;
      plot3d(surface, h, z=-2..2, grid=[12,12]);
    fi;
  fi;
end:

end module: