(*
  1  p=1/sqrt(1-e^2)       1  v
  2  e                     2  E
  3  s = sin(i)            3  g
  4  c = cos(i)            4  h
  5  p'=1/sqrt(1-e'^2)     5  v'
  6  e'                    6  E'
  7  s' = sin(i')          7  g'
  8  1-c' = 1-cos(i')      8  h'
  9  a
 10  a'
 11  r
 12  r'
 13  n
 14  n'
 to obtain numerical values
 of partial derivations from F with L G H l g h
*)
unit unthelem ; { correction to canonical elements }
interface
uses
  unsagcot , { types tsagelm tsagano }
  untypvar ; { types trowvar trowptr }
procedure corelem ( xo  : tsagano ;        { L G H l g h canonical }
                    me  : word ;           { number of terms }
                    ro  : trowptr ;        { dynamical array for terms }
                    var  ys  : tvarsix ) ; { corrections to elements }
implementation                             { or partial derivations }
uses
  untheled , { procedure simpartial }
  unthelep , { procedure simparamon }
  untheleg , { procedure giparargum }
  unthelpo ; { procedure allparapow }
{ to sum angles with coefficients to argument value } 
procedure cosargu ( r  : trowvar ;             { current term }
                    var  s , c  : extended ) ; { co sin argument }
var
  i  : byte ;     { count }
  g  : extended ; { argument }
begin             { types and global untypvar }
  g:=0.0;  { initial nullo to argument }
  for i:=1 to 8 do  { for all angles v E g h v' E' g' h' }
    g:=g+r.mar[i]*argall[i]; { our argument in radian }
  s:=Sin(g); {   sinus of our current argument }
  c:=Cos(g); { cosinus of our current argument }
end;
{ for current term
  numerical value of amplitude and derivations }
procedure difacto ( r  : TRowVar ;        { term types untypvar }
                    var  a  : extended ;  { amplitude with power }
                    var  d  : tampdif ) ; { derivations with L G H }
var
  i,j  : byte ; { simple counts }
begin           { all global variables from untypvar }
  { to multiply parameters in power and derivations }
  a:=r.amp[1];                  { initial amplitude value }
  for j:=1 to 3 do  d[j]:=0.0;  { initial nullo to sum derivations }
  for i:=1 to numsel do         { for all parameters }
    if  r.mec[i] <> 0  then     { power of parameter not nullo }
      begin
        a:=a*powall[i][r.mec[i]]; { to multiply with power }
        case  i  of               { 1 , 2 , 3 , 4 , 9 , 13 }
          1 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[1]*amplip[j]; { d(p)/d(L) }
          2 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[2]*amplie[j]; { d(e)/d(L) }
          3 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[3]*amplis[j]; { d(s)/d(L) }
          4 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[4]*amplic[j]; { d(c)/d(L) }
          9 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[9]*amplia[j]; { d(a)/d(L) }
         13 : for j:=1 to 3 do                 { for L G H }
                d[j]:=d[j]+r.mec[i]*amplin[j]; { d(n)/d(L) }
        end; { case i = 1 2 3 4 9 13 parameter p e s c a n }
      end;
end;
{ to add d/d(L G H) to array ys[4,5,6] }
procedure difampo ( z  : extended ;   { sinus or cosinus }
                    a  : extended ;   { amplitude value }
                    d  : tampdif ;    { derivations value }
                    var  y  : tvarsix ) ; { d(r)/d(L G H) }
begin
  y[4]:=y[4]+a*d[1]*z; { to add with current term d/dL }
  y[5]:=y[5]+a*d[2]*z; { to add with current term d/dL }
  y[6]:=y[6]+a*d[3]*z; { to add with current term d/dL }
end;
{ angle E may be in argument
  to take into account derivations
     dE/dl dE/de*de/dL dE/de*de/dG
     dE/dl=1/(1-e*cos(E))  dE/de=sin(E)/(1-e*cos(E))
  values:  E argall[2]  e parall[2]  de/d(L G) e*amplie[1 2] }
procedure decan ( s , c  : extended ; { sunus cosinus argument }
                  a  : extended ;     { amplitude with powers }
                  r  : trowvar ;      { current term }
                  var  y  : tvarsix ) ; { to add derivations }
