Regression fit for ENSO

gem-lte.ads

      Doodson_Set : Doodson_List :=
     (
       (1, 0, 0, 0.0, 1.0), -- T
       (1, 0, 0, 1.0, 1.0), -- D
--       (0, 0, 2, 2.0, 1.0), -- 3Y
       (0, 0,-2, 2.0, 1.0), -- 5.7 
--       (0, 0, 0, 2.0, 1.0), -- N/2
       (2, 0, 0, 1.0, 1.0), -- DT
       (2, 0, 0, 0.0, 1.0), -- t
       (2, 0, 0, 2.0, 1.0), -- d
--       (1, 0,-1, 0.0, 1.0), -- A
--       (2, 0,-2, 0.0, 1.0), -- a
       (0, 0, 0, 1.0, 1.0), -- N
       (0, 0, 2, 0.0, 1.0), -- p 
--       (1, 0, 1, 0.0, 1.0), -- E
       (2, 0,-1, 1.0, 1.0), -- AD
       (0, 0, 2, 1.0, 1.0), -- 3.56
--       (1, 0,-1, 1.0, 1.0), -- A+
--       (1, 0,-1,-1.0, 1.0), -- A-
       (0, 0, 1, 1.0, 1.0), -- 6Y
--       (1, 0, 1, 1.0, 1.0), -- 26.985
--       (0, 0, 1,-1.0, 1.0),-- 16.8y
       (0, 0,-1, 0.0, 1.0) -- P
       --(1, 0, 1, -1.0, 1.0) -- Fake
--       (3, 0,-1, 1.0, 1.0), -- Mt'
--       (3, 0,-1, 0.0, 1.0),  -- Mt
--       (3, 0,-1, 2.0, 1.0)  -- Mtd
     );

      Doodson_RSet : Doodson_List :=
     (
       (1, 0, 0, 0.0, 1.0), -- T
       (1, 0, 0, 1.0, 1.0), -- D
--       (0, 0, 2, 2.0, 1.0), -- 3Y
--       (0, 0,-2, 2.0, 1.0), -- 5.7 
--       (0, 0, 0, 2.0, 1.0), -- N/2
       (2, 0, 0, 1.0, 1.0), -- DT
       (2, 0, 0, 0.0, 1.0), -- t
       (2, 0, 0, 2.0, 1.0), -- d
       (1, 0,-1, 0.0, 1.0), -- A
       (1, 0, 1, 0.0, 1.0), -- evect
--       (0, 0, 0, 1.0, 1.0), -- N
--       (0, 0, 2, 0.0, 1.0), -- p 
--       (1, 0, 1, 0.0, 1.0), -- E
       (2, 0,-1, 1.0, 1.0) -- AD
--       (0, 0, 2, 1.0, 1.0), -- 3.56
--       (1, 0,-1, 1.0, 1.0), -- A+
--       (1, 0,-1,-1.0, 1.0), -- A-
--       (0, 0, 1, 1.0, 1.0), -- 6Y
--       (1, 0, 1, 1.0, 1.0), -- 26.985
--       (0, 0, 1,-1.0, 1.0),-- 16.8y
--       (0, 0,-1, 0.0, 1.0) -- P
     );
   
     Annual_Set : Doodson_List :=
     (
       (0, 1, 0, 0.0, 365.242), -- annual
       (0, 2, 0, 0.0, 365.242) -- semi
     );
   

GEM-Mix-Regression,adb

with Text_IO;
with GEM.LTE.Primitives;

function GEM.Mix_Regression (File_Name : in String) return Gem.LTE.Period_Set is
   use GEM.LTE, GEM.LTE.Primitives;
   D : Data_Pairs := Make_Data(File_Name);
   First, Last : Integer;
   Singular : Boolean;
   Forcing : Data_Pairs := D;
   Model : Data_Pairs := D;
   P_A1, P_A2, P_A3, P_A4, P_A5 : Gem.LTE.Long_Periods := Gem.LTE.LP_Set;
   P_B1                         : Gem.LTE.Long_Periods := Gem.LTE.LP_RSet;
   P_Annual                     : Gem.LTE.Long_Periods := Gem.LTE.LP_Annual;
   A_A1, A_A2, A_A3, A_A4, A_A5 : Gem.LTE.Long_Periods_Amp_Phase := Gem.LTE.LPAP_Set; 
   A_B1                       : Gem.LTE.Long_Periods_Amp_Phase := Gem.LTE.LPAP_RSet;
   A_Annual                     : Gem.LTE.Long_Periods_Amp_Phase := Gem.LTE.LPAP_Annual;
   Level, K0, Trend, Accel : Long_Float := 0.0;
   Last_Time : Long_Float;
   Freq, Freq2 : Long_Float;
   One : constant Long_Float := 1.0;
