%---------------------------------------------------------------------------------------------
%
%---------------------------------------------------------------------------------------------
clear all; format short g; format compact
isave = 1;
if isave == 0
    %---------------------------------------------------------------------------------------------
    Zz = 92; Nn = 146; Aa = Nn + Zz;
    %---------------------------------------------------------------------------------------------
    % units fm (space, momentum^-1) and MeV (energy, mass, time^-1)
    % Coulomb now correct 
    %---------------------------------------------------------------------------------------------
    Lx = 64; Ly = 64; Lz = 64; 
    Nx = 64; Ny = 64; Nz = 64; dx = Lx/Nx; dy = Ly/Ny; dz = Lz/Nz; dV = dx*dy*dz;
    N3 = 3*max([Nx,Ny,Nz]); Ll = max([Lx,Ly,Lz]); L = sqrt(3)*Ll;
    x  = dx*(-Nx/2:Nx/2-1); y = dy*(-Ny/2:Ny/2-1); z = dz*(-Nz/2:Nz/2-1);
    kx = 2*pi/Lx*[0:Nx/2-1,-Nx/2:-1];
    ky = 2*pi/Ly*[0:Ny/2-1,-Ny/2:-1];
    kz = 2*pi/Lz*[0:Nz/2-1,-Nz/2:-1];
    [Kx,Ky,Kz] = ndgrid(kx,ky,kz); K2 = Kx.^2 + Ky.^2 + Kz.^2;
    Kx(Nx/2+1,:,:) = 0; Ky(:,Ny/2+1,:) = 0; Kz(:,:,Nz/2+1) = 0;  % required for calculation of first order derivatives
    %---------------------------------------------------------------------------------------------
    hbarc  = 197.326;
    m_n = 939.565; h2_2m_n  = hbarc^2/(2*m_n); aN_n = h2_2m_n*(3*pi^2)^(2/3);
    m_p = 938.272; h2_2m_p  = hbarc^2/(2*m_p); aN_p = h2_2m_p*(3*pi^2)^(2/3);
    %---------------------------------------------------------------------------------------------
    [X,Y,Z] = ndgrid(x,y,z); rr = X.^2 + Y.^2 + Z.^2; %clear X Y Z x y z kx ky kz
    r     = sqrt(rr); R = 1.12*Aa^(1/3)-0.86/Aa^(1/3); a = 0.54;
    phi_p = 1./(1+cosh(r/a)/cosh(R/a)); norm = 1/sqrt( sum(sum(sum( phi_p.^2 )))*dV );
    phi_p = phi_p*norm; phi_n = phi_p; Vc = phi_p*0;
    %---------------------------------------------------------------------------------------------
     
        % nuclear_binding_energy.txt, only energies, no charge formfacxtors
        at = -741.218455246; bt = 939.925985887; ct = 138.763424223; ct1 = -160.582591473; eta_t = 0.470292092182;
    
    gam = 7/3; e2  =  hbarc/137.036; hbarc_t = hbarc*sqrt(eta_t);
    h2_2m_n_t = h2_2m_n*eta_t;  h2_2m_p_t = h2_2m_p*eta_t;
    Lambda = hbarc_t*max([pi/dx,pi/dy,pi/dz]); m_N = 2*m_n*m_p/(m_n + m_p); Ec = Lambda^2/(2*m_N); tau = 0.25/Ec; % [MeV^-1]
    %---------------------------------------------------------------------------------------------
    qx = 2*pi/Ll*[0:N3/2-1,-N3/2:-1]/3;
    qy = 2*pi/Ll*[0:N3/2-1,-N3/2:-1]/3;
    qz = 2*pi/Ll*[0:N3/2-1,-N3/2:-1]/3;
    [Qx,Qy,Qz]= ndgrid(qx,qy,qz); Q2 = Qx.^2 + Qy.^2 + Qz.^2;
    Vq = 4*pi*e2*(1-cos(L*sqrt(Q2)))./(Q2+eps); Vq(1,1,1) = 4*pi*e2*L^2/2; clear qx qy qz Qx Qy Qz Q2
    %-------------------------------------------------------------------------------------------------
    % Now we find the ground state of the nucleus
    %-------------------------------------------------------------------------------------------------
    tfin   = 0.6; nwrite = 30; nsteps = ceil(tfin/nwrite/tau); itime = 0;  fftw('planner','measure');
    %-------------------------------------------------------------------------------------------------
    tic
    for it = 1:nwrite
        fn = fftn(phi_n); fp = fftn(phi_p);
        ek_n = exp( -tau*h2_2m_n*K2/2*eta_t ); phi_n = ifftn(ek_n.*fn); % half step
        ek_p = exp( -tau*h2_2m_p*K2/2*eta_t ); phi_p = ifftn(ek_p.*fp);
        for kt = 1:nsteps
            itime = itime + 1; tcur = itime*tau;
            %-------------- Potentials -----------------------------------------
            rho_p  = Zz*abs(phi_p).^2; rho_n  = Nn*abs(phi_n).^2;
            rho_0  = rho_n + rho_p; rho_1  = rho_n - rho_p;
            U_p    = 2*at*rho_0 + gam*bt*rho_0.^(gam-1) - ct/3*rho_1.^2./rho_0.^(4/3);
            U_n    = U_p + 2*ct*rho_1./rho_0.^(1/3) + 2*ct1*rho_1 + aN_n*rho_n.^(2/3);
            U_p    = U_p - 2*ct*rho_1./rho_0.^(1/3) - 2*ct1*rho_1 + aN_p*rho_p.^(2/3);
            if e2 > 0
                %-------------   Coulomb potential ---------------------------------
                rho3   = zeros(N3,N3,N3); rho3(1:Nx,1:Ny,1:Nz) = rho_p;
                rho3   = ifftn(Vq.*fftn(rho3)); Vc = rho3(1:Nx,1:Ny,1:Nz);
                U_p    = U_p + Vc - e2*(3/pi)^(1/3)*rho_p.^(1/3);
            end
            %-------------------------------------------------------------------
            phi_n = exp( -tau*U_n ).*phi_n; phi_p = exp( -tau*U_p ).*phi_p;
            %-------------------------------------------------------------------
            fn   = fftn(phi_n); fp   = fftn(phi_p);
            if     kt  < nsteps;  ek_n = exp( -eta_t*tau*h2_2m_n*K2   );  ek_p = exp( -eta_t*tau*h2_2m_p*K2   ); % full step
            elseif kt == nsteps ; ek_n = exp( -eta_t*tau*h2_2m_n*K2/2 );  ek_p = exp( -eta_t*tau*h2_2m_p*K2/2 ); % half step
            end
            phi_n = ifftn(ek_n.*fn); norm = 1/sqrt(sum(sum(sum(abs(phi_n).^2)))*dV); phi_n  = phi_n*norm;
            phi_p = ifftn(ek_p.*fp); norm = 1/sqrt(sum(sum(sum(abs(phi_p).^2)))*dV); phi_p  = phi_p*norm;
            mu_n = sum(sum(sum( h2_2m_n*eta_t*phi_n.*ifftn(K2.*fn) + U_n.*phi_n.^2 )))*dV;
            mu_p = sum(sum(sum( h2_2m_p*eta_t*phi_p.*ifftn(K2.*fp) + U_p.*phi_p.^2 )))*dV;
        end
        rho_p  = Zz*abs(phi_p).^2; rho_n  = Nn*abs(phi_n).^2; rho_0  = rho_n + rho_p; rho_1  = rho_n - rho_p;
        if e2 > 0
            rho3   = zeros(N3,N3,N3); rho3(1:Nx,1:Ny,1:Nz) = rho_p;
            rho3   = ifftn(Vq.*fftn(rho3)); Vc = rho3(1:Nx,1:Ny,1:Nz);
        end
        energy_vol = real(sum(sum(sum( 0.6*aN_n*rho_n.^(5/3) + 0.6*aN_p*rho_p.^(5/3) ...
            + at*rho_0.^2 + bt*rho_0.^gam + ct*rho_1.^2./(rho_0.^(1/3)+eps) + ct1*rho_1.^2 )))*dV);
        energy_Coul = real(sum(sum(sum(  Vc.*rho_p/2 - e2*0.75*(3/pi)^(1/3)*rho_p.^(4/3) )))*dV);
        fn   = fftn(phi_n); fp   = fftn(phi_p);
        energy_surf = eta_t*real(sum(sum(sum( ( h2_2m_n*Nn*conj(phi_n).*ifftn(K2.*fn) + h2_2m_p*Zz*conj(phi_p).*ifftn(K2.*fp) ) )))*dV );
        energy = energy_vol + energy_surf + energy_Coul;
        disp([toc, itime, hbarc_t*tcur, energy/Aa, energy_Coul/Aa, energy_surf/Aa, energy_vol/Aa, mu_n, mu_p])
    end
    % Liquid drop energy
    e_vol = -15.8*Aa; e_surf = 18.3*Aa^(2/3); e_Coul = 0.714*Zz^2/Aa^(1/3); e_sym = 23.3*(Nn-Zz)^2/Aa;
    e_mother = e_vol + e_sym + e_Coul + e_surf;
    e_dauthers = 2*( -15.8*Aa/2 + 18.3*(Aa/2)^(2/3) + 0.741*(Zz/2)^2/(Aa/2)^(1/3) + 23.3*(Nn/2-Zz/2)^2/(Aa/2) );
    e_barrier = (e_dauthers - e_mother) + e2*(Zz/2)^2/(2.4*(Aa/2)^(1/3));
    disp( [(e_vol + e_sym), e_surf, e_Coul, (e_vol + e_sym + e_Coul + e_surf), e_barrier, (e_dauthers - e_mother)] )
    disp([ energy_vol, energy_surf, energy_Coul, energy ])
    %===================================================================================
    Energy_gs = energy;
    save U238_gs.mat
