I'm not sure what you're trying to do with this code. I think you want your loop variable to be i, not date (which leaves i undefined). I don't think there's any advantage to making nested "display"'s, as you seem to be doing. Instead, you could construct lists of circles, and then each time you want to make a graphics file you call "display" once with the WorldMap_Complete and the appropriate circles. You can indeed save animations as .gif files.

I don't know what you mean by "traces" in this context. What is k? If you want to plot a curve given by the intersection of two surfaces, you could use intersectplot in the plots package. Or if you want to use spacecurve you'd have to rewrite the curve in parametric form. Thus:

> with(plots):
  intersectplot(x^2 + y^2 = 1, z = 0, x=-1..1,y=-1..1,z=-1..1);
  spacecurve([cos(t),sin(t),0],t=0..2*Pi);

You might try something like this.

> S := [solve(whatever)];
  S1 := map(allvalues, S);  # in case there are RootOf's involved
  seemsreal:= x -> not has(simplify(fnormal(evalf(x)),zero),I);
  select(seemsreal, S1);

Of course, this can be fooled by roundoff error.

The folder containing these files should be included in libname. For example:

> libname:= libname,"c:\\maple\\FourierTrigSeries";

Irrespective of the "assume" versus "assuming" issue, if you want to plug in a specific number for d in your expression, it's best to use "eval" rather than assigning a value to d. Thus:

> eval(I1, d = 3.0);

-3.0*alpha*lambda^3*(-(b^2+9.00)^(3/2)+(a^2+9.00)^(3/2)) /(a^2+9.00)^(3/2)/(b^2+9.00)^(3/2) If you really do want to remove assumptions on a variable, you can do that as follows:

> assign(d,'d'): unassign('d'):

With A as before:

> deg:= Pi/180;
  Globe:= (lon,lat) ->(cos(lon)*cos(lat), sin(lon)*cos(lat),
      sin(lat));
  with(plots):
  n:= ArrayDims(A)[1];
  dates := {seq(A[i,1],i=n)};
  for date in dates do
    plotsetup(gif, plotoutput=cat("globe",
      StringTools[SubstituteAll](date,"/","-"),".gif"),
      plotoptions="height=300,width=300");
    display([plottools[sphere]([0,0,0],1, style=patchnogrid),
      seq(circonsphere(0.05, A[i,2]*deg, A[i,3]*deg,
      colour=blue), i = select(j -> (A[j,1]=date),[$n]))],
      scaling=constrained);
  end do;
  plotsetup(default);

p is sqrt(x0^2 + y0^2), which is used in calculating C. C is a parametric equation for the circle in cartesian coordinates. Let P = [x0, y0, z0] be the point on the sphere with longitude lon and latitude lat. In the case where P is not one of the poles, the three vectors P, U = [-y0/p, x0/p, 0] and V = [x0*z0/p, y0*z0/p, p] form an orthonormal set. I then take C = sqrt(1-r^2)*P + r*cos(t)*U + r*sin(t)*V. Do you mean you want one plot with all the points, or separate plots for each? I'll assume it's the former. If you have used ImportMatrix to get an n x 3 Matrix A whose rows are the [date, longitude, latitude] values, for each of these you just need to call circonsphere with the longitude and latitude values, and use display to combine all of the results with a globe. Thus:

> deg:= Pi/180;
  Globe:= (lon,lat) ->(cos(lon)*cos(lat), sin(lon)*cos(lat),
      sin(lat));
  with(plots):
  display([plottools[sphere]([0,0,0],1, style=patchnogrid), 
    seq(circonsphere(0.05, A[i,2]*deg, A[i,3]*deg,
      colour=blue), i=ArrayDims(A)[1])]);

If you only want to plot some of the circles, use some of the rows...

Here's a procedure for plotting a circle of a given radius on the surface of the unit sphere, centred at a point given in terms of longitude and latitude.

> circonsphere:= proc(r, lon, lat)
  local x0,y0,z0,p,C;
  x0, y0, z0 := Globe(lon, lat); # as you defined Globe
  p := sqrt(x0^2 + y0^2);
  if p = 0 then C := [r*cos(t),r*sin(t),sqrt(1-r^2)*z0]
  else C:= [sqrt(1-r^2)*x0 - r*(y0*cos(t)+x0*z0*sin(t))/p,
            sqrt(1-r^2)*y0 + r*(x0*cos(t)-y0*z0*sin(t))/p,
            sqrt(1-r^2)*z0 + r*sin(t)*p]
  end if;
  plots[spacecurve](C,t=0..2*Pi,args[4..nargs]);
  end proc;

Your longitude and latitude seem to be in degrees, so these must be converted to radians. Thus:

> deg:= Pi/180;
  Globe:= (lon,lat) ->(cos(lon)*cos(lat), sin(lon)*cos(lat),
      sin(lat));
  with(plots):
  display(plottools[sphere]([0,0,0],1, style=patchnogrid), 
   circonsphere(0.1, -21.5750*deg, -60.9340*deg, colour=blue),
   orientation=[30,120]);

I answered this in the other thread: the error message shows that the problem is not the undefined elements, but the fact that A is not an Array.

It seems your plot3d structure doesn't use an Array, but rather a list of lists. I'm not sure what circumstances might cause this (just having undefined entries is not sufficient). Perhaps you're using an earlier version of Maple. In any case, the fix is easy enough. If you are using an earlier version (10 or less) of Maple, max['defined'] won't work (this causes max to omit any undefined entries), so that will have to be changed too. Try this:

