How to write a function in the following form f[a_, b_, c_,...] so that it is flexible with respect to the number of arguments?

Question

@lericr have written under several of my questions that the way I try to modernize the code could be made much simpler. Below is the original part of the code. The remaining parts have similar designs, so if what I would like to get can be simplified, then this can be applied to other parts.

The goal I would like to achieve is to have the function EER[a1_, b1_, qa_, qb_, c11_, c12_, c13_,..., c21_, c22_, c23_,...] that depends on an arbitrary number of arguments of this type cij. The number of variables cij depends on the imax and jmax.

The function EER is the sum of the Kxx and Pxx each of which depends on variables a1, b1, qa, qb, c11, c12, c13,..., c21, c22, c23,.... I managed to present the function Kxx in such a way that it is flexible in relation to the number of variables (that is, when changing imax and jmax, a specific value of this function can be easily obtained by substituting the numerical values of all variables). But I can't do the same procedure with the function Pxx.

Once again, briefly goal: to make the function Pxx[a1_, b1_, qa_, qb_, c11_, c12_, c13_,..., c21_, c22_, c23_,...] universal in relation to the number of variables, and also, if possible, to reasonably simplify the code.

The code is written for imax=1 and jmax=2. If you change them to imax=2 and jmax=3 for example, you can see that the function Kxx is flexible and you can immediately get its values (in:Kxx[1.22, 0.44, 1.40, 1.81, 1, 2, 1, 3, 2, 1] out: 281.248), but this is not done with the function Pxx (in:Pxx[1.22, 0.44, 1.40, 1.81, 1, 2, 1, 3, 2, 1] out: Pxx[1.22, 0.44, 1.4, 1.81, 1, 2, 1, 3, 2, 1]).

ClearAll["Global`*"]
imax = 1; jmax = 2;

Psi[r_, z_, i_, j_] := Exp[-b[j]*z^2]*Exp[-a[i]*r^2];

(*variables cij*)
c[i_, j_] := Symbol["c" <> ToString[i] <> ToString[j]]
var = Join[{a1, b1, qa, qb}, 
   Flatten[Table[c[i, j], {i, 1, imax}, {j, 1, jmax}]]];
var1 = Pattern[#, Blank[]] & /@ var;

(*a[i] and b[i] are geometric progressions*)
a[1] = a1;
b[1] = b1;
Do[a[i] = a[i - 1] qa, {i, 2, imax}];
Do[b[j] = b[j - 1] qb, {j, 2, jmax}];

Kk = FullSimplify[
   Psi[r, z, i2, j2] *
    Laplacian[Psi[r, z, i1, j1] , {r, \[Theta], z}, "Cylindrical"]];
Kk1 = Integrate[Kk r, {r, 0, \[Infinity]}, {z, -Infinity, Infinity}, 
   Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, 
     a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}];
Kx = -1/2 2 Pi Sum[
    c[i1, j1] c[i2, j2] Kk1, {i1, 1, imax}, {i2, 1, imax}, {j1, 1, 
     jmax}, {j2, 1, jmax}];
Kxx[Sequence @@ var1] = Kx;

VB1[r_, z_] := -(1/(2*Sqrt[r^2 + z^2]))*Exp[-Sqrt[r^2 + z^2]*1/2];
Px1[a10_, b10_, qa0_, qb0_, i1_, j1_, i2_, j2_] := 
  Block[{a1 = a10, b1 = b10, qa = qa0, qb = qb0}, 
   NIntegrate[
    Psi[r, z, i2, j2]*VB1[r, z]*Psi[r, z, i1, j1]*r, {r, 
     0, \[Infinity]}, {z, -Infinity, Infinity}]];
Pxx[a10_, b10_, qa0_, qb0_, c111_, c121_] := 
  Block[{c11 = c111, c12 = c121}, 
   2 Pi Sum[
     c[i1, j1] c[i2, j2] Px1[a10, b10, qa0, qb0, i1, j1, i2, j2], {i1,
       1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}]];

With[{var = var}, 
  EER[Sequence @@ var1] := 
   N[Kxx[Sequence @@ var] + Pxx[Sequence @@ var]]];

(*imax=1, jmax=2*)
(*Kxx[1.22,0.44,1.40,1.81,1,2]*)
(*27.549720475849007` *)
(*Pxx[1.22,0.44,1.40,1.81,1,2]*)
(*-10.388782043150181` *)

(*imax=2, jmax=3*)
(*Kxx[1.22,0.44,1.40,1.81,1,2,1,3,2,1]*)
(*281.2480800909348` *)
(*Pxx[1.22,0.44,1.40,1.81,1,2,1,3,2,1]*)
(*Pxx[1.22`,0.44`,1.4`,1.81`,1,2,1,3,2,1]*)

(*imax=1, jmax=2*)
(*EER[1.22,0.44,1.40,1.81,1,2]*)
(*17.160938432698828` *)

(*imax=2, jmax=3*)
(*EER[1.22,0.44,1.40,1.81,1,2,1,3,2,1]*)
(*EER[1.22`,0.44`,1.4`,1.81`,1,2,1,3,2,1]*)

Answer 1

Okay, this will be long and involved. I'm basing this as much as I can on your code so that you don't have to try to translate too much. But I've only gotten halfway to your final objective. I hope you can follow this pattern to get the rest of the way (ask further questions if you can't). There are many other improvements we could make, but I wanted to stay focused as much as possible on how to write an extensible function.

Our first objective is to rewrite this:

Kxx[Sequence @@ var1] = Kx

I'm deferring tackling the Pxx and EER and so forth.