elseif isave == 1
    load U238_gs.mat;
end

%=====================================================================================================================================

%-----------------------------------------------------------
% give the nucleus an initial kick and follow in time
%-----------------------------------------------------------
alpha_2 = 0.35*min([kz(2),kx(2),ky(2)].^2);  %alpha_2 = 0.00342694597260047;
alpha_3 = 0.75*alpha_2*0;
theta_n = alpha_2*(2*Z.^2 - X.^2 - Y.^2) + alpha_3/15*(5*Z.^2-3).*Z; theta_p = theta_n;
phi_n = phi_n.*exp(1i*theta_n); phi_p = phi_p.*exp(1i*theta_p);
%------------------------------------------------------------------------------------------------------------
% put CM velocity to zero
%------------------------------------------------------------------------------------------------------------
fn = fftn(phi_n); fp = fftn(phi_p);
vx = real( sum(sum(sum( Nn*conj(phi_n).*ifftn(Kx.*fn) + Zz*conj(phi_p).*ifftn(Kx.*fp) )))*dV)*hbarc_t/m_N/Aa;
vy = real( sum(sum(sum( Nn*conj(phi_n).*ifftn(Ky.*fn) + Zz*conj(phi_p).*ifftn(Ky.*fp) )))*dV)*hbarc_t/m_N/Aa;
vz = real( sum(sum(sum( Nn*conj(phi_n).*ifftn(Kz.*fn) + Zz*conj(phi_p).*ifftn(Kz.*fp) )))*dV)*hbarc_t/m_N/Aa;
theta_cm = (vx*X+vy*Y+vz*Z)*m_N/hbarc_t; phi_n = phi_n.*exp(-1i*theta_cm); phi_p = phi_p.*exp(-1i*theta_cm);
%-------------------------------------------------------------------------------------------------------------
tau = tau/1;  itau = 1i*tau; % time-step should be smaller for real-time evolution
tfin   = 1250/hbarc_t; nwrite = 1000; nsteps = ceil(tfin/nwrite/tau); itime = 0;  fftw('planner','measure');
icur   = 0; tcur = tau*icur;
Time = zeros(nwrite,1); Q20 = Time; q20 = Time;
E_vol = Time; E_W = Time; E_Coul = Time;  E_def = Time; E_flow = Time;
x1 = Time; y1 = Time; z1 = Time; x2 = Time; y2 = Time; z2 = Time;
%-------------------------------------------------------------------------------------------------
tic
energy_vol = real(sum(sum(sum( (0.6*aN_n*rho_n.^(5/3) + 0.6*aN_p*rho_p.^(5/3)) ...
    + at*rho_0.^2 + bt*rho_0.^gam + ct*rho_1.^2./(rho_0.^(1/3)+eps) + ct1*rho_1.^2 )))*dV);
