%                   coefsRKNmod_proc.m
% Coefficients of selected RKN methods with modified potentials
% S_h= P_h^(-1) K_h P_h 
%  = (-ap_{1}, -(bp_1,cp_1), ..., -ap_k, -(bp_k,cp_k)
%    (a_{s+1}, (b_s,c_s), a_s, ..., (b_1,c_1), a_1)
%    ((bp_k,cp_k), ap_k, ..., (bp_1,cp_1), ap_1)
% Methods with number of stages s={1 2,3}
% There aresymmetric compositions ABA and BAB 
%  compositions: sa=1 calls for a symmetric ABA scheme and sa=2 for the BAB
% S. Blanes, F. Casas and A. Murua, Splitting methods for differential equations,
%   Acta Numerica (2024)     (Section 8, Table 8.4)


function [a b c ap bp cp] = coefsRKNmod_proc(s,sa)
% s = Number of stages 
if (s==1)  % 1-stage 4th-order: TI84-ROW91-WHT96-BCR99
        a = [1/2  , 1/2 ];
        b = [1 , 0];
        c = [1/24 , 0];   
        
        z1 = 0 ;
        z2 =-0.8749306155955435 ;
        z3 = 0.2371066801510219 ;
        z4 = 0.5363539829039128 ;
        
        y1 =-0.1859353996846055 ;
        y2 =-0.0731969797858114; 
        y3 = 0.1576624269298081; 
        y4 = 0 ;
    
    ap=-[z1 z2 z3 z4];
    bp=-[y1 y2 y3 y4];
    cp= [0 0 0 0 ];
        
    
elseif (s==2)  % 2-stage 5th-order non-symmetric kernel (BCR99)

    if sa==1
        a1=-2/5;
        b1=(42+sqrt(1389))/84;
        c1=(26391+868*sqrt(1389))/933408;
        c2=(26391-868*sqrt(1389))/933408;
        a = [1-a1 , a1];
        b = [1-b1 , b1];
        c = [c2 , c1];
    
        z1 =-0.2598277426708327 ;
        z2 =-0.8789034357785656;
        y1 =-0.6249287013618940;
        y2 = 0.0087247119259793 ;
        v1 =-0.0156990164295207;
        v2 =-0.0008843356433671;
        ap=[z1 z2];
        bp=[y1 y2];
        cp=[v1 v2];
    else
        
    end

elseif (s==3)   % 3-stage 6th-order symmetric ABA kernel (BCR99)

    if sa==1
          b1 = 0.2621129352517028; 
          c1 = 0;
          a1 = -0.0682610383918630;
          b2 = 1 - 2*b1;
          c2 = 0.0164011128160783;
          a2 = 1/2 - a1;

         a = [a1 , a2 , a2 , a1];
         b = [b1 , b2 , b1 ,  0];
         c = [c1 , c2 , c1 ,  0];
    
        z1 = 0.1604630501234888 ;
        z2 =-0.1222126706298830 ;
        z3 = 0.1916801124727711;
        z4 = 0.5630722377955035 ;
        z5 =-0.7612758792358986 ;
        z6 =-(z1+z2+z3+z4+z5);
        y1 =-0.012334538446142270 ;
        y2 =-0.6610294848488182 ;
        y3 =-0.023112349678219939;
        y4 = 0.000181521815949959 ;
        y5 = 2.3768244683666757 ;
        y6 =-(y1+y2+y3+y4+y5);
        v1 = 0.013816178183636998 ;
        v2 =-0.050288359617427786 ;
        v3 =-0.013462400168471472;
        v4 = 0.000603819193361427 ;
        v5 =-0.01 ;
        v6 = 0.01;
            
        ap=-[z1 z2 z3 z4 z5 z6];
        bp=-[y1 y2 y3 y4 y5 y6];
        cp=-[v1 v2 v3 v4 v5 v6];
    else
        
    end
end

