"""Design calculations only. No experimental measurements are generated."""
from pathlib import Path
import csv, json, math
from statistics import NormalDist

HERE = Path(__file__).resolve().parent

def sound_speed(t_c, rh, pressure_pa, co2_ppm):
    """Cramer 1993 polynomial; RH is a fraction, actual local pressure in Pa.

    Equation independently transcribed from published coefficients. Applicability
    guards follow the documented Cramer range; no uncertainty is implied.
    """
    if not all(math.isfinite(v) for v in (t_c, rh, pressure_pa, co2_ppm)):
        raise ValueError('Finite environmental inputs required')
    if not (0 <= t_c <= 30 and 0 <= rh <= 1 and 75000 <= pressure_pa <= 102000 and 0 <= co2_ppm < 10000):
        raise ValueError('Outside declared Cramer environmental range')
    t=t_c; T=t+273.15; p=pressure_pa; xc=co2_ppm*1e-6
    enhancement=1.00062+3.14e-8*p+5.6e-7*t*t
    saturated=math.exp(1.2811805e-5*T*T-1.9509874e-2*T+34.04926034-6353.6311/T)
    xw=rh*enhancement*saturated/p
    if xw>.06: raise ValueError('Water-vapor mole fraction outside range')
    return (331.5024+.603055*t-.000528*t*t
      +(51.471935+.1495874*t-.000782*t*t)*xw
      +(-1.82e-7+3.73e-8*t-2.93e-10*t*t)*p
      +(-85.20931-.228525*t+5.91e-5*t*t)*xc
      -2.835149*xw*xw-2.15e-13*p*p+29.179762*xc*xc+.000486*xw*p*xc)

def pipe_mode(length_m, radius_m, mode, speed=343.0, eta=.6133):
    if not all(math.isfinite(v) for v in (length_m,radius_m,speed,eta)):
        raise ValueError('Finite geometry required')
    if min(length_m,radius_m,speed)<=0 or eta<0 or isinstance(mode,bool) or int(mode)!=mode or mode<1:
        raise ValueError('Positive dimensions and integer mode required')
    effective=length_m+2*eta*radius_m
    return dict(frequency_hz=mode*speed/(2*effective), ka=mode*math.pi*radius_m/effective,
                primary_low_ka_eligible=mode*math.pi*radius_m/effective<=.3)

def softening_from_frequency(f0, f, zeta0=0.0, zeta=0.0):
    """Positive D means reduced effective stiffness. f inputs are ringdown fd.

    Requires unchanged effective mass and calibrated fixture, with underdamped
    single-mode behavior. Not valid for a driven resonance peak without a model.
    """
    if not all(math.isfinite(v) for v in (f0,f,zeta0,zeta)) or min(f0,f)<=0 or not (0<=zeta0<1 and 0<=zeta<1):
        raise ValueError('Positive ringdown frequencies and 0 <= damping < 1 required')
    return 1-(f/f0)**2*(1-zeta0*zeta0)/(1-zeta*zeta)

def pilot_sample_size(sd, effect, margin=.02, power=.9, alpha=.05):
    """Normal approximation, independent equal-size groups, positive contrast.

    Targets the lower two-sided CI bound above the practical margin. Excludes
    attrition, multiple comparisons, repeated measures and batch effects.
    """
    if not all(math.isfinite(v) for v in (sd,effect,margin,power,alpha)) or not (sd>0 and effect>margin>=0 and 0<power<1 and 0<alpha<1):
        raise ValueError('Effect must exceed the practical margin; check inputs')
    z=NormalDist().inv_cdf(1-alpha/2)+NormalDist().inv_cdf(power)
    return math.ceil(2*sd*sd*z*z/(effect-margin)**2)

def pulse_sequence(tp, tr, tau_on, tau_off, pulses):
    """Exact scalar on/off relaxation, x=0 initially and equilibrium x=1 on.

    Values are normalized state at pulse ends, not measured current or a
    reproduction of a specified paper figure. Off equilibrium is zero.
    """
    if not all(math.isfinite(v) for v in (tp,tr,tau_on,tau_off)) or tp<=0 or tr<0 or min(tau_on,tau_off)<=0 or int(pulses)!=pulses or pulses<1:
        raise ValueError('Positive pulse/relaxation times and pulse count required')
    alpha=math.exp(-tp/tau_on); beta=math.exp(-tr/tau_off)
    a_star=alpha*(-math.expm1(-tr/tau_off))/(-math.expm1(-tp/tau_on-tr/tau_off))
    return [1-(a_star+(alpha-a_star)*math.exp((n-1)*(-tp/tau_on-tr/tau_off))) for n in range(1,pulses+1)]

def main():
    out=HERE/'design-output';out.mkdir(exist_ok=True)
    v=sound_speed(20,.5,101325,420); rows=[]
    for i,L in enumerate([.3,.4,.5,.6,.7,.8,.9,1.0],1):
        tube=f'T{i:02}'; split='test' if tube in ('T02','T05','T07') else 'calibration'
        for n in range(1,7):
            rows.append(dict(tube_id=tube,split=split,length_m=L,radius_m=.025,mode=n,**pipe_mode(L,.025,n,v)))
    with (out/'acoustic-design.csv').open('w',newline='') as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    report=dict(status='Design calculations; no collected measurements',sound_speed_example_m_s=v,
        environment=dict(temperature_c=20,relative_humidity_fraction=.5,pressure_pa=101325,co2_ppm=420),
        low_ka_limit=.3,primary_modes=sum(r['primary_low_ka_eligible'] for r in rows),
        candidate_modes=len(rows),independent_tubes=8,
        softening_example=dict(f0_hz=100,f_hz=95,damping_unchanged_zero=True,D=softening_from_frequency(100,95)),
        power_sensitivity=[dict(sd_fraction=sd,effect_fraction=effect,practical_margin=.02,n_per_group_normal_approx=pilot_sample_size(sd,effect)) for sd in (.01,.02,.04) for effect in (.03,.05)],
        pulse_example=dict(tp=.1,tr=.2,tau_on=.5,tau_off=1.0,end_state=pulse_sequence(.1,.2,.5,1.0,20),scope='Illustrative dimensionless state, times in seconds; not a source-paper parameter set'))
    (out/'design-results.json').write_text(json.dumps(report,indent=2)+'\n')
    print(json.dumps(report,indent=2))

if __name__=='__main__':main()