This post is related to the question. There were  proposed two ways of finding the volume of the cutted part of a sphere in the form of a wedge.  Here the procedure is presented that shows the rest of the sphere. Parameters procedure: R - radius of the sphere, H1 - the distance the first cutting plane to the plane  xOy,  H2 -  the distance the second cutting plane to the plane  zOy. Necessary conditions:  R>0,  H1>=0,  H2>=0,  H1^2+H2^2<R^2 . For clarity, different surfaces are painted in different colors.

restart;

Pic := proc (R::positive, H1::nonnegative, H2::nonnegative)

local A, B, C, E, F;

if R^2 <= H1^2+H2^2 then error "Should be H1^(2)+H2^(2)<R^(2)" end if;

A := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = arctan(sqrt(-H1^2-H2^2+R^2), H2) .. 2*Pi-arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. Pi, color = green);

B := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = 0 .. arccos(sqrt(R^2-H2^2-H2^2*tan(phi)^2)/R), color = green);

C := plot3d([R*sin(theta)*cos(phi), R*sin(theta)*sin(phi), R*cos(theta)], phi = -arctan(sqrt(-H1^2-H2^2+R^2), H2) .. arctan(sqrt(-H1^2-H2^2+R^2), H2), theta = arccos(H1/R) .. Pi, color = green);

E := plot3d([r*cos(phi), r*sin(phi), H1], phi = -arccos(H2/sqrt(R^2-H1^2)) .. arccos(H2/sqrt(R^2-H1^2)), r = H2/cos(phi) .. sqrt(R^2-H1^2), color = blue);

F := plot3d([H2, r*cos(phi), r*sin(phi)], phi = arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)) .. Pi-arccos(sqrt(-H1^2-H2^2+R^2)/sqrt(R^2-H2^2)), r = H1/sin(phi) .. sqrt(R^2-H2^2), color = gold);

plots[display](A, B, C, E, F, axes = none, view = [-1.5 .. 1.5, -1.5 .. 1.5, -1.5 .. 1.5], scaling = constrained, lightmodel = light4, orientation = [60, 80]);

end proc:

Example of use:

Pic(1,  0.5,  0.3);

The work contains two procedures called  SetPartition  and  NumbPart .

The first procedure  SetPartition  generates all the partitions of a set  S  into disjoint subsets of specific sizes defined by a list  Q. The third optional parameter is the symbol  `ordered`  or  `nonordered`  (by default) . It determines whether the order of subsets in a partition is significant.

The second procedure  NumbPart  returns the number of all partitions of the sizes indicated by  Q .

Codes of these procedures:

SetPartition:=proc(S::set, Q::list(posint), K::symbol:=nonordered)  # Procedure finds all partitions of the set  S  into subsets of the specific size given Q.

local L, P, T, S1, S2, n, i, j, m, k, M;

uses ListTools,combinat;

if `+`(op(Q))<>nops(S) then error "Should be `+`(op(Q))=nops(S)" fi;

L:=convert(Q,set);

T:=[seq([L[i],Occurrences(L[i],Q)], i=1..nops(L))];

P:=`if`(K=ordered,convert(T,list),convert(T,set));

S1:=S;  S2:={`if`(K=ordered,[],{})}; n:=nops(P);

for i to n do

m:=P[i,1]; k:=P[i,2];

for j to k do

S2:={seq(seq(`if`(K=ordered,[op(s),t],{op(s),t}), t=choose(S1 minus `union`(op(s)),m)), s=S2)};

od; od;

if K=ordered then {map(op@permute,S2)[]} else S2 fi;

end proc:

NumbPart:=proc(Q::list(posint), K::symbol:=nonordered)  # Procedure finds the number of all partitions of a set into subsets of the specific size given  Q

local L, T, P, n, S, N;

uses ListTools;

L:=convert(Q,set);

T:=[seq([L[i],Occurrences(L[i],Q)], i=1..nops(L))];

P:=`if`(K=ordered,convert(T,list),convert(T,set));

n:=nops(P);  N:=add(P[i,2], i=1..n);

S:=add(P[i,1]*P[i,2],i=1..n)!/mul(P[i,1]!^P[i,2],i=1..n)/mul(P[i,2]!,i=1..n);

if K=nonordered then return S else  S*N! fi;

end proc:

Examples of use:

SetPartition({a,b,c,d,e}, [1,1,3]);  nops(%);  # Nonordered partitions and their number

SetPartition({a,b,c,d,e}, [1,1,3], ordered);  nops(%);  # Ordered partitions and their number

Here's a more interesting example. 5 fruits  {apple, pear, orange, mango, peach}  must be put on three plates so that  on each of two plates there are 2  fruits, and there is one fruit  on one plate. Two variants to solve: 1) plates are indistinguishable and 2) different plates. In how many ways can this be done?

SetPartition({apple,pear,orange,mango,peach}, [1,2,2]);  nops(%);  # plates are indistinguishable

NumbPart([1,2,2], ordered);  # Number of ways for  different plates

                                                              90

Another example - how many ways can be divided  100 objects into 10 equal parts?

NumbPart([10$10]);

    64954656894649578274066349293466217242333450230560675312538868633528911487364888307200

 SetPartition.mws

The procedure  Partition  significantly generalizes the standard procedure  combinat[partition]  in several ways. The user specifies the number of parts of the partition, and can also set different limitations on parts partition.

