FRCMFD-v2: TEST 1C-N — NARROW-BAND RESONANCE DIAGNOSTIC (0.28–0.32 v, +0.35)

""" FRCMFD-v2: TEST 1C-N — NARROW-BAND RESONANCE DIAGNOSTIC (0.28–0.32 v, +0.35) Focus: Clean COM, dv/dt, canonical momentum P_z, and momentum drift With periodic z-boundaries, consistent operators, and refined timestep. """ import numpy as np import scipy.sparse as sp from datetime import datetime import json import glob import os from scipy.signal import savgol_filter # ================================================================ # LOAD SOLITON # ================================================================ preferred = "/content/test_0A_soliton_20260523_195958.npz" if os.path.exists(preferred): soliton_file = preferred else: soliton_files = glob.glob("test_0A_*_soliton.npz") if not soliton_files: soliton_files = glob.glob("/content/drive/MyDrive/FRCMFD_v2_Backups/*/test_0A_*_soliton.npz") soliton_file = sorted(soliton_files)[-1] data = np.load(soliton_file) Psi_soliton = data["Psi_soliton"] r_grid = data["r_grid"] z_grid = data["z_grid"] dr = float(data["dr"]) dz = float(data["dz"]) v = float(data["v"]) mu = float(data["mu"]) lam = float(data["lam"]) kappa = float(data["kappa"]) m = int(data["m"]) S_max = float(data["S_max"]) Psi_sat = float(data["Psi_sat"]) # Shape handling: enforce (nz, nr) ordering if Psi_soliton.ndim == 2: s0, s1 = Psi_soliton.shape nr_from_r = len(r_grid) nz_from_z = len(z_grid) if (s0, s1) == (nz_from_z, nr_from_r): nz, nr = s0, s1 Psi_base_2d = Psi_soliton.copy() else: nr, nz = s0, s1 Psi_base_2d = Psi_soliton.T.copy() else: nr = len(r_grid) nz = len(z_grid) Psi_base_2d = Psi_soliton.reshape((nz, nr)) Psi_flat = Psi_base_2d.reshape(-1) Lz = z_grid[-1] - z_grid[0] # ================================================================ # BUILD OPERATORS # ================================================================ def build_radial_operator(r_grid, dr): nr = len(r_grid) r_face = np.zeros(nr + 1) r_face[0] = r_grid[0] - dr/2 for i in range(1, nr + 1): r_face[i] = r_grid[i-1] + dr/2 flux_right = r_face[1:] / dr flux_left = r_face[:-1] / dr main_diag = -(flux_left + flux_right) lower_diag = flux_left[1:] upper_diag = flux_right[:-1] M = sp.diags([lower_diag, main_diag, upper_diag], [-1,0,1], format="csr") w_r = r_grid * dr W_r = sp.diags(w_r, format="csr") W_r_inv = sp.diags(1.0/w_r, format="csr") return W_r_inv @ M, W_r def build_axial_second_derivative(nz, dz): main_diag = np.full(nz, -2.0/dz**2) off_diag = np.full(nz-1, 1.0/dz**2) L_z = sp.diags([off_diag, main_diag, off_diag], [-1, 0, 1], format="lil") # periodic wrap: top <-> bottom L_z[0, nz-1] = 1.0/dz**2 L_z[nz-1, 0] = 1.0/dz**2 L_z = L_z.tocsr() W_z = sp.diags(np.ones(nz)*dz, format="csr") return L_z, W_z def build_axial_first_derivative(nz, dz): off_diag = np.full(nz-1, 1.0/(2.0*dz)) D_z = sp.diags([-off_diag, off_diag], [-1, 1], format="lil") # periodic wrap: top <-> bottom D_z[0, nz-1] = -1.0/(2.0*dz) D_z[nz-1, 0] = 1.0/(2.0*dz) return D_z.tocsr() L_r, W_r = build_radial_operator(r_grid, dr) L_z, W_z = build_axial_second_derivative(nz, dz) D_z = build_axial_first_derivative(nz, dz) I_r = sp.eye(nr, format="csr") I_z = sp.eye(nz, format="csr") L_2D = sp.kron(I_z, L_r) + sp.kron(L_z, I_r) W_2D = sp.kron(W_z, W_r) D_2D = sp.kron(D_z, I_r) dV = W_2D.diagonal() * 2*np.pi r_mesh_2d = np.tile(r_grid, nz) # ================================================================ # PHYSICS FUNCTIONS # ================================================================ def compute_energy(Psi): psi_sq = np.abs(Psi)**2 kin = -0.5*v**2 * np.real(np.sum(np.conj(Psi)*(L_2D@Psi)*dV)) mass = -0.5*mu*np.sum(psi_sq*dV) nl = 0.25*lam*np.sum(psi_sq*psi_sq*dV) S = S_max*np.tanh(psi_sq/(Psi_sat**2)) tension = 0.5*kappa*np.sum(S*psi_sq*dV) cent = 0.5*v**2*m**2*np.sum(psi_sq/(r_mesh_2d**2+1e-12)*dV) return (kin+mass+nl+tension+cent).real def acceleration(Psi): psi_sq = np.abs(Psi)**2 S = S_max*np.tanh(psi_sq/(Psi_sat**2)) dS = (S_max/(Psi_sat**2))*(1/np.cosh(psi_sq/(Psi_sat**2))**2) term_kin = -v**2*(L_2D@Psi) term_mass = mu*Psi term_nl = lam*psi_sq*Psi term_tension = kappa*(S + psi_sq*dS)*Psi term_cent = v**2*m**2*Psi/(r_mesh_2d**2+1e-12) return -(term_kin + term_mass + term_nl + term_tension + term_cent) def compute_z_com(Psi): Psi_2d = Psi.reshape((nz,nr)) dens = np.abs(Psi_2d)**2 dV_r = 2*np.pi*r_grid*dr*dz z_profile = np.sum(dens*dV_r, axis=1) return np.sum(z_grid*z_profile)/np.sum(z_profile) def compute_momentum_z(Psi): # Use the same D_2D operator family for consistency Psi_2d = Psi.reshape((nz,nr)) dPsi_flat = D_2D @ Psi dPsi_2d = dPsi_flat.reshape((nz,nr)) integrand = np.imag(np.conj(Psi_2d)*dPsi_2d) R = np.tile(r_grid, (nz,1)) dV_loc = 2*np.pi*R*dr*dz return float(np.sum(integrand*dV_loc)) def compute_asymmetry(Psi): Psi_2d = Psi.reshape((nz,nr)) mid = nr//2 phase = np.unwrap(np.angle(Psi_2d[:,mid])) grad = np.abs(np.gradient(phase, dz)) amp = np.abs(Psi_2d[:,mid]) c = np.argmax(amp) m = min(20, nz//10) if cnz-m: return 0.0 f = np.mean(grad[c:c+m]) r = np.mean(grad[c-m:c]) return f/max(r,1e-10) def compute_wake_power(Psi): Psi_2d = Psi.reshape((nz,nr)) R = np.tile(r_grid, (nz,1)) mask = R>3.0 phase = np.unwrap(np.angle(Psi_2d), axis=0) return float(np.var(phase[mask])) def unwrap(z): out = z.copy() off = 0 for i in range(1,len(z)): dz_loc = z[i]-z[i-1] if dz_loc>0.5*Lz: off -= Lz elif dz_loc<-0.5*Lz: off += Lz out[i] += off return out def vel_accel(t,z): dt_loc = np.mean(np.diff(t)) v_loc = np.gradient(z,dt_loc) v_s = savgol_filter(v_loc,9,2) if len(v_loc)>=9 else v_loc a_loc = np.gradient(v_s,dt_loc) return v_s,a_loc # ================================================================ # NARROW-BAND VELOCITIES # ================================================================ vels = [0.28, 0.30, 0.32, 0.35] results = {} dt = 0.0005 t_max = 100.0 n_steps = int(t_max/dt) n_save = 400 # ~0.2 time units between diagnostics for vf in vels: vs = vf*v # Start from base soliton in (nz, nr) Psi_2d = Psi_base_2d.copy() z_mesh = np.tile(z_grid.reshape(-1,1),(1,nr)) Psi_2d = Psi_2d * np.exp(1j*(vs/v)*z_mesh) Psi = Psi_2d.reshape(-1) # Initialize velocity using consistent D_2D operator Psi_dot = np.real(-vs * (D_2D @ Psi)) P0 = compute_momentum_z(Psi) t_list = [] z_list = [] P_list = [] asym = [] wake = [] amps = [] E_list = [] for step in range(n_steps): acc = acceleration(Psi) P_half = Psi_dot + 0.5*dt*acc Psi_new = Psi + dt*P_half acc2 = acceleration(Psi_new) Psi_dot_new = P_half + 0.5*dt*acc2 Psi, Psi_dot = Psi_new, Psi_dot_new if step % n_save == 0: t = step*dt t_list.append(t) z_list.append(compute_z_com(Psi)) P_list.append(compute_momentum_z(Psi)) asym.append(compute_asymmetry(Psi)) wake.append(compute_wake_power(Psi)) amps.append(np.max(np.abs(Psi))) E_list.append(compute_energy(Psi)) t_arr = np.array(t_list) z_arr = np.array(z_list) P_arr = np.array(P_list) z_un = unwrap(z_arr) v_inst, a_inst = vel_accel(t_arr, z_un) mask = (t_arr>=15) & (t_arr<=80) mean_v = float(np.mean(v_inst[mask])) mean_a = float(np.mean(a_inst[mask])) mean_P = float(np.mean(P_arr[mask])) dPdt = np.gradient(P_arr[mask], np.mean(np.diff(t_arr[mask]))) mean_dP = float(np.mean(dPdt)) Pf = P_arr[-1] drift = (Pf-P0)/abs(P0)*100 if abs(P0)>1e-15 else 0.0 results[vf] = { "mean_asymmetry": float(np.mean(asym)), "mean_wake_power": float(np.mean(wake)), "mean_COM_velocity": mean_v, "mean_dvdt": mean_a, "mean_Pz": mean_P, "mean_dPz_dt": mean_dP, "momentum_initial": float(P0), "momentum_final": float(Pf), "momentum_drift_percent": float(drift), "energy_drift_percent": float((E_list[-1]-E_list[0])/abs(E_list[0])*100), "amplitude_preservation_percent": float(amps[-1]/amps[0]*100), } # ================================================================ # SAVE RESULTS # ================================================================ ts = datetime.now().strftime("%Y%m%d_%H%M%S") name = f"test_1C_NarrowBand_periodic_consistent_{ts}" with open(f"{name}.json","w") as f: json.dump(results,f,indent=2) print("✓ Narrow-band Test 1C-N (periodic z, consistent operators, dt=0.0005) complete.") print(f"✓ Results saved: {name}.json")

