File : run_rminusr1.adb

with Reps, Messages, Fouriers_Init, Lin, Lin.IO;
use Reps, 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 Run_RminusR1 is

  subtype Scalar is Interval;

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

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

  package Run_RminusR1_Ops is new RG_Ops (Lmax => Lmax, Wmax => Wmax, Pdeg => Pdeg, Dho => Dho,
                                          Sigma => Sigma, Kappa => Kappa, Scalar => Scalar);
  use Run_RminusR1_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
  HDec: constant Boolean := False;         -- use decimal IO for Hamiltonians

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

  Rs: Weights;
  Hb1,Hb2: HBall;
  H1,Psi0,Hq,Hzp: Hamilton;
  W: Sparse(0 .. Nyy);

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("W-" & Yname & ".tdr",W);

  Read("Fix-" & Iname & ".nh",0,(Two,One),H1,HDec);
  Read("Psi-" & Iname & ".nh",0,(Two,One),Psi0,HDec);

  UpsiAux(Err,Rpp,Psi0,Hq,Hzp);            -- Rpp describes domain after HoToU
  Write("Hzpp-" & Oname & ".nh",Hzp,HDec); -- save Hzp, just in case ...

  for Fac in 1 .. 2 loop
    Hb1.R := R;
    Hb1.E := Scal(Rep(Fac)*RBall);
    Show("Ball:",Sup(Hb1.E));

    Show("Computing (R-R1)(H1+M*Ball)");
    RminusR1oM(Rs,W,Psi0,Hq,Hzp,H1,Hb1,Hb2); -- NminusN1 eats domain after HoToU
    Show("Norm:",Sup(Hb2.E));
  end loop;

end Run_RminusR1;