File : checktorus.adb

with Reps, Ada.Text_IO, Messages;
use Reps, Ada.Text_IO, Messages;

with Intervals, Intervals.Ops, Intervals.IO;
use Intervals, Intervals.Ops, Intervals.IO;
pragma Elaborate_All(Intervals, Intervals.Ops, Intervals.IO);
-- for purely numeric computations replace "Interval" by "Numeric" in this program


with Fouriers_Init;
use Fouriers_Init;
with Torus_Ops;

procedure CheckTorus is

  subtype Scalar is Interval;

  Lmax: constant Positive := 40;
  Wmax: constant Positive := 40;
  Pdeg: constant Positive := 16;
  Dho:  constant Natural  :=  8;

  Sigma: constant Rep := 0.85001;
  Kappa: constant Rep := Sigma/0.4;

  package CheckTorus_Ops is new Torus_Ops (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho,
                                          Sigma => Sigma, Kappa => Kappa, Scalar => Scalar);
  use CheckTorus_Ops;
  use Scal_Fou, Scal_Fou_Ops, Scal_Fou_IO;

  HDec: constant Boolean := False;         -- decimal IO

  Rt:  constant Radius := 0.0001;          -- r, for norm of tori
  B:   constant Radius := 2.0;             -- b, for norm of tori
  Rad: constant Radius := 0.002;           -- radius of ball
  RQ:  constant Radius := 0.7;             -- domain for q-derivatives

  function Get_Mu return Scalar is
    InvTheta: constant Scalar := Lambda(0,-1);
    Ch: Scalar;
    F: File_Type;
  begin
    Open(F,In_File,"c_h");  -- saved by UComp
    Get(F,Ch,False);
    Close(F);
    return InvTheta*InvTheta*InvTheta/Ch;
  end Get_Mu;

  Mu,Ex,Ez,K: Scalar;
  R: Weights;
  T1,T2: Torus;
  U0,UR,U,Ut: Canonical;

begin
  Trace;
  Mu := Get_Mu;

  Read("uzero",(Two,One),U0,HDec);
  Read("urest",(Two,One),UR,HDec);
  Show("Computing U = U0 o UR");
  Comp(U0,UR,U);
  Write("u",U,HDec);

  Read("tapprox-40-16-8",Rt,T1,HDec);

  Ex := Norm(Rt,T1.Q)+Scal(Rad);
  Ez := Norm(Rt,T1.P)+Scal(Rad)/B;
  Show("E         =",Sup(Ex),Sup(Ez));

  R := (RQ,Sup(Ez));
  Show("rho       =",R.Q,R.P);
  TmuComp(Mu,U,Ut);
  K := NormDer(R,B,Ut);
  Show("|T_mu Du| =",Sup(K));

  K := (InvTheta+K*BetaMax(RQ,Rt,Ex))*Exp(Rt*Log(Theta));
  Show("K         =",Sup(K));
  Show("(1-K)*R   =",Inf(Rad*(ScalOne-K)));

  Show("Computing NN(T)-T");
  NN(Mu,U,T1,T2);
  Show("|NN(T)-T| =",SupAbs(Norm(Rt,B,T2-T1)));

end CheckTorus;