File : rg.adb
with Rounding, Ints.IO, Messages;
use Rounding, Ints.IO, Messages;
package body RG is
use RG_Funs, RG_Taylors, RG_Mini, RGV, RGV_ScalVecs, RG_Funs_Num, Num_Funs;
----------------- first conversions Numeric <--> Scalar
function N2S(N: Numeric) return Scalar is
begin
return Scal(N2R(N));
end N2S;
function S2N(S: Scalar) return Numeric is
begin
return R2N(Approx(S));
end S2N;
procedure N2S(F: in NTay; G: out STay) is
use NFT;
begin
G.R := F.R;
for I in Power loop
for J in Power loop
G.C(I,J) := N2S(F.C(I,J));
end loop;
end loop;
end N2S;
procedure S2N(F: in STay; G: out NTay) is
begin
G.R := F.R;
for I in Power loop
for J in Power loop
G.C(I,J) := S2N(F.C(I,J));
end loop;
end loop;
end S2N;
procedure N2S(F: in NFun; G: out SFun) is
begin
G.A := F.A;
G.E := True;
N2S(F.P,G.P);
N2S(F.Q,G.Q);
end N2S;
function N2S(F: NFun) return SFun is
G: SFun;
begin
N2S(F,G);
return G;
end N2S;
procedure S2N(F: in SFun; G: out NFun) is
begin
G.A := F.A;
G.E := True;
S2N(F.P,G.P);
S2N(F.Q,G.Q);
end S2N;
function S2N(F: SFun) return NFun is
G: NFun;
begin
S2N(F,G);
return G;
end S2N;
------------------ some procedures that use numeric approximations
procedure NumRead(Name: in String; F: out SFun; Decimal: in Boolean := False) is
--- just truncates powers higher than PDeg
begin
if Trunc then Read(Name,F,Decimal); return; end if;
declare
use NFF;
G: NFun;
begin
Relax_Checks;
Read(Name,G,Decimal);
N2S(G,F);
Restore_Checks;
end;
end NumRead;
procedure NumChangeVar(A: in Args; F: in out SFun) is
--- change variables from F.A to A
begin
if Trunc then ChangeVar(A,F); return; end if;
declare
use NFF;
G: NFun;
begin
Relax_Checks;
S2N(F,G);
ChangeVar(A,G);
N2S(G,F);
Restore_Checks;
end;
end NumChangeVar;
procedure NumRVal(G: in SFun; X,C: in Scalar; S: in out Scalar) is
--- approximation for RValue
use NFF;
N: Numeric;
begin
Relax_Checks;
N := S2N(S);
RVal(S2N(G),S2N(X),S2N(C),N);
S := N2S(N);
Restore_Checks;
end NumRVal;
procedure RValue(G: in SFun; X,C: in Scalar; S: in out Scalar) is
--- solves G(x,s)=c
begin
if not Trunc then
NumRVal(G,X,C,S);
end if;
RVal(G,X,C,S);
end RValue;
procedure NumLVal(F: in SFun; Y,C: in Scalar; S: in out Scalar) is
--- approximation for LValue
use NFF;
N: Numeric;
begin
Relax_Checks;
N := S2N(S);
LVal(S2N(F),S2N(Y),S2N(C),N);
S := N2S(N);
Restore_Checks;
end NumLVal;
procedure LValue(F: in SFun; Y,C: in Scalar; S: in out Scalar) is
--- solves F(s,y)=c
begin
if not Trunc then
NumLVal(F,Y,C,S);
end if;
LVal(F,Y,C,S);
end LValue;
procedure ZValue(F,G: in SFun; X,Y: in Scalar; S: in out Scalar) is
--- computes midpoint z: F(x,z)+G(z,y)=0
use NFF;
N: Numeric;
begin
if not Trunc then
Relax_Checks;
N := S2N(S);
ZVal(S2N(F),S2N(G),S2N(X),S2N(Y),N);
S := N2S(N);
Restore_Checks;
end if;
ZVal(F,G,X,Y,S);
end ZValue;
procedure VValue(F,G: in SFun; X: in Scalar; S: in out Scalar) is
--- computes midpoint s: G(x,s)+F(s,-s)=0
use NFF;
N: Numeric;
begin
if not Trunc then
Relax_Checks;
N := S2N(S);
VVal(S2N(F),S2N(G),S2N(X),N);
S := N2S(N);
Restore_Checks;
end if;
VVal(F,G,X,S);
end VValue;
procedure YValue(G: in SFun; S: in out Scalar) is
--- used to compute Lambda: Gx(s^2,s)=[Gy(s^2,s)]^2
use NFF;
N: Numeric;
begin
if not Trunc then
Relax_Checks;
N := S2N(S);
YVal(S2N(G),N);
S := N2S(N);
Restore_Checks;
end if;
YVal(G,S);
end YValue;
---------------------------------------
procedure MakePackIndex(Verbose: Boolean := False) is
--- enumerate low modes, independently of PDeg
I,K: Integer := 0;
begin
for L in Power loop JH(L) := 0; end loop;
for D in 1 .. PDeg loop
for J in 0 .. D loop
I := D-J;
JH(I) := IMax(J+1,JH(I));
K := K+1;
II(K) := I;
JJ(K) := J;
if Verbose then
Show("I J K " & Strng(I,4) & Strng(J,4) & Strng(K,6));
end if;
if K=SVDim then return; end if;
end loop;
end loop;
end MakePackIndex;
procedure ShowPackIndex is
begin
MakePackIndex(True);
end ShowPackIndex;
procedure Pack(F: in STay; V: out SVec) is
--- copy low coefficients from F into V
C: Scalar;
begin
for K in SVRange loop
C := F.C(II(K),JJ(K));
if not IsNumbers(C) then
Error("RG.Pack error: not implemented");
end if;
V(K) := C;
end loop;
end Pack;
procedure UnPack(V: in SVec; F: out STay) is
begin
F := TayZero;
for K in reverse SVRange loop
F.C(II(K),JJ(K)) := V(K);
end loop;
end UnPack;
procedure Prod(A: in SMat; F: in STay; G: out STay) is
--- rank VDim operator
V: SVec;
begin
Pack(F,V);
Mult(A,V);
UnPack(V,G);
G.R := F.R;
end Prod;
function Norm(R: RadPair; A: SMat) return Scalar is
--- operator norm
E: Scalar := ScalZero;
F,G: STay;
begin
for K in SVRange loop
UnitVector(II(K),JJ(K),R,F);
Prod(A,F,G);
E := Max(Norm(R,G),E);
end loop;
return E;
end Norm;
---------------------------------------
function Inv(F: SFun) return SFun is
--- first 1/F numerically, then error bound
begin
if Trunc then return NumInv(F); end if;
declare
use NFN;
N: constant Integer := 256;
Del,Eps,Rad,K: Scalar;
H,U: SFun;
begin
Relax_Checks;
H := N2S(NumInv(S2N(F)));
Restore_Checks;
Prod(F,H,U);
AddConst(-ScalOne,U);
Del := Norm(U);
Eps := Norm(H)*Del;
Rad := Rep(N+1)*Eps/Rep(N);
AddBall(Rad,H);
K := Two*(Del+Rad*Norm(F));
if Sup(Rep(N+1)*K)>One then
Error("RG.Inv error: (N+1)*K>1",Sup(Rep(N+1)*K),T3);
end if;
return H;
end;
end Inv;
function "/"(F,G: SFun) return SFun is
begin
return F*Inv(G);
end "/";
function Inv(L: Flop) return Flop is
--- inverse of operator L
H: SFun;
M: Flop;
begin
H := Inv(L(1)*Reflect(L(1))-L(2)*Reflect(L(2)));
M(1) := H*Reflect(L(1));
M(2) := -H*L(2);
return M;
end Inv;
--------------------------- now RG related procedures
procedure BasicData(G: in SFun) is
--- show some info, V can be zero
GX,GY,La,Mu: Scalar;
begin
ValDer(ScalOne,ScalZero,G,GX,GY);
LambdaMu(G,La,Mu);
Show1("G:",G.A);
Show("G: Rho = ",Rho1(G),Rho2(G));
Show1("G: Gx(1,0) = ",GX);
Show1("G: Gy(1,0) = ",GY);
Show1("G: Lambda = ",La);
Show1("G: Mu = ",Mu);
New_Line;
end BasicData;
procedure CompInfo(G,V: in SFun) is
--- additional info, was used for determining appropriate arguments
--- V can be zero
R: RadPair;
La,Mu: Scalar;
Vt: SFun := V;
F: SFun;
begin
MakeLaMuF(G,La,Mu,F);
GuessV(F,G,Vt);
R := CompYRadii(Vt,G.A);
RhoInfo("G:",G,R);
R := CompVWRadii(Vt,F.A);
RhoInfo("F:",F,R);
end CompInfo;
procedure PreNormalize(R: in RadPair; G: in out SFun) is
--- normalization via scaling
Rt: RadPair;
S,St: Scalar;
Gx,Gy: SFun;
begin
Symm(G);
Rt := Reps.Between(ValRadii(ScalOne,ScalZero,G.A),Rho(G));
Der(Rt,G,Gx,Gy);
S := ScalOne;
LValue(Gy,ScalZero,ScalZero,S);
Scale(S,G.A,G);
ValDer(ScalOne,ScalZero,G,S,St);
Mult(Inv(S),G);
SetRho(R,G);
ZapConst(G.P);
end PreNormalize;
procedure Normalize(S1,S2: in Scalar; G: in out SFun) is
--- ad hoc normalization, avoid scaling for accuracy
GX,GY,A11,A12,A21,A22,D: Scalar;
begin
Symm(G);
ValDer(ScalOne,ScalZero,G,GX,GY);
GX := GX-S1;
GY := GY-S2;
A11 := ScalTwo+Scal(G.A.B);
A21 := ScalTwo-Scal(G.A.B);
A12 := Two*Scal(G.A.A)+ScalOne;
A22 := Two*Scal(G.A.A)-ScalOne;
D := Four*(ScalOne-G.A.A*Scal(G.A.B));
Add((A22*GX-A12*GY)/D,G.P.C(1,0));
Add((A11*GY-A21*GX)/D,G.P.C(0,1));
ZapConst(G.P);
end Normalize;
procedure Normalize(G: in out SFun) is
begin
Normalize(ScalOne,ScalZero,G);
end Normalize;
procedure LambdaMu(G: in SFun; La,Mu: out Scalar) is
--- compute Lambda and Mu
--- good accuracy of La comes from large 2nd derivatives
GX,GY: Scalar;
begin
La := Scal(La0);
YValue(G,La);
ValDer(Sqr(La),La,G,GX,GY);
Mu := -GY;
end LambdaMu;
procedure Scale(G: in SFun; La,Mu: in Scalar; F: in out SFun) is
A: constant Args := F.A; --- needs to be defined
begin
F := G;
Scale(La,A,F);
ZapConst(F.P);
Mult(Inv(La*Mu),F);
end Scale;
procedure MakeLaMuF(G: in SFun; La,Mu: out Scalar; F: out SFun) is
--- compute La,Mu,F from G
begin
LambdaMu(G,La,Mu);
SetZero(FArgs,F);
Scale(G,La,Mu,F);
end MakeLaMuF;
function VArgs(A: Args) return Args is
--- arguments for midpoint function V
L2: constant Rep := 0.4995601179710913; -- approximate La^2
VA: Args;
begin
VA.T0 := A.T0*L2;
VA.S0 := A.S0*L2;
VA.B := A.B*L2;
VA.A := A.A/L2;
return VA;
end VArgs;
function VRadii(R: RadPair) return RadPair is
--- radii for midpoint function V
L2: constant Rep := 0.4996; -- slightly larger than La^2
begin
return (L2*R(1),L2*R(2));
end VRadii;
procedure MidPointGuess(F,G: in SFun; R: in RadPair; V: in out SFun) is
--- approximate fixed point function V
begin
NumChangeVar(VArgs(G.A),V);
if Trunc then NumCompZeroD(F,G,R,V); return; end if;
declare
use NFN;
NV: NFun;
begin
Relax_Checks;
S2N(V,NV);
NumCompZeroD(S2N(F),S2N(G),R,NV);
N2S(NV,V);
Restore_Checks;
end;
end MidPointGuess;
procedure GuessV(F,G: in SFun; V: in out SFun) is
--- approximate fixed point function V
Rv: constant RadPair := VRadii(Rho(G));
begin
NumChangeVar(VArgs(G.A),V);
if IsZero(V) then
MidPointGuess(F,G,Rv,V);
end if;
SetRho(Rv,V);
end GuessV;
procedure FGDers(F,G,V: in SFun; Fx,Gy,Fxx,Fxy,Gyy: out SFun) is
--- derivatives needed in MidPoint
RF: constant RadPair := CompVWRadii(V,F.A);
RG: constant RadPair := CompYRadii(V,G.A);
Hx,Hy,Hxx,Hxy,Hyy: SFun;
begin
Der(RF,F,Fx,Hy,Fxx,Fxy,Hyy);
Der(RG,G,Hx,Gy,Hxx,Hxy,Gyy);
end FGDers;
procedure MidPoint(F,G: in SFun; Rv: in RadPair; V: in out SFun) is
--- solves Gy(x,V)+Fx(V,W)=0
--- initial V can be zero
N: constant Integer := 256;
Eps,Rad,K: Scalar;
Fx,Fxx,Fxy,Gy,Gyy: SFun;
L,Li: Flop;
begin
TraceEnter("RG.MidPoint");
MidPointGuess(F,G,Rv,V);
if Trunc then TraceLeave; return; end if;
FGDers(F,G,V,Fx,Gy,Fxx,Fxy,Gyy);
L := DFF(Fxx,Fxy,Gyy,V);
Li := Inv(L);
Eps := Norm(Rv,Li*FF(Fx,Gy,V));
Rad := Rep(N+1)*Eps/Rep(N);
AddBall(Rad,V);
FGDers(F,G,V,Fx,Gy,Fxx,Fxy,Gyy);
K := Norm(Rv,Li*(DFF(Fxx,Fxy,Gyy,V)-L));
if Sup(Rep(N+1)*K)>One then
Error("RG.MidPoint error: (N+1)*K>1",Sup(Rep(N+1)*K));
end if;
TraceLeave;
end MidPoint;
procedure Compose(F,G,V: in SFun; H: out SFun) is
--- assuming F and G are symmetric
begin
H := ScalTwo*Symm(CompY(V,G))+Symm(CompVW(V,F));
end Compose;
procedure Renorm(R: in RadPair; G,V: in out SFun; La,Mu: out Scalar; F: out SFun) is
--- G is assumed to be normalized
--- initial V can be zero
Rv: constant RadPair := VRadii(R);
H: SFun;
begin
TraceEnter("RG.Renorm");
MakeLaMuF(G,La,Mu,F);
MidPoint(F,G,Rv,V);
Compose(F,G,V,H);
Scale(H,Sqr(La),Sqr(Mu),G); -- now G should have (at least) domain R
SetRho(R,G);
TraceLeave;
end Renorm;
procedure Renorm(R: in RadPair; G,V: in out SFun) is
La,Mu: Scalar;
F: SFun;
begin
Renorm(R,G,V,La,Mu,F);
end Renorm;
procedure Renorm(R: in RadPair; G: in out SFun) is
V: SFun;
begin
SetZero(VArgs(G.A),V);
Renorm(R,G,V);
end Renorm;
procedure Renorm(Rt: in RadPair; G: in SFun; B: out RGData) is
--- need Rt larger than Rho(G) for DRenorm
begin
B.G := G;
B.Gt := G;
SetZero(VArgs(G.A),B.V);
Renorm(Rt,B.Gt,B.V,B.La,B.Mu,B.F);
end Renorm;
-------------------------------
procedure Write(Name: in String; B: in RGData; Decimal: in Boolean := False) is
F: File_Type;
V: SFun := B.V;
begin
Show("Writing " & Name,1);
Create(F,Out_File,Name);
Put(F,B.La,Decimal);
Put(F,B.Mu,Decimal);
Put(F,B.G,Decimal);
Put(F,B.F,Decimal);
Put(F,V,Decimal);
Put(F,B.Gt,Decimal);
Close(F);
end Write;
procedure Read(Name: in String; B: out RGData; Decimal: in Boolean := False) is
F: File_Type;
begin
Show("Reading " & Name,1);
Open(F,In_File,Name);
Get(F,B.La,Decimal);
Get(F,B.Mu,Decimal);
Get(F,B.G,Decimal);
Get(F,B.F,Decimal);
Get(F,B.V,Decimal);
Get(F,B.Gt,Decimal);
Close(F);
end Read;
procedure DNormalize(G: in out SFun) is
--- derivative of Normalize
begin
Normalize(ScalZero,ScalZero,G);
end DNormalize;
procedure DLambdaMu(G,DG: in SFun; La: in Scalar; DLa,DMu: out Scalar) is
--- derivative of LambdaMu
La2: constant Scalar := Sqr(La);
DGX,DGY,GX,GY,GXX,GXY,GYY: Scalar;
begin
ValDer(La2,La,DG,DGX,DGY);
ValDer(La2,La,G,GX,GY,GXX,GXY,GYY);
DLa := -(DGX-Two*GY*DGY)/(Two*La*GXX+GXY-Two*GY*(Two*La*GXY+GYY));
DMu := -DGY-(Two*La*GXY+GYY)*DLa;
end DLambdaMu;
procedure DScale(F,DG: in SFun; La,Mu,QLa,QMu: in Scalar; R: in RadPair; DF: out SFun) is
--- derivative of Scale, QLa is DLa/La, QMu is DMu/Mu
H: SFun;
begin
SetZero(F.A,DF);
Scale(DG,La,Mu,DF);
MultAdd(-QLa-QMu,F,DF);
MDer(R,F,H);
MultAdd(QLa,H,DF);
ZapConst(DF.P);
end DScale;
procedure DRenorm(B: in RGData; DG: in out SFun) is
--- derivative of Renorm
--- G is assumed to be normalized
RDG: constant RadPair := Rho(B.G);
RDF: RadPair;
DLa,DMu: Scalar;
DF,DH: SFun;
begin
TraceEnter("RG.DRenorm");
if IsZero(DG) then TraceLeave; return; end if;
DNormalize(DG);
DLambdaMu(B.G,DG,B.La,DLa,DMu);
SetZero(FArgs,DF);
RDF := Reps.Between(CompVWRadii(B.V,B.F.A),Rho(B.F));
DScale(B.F,DG,B.La,B.Mu,DLa/B.La,DMu/B.Mu,RDF,DF); --- F,DG -> DF
Compose(DF,DG,B.V,DH);
DScale(B.Gt,DH,Sqr(B.La),Sqr(B.Mu),Two*DLa/B.La,Two*DMu/B.Mu,RDG,DG); --- Gt,DH -> DG
TraceLeave;
end DRenorm;
-------------------------------
function Pack(F: SFun) return SVec is
--- copy low coefficients from F.P into V
V: SVec;
begin
Pack(F.P,V);
return V;
end Pack;
function UnPack(D: Domain; V: SVec) return SFun is
--- copy V into low coefficients of F.P
F: SFun;
begin
F.A := D.A;
UnPack(V,F.P);
F.E := IsNumbers(F.P); -- should always be true
F.Q := TayZero;
SetRho(D.R,F);
return F;
end UnPack;
procedure Sharpen(A: in out SMat) is
begin
if Trunc then return; end if;
for I in SVRange loop
for J in SVRange loop
A(I,J) := Scal(Approx(A(I,J)));
end loop;
end loop;
end Sharpen;
procedure Normalize(A: in Args; M: in out SMat; R: RadPair := Some_RadPair) is
D: constant Domain := (A,R);
H0,H: SFun;
begin
for K in SVRange loop
H0 := UnPack(D,GetColumn(K,M));
H := H0;
DNormalize(H);
if R /= Some_RadPair then
Show(Strng(K,2) & " ||H-H0|| =",Sup(Norm(R,H-H0)),1);
end if;
SetColumn(K,Pack(H),M);
end loop;
end Normalize;
function NumDRenorm(B: RGData) return SMat is
--- approx derivative of Renorm
U: SVec;
A: SMat;
DG: SFun;
begin
TraceEnter("RG.NumDRenorm");
if not Trunc then
Show("RG.NumDRenorm: run this in Numeric mode");
end if;
for K in SVRange loop
Show1("column" & Strng(K,4),1);
SetZero(U);
U(K) := ScalOne;
DG := UnPack(Dom(B.G),U);
DRenorm(B,DG);
DNormalize(DG);
SetColumn(K,Pack(DG),A);
end loop;
Sharpen(A); -- just in case ...
TraceLeave;
return A;
end NumDRenorm;
function ContractionMatrix(B: RGData) return SMat is
--- compute matrix used in contraction
--- should be run in Numeric mode
use RGV_Ops;
A: SMat;
begin
A := NumDRenorm(B);
AddConst(-One,A);
Invert(A);
AddConst(One,A);
Neg(A);
Sharpen(A); -- just in case ...
return A;
end ContractionMatrix;
function ContractionMatrix(G: SFun) return SMat is
Fac: constant Rep := Rep(15)/Rep(16);
R: constant RadPair := (Fac*Rho1(G),Fac*Rho2(G));
B: RGData;
begin
Renorm(R,G,B);
return ContractionMatrix(B);
end ContractionMatrix;
procedure Prod(A: in SMat; H: in SFun; G: out SFun) is
--- rank SVDim operator
--- avoid if H has low order error terms
begin
SetZero(H,G);
Prod(A,H.P,G.P);
SetRho(Rho(H),G);
end Prod;
procedure MProd(A: in SMat; H: in SFun; G: out SFun) is
--- avoid if H has low order error terms
begin
Prod(A,H,G);
Add(H,G);
end MProd;
procedure Contract(R: in RadPair; A: in SMat; G0: in SFun; V,H: in out SFun) is
--- the contraction
--- H can be zero but need valid Args
AH,G: SFun;
begin
TraceEnter("RG.Contract");
Prod(A,H,AH);
Sum(H,AH,G);
Add(G0,G);
Normalize(G);
Renorm(R,G,V);
Sub(G0,G);
Diff(G,AH,H);
TraceLeave;
end Contract;
procedure DContract(B: in RGData; A: in SMat; DH: in out SFun) is
--- derivative of Contract
ADH: SFun;
begin
TraceEnter("RG.DContract");
Prod(A,DH,ADH);
Add(ADH,DH);
DRenorm(B,DH);
Sub(ADH,DH);
TraceLeave;
end DContract;
function DContractNorm(B: RGData; A: SMat) return Scalar is
--- bound on the norm of DContract on the set B.G
MaxI: constant Integer := IMin(30,PDeg);
MaxJ: Integer := MaxI;
Eps: constant Radius := 0.1;
R: constant RadPair := Rho(B.G);
I: Integer := 0;
J: Integer := 1;
E,Ei,Eij: Radius := Zero;
DH: SFun;
begin
TraceEnter("RG.DContractNorm");
DH.A := B.G.A;
loop
UnitVector(I,J,R,DH.P);
DH.Q := TayZero;
if I >= MaxI then MaxJ := 0; end if;
if J >= MaxJ and not Trunc then ToBalls(DH); end if;
DContract(B,A,DH);
Eij := Sup(Norm(R,DH));
Show("norm(" & Strng(I,2) & "," & Strng(J,2) & ") <=",Eij,1);
Ei := RMax(Ei,Eij);
if J >= MaxJ then
E := RMax(E,Ei);
exit when MaxJ=0;
if JH(I)=0 and then Ei<Eps then MaxJ := 0; end if;
I := I+1;
J := 0;
MaxJ := IMax(JH(I),MaxJ);
Ei := Zero;
else
J := J+1;
if J>JH(I) and then Eij<Eps then MaxJ := J; end if;
end if;
end loop;
Show("norm(**,**) <=",E,1);
TraceLeave;
return Scal(E);
end DContractNorm;
---------------------------
procedure MakeMapData(A: in RGData; R: in Rectangle; B: out MapData) is
begin
B.G := A.G;
B.La := A.La;
B.Mu := A.Mu;
B.R := R;
end MakeMapData;
function DG0(G: SFun; X,Y: Scalar) return MapDer is
--- derivative of G0
GX,GY,GXX,GXY,GYY: Scalar;
D: MapDer;
begin
ValDer(X,Y,G,GX,GY,GXX,GXY,GYY);
D(1,1) := -GXX/GXY;
D(1,2) := -ScalOne/GXY;
D(2,1) := GXY-GXX*GYY/GXY;
D(2,2) := -GYY/GXY;
return D;
end DG0;
function DG0(B: MapData; X,Y: Scalar) return MapDer is
begin
return DG0(B.G,X,Y);
end DG0;
function DF0(B: MapData; X,Y: Scalar) return MapDer is
D: MapDer;
begin
D := DG0(B.G,X*B.La,Y*B.La);
D(1,2) := D(1,2)*B.Mu/B.La;
D(2,1) := D(2,1)*B.La/B.Mu;
return D;
end DF0;
function DH0(B: MapData; X,Y: Scalar) return MapDer is
D: MapDer;
begin
D := DG0(B.G,X/B.La,Y/B.La);
D(1,2) := D(1,2)*B.La/B.Mu;
D(2,1) := D(2,1)*B.Mu/B.La;
return D;
end DH0;
procedure CheckJ(B: in MapData) is
--- test used to show that J=Id
DF,DG,DH,DJ: MapDer;
begin
DH := DH0(B,B.La,ScalZero);
DF := DF0(B,B.La,ScalOne);
DG := DG0(B,ScalOne,ScalZero);
DJ := DG*DF*Inv2(DH);
Show1("DJ",DJ);
Check(DJ(1,1)>ScalZero);
end CheckJ;
function DMap(G: SFun; X,Y: Scalar) return MapDer is
--- derivative of Map(G)
GX,GY,GXX,GXY,GYY: Scalar;
D: MapDer;
begin
ValDer(X,Y,G,GX,GY,GXX,GXY,GYY);
D(1,1) := -GXX/GXY;
D(1,2) := -ScalOne/GXY;
D(2,1) := GXY-GXX*GYY/GXY;
D(2,2) := -GYY/GXY;
return D;
end DMap;
procedure CheckJ(G: in SFun) is
--- test used to show that J1=Id
La,Mu: Scalar;
DF,DG,DH,DJ: MapDer;
F,H: SFun;
begin
LambdaMu(G,La,Mu);
F := G;
Scale(La,F);
Mult(Inv(La*Mu),F);
H := G;
Scale(Inv(La),H);
Mult(La*Mu,H);
DH := DMap(H,La,ScalZero);
DF := DMap(F,La,ScalOne);
DG := DMap(G,ScalOne,ScalZero);
DJ := DG*DF*Inv2(DH);
Show1("DJ",DJ);
Check(DJ(1,1)>ScalZero);
end CheckJ;
--------------------------
procedure CheckRectangle(R: in Rectangle) is
begin
if R(1)>R(2) or R(3)>R(4) then
Error("RG.CheckRectangle error");
end if;
end CheckRectangle;
function WellInside(P: Point; R: Rectangle) return Boolean is
X: constant Scalar := P(1)+DiskAt0(Small); -- leave room for closure
Z: constant Scalar := P(2)+DiskAt0(Small);
begin
return Inf(X)>R(1) and Sup(X)<R(2) and Inf(Z)>R(3) and Sup(Z)<R(4);
end WellInside;
procedure Scale(SX,SZ: in Scalar; R: in Rectangle; Ri,Ro: out Rectangle) is
--- usable only for SX,SZ negative
TX: constant Scalar := SX+DiskAt0(Small);
TZ: constant Scalar := SZ+DiskAt0(Small);
begin
Ri(1) := Sup(R(2)*TX);
Ri(2) := Inf(R(1)*TX);
Ri(3) := Sup(R(4)*TZ);
Ri(4) := Inf(R(3)*TZ);
CheckRectangle(Ri);
Ro(1) := Inf(R(2)*TX);
Ro(2) := Sup(R(1)*TX);
Ro(3) := Inf(R(4)*TZ);
Ro(4) := Sup(R(3)*TZ);
CheckRectangle(Ro);
end Scale;
function Lambda(B: in MapData; P: Point) return Point is
begin
return (P(1)*B.La,P(2)*B.Mu);
end Lambda;
function InvLambda(B: in MapData; P: Point) return Point is
begin
return (P(1)/B.La,P(2)/B.Mu);
end InvLambda;
function GMap0(G: in SFun; P: Point) return Point is
--- image (y,w) of (x,z) under map G: -Gx(x,y)=z, w=Gy(x,y)
X: constant Scalar := P(1);
Z: constant Scalar := P(2);
R: RadPair := (0.1,0.1);
Y,T,W: Scalar := ScalZero;
Gx,Gy,Gxx,Gxy,Gyy: SFun; -- make global if need efficiency
begin
if not Trunc then
Der(R,G,Gx,Gy);
NumRVal(Gx,X,-Z,Y);
R := Reps.Between(ValRadii(X,Y,G.A),Rho(G));
end if;
Der(R,G,Gx,Gy,Gxx,Gxy,Gyy);
RVal(Gx,X,-Z,Y);
T := Val(X,Y,Gxy);
if Inf(T) <= One then
Error("RG.GMap0 error: twist not as promised");
end if;
W := Val(X,Y,Gy);
return (Y,W);
end GMap0;
function GMap0(B: in MapData; P: Point) return Point is
begin
return GMap0(B.G,P);
end GMap0;
function InvGMap0(G: in SFun; P: Point) return Point is
--- image (y,w) of (x,z) under map G: -Gx(x,y)=z, w=Gy(x,y)
Y: constant Scalar := P(1);
W: constant Scalar := P(2);
R: RadPair := (0.1,0.1);
X,T,Z: Scalar := ScalZero;
Gx,Gy,Gxx,Gxy,Gyy: SFun; -- make global if need efficiency
begin
if not Trunc then
Der(R,G,Gx,Gy);
NumLVal(Gy,Y,W,X);
R := Reps.Between(ValRadii(X,Y,G.A),Rho(G));
end if;
Der(R,G,Gx,Gy,Gxx,Gxy,Gyy);
LVal(Gy,Y,W,X);
T := Val(X,Y,Gxy);
if Inf(T) <= Zero then
Error("RG.InvGMap0 error: cannot verify twist");
end if;
Z := -Val(X,Y,Gx);
return (X,Z);
end InvGMap0;
function InvGMap0(B: in MapData; P: Point) return Point is
--- inverse of GMap0
begin
return InvGMap0(B.G,P);
end InvGMap0;
function FMap0(B: in MapData; P: Point) return Point is
--- image (y,w) of (x,z) under map F: -Fx(x,y)=z, w=Fy(x,y)
begin
return InvLambda(B,GMap0(B,Lambda(B,P)));
end FMap0;
function GMap1(B: in MapData; P: Point) return Point is
begin
return InvLambda(B,FMap0(B,GMap0(B,Lambda(B,P))));
end GMap1;
function GMap(B: in MapData; P: Point) return Point is
--- combine G0 and G1, not used
begin
if WellInside(P,B.R) then
return GMap0(B,P);
else
return GMap1(B,P);
end if;
end GMap;
function FMap(B: in MapData; P: Point) return Point is
--- not used
begin
return InvLambda(B,GMap(B,Lambda(B,P)));
end FMap;
function GLambda0(B: in MapData; P: Point) return Point is
begin
return GMap0(B,Lambda(B,P));
end GLambda0;
function GLambda1(B: in MapData; P: Point) return Point is
begin
return GMap1(B,Lambda(B,P));
end GLambda1;
procedure CheckMap(NX,NZ: in Positive; B: in MapData; D,R: in Rectangle) is
DX: constant Rep := (D(2)-D(1))/Rep(NX);
DZ: constant Rep := (D(4)-D(3))/Rep(NZ);
E: constant Scalar := DiskAt0(Small);
XL,ZL: Rep;
P,Q: Point;
begin
XL := D(1);
loop
Show("XL =",XL,3);
P(1) := Disk(XL,XL+DX)+E; -- ends up covering [D(1),D(2)]
ZL := D(3);
loop
P(2) := Disk(ZL,ZL+DZ)+E; -- ends up covering [D(3),D(4)]
Q := Map(B,P); -- also serves as domain check
if not WellInside(Q,R) then
Show1("Q =",Q);
Error("RG.CheckMap: range error");
end if;
ZL := ZL+DZ;
exit when ZL>D(4);
end loop;
XL := XL+DX;
exit when XL>D(2);
end loop;
end CheckMap;
procedure CheckG0 is new CheckMap(Map => GMap0);
procedure CheckGLambda0 is new CheckMap(Map => GLambda0);
procedure CheckGLambda1 is new CheckMap(Map => GLambda1);
procedure CheckMaps(B: in MapData) is
--- notice that Scale and CheckMap allow margins for closure
NX: constant Integer := 500;
NZ: constant Integer := 500;
RBig: constant Rectangle := (-Big,Big,-Big,Big);
D,Df,Rt: Rectangle;
begin
Show("Checking domain of G0");
CheckG0(NX,NZ,B,B.R,RBig); --- R is in the domain of G0
Show("Checking domain and range of G1 Lambda");
Scale(Inv(B.La),Inv(B.Mu),B.R,Df,Rt); --- Df is in the domain of F0
CheckGLambda1(NX,NZ,B,B.R,Df); --- G1*Lambda maps R inside Df
Show("Checking domain and range of G0 Lambda");
Scale(B.La,B.Mu,B.R,Rt,D); --- D contains Lambda*R
CheckGLambda0(NX,NZ,B,D,Df); --- G0*Lambda maps D inside Df
Show("Checks passed");
end CheckMaps;
begin
ShowType(ScalZero);
Show1(" PDeg =",PDeg);
Show1(" SVDim =",SVDim);
MakePackIndex;
end RG;