Popular posts from this blog

THE GOLDEN BALLROOM/BUNKER

Image
Ben Meiselas reports on the shocking admission by Donald Trump’s DOJ in a court case where the DOJ admits to a secret project underneath the ballroom which they claim is needed to protect Donald Trump’s life for “national security purposes.” "You unlock this door with the key of complicity. Beyond it is another dimension — a dimension of betrayal, of indulgence, of fear. You’re moving into a land of both shadow and substance, of politics and paranoia. You’ve just crossed into… the MAGA Zone." "Tonight’s story: A leader sworn to protect his nation makes a bargain with its enemies. The deal? Silence in the face of nuclear annihilation. No retaliation, no defense — only surrender dressed in secrecy. While citizens live unaware, their president builds a palace beneath the earth, a ballroom of gold, of marble and chandeliers, a masquerade hall for billionaires. But behind the gilded doors lies not music and laughter, but a bomb shelter — a sanctuary for the few, pur...
Read more

Conceptual Summary #2: (∂t2​S−c2∇2S+βS3)=σ(x,t)⋅FR​(C[Ψ])

Image
The tests designed to find the Aether—principally the Michelson-Morley experiment (M-M)—would be largely ineffective at detecting the fundamental S field (Substrate) in the Reactive Substrate Theory (RST) framework, though their null results are highly compatible with RST's structure. The reason lies in the distinct conceptual definition and behavior of the two media. The Test: Michelson-Morley and the Aether. The M-M experiment was designed to detect the "aether wind," a predicted change in the speed of light caused by the Earth's motion relative to a stationary, rigid medium—the classical Luminiferous Aether. Classical Aether: A static, absolute reference frame; an elastic, passive, material-like medium that fills space and serves as a carrier for light waves. Predicted Result: Light traveling parallel to Earth's motion through the Aether should be slightly slower than light traveling perpendicular to it. Actual Result (Null Result): No significant differenc...
Read more

ICE PROUDLY ANNOUNCES NEW “ELITE” TASK FORCE COMMANDER JEREMY DEWITTE

Image
For Immediate Release OFFICE OF STRATEGIC SUPPRESSION & DESPERATE RECRUITMENT U.S. IMMIGRATION AND CUSTOMS ENFORCEMENT (ICE) SUBJECT: ICE PROUDLY ANNOUNCES NEW “ELITE” TASK FORCE COMMANDER JEREMY DEWITTE; AGENCY TO ADOPT FULL “MOTOR ESCORT” DOCTRINE BY Q3 WASHINGTON, D.C. — In a historic moment that experts are calling “a cry for help” and insiders are calling “Tuesday,” U.S. Immigration and Customs Enforcement (ICE) has officially appointed Jeremy Dewitte as National Coordinator of the newly created Tactical High-Visibility Fallujah Squadron & Escort Wing (T.H.V.F.S.E.W.) . This appointment follows a PBS NewsHour report highlighting ICE’s bold new hiring strategy: removing every meaningful hiring standard. After eliminating outdated requirements such as "must not be a professional police impersonator," "must possess a valid driver's license," and "must not be a registered sex ...
Read more