"""High-speed addendum anchored to published Chiron stopping distances.

This is a distance-equivalent force-envelope model, not a Bugatti digital twin.
One aggregate parameter is fitted separately to each published starting-speed
test. The shared air-brake envelope is retained for both reference/prototype;
the new fan, seal, mass and tyre-load effects remain conditional assumptions.
"""
import copy,csv,hashlib,json,math
from pathlib import Path
import numpy as np
from scipy.optimize import brentq
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import vehicle
from ground_module import config as ground_config,steady

R=Path(__file__).resolve().parents[1]
BENCHMARKS={100:31.4,200:114.,300:247.,400:491.}
SPEC='https://bugatti-newsroom.imgix.net/66703700d9bf8f4b7ce9211c/211122_BU_Chiron%20ENG.pdf'
RECORD='https://newsroom.bugatti.com/en/press-releases/0-400-0-km-h-in-42-seconds-bugatti-chiron-sets-world-record'


def config(speed,area=2.):
    c=copy.deepcopy(vehicle.config())
    c.update(base_mass_kg=1995.,wheelbase_m=2.711,track_m=(1.749+1.661)/2,
             yaw_inertia_base_kgm2=2500.*1995/1500,
             initial_fan_spin_energy_J_assumed=8000.)
    c['airbrake_reference']={
        'initial_Cd':.36 if speed>=380 else .40,'airbrake_Cd':.60,
        'projected_frontal_area_m2_assumed':area,'deployment_duration_s_assumed':.8,
        'rear_downforce_increment_kg_equivalent_at_400kmh':900.,
        'downforce_application_x_m_assumed':-c['wheelbase_m']/2,
        'source':SPEC,'rear_load_and_deployment_source':RECORD,
        'scope':'Published Cd endpoints and nominal rear-load increment. Frontal area, linear deployment interpolation, force location, CG, axle loads and inertia are assumptions. Other aero/tyre/brake effects are absorbed into a separately fitted equivalent parameter.'}
    return c


def fit(c,speed,target):
    s={'name':'published-distance-calibration','surface_factors':[1,1]}
    def residual(mu):return vehicle.run(c,mu,s,'reference',speed,record=False)[0]['path_distance_m']-target
    mu=brentq(residual,.25,3.,xtol=3e-7)
    z,_=vehicle.run(c,mu,s,'reference',speed,record=False)
    assert abs(z['path_distance_m']-target)<.01
    return mu,z


def budget(z,h,c):
    mass=c['base_mass_kg']+(c['installed_hardware_mass_target_kg'] if z['mode']!='reference' else 0)
    ke=.5*mass*(z['speed_kmh']/3.6)**2
    drag=float(np.trapezoid([r['aerodynamic_drag_power_W'] for r in h],[r['time_s'] for r in h]))
    peak=max(max(0.,-mass*r['body_ax_ms2']*r['speed_ms']-r['aerodynamic_drag_power_W']) for r in h)
    return {'initial_kinetic_energy_MJ':ke/1e6,'aerodynamic_energy_dissipated_MJ':drag/1e6,
            'remaining_brake_tyre_energy_reference_MJ':(ke-drag)/1e6,
            'peak_combined_brake_tyre_power_reference_MW':peak/1e6,
            'scope':'Straight-line body-energy partition. No disc/pad/tyre temperature, heat split, fade or torque ceiling is predicted.'}


def save_rows(name,rows):
    with (R/'output'/name).open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)


