# How do I ensure the substitution happens after I take the first derivative before taking my second 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 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_] :=
J2*rE^2*D[
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

Edit: I believe a better way to ask this is the following: How do I ensure the substitution of the delaunay variables happens after I take the first derivative for the lagrange brackets and then take the second derivative with respect to my orbital elements.

Edit 2: I've continued to update the code on github with my latest attempt

Edit 3: As some of yall have upvoted the response to post the code I decided to put it into the main area as well: https://github.com/pjdog/mathematicaStuff/blob/main/Just%20the%20code%20copied%20as%20github%20doesnt%20display%20mathematica%20well .