var
  e,u,q,p,z  : extended ; { temporary variables }
begin
  e:=parall[2]; { value e eccentricity global untypvar }
  u:=argall[2]; { value E eccentric anomaly in radian }
  q:=1.0/(1.0-e*cos(u)); { 1/(1-e*cos(E)) dE/dl partial }
  p:=sin(u)*q;      { sin(E)/(1-e*cos(E)) dE/de partial }
  if  r.fcs = 'C'   { for d(r)/dE our function is cosinus }
    then  z:=-r.mar[2]*a*s  { - sinus after differentiating  }
    else  z:=+r.mar[2]*a*c; { + cosinus for sinus d(r)/dE }
  y[1]:=y[1]+z*q; { d(r)/dl=d(r)/dE*dE/dl global untypvar }
  y[4]:=y[4]+z*p*e*amplie[1]; { d(r)/dE*dE/de*de/dL }
  y[5]:=y[5]+z*p*e*amplie[2]; { d(r)/dE*dE/de*de/dG }
end;
{ for simple derivations  d(r)/dg to y[2]  d(r)/dh to y[3] }
procedure difanglo ( s , c  : extended ; { sunus cosinus argument }
                     a  : extended ;     { amplitude with powers }
                     r  : trowvar ;      { current term }
                     var  y  : tvarsix ) ; { to add derivations }
var
  z  : extended ; { temporary variable }
begin
  if  r.fcs = 'C'   { for d(r)/d(g h) but our function is cosinus }
    then  z:=-a*s   { - sinus after differentiating  }
    else  z:=+a*c;  { + cosinus for sinus d(r)/d(g h) }
  y[2]:=y[2]+r.mar[3]*z;   { d(r)/dg global untypvar }
  y[3]:=y[3]+r.mar[4]*z;   { d(r)/dh global untypvar }
end;
{ for any row R to calculate
  \partial d(R)/d(l) to ys[1]
  \partial d(R)/d(g) to ys[2]
  \partial d(R)/d(h) to ys[3]
  \partial d(R)/d(L) to ys[4]
  \partial d(R)/d(G) to ys[5]
  \partial d(R)/d(H) to ys[6] }
procedure corelem ( xo  : tsagano ;        { L G H l g h canonical }
                    me  : word ;           { number of terms }
                    ro  : trowptr ;        { dynamical array for terms }
                    var  ys  : tvarsix ) ; { corrections to elements }
var                                        { or partial derivations }
  m  : word ;     { count for terms }
  r  : trowvar ;  { untypvar current term types untypvar }
  a  : extended ; { current amplitude value }
  d  : tampdif ;  { to sum partial derivations with L G H }
  s  : extended ; { sinus argument }
  c  : extended ; { cosinus argument }
begin
  { numerical values of partial derivations  with L G H }
  simpartial(xo.e); { untheled to global arrays amplia p e s c n }
  simparamon(xo.e); { p e s c a n  num.value to global parall unthelep }
  simargumon(xo.e); { numerical values for angles v E g h  unthelep }
  giparargum(xo.t); { parameter and angles of giant planet untheleg }
  allparapow(me,ro); { unthelpo our parameters to powers }
                 { to sum terms in our row for corrections }
  m:=0;          { for all terms }
repeat           { term  by term }
  m:=m+1;        { the next number }
  r:=ro^[m];     { for the next term }
  { to sum angles with coefficients to argument value } 
  cosargu(r,s,c); { term sinus cosinus argument }
  { numerical value of amplitude and derivations }
  difacto(r,a,d);   { term amplitude derivations }
  if  r.fcs = 'C'            { may be cosinus or sinus }
    then  difampo(c,a,d,ys)  { to add d/d(L G H) to ys }
    else  difampo(s,a,d,ys); { to add d/d(L G H) to ys }
  if  r.mar[2] <> 0          { there is angle E in argument }
    then  decan(s,c,a,r,ys); { dE/dl dE/de*de/dL dE/de*de/dG }
  difanglo(s,c,a,r,ys);      { d(r)/dg d(r)/dh }
until  m = me ;  { for all terms }
end;
end.