VMS_SPARK_README.txt

VMS + NFW + MOND analysis for SPARC rotation curves
(Nov 2025 fixed version Ð units clean, fits honest, no magic constants hiding in the walls)

0. What this script is

This is a single-file, no-nonsense comparison between three rotation-curve models on SPARC galaxies:
* VMS: my geometric residue profile (effectively a Zhao n=0 spherical halo),
with one free parameter per galaxy: ?? (a density scale in M?/kpc?).
* NFW: standard cold dark matter halo, but constrained to one free parameter per galaxy: ?_s.
The scale radius r_s is fixed per galaxy from a simple heuristic, on purpose, so VMS and NFW are on equal footing.
* MOND: classic Milgrom MOND with fixed
a? = 1.2 ? 10??? m/s?,
so zero free parameters per galaxy in this implementation.

All three models use the same baryonic baseline from SPARC:
* Vgas, Vdisk, Vbul are combined as
Vbar = sqrt(Vgas? + Vdisk? + Vbul?).
* MOND modifies that baryonic curve directly.
* VMS and NFW add a dark/geometric halo on top of the same Vbar:

v_tot?(r) = v_bar?(r) + v_DM?(r)

so nobody gets a special advantage or handicap on the baryons.

The code uses standard SPARC units:
* Radius r in kpc
* Velocity v in km/s
* Densities in M?/kpc?
* G in kpc (km/s)? / M? (and we do not Òre-convertÓ it into some Frankenstein unit system)

If youÕve ever looked at SPARC, this will feel familiar.

1. What the script does, in plain terms