After inspecting your var1, we see that has four arguments at the front and then a variable number of arguments after that. But those variable arguments can always be made with a Table and iteration. That means they're really just a matrix. So, let's try that:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] := Kx

Note the :=. That suspends evaluation until we need it. Now let's plug in your Kx:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
  -1/2 2 Pi Sum[
    c[[i1, j1]] c[[i2, j2]] Kk1,
    {i1, 1, imax}, {i2, 1, imax}, {j1, 1, jmax}, {j2, 1, jmax}]

Notice I changed [] to [[]], because now we're accessing a matrix. Now consider that imax and jmax can be derived from c. Let's use that fact:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
  With[
    {dims = Dimensions[c]},
    -1/2 2 Pi Sum[
      c[[i1, j1]] c[[i2, j2]] Kk1,
      {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]]

Time to plug in your Kk1, but first let's inspect it on its own:

Integrate[Kk r, 
  {r, 0, \[Infinity]}, {z, -Infinity, Infinity},
  Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}]

We need definitions for a and b. Your definitions look like this:

a[1] = a1;
b[1] = b1;
Do[a[i] = a[i - 1] qa, {i, 2, imax}];
Do[b[j] = b[j - 1] qb, {j, 2, jmax}];

There is a problem here and an opportunity. The problem is that these definitions refer to a1, qa, etc. But as you see above we're passing in those arguments to Kxx, so we need to feed them to a and b. The opportunity is that a and b are simple functions, so we don't need to use a loop to pre-define values for them. In fact, they're the same function. As you noted, they are geometric sequences. So, let's do this:

geoseq[init_, r_, n_Integer?Positive] := init*r^(n - 1)

Let's plug this into our Kxx, but let's keep it slightly simpler by redefining a and b locally:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
  With[
    {dims = Dimensions[c],
     a = geoseq[a1, qa, #] &,
     b = geoseq[b1, qb, #] &},
    -1/2 2 Pi Sum[
      c[[i1, j1]] c[[i2, j2]] Integrate[Kk r, {r, 0, Infinity}, {z, -Infinity, Infinity}, 
        Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}],
      {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]]

Okay, now for Kk. It looks like this:

FullSimplify[Psi[r, z, i2, j2]*Laplacian[Psi[r, z, i1, j1], {r, \[Theta], z}, "Cylindrical"]]

So, we need Psi. Your Psi depended on a and b. But we've replaced those with geoseq. And now this should be déjà vu... we need to add arguments:

Psi[a1_, b1_, qa_, qb_, r_, z_, i_, j_] := Exp[-geoseq[b1, qb, j]*z^2]*Exp[-geoseq[a1, qa, i]*r^2]

Okay, putting it all together, and hoping I can do that without a copy-paste error:

Kxx[a1_, b1_, qa_, qb_, c_List?MatrixQ] :=
  With[
    {dims = Dimensions[c],
     a = geoseq[a1, qa, #] &,
     b = geoseq[b1, qb, #] &},
    -1/2 2 Pi 
    Sum[
      c[[i1, j1]] c[[i2, j2]] 
      Integrate[
        FullSimplify[
          Psi[a1, b1, qa, qb, r, z, i2, j2]*
          Laplacian[
            Psi[a1, b1, qa, qb, r, z, i1, j1], {r, \[Theta], z}, "Cylindrical"]] r, 
        {r, 0, Infinity}, {z, -Infinity, Infinity}, 
        Assumptions -> {a[i1] > 0, b[i1] > 0, a[i2] > 0, b[i2] > 0, a[j1] > 0, b[j1] > 0, a[j2] > 0, b[j2] > 0}],
      {i1, 1, dims[[1]]}, {i2, 1, dims[[1]]}, {j1, 1, dims[[2]]}, {j2, 1, dims[[2]]}]]

Let's take it for a spin:

Kxx[1.22, 0.44, 1.40, 1.81, {{1, 2}}]
(* 27.5497 *)

Notice the last argument is a matrix, we needed to reformulate from what you had. Let's try another:

Kxx[1.22, 0.44, 1.40, 1.81, {{1, 2, 1}, {3, 2, 1}}]
(* 281.248 *)

At this point there are lots of opportunities for clean up and maybe some performance improvement, but I'm going to stop there for now.

Answer 2

Here's an example using a Module construct:

 myFun[varList_List] := Module[{i, j, varLen},
      varLen = Length@varList;
      Print["Total variables: ", varLen];
      Print["List of variables: ", varList];
      If[varLen > 1,
       Print["Last variable value: ", varList[[-1]]];
       ];
      ]
    
   myFun[{1, 2, 3, 4}]

(* Total variables: 4
   List of variables: {1,2,3,4}
   Last variable value: 4 *)

Answer 3

Following josh's answer:

c1[i_, j_] := Symbol["c" <> ToString[i] <> ToString[j] <> "1"]
varc1 = Join[{a1, b1, qa, qb}, 
   Flatten[Table[c1[i, j], {i, 1, imax}, {j, 1, jmax}]]];
varc11 = Pattern[#, Blank[]] & /@ var;
FF = Table[c[i, j] -> c1[i, j], {i, 1, imax}, {j, 1, jmax}]

Pxx[Sequence @@ varc11] := 
 Module[{Sequence @@ varc1}, PxB1[a10, b10, qa0, qb0] /. FF]

(*Pxx[1.22, 0.44, 1.40, 1.81, 1, 2]*)

(*Module[{Sequence @@ varc1}, PxB1[a10, b10, qa0, qb0] /. FF]*)