%                   coefsRKNmod.m
% Coefficients of selected RKN methods with modified potentials
%  S_h = (a_{s+1} (b_s,c_s) a_s ... (b_1,c_1) a_1)
% Methods with number of stages s={1 2,3,4,5,11}
% 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] = coefsRKNmod(n,ab)
% Number of stages 
nk=n;
if (n==1)  %Strang/leap-frog (effective order 4)
    if (ab==1)
        a = [1/2  , 1/2 ];
        b = [1 , 0];
        c = [1/24 , 0];       
    elseif (ab==2)
          nk=1;
        a = [0  , 1 ];
        b = [1/2, 1/2];
        c = [1/48 , 1/48];
    end
elseif (n==2)  % 2-stage 4th-order Koseleff--Chin (KOS93,CHI97)
    if (ab==1) %ABA
         a1=1/2 - 1/(2*sqrt(3));
        a = [a1  , 1-2*a1 , a1  ];
        b = [1/2, 1/2 , 0];
        c = [1/24-sqrt(3)/48 , 1/24-sqrt(3)/48 , 0];
        
    elseif (ab==2) %BAB
        a = [0  , 1/2 , 1/2  ];
        b = [1/6, 2/3 , 1/6];
        c = [0 , 1/72 , 0];

    end
elseif (n==3)  %3-satage 4th-order methods (OMF03)
    if (ab==1) %ABA
          a1 = 0.1159953608486416; 
          b1 =	0.2825633404177051;
          c1 = 0.001226088989536361;
          a2 =1/2 - a1;
          b2 = 1 - 2*b1;
          c2 = 0.003035236056708454;

         a = [a1 , a2 , a2 , a1 ];
         b = [b1 , b2 , b1 , 0];
         c = [c1 , c2 , c1 , 0];
    elseif (ab==2) %BAB

          b1 = 0.08002565306418866; 
          c1 = 0.0002725753410753895;
          a1 = 0.2728983001988755;
          b2 = 1/2 - b1;
          c2 = 0.002960781208329478;
          a2 = 1 - 2*a1;

         a = [0  , a1 , a2 , a1 ];
         b = [b1 , b2 , b2 , b1];
         c = [c1 , c2 , c2 , c1];
    elseif (ab==3) %ABA (witho only 1 modified potential)

          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];
    end
elseif (n==4)  % 4-satage 4th-order methods (OMF03)
    if (ab==1) % ABA
          a1 = 0.09325912861071900;  
          b1 = 0.2247800288685984;
          c1 = 0.0005562281089130940;
          a2 = 0.2791634819768266;
          b2 = 1/2 - b1;
          c2 = 0.0008405927247441154;
          a3 = 1 - 2*(a1 + a2);

        a = [a1 , a2 , a3 , a2 , a1 ];
        b = [b1 , b2 , b2 , b1 , 0 ];
        c = [c1 , c2 , c2 , c1 , 0 ];
    elseif (ab==2) % BAB
        % The 4-stage 4th-order CABACABAC method of Omelyan et al. eq. (53)
        b1 = 0.05851872613455621; 
        c1 =0.4339598806816256e-03;
        a2 = 0.1921125277429464;
        b2 = 0.2852162240687091;
        c2=0;
        a3 = (1 - 2*a2)/2;
        b3 = 1 - 2*(b1 + b2); 
        c3 =0.2427475259663050e-02;

        a = [0  , a2 , a3 , a3 , a2 ];
        b = [b1 , b2 , b3 , b2 , b1];
        c = [c1 , c2 , c3 , c2 , c1];
    end
elseif (n==5)  % 5-satage 4th-order methods (OMF03)
    if (ab==1)