energy_Coul = real(sum(sum(sum(  Vc.*rho_p/2 - e2*0.75*(3/pi)^(1/3)*rho_p.^(4/3) )))*dV);

fn   = fftn(phi_n); fp   = fftn(phi_p);
energy_kin = real(sum(sum(sum( ...
    ( h2_2m_n_t*Nn*conj(phi_n).*ifftn(K2.*fn) + h2_2m_p_t*Zz*conj(phi_p).*ifftn(K2.*fp) ) )))*dV );
energy = energy_vol + energy_Coul + energy_kin;

fn = sqrt(rho_n); fp = sqrt(rho_p);
energy_W    = real( sum(sum(sum( h2_2m_n_t*fn.*ifftn(K2.*fftn(fn)) + h2_2m_p_t*fp.*ifftn(K2.*fftn(fp)) )))*dV );
energy_flow = energy_kin - energy_W;
energy_def  = energy - energy_flow - Energy_gs;
im = 0;
disp([toc, toc, itime, hbarc_t*tcur, energy, energy_Coul, energy_W, energy - Energy_gs, energy_flow, energy_def])
dt = toc;
Frames = zeros(nwrite/5,Ny,Nz);
for it = 1:nwrite
    
    fn = fftn(phi_n); fp = fftn(phi_p);
    ek_n = exp( -itau*h2_2m_n_t*K2/2 );   ek_p = exp( -itau*h2_2m_p_t*K2/2 ); % half step
    phi_n = ifftn(ek_n.*fn);  phi_p = ifftn(ek_p.*fp);
    
    for kt = 1:nsteps
        itime = itime + 1; tcur = itime*tau;
        %-------------- Potentials -----------------------------------------
        rho_p  = Zz*abs(phi_p).^2; rho_n  = Nn*abs(phi_n).^2;
        rho_0  = rho_n + rho_p; rho_1  = rho_n - rho_p;
        U_p    = 2*at*rho_0 + gam*bt*rho_0.^(gam-1) - ct/3*rho_1.^2./rho_0.^(4/3);
        U_n    = U_p + 2*ct*rho_1./rho_0.^(1/3) + 2*ct1*rho_1 + aN_n*rho_n.^(2/3);
        U_p    = U_p - 2*ct*rho_1./rho_0.^(1/3) - 2*ct1*rho_1 + aN_p*rho_p.^(2/3);
        if e2 > 0
            %-------------   Coulomb potential ---------------------------------
            rho3   = zeros(N3,N3,N3); rho3(1:Nx,1:Ny,1:Nz) = rho_p;
            rho3   = ifftn(Vq.*fftn(rho3)); Vc = rho3(1:Nx,1:Ny,1:Nz);
            U_p    = U_p + Vc - e2*(3/pi)^(1/3)*rho_p.^(1/3);
        end
        %-------------------------------------------------------------------
        phi_n = exp( -itau*U_n ).*phi_n; phi_p = exp( -itau*U_p ).*phi_p;
        %-------------------------------------------------------------------
        fn = fftn(phi_n); fp = fftn(phi_p);
        if     kt < nsteps;  ek_n = exp( -itau*h2_2m_n_t*K2   );   ek_p = exp( -itau*h2_2m_p_t*K2   ); % full step
        elseif kt == nsteps; ek_n = exp( -itau*h2_2m_n_t*K2/2 );   ek_p = exp( -itau*h2_2m_p_t*K2/2 );% half step
        end
        phi_n = ifftn(ek_n.*fn); norm = 1/sqrt(sum(sum(sum(abs(phi_n).^2)))*dV); phi_n  = phi_n*norm;
        phi_p = ifftn(ek_p.*fp); norm = 1/sqrt(sum(sum(sum(abs(phi_p).^2)))*dV); phi_p  = phi_p*norm;
    end
    rho_p  = Zz*abs(phi_p).^2; rho_n  = Nn*abs(phi_n).^2;
    rho_0  = rho_n + rho_p; rho_1  = rho_n - rho_p;
    if e2 > 0
        rho3   = zeros(N3,N3,N3); rho3(1:Nx,1:Ny,1:Nz) = rho_p;
        rho3   = ifftn(Vq.*fftn(rho3)); Vc = rho3(1:Nx,1:Ny,1:Nz);
    end
    
    energy_vol = real(sum(sum(sum( (0.6*aN_n*rho_n.^(5/3) + 0.6*aN_p*rho_p.^(5/3)) ...
        + at*rho_0.^2 + bt*rho_0.^gam + ct*rho_1.^2./(rho_0.^(1/3)+eps) + ct1*rho_1.^2 )))*dV);
    energy_Coul = real(sum(sum(sum(  Vc.*rho_p/2 - e2*0.75*(3/pi)^(1/3)*rho_p.^(4/3) )))*dV);
    fn   = fftn(phi_n); fp   = fftn(phi_p);
    energy_kin = real(sum(sum(sum( ...
        ( h2_2m_n_t*Nn*conj(phi_n).*ifftn(K2.*fn) + h2_2m_p_t*Zz*conj(phi_p).*ifftn(K2.*fp) ) )))*dV );
    energy = energy_vol + energy_Coul + energy_kin;
    
    fn = sqrt(rho_n); dd_srho =           h2_2m_n_t*fn.*ifftn(K2.*fftn(fn));
    fp = sqrt(rho_p); dd_srho = dd_srho + h2_2m_p_t*fp.*ifftn(K2.*fftn(fp));
    energy_W    = real( sum(sum(sum( dd_srho )))*dV );
    energy_flow = energy_kin - energy_W;
    energy_def  = energy - energy_flow - Energy_gs;
    
    disp([toc, toc - dt, itime, hbarc_t*tcur, energy, energy_Coul, energy_W, energy - Energy_gs, energy_flow, energy_def])
    dt = toc;
    
    icur = icur + 1;
    Q20(icur)  = sum(sum(sum(rho_0.*(2*Z.^2-X.^2-Y.^2))))*dV/Aa; q20(icur)  = sum(sum(sum(rho_1.*(2*Z.^2-X.^2-Y.^2))))*dV/Aa;
    Time(icur) = hbarc_t*tcur; E_def(icur) = energy_def; E_flow(icur) = energy_flow;
    E_vol(icur) = energy_vol; E_def(icur) = energy_def; E_W(icur) = energy_W; E_Coul(icur) = energy_Coul;
    if mod(icur,5) == 0
    %-----------  plotting only ----------------------
    figure(11); set(gca,'FontSize',20); colormap('hot')
    pcolor(squeeze(Y(Nx/2,:,:)),squeeze(Z(Nx/2,:,:)),squeeze(rho_0(Nx/2,:,:)) )
    shading interp; colorbar; axis equal; xlim([-Ly/2,Ly/2]); ylim([-Lz/2,Lz/2]); 
    xlabel('Y [fm]'); ylabel('Z [fm]'); title('isoscalar density','FontSize',27); drawnow
    im = im + 1; Frames(im,:,:) = squeeze(rho_0(Nx/2,:,:));
    figure(12); set(gca,'FontSize',20); colormap('copper')
    surf(squeeze(Y(Nx/2,:,:)),squeeze(Z(Nx/2,:,:)),squeeze(rho_0(Nx/2,:,:)) )
    shading interp; colorbar;  xlim([-Ll/2,Ll/2]); ylim([-Ll/2,Ll/2]); zlim([0,0.18]);  
    xlabel('Y [fm]'); ylabel('Z [fm]'); title('isoscalar density','FontSize',27); drawnow
   
    figure(31); set(gca,'FontSize',20); plot(Time(1:icur),Q20(1:icur),'b','LineWidth',2); 
    hold on; xlabel('Time [fm/c]'); ylabel('Q_{20} [fm^2]'); grid on; drawnow

    figure(33); set(gca,'FontSize',20); plot(Time(1:icur),E_def(1:icur),'b','LineWidth',2); 
    hold on; xlabel('Time [fm/c]'); ylabel('E_{def} [MeV]'); grid on; drawnow
    
    figure(35); set(gca,'FontSize',20); plot(Q20(1:icur),E_def(1:icur),'b','LineWidth',2); 
    hold on; xlabel('Q_{20} [fm^2]'); ylabel('E_{def} [MeV]');grid on; drawnow
 
    end
    
end
%===================================================================================
save teqf.mat Time E_def Q20 Frames
for im = 1:200
    figure(11); set(gca,'FontSize',20); colormap('hot')
    RHO =   squeeze(Frames(im,:,:))';
    pcolor(  RHO')
    shading interp;  zlim([0, 0.14]); axis image off
    % colorbar; axis equal; xlim([-Ly/2,Ly/2]); ylim([-Lz/2,Lz/2]); drawnow
    % xlabel('Y [fm]'); ylabel('Z [fm]'); title('Fissioning ^{238}U - isoscalar density','FontSize',27)
    Mm(im) = getframe(11);

end