begin
   First := D'First;
   Last := D'Last;
   Text_IO.Put_Line("records from " & First'Img & " ... " & Last'Img);
   Text_IO.Put_Line("factors from " & P_A1'First'Img & " ... " & P_A1'Last'Img);
   for I in Forcing'Range loop
      Forcing(I).Value := Forcing(I).Date;
      Last_Time := Forcing(I).Value;
   end loop;
   Text_IO.Put("updated forcing  ");
   Put(Last_Time);
   Text_IO.New_Line;
   for I in P_A1'Range loop -- to aliased frequency
      Freq := GEM.LTE.Year_Length/P_A1(I);
      Freq := abs Long_Float'Remainder(Freq,One);
      P_A1(I) := Freq;
      P_A2(I) := One - Freq;
      P_A3(I) := One + Freq;
      P_A4(I) := 2.0*One - Freq;
      P_A5(I) := 2.0*One + Freq;
   end loop;

   for I in P_B1'Range loop -- to aliased frequency
      Freq := 2.0*GEM.LTE.Year_Length/P_B1(I);
      Freq2 := 4.0*GEM.LTE.Year_Length/P_B1(I);
      Text_IO.Put_Line("Q=" & Freq2'Image);
      --Text_IO.Put_Line("F=" & Freq'Image);
      if Integer(Freq) mod 2 = 1 then
         Freq := abs Long_Float'Remainder(Freq,One);
      else
         Freq := One - abs Long_Float'Remainder(Freq,One);
      end if;
      P_B1(I) := Freq*0.5;
      -- P_B2(I) := (One-Freq)*0.5;
      
--      if Integer(Freq2) mod 2 = 1 then
--         Freq2 := abs Long_Float'Remainder(Freq2,One);
--      else
--         Freq2 := One - abs Long_Float'Remainder(Freq2,One);
--      end if;
--      P_Q1(I) := Freq2*0.25;
   end loop;

   declare
      P : Gem.LTE.Long_Periods := P_A1 & P_A2 & P_A3 & P_A4 & P_A5 & P_B1 
        & P_Annual; -- & P_Q1;
      A : Gem.LTE.Long_Periods_Amp_Phase := A_A1 & A_A2 & A_A3 & A_A4 & A_A5 & A_B1 
        & A_Annual; -- & A_Q1;
   begin
   Trend := 0.0;
   Regression_Factors (Data_Records => D(800..1700), -- Time series
                       Forcing => Forcing(800..1700),  -- Value @ Time
                       NM => P'Last, -- # modulations
                       DBLT => P, --D.B.LT,
                       DALTAP => A, --D.A.LTAP,
                       DALEVEL => LEVEL,
                       DAK0 => K0,
                       Secular_Trend => Trend,
                       Accel => Accel,
                       Singular => Singular
                       );

   Text_IO.Put_Line("Singular? " & Singular'Img);
   for I in P'Range loop
      Put(P(I)); Text_IO.Put("   ");
      -- INTEGRATE dLOD -- DBLTAP(I).Amplitude := DBLTAP(I).Amplitude * 26.736 / DBLT(I);
      Put(A(I).Amplitude); Text_IO.Put("   ");
      -- INTEGRATE dLOD -- DBLTAP(I).Phase := DBLTAP(I).Phase + Ada.Numerics.Pi/2.0;
      Put(A(I).Phase); Text_IO.Put("   ");
      Text_IO.New_Line;
   end loop;

   Model := Gem.LTE.Primitives.LTE(Forcing => Forcing,
                Wave_Numbers => P,
                Amp_Phase => A,
                Offset => Level,
                K0 => K0,
                Trend => 0.0,
                NonLin => 1.0);
   Put(CC(D,Model), "=CC " );
      Put(Year_Length, "=Yr", True);
      
   GEM.LTE.Primitives.Save(Model, D, Forcing);

      declare
         PS : constant Gem.LTE.Period_Set := (N => P'Length,
                                  LP => P,
                                  AP => A);
      begin
         return PS;
      end;
      --return A;
   end;

end;

CW fingerprint. Note how the draconic CW frequency at 0.845 is replicated at 1.845 -- modulation of an annual impulse train.
Note the minor peaks at 0.5, 1.5, 2.5 -- that's a biennial modulation of the annual carrier
Also minor peaks at 1.678, 2.678 -- related to the 1/D+1/A beat of the 1/D-1/A 6 year frequency.

Re-evaluation of the experiment described here https://geoenergymath.com/2021/01/07/chandler-wobble-forcing/