% order 4 ABA
    a1=0.6419108866816235e-1;
    a2=0.1919807940455741;
    b1=0.1518179640276466;
    b2=0.2158369476787619;
    c1=0;
    c2=0.9628905212024874e-3;
    c3=0;
          b3 = 1 - 2*(b1+b2);
          a3 =1/2 - (a1 + a2);

        a = [a1 , a2 , a3 , a3 , a2 , a1 ];
        b = [b1 , b2 , b3 , b2 , b1 , 0 ];
        c = [c1 , c2 , c3 , c2 , c1 , 0 ];
    elseif (ab==2) % BAB
    a1=0.2029270564692829;
    a2=0.1926052063353027;
    b1=0.6668876199434440e-1;
    b2=0.2620356629687677;
    c1=0;
    c2=0.1042387551227681e-02;
    c3=0;
          b3 = 1/2 - (b1+b2);
          a3 = 1 - 2*(a1 + a2);

        a = [0  , a1 , a2 , a3 , a2 , a1 ];
        b = [b1 , b2 , b3 , b3 , b2 , b1];
        c = [c1 , c2 , c3 , c3 , c2 , c1];
    elseif (ab==3) % order 6 ABA
    a1=0.1097059723948682;
    a2=0.4140632267310831;
    b1=0.2693315848935301;
    b2=1.1319803486515564;
    c1=0.8642161339706166e-3;
    c2=0;
    c3=-0.1324638643416052e-1;
          b3 = 1 - 2*(b1+b2);
          a3 =1/2 - (a1 + a2);

        a = [a1 , a2 , a3 , a3 , a2 , a1 ];
        b = [b1 , b2 , b3 , b2 , b1 , 0 ];
        c = [c1 , c2 , c3 , c2 , c1 , 0 ];
    elseif (ab==4) % order 6 BAB
    a1=(1+1/sqrt(5))/2;
    a2=-1/sqrt(5);
    b1=1/12;
    b2=5/12-sqrt(50+22*sqrt(5))/24;
    c1=0;
    c2=(15+5*sqrt(5))/1152-sqrt(50+22*sqrt(5))*(1/2880+sqrt(5)/1152);
    c3=-(11+5*sqrt(5))/1152+sqrt(50+22*sqrt(5))*(1/2880+sqrt(5)/1152);
          b3 = 1/2 - (b1+b2);
          a3 = 1 - 2*(a1 + a2);

        a = [0  , a1 , a2 , a3 , a2 , a1 ];
        b = [b1 , b2 , b3 , b3 , b2 , b1];
        c = [c1 , c2 , c3 , c3 , c2 , c1];
    end
elseif (n==11) %Omelyan 11-satage order 8 (OMF03)
    if (ab==1) % order 8 ABA
        a1 = 0.41009674738801111928784693005080;
        b1 =0.48249309817414952912695842664785e-2;
        c1 =0.14743936907797528364717244760736e-3;
        a2 =-0.34123345756052780489101697378499;
        b2 =0.17492394861090375603419001374207;
        c2 =0.23288450531932545357194967600155e-3;
        a3 =0.25644714021068150492361761631743;
        b3 =0.29304366370957066164364546204288;
        c3 =0.61648659635535962497705619884752e-2;
        a4 =0.27765273975812438394100476242641;
        b4 =0.47448940168459770284238136482511e-1;
        c4 =-0.12307516860831240716732016960034e-1;
        a5 =-0.56926266869753773902939657321159;
        b5 =-0.15299863411743974499219652320477e-2;
        c5 =-0.73296648559126385387017161643798e-4;
        %a6 = 0.46629949890124853576794423820194;
        %b6=-0.37422994259002571606842462603791e-1;
        c6 = 0.15295860994523744731993293847001e-1;

          b6 =1 - 2*(b1+b2+b3+b4+b5);
          a6 =1/2 - (a1+a2+a3+a4+a5);

        a = [a1, a2, a3, a4, a5, a6, a6, a5, a4, a3, a2, a1];
        b = [b1, b2, b3, b4, b5, b6, b5, b4, b3, b2, b1,  0];
        c = [c1, c2, c3, c4, c5, c6, c5, c4, c3, c2, c1,  0];
    elseif (ab==2) % order 8 BAB
            b1 = 0.1839699354244402;
            c1 = 0;
            a1 = 0.6922517172738832 ;
            b2 = 0.7084389757230299;
            c2 = 0.3976209968238716E-1;
            a2 = -0.3183450347119991 ;
            b3 = 0.1981440445033534;
            c3 = 0.2245403440322733E-1;
            a3 = 0.6766724088765565;
            b4 = -0.6409380745116974E-1;
            c4 = 0.9405266232181224E-3;
            a4 = -0.7207972470858706 ;
            b5 = -0.6887429532761409;
            c5 = -0.7336500519635302E-1;
            a5 = 0.3580316862350045; 
            b6 = 0.1622838050764871;
            c6 = 0.2225664796363730E-1;
            a6=-0.3756270611751488;

          b6 =1/2 - (b1+b2+b3+b4+b5);
          a6 = 1 - 2*(a1+a2+a3+a4+a5);

        a = [ 0, a1, a2, a3, a4, a5, a6, a5, a4, a3, a2, a1];
        b = [b1, b2, b3, b4, b5, b6, b6, b5, b4, b3, b2, b1];
        c = [c1, c2, c3, c4, c5, c6, c6, c5, c4, c3, c2, c1];
    end
end

