File : scaling.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 Scaling 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 Scaling_Ops is new Torus_Ops (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho,
                                          Sigma => Sigma, Kappa => Kappa, Scalar => Scalar);
  use Scaling_Ops;
  use Scal_Fou, Scal_Fou_Ops, Scal_Fou_IO;
  use Scal_Tay, Scal_Tay_Ops;

  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

  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;

  function Contains(Z1,Z2: Scalar) return Boolean is
  begin
    if IsNumeric(ScalZero) then return True; end if;
    return (Inf(Z1)<Inf(Z2)) and (Sup(Z1)>Sup(Z2));
  end Contains;

  Rp: Radius;
  RpHalf,Mu,ZR,J0,J1,J2,J3,DU: Scalar;
  U: Canonical;
  T: Torus;
  P: Taylor;

begin
  Trace;
  Mu := Get_Mu;
  Read("u",(Two,One),U,HDec);
  P := ValQ0(U);
  Rp := Rho(U).P;

  RpHalf := Half*Scal(Rp);
  Read("tapprox-40-16-8",Rt,T,HDec);
  if Sup(Norm(Rt,T.P)+Scal(Rad)/B) >= Inf(RpHalf) then
    Message("Scaling",Domain_Error);
  end if;

  ZR := RpHalf/Rep(200);
  J0 := Ball(RpHalf);
  J3 := Scal(-0.014675)+Ball(Rep(0.000001));

  for M in -100 .. 100 loop
    J2 := (Rep(M)*RpHalf)/Rep(100)+Ball(ZR);
    if Inf(DerVal(Rp,J2,P)) <= Zero then
      Message("Scaling",Domain_Error);   -- P not strictly monotone
    end if;
    for N in 1 .. 15 loop                -- check contraction to J3
      J1 := J2;
      J2 := -Mu*Theta*Val(Rp,J1,P);
      if not Contains(J0,J2) then
        Message("Scaling",Domain_Error); -- got mapped outside J0
      end if;
    end loop;
    if not Contains(J3,J2) then
      Message("Scaling",Domain_Error);
    end if;
  end loop;

  J2 := J3;
  for N in 1 .. 11 loop                 -- contract J3 down to J
    J1 := J2;
    J2 := -Mu*Theta*Val(Rp,J1,P);
    if not Contains(J1,J2) then
      Message("Scaling",Domain_Error);
    end if;
  end loop;

  DU := -Mu*Theta*DerVal(Rp,J1,P);

  Show("  -1/theta : ",Inf(-InvTheta),Sup(-InvTheta),17);
  Show("        mu : ",Inf(Mu),Sup(Mu),17);
  Show("         J : ",Inf(J1),Sup(J1),17);       -- contains z_infinity
  Show("lambda_tau : ",Inf(Theta),Sup(Theta),17);
  Show("  lambda_x : ",Inf(Mu/DU),Sup(Mu/DU),17);
  Show("  lambda_y : ",Inf(Mu/Theta),Sup(Mu/Theta),17);
  Show("  lambda_z : ",Inf(DU),Sup(DU),17);

end Scaling;