function largeAmplitudeB_1D

set(0,'DefaultTextInterpreter','latex');
set(0,'DefaultLineLineWidth',1.0);
set(0,'DefaultAxesFontSize',15);
set(0,'DefaultTextFontSize',15);
set(0,'DefaultAxesLineWidth',1.0);
myblue   = [16 116 176]/256;
myred    = [227 35 31]/256;
myorange = [228 157 37]/256;

%
% Finite difference 1D code to evolve fields to become large amplitude
% Author: Jonathan Squire (modified for GPP2 HW03 by MWK)
%

amp    = 0.01;          % initial amplitude of fluctuation
theta0 = 85;            % angle (in degrees) between k and <B>
ifield = 2;             % choice of initial B-field fluctuation (see options below)
noplot = 0;             % set = 1 if you'd rather skip the plot (for speed)
pauseg = 0;             % set = 1 if you'd like to click to frame-advance the movie

tmax   = 10;            % final integration time
N      = 1024;          % number of cells
fdorder= 3;             % 3 = 6th-order finite differencing; 2 = 4th-order
L = 1;                  % length of box

tout = [0:0.1:(tmax-1)];
nt = length(tout);
ampmax = 100;
n2s= @(s) num2str(s);

p0 = [1 sind(theta0)/cosd(theta0) 0];    % wavenumber direction, assuming va=[1,0,0];

lam = (0:L/N:(L*(1-L/N))).';
k = 2i*pi/L * [0:N/2-1 -N/2:-1].'; % For ICs

%
% THREE DIFFERENT INITIAL CONDITIONS
%
% ifield=1: k = 1 monochromatic
% ifield=2: k = -3, -2, -1, 1, 2, 3 with random amplitudes and phases
% ifield=3: k = 3, 4, 5 with pre-set amplitudes and phases
%
%
if ifield==1
    bninit = cos(2*pi*lam);
end
if ifield==2
    rng("shuffle");
    filt = abs(k)<=6*pi/L; filt(1)=0;
    bninit = 1*real(ifft(filt.*fft(randn(size(lam)))));
end
if ifield==3
    bninit = - 0.15*sin(6*pi*lam+2) - 0.1*cos(8*pi*lam-1) +0.1*sin(10*pi*lam);
end


% normalize appropriately
bninit = bninit.*amp/sqrt(mean(bninit.^2));
binit = formICs(p0,bninit);

[dmat,ddmat,ddmatP]=finite_diff_mat(lam,fdorder);

options = odeset('RelTol',1e-6,'AbsTol',1e-10);
tic
[t,yout] = ode45(@(t,y) evalRHS(t,y,N,p0,dmat,ddmatP),tout, binit(:), options);
toc

