import numpy as np
k=np.array([0.64,1.64,2.64])     # 2k2-k1=k3 exact
# density amplitudes drho ~ k^2*Phi_weyl, ratios [1,0.505,0.447]; acoustic signs compression/rarefaction
dr=k**2*np.array([1.0,0.505,0.447]); dr=dr/dr[0]
N=2**19; L=2*np.pi*60/k[0]; x=np.linspace(0,L,N,endpoint=False); dx=x[1]-x[0]
def dmond_phi_k3(amps,signs):
    src=sum(signs[i]*amps[i]*np.cos(k[i]*x) for i in range(len(amps)) if amps[i]!=0)
    G1=np.cumsum(src)*dx; G1-=G1.mean()                  # int rho dx
    Pp=np.sign(G1)*np.sqrt(np.abs(G1)); PM=np.cumsum(Pp)*dx; PM-=PM.mean()
    c=2*np.sum(PM*np.cos(k[2]*x))*dx/L; s=2*np.sum(PM*np.sin(k[2]*x))*dx/L
    return np.hypot(c,s)
for sgnlabel,signs in [("[+,-,+]",[1,-1,1]),("[+,+,+]",[1,1,1]),("[+,-,-]",[1,-1,-1])]:
    full=dmond_phi_k3([dr[0],dr[1],dr[2]],signs)
    k3only=dmond_phi_k3([0,0,dr[2]],[0,0,signs[2]])
    print(f"  signs {sgnlabel}: Phi_M(k3) full(k1,k2,k3)={full:.4f}  k3-only={k3only:.4f}  FEEDING factor={full/k3only:.3f}")
print("\n  (>1 = peaks1,2 constructively boost peak3 ; ~1 = no feeding ; <1 = destructive)")
print("  density amps used:",np.round(dr,2))
