File : torus_ops.adb

with Messages, Reps.Ops, Lin.Ops, Numerics.Ops.Fou;
use Messages;

pragma Elaborate_All (Messages, Lin.Ops);

package body Torus_Ops is

  function IMin(I1,I2: Integer) return Integer renames Integer'Min;
  function IMax(I1,I2: Integer) return Integer renames Integer'Max;

  pragma Inline (Imin,Imax);

  procedure ShowParam is
  begin
    if Trunc then
      Show("       Scalar =    Numeric");
    else
      Show("       Scalar =   Interval");
    end if;
    Show("         Lmax =",Lmax);
    Show("         Wmax =",Wmax);
    Show("        Size1 =",Size1);
    Show("        Size  =",Size2);
    Show("        SizeE =",SizeE);
    Show("         Pdeg =",Pdeg);
    Show("          Dho =",Dho);
    Show("    DmaxCount =",DmaxCount);
    Show("     ErrCount =",ErrCount);
    Show(" Fourier'Size =",Fourier'Size);
    Show(" Sobolev'Size =",Sobolev'Size);
    Show("        Sigma =",Sigma);
    Show("        Kappa =",Kappa);
  end ShowParam;

  procedure Show(U: in Canonical; Cut: in Rep := Zero) is
  begin
    Show(U.Q,Cut);
    Show(U.P,Cut);
  end Show;

  procedure Read(Name: in String; R: in Weights; U: out Canonical; Decimal: in Boolean := True) is
  begin
    Read(Name & "_x",1,R,U.Q,Decimal);
    Read(Name & "_z",0,R,U.P,Decimal);
  end Read;

  procedure Write(Name: in String; U: in Canonical; Decimal: in Boolean := True) is
  begin
    Write(Name & "_x",U.Q,Decimal);
    Write(Name & "_z",U.P,Decimal);
  end Write;

  procedure Show(T: in Torus; Cut: in Rep := Zero) is
  begin
    Show(T.Q,Cut);
    Show(T.P,Cut);
  end Show;

  procedure Read(Name: in String; R: in Radius; T: out Torus; Decimal: in Boolean := True) is
  begin
    Read(Name & "_x",1,R,T.Q,Decimal);
    Read(Name & "_z",0,R,T.P,Decimal);
  end Read;

  procedure Write(Name: in String; T: in Torus; Decimal: in Boolean := True) is
  begin
    Write(Name & "_x",T.Q,Decimal);
    Write(Name & "_z",T.P,Decimal);
  end Write;

  function Rho(T: Torus) return Radius is
  begin
    return Radius'Min(Rho(T.Q),Rho(T.P));
  end Rho;

  function Rho(U: Canonical) return Weights is
  begin
    return Min(Rho(U.Q),Rho(U.P));
  end Rho;

  procedure SetZero(R: in Radius; T: out Torus) is
  begin
    SetZero(1,R,T.Q);
    SetZero(0,R,T.P);
  end SetZero;

  function "-"(T1,T2: Torus) return Torus is
    T: Torus;
  begin
    T.Q := T1.Q-T2.Q;
    T.P := T1.P-T2.P;
    return T;
  end "-";

  function Norm(R,B: Radius; T: Torus) return Scalar is
  begin
    return Max(Norm(R,T.Q),B*Norm(R,T.P));
  end Norm;

  function NormDer(R: Weights; B: Radius; U: Canonical) return Scalar is
    Rqq,Rqp,Rpq,Rpp: Scalar;
    Ht: Fourier;
  begin
    OmegaDq(R.Q,U.Q,Ht);
    Rqq := Norm(R,Ht);   Show("Rqq =", Sup(Rqq));
    Der3(R.P,U.Q,Ht);
    Rqp := Norm(R,Ht);   Show("Rqp =", Sup(Rqp));
    OmegaDq(R.Q,U.P,Ht);
    Rpq := Norm(R,Ht);   Show("Rpq =", Sup(Rpq));
    Der3(R.P,U.P,Ht);
    Rpp := Norm(R,Ht);   Show("Rpp =", Sup(Rpp));
    return Max(Rqq+Rqp/B,B*Rpq+Rpp);
  end NormDer;

  function Val(Z: Scalar; U: Canonical) return Scalar is
  begin
    return Z+Val(Z,U.P);
  end Val;

  function DerVal3(Z: Scalar; U: Canonical) return Scalar is
  begin
    return ScalOne+DerVal3(Z,U.P);
  end DerVal3;

  function ValQ0(U: Canonical) return Taylor is
    P: Taylor;
  begin
    P := ValQ0(U.P);
    AddComponent(1,ScalOne,P);
    return P;
  end ValQ0;

  procedure TmuComp(Mu: in Scalar; U1: in Canonical; U2: out Canonical) is
  begin
    Prod(-InvTheta,U1.Q,U2.Q);
    Prod(-Theta*Mu,U1.P,U2.P);
  end TmuComp;

  procedure Comp(U1,U2: in Canonical; U3: out Canonical) is
    Hp: Fourier;
  begin
    TraceEnter("Comp");
    Copy(U2.P,Hp);
    AddComponent(0,0,1,ScalOne,Hp);
    CompOmega3(U1.Q,U1.P,U2.Q,Hp,U3.Q,U3.P);
    Add(U2.Q,U3.Q);
    Add(U2.P,U3.P);
    TraceLeave;
  end Comp;

  procedure Comp(U: in Canonical; T1: in Torus; T2: out Torus) is
  begin
    TraceEnter("Comp");
    CompOmega3(U.Q,U.P,T1.Q,T1.P,T2.Q,T2.P);
    Add(T1.Q,T2.Q);
    Add(T1.P,T2.P);
    TraceLeave;
  end Comp;

  procedure TmuComp(Mu: in Scalar; T1: in Torus; T2: out Torus) is
  begin
    Prod(-InvTheta,T1.Q,T2.Q);
    Prod(-Theta*Mu,T1.P,T2.P);
  end TmuComp;

  procedure TmuComp(Mu: in Scalar; T: in out Torus) is
  begin
    Mult(-InvTheta,T.Q);
    Mult(-Theta*Mu,T.P);
  end TmuComp;

  procedure CompTinv(T1: in Torus; T2: out Torus) is
  begin
    CompTinv(T1.Q,T2.Q);
    CompTinv(T1.P,T2.P);
  end CompTinv;

  procedure NN(Mu: in Scalar; U: in Canonical; T1: in Torus; T2: out Torus) is
    T: Torus;
  begin
    TraceEnter("NN");
    Comp(U,T1,T2);
    CompTinv(T2,T);
    TmuComp(Mu,T,T2);
    TraceLeave;
  end NN;


end Torus_Ops;