Операции над элементарными слагаемыми
В исходных текстах вычислительных процедур определим элементарное слагаемое в виде структуры:
{ for example case expansion with Bessel function }
unit unalltyp ; { for different types for some case }
interface
const
maxarr = 74967 ; // dimension of arrays
numsel = 8 ; // selected elements but only 8 in use
maxsel : array[1..numsel] of ShortInt = // maximum power
( 10 , // e in use
10 , // s=sin(i) in use
10 , // e' in use
10 , // s'=sin(i') in use
30 , // a semimajor axis in use
30 , // a' semimajor axis in use
8 , // n mean motion im use
8 ) ; // n' mean motion in us
maxarg = 6 ; // number of arguments l g h l' g' h'
type
TDimSel = array[1..numsel] of ShortInt ; { for power of var. 8 bytes }
TDimArg = array[1..maxarg] of ShortInt ; { for arguments 6 bytes }
{ elementary term has a form
A amplitude numeric value
*e^i_1*(sin i)^i_3
*e'^i_2*(sin i')^i_4
*a^i_5*a'^i_6
*cos(k_1*l+k_2*g+k_3*h+k_4*l'+k_5*g'+k_6*h') }
TRowVar = record
npa : Byte ; { the first planet 1 byte }
npb : Byte ; { the second planet 1 byte }
amp : Extended ; { amplitude simple }
mec : TDimSel ; { power of variables 8 bytes }
mar : TDimArg ; { coef. for argument 6 bytes }
fcs : Char ; { function sin or cos 1 byte }
end;
type
{ for the array }
TRowDim = array[1..maxarr] of TRowVar ;
TOrdNum = array[0..27] of longint ;
type
{ for the pointer to array }
TRowPtr = ^TRowDim ;
implementation
end.
Определим также несколько динамических массивов и выделим для них место в оперативной памяти компьютера:
{ for global variables }
unit unallvar ;
interface
uses
unalltyp ; { for all types }
var
maxsum : ShortInt ;{ maximum for sum of powers }
memrec : Word ; { size of record TRowVar UnAllTyp }
recnul : TRowVar ; { for nullo record type from UnAllTyp }
mecone : longint ; { number of terms in rowone array }
mectwo : longint ; { number of terms in rowtwo array }
mectre : longint ; { number of terms in rowtre array }
mecfur : longint ; { number of terms in rowfur array }
mecfiv : longint ; { number of terms in rowfiv array }
mecsix : longint ; { number of terms in rowsix array }
mecaux : longint ; { number of terms in rowaux array }
mectmp : longint ; { number of terms in tempor array }
meccoh : longint ; { number of terms in rowcoh array }
rowone : TRowPtr ; { to the pointer for array unalltyp }
rowtwo : TRowPtr ; { to the pointer for array }
rowtre : TRowPtr ; { to the pointer for array }
rowfur : TRowPtr ; { to the pointer for array }
rowfiv : TRowPtr ; { to the pointer for array }
rowsix : TRowPtr ; { to the pointer for array }
rowaux : TRowPtr ; { to the pointer for array }
rowtmp : TRowPtr ; { pointer for tempor array }
rowcoh : TRowPtr ; { pointer for rowcoh array }
implementation
end.
unit ungetmem ; { for dynamic memory }
interface
procedure GetRowMem ;
procedure FreeRowMem ;
implementation
uses
unalltyp , { for global types }
unallvar ; { for global variables }
procedure GetRowMem ; { array to heap by pointer }
begin
memrec:=SizeOf(TRowVar);
GetMem(rowone,SizeOf(TRowDim));
GetMem(rowtwo,SizeOf(TRowDim));
GetMem(rowtre,SizeOf(TRowDim));
GetMem(rowfur,SizeOf(TRowDim));
GetMem(rowfiv,SizeOf(TRowDim));
GetMem(rowsix,SizeOf(TRowDim));
GetMem(rowaux,SizeOf(TRowDim));
GetMem(rowtmp,SizeOf(TRowDim));
GetMem(rowcoh,SizeOf(TRowDim));
end;
procedure FreeRowMem ; { array to heap by pointer }
begin
FreeMem(rowcoh,SizeOf(TRowDim));
FreeMem(rowtmp,SizeOf(TRowDim));
FreeMem(rowaux,SizeOf(TRowDim));
FreeMem(rowsix,SizeOf(TRowDim));
FreeMem(rowfiv,SizeOf(TRowDim));
FreeMem(rowfur,SizeOf(TRowDim));
FreeMem(rowtre,SizeOf(TRowDim));
FreeMem(rowtwo,SizeOf(TRowDim));
FreeMem(rowone,SizeOf(TRowDim));
end;
end.
Выполним присвоение начальных значений:
unit UnRowIni ; { for initial data for rows }
interface
procedure RecNullo ; { without planets numbers }
procedure RowNullo ; { initial nullo for six arrays }
procedure RowTreNul ; { initial nullo for rowtre^ rowaux^ }
implementation
uses
UnAllTyp , { for type TRowVar var row one two tre aux coh }
UnAllVar ; { for global variables }
procedure RecNullo ; { without planets numbers }
var
j : Byte ;
begin
recnul.npa:=0;
recnul.npb:=0;
recnul.amp:=0.0; { amplitide nullo initial }
for j:=1 to numsel do { numsel = 4 e s e' s' }
recnul.mec[j]:=0; { for all power nullo }
for j:=1 to maxarg do
recnul.mar[j]:=0; { for all argument nullo }
recnul.fcs:='C'; { let it be cosinus but nullo amplitude }
end;
procedure RowNullo ; { initial nullo for six arrays }
var
i : longint ; { simple count }
begin
mecone:=0; { initial nullo terms in six array }
mectwo:=0;
mectre:=0;
mecfur:=0; { initial nullo terms in six array }
mecfiv:=0;
mecsix:=0;
mecaux:=0;
mectmp:=0;
meccoh:=0;
for i:=1 to maxarr do { const maxe from UnTypEcc }
begin
rowone^[i]:=recnul;
rowtwo^[i]:=recnul;
rowtre^[i]:=recnul;
rowfur^[i]:=recnul;
rowfiv^[i]:=recnul;
rowsix^[i]:=recnul;
rowaux^[i]:=recnul;
rowtmp^[i]:=recnul;
rowcoh^[i]:=recnul;
end;
end;
procedure RowTreNul ; { initial nullo values for rowtre array }
var
i : longint ; { simple count }
begin
mectre:=0;
mecaux:=0;
for i:=1 to maxarr do { const maxe from UnTypEcc }
begin
rowtre^[i]:=recnul;
rowaux^[i]:=recnul;
end;
end;
end.
Для элементарных слагаемых опеределены операции сложения, умножения и приведения подобных членов:
unit UnRowMul ; { to multiply two rows result to rowtre }
interface
uses
Windows , { for function GetTickCount }
unalltyp ; { for global types }
procedure AmpFactor ( c : Extended ; var a : extended ); // UnTypVar
procedure AmplMults ( a , b : extended ; var c : extended ) ; // UnTypVar
procedure AmplPower ( x : extended ; var y : extended ) ; { UnTypVar }
function BooMaxSum ( a , b : TRowVar ) : Boolean ;
procedure PowerParm ( a , b : TRowVar ; var c : TRowVar ) ;
function BooComPair ( x , e : TRowVar ) : Boolean ;
function BooCompare ( x , e : TRowVar ) : Boolean ;
procedure ToPositiveA ( var x : TRowVar ) ;
procedure ForSinOrCos ( x , y : TRowVar ;
var a , b : extended ;
var c : Char ) ;
procedure ToMultRow ( ua , ub : TRowVar ) ; { simple records }
procedure ForAmpNul ; { rowtre may has nullo amplitudes }
procedure MultForRow ( am : Extended ) ; { am * rowone * rowtwo = rowtre }
implementation
uses
UnAllVar , { for global variables }
UnRowLow , { for procedure LowerAux }
ungetida ; { for current date and time function strCurtida }
procedure AmpFactor ( c : Extended ; var a : extended );
begin
a:=c*a; { amplitude with the factor c }
end;
procedure AmplMults ( a , b : extended ; { two amplitudes }
var c : extended ) ; { result amplitude }
begin { to multiple two amplitudes }
c:=a*b; { only amplitudes }
end;
{ amplitude in power -1 : y = 1/x = x^(-1) }
procedure AmplPower ( x : extended ; var y : extended ) ;
begin
if ( Abs(x) < 1.0e-15 ) { too small }
then y:=0.0 { nullo simple }
else y:=1.0/x; { +1/x coefficients }
end;
function BooMaxSum ( a , b : TRowVar ) : Boolean ;
var
m : ShortInt ;
begin
m:=a.mec[1]+b.mec[1]+a.mec[2]+b.mec[2] // e s
+a.mec[3]+b.mec[3]+a.mec[4]+b.mec[4]; // e' s'
if m > maxsum // const maxsum from UnTypVar
then // sum of powers is too large
BooMaxSum:=False // exit no term to row
else
BooMaxSum:=True; // there is our term
end;
procedure PowerParm ( a , b : TRowVar ; var c : TRowVar ) ;
var // to analyse power of variables
i : Byte ;
begin
c.npa:=a.npa; { only for selected pair }
c.npb:=a.npb; { a.npa npb = b.npa npb }
for i:=1 to numsel do { for each parameter numsel = 8 e s e' s' }
c.mec[i]:=a.mec[i]+b.mec[i]; { to sum of two powers a a' n n' }
end;
function BooComPair ( x , e : TRowVar ) : Boolean ;
begin
BooComPair:= ( ( x.npa = e.npa ) and ( x.npb = e.npb ) ) ;
end;
function BooCompare ( x , e : TRowVar ) : Boolean ;
var
j : Byte ;
b : Boolean ;
begin
b:= ( x.fcs = e.fcs ) ; { characters for functions 'C' or 'S' }
if b
then { to compare numbers for two planets }
b:=BooComPair(x,e);
j:=0;
while ( b and ( j < numsel ) ) do
begin
j:=j+1;
b:= ( x.mec[j] = e.mec[j] ) ; { the powers of variables }
end;
j:=0;
while ( b and ( j < maxarg ) ) do
begin
j:=j+1;
b:= ( x.mar[j] = e.mar[j] ) ; { the arguments of variables }
end;
BooCompare:=b; { true if it is similar term or false if it is not }
end;
procedure TheNextTerm ( x : TRowVar ) ; { to find similar if it exists }
var
m : longint ;
b : Boolean ;
e : TRowVar ; { current record for compare }
begin { to find similar terms for current term in array rowtre^[mectre] }
b:=False;
m:=0; { nullo }
while ( ( not b ) and ( m < mectre ) ) do
begin
m:=m+1; { to try the next term }
e:=rowtre^[m]; { current record for compare }
b:=BooCompare(x,e);
end;
if b { there is similar term }
then { to sum only amplitudes }
rowtre^[m].amp:=rowtre^[m].amp+x.amp
else { there is no similar term in current array }
begin { it is our new term in result array }
if mectre = maxarr then Exit ; { too many terms }
mectre:=mectre+1;
rowtre^[mectre]:=x;
end;
end;
procedure ToPositiveA ( var x : TRowVar ) ;
var
i,j : Byte ;
begin { to positive for the first nonnullo argument }
for j:=1 to maxarg do
begin
if x.mar[j] > 0
then
Exit ; { well the argument is positive }
if x.mar[j] < 0 { the argument may be negative }
then
begin
for i:=j to maxarg do { to positive for the current argument }
x.mar[i]:=-x.mar[i];
if x.fcs = 'S' { for sinus function minus for amplitude }
then
AmpFactor(-1.0,x.amp); // minus for amplitude
Exit ; { well the argument is positive }
end;
end;
end;
procedure ForSinOrCos ( x , y : TRowVar ;
var a , b : extended ;
var c : Char ) ;
begin { 2*cos(x)*sin(y) =-sin(x-y) + sin(x+y) }
b:=a; { amlpitude may be minus }
if x.fcs = y.fcs { 2*cos(x)*cos(y) = cos(x-y) + cos(x+y) }
then { cos*cos or sin*sin }
begin { +cos(x-y) }
c:='C'; { cosinus function as result }
if x.fcs = 'S' { sin*sin case }
then { 2*sin(x)*sin(y) = cos(x-y) - cos(x+y) }
AmpFactor(-1.0,b); { -cos(x+y) }
end { 2*sin(x)*cos(y) = sin(x-y) + sin(x+y) }
else { sin*cos or cos*sin }
begin { +sin(x+y) }
c:='S'; { sinus function as result }
if x.fcs = 'C' { cos*sin case }
then { 2*cos(x)*sin(y) =-sin(x-y) + sin(x+y) }
AmpFactor(-1.0,a); { -sin(x-y) }
end;
end;
procedure ToMultRow ( ua , ub : TRowVar ) ; { simple records }
var { 2*cos(x)*cos(y) = cos(x-y) + cos(x+y) }
j : Byte ; { 2*sin(x)*sin(y) = cos(x-y) - cos(x+y) }
b : extended ; { 2*sin(x)*cos(y) = sin(x-y) + sin(x+y) }
x : TRowVar ; { 2*cos(x)*sin(y) =-sin(x-y) + sin(x+y) }
begin { from two record may be two new records }
AmplMults(ua.amp,ub.amp,x.amp); { to multiply amplitudes }
if ( ( Abs(x.amp) < 1.0e-19 ) { may be nearly nullo term }
or ( not BooMaxSum(ua,ub) ) )
then Exit ; { sum of powers is too large }
AmpFactor(0.5,x.amp); { factor 0.5 for amplitude }
PowerParm(ua,ub,x); { to sum all powers }
ForSinOrCos(ua,ub,x.amp,b,x.fcs); { amplitude function indentify }
for j:=1 to maxarg do
x.mar[j]:=ua.mar[j]-ub.mar[j]; { for the first all arguments minus }
ToPositiveA(x);
TheNextTerm(x); { for minus in argument x-y }
x.amp:=b; { for plus in argument x+y }
for j:=1 to maxarg do
x.mar[j]:=ua.mar[j]+ub.mar[j]; { for the second all arguments plus }
ToPositiveA(x);
TheNextTerm(x);
end;
procedure ForAmpNul ; { rowtre may has nullo amplitudes }
var
m : longint ;
begin
mecaux:=0;
for m:=1 to mectre do
if ( Abs(rowtre^[m].amp) > 1.0e-17 )
then
begin
mecaux:=mecaux+1; { non nullo term }
rowaux^[mecaux]:=rowtre^[m];
end;
if mecaux < mectre
then
LowerAux(3) ; { rowaux to rowtre from UnRowLow }
end;
{ to multiply two rows in arrays
for current date and time there is
function strCurtida of shortstring from ungetida }
procedure MultForRow ( am : Extended ) ; { am * rowone * rowtwo = rowtre }
var { with the factor am }
ma,mb : longint ;
ua,ub : TRowVar ; { type from UnTypVar }
la,lb : longint ; { current time in millisecond }
begin
la:=GetTickCount ; { function from unit Windows in millisecond }
for ma:=1 to mecone do
begin
ua:=rowone^[ma]; { current record from the first row }
AmpFactor(am,ua.amp); { to multiply by amplitude }
for mb:=1 to mectwo do
begin
ub:=rowtwo^[mb]; { current record from the second row }
if ( ( ua.npa = ub.npa ) and ( ua.npb = ub.npb ) )
then { result to rowtre with the count mectre }
ToMultRow(ua,ub); { only pair numbers na nb }
lb:=GetTickCount ; { current time value in millisecond }
if ( (lb-la) > 30000 ) then { too long interval }
begin { in millisecond }
writeln(' ',ma,' * ',mb,' = ',mectre,' term by term');
la:=lb; { to replace the moment }
end;
end;
end;
ForAmpNul ; { in current array rowtre may be nullo amplitude }
WriteLn('to multiply ',
mecone,' terms * ',mectwo,' terms = ',mectre,' terms',
' power = ',maxsum);
end;
end.
Здесь можно вернуться к возмущающей функции
Здесь можно вернуться на страницу Программные приложения