Required parameters:  n - a nonnegative integer, k - a positive integer or a range (k  specifies the number of parts of the partition). The parameter  res  is the optional parameter (by default  res is  1 ). If  res  is a number, all elements of  k-tuples must be greater than or equal  res .  If  res  is a range  a .. b ,   all elements of  k-tuples must be greater than or equal  a  and  less than or equal  b . The optional parameter  S  - set, which includes elements of the partition. By default  S = {$ 0.. n} .

The code of the procedure:

Partition:=proc(n::nonnegint, k::{posint,range}, res::{range, nonnegint} := 1, S::set:={$0..n})  # Generates a list of all partitions of an integer n into k parts

local k_Partition, n1, k1, L;

k_Partition := proc (n, k::posint, res, S)

local m, M, a, b, S1, It, L0;

m:=S[1]; M:=S[-1];

if res::nonnegint then a := max(res,m); b := min(n-(k-1)*a,M)  else a := max(lhs(res),m); b := min(rhs(res),M) fi;

S1:={$a..b} intersect S;

if b < a or b*k < n or a*k > n  then return [ ] fi;

It := proc (L)

local m, j, P, R, i, N;

m := nops(L[1]); j := k-m; N := 0;

for i to nops(L) do

R := n-`+`(op(L[i]));

if R <= b*j and a*j <= R then N := N+1;

P[N] := [seq([op(L[i]), s], s = {$ max(a, R-b*(j-1)) .. min(R, b)} intersect select(t->t>=L[i,-1],S1) )] fi;

od;

[seq(op(P[s]), s = 1 .. N)];

end proc;

if k=1 then [[b]] else (It@@(k-1))(map(t->[t],S1))  fi;

end proc;

if k::posint then return k_Partition(n,k,res,S) else n1:=0;

for k1 from lhs(k) to rhs(k) do

n1:=n1+1; L[n1]:=k_Partition(n,k1,res,S)

od;

L:=convert(L,list);

[seq(op(L[i]), i=1..n1)] fi;

end proc:

Examples of use:

Partition(15, 3);

Partition(15, 3..5, 1..5);  # The number of parts from 3 to 5, and each summand from 1 to 5

Partition(15, 5, {seq(2*n-1, n=1..8)});  # 5 summands and all are odd numbers 

A more interesting example.
There are  k banknotes in possible denominations of 5, 10, 20, 50, 100 dollars. At what number of banknotes  k  the number of variants of exchange  $140  will be maximum?

n:=0:

for k from 3 to 28 do

n:=n+1: V[n]:=[k, nops(Partition(140, k, {5,10,20,50,100}))];

od:

V:=convert(V, list);

max(seq(V[i,2], i=1..nops(V)));

select(t->t[2]=8, V);

Here are these variants:

Partition(140, 10, {5,10,20,50,100});

Partition(140, 13, {5,10,20,50,100});

 Partition.mws

A heart shape in 3d:

The code of the animation:

A := plots[animate](plot3d, [[16*sin(t)^3*cos(s), 16*sin(t)^3*sin(s), 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)], t = 0 .. u, s = 0 .. 2*Pi, color = red, style = surface, axes = none], u = 0 .. Pi, frames = 100):

B := plots[animate](plot3d, [[16*sin(t)^3*cos(s), 16*sin(t)^3*sin(s), 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t)], t = u .. Pi, s = 0 .. 2*Pi, color = "LightBlue", style = surface, axes = none], u = 0 .. Pi, frames = 100):

plots[display](A, B);

Edited. The direction of painting changed.

The well-known  combinat[composition]  command computes and returns a list containing all distinct ordered  k-tuples of positive integers whose elements sum equals  n . These are known as the compositions of  n .  For some applications, additional constraints are required for the elements of these k-tuples, for example, that they are within a certain range.

The  Composition  procedure solves this problem. Required parameters:  n - a nonnegative integer, k - a positive integer. The parameter  res  is the optional parameter (by default  res is  0 ). If  res  is a number, all elements of  k-tuples must be greater than or equal  res .  If  res  is a range  a .. b ,   all elements of  k-tuples must be greater than or equal  a  and  less than or equal  b .  Composition(n,k,1)  is equivalent to  combinat[composition](n,k) .

The code of the procedure:

Composition := proc (n::nonnegint, k::posint, res::{range, nonnegint} := 0)

local a, b, It, L0; 

if res::nonnegint then a := res; b := n-(k-1)*a  else a := lhs(res); b := rhs(res) fi;

if b < a or b*k < n then return `No solutions` fi; 

It := proc (L)

local m, j, P, R, i, N;

m := nops(L[1]); j := k-m; N := 0;

for i to nops(L) do

R := n-`+`(op(L[i]));

if R <= b*j and a*j <= R then N := N+1;

P[N] := [seq([op(L[i]), s], s = max(a, R-b*(j-1)) .. min(R, b))] fi;

od;

[seq(op(P[s]), s = 1 .. N)];

end proc;

L0 := [[]];

(It@@k)(L0); 

end proc:

Three simple examples:

Composition(10,3); ``;   # All terms greater than or equal 0

Composition(10,3, 2);   # All terms greater than or equal 2

Composition(10,3, 2..4);   # All terms greater than or equal 2 and less than or equal to 4 

A more complex example. The problem - to find all the numbers in the range  1 .. 99999999  whose digits sum is equal to 21 .

Each number is represented by a list of digits from left to right, replacing missing digits at the left with zeros.

M:=Composition(21,8, 0..9):  

nops(M);  # The number of solutions

[seq(M[1+100000*i], i=0..9)]; # 10 solutions from the list M starting the first one

seq(add(%[i,k]*10^(8-k), k=1..8),i=1..nops(%));  # Conversion into numbers

Composition.mws