calculus and analysis – How do I take a derivative, substitute a function into the initial function, and take a subsequent derivative

The title may be a little confusing so let me explain what I mean:

I am using lagrange brackets to produce some equations of motion for orbital mechanics, where the motion is split between oscillating and mean elements. I’m doing the same thing as brouwer theory if you’re familiar with that.

Basically motion is decomposed into a mean element which varies linearly in time and a quickly var oscillating portion.

Anyways, to make the oscillating state elements I need to take derivatives of a generating function. I also further would like to take the derivative of the oscillating elements with respect to the mean elements for least squares routines to estimate mean elements.

So I have been able to successfully produce the oscillating elements, but when attempting to take the derivative of those I run into a problem.

To get the oscillating elements I have to take the lagrange brackets of my elements with the generaging function ie <ei,W> where the derivatives are in denaulay variables so

Where ei are the state vector I am using, W is the generating function and l, g, h,L,G,H are a seperate set of state variables other than those in ei called delaunay variables. So far so good

Next I want to use least squares to find the mean variables. To do so as stated above I need the derivative of the oscilating to mean elements. To do this I need to define the l, g, h, L, G, and H in terms of my state variables, and then take the derivative.

Unfortunately since I’m using DelaySet, mathematica attempts to take the derivatives all at once at the end, which cause the derivative to fail, as it cannot take the derivative with respect to l, g, h, L, G, and H as they are now functions of my variables ei.

How do I fix this without hand transcription? I have tried using evaluate around the derivatives with respect to l, g, h, L, G, and H, and using set rather than delay set (as in = instead of :=). I have not been able to get the desired behavior

Here is a portion of the code cut out to illustrate whats not working:

L(n_, (Mu)_) := Sqrt((Mu)*a(n, (Mu)));

e(h_, k_) := Sqrt(h^2 + k^2);

G(n_, (Mu)_, h_, k_) := L(n, (Mu))*(1 - e(h, k)^2)^(1/2);

i(p_, q_) := 2*ArcTan(p^2 + q^2);
H (n_, (Mu)_, h_, k_, p_, q_) := G(n, (Mu), h, k)*Cos(i(p, q))
xsi(h_, k_) := ArcTan(h/k);
l(h_, k_, lm_) := lm - xsi(h, k);
g(q_, p_) := ArcTan(p/q);

hDel(h_, k_, q_, p_) := Simplify(xsi(h, k) - g(p, q))

W1lp1(l1_, g1_, hDel1_, L1_, G1_, H1_, n_, h_, k_, q_, p_, 
 lm_) := -(1/(32*G1^3))*(1 - G1^2/L1^2)*(1 - 5*H1^2/G1^2)^-1*(1 - 
   16*H1^2/G1^2 + 15*H1^4/G1^4)*Sin(2*g1)
Mean to Oscillating Transformation Long period  portion

nW1lp1(l1_, g1_, hDel1_, L1_, G1_, H1_, n_, h_, k_, q_, p_, lm_, 
  mu_) := -dndL(n, mu)*
   D(W1lp1 (l1, g1, hDel1, L1, G1, H1, n, h, k, q, p, lm), l1) - 
  dndG(n, mu, k, h)*
   D(W1lp1(l1, g1, hDel1, L1, G1, H1, n, h, k, q, p, lm), g1) - 
  dndH (n, mu, h, k, p, q)*
   D(W1lp1 (l1, g1, hDel1, L1, G1, H1, n, h, k, q, p, lm), hDel1);

DnOscDnMeanlp1(n_, h_, k_, q_, p_, lm_, mu_, J2_, rE_) := 
   nW1lp1(l(h, k, lm), g(q, p), hDel(h, k, q, p), L(n, mu), 
    G(n, mu, h, k), H(n, mu, h, k, p, q), n, h, k, q, p, lm, mu, J2, 
    rE), n);

The problem is with the last function as I can get the ones above to work. I’ve also posted the full code in a comment below