"""Hypothetical engineering-stress measurement study, not a material test."""
from pathlib import Path
import csv,json,math,random
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
OUT=Path(__file__).resolve().parent
SEED=20260924;N=10000;WIDTH=5.;THICKNESS=.5;LENGTH=25.;E=1200.;KINK=.025;TANGENT=40.;NOISE=.05

def imposed_stress(epsilon):return E*epsilon if epsilon<=KINK else E*KINK+TANGENT*(epsilon-KINK)
def write(name,rows):
    with (OUT/name).open('w',newline='',encoding='utf-8') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
rng=random.Random(SEED);area=WIDTH*THICKNESS;curve=[]
for i in range(101):
    strain=i/1000;stress=imposed_stress(strain);force=stress*area;reading=force+rng.gauss(0,NOISE)
    curve.append(dict(engineering_strain=strain,engineering_strain_percent=100*strain,nominal_gauge_extension_mm=strain*LENGTH,imposed_engineering_stress_MPa=stress,imposed_force_N=force,simulated_load_reading_N=reading,reconstructed_engineering_stress_MPa=reading/area))
errors=[];draws=[];reference=imposed_stress(.1)
for i in range(N):
    width=WIDTH*(1+rng.gauss(0,.02));thickness=THICKNESS*(1+rng.gauss(0,.02));assert min(width,thickness)>0
    estimated=reference*area/(width*thickness);err=estimated/reference-1
    errors.append(err);draws.append(dict(draw=i+1,assumed_measured_width_mm=width,assumed_measured_thickness_mm=thickness,inferred_stress_MPa=estimated,relative_stress_error=err))
