"""Independent smooth conservative Fourier–Galerkin Euler–Poisson reference.

This has no finite-volume reconstruction or Riemann flux. The2/3cutoff suppresses
aliased products; before shocks it checks the finite-volume solution. It is not
used as a nonlinear-shock reference or a mass-spectrum prediction.
"""
import numpy as np


def integrate_spectral(initial,sound,end_time,step=.005):
    state=initial.copy();n=state.shape[-1]
    modes=np.fft.fftfreq(n)*n
    kx,ky=np.meshgrid(2*np.pi*modes,2*np.pi*modes,indexing='ij')
    mask=(abs(kx)<2*np.pi*n/3)&(abs(ky)<2*np.pi*n/3)
    squared=kx**2+ky**2;inverse=np.divide(-1.,squared,out=np.zeros_like(squared),where=squared!=0)
    def rhs(value):
        density=value[0]
        if density.min()<=0:raise FloatingPointError('Spectral reference lost positive density')
        vx,vy=value[1]/density,value[2]/density
        flux_x=np.array([value[1],value[1]*vx+sound**2*density,value[2]*vx])
        flux_y=np.array([value[2],value[1]*vy,value[2]*vy+sound**2*density])
        derivative=-np.fft.ifft2((1j*kx*np.fft.fft2(flux_x)+1j*ky*np.fft.fft2(flux_y))*mask).real
        potential=np.fft.fft2(density-density.mean())*inverse
        gx=np.fft.ifft2(-1j*kx*potential).real;gy=np.fft.ifft2(-1j*ky*potential).real
        derivative[1]+=np.fft.ifft2(np.fft.fft2(density*gx)*mask).real
        derivative[2]+=np.fft.ifft2(np.fft.fft2(density*gy)*mask).real
        return derivative
    time=0.
    while time<end_time-1e-14:
        dt=min(step,end_time-time)
        a=rhs(state);b=rhs(state+dt*a/2);c=rhs(state+dt*b/2);d=rhs(state+dt*c)
        state+=dt*(a+2*b+2*c+d)/6;time+=dt
    return state


def spectral_cell_averages(point_values):
    n=point_values.shape[-1];modes=np.fft.fftfreq(n)*n
    window=np.sinc(modes[:,None]/n)*np.sinc(modes[None,:]/n)
    return np.fft.ifft2(np.fft.fft2(point_values)*window).real