def main():
    calibrations=[];cases=[];histories={};mus={}
    straight={'name':'uniform-high-speed','surface_factors':[1,1]}
    for speed,target in BENCHMARKS.items():
        c=config(speed);mu,ref=fit(c,speed,target);mus[speed]=mu
        calibrations.append({'speed_kmh':speed,'published_distance_m':target,'calibrated_distance_m':ref['path_distance_m'],
                             'equivalent_parameter_not_measured_tyre_mu':mu,'calibrated_model_stop_time_s':ref['stop_time_s']})
        print('Calibrated',speed,'km/h:',mu,ref['path_distance_m'],flush=True)
        for mode in ['reference','mass-only','body-suction']:
            z,h=vehicle.run(c,mu,straight,mode,speed);z['energy_partition']=budget(z,h,c);cases.append(z)
            name=f'highspeed-{speed}-{mode}';histories[name]=h;save_rows(name+'.csv',h)
            print('Case',speed,mode,z['path_distance_m'],'m',flush=True)
    # Keep steering. No forced wheel deployment is credited at high speed.
    stress=[];speed=400;c=config(speed);mu=mus[speed]
    scenarios=[{'name':'400-crosswind-pulse','surface_factors':[1,1],'lateral_disturbance_N':1200},
               {'name':'400-curve-0p4g-initial','surface_factors':[1,1],'curvature_m_inv':.4*9.81/(400/3.6)**2,'initial_yaw_rate_rads':.4*9.81/(400/3.6)},
               {'name':'400-one-channel-loss','surface_factors':[1,1],'fan_channel_failure_s':1.1}]
    for s in scenarios:
        for mode in ['reference','body-suction']:
            z,h=vehicle.run(c,mu,s,mode,speed);stress.append(z);histories[s['name']+'-'+mode]=h
    sensitivities=[]
    for label,gap,start,extra_mass in [('half-mm-effective-gap',.00005,.15,170),
                                     ('three-mm-effective-gap',.0003,.15,170),
                                     ('one-second-fan-delay',.0001,1.,170),
                                     ('two-second-fan-delay',.0001,2.,170),
                                     ('340kg-hardware',.0001,.15,340)]:
        cc=copy.deepcopy(c);cc.update(deployment_start_s_assumed=start,installed_hardware_mass_target_kg=extra_mass)
        s={'name':label,'surface_factors':[1,1],'ground_map_overrides':{'nominal_seal_gap_m':gap}}
        z,_=vehicle.run(cc,mu,s,'body-suction',400,record=False);sensitivities.append(z)
    # Area is not given in the source sheet. Preserve every reference distance
    # under each area assumption, then show how parameter ambiguity affects gain.
    area_checks=[]
    for area in [1.6,2.4]:
        cc=config(400,area);mm,zz=fit(cc,400,491.)
        z,_=vehicle.run(cc,mm,straight,'body-suction',400,record=False)
        area_checks.append({'assumed_projected_area_m2':area,'refit_equivalent_parameter':mm,'reference_distance_m':zz['path_distance_m'],'prototype_projection_m':z['path_distance_m']})
    dt_checks=[]
    for mode in ['reference','body-suction']:
        coarse=next(z for z in cases if z['speed_kmh']==400 and z['mode']==mode)
        fine,_=vehicle.run(c,mu,straight,mode,400,dt=c['step_s']/2,record=False)
        dt_checks.append({'mode':mode,'distance_relative_change':fine['path_distance_m']/coarse['path_distance_m']-1,
                          'normal_load_residual_N':fine['max_normal_load_residual_N']})
    g=ground_config();uf=c['underfloor'];g.update(effective_total_area_m2=uf['bench_area_m2'],effective_total_perimeter_m=uf['bench_total_perimeter_m'],cavity_height_m=uf['bench_cavity_height_m'],feed_duct_diameter_m=uf['bench_feed_duct_diameter_m'],duct_loss_coefficient_assumed=uf['duct_loss_coefficient_assumed'])
    maps=[steady(x,g,10) for x in [.0005,.001,.003]]
    pressure_scales=[]
    for speed in [100,200,300,400]:
        v=speed/3.6;q=.5*c['air_density_kgm3']*v*v
        pressure_scales.append({'speed_kmh':speed,'ambient_dynamic_pressure_Pa':q,'illustrative_delta_Cp_0p1_pressure_Pa':.1*q,'Mach_at_343ms_assumed':v/343.,'kinetic_energy_relative_to_100kmh':(speed/100)**2})
    allcases=cases+stress+sensitivities
    assert all(z['stopped'] and z['max_force_circle_utilization']<=1+1e-8 and z['max_normal_load_residual_N']<1e-6 for z in allcases)
    assert all(abs(x['distance_relative_change'])<.002 for x in dt_checks)
    comparisons=[]
    for speed,target in BENCHMARKS.items():
        z=next(z for z in cases if z['speed_kmh']==speed and z['mode']=='body-suction')
        comparisons.append({'speed_kmh':speed,'published_Chiron_distance_m':target,
                            'conditional_prototype_distance_m':z['path_distance_m'],
                            'conditional_distance_reduction_m':target-z['path_distance_m'],
                            'conditional_distance_reduction_fraction':1-z['path_distance_m']/target,
                            'fan_energy_used_kJ_including_assumed_prespinnning':z['fan_energy_used_J']/1000,
                            'example_stability_limits_pass':z['meets_illustrative_stability_target']})
    result={'evidence':'PUBLISHED_BENCHMARKS_PLUS_CALIBRATED_HIGH_SPEED_FORCE_ENVELOPE_PROJECTIONS_NO_HARDWARE_TEST',
            'source_specifications':SPEC,'source_400kmh_record':RECORD,'published_distances_m':BENCHMARKS,
            'published_DIN_empty_mass_kg':1995,'actual_record_test_mass_known':False,
            'reference_configuration':config(400),'calibrations':calibrations,'comparisons':comparisons,
            'cases':cases,'stress_cases':stress,'sensitivity_cases':sensitivities,'projected_area_sensitivity_with_refit':area_checks,
            'time_step_checks':dt_checks,'installed_static_fan_map_assumptions':maps,'high_speed_pressure_scales':pressure_scales,
            'limits':[
                'Each starting-speed test has its own fitted aggregate parameter. These independent calibration points do not identify a common measured tyre, brake, fade or aero map, and are not independent model validation.',
                'The shared reference includes the published air-brake Cd and the stated rear-load increment. Frontal area, force location, deployment interpolation, CG and inertia are assumptions.',
                'Additional normal load is assumed to be exploitable by tyre/brake control. A real torque or thermal ceiling can remove that benefit; disc/pad/tyre temperatures and fade are not modelled.',
                'The assumed 170 kg installed mass, 1 mm effective seal gap, fan map and rapid pressure availability are not achieved hardware measurements.',
                'An 8 kJ rotor pre-spin allowance is charged against the 100 kJ store. It is an illustrative energy allowance; maintaining standby speed, motor peak power, inertia and valve operation remain unverified.',
                'At 400 km/h ambient dynamic pressure exceeds the assumed fan stall head. Vehicle-specific underfloor Cp, moving-road leakage, compressibility, wheel speed and flexible skirts require coupled high-speed analysis and tests.',
                'No 16.4% low-speed CFD drag benefit is transferred. Reference and prototype retain the same assumed external-airbrake map.',
                'The direct cold-only bench model does not supply this fan power. The high-speed gain is a powered-suction hypothesis, not proof of net cryogenic energy generation.',
                'The original high-speed road-wheel deployment concept is not credited. Tyres retain rolling steering contact in these comparisons.'
            ],'physical_measurements_for_CryoPulse':0}
    (R/'output/high-speed-results.json').write_text(json.dumps(result,indent=2,allow_nan=False)+'\n')
    save_rows('high-speed-summary.csv',comparisons)
    compact={name:[[round(float(r[k]),6) for k in ['time_s','X_m','Y_m','heading_rad','speed_ms','fan_force_N','fan_energy_remaining_J','aerodynamic_downforce_N']] for r in rows[::max(1,len(rows)//350)]] for name,rows in histories.items()}
    (R/'output/high-speed-trajectories.json').write_text(json.dumps(compact,separators=(',',':')))
    fig,axs=plt.subplots(1,2,figsize=(12,5),layout='constrained')
    speeds=list(BENCHMARKS);axs[0].plot(speeds,list(BENCHMARKS.values()),'o-',color='#203f4d',label='Published Chiron distances')
    axs[0].plot(speeds,[z['conditional_prototype_distance_m'] for z in comparisons],'o-',color='#187f92',label='Conditional powered-suction projection')
    axs[0].set(xlabel='Initial speed (km/h)',ylabel='Stopping path (m)',title='Published benchmark / conditional fitted-model comparison');axs[0].legend(fontsize=8);axs[0].grid(alpha=.2)
    for mode,col in [('reference','#203f4d'),('body-suction','#187f92')]:
        rows=histories['highspeed-400-'+mode];axs[1].plot([r['time_s'] for r in rows],[r['speed_ms']*3.6 for r in rows],label=mode,color=col)
    axs[1].set(xlabel='Time after braking command (s)',ylabel='Speed (km/h)',title='400 km/h speed trace in the force-envelope model');axs[1].legend(fontsize=8);axs[1].grid(alpha=.2)
    fig.suptitle('CryoPulse high-speed addendum / powered suction with rolling steering\nConditional model: actual Chiron stopping distances anchor each reference',fontsize=14)
    fig.savefig(R/'output/high-speed-comparison.png',dpi=170);plt.close(fig)
    print(json.dumps({'comparisons':comparisons,'sensitivities':[{k:z[k] for k in ['scenario','path_distance_m','fan_energy_used_J']} for z in sensitivities],'time_step_checks':dt_checks,'area_sensitivity':area_checks},indent=2),flush=True)


if __name__=='__main__':main()