{ legend for partial differentiation:
  L=sqrt(f*m*a)
  n^2*a^3=fm
  L=sqrt(f*m*a)=n*a^2
  G=L*sqrt(1-e^2)=n*a^2/p
  H=G*cos(i)=G*c
  a=L^2/(f*m)
  p=1/sqrt(1-e^2)=L/G
  e^2=1-G^2/L^2
  c=H/G
  s^2=1-H^2/G^2
  n=(f*m)^2/L^3
  da/dL=2*L/fm=2/(n*a)   ! factor +2
  dn/dL=-3/a^2           ! factor -3
  dp/dL=1/G=p/(n*a^2)
  dp/dG=-L/G^2=-p^2/(n*a^2)
  ds/dG=(1/s)*(H^2/G^3)=p*c^2/(s*n*a^2)
  ds/dH=-1/s*(H/G^2)=-p*c/(s*n*a^2)
  dc/dG=-H/G^2=-p*c/(n*a^2)
  dc/dH=1/G=p/(n*a^2)
  de/dL=((1-e^2)/e)*(1/L)=1/(p^2*e*n*a^2)
  de/dG=-(1/e)*(sqrt(1-e^2)/L)=-1/(p*e*n*a^2) }
unit untheled ; { for partial derivations }
interface       { numerical values }
uses
  unsagcot ; { type tvarsix }
{ we use numerical values
  of canonical elements L G H 
  to obtain numerical values of partial derivations
  from parameters  p a e s c  with respect to L G H }
procedure simpartial ( e  : tvarsix ) ; 
implementation
uses
  untypvar , { global amplia amplip ..e ..s ..c ..n[..] }
  ungiamod ; { global value PlanetFMain instead fm for giant planet }
procedure allparmnul ; { initial nullo to arrays }
var
  i  : byte ;  { count for three elements L G H }
begin            { e eccentricity i inclination }
  for i:=1 to 3 do { nullo to amplitudes }
    begin
      amplia[i]:=0.0; { arrays global from untypvar }
      amplin[i]:=0.0;
      amplip[i]:=0.0;
      amplie[i]:=0.0;
      amplis[i]:=0.0;
      amplic[i]:=0.0;
    end;
end;
{ we use numerical values
  of canonical elements e[1]=L e[2]=G e[3]=H 
  to obtain numerical values of partial derivations
  from parameters  p a e s c  with respect to L G H }
procedure simpartial ( e  : tvarsix ) ; { type unsagcot L G H l g h }
var                  { e[1]=L e[2]=G e[3]=H e[4]=l e[5]=g e[6]=h }
  p,q  : extended ; { temporary values }
begin          { value for fm in global PlanetFMain }
  allparmnul ; { initial nullo to arrays amplia ..p ..[1..3] }
  { da/dL=2*L/fm  da/dG=0  da/dH=0 }
  amplia[1]:=+2.0*e[1]/PlanetFMain; { da/dL }
  { dp/dL=1/G  dp/dG=-L/G^2  dp/dH=0 }
  amplip[1]:=+1.0/e[2];       { dp/dL=1/G }
  amplip[2]:=-e[1]/sqr(e[2]); { dp/dG=-L/G^2 }
  { de/dL=((1-e^2)/e)*(1/L)  de/dG=-(1/e)*(sqrt(1-e^2)/L)  de/dH=0 }
  p:=sqr(e[2]/e[1]); { 1-e^2=G^2/L^2 }
  q:=sqrt(1.0-p);      { e=sqrt(1-G^2/L^2) }
  amplie[1]:=+(p/q)*(1.0/e[1]);       { de/dL=((1-e^2)/e)*(1/L) }
  amplie[2]:=-(sqrt(p)/q)*(1.0/e[1]); { de/dG=-(1/e)*(sqrt(1-e^2)/L) }
  { n=(f*m)^2/L^3  dn/dL=-3*(fm)^2/L^4  dn/dG=0  dn/dH=0 }
  amplin[1]:=-3.0*sqr(PlanetFMain)/sqr(e[1]*e[1]); { dn/dL=-3*(fm)^2/L^4 }
  { s^2=1-H^2/G^2  ds/dL=0  ds/dG=(1/s)*(H^2/G^3)  ds/dH=-1/s*(H/G^2) }
  p:=sqr(e[3]/e[2]); { 1-s^2=H^2/G^2 }
  q:=sqrt(1.0-p);    { sin(i)=sqrt(1-H^2/G^2) }
  amplis[2]:=+(1.0/q)*(p/e[2]); { ds/dG=(1/s)*(H^2/G^3) }
  amplis[3]:=-(1.0/q)*(p/e[3]); { ds/dH=-(1/s)*(H/G^2) }
  { c=H/G  dc/dL=0  dc/dG=-H/G^2  dc/dH=1/G }
  amplic[2]:=-p/e[3];   { dc/dG=-H/G^2 }
  amplic[3]:=+1.0/e[2]; { dc/dH=1/G }
end;
end.