errors.sort();quant=lambda p:errors[int(p*(N-1))]
write('engineering-curve.csv',curve);write('geometry-uncertainty.csv',draws)
# Direct force/area, extension/length identities; SI units provide independent scale checks.
max_si=max(abs(r['imposed_force_N']/(area*1e-6)/1e6-r['imposed_engineering_stress_MPa']) for r in curve)
assert max_si<1e-12
assert all(math.isclose(r['nominal_gauge_extension_mm']/LENGTH,r['engineering_strain'],abs_tol=1e-14) for r in curve)
assert imposed_stress(0)==0 and math.isclose(imposed_stress(KINK),30)
assert math.isclose(imposed_stress(.1),33)
assert all(math.isclose(d['relative_stress_error'],area/(d['assumed_measured_width_mm']*d['assumed_measured_thickness_mm'])-1,abs_tol=1e-14) for d in draws)
# Deterministic geometry bias: both dimensions 2% too large imply stress/(1.02²).
assert math.isclose((reference*area)/((WIDTH*1.02)*(THICKNESS*1.02))/reference,1/(1.02**2),rel_tol=1e-14)
# Independent closed-form first-order propagation predicts sqrt(.02²+.02²).
mean=sum(errors)/N;sd=math.sqrt(sum((x-mean)**2 for x in errors)/(N-1));linear_sd=math.sqrt(2)*.02
assert abs(sd/linear_sd-1)<.05
meta={'title':'Reading a tensile curve: geometry changes the reported stress','kind':'Seeded synthetic measurement demonstration; no physical specimen or measured material performance','seed':SEED,'random_generator':'Python random.Random MT19937, Gaussian gauss method','parameters':{'nominal_width_mm':WIDTH,'nominal_thickness_mm':THICKNESS,'nominal_gauge_length_mm':LENGTH,'assigned_initial_slope_MPa':E,'assigned_transition_strain':KINK,'assigned_post_transition_slope_MPa':TANGENT,'simulated_force_noise_SD_N':NOISE,'independent_dimension_relative_SD':.02},'formulas':{'engineering_stress':'sigma=F/A0 using original cross-sectional area','engineering_strain':'epsilon=DeltaL/L0','hypothetical_curve':'sigma=1200*epsilon up to 0.025; sigma=30+40*(epsilon-0.025) above','geometry_error':'sigma_inferred/sigma_imposed-1=A_nominal/A_assumed_measured-1'},'row_counts':{'curve':len(curve),'geometry_draws':len(draws)},'summary':{'chosen_curve_endpoint_MPa':reference,'endpoint_force_N':reference*area,'geometry_error_empirical_p2_5':quant(.025),'geometry_error_empirical_median':quant(.5),'geometry_error_empirical_p97_5':quant(.975),'geometry_error_sample_SD':sd,'first_order_predicted_SD':linear_sd},'validation':{'independent_SI_stress_conversion_max_error_MPa':max_si,'extension_strain_identity':True,'transition_continuity':True,'geometry_inverse_area_identity_all_draws':True,'first_order_SD_agreement_within5percent':True},'limitations':['All material-curve and noise parameters are assigned hypothetical inputs, not identified polymer properties.','Engineering stress uses original area. No instantaneous-area true-stress model is present.','33 MPa is the selected endpoint maximum, not ultimate tensile strength or fracture stress; no necking or failure is modelled.','The 2.5 percent kink is an assigned slope change, not an experimentally measured yield point. Unloading and permanent deformation are not defined.','Gaussian load noise and independent 2 percent dimensional errors are illustrative; real calibration, bias and correlation can differ.','The geometry histogram is a propagated assumption distribution, not a confidence interval for a real material.','A small Gaussian dimensional-error model can mathematically yield nonpositive dimensions; all 10,000 generated dimensions are explicitly checked positive.','The model does not determine toughness to fracture, ductility at failure, fatigue, rate/temperature dependence, bending properties or design allowables.'],'sources':[{'title':'NIST, A Tensile-Testing Technique for Micrometer-Sized Free-Standing Thin Films','url':'https://www.nist.gov/publications/tensile-testing-technique-micrometer-sized-free-standing-thin-films'},{'title':'NIST, Diffraction Stress Measurement Under Applied Load','url':'https://www.nist.gov/programs-projects/ncal-diffraction-stress-measurement-under-applied-load'}]}
(OUT/'study.json').write_text(json.dumps(meta,indent=2)+'\n',encoding='utf-8')
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False})
fig,ax=plt.subplots(1,2,figsize=(11.5,4.8),constrained_layout=True)
ax[0].plot([r['engineering_strain_percent'] for r in curve],[r['imposed_engineering_stress_MPa'] for r in curve],color='#393f94',lw=2,label='Assigned noiseless engineering curve')
ax[0].scatter([r['engineering_strain_percent'] for r in curve],[r['reconstructed_engineering_stress_MPa'] for r in curve],s=9,color='#b42e8d',alpha=.65,label='Simulated load reading / original area')
ax[0].axvline(2.5,color='#777',ls=':');ax[0].set(xlabel='Engineering strain (%)',ylabel='Engineering stress (MPa)',title='A  |  An imposed curve, not a tested material');ax[0].legend(fontsize=8,loc='lower right');ax[0].text(4,17,'Assigned slope change at 2.5%\nNo fracture or unloading model',fontsize=9)
ax[1].hist([x*100 for x in errors],bins=55,density=True,color='#087f8c',alpha=.8)
for v in [quant(.025),quant(.975)]:ax[1].axvline(v*100,color='#b42e8d',ls='--')
ax[1].set(xlabel='Stress error from assumed measured dimensions (%)',ylabel='Probability density per percentage point',title='B  |  Small dimension errors affect stress');ax[1].text(.04,.94,f'10,000 seeded geometry draws\nCentral 95%: {quant(.025)*100:.2f}% to {quant(.975)*100:.2f}%',transform=ax[1].transAxes,va='top',fontsize=9,bbox=dict(facecolor='white',alpha=.9,edgecolor='none'))
for a in ax:a.grid(alpha=.15)
fig.suptitle('Force is measured; stress also depends on specimen geometry',fontsize=15,weight='bold')
fig.supxlabel('All curves and error distributions are synthetic. Width 5 mm, thickness 0.50 mm, gauge 25 mm; dimension errors are independent 2% SD.',fontsize=9)
for ext in ['png','svg','pdf']:fig.savefig(OUT/f'tensile-measurement.{ext}',dpi=180)
print(json.dumps(meta['summary']))