File : ucomp.adb

with Reps, Ada.Text_IO, Messages, Fouriers_Init, Lin, Lin.IO;
use Reps, Ada.Text_IO, Messages, Fouriers_Init, Lin, Lin.IO;
with RG_Ops;

with Intervals, Intervals.Ops, Intervals.IO;
use Intervals, Intervals.Ops, Intervals.IO;
pragma Elaborate_All(Intervals.Ops, Intervals.IO);

procedure UComp 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 UComp_Ops is new RG_Ops (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho,
                                   Sigma => Sigma, Kappa => Kappa, Scalar => Scalar);
  use UComp_Ops;
  use Scal_Fou, Scal_Fou_Ops, Scal_Fou_IO;

  R:   constant Weights := (0.85,0.150);   -- for domain of RG
  Rpp: constant Weights := (0.90,0.165);   -- for domain of NN

  RBall: constant Radius := 3.0E-12;       -- for contraction MM

  Iname: constant String := "40-16-8";     -- name stem for input Hamiltonians
  Oname: constant String := "40-16-8";     -- name stem for output
  Yname: constant String := "85B04";       -- name stem for input matrices
  DecIO: constant Boolean := False;        -- use decimal IO

  Nyy: constant Positive := 294911;        -- maximum number of entries in sparse matrix W
  Err:   constant Radius := 1.0E-14;       -- precision for NumFixComp3 etc.

  F: File_Type;
  Rs: Weights;
  Hb1,Hb2: HBall;
  H1,Psi0,Hq,Hzp: Hamilton;
  W: Sparse(0 .. Nyy);
  C: Scalar;
  U: UBall;

begin
  Trace;                                   -- activate TraceEnter/TraceLeave
  ShowParam;
  Show("  R =",R);
  Show("Rpp =",Rpp);
  Rs := (Up(R.Q),Up(R.P));                 -- for range of RG
  Show(" Rs =",Rs);
  Show("Err =",Err);

  Read("Psi-" & Iname & ".nh",0,(Two,One),Psi0,DecIO);
  UpsiAux(Err,Rpp,Psi0,Hq,Hzp);            -- Rpp describes domain after HoToU
  Comp3(Hq,Hzp,H1);
  C := Component(0,0,1,Hzp);
  AddComponent(0,0,1,-ScalOne,Hzp);        -- remove contribution from Id
  Write("uzero_x",H1,DecIO);
  Write("uzero_z",Hzp,DecIO);
  SetComponent(0,0,1,C,Hzp);               -- restore contribution from Id

  Read("Fix-" & Iname & ".nh",0,(Two,One),H1,DecIO);
  Read("W-" & Yname & ".tdr",W);

  Hb1.R := R;
  Hb1.E := Scal(RBall);
  Show("Ball:",Sup(Hb1.E));

  Show("Computing (R-R1)(H1+M*Ball)");
  Show("Computing also Scaling and 2nd canonical transformation");
  RminusR1oM(Rs,W,Psi0,Hq,Hzp,H1,Hb1,Hb2,C,U);
  Show("Norm:",Sup(Hb2.E));
  Write("urest",U,DecIO);

  Show("Scaling interval, center and radius:");
  Put(Current_Output,C,DecIO);
  New_Line;
  Show("Writing scaling interval to c_h");
  Create(F,Out_File,"c_h");
  Put(F,C,DecIO);
  New_Line(F);
  Close(F);

end UComp;