Composite function derivatives

Question

I would like to calculate the numerical partial derivatives (numerical gradient) of this function:

R[M_, d_, \[Theta]_, \[Phi]_, \[Psi]_, \[Iota]_, f_, \[CapitalPhi]_] := 1/\[CapitalOmega] K ((1 + Cos[\[Iota]]^2) P[
       Sin[\[CapitalOmega] t + \[Eta] + \[CapitalPhi]] - 
        Sin[\[Eta] + \[CapitalPhi]]] + 
     2 Cos[\[Iota]] X[
       Cos[\[CapitalOmega] t + \[Eta] + \[CapitalPhi]] - 
        Cos[\[Eta] + \[CapitalPhi]]]) - 
 1/\[Omega] A ((1 + Cos[\[Iota]]^2) P[
       Sin[\[Omega] t + \[CapitalPhi]] - Sin[\[CapitalPhi]]] + 
     2 Cos[\[Iota]] X[
       Cos[\[Omega] t + \[CapitalPhi]] - Cos[\[CapitalPhi]]])

in which K, A, P, X, [Eta], [Tau], [Omega] are functions of the same variables of R, for example:

\[Omega] = N[(2*f*\[Pi]];
\[Eta] = N[(-2*\[Pi]*
    f*(L/c)*(1 +  {-Sin[\[Theta]] Cos[\[Phi]], -Sin[\[Theta]] Sin[\
\[Phi]], -Cos[\[Theta]]}.{-Sin[zen] Cos[az], -Sin[zen] Sin[az], -Cos[
          zen]}))];

A = N[1.3*10^-15*(f/10^-7)^(2/3)*(M/(s*10^9))^(5/3)*(d/G)^-1];

P = N[(((Cos[zen] Cos[\[Psi]] Sin[\[Theta]] - 
         Sin[az] Sin[
           zen] (Cos[\[Theta]] Cos[\[Psi]] Sin[\[Phi]] + 
            Cos[\[Phi]] Sin[\[Psi]]) - 
         Cos[az] Sin[
           zen] (Cos[\[Theta]] Cos[\[Phi]] Cos[\[Psi]] - 
            Sin[\[Phi]] Sin[\[Psi]]))^2 - (Cos[
           zen] Sin[\[Theta]] Sin[\[Psi]] - 
         Cos[az] Sin[
           zen] (Cos[\[Psi]] Sin[\[Phi]] + 
            Cos[\[Theta]] Cos[\[Phi]] Sin[\[Psi]]) - 
         Sin[az] Sin[
           zen] (-Cos[\[Phi]] Cos[\[Psi]] + 
            Cos[\[Theta]] Sin[\[Phi]] Sin[\[Psi]]))^2)/(2 (1 + 
         Cos[zen] Cos[\[Theta]] + 
         Cos[az] Cos[\[Phi]] Sin[zen] Sin[\[Theta]] + 
         Sin[az] Sin[zen] Sin[\[Theta]] Sin[\[Phi]])))];

X = N[(((Cos[zen] Cos[\[Psi]] Sin[\[Theta]] - 
         Sin[az] Sin[
           zen] (Cos[\[Theta]] Cos[\[Psi]] Sin[\[Phi]] + 
            Cos[\[Phi]] Sin[\[Psi]]) - 
         Cos[az] Sin[
           zen] (Cos[\[Theta]] Cos[\[Phi]] Cos[\[Psi]] - 
            Sin[\[Phi]] Sin[\[Psi]])) (Cos[
           zen] Sin[\[Theta]] Sin[\[Psi]] - 
         Cos[az] Sin[
           zen] (Cos[\[Psi]] Sin[\[Phi]] + 
            Cos[\[Theta]] Cos[\[Phi]] Sin[\[Psi]]) - 
         Sin[az] Sin[
           zen] (-Cos[\[Phi]] Cos[\[Psi]] + 
            Cos[\[Theta]] Sin[\[Phi]] Sin[\[Psi]])))/(1 + 
       Cos[zen] Cos[\[Theta]] + 
       Cos[az] Cos[\[Phi]] Sin[zen] Sin[\[Theta]] + 
       Sin[az] Sin[zen] Sin[\[Theta]] Sin[\[Phi]]))];
\[Tau] = N[
   1/c L (1 + Cos[zen] Cos[\[Theta]] + 
      Cos[az] Cos[\[Phi]] Sin[zen] Sin[\[Theta]] + 
      Sin[az] Sin[zen] Sin[\[Theta]] Sin[\[Phi]])];
K = N[1.3*10^-15*(f - 
      15*(M/(s*10^(8.5)))^(5/3)*(f/n)^(11/3)*\[Tau]*10^(-9)/10^-7)^(
    2/3)*(M/(s*10^9))^(5/3)*(d/G)^-1];

I need to calculate the numerical gradient and evaluate it at a given point, but I can't substitute the parameters in the function or calculate composite derivatives. Thank you in advance for you help