From top to bottom, the script:
1. Defines constants and unit system
o G_KPC = 4.302e-6 in kpc (km/s)? / M? Ð the standard astrophysical G.
o A0 = 1.2e-10 m/s? Ð the standard MOND acceleration scale.
2. Loads SPARC rotation curve files
o Looks for a file matching ./sparc_data/*{GALAXY}*.dat.
o Handles both SPARC formats:
* 5 columns: Rad, Vobs, errV, Vgas, Vdisk ? sets Vbul = 0.0.
* 6 columns: Rad, Vobs, errV, Vgas, Vdisk, Vbul.
o Drops NaNs and resets the index.
o Prints how many data points it found for that galaxy.
3. Implements the three models
o VMS:
v?(r) = 4 ¹ G ?? r / (r + r_s)
vms_velocity(r, rho0, rs) returns v_DM(r) in km/s.
o NFW:
* Mass:
M(<r) = 4 ¹ ?_s r_s? [ln(1 + x) ? x/(1 + x)], where x = r/r_s
* Velocity:
v_DM?(r) = G M(<r) / r
Implemented as nfw_mass + nfw_velocity.
A tiny floor on r (r_safe = max(r, 1e-10)) avoids divide-by-zero if a radius is exactly zero.
o MOND:
* Uses the simple relation v? = Vbar? a? r.
* Converts r to metres, Vbar to m/s, applies the formula, and converts back to km/s.
* Implemented in mond_velocity(r, Vbar).
4. Computes statistics
o Residuals: residuals = v_obs ? v_model.
o ??: ?[(residual/?)?].
o DOF: N ? n_params.
o Reduced ??: ??_red = ??/DOF (or NaN if DOF ² 0).
o RMS: square root of the mean squared residual, in km/s.
o This is all wrapped in calculate_statistics.
5. Fits VMS and NFW (one parameter each)
o For VMS:
* Fix r_s per galaxy (rs), based on the galaxyÕs radial extent.
* Use the asymptotic relation v_flat? = 4 ¹ G ?? r_s to guess ??.
* Fit only ?? with curve_fit, using v_tot(r) = sqrt(v_DM? + v_bar?).
* absolute_sigma=True so the parameter error is a proper 1? uncertainty based on the actual data errors.
o For NFW:
* Fix r_s per galaxy again (same rs), so weÕre only fitting ?_s.
* Use a heuristic asymptotic relation to guess ?_s (rough but good enough).
* Fit only ?_s with curve_fit, again using v_tot? = v_bar? + v_DM?.
6. Runs a single-galaxy test at the bottom
o You specify a galaxy label (e.g. ÔUGC02259Õ).
o It:
* Loads the data.
* Builds v_bar.
* Floors errV at 0.5 km/s so we donÕt let one ÒperfectÓ point dominate the fit.
* Computes MOND prediction (0 parameters).
* Sets v_flat from the last few observed velocities.
* Sets rs_guess = max(r)/3 as a reasonable VMS/NFW scale radius.
* Fits VMS (??) and NFW (?_s).
* Prints ??_red and RMS for all three models, plus the fitted densities.

2. How to run it

2.1 Install dependencies

You need Python 3.8+ with:
* numpy
* pandas
* scipy

If youÕre in a typical scientific environment, you probably already have these.

2.2 Get the SPARC data
1. Go to:
http://astroweb.cwru.edu/SPARC/
2. Download the rotation curve / Rotmod data.
3. Extract the .dat files into a folder called ./sparc_data/ right next to this script.

You should end up with paths like:
* ./sparc_data/UGC02259_Rotmod.dat
* ./sparc_data/NGC3198_Rotmod.dat
* etc.

2.3 Run a single galaxy

Open the script and look at the bottom:
if __name__ == '__main__':
    galaxy = 'UGC02259'          # change as needed
    df = load_sparc_galaxy(galaxy)
    ...
Change 'UGC02259' to whatever SPARC galaxy you want, e.g.:
galaxy = 'NGC3198'
Then from a terminal:
python vms_analysis.py
YouÕll see output like:
* MOND ??_red and RMS
* VMS ??_red, RMS, and ?? ± error
* NFW ??_red, RMS, and ?_s ± error

If youÕve ever fit rotation curves before, youÕll know how to read those immediately.

3. Caveats, gotchas, and deliberate choices

IÕm not hiding the compromises here; theyÕre baked in on purpose:
1. One-parameter halos for VMS and NFW
o In this particular script, r_s is not fitted. ItÕs fixed from the galaxyÕs size (max radius / 3).
o That means VMS and NFW each get only 1 free parameter per galaxy: ?? or ?_s.
o MOND gets 0 (a? fixed globally).
o This is about putting everyone on a minimalist, symmetric footing, not about squeezing every last bit of ?? performance out of NFW.
2. Error floor on v_err
o v_err = max(errV, 0.5 km/s) is there so that one datapoint with a vanishingly small formal error doesnÕt blow up ?? or dominate the fit.
o If you hate that choice, change 0.5 to something else. The code wonÕt fight you.
3. Heuristic NFW initial guess
o The ?_s initial guess is based on a rough asymptotic scaling with v_flat and r_s.
o ItÕs not trying to solve the full NFW mass profile analytically; itÕs just good enough to let curve_fit converge fast.
4. Guard rails for pathological data
o NFW uses r_safe = max(r, 1e-10) to avoid dividing by zero.
o If your SPARC subset has some weird points, the code will try to behave instead of detonating.
5. MOND implementation is simple on purpose
o No fancy interpolating function, no ?(x) logistics Ñ just the simple v? = v_bar? a? r relation.
o If you want a more sophisticated MOND flavor, you can absolutely extend it.

Overall: this is meant to be a clean, honest, low-parameter comparison, not a maximal-parameter, galaxy-by-galaxy tuning exercise.

4. Optional: go crazy on many galaxies (batch wrapper)

If you want to run this on dozens or hundreds of SPARC galaxies, you donÕt need to rewrite everything. You just need a thin wrapper that loops over a list of galaxy names, runs the same logic, and writes a CSV.

HereÕs a simple multi-galaxy runner you can bolt onto the bottom of the file.

How to use this wrapper
* Option A (replace):
Replace the existing if __name__ == '__main__': block with the one below.
* Option B (switch modes):
Keep both, but comment out the single-galaxy block when you want to run the batch.

Either way, hereÕs the wrapper:
# =============================================================================
# Optional: multi-galaxy batch runner
# =============================================================================
if __name__ == '__main__':
    # Edit this list however you like. Use SPARC galaxy names that match your .dat files.
    galaxy_list = [
        'UGC02259',
        'NGC3198',
        'NGC2403',
        'NGC5055',
        'M33',
        'UGC2885'
        # ...add more SPARC IDs here...
    ]

    rows = []

    for galaxy in galaxy_list:
        print(f"\n=== {galaxy} ===")
        try:
            df = load_sparc_galaxy(galaxy)
            r = df['Rad'].values
            v_obs = df['Vobs'].values
            v_err = np.maximum(df['errV'].values, 0.5)  # floor on errors
            v_bar = np.sqrt(df['Vgas']**2 + df['Vdisk']**2 + df['Vbul']**2)

            # MOND (0 parameters)
            v_mond = mond_velocity(r, v_bar)
            mond_stats = calculate_statistics(v_obs, v_mond, v_err, n_params=0)

            # Flat velocity & scale radius guesses
            v_flat = np.median(v_obs[-5:]) if len(v_obs) >= 5 else np.max(v_obs)
            rs_guess = np.max(r) / 3.0

            # VMS (1 parameter: ??)
            rho0, rho0_err = fit_vms(r, v_obs, v_err, rs_guess, v_flat, v_bar)
            v_dm_vms = vms_velocity(r, rho0, rs_guess)
            v_vms = np.sqrt(v_dm_vms**2 + v_bar**2)
            vms_stats = calculate_statistics(v_obs, v_vms, v_err, n_params=1)

            # NFW (1 parameter: ?_s)
            rho_s, rho_s_err = fit_nfw(r, v_obs, v_err, rs_guess, v_flat, v_bar)
            v_dm_nfw = nfw_velocity(r, rho_s, rs_guess)
            v_nfw = np.sqrt(v_dm_nfw**2 + v_bar**2)
            nfw_stats = calculate_statistics(v_obs, v_nfw, v_err, n_params=1)

            rows.append({
                'Galaxy': galaxy,
                'N_points': len(r),

                'MOND_chi2_red': mond_stats['chi2_red'],
                'MOND_rms': mond_stats['rms'],

                'VMS_chi2_red': vms_stats['chi2_red'],
                'VMS_rms': vms_stats['rms'],
                'VMS_rho0': rho0,
                'VMS_rho0_err': rho0_err,

                'NFW_chi2_red': nfw_stats['chi2_red'],
                'NFW_rms': nfw_stats['rms'],
                'NFW_rho_s': rho_s,
                'NFW_rho_s_err': rho_s_err,
            })

        except Exception as e:
            print(f"Skipping {galaxy}: {e}")

    if rows:
        summary = pd.DataFrame(rows)
        summary.to_csv('sparc_vms_mond_nfw_summary.csv', index=False)
        print("\nSaved summary to sparc_vms_mond_nfw_summary.csv")
        print(summary)
    else:
        print("\nNo successful fits Ð check your galaxy names and data files.")

What this gives you
* One CSV: sparc_vms_mond_nfw_summary.csv
* Each row: one galaxy, with:
o Number of points.
o MOND ??_red, RMS.
o VMS ??_red, RMS, ?? ± err.
o NFW ??_red, RMS, ?_s ± err.

From there, you can:
* Make histograms of ??_red for each model.
* Look for outliers where VMS crushes NFW (or vice versa).
* Check correlations between ??, ?_s, and galaxy properties.

5. Sharing and publishing

If you do run this on a large SPARC sample and get interesting trends:
* YouÕre absolutely free to:
o Share the results anywhere.
o Publish your own analysis.
o Fork the script and modify it.

If you want to share your summary with me for inclusion on the VMS Institute site (or for people to compare their own runs against), feel free to send:
* Your sparc_vms_mond_nfw_summary.csv.
* A short note on:
o Which SPARC galaxies you used.
o Any changes you made to the script (error floors, r_s strategy, MOND variant, etc.).

The whole point of this is to make it easy and transparent for anyone to test VMS against the usual suspects, on public data, with code they can read in one sitting.