# Composite function derivatives

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