>  maxpt:= proc(P :: specfunc(anything, PLOT3D))
     local grid, a,b,c,d,A,Adims,zmax,maxpts;
     grid := indets(P, specfunc(anything, GRID));
     if nops(grid) <> 1 then 
       error "There were %1 GRID objects",nops(grid)
     end if;
     grid := op(grid);
     a:= op([1,1],grid);
     b:= op([1,2],grid);
     c:= op([2,1],grid);
     d:= op([2,2],grid);
     A:= op(3,grid);
     if type(A,Array) then Adims:= ArrayDims(A)
     else Adims:= 1..nops(A), 1..nops(A[1])
     end if;
     zmax :=  
       max(seq(seq(`if`(type(A[i,j],undefined),NULL,A[i,j]),
         i=Adims[1]), j=Adims[2]));
     maxpts:= select((p -> A[p[1],p[2]] = zmax),  
             [seq(seq([i,j],i=Adims[1]),j=Adims[2])]);     
     map(p -> [a + (p[1]-1)/(rhs(Adims[1])-1)*(b-a),
               c + (p[2]-1)/(rhs(Adims[2])-1)*(d-c),
               zmax], maxpts);
     end proc;

The problem is that the 2D and 3D plot structures are arranged differently. In the 2D case you have CURVES(...) where the first operand is a list of points. In the 3D case you have GRID(...) where the first two operands give the intervals for the x and y axes and the third is a two-dimensional Array of the z values. Here's a procedure that will work to get the maximum point(s) from a 3D plot containing a GRID object.

> maxpt:= proc(P :: specfunc(anything, PLOT3D))
    local grid, a,b,c,d,A,Adims,zmax,maxpts;
    grid := indets(P, specfunc(anything, GRID));
    if nops(grid) <> 1 then 
      error "There were %1 GRID objects",nops(grid)
    end if;
    grid := op(grid);
    a:= op([1,1],grid);
    b:= op([1,2],grid);
    c:= op([2,1],grid);
    d:= op([2,2],grid);
    A:= op(3,grid);
    Adims:= ArrayDims(A);
    zmax := max['defined'](seq(seq(A[i,j],i=Adims[1]),
            j=Adims[2]));
    maxpts:= select((p -> A[p[1],p[2]] = zmax),  
            [seq(seq([i,j],i=Adims[1]),j=Adims[2])]);     
    map(p -> [a + (p[1]-1)/(rhs(Adims[1])-1)*(b-a),
              c + (p[2]-1)/(rhs(Adims[2])-1)*(d-c),
              zmax], maxpts);
    end proc;

So for example:

> maxpt(plot3d(sin(x) - sin(y), x=0 .. 2, y = 0 .. 2));

[[1.583333333, 0., .999921412762874716]] Caution: the result is of course only approximate, because the function is only calculated on a certain grid of points.

I don't know what your FOC is, but I assume it's an expression containing variables p, s, pa and Q. So Q_opt is a procedure that should work for numerical values of its three arguments. It won't work if you call Q_opt with symbolic arguments, though. Now you want to plot Q_opt(s_val, 0, 0.6) for s_val from 0.5 to 0.59 without actually calling Q_opt with that symbolic argument s_val. This can be done if you give plot a procedure rather than an expression as its first argument. You could write that procedure as s -> Q_opt(s, 0, 0.6) or you could use rcurry(Q_opt, 0, 0.6). The second argument must then be the interval 0.5 .. 0.59 rather than an equation s_val = 0.5 .. 0.59. If you want to label the x axis, you can do that using the labels option. Thus:

 plot(rcurry(Q_opt, 0, 0.6), 0.5 .. 0.59, axes = box,
   labels = [s_val, Q]);

ought to work, assuming fsolve returns one solution for each s_val in this interval.

The result of your crossprod command is a vector. Now solve won't work on a vector, it wants a set or list of expressions or equations. So you could say

> A := crossprod([5,1,-2],[x-0,y-1,z-(-1)]);
  solve(convert(A,set));

But this isn't really what you want either. I suggest you think again about the mathematics behind this, before turning the Maple. What can you say about the normal to the plane if the plane is parallel to (5,1,-2)?

Caution: RealDomain doesn't always work as you might expect. For example, the cubic polynomial x^3 - 3*x + 1 has three real roots, but the closed-form expressions for them in terms of radicals use complex numbers. Therefore RealDomain does not return any roots:

> with(RealDomain):
  solve(x^3 - 3*x + 1 = 0);

The fixed points of f are the solutions of f(x) = x. For a rational function, solve should do it (although Groebner[Solve] may be better in some cases). Otherwise there is fsolve for numerical approximation of solutions. Here's an example with x in R^4:

> F:= [x4 - x1*x3, x1 - x2*x4, x2 - x3*x1, x3 - x4*x2];
  S:= solve(F - [x1,x2,x3,x4]);

S := {x2 = 0, x3 = 0, x1 = 0, x4 = 0}, {x1 = RootOf(_Z^2+2*_Z+2), x4 = -RootOf(_Z^2+2*_Z+2)-2, x2 = -RootOf(_Z^2+2*_Z+2)-2, x3 = RootOf(_Z^2+2*_Z+2)}

> allvalues(S[2]);

{x1 = -1+I, x4 = -1-I, x2 = -1-I, x3 = -1+I}, {x2 = -1+I, x3 = -1-I, x1 = -1-I, x4 = -1+I} So in this case the only real solution is [0,0,0,0]. BTW: why can't the Maple Math tag handle expression sequences correctly? I had to insert tags separately around the two allvalues results, otherwise just one shows.