avec = 1+t;
bx = yout(:,1:N).';
by = yout(:,N+1:2*N).';
bz = yout(:,2*N+1:3*N).';
bmag = sqrt((bx+1.').^2 + (by).^2 + bz.^2);
pdb = p0(1)*bx + p0(2)*by + p0(3)*bz;
bdb = bx./bmag;

if ~noplot 
    figure(1); clf;
    for ttt=1:nt
        g= tiledlayout(2, 1, 'TileSpacing', 'compact'); 
        nexttile;
        plot(lam,bmag(:,ttt),'k','Linewidth',1.2)
        hold on
        plot(lam,bx(:,ttt)+1,'Linewidth',1.2,'Color',myorange);
        plot(lam,by(:,ttt),'Linewidth',1.2,'Color',myblue);
        plot(lam,bz(:,ttt),'Linewidth',1.2,'Color',myred);
        plot(lam,0*lam,'k:','Linewidth',1.2);
        hold off
        ylim(ampmax*[-1 1]); %ylim(1.1*mean(bmag(:,ttt))*[-1 1])
        ax = gca; ax.TickLabelInterpreter = 'latex';
        title(['$t = ' n2s(avec(ttt)) '$'])
        xticklabels({})
        ylabel('$\textbf{\textit{B}}\,/\,|\langle\textbf{\textit{B}}\rangle|$')
        leg = legend('$~|B|$','$~B_x$','$~B_y$','$~B_z$', ...
                     'Location','SouthEast','Interpreter','latex','FontSize',16,'Orientation','horizontal');
        leg.ItemTokenSize = [32,20];
        legend('boxoff');
        nexttile;
        plot(lam,bdb(:,ttt),'k','Linewidth',1.2)
        ylim([-1.2 1.2]);
        ax = gca; ax.TickLabelInterpreter = 'latex';
        xlabel('$\lambda$'); ylabel('$\delta B_\parallel\,/\,|\textbf{\textit{B}}|$');
    
        drawnow
        %pause(0.01);
        if pauseg;ginput(1);end
    end
end

fprintf('maximum value of |dBprl/B| = %f\n',max(max(abs(bdb))));

end


function dydt = evalRHS(t,y,N,p,dmat,ddmat)
% Timescale defined by tA(t=0)=1, absorb alpha into d/dt
% Setup
p2 = vecnorm(p)^2;
pdva = p(1);
vA = [1 0 0 ];
bvec = reshape(y,[N,3]);

todia = @(v) spdiags(v,0,N,N);

b2tot = vecnorm(vA + bvec,2,2).^2;
phiOp = todia(b2tot*p2 - pdva^2)*ddmat;

phi = - phiOp\bvec(:,1);

dbdtau =  bvec - ...
    dmat*((dmat*phi).*( p2*(vA+bvec) - p*pdva));

dydt = dbdtau(:);

end


function S = solutionStats(S,t,l,bvec,dlmat)
% Various statistics from the evolution. 
b2 = (1+bvec(:,1)).^2 + bvec(:,2).^2 + bvec(:,3).^2;

% Statistics – steepness, CB2 
S.CB2 = [S.CB2 std(b2)./(std(bvec(:,1)).^2+std(bvec(:,2)).^2+std(bvec(:,3)).^2)];
S.Ql = [S.Ql mean((dlmat*bvec(:,1)).^2 + (dlmat*bvec(:,2)).^2 + (dlmat*bvec(:,3)).^2)...
            / mean(bvec(:,1).^2 + bvec(:,2).^2 + bvec(:,3).^2)];
% Measures of the gradients in different norms – do we form a discontinuity? 
for iii=1:3
    Q1(iii) = mean(abs(dlmat*bvec(:,iii)) )/mean(abs(bvec(:,iii)));
    Q2(iii) = (mean(abs(dlmat*bvec(:,iii)).^2 )/mean(abs(bvec(:,iii)).^2))^0.5 ;
    Q4(iii) =  (mean(abs(dlmat*bvec(:,iii)).^4 )/mean(abs(bvec(:,iii)).^4))^0.25 ;
    Q12(iii) =  (mean(abs(dlmat*bvec(:,iii)).^12 )/mean(abs(bvec(:,iii)).^12))^(1/12) ;
    Qinf(iii) =  (max(abs(dlmat*bvec(:,iii))) )/max(abs(bvec(:,iii)));
end
S.Q1 = [S.Q1 Q1.'];
S.Q2 = [S.Q2 Q2.'];
S.Q4 = [S.Q4 Q4.'];
S.Q12 = [S.Q12 Q12.'];
S.Qinf = [S.Qinf Qinf.'];        
        
S.A = [S.A sqrt(mean(bvec(:,1).^2 + bvec(:,2).^2 + bvec(:,3).^2))];
S.dBBx = [S.dBBx std(b2(:))/std(bvec(:,1).^2)];
mbx = mean(bvec(:,1));
S.dBprl = [S.dBprl sqrt(mean((bvec(:,1) - mbx).^2) )];

S.t = [S.t t];

end


function binit = formICs(p0,bninit)
vA0 = [1 0 0 ];
% Form ICs by constructing spherically polarized wave given bz
nh = cross(p0,vA0);
nh = nh/vecnorm(nh);
mh = cross(cross(p0,vA0),p0); 
mh = mh/vecnorm(mh);
mhvA = dot(vA0,mh);

C = fzero(@(C) -mhvA + mean(sqrt(C^2 - bninit.^2)), [1.00000001 1e6]*max(abs(bninit)) );

bminit =  1*(-mhvA + sqrt(C^2 - bninit.^2));
% bminit = -0.5*(bninit.^2 - mean(bninit.^2));

binit = [bminit*mh(1)+bninit*nh(1) bminit*mh(2)+bninit*nh(2) bminit*mh(3)+bninit*nh(3)];

end


function [dmat,ddmat,ddmatP]=finite_diff_mat(x,order)
% 1D finite difference matrices. dmat and ddmat standard periodic matrixes
% ddmatP has a point fixed in order to solve Poisson's equation
% order is actually twice the input number (e.g. order=2 is a 4th order solver)

% FD coefficients, from dfd https://au.mathworks.com/matlabcentral/fileexchange/76447-finite-difference-derivatives
% order=1 actually 2, 2 actually 4 etc.
fd_coeff = [      1/2     0       0       0;
                   2/3     -1/12   0       0; 
                  3/4     -3/20   1/60    0;
                   4/5     -1/5    4/105   -1/280];
% Second derivative    
fd_coeff_dd = [-2      1       0       0       0;
            -5/2    4/3     -1/12   0       0;
            -49/18  3/2     -3/20   1/90    0;
            -205/72 8/5     -1/5    8/315   -1/560];

dx = x(2)-x(1);
nx = length(x);

% First order
on = ones(nx,1)/dx;
matdiags = [];
for jjj = 1:order
    matdiags = [-fd_coeff(order,jjj)*on matdiags fd_coeff(order,jjj)*on];
end
% Go backwards to make periodic
for jjj = order:-1:1
    matdiags = [fd_coeff(order,jjj)*on matdiags -fd_coeff(order,jjj)*on];
end
dlevs = [1:order nx-order:nx-1];
dmat = spdiags(matdiags,[-fliplr(dlevs) dlevs],nx,nx);

% Second order
on = ones(nx,1)/dx.^2;
matdiags = fd_coeff_dd(order,1)*on;
for jjj = 2:order+1
    matdiags = [fd_coeff_dd(order,jjj)*on matdiags fd_coeff_dd(order,jjj)*on];
end
% Go backwards to make periodic
for jjj = order+1:-1:2
    matdiags = [fd_coeff_dd(order,jjj)*on matdiags fd_coeff_dd(order,jjj)*on];
end
dlevs = [1:order nx-order:nx-1];
ddmat = spdiags(matdiags,[-fliplr(dlevs) 0 dlevs],nx,nx);


% Make well posed for inversion by fixing value at one point
% https://scicomp.stackexchange.com/questions/27370/linear-solvers-how-to-deal-with-a-singular-system-poisson-equation-with-neuma
% Make sure sum(f)=0 in solution A.phi = f, then subtract mean(phi) after
ddmatP = ddmat;
idx = floor(nx/2);
% ddmatP(idx,:) = 0;
% ddmatP(idx,idx) = 1;
% idx = nx-1;
ddmatP(idx,:) = ddmatP(idx,:) + 1;
% Seems like a cluge, but works well for solving Poissons eqn
% Sometimes works fine with normal ddmat too – should try both

end
