"""Independent SciPy forward integral; no browser code or interpolation tables."""
import numpy as np
from scipy.integrate import quad
from scipy.special import wofz
from numpy.polynomial.legendre import leggauss

MPC_CM=3.0856775814913673e24

def normalization(configuration):
    factors=[]
    for amplitude,harmonic,shift in zip(configuration['fundamentalAmplitudes'],configuration['secondHarmonicAmplitudes'],configuration['secondHarmonicPhases']):
        factor,error=quad(lambda phase:np.exp(amplitude*np.cos(phase)+harmonic*np.cos(2*phase+shift)),0,2*np.pi,epsabs=1e-12,epsrel=1e-12)
        assert error<1e-9
        factors.append(factor/(2*np.pi))
    return np.prod(factors)

def neutral_and_velocity(configuration,settings,positions,density_normalization):
    basis=np.asarray(configuration['basis'],dtype=float)
    phase=2*np.pi/configuration['boxComovingMpc']*(basis@positions)+np.asarray(configuration['fundamentalPhases'])[:,None]
    fundamental=np.asarray(configuration['fundamentalAmplitudes'])[:,None]
    harmonic=np.asarray(configuration['secondHarmonicAmplitudes'])[:,None]
    shift=np.asarray(configuration['secondHarmonicPhases'])[:,None]
    density=configuration['meanHydrogenPerCm3']*np.exp(np.sum(fundamental*np.cos(phase)+harmonic*np.cos(2*phase+shift),axis=0))/density_normalization
    temperature=settings['temperatureKelvin']
    lower=np.sqrt(temperature/3.148);upper=np.sqrt(temperature/7.036e5)
    recombination=7.982e-11/(lower*(1+lower)**.252*(1+upper)**1.748)
    ratio=recombination*density/settings['photoionizationPerSecond']
    fraction=2*ratio/(1+2*ratio+np.sqrt(1+4*ratio))
    expansion=configuration['hubbleTodayKmPerSecondPerMpc']*np.sqrt(configuration['omegaMatter']*(1+configuration['redshift'])**3+configuration['omegaLambda'])/(1+configuration['redshift'])
    wave=2*np.pi/configuration['boxComovingMpc']
    peculiar=-settings['velocityScale']*configuration['growthRateApproximation']*expansion/wave*np.sum((basis[:,0]/np.sum(basis*basis,axis=1))[:,None]*(fundamental*np.sin(phase)+harmonic/2*np.sin(2*phase+shift)),axis=0)
    return density*fraction,expansion*positions[0]+peculiar

def spectrum(configuration,settings,cells,order):
    if not (5000<=settings['temperatureKelvin']<=30000 and .5e-12<=settings['photoionizationPerSecond']<=2e-12 and 0<=settings['velocityScale']<=2):raise ValueError('Outside declared screen domain')
    if any(abs(settings[key])>4 for key in ['transverseYMpc','transverseZMpc']):raise ValueError('Outside finite transverse domain')
    nodes,weights=leggauss(order);width=configuration['boxComovingMpc']/cells
    positions=(-4+(np.arange(cells)[:,None]+.5)*width+width/2*nodes).ravel()
    proper_weights=np.tile(weights,cells)*width/2*MPC_CM/(1+configuration['redshift'])
    points=np.array([positions,np.full(len(positions),settings['transverseYMpc']),np.full(len(positions),settings['transverseZMpc'])])
    neutral,velocity=neutral_and_velocity(configuration,settings,points,normalization(configuration))
    column=np.dot(neutral,proper_weights)
    assert column*6.3e-18<.1,'Optically thin ionizing-screen assumption failed'
    thermal=np.sqrt(2*1.380649e-16*settings['temperatureKelvin']/1.67262192595e-24)
    damping=6.265e8*1215.67e-8/(4*np.pi*thermal)
    strength=np.pi*2.8179403205e-13*2.99792458e10*.4162*1215.67e-8
    grid=np.linspace(-500,500,601);depth=[]
    for observed in grid:
        cross_section=strength/(np.sqrt(np.pi)*thermal)*wofz((observed-velocity)*1e5/thermal+1j*damping).real
        depth.append(np.dot(neutral*proper_weights,cross_section))
    depth=np.asarray(depth)
    return grid,depth,np.exp(